edgeR 包——学习

edgeR：differential analysis of sequence read count data

数据来源：RNA测序，基因表达系列分析（serial analysis of gene expression， SAGE）
CHIP测序等的read count 数据
分析范围：基因（gene），外显子（exon），转录本的差异表达
特点：使用经验贝叶斯算法，适用于基因特异性生物学变异的估计。
功能：对read count数据进行差异分析（differential expression analysis）。

image.png

官方介绍说明
两种检验方法：quasi-likelihood method，用于bulk RNA-Seq数据差异表达分析；
likelihood ratio test，用于单细胞转录组测序（single cell RNA-seq）以及没有重复的数据集。

分析过程

1. 读入数据

2. 创建EDGList——edgeR存储数据形式
edge R存储数据在一个叫做DEGList的类似于list的对象中。

 group <- c(1,1,2,2) #分组信息
 y <- DGEList(counts=x, group=group) #创建一个DEGlist，参数：counts数据（数据框或者矩阵），分组信息（一般为因子）

3. 过滤
在所有文库中，计数非常低的基因几乎不能提供差异表达的证据。
从生物学的角度来看，一个基因必须在一定程度上表达后才有可能被表达 ，从而被翻译成蛋白质或在生物学上具有重要意义。在进一步分析之前，这些基因应该被过滤掉。
根据经验法则，如果基因不能在所有样本中表达，应该被去除。我们可以设置自己的基因表达定义。 通常
一个基因需要在一个文库中有5-10个计数，才能被认为在该文库中表达 。我们需要使用CPM（counts/1000000）过滤，而不是在直接过滤（因为直接过滤不考虑样本之间文库大小的差异）。

 y$samples
group lib.size norm.factors
Sample1 1 10880519 1
Sample2 1 9314747 1
Sample3 1 11959792 1
Sample4 2 7460595 1
Sample5 2 6714958 1
#lib.size表示每个样本的文库大小或者测序深度

###我们可以使用下面的代码进行低表达基因过滤
design <- model.matrix(~ Patient + Treatment)  #分组矩阵建立
> keep <- filterByExpr(y, design=Treatment)
> y <- y[keep, , keep.lib.sizes=FALSE]

4. 数据标准化

• 在对数据进行标准化的时候，我们不需要考虑校正基因长度问题，因为对于每个RNA样本而言，基因长度对read counts有着相同的相对影响。
  对read counts 影响最显著的技术因素是测序深度（测序得到的碱基总量（bp）与基因组大小（Genome）的比值，它是评价测序量的指标之一）。
  calcNormFactors函数通过查找一组缩放因子来标准化文库的大小
  对于文库的大小，使大多数基因的样本之间的logFC最小。
  计算这些比例因子的默认方法使用修剪后的m值的平均值
  (TMM， trimmed mean of M-values)每对样品之间。 我们称其为原始文库大小的乘积
  和缩放因子的有效文库大小。**在所有下游分析有效文库大小将取代原始文库大小 **
  建议大多数RNA-Seq数据使用TMM，在这些数据中，大多数（超过一半）的基因被认为在任何一对样本之间没有差异表达。

y <- calcNormFactors(y)  #数据标准化
y$samples
group lib.size norm.factors
Sample1 1 10880519 1.17
Sample2 1 9314747 0.86
Sample3 1 11959792 1.32
Sample4 2 7460595 0.91
Sample5 2 6714958 0.83

DGEList的所有归一化系数都乘以1，确保有效库大小的几何平均值与原始库大小的几何平均值相同。
大小的几何平均数。归一化系数低于1，表明少数高计数的基因垄断了测序，导致其他基因的计数低于通常的文库大小。因此，库的大小将被缩小，类似于缩减该库中的计数。反之，一个高于1的系数会扩大规模，类似于降低计数的规模。
5.估计离散程度（Estimating dispersions）

• 单因子

y <- estimateDisp(y)
y <- estimateCommonDisp(y)
y <- estimateTagwiseDisp(y)

• 多因子
  输入数据，DEGList and design分组矩阵
  一个函数估计共同离散度，趋势离散度，标签离散度

y <- estimateDisp(y,design)

分别计算用以下函数

 y <- estimateGLMCommonDisp(y, design)

> y <- estimateGLMTrendedDisp(y, design)

> y <- estimateGLMTagwiseDisp(y, design)

**拟合采用的是广义线性模型（GLM）
6. 差异基因检验

 group <- factor(c(1,1,2,2,3,3))
 design <- model.matrix(~group) #构建分组矩阵
fit <- glmQLFit(y, design)  #广义线性模型 quasi-likelihood (QL) F-test拟合
qlf.2vs1 <- glmQLFTest(fit, coef=2)
 topTags(qlf.2vs1)

根据FC阈值筛选差异表达基因

• 结合logFC与pvalue共同选择差异表达基因

 fit <- glmQLFit(y, design)
 tr <- glmTreat(fit, coef=2, lfc=1)
 topTags(tr)

其他
DEseq还可以进行GO以及KEGG通路富集分析，and基因集检测，热图，聚类，CRISPR-CAS8,shRNA-seq screen analysis，可变性剪切以及差异甲基化分析

参考资料：edgeR: differential analysis of sequence read count data User's Guide (bioconductor.org)