使用limma进行两组间的差异分析

生信修炼手册

发布于 2020-05-08 08:46:51

6.9K015

代码可运行

文章被收录于专栏：生信修炼手册生信修炼手册

运行总次数：15

代码可运行

欢迎关注”生信修炼手册”!

limma这个R包可以用于分析芯片数据，也可以分析NGS测序的数据，其核心是通过线性模型去估算不同分组中基因表达量的均值和方差，从而进行差异分析。

limma也是基于raw count的定量方式，但是它并不提供归一化的算法。在官方手册中，推荐采用edgeR的TMM归一化算法。完整代码如下

1. 读取文件

读取基因在所有样本中的表达量文件，示例如下

gene_id ctrl-1 ctrl-2 ctrl-3 case-1 case-2 case-3
geneA 14  0  11  4  0  12
geneB 125 401 442 175 59 200

每一行为一个基因，每一列代表一个样本。读取数据的代码如下

# 读取表达量的表格
counts <- read.table(
 "gene.counts.tsv",
 header=T,
 sep="\t",
 row.names=1,
 comment.char="",
 check.names=F)# 设置样本分组
group <- factor(rep(c("control", "case"), each = 3))
design <- model.matrix(~group)# 构建edgeR中的对象
library(edgeR)
y <- DGEList(counts=count)

之所以采用edgeR来读取数据，是为了方便后续的预处理和归一化。

2. 过滤count数很低的基因

和edgeR中的预处理过程类似，根据CPM表达量对基因进行过滤，代码如下

keep <- rowSums(cpm(y)>1) >= 2
y <- y[keep, , keep.lib.sizes=FALSE]

3. 归一化

默认采用TMM归一化算法，计算每个样本的 sizefactor, 代码如下

y <- calcNormFactors(y)

4. 表达量转换

在进行差异分析前，需要对表达量进行转换，有以下两种选择

logCPM
voom

第一种转换就是计算logCPM值，第二种转换适用于样本间sizaFactors差异较大的情况。转换的代码如下

# logCPM
logCPM <- cpm(dge, log=TRUE, prior.count=3)
# voom
v <- voom(dge, design, plot=TRUE)

5. 差异分析

转换之后的表达量就可以进行差异分析了，代码如下

fit <- lmFit(logCPM, design)
fit <- eBayes(fit, trend=TRUE)
res<- topTable(fit, coef=ncol(design))

上述代码采用的是logCPM值，当然也可以采用voom转换后的值，当采用voom转换时，注意trend参数为FALSE。

这里只是介绍了最简单的用法，更多复杂案例，比如多个分组，时间序列的差异分析等，请参考官方文档。

·end·

—如果喜欢，快分享给你的朋友们吧—

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2018-10-03，如有侵权请联系 cloudcommunity@tencent.com 删除

黑石物理服务器

编程算法

本文分享自生信修炼手册微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体同步曝光计划，欢迎热爱写作的你一起参与！

黑石物理服务器

编程算法

登录后参与评论

暂无评论