limma：RNA-Seq Data

生成计数矩阵

标准化和过滤

1. 第一步先创建DGEList对象。
   dge = DGEList(counts = counts)
2. 移除0和低计数值得行。

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

3. 标准化。
   dge = calcNormFactors(dge)

差异表达：limma-trend

适用： 测序深度在RNA样本是一致的，最大文库和最小文库的大小比率不超过3倍。

1. 将计数转换为cpm值。使用prior.count参数降低计数对数的方差。
   logCPM = cpm(dge,log = TRUE,prior.count = 3)
2. 运行eBayes或treat时使用trend=TRUE参数。

fit = lmFit(logCPM,design)
fit = eBayes(fit,trend = TRUE)
topTable(fit,coef = ncol(design))

3. 如果想让gene ranking的fold changes更明显。

fit = lmFit(logCPM,design)
fit = treat(fit,lfc = log2(1.2))
topTreat(fit,coef = ncol(design))

差异表达：limma-voom

适用： 样品间文库变化很大的时候。

1. voom()应用于标准化和过滤后的DGEList
   v = voom(dge,design,plot = TRUE)
   也可以直接对计数矩阵进行voom转换。
   v = voom(counts,design,plot = TRUE)
   如果数据噪音很多：
   v = voom(counts,design,plot = TRUE,normalize = "quantile")

2. 常规limma pipeline

fit = lmFit(v,design)
fit = eBayes(fit)
topTable(fit,coef = ncol(design))

3. 样本质量加权
   plotMDS

差异剪切

limma检测数基因在不同条件下剪切的区别。

1. 这个情况下，计数矩阵必须是外显子单位计数，每行一个外显子。
   DGEList的genes矩阵中的GeneID标示每个外显子属于哪个基因。

dge = DGEList(counts = counts)
dge$genes$GeneID = GeneID

2. 过滤，标准化。

A  = rowSums(dge$counts)
dge = dge[A>10,,keep.lib.sizes = FALSE]
dge = calcNormFactors(dge)

3. voom转换，拟合线性模型。

v = voom(dge,design,plot = TRUE)
fit = lmFit(v,design)

4. 测试不同剪切和线性模型系数之间的联系。运行diffSplice（）

ex = diffSplice(fit,geneid = "GeneID")
topSplice(ex, coef=2, test="simes") #指定与线性模型第二个系数的关系

5. 展示单个外显子相对于同一个基因中其他外显子富集情况。
   plotSplice(ex)