我对 R 相当陌生，最近开始使用它来分析一些微阵列数据。分析的总体目标是采用 DC2 并比较该人群中的 WT 与 KO 组。但是我在 limma 处理方面遇到了一些问题。使用 oligo 包处理数据后，我尝试使用 limma 创建用于分析的设计矩阵。这是我对 DC2 的 ExpressionSet 的工作流程:
pData(DC2)
index filename genotype cell_type
1 KO DC2 2 HP10.CEL KO DC2
2 KO DC2 3 HP11.CEL KO DC2
3 KO DC2 4 HP12.CEL KO DC2
1 WT DC2 10 HP7.CEL WT DC2
2 WT DC2 11 HP8.CEL WT DC2
3 WT DC2 12 HP9.CEL WT DC2
design <- model.matrix(~DC2$genotype)
design
(Intercept) DC2$genotypeWT
1 1 0
2 1 0
3 1 0
4 1 1
5 1 1
6 1 1
fit <- lmFit(DC2, design)
fit <- eBayes(fit)
toptable(fit)
这将输出基因列表，如下所示:
logFC t P.Value adj.P.Val B
17551163 14.09722 208.2627 2.990326e-13 2.700912e-10 17.14467
17511316 13.91167 205.0811 3.292503e-13 2.700912e-10 17.12716
17551167 13.92093 204.5801 3.343243e-13 2.700912e-10 17.12434
17375373 13.76320 202.1271 3.605170e-13 2.700912e-10 17.11025
17550685 13.74022 201.5428 3.671032e-13 2.700912e-10 17.10682
但是，当我使用以下代码手动检查(仅采用第一个功能)时:
toptable(fit, n=1)
genename <- rownames(toptable(fit, n=1))
typeMean <- tapply(exprs(DC2)[genename,], DC2$genotype, mean)
typeMean["KO"] - typeMean["WT"]
相同特征“17551163”的输出不同
KO
0.04538317
我试图四处寻找答案，但没有运气。我假设这可能与矩阵设计有关？任何帮助将不胜感激。
谢谢
请您参考如下方法:
对于那些跳过问题下方评论中的讨论的人来说，这是一个答案。
使用 lmFit 进行估算后和 eBayes我们可以质疑我们在 model.matrix 中提供的所有对比之间的最高区分基因。步。
在这里，作者设计如下:design <- model.matrix(~DC2$genotype) .请记住，(Intercept)如果我们想明确说明我们想要与 DC2$genotype 相关的对比度，则是第一个系数。 ，所以调用应该是:
toptable(fit, coef = 2)
自然地，如果设计包含更多对比，它们将被分配连续的自然数。
备注
如果我们想从设计中移除截距
design <- model.matrix(~ -1 + DC2$genotype) ;第一个系数现在是
DC2$genotype .
