limma包是2015年发表在Nucleic Acids Resarch一个做差异分析的工具，目前引用次数高达七千多次。最流行的差异分析软件就是limma了，它现在更新了一个voom的算法，所以既可以对芯片数据，也可以对转录组高通量测序数据进行分析，其它所有的差异分析软件其实都是模仿这个的。

做差异分析，需要三个数据：

• 表达矩阵
• 分组矩阵
• 差异比较矩阵

前面两个肯定是必须的，有表达矩阵，样本必须进行分组，才能分析，但是我看到过好几种例子，有的有差异比较矩阵，有的没有。

后来我仔细研究了一下limma包的说明书，发现这其实是一个很简单的问题。

大家仔细观察下面的两个代码

首先是不需要差异比较矩阵的

<pre style="border: 0px; font: 13px / 1.5 "Courier 10 Pitch", Courier, monospace; margin: 0px 0px 1.625em; outline: 0px; padding: 0.75em 1.625em; vertical-align: baseline; background: rgb(244, 244, 244); overflow: auto;"> library(CLL)
data(sCLLex)
library(limma)
design=model.matrix(~factor(sCLLex$Disease))
fit=lmFit(sCLLex,design)
fit=eBayes(fit)
options(digits = 4)
#topTable(fit,coef=2,adjust='BH')
> topTable(fit,coef=2,adjust='BH')
logFC AveExpr t P.Value adj.P.Val B
39400_at 1.0285 5.621 5.836 8.341e-06 0.03344 3.234
36131_at -0.9888 9.954 -5.772 9.668e-06 0.03344 3.117
33791_at -1.8302 6.951 -5.736 1.049e-05 0.03344 3.052
1303_at 1.3836 4.463 5.732 1.060e-05 0.03344 3.044
36122_at -0.7801 7.260 -5.141 4.206e-05 0.10619 1.935
36939_at -2.5472 6.915 -5.038 5.362e-05 0.11283 1.737
41398_at 0.5187 7.602 4.879 7.824e-05 0.11520 1.428
32599_at 0.8544 5.746 4.859 8.207e-05 0.11520 1.389
36129_at 0.9161 8.209 4.859 8.212e-05 0.11520 1.389
37636_at -1.6868 5.697 -4.804 9.355e-05 0.11811 1.282
</pre>

然后是需要差异比较矩阵的

<pre style="border: 0px; font: 13px / 1.5 "Courier 10 Pitch", Courier, monospace; margin: 0px 0px 1.625em; outline: 0px; padding: 0.75em 1.625em; vertical-align: baseline; background: rgb(244, 244, 244); overflow: auto;"> library(CLL)
data(sCLLex)
library(limma)
design=model.matrix(~0+factor(sCLLex$Disease))
colnames(design)=c('progres','stable')
fit=lmFit(sCLLex,design)
cont.matrix=makeContrasts('progres-stable',levels = design)
fit2=contrasts.fit(fit,cont.matrix)
fit2=eBayes(fit2)
options(digits = 4)
topTable(fit2,adjust='BH')

           logFC AveExpr      t   P.Value adj.P.Val     B
39400_at -1.0285   5.621 -5.836 8.341e-06   0.03344 3.234
36131_at  0.9888   9.954  5.772 9.668e-06   0.03344 3.117
33791_at  1.8302   6.951  5.736 1.049e-05   0.03344 3.052
1303_at  -1.3836   4.463 -5.732 1.060e-05   0.03344 3.044
36122_at  0.7801   7.260  5.141 4.206e-05   0.10619 1.935
36939_at  2.5472   6.915  5.038 5.362e-05   0.11283 1.737
41398_at -0.5187   7.602 -4.879 7.824e-05   0.11520 1.428
32599_at -0.8544   5.746 -4.859 8.207e-05   0.11520 1.389
36129_at -0.9161   8.209 -4.859 8.212e-05   0.11520 1.389
37636_at  1.6868   5.697  4.804 9.355e-05   0.11811 1.282</pre>

大家运行一下这些代码就知道，两者结果是一模一样的。

而差异比较矩阵的需要与否，主要看分组矩阵如何制作的！

design=model.matrix(~factor(sCLLex$Disease))

design=model.matrix(~0+factor(sCLLex$Disease))

有本质的区别！！！

前面那种方法已经把需要比较的组做出到了一列，需要比较多次，就有多少列，第一列是截距不需要考虑，第二列开始往后用coef这个参数可以把差异分析结果一个个提取出来。

