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

什么是limma？

首先要明白，不管哪种差异分析，其本质都是广义线性模型。limma也是广义线性模型的一种，其对每个gene的表达量拟合一个线性方程。limma的分析过程包括ANOVA分析、线性回归等。
$Y=β0+β1X1+β2X2+⋯+βpXp+ϵ$
limma对每个gene拟合出这样一个方程，其中：
$X$ 可以是：

一个连续变量：如pH，RIN值，年龄，体重，身高...
一个分类变量：如性别、种族、与中位数比较的gene高低表达...
$β$ 是limma将要求出的值
$ϵ$ 是假定在整个数据集中正态分布的残差(residual)

数据解释

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

品种：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.png

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.png

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

使用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分析的全部内容，注释写的很清楚，可以用这个流程分析任何转录组数据，进行差异表达分析。

limma 差异分析透彻讲解