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