用makeContrasts和eBayes比较大体积核糖核酸的比较

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.Cycle1 -> B.Cycle1 -> B.Cycle1 en28 20 en23 en25 en27

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

这是我最小的工作例子。

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))

到目前为止，我没有剩余的自由度错误。

我甚至不确定这是否是在统计学上正确的方法，以解决我的问题，获得不同的基因之间的第三周期治疗之间的治疗前，在所有患者。

任何帮助都将不胜感激。

谢谢!

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