运用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

获取代码

以下是我的个人公众号,该篇的数据及代码可在公众号中回复“差异基因分析”即可获得,谢谢大家支持。
在这里插入图片描述

  • 25
    点赞
  • 196
    收藏
    觉得还不错? 一键收藏
  • 1
    评论
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值