多分组的差异分析只需要合理设置design矩阵即可

通常情况下，我们是做正常组和对照组差异分析，或者药物处理前后，这样的分组都是超级简单的。

如果加入药物加上浓度，梯度会稍微复杂一点。比如针对cisplatin这个药物的数据集：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE116439   ，几百个转录组芯片样品数据，其中一个细胞系：

GSM3232610 A498_cisplatin_0nM_24h
GSM3232611 A498_cisplatin_0nM_2h
GSM3232612 A498_cisplatin_0nM_6h
GSM3232613 A498_cisplatin_15000nM_24h
GSM3232614 A498_cisplatin_15000nM_2h
GSM3232615 A498_cisplatin_15000nM_6h
GSM3232616 A498_cisplatin_3000nM_24h
GSM3232617 A498_cisplatin_3000nM_2h
GSM3232618 A498_cisplatin_3000nM_6h

就是不同浓度的cisplatin药物处理不同时间，两个变量的不同梯度。这样的实验设计的差异分析就会比较复杂。

尤其是数据挖掘领域，有一个强行找差异的套路，就是按照基因表达量高低把肿瘤病人分组后，然后做差异分析这样的基因表达量高低分组，极端情况下其实是可以把两万个基因都做一遍。不过，大概率大家应该是并不会这样弄，比较喜欢的是，挑选某个生物学意义的基因集，比如铁死亡，缺氧等等。每个基因集里面都是几十甚至上百个基因，他们都可以根据表达量高低对病人进行分组，然后做差异分析。

我们来模拟一下这样的情况，然后给出一个合理的差异分析解决方案。

构造 design矩阵

如下所示的代码，模拟出来20个基因在50个样品的表达量随机数值，然后仅仅是在各个基因内部根据表达量高低进行  0和1的二值化，构成design矩阵：

tmp=as.data.frame(matrix(rnorm(1000),50))
cg='SF3B1|TET2|SRSF2|ASXL1|DNMT3A|RUNX1|U2AF1|TP53|EZH2|IDH2|NPM1|MLL2|PTPN11|CREBBP|KIT|MPL|NF1|WT1|IRF1|RAD21'
colnames(tmp)=strsplit(cg,'\\|')[[1]]
tmp[tmp>0]=1
tmp[tmp <0]=0
design= cbind(offset=1,tmp)
minF=5 ## Minimal number of alterations
design = design[,colSums(design)>=minF]
head(design[,1:5])

可以看到，design矩阵 如下所示：

> head(design[,1:5])
  offset SF3B1 TET2 SRSF2 ASXL1
1      1     1    0     1     0
2      1     0    0     1     0
3      1     0    0     0     1
4      1     1    1     1     0
5      1     1    1     1     0
6      1     1    0     0     0

因为design里面有offset，所以无需使用 makeContrasts 函数去 设置contrast.matrix矩阵。

可以看到，我们的20个基因，都是可以根据其在对应的50个样品里面的表达量差异分成高低两个分组， 后面就可以做20次差异分析。

创造模拟的表达量矩阵

这里面的表达量矩阵我们仍然是模拟的，还是选择rnorm函数：

geneExpr=as.data.frame(matrix(rnorm(1000),20))
head(geneExpr[,1:5])
dim(geneExpr)

如下所示，可以看到，是 50个样品的表达量矩阵， 每个样品里面都是20个基因的表达量。

> head(geneExpr[,1:5])
          V1         V2         V3         V4         V5
1 -0.6973356 -0.3119271  0.1975896 -0.9088014  1.8884890
2  0.8279720  1.8873989  1.1911934  0.3008998  1.2488066
3  0.9759161  0.3792309  0.2255550 -0.6953532 -0.1491860
4  1.0270041  1.3005796  1.7677480  1.0649399  0.3224370
5 -0.5429206  0.2575444 -0.8916706 -0.4578587 -1.8454486
6 -1.3586660  0.9819655  0.4961708  0.6441539 -0.6976303
> dim(geneExpr)
[1] 20 50

###整体的差异分析

使用 limma 包做差异分析，其实就  lmFit，eBayes，topTable 即可

前面的 design矩阵可以是无限制的列，每一列都是一个分组的可能性，目前仅仅是支持二分组哦，所以是0和1的二值化。

library(limma)
glm = lmFit(geneExpr , design = design ) 
glm = eBayes(glm)

diff1=topTable(glm, coef=2, n=Inf)
diff2=topTable(glm, coef=3, n=Inf)
head(topTable(glm, coef=ncol(design), n=Inf))

可以使用 topTable 函数取出任意一次差异分析结果表格，每个结果都可以去判断统计学显著的上下调基因，去富集分析，去做火山图，热图。这样的分析看我六年前的表达芯片的公共数据库挖掘系列推文即可哈 ：

一个思考题

这样的多次差异分析，每个都有自己的统计学显著的上下调基因，这样的话火山图，热图太多了，并不是一个很好的展现方式。

怎么样把前面的 topTable 函数得到 得到 diff1，diff2等数据框整合起来进行展示呢？

比如如下所示的两个火山图，就展示起来很耗费PPT了，如果20个基因的分组后的 20次差异分析结果，展现就很困难。
（火山图来源于昨天的学徒作业: 上下调基因各自独立进行GO数据库的3分类富集(求美图代码)）

GSE97251复现，陈杭 2021-7-15 

如果你确实想不到，明天我们的教程，你一定很需要！

文末友情推荐