r-Limma使用makeContrasts和eBayes比较批量RNA序列

Question

经过一天的谷歌搜索，我决定最好在这里提问。

因此，实验是我有来自3名患者的大量RNA-seq数据：A、B、C。他们的RNA-seq数据是针对预处理、治疗周期1、治疗周期2、治疗周期3获得的。

所以我总共有12个批量RNA-seq:样本

• A.PreTreat->A.Cycle1->A.Cycle2->A.Cycle3

• B.PreTreat->B.Cycle1->B.Cycle2->B.Cycle3

• C.预处理->C.周期1->C.周期2->C.周期3

我想使用model.matrix(), lmFit(), makeContrasts(), contrasts.fit(), eBayes()获得不同周期(即周期3至预处理，周期3至周期2(之间的差异基因列表，所有这些都在limma包中。

这是我最简单的工作示例。

library(limma)
# Already normalized expression set: rows are genes, columns are the 12 samples
normalized_expression <- matrix(data=sample(1:100), nrow=10, ncol=12)
colnames(normalized_expression) <- c("A.PreTreat", "A.Cycle1", "A.Cycle2", "A.Cycle3", "B.PreTreat", "B.Cycle1", "B.Cycle2", "B.Cycle3", "C.PreTreat", "C.Cycle1", "C.Cycle2", "C.Cycle3")
patient_and_treatment <- factor(colnames(normalized_expression), levels = colnames(normalized_expression))
design.matrix <- model.matrix(~0 + patient_and_treatment)
colnames(design.matrix) <- patient_and_treatment
fit <- lmFit(normalized_expression, design.matrix)
# I want to get a contrast matrix to get differential genes between cycle 3 treatment and pre-treatment in all patients
contrast.matrix <- makeContrasts("A.Cycle3+B.Cycle3+C.Cycle3-A.PreTreat-B.PreTreat-C.PreTreat",
levels = levels(patient_and_treatment))
# Outputs Error of no residual degree of freedom
fit2 <- eBayes( contrasts.fit( fit, contrast.matrix ) )
# Want to run but cannot
summary(decideTests(fit2))

到目前为止，我还停留在没有剩余自由度的错误上。

我甚至不确定这是否是limma中统计学上正确的方法来解决我的问题，即在所有患者的第3周期治疗和预治疗之间获得差异基因列表。

任何帮助都将不胜感激。

谢谢！

Answer 1

每组不能有一个观察结果，这使得回归变得毫无意义，因为你要将每个数据点拟合到自己身上。

简单地说，你要寻找的是在所有患者中观察到的常见影响，例如Cycle3与PreTreat等，建立这样的模型：

library(limma)
metadata = data.frame(
Patient=gsub("[.][^ ]*","",colnames(normalized_expression)),
Treatment=gsub("^[A-Z][.]*","",colnames(normalized_expression))
)
Patient    Treatment
1        A PreTreat
2        A   Cycle1
3        A   Cycle2
4        A   Cycle3
5        B PreTreat
6        B   Cycle1
7        B   Cycle2
8        B   Cycle3
9        C PreTreat
10       C   Cycle1
11       C   Cycle2
12       C   Cycle3

现在指定模型矩阵，患者术语是为了说明患者之间起始水平的差异：

design.matrix <- model.matrix(~0 + Treatment+Patient,data=metadata)
fit <- lmFit(normalized_expression, design.matrix)
contrast.matrix <- makeContrasts(TreatmentCycle3-TreatmentPreTreat,
TreatmentCycle1-TreatmentPreTreat,levels=design.matrix)
fit2 = contrasts.fit(fit, contrast.matrix)
fit2 = eBayes(fit2)

你可以检查系数是否满足你的要求：

fit2$coefficients
Contrasts
TreatmentCycle3 - TreatmentPreTreat
[1,]                           -3.666667
[2,]                          -13.666667
[3,]                            1.666667
[4,]                          -40.666667
[5,]                           12.000000
[6,]                          -46.000000
[7,]                          -32.000000
[8,]                            4.666667
[9,]                           11.333333
[10,]                            5.666667
Contrasts
TreatmentCycle1 - TreatmentPreTreat
[1,]                           -11.33333
[2,]                           -19.33333
[3,]                           -27.33333
[4,]                           -42.33333
[5,]                            27.33333
[6,]                           -32.66667
[7,]                           -33.00000
[8,]                           -30.66667
[9,]                            46.00000
[10,]                            17.33333