而后面那种方法，仅仅是分组而已，组之间需要如何比较，需要自己再制作差异比较矩阵，通过makeContrasts函数来控制如何比较！

实战实战实战实战

基因表达差异分析是我们做转录组最关键根本的一步，edgeR+limma是目前最为推荐的方式。本文结合示例数据，将对这个过程进行梳理，让你明白limma包的why，what，how。
本文示例数据下载

什么是limma？

首先要明白，不管哪种差异分析，其本质都是广义线性模型。limma也是广义线性模型的一种，其对每个gene的表达量拟合一个线性方程。limma的分析过程包括ANOVA分析、线性回归等。

image.png

数据解释

本文数据有两个因素，均为分类变量

• 品种：cultivar（C,I5/I8）
• 时间：time(6,9)
  cols为样本编号，rows为基因表达。和我们平时用的数据一致。

开始分析

1. 准备数据

library(edgeR) #edgeR将同时引入limma
counts <- read.delim("all_counts.txt", row.names = 1)
head(counts)

d0 <- DGEList(counts)
# 注意： calcNormFactors并不会标准化数据，只是计算标准化因子
d0 <- calcNormFactors(d0)
d0

# 过滤低表达
cutoff <- 1
drop <- which(apply(cpm(d0), 1, max) < cutoff)
d <- d0[-drop,] 
dim(d) # number of genes left

# sample names
snames <- colnames(counts)  
snames
# 此数据有两个因素：cultivar（C,I5/I8）和time(6,9)
cultivar <- substr(snames, 1, nchar(snames) - 2) 
time <- substr(snames, nchar(snames) - 1, nchar(snames) - 1)
cultivar
time

# Create a new variable “group” that combines cultivar and time
group <- interaction(cultivar, time)
group

# Multidimensional scaling (MDS) plot
plotMDS(d, col = as.numeric(group))

我们首先构建了edgeR的DGEList对象，这个对象将来将会转化成limma中的EList对象。然后计算了标准化因子，过滤低表达基因。然后按照分组整理了列名，并进行初步的MDS plot，以便看到样品的大概分布

image

2. limma

mm <- model.matrix(~0 + group)
par(mfrow = c(1, 3))
y <- voom(d, mm, plot = T)
#voom的曲线应该很光滑，比较一下过滤低表达gene之前的图形：
voom(d0, mm, plot = T)
# lmFit fits a linear model using weighted least squares for each gene:
fit <- lmFit(y, mm)
head(coef(fit))

#Comparisons between groups (log fold-changes) are obtained as contrasts of these fitted linear models:
# Comparison between times 6 and 9 for cultivar I5
# makeContrasts实际就是定义比较分组信息
contr <- makeContrasts(groupI5.9 - groupI5.6, levels = colnames(coef(fit)))
# 比较每个基因
tmp <- contrasts.fit(fit, contr)
# Empirical Bayes smoothing of standard errors ， (shrinks standard errors that are much larger or smaller than those from other genes towards the average standard error) 
# (see https://www.degruyter.com/doi/10.2202/1544-6115.1027)
tmp <- eBayes(tmp)
# 使用plotSA 绘制了log2残差标准差与log-CPM均值的关系。平均log2残差标准差由水平蓝线标出
plotSA(tmp, main="Final model: Mean-variance trend")
# topTable 列出差异显著基因
top.table <- topTable(tmp, sort.by = "P", n = Inf)
# logFC: log2 fold change of I5.9/I5.6
# AveExpr: Average expression across all samples, in log2 CPM
# t: logFC divided by its standard error
# P.Value: Raw p-value (based on t) from test that logFC differs from 0
# adj.P.Val: Benjamini-Hochberg false discovery rate adjusted p-value
# B: log-odds that gene is DE (arguably less useful than the other columns)
head(top.table, 20)

# p值<0.05的基因有多少个？
length(which(top.table$adj.P.Val < 0.05))

#Write top.table to a file
top.table$Gene <- rownames(top.table)
top.table <- top.table[,c("Gene", names(top.table)[1:6])]
write.table(top.table, file = "time9_v_time6_I5.txt", row.names = F, sep = "\t", quote = F)

limma的核心步骤包括voom、fit、eBays等步骤，注释里都有详细说明。最后我们用topTable方法按照p值排序输出结果。
下图是我为了说明绘制的，顺序反了。。。中间的是没有过滤低表达基因之前的，左边是过滤后的，最后是fit后的，可以明显的看出区别。

