专栏首页小白鱼的生统笔记R包limma作差异基因分析

R包limma作差异基因分析

R包limma差异基因分析

limma包最初设计用于分析芯片数据,后经扩展也可适用于RNA-seq数据。通过将read count数据转换为log2-counts per million(logCPM),估计均值-方差(mean-variance)关系以确定在线性建模之前每次观察的权重,之后将数据应用于线性建模,并通过经验贝叶斯统计差异基因。

本篇简介R语言limma包的差异基因分析流程。

下文中使用的示例数据及相关R代码的百度盘链接:

https://pan.baidu.com/s/1WWuJ6aE25QGXpfFzbU91VA

limma差异基因分析

R包安装

limma直接使用Bioconductor安装。

#Bioconductor 安装 limma
BiocManager::install('limma')
 
library(limma)

准备数据

网盘示例数据“gene.txt”,是一个基因表达数据矩阵,包含10列样本(可分为control和treat两组,每组各5个样本)共计1000行基因。

#基因表达矩阵
exprSet <- read.delim('gene.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
head(exprSet)

limma分析需要指定一试验设计矩阵,即样本分组,要保证基因表达矩阵中的样本名称顺序和分组顺序是一一对应的。

#试验设计矩阵
#注意要保证表达矩阵中的样本顺序和这里的分组顺序是一一对应的
group <- factor(rep(c('control', 'treat'), each = 5), levels = c('control', 'treat'))
design <- model.matrix(~0+group)
colnames(design) = levels(factor(group))
rownames(design) = colnames(group)
design

基因表达值标准化(voom标准化)

计数数据(read count)存在不可忽视的均值-方差关系,原始计数表现出随着计数大小的增加而增加的方差,而对数计数则通常表现为均值-方差趋势减小。

limma包中voom()函数用于标准化RNA-Seq或芯片数据。该函数将计数(所有计数加0.5以避免对数取零)转换为logCPM值,之后将logCPM值矩阵进行标准化。通过估计对数计数的均值-方差趋势,然后根据其预测方差为每个观察值分配权重,用于在线性建模过程中使用权重来调整异方差。

#voom 标准化,详见 ?norm
norm <- voom(exprSet, design, plot = TRUE)

数据标准化后,根据每个基因的平均log2CPM构建线性模型,计算残差,并拟合log2CPM与残差标准差平方根的关系,得到的趋势线(红色平滑曲线)可为样本和基因分配权重。

对于该趋势线,若左侧0起点处所示残差标准差明显偏高,或者0起点处出现上式趋势,则表明数据中存在较多的低表达(low counts数)基因。若觉得它们可能会对后续的计算带来较大的误差,不妨在标准化前手动过滤下。

数据执行标准化后,原本离散程度较大的reads counts矩阵的离散程度降低。

#标准化前后数据比较
par(mfrow = c(2, 2))
boxplot(exprSet)
plotDensities(exprSet)
boxplot(norm$E)
plotDensities(norm$E)

此外,对于voom()标准化功能,除了像上述这样指定原始的基因表达矩阵外,也可将其应用于edgeR的DGEList对象中实现数据标准化,因此limma也常和edgeR结合使用。前文在简介edgeR时,在最后就展示了一例voom()标准化DGEList的例子。

差异分析

接下来就是使用标准化后的基因表达数据,通过加权或广义最小二乘法对每个基因拟合线性模型,并通过经验贝叶斯计算出适度的t统计值、F统计值和差异表达值,最后获得差异基因分析结果。

#线性拟合,详见 ?lmFit
fit <- lmFit(norm, design, method = 'ls')
 
#确定比较的两组
#后续将计算标记为 1 的组相对于 -1 的组,基因表达值的上调/下调状态
contrast <- makeContrasts('treat-control', levels = design)
contrast
 
#使用经验贝叶斯模型拟合标准误差,详见 ?eBayes
fit2 <- contrasts.fit(fit, contrast)
fit2 <- eBayes(fit2)
 
qqt(fit2$t, df = fit2$df.prior+fit2$df.residual, pch = 16, cex = 0.2)
abline(0,1)
 
#p 值校正、提取差异分析结果,详见 ?topTable
diff_gene <- topTable(fit2, number = Inf, adjust.method = 'fdr')
head(diff_gene, 10)
 
write.table(diff_gene, 'gene_diff.txt', col.names = NA, sep = '\t', quote = FALSE)

后续自定义根据log2FC值以及校正后的p值等,自定义筛选差异基因就可以了。

火山图示例。

#例如这里根据 |log2FC| >= 1 & FDR p-value < 0.01 定义“差异”
diff_gene[which(diff_gene$logFC >= 1 & diff_gene$adj.P.Val < 0.01),'sig'] <- 'red'
diff_gene[which(diff_gene$logFC <= -1 & diff_gene$adj.P.Val < 0.01),'sig'] <- 'blue'
diff_gene[which(abs(diff_gene$logFC) < 1 | diff_gene$adj.P.Val >= 0.01),'sig'] <- 'gray'
 
log2FoldChange <- diff_gene$logFC
FDR <- diff_gene$adj.P.Val
plot(log2FoldChange, -log10(FDR), pch = 20, col = diff_gene$sig)
abline(v = 1, lty = 2)
abline(v = -1, lty = 2)
abline(h = -log(0.01, 10), lty = 2)

本文分享自微信公众号 - 小白鱼的生统笔记(gh_5f751e893315),作者:生信小白鱼 鲤小白

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2019-11-03

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 突破韦恩图数量限制,R包UpSetR集合可视化

    传统的维恩图通常很难展示>5组以上的集合情形,而花瓣图则丢失了很多信息。与之相比,UpSet绘图提供了一种有效的方法来可视化多个集合的交集,它通过R包UpSet...

    用户7585161
  • R语言绘制差异火山图示例

    前几篇分别简介了limma、edgeR、DESeq2、EBSeq的差异基因分析方法,对于两组数据间的差异表达基因的可视化,我们一般就以火山图的样式呈现。

    用户7585161
  • 多元因子分析(MFA)及其在R中实现

    多元因子分析(multiple factor analysis,MFA)是主成分分析(PCA)的扩展,旨在描述多组变量集的相互关系(Escofier and P...

    用户7585161
  • 使用Mfuzz进行转录组表达模式聚类分析

    Mfuzz是用来进行不同时间点转录组数据表达模式聚类分析的R包,使用起来非常方便,直接输入不同样本归一化后的counts或者FPKM及TPM值就可进行聚类。

    生信小王子
  • GSEA结果解读

    ES是GSEA最初的结果,反应全部杂交data排序后,在此序列top或bottom富集的程度。 ES原理:扫描排序序列,当出现一个功能集中的gene时,增加E...

    Y大宽
  • UPS宣布遭黑客入侵 客户数据存在泄露风险

    选择了UPS的消费者,你们的包裹被准时送达了么?但如果你的数据遭到了泄露,或许就不会感谢他们了。一份新出的报告表明,继多家零售商的系统被入侵之后,UPS也遭遇了...

    安恒信息
  • 给你的女神洗洗脸

    请大家想一下,当自己存储了好久女神的图片被噪声污染了,那是一种怎样的伤心欲绝的事情啊。但是,有matlab在。你什么都不用担心,matlab会用滤波的办法给图片...

    matlab爱好者
  • 【静态延迟加载】self关键字和static关键字的区别

    上面是一个很经典很普通的工厂模式代码,我们期望的是输出各自手机的品牌名称,但是结果显示的是父类中的品牌名称。这说明我们调用的 self 关键代表的是代码中它所在...

    码缘
  • 如何取消ajax请求的回调

    我们在开发过程中有时候会碰到这样的需求,连续发送多个ajax请求,请求个数大于等于2,后面的ajax请求发送时,如果前面的ajax请求还没有返回,就取消前面aj...

    挥刀北上
  • Java之static作用的全方位总结

     引用一位网友的话,说的非常好,如果别人问你static的作用;如果你说静态修饰 类的属性 和 类的方法 别人认为你是合格的;如果是说 可以构成 静态代码块,...

    小勇DW3

扫码关注云+社区

领取腾讯云代金券