052-limma前瞻

刘小泽写于18.11.15

limma是做差异分析的好工具，真的需要好好阅读原版工具书 我需要静下心来好好看看limma内涵了，不想就这么学一个流程

什么是limma？

limma是芯片、转录组分析的金标准,从2005年被开发出来，已经被引用了快5000次，之所以厉害是因为它包含了数据导入、预处理、质控、标准化、构建线性模型、差异基因分析等一整套的流程

它需要的格式是行为基因，列为样本，中间数据为原始表达量（raw count）的表达矩阵。我们最关注的就是差异分析操作，核心就是利用线性模型去估计各个分组的基因表达量的均值与方差。主要利用两个函数，lmFit （用于线性拟合）与 eBayes（根据拟合结果进行推断）。

流程

首先需要三大矩阵：

• 表达矩阵：就是熟知的matrix形式
• 分组矩阵：也就是实验设计（design）矩阵
• 差异比较矩阵：就是说那几个组之间进行比较

然后就是线性模型拟合、贝叶斯检验、得到结果

因此我们看到最关键的部分就是如何把我们的数据构建成软件要求的样子

构建分组矩阵

design <- model.matrix(~0+group+lane)
colnames(design) <- gsub("group", "", colnames(design)) 
design

构建分组矩阵一般有两种：

• 一种是这里演示的~0+group+lane，它不需要考虑第一列group的截距，第二列lane保留了截距，那么接下来怎么比较，需要自己再用makeContrscts定义差异比较矩阵；

  关于截距intercept，我的理解是这样的：在有截距的模型中，截距是计算某种表型的平均值，所有其他的表型系数都是该种表型与截距之间的差值。因此，这些表达值的差异由大到小，取决于两种表型的差异

• 第二种是~group+lane，在group和lane中都保留了截距，因此会自己比较多次，再将差异结果一一提取

这里选择第一种构建方式，原因是不考虑group的截距时，设置的比较模型更加直观

构建差异比较矩阵

contr.matrix <- makeContrasts(
   BasalvsLP = Basal-LP,
   BasalvsML = Basal - ML,
   LPvsML = LP - ML,
   levels = colnames(design))

得到

contr.matrix
## Contrasts
## Levels BasalvsLP BasalvsML LPvsML 
##	Basal 	1			1		0 
##	LP 		-1			0		1 
##	ML 		0			-1		-1 
## laneL006 0 			0 		0 
## laneL008 0 			0 		0

limma包的优势在这里体现出来了，它可以构建简单的比较矩阵（这里的celltype和lane），还能hold住更多组比较

正式的limma流程

# 第一步
vfit <- lmFit(v, design)
# 第二步
vfit <- contrasts.fit(vfit, contrasts=contr.matrix)
efit <- eBayes(vfit)
# 第三步
ofit <- topTable(efit, coef=1, n=Inf)
DEG <- na.omit(ofit)
head(DEG)