edgeR差异分析

查看当前工作路径
> getwd()
[1] "C:/Users/YLAB/Documents"

重新设置工作路径
> setwd("E:/ZYH/R.project/rna-seq/lianxi1/exon_level/edgeR/")

载入数据
> raw_count <- read.table("E:/ZYH/R.project/rna-seq/lianxi1/exon_level/edgeR/merge.txt",sep="\t", header = T)

#先把第一列当作行名来处理
row.names(raw_count)<-make.names(raw_count[,1],TRUE)

#删除第一列
raw_count<-raw_count[,-1]

初步过滤，删除表达矩阵中含有0的行
raw_count= raw_count[which(rowSums(raw_count==0)==0),]

#指定分组，注意要保证表达矩阵中的样本顺序和这里的分组顺序是一一对应的
#对照组在前，处理组在后
group <- rep(c('control', 'treat'), each = 4)

载入R包
> library(limma)
> library(edgeR)
> library(statmod)
#数据预处理
#（1）构建 DGEList 对象
> dgelist <- DGEList(counts = raw_count, group = group)
#（2）过滤 low count 数据，例如 CPM 标准化（推荐）
keep <- rowSums(cpm(dgelist) > 1 ) >= 2
dgelist <- dgelist[keep, , keep.lib.sizes = FALSE]

#（3）标准化，以 TMM 标准化为例
dgelist_norm <- calcNormFactors(dgelist, method = 'TMM')

计算差异倍数列表：整体过程分为两步。
1.一是估计离散度。根据试验设计的样本分组，通过加权似然经验贝叶斯模型，估算每个基因表达的离散度，其代表了基因表达值在个体间的差异。
2.二是模型拟合。edgeR包中提供了多种计算差异基因的算法模型，之间略有不同，如下展示了两种。
#差异表达基因分析
#首先根据分组信息构建试验设计矩阵，分组信息中一定要是对照组在前，处理组在后
design <- model.matrix(~group)

#（1）估算基因表达值的离散度
dge <- estimateDisp(dgelist_norm, design, robust = TRUE)

#（2）模型拟合，edgeR 提供了多种拟合算法
#负二项广义对数线性模型
fit <- glmFit(dge, design, robust = TRUE)
lrt <- topTags(glmLRT(fit), n = nrow(dgelist$counts))

write.table(lrt, 'control_treat.glmLRT.txt', sep = '\t', col.names = NA, quote = FALSE)

#拟似然负二项广义对数线性模型
fit <- glmQLFit(dge, design, robust = TRUE)
lrt <- topTags(glmQLFTest(fit), n = nrow(dgelist$counts))

write.table(lrt, 'control_treat.glmQLFit.txt', sep = '\t', col.names = NA, quote = FALSE)

筛选差异表达基因：通过输出的差异表达倍数表，即可自定义阈值筛选差异表达基因了。例如这里以上述输出的负二项广义对数线性模型的结果为例，将它再次读入到R中，并根据|log2FC|≥2以及FDR<0.05定义显著变化的基因，以及通过“up”和“down”分别区分上、下调的基因。

#读取上述输出的差异倍数计算结果
gene_diff <- read.delim('control_treat.glmLRT.txt', row.names = 1, sep = '\t', check.names = FALSE)

#首先对表格排个序，按 FDR 值升序排序，相同 FDR 值下继续按 log2FC 降序排序
gene_diff <- gene_diff[order(gene_diff$FDR, gene_diff$logFC, decreasing = c(FALSE, TRUE)), ]

#log2FC≥2 & FDR<0.05 标识 up，代表显著上调的基因
#log2FC≤-2 & FDR<0.02 标识 down，代表显著下调的基因
#其余标识 none，代表非差异的基因
gene_diff[which(gene_diff$logFC >= 2 & gene_diff$FDR < 0.05),'sig'] <- 'up'
gene_diff[which(gene_diff$logFC <= -2 & gene_diff$FDR < 0.05),'sig'] <- 'down'
gene_diff[which(abs(gene_diff$logFC) <= 2 | gene_diff$FDR >= 0.05),'sig'] <- 'none'

#输出选择的差异基因总表
gene_diff_select <- subset(gene_diff, sig %in% c('up', 'down'))
write.table(gene_diff_select, file = 'control_treat.glmQLFit.select.txt', sep = '\t', col.names = NA, quote = FALSE)

#根据 up 和 down 分开输出
gene_diff_up <- subset(gene_diff, sig == 'up')
gene_diff_down <- subset(gene_diff, sig == 'down')

write.table(gene_diff_up, file = 'control_treat.glmQLFit.up.txt', sep = '\t', col.names = NA, quote = FALSE)
write.table(gene_diff_down, file = 'control_treat.glmQLFit.down.txt', sep = '\t', col.names = NA, quote = FALSE)