基于转录序列定量表达值，识别差异表达变化的基因或其他非编码RNA分子的方法实际上非常多。 然而，目前DESeq2和edgeR是出现频率最高的两种方法。
DESeq2在上一篇文章中有介绍，本文将继续展示r包装edgeR的差异基因分析过程。 像DESeq2一样，edgeR作为一种广泛使用的r包，其身影在文献中很常见。 例如以下所示。
相关文献的说明
edgeR安装
edgeR可以直接使用Bioconductor安装，但还是非常简单。
#Bioconductor安装编辑器
要安装#install.packages('BiocManager”)，必须先安装biocmanager。 如果尚未安装，请执行此步骤
biocmanager :3360安装(edger ) )。
应用edgeR鉴定差异表达基因
edgeR利用经验贝叶斯估计和基于负二项模型的精确检验来确定差异基因，通过调节基因间的过度离散度，像费休精确检验一样，但利用适应过度方差数据的精确检验来评估各基因的差异表达。
以下是edgeR分析差异表达基因的一般过程。
2.1准备读取数据和文件
首先准备基因表达量矩阵。
本文的所有测试数据和r代码都可以在文末获得
“control_treat.count.txt”是6个测序样本的基因表达值矩阵，包括3个处理组(treat )和3个对照组(control )。 需注意在对照组前，处理组后。 第一列是基因名称，注意不能有重复值。
样本文件格式，行为基因，列为样本，内容为基因表达值
然后，将准备好的基因表达量矩阵读入r，预先定义样品组信息。
#基因表达矩阵的读取
targets
#指定组，使表现矩阵中的样本顺序和这里的组顺序一一对应
对照组为前，处理组为后
组
2.2表达值预处理与标准化
在进行差异表达基因分析之前，将输入的基因表达当前序列和组信息构建到edgeR的DGEList对象中，便于数据的保存和中间运算； 并对表达进行标准化，消除样品制备或文库测序过程的影响。
同时，需要手动过滤低表达量的基因。 为什么不需要在前面的DESeq2中手动过滤？ 你可能会怀疑。 这是因为，从设计DESeq2之初就将自动过滤步骤封装在函数中，可以进行低表达基因的自动处理； edgeR中没有，需要额外的操作。 建议根据CPM(count-per-million，每百万个碱基目标基因的read count值)的值进行过滤。 例如，以CPM值为1为基准。 也就是说，某个基因的read count最低的样本(文库)中的count值)，) read count大于最低的样本的count总数
库(edger )
#数据预处理
#(1)构建DGEList对象
dgelist
(2)过滤CPM标准化(推荐)等低计数数据
keep 1)=2
dgelist
#3)标准化，以TMM标准化为例
dgelist_norm
2.3计算差异倍数列表
整个过程分为两个阶段。
一个是估计离散度。 根据实验设计的样本组，通过加权似然经验贝叶斯模型，估算了每个基因表达的离散度，揭示了个体间基因表达值的差异。
二是模型拟合。 edgeR包提供了几种计算差异基因的算法模型，之间略有不同。 以下是两种类型。
#差异表达基因分析
#首先，基于组信息构建试验设计矩阵，组信息中必有对照组在前，处理组在后
设计
(1)估算基因表达值的偏差
dge
(2)模型拟合，edgeR给出了多种拟合算法
#负二项广义对数线性模型
fit
lrt
write.table(lrt，' control_treat.glmLRT.txt '，sep='t '，col.names=NA，quote=FALSE ) ) )。
拟似然负二项广义对数线性模型
fit
lrt
write.table(lrt，' control_treat.glmQLFit.txt '，sep='t '，col.names=NA，quote=FALSE ) ) )。
最后输出了各基因的差异表达倍数表。
以负二项广义对数线性模型的结果为例，包括基因id、log2变换后的差异倍数(Fold Change )值、显著的p值和校正后的p值(缺省FDR校正)等主要信息。
差异表达基因分析
结果2.4 筛选差异表达基因
通过输出的差异表达倍数表，即可自定义阈值筛选差异表达基因了。
例如这里以上述输出的负二项广义对数线性模型的结果为例，将它再次读入到R中，并根据|log2FC|≥1以及FDR<0.01定义显著变化的基因，以及通过“up”和“down”分别区分上、下调的基因。
##筛选差异表达基因
#读取上述输出的差异倍数计算结果
gene_diff
#首先对表格排个序，按 FDR 值升序排序，相同 FDR 值下继续按 log2FC 降序排序
gene_diff
#log2FC≥1 & FDR<0.01 标识 up，代表显著上调的基因
#log2FC≤-1 & FDR<0.01 标识 down，代表显著下调的基因
#其余标识 none，代表非差异的基因
gene_diff[which(gene_diff$logFC >= 1 & gene_diff$FDR < 0.01),'sig']
gene_diff[which(gene_diff$logFC <= -1 & gene_diff$FDR < 0.01),'sig']
gene_diff[which(abs(gene_diff$logFC) <= 1 | gene_diff$FDR >= 0.01),'sig']
#输出选择的差异基因总表
gene_diff_select
write.table(gene_diff_select, file = 'control_treat.glmQLFit.select.txt', sep = 't', col.names = NA, quote = FALSE)
#根据 up 和 down 分开输出
gene_diff_up
gene_diff_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)
标识显著上下调的基因，添加在最后一列
新获得的3个表格中，分别包含了根据自定义阈值筛选的所有差异基因、上调基因或下调基因列表，最后一列标注了上下调概况。
2.5 差异表达火山图示例
已经识别了显著差异表达的基因，最后期望通过火山图将它们表示出来，火山图是文献中常见统计图表之一。
ggplot2为此提供了优秀的作图方案，参考以下示例。
##ggplot2 差异火山图
library(ggplot2)
#默认情况下，横轴展示 log2FoldChange，纵轴展示 -log10 转化后的 FDR
p
geom_point(size = 1) + #绘制散点图
scale_color_manual(values = c('red', 'gray', 'green'), limits = c('up', 'none', 'down')) + #自定义点的颜色
labs(x = 'log2 Fold Change', y = '-log10 FDR', title = 'control vs treat', color = '') + #坐标轴标题
theme(plot.title = element_text(hjust = 0.5, size = 14), panel.grid = element_blank(), #背景色、网格线、图例等主题修改
panel.background = element_rect(color = 'black', fill = 'transparent'),
legend.key = element_rect(fill = 'transparent')) +
geom_vline(xintercept = c(-1, 1), lty = 3, color = 'black') + #添加阈值线
geom_hline(yintercept = 2, lty = 3, color = 'black') +
xlim(-12, 12) + ylim(0, 35) #定义刻度边界
p
差异表达基因火山图