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

2019年7月28日 196次阅读

我是R的新手,最近开始用它来分析一些微阵列数据.该分析的总体目标是采用DC2并比较该群体中的WT与KO组.但是我遇到了一些问题处理问题.在使用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进行估计之后,我们可以质疑我们在模型中提供的所有对比之间的最高分化基因.矩阵步骤.

在这里,作者创建了如下设计：design< – model.matrix(~DC2 $genotype).请记住,(拦截)是第一个系数,如果我们需要明确说我们想要与DC2 $基因型相关的对比,那么调用应该是：

toptable(fit, coef = 2)

当然,如果设计包含更多对比,则会为它们分配连续的自然数.

备注

如果我们想从设计设计中删除截距< – model.matrix(〜-1 DC2 $genotype);第一个系数现在是DC2 $基因型.

点赞