edgeR：差异表达分析

edgeR的适用性

适用于RNA-Seq，SAGE-Seq，Chip-Seq，CRISPR-Cas9，DNA methylation研究。

快速入门

glm approach 相比经典方法更灵活。旗下包含quasi-likelihood F-test method 和 likelihood ratio。
quasi-likelihood： 建议用于大量RNA-seq数据的差异表达分析。
likelihood ratio： 建议用于单细胞RNA-seq分析和没有重复的实验分析中。

x = read.delim("TableOfCounts.txt",row.name = "Symbol")
group = factor(c(1,1,2,2))
y = DGEList(counts = x,group = group)
y = calcNormFactors(y)
design = model.matrix(~group)
y = estimateDisp(y,design)

运行quasi-likelihood F test：
fit = glmQLFit(y,design)
qlf = glmQLFTest(fit,coef = 2)
topTags(qlf)

运行likelihood ratio tests：
fit = glmFit(y,design)
lrt = glmLRT(fit,coef = 2)
topTags(lrt)

edgeR 介绍

• DGEList Object
  DGEList(counts = x,group = group)
  DGEList有一个count matrix记录数值型计数、df sample包含样本和文库（lib.size列）信息、df gene记录了基因和基因区间注释信息。

• 过滤
  通常一个基因必须在一个文库中具有5-10的计数才能被认为在该基因中表达。
  应以CPM进行过滤。

keep = filterByExpr(y)
y = y[keep,,keep.lib.sizes = FALSE]

• 标准归一化

1. 测序深度
2. RNA比例：calNormFactors函数通过缩放比例来标准化RNA。
3. GC 含量
4. 基因长度
5. 基于模型的标准化

• 负二项模型
  BCV:生物变异系数

• 两个或多个组之间的成对比较（classic）

1. 估计离散度

估计qCML共有离散和标记离散
y = estimateDisp(y)
估计qCML共有离散
y = estimateCommonDisp(y)
估计qCML标记离散
y = estimateTagwiseDisp(y)

2. 检测DE基因

et = exactTest(y)
togTags(et)

• 更复杂的实验（glm）

3. 估计离散度（公共、趋势、标记）
   y = estimateDisp(y,design)
4. 检测DE基因

group = factor(c(1,1,2,2,3,3))
design = model.amtrix(~group)
fit = glmQLFit(y,design)

coef = 2:2vs1 coef = 3:3vs1 
contrast = c(0,-1,1) : 3vs2
qlf.2vs1 = glmQLFTest(fit,coef = 2 )

• 如果没有实验重复怎么办？

1. 根据经验选择一个合理的离散值，将其运用到exactTest或glmFit中。
   常见BCV的典型值为：人类数据为0.4，基因相同的模式生物数据为0.1，技术复制为0.01。

bcv = 0.2
y = DEGList(counts = counts,group = 1:2)
y = exactTest(y,dispersion = bcv ** 2)

2. 如果存在相当数量的control组转录本不是DE，可利用他们估计离散值。

创建只有一个治疗组的副本
y1 = y
y1$sample$group = 1

利用housekeeping基因估计
y0 = estimateDisp(y1[housekeeping,],trend = "none",tagwise = FALSE)
y$common.dispersion = y0$common.dispersion
fit = glmFit(y,design)
lrt = glmLRT(fit)

• 富集分析和通路分析
  goana()函数进行GO分析， kegga()函数进行KEGG路径分析。
  两个函数都接受一个DGELRT或DGEExact对象。
  输入对象的行名为Entrez ID。
  topGO函数可以查看最丰富的顶部集合，topKEGG函数可以查看路径的顶部集合。

qlf = glmQLFTest(fit,coef = 2)
go = goana(qlf,species = "Mm"，ont = "BP")
topGo(go,sort = "up")

keg = kegga(qlf, species="Mm")
topKEGG(keg, sort="up")

• 聚类和热图
  使用cpm()函数：y是标准化的DGEList对象
  logcpm = cpm(y,log = TRUE)

数据挖掘

• MDS图
  X轴为样本之间的距离，Y轴为相对应的BCV。
  plotMDS(y)
  

• MD图
  绘制 log-FC 与 log-cpm 均值 - 差异图。
  图中的蓝色横线代表2倍差异 (2-fold changes)。

plotMD(lrt)
abline(h = c(-1,1),col = "blue")