运用limma对基因进行差异分析

此次案例我们是要比较ARID1A这个基因突变与不突变两亚类肝癌病人哪些基因的表达存在差异（寻找差异基因）。
在之前的三篇TCGA的数据采集文章中，我们已经采集到了所需要的所有数据，具体见下：

那么接下来我们就要对这些数据进行差异分析，采用的工具为R语言，所需要的包有limma、edgeR。

第一部分：安装及加载edgeR、limma包

edgeR的安装我们可以用先前的方法用Bioconductor进行安装。由于edgeR和limma互为依赖关系，所以limma也会在同时安装：

BiocManager::install(“edgeR”)
library(limma)
library(edgeR)
library(statmod)

第二部分：导入数据，设置分组，构建DGEList对象

此次样本数据中我们选择了83个肝癌病人的数据，其中56个属于未突变组，27个属于突变组。

targets <- as.matrix(read.csv("F:/公众号/图文素材/运用limma对基因进行差异分析/limma_data.csv", header = TRUE, row.names = 1))

group <- rep(c('C1', 'C2'), each = 56, length.out = 83) #C1是未突变组，C2是突变组

再根据基因表达量矩阵以及样本分组信息，构建DGEList对象。

dgelist <- DGEList(counts = targets, group = group)

第三部分：过滤低表达的基因，进行标准化处理

输入的基因表达量矩阵文件中，可能会含有一些低表达量的基因（甚至还有一些被忽略的全部为0值的行），需要在执行差异分析前将它们剔除。原因在于，这些基因未表达到具有生物学意义的程度（从生物学角度，有生物学意义的基因的表达量必须高于某一个阈值）；并且低表达量的基因受到随机因素影响比较大，故其统计结果也不可靠，还会影响p值校正过程；此外，过滤后的数据量降低，也能加快运行速率。
过滤的方法有很多种，在这里我们根据CPM值进行过滤，如果有对标准化方法有兴趣的小伙伴可以阅读“Comparing the normalization methods for the differential analysis of Illumina high-throughput RNA-Seq data”, 了解不同标准化方法之间的差异。

keep <- rowSums(cpm(dgelist) > 1 ) >= 2 #过滤
dgelist <- dgelist[keep, ,keep.lib.sizes = FALSE]

使用edgeR中的calcNormFactors()函数对数据标准化，以消除由于样品制备或建库测序过程中带来的影响。这里以TMM标准化方法为例：

dgelist_norm <- calcNormFactors(dgelist, method = 'TMM')  #TMM标准化

说明：
calcNormFactors()的方法有"TMM"、“RLE”、"upperquartile"三个，其它参数设置可参考：http://www.biostatistic.net/thread-10332-1-1.html

计算得到的归一化系数被用作样本大小的缩放系数。标准化后的DGEList对象中，增添了标准化因子信息。

然后，可以利用limma包中plotMDS()，绘制MDS图，使用无监督聚类方法展示出了样品间的相似性（或差异）。可据此查看各样本是否能够很好地按照分组聚类，评估试验效果，判别离群点，追踪误差的来源等。

plotMDS(dgelist_norm, col = rep(c('red', 'blue'), each = 5), dim = c(1, 2))

第四部分：估算离散值

负二项分布（negative binomial，NB）模型需要均值和离散值两个参数。首先需要根据样本分组，或者说根据实验设计、目的，构建一个分组矩阵。然后使用该分组矩阵，通过加权似然经验贝叶斯（weighted likelihood empirical Bayes）模型，估算每个基因的负二项离散Tagwise（即经验贝叶斯稳健离散值）、Common（即经验贝叶斯稳健离散值的均值）、Trended（即经验贝叶斯稳健离散值的拟合值）。
edgeR中，分组矩阵使用model.matrix()获得，并可以通过estimateDisp()估算离散值。edgeR具体说明：

design <- model.matrix(~group)    #构建分组矩阵
dge <- estimateDisp(dgelist_norm, design, robust = TRUE) #估算离散值

对于两个分组，分别标识为0和1，比方说上图中，对于各样本是否属于“C2分组”（groupC2），前56个样本标识为0意为不属于“C2分组”，后27个标识为1意为属于“C2分组”。

需要注意，标识好0、1类型后，后面的差异分析将以“分组1”的基因表达量相较于“分组0”是上调还是下调为准进行统计。因此在本示例中，后续分析获得的基因表达量上调/下调均为“分组C2”相较于“分组C1”而言的。实际的分析中，切记不要搞反了。

使用plotBCV()展示估算的离散值

plotBCV(dge) #作图查看

第五部分：差异分析

edgeR和limma包中提供了多种计算差异基因的方法，建立在不同模型的基础上，本文采用的为负二项式广义对数线性模型（edgeR）。

首先拟合负二项式广义对数线性模型（negative binomial generalized log-linear model），获取差异基因。这种方法大致可以这样理解，如果某个基因的表达值偏离这个分布模型，那么该基因即为差异表达基因。

使用edgeR包中的函数glmFit()和glmLRT()实现，其中glmFit()用于将每个基因的read count值拟合到模型中，glmLRT()用于对给定系数进行统计检验。

fit <- glmFit(dge, design, robust = TRUE)     #拟合模型
lrt <- glmLRT(fit)   #统计检验
topTags(lrt)

topTags(lrt)简要展示了分析结果，其中：logFC即log2转化后的 Fold Change值，但是要注意，这里不是简单的将基因的read count值直接对比，而是分别计算了基因在两组中的CPM值，然后据此计算的logFC；logCPM是log2转化后的CPM值；LR，似然比统计；PValue，差异表达的p值；FDR，FDR校正后的p值。

write.csv(topTags(lrt, n = nrow(dgelist$counts)), ' F:/公众号/图文素材/运用limma对基因进行差异分析/glmLRT.csv', quote = FALSE) #输出主要结果

dge_de <- decideTestsDGE(lrt, adjust.method = 'fdr', p.value = 0.05)  #查看默认方法获得的差异基因
summary(dge_de)

decideTestsDGE()可用于统计差异基因数量，屏幕输出了其默认值（供参考，大多数情况下我们还是优先根据Fold Change值以及p值等手动去筛选，而不会在意这个程序自己判断的数值），-1表示下调基因数量，1表示上调基因数量，0表示无差异基因数量。注意，对于这里的示例数据，基因表达量上调/下调均为“分组C2”相较于“分组C1”而言的。

plotMD()是limma包中的方法，可以初步绘制火山图观测差异基因分析结果。下图为程序默认的差异分析结果，对应了decideTestsDGE()统计的差异基因数量。纵轴为log2 Fold Change值；横轴为log2 CPM值，反映了基因表达量信息；蓝色的点表示上调基因，红色的点表示下调基因，黑色的点为无差异基因。我们也可以据此绘制火山图和热图

plotMD(lrt, status = dge_de, values = c(1, -1), col = c('blue', 'red'))#作图观测
abline(h = c(-1, 1), col = 'gray', lty = 2)

参考文章：
1、代码及解释主要参考于刘尧老师的文章：http://blog.sciencenet.cn/blog-3406804-1188480.html
2、edgeR包详细说明：http://www.bioconductor.org/packages/release/bioc/html/edgeR.html
3、Comparing the normalization methods for the differential analysis of Illumina high-throughput RNA-Seq data ：https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-015-0778-7

获取代码

以下是我的个人公众号，该篇的数据及代码可在公众号中回复“差异基因分析”即可获得，谢谢大家支持。