image

这时差异分析就有已经完成了，怎样，是不是很简单？

使用limma进行双变量、多变量、连续变量分析

###################双变量分析（cultivar+time）########################
mm <- model.matrix(~cultivar*time)
y <- voom(d, mm, plot = F)
fit <- lmFit(y, mm)
head(coef(fit))
# (Intercept)  cultivarI5 cultivarI8      time9 cultivarI5:time9 cultivarI8:time9
# AT1G01010    4.837410  0.53644370  0.2279446 0.20580445      -0.05565729       0.09265044
# AT1G01020    3.530869 -0.03152318 -0.3180096 0.15875297       0.06289715       0.36468449
# AT1G01030    1.250817 -0.32143420  0.3084243 0.03477863      -0.48099113      -0.37842909
# AT1G01040    5.676015  0.27097286  0.1028739 0.50635951      -0.58923660      -0.46975071
# AT1G01050    6.598712 -0.09734846 -0.1347759 0.02052702       0.23139851       0.22730960
# AT1G01060    7.807988 -0.34550979 -0.4172467 1.15805850      -0.34989810      -0.17267051
# 这个表中显示的是coefficient（相关系数） 
# cultivarI5 这一列表示cultivar I5 组均值 vs cultivar C（参考cultivar）的差异, for time 6 (the reference level for time)
#  time9 这一列表示time9 组均值 vs time6 ,forcultivar C的差异
# cultivarI5:time9 : the difference between times 9 and 6 of the differences between cultivars I5 and C (interaction effect)

# 接下来我们可以定义fit中的coef参数，来进行组间fit
# Let’s estimate the difference between cultivars I5 and C at time 6
tmp <- contrasts.fit(fit, coef = 2) #  the difference in mean expression between cultivar I5 and the reference cultivar (cultivar C), for time 6 (the reference level for time)
tmp <- eBayes(tmp)
top.table <- topTable(tmp, sort.by = "P", n = Inf)
head(top.table, 20)

tmp <- contrasts.fit(fit, coef = 5) # Test cultivarI5:time9
tmp <- eBayes(tmp)
top.table <- topTable(tmp, sort.by = "P", n = Inf)
head(top.table, 20)

####################多变量分析########################
#让事情更复杂一点，我们加入批次信息
batch <- factor(rep(rep(1:2, each = 2), 6))
# 只需要重新定义model matrix，其余都一样
mm <- model.matrix(~0 + group + batch)
y <- voom(d, mm, plot = F)
fit <- lmFit(y, mm)
contr <- makeContrasts(groupI5.6 - groupC.6, levels = colnames(coef(fit)))
tmp <- contrasts.fitit(fit, contr)
tmp <- eBayes(tmp)
top.table <- topTable(tmp, sort.by = "P", n = Inf)
head(top.table, 20)

# 加入连续变量
# Generate example RIN data
set.seed(99)
RIN <- rnorm(n = 24, mean = 7.5, sd = 1)
RIN
mm <- model.matrix(~0 + group + RIN)
y <- voom(d, mm, plot = F)
fit <- lmFit(y, mm)
contr <- makeContrasts(groupI5.6 - groupC.6, levels = colnames(coef(fit)))
tmp <- contrasts.fit(fit, contr)
tmp <- eBayes(tmp)
top.table <- topTable(tmp, sort.by = "P", n = Inf)
head(top.table, 20)

# What if we want to look at the correlation of gene expression with a continuous variable like pH?
# Generate example pH data
set.seed(99)
pH <- rnorm(n = 24, mean = 8, sd = 1.5)
pH
mm <- model.matrix(~pH)
head(mm)
y <- voom(d, mm, plot = F)
fit <- lmFit(y, mm)
tmp <- contrasts.fit(fit, coef = 2) # test "pH" coefficient
tmp <- eBayes(tmp)
top.table <- topTable(tmp, sort.by = "P", n = Inf)
head(top.table, 20)

上面，我们分别加入了额外的二分变量、连续变量进行limma分析，结果都很好。
这就是有关limma分析的全部内容，注释写的很清楚，可以用这个流程分析任何转录组数据，进行差异表达分析。

参考：
http://www.bio-info-trainee.com/tag/limma
https://www.jianshu.com/p/f009bea514af?utm_campaign=haruki