edgeR的简单实用流程

前些天，实验室里的同事，为了补充计算几个比较组合的差异分析，花了几千元，真够冤枉的。但碍于我们实验室对生信技术的鄙夷，我也懒得理睬。这里简单介绍一下如何使用edgeR进行差异基因的计算。文末附上完整代码，直接复制粘贴更改名称即可使用。

1. 安装

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("edgeR")

2. 读取数据

使用edgeR计算差异基因，使用的是readcount数据。我们这里使用最简单的双处理三重复数据进行举例。

OS = read.table("readcounts.txt",sep = "\t",header = T)
rownames(OS)=OS$GeneID
OS = OS[,2:7]
head(OS)

共两个处理：MM和CM，查看输入格式如下：

数据格式

3.构建DGEList

Treat = factor(substring(colnames(OS),1,2))
Treat = relevel(Treat, ref = "MM")
Time = factor(substring(colnames(OS),3,3))
#构建DEGList
y = DGEList(counts = OS, group = Treat)

Treat是指处理，利用substring(object, start, end),来指定处理名称，及去掉数字的部分，这里使用start和end确定处理的名字的位置，两位字节，即MM和CM。ref来确定对照。Time是指重复，这里的start和end则定位重复的标记，即1,2和3。这里MM1, CM1等，数字出现在第三位，因此start=3, end=3。
使用DEGList直接构建。

4. 过滤和标准化

部分readcount太低的需要去除。

keep = filterByExpr(y)
table(keep)
y = y[keep,,keep.lib.sizes = FALSE]
y = calcNormFactors(y)
y$samples

使用filterByExpr()来过滤，过滤掉372个基因。结果如下：

过滤

使用calcNormFactors()进行标准化，结果如下：

标准化

5. 设计矩阵和估计离散值

我们直接用GLM方法来计算

design = model.matrix(~Time+Treat)
rownames(design) = colnames(y)

设计的矩阵如下：
接着计算离散值：

y = estimateDisp(y, design, robust=TRUE)
y$common.dispersion
plotBCV(y)

BCV

fit = glmQLFit(y, design, robust = TRUE)
plotQLDisp(fit)

image.png

6. 计算差异值

qlf = glmQLFTest(fit)
topTags(qlf)

glmQLFTest()计算差异值，topTags()默认显示最显著的10个基因。TopTags()结果如下：

展示FDR最低的十个基因

top = rownames(topTags(qlf))
cpm(y)[top,]

可以得到各个样本的readcount中间值：

readcount中间值

summary(decideTests(qlf))
plotMD(qlf)
abline(h = c(-1,1), col = "blue")

最后展示上下调基因数并作图：

上下调基因概述

FDR<0.05的基因

7. 提取结果

output = topTags(qlf, n =9000)     n大于基因总数就行
top2 = rownames(output)
cpm2 = cpm(y)[top2,]
output2 = as.data.frame(output)
cpm3 = as.data.frame(cpm2)
output_all = cbind(cpm3,output2)
output_DEG = output_all[(output_all$FDR < 0.05 & output_all$logFC > 1) | (output_all$FDR < 0.05 & output_all$logFC < -1),]
output_Up = output_all[output_all$FDR < 0.05 & output_all$logFC > 1,]
output_Down = output_all[output_all$FDR < 0.05 & output_all$logFC < -1,]

write.table(output_all, file = "All.txt", row.names = T, sep = "\t", quote = F)
write.table(output_DEG, file = "DEG.txt", row.names = T, sep = "\t", quote = F)
write.table(output_Up, file = "Up.txt", row.names = T, sep = "\t", quote = F)
write.table(output_Down, file = "Down.txt",, row.names = T, sep = "\t", quote = F)

这里分别提取所有基因，差异基因，上调基因和下调基因的readcount中间值，FC和FDR值。

完整代码

最后附上所有代码：

setwd("...")
library(edgeR)

OS = read.table("readcounts.txt",sep = "\t",header = T)
rownames(OS)=OS$GeneID
OS = OS[,2:7]
head(OS)

Treat = factor(substring(colnames(OS),1,2))
Treat = relevel(Treat, ref = "MM")
Time = factor(substring(colnames(OS),3,3))
y = DGEList(counts = OS, group = Treat)

keep = filterByExpr(y)
table(keep)
y = y[keep,,keep.lib.sizes = FALSE]
y = calcNormFactors(y)
y$samples

design = model.matrix(~Time+Treat)
rownames(design) = colnames(y)
design

y = estimateDisp(y, design, robust=TRUE)
y$common.dispersion
plotBCV(y)
fit = glmQLFit(y, design, robust = TRUE)
plotQLDisp(fit)

qlf = glmQLFTest(fit)
topTags(qlf)
top = rownames(topTags(qlf))
cpm(y)[top,]
summary(decideTests(qlf))
plotMD(qlf)
abline(h = c(-1,1), col = "blue")

output = topTags(qlf, n =100000)   
top2 = rownames(output)
cpm2 = cpm(y)[top2,]
output2 = as.data.frame(output)
cpm3 = as.data.frame(cpm2)
output_all = cbind(cpm3,output2)
output_DEG = output_all[(output_all$FDR < 0.05 & output_all$logFC > 1) | (output_all$FDR < 0.05 & output_all$logFC < -1),]
output_Up = output_all[output_all$FDR < 0.05 & output_all$logFC > 1,]
output_Down = output_all[output_all$FDR < 0.05 & output_all$logFC < -1,]
write.table(output_all, file = "All.txt", row.names = T, sep = "\t", quote = F)
write.table(output_DEG, file = "DEG.txt", row.names = T, sep = "\t", quote = F)
write.table(output_Up, file = "Up.txt", row.names = T, sep = "\t", quote = F)
write.table(output_Down, file = "Down.txt",, row.names = T, sep = "\t", quote = F)

参考：Bioconductor-edgeR