用limma做差异基因/蛋白分析

limma（Linear Models for Microarray and RNA-Seq Data）早期主要用于生物芯片数据的输入、前处理（归一化）、批次效应、背景去除、等数据处理功能和对差异化基因分析的算法，limma设计的线性回归和经验贝叶斯方法找差异基因非常值得学习，被迁移到很多生信分析场景，如：NGS差异数据分析、蛋白数据分析等等。

limma包的文档也非常的全面，总结了各种条件下的数据方法，是生物信息领域中最经典的数据分析包。

核心函数：

# 线性模型拟合给定数组的每个基因
lmFit(object, design=NULL, ndups=1, spacing=1, block=NULL, correlation, weights=NULL,method="ls", ...)

# 计算给定分组的标准误差和估计系数
contrasts.fit(fit, contrasts=NULL, coefficients=NULL)

# 计算T检验、F检验，并通过标准误差的经验Bayes来计算差异表达的概率
eBayes(fit, proportion = 0.01, stdev.coef.lim = c(0.1,4),
       trend = FALSE, robust = FALSE, winsor.tail.p = c(0.05,0.1))

# 提取差异信息
topTable(fit, coef=NULL, number=10, genelist=fit$genes, adjust.method="BH", sort.by="B", resort.by=NULL, p.value=1, lfc=0, confint=FALSE)

经典分析套路：

# 定义模型、标签，实验设计矩阵
design <- model.matrix(~0+group)
colnames(design) <- levels(group)
# 建立比较分组
contrast <- makeContrasts(contrasts=paste0(unique(group), collapse = "-"), levels = design)
# 线性回归方法
fit <- lmFit(logCPM, design)
fit <- contrasts.fit(fit, contrast) # 可选
# Bayes估计
fit <- eBayes(fit, trend=TRUE)
# 过滤差异基因信息
# adjust.method包含"none", "BH", "BY" 和 "holm"方法
topTable(fit, adjust.method = "BH", coef=ncol(design), number = Inf)

提取差异基因用topTable即可，具体可以利用R帮助文档调整过滤条件。

limma的文档中提到很多案例，建议手动试试。

参考资料：

1.http://www.bioconductor.org/packages/release/bioc/vignettes/limma/inst/doc/usersguide.pdf