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")