用于分析微阵列基因表达数据的 limma 问题 - 可能与设计矩阵有关答案

【问题标题】：Issues with limma for analysis of microarray gene expression data - possibly related to design matrix用于分析微阵列基因表达数据的 limma 问题 - 可能与设计矩阵有关
【发布时间】：2018-04-08 21:13:56
【问题描述】：

我对 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

我试图四处寻找答案，但没有运气。我假设这可能与矩阵设计有关？任何帮助将不胜感激。

谢谢

【问题讨论】：

您能否发布“17551163”的表达式值？
@PoGibas 感谢您的回复。表达值为：“1 KO DC2”14.0912151923856“2 KO DC2”14.0912151923856“3 KO DC2”14.14.092351251855“1 WT DC2”14.0592734145348“2 WT DC2”14.03271978068“3 WT DC2”14.0635228171713 span>
尝试将coef 参数放入topTable，以声明您感兴趣的对比度：topTable(fit, coef = 2)。这会改变输出吗？
@storaged 啊，太好了。这似乎奏效了——非常感谢。只是从概念上讲，coef = 2 位是否指定了DC2$genotypeWT 列中 0 和 1 之间的比较？那么当我设置coef = 1 时，会将(Intercept) 中的1 与DC2$genotypeWT 列中的0 和1 进行比较？
@Jack 据我所知limma 包它就是你描述它的方式。 coef = 2 明确表示我们希望在 DC2$genotypeWT 内进行比较。很高兴它对你有用！

标签： r bioinformatics bioconductor limma

【解决方案1】：

答案，对于那些跳过阅读问题下方 cmets 讨论的人。

在使用lmFit 和eBayes 进行估计之后，我们可以质疑我们在model.matrix 步骤中提供的所有对比之间的最高区分基因。

在这里，作者创建了如下设计：design <- model.matrix(~DC2$genotype)。请记住，(Intercept) 是第一个系数，如果我们需要明确表示我们想要与DC2$genotype 相关的对比度，那么调用应该是：

toptable(fit, coef = 2)

当然，如果设计包含更多对比，它们会被分配连续的自然数。

备注

如果我们想从设计中删除截距design <- model.matrix(~ -1 + DC2$genotype)；第一个系数现在是DC2$genotype。

【讨论】：