limma包差异分析

差异分析的原理就是判断组间（处理）差异是否显著大于组内（误差）差异。比如两组样本，高浓度镉离子处理的拟南芥，和正常对照（没有处理）的拟南芥，每组三个植株。由于不同个体间会有个体差异，那么如何判断基因表达差异是由个体间差异导致的还是处理导致的？这就涉及到组间差异和组内差异的比较了。一般差异分析软件如edgeR、DESeq的原理就是比较组间差异和组内差异。所以说通常情况下，如果没有重复，则组内差异无法计算，这就是为什么我们一直强调要有生物学重复样本的原因。但是软件在没有重复样本的时候也可以计算组内差异，是通过公式预估或者人为预设误差差异的。所以说现在不管有没有重复样本，都是可以进行差异分析的，只是有重复样本的话，才能准确估算个体差异带来的误差，差异分析结果才会更准确。
基于R语言进行差异分析的包有很多个，本次详细介绍使用limma（Linear Models for Microarray Data）进行差异分析。

什么是limma？

原理： limma 最初是针对基因芯片数据开发的，但后来也被应用于 RNA 测序数据。limma 基于线性模型，使用加权最小二乘法来估计基因表达的差异，并通过贝叶斯方法来校正多重检验问题。limma 通常假设基因表达数据服从正态分布，在处理大规模数据时表现出色，适用于高通量数据分析，如芯片和大规模RNA测序数据，能够很好地控制假阳性率。

limma也是广义线性模型的一种，其对每个gene的表达量拟合一个线性方程。limma的分析过程包括ANOVA分析、线性回归等。
limma对每个gene拟合出这样一个方程，其中：
X可以是：

• 一个连续变量：如pH，RIN值，年龄，体重，身高…
• 一个分类变量：如性别、种族、与中位数比较的gene高低表达…

$\beta$ 是limma将要求出的值
$ϵ$ 是假定在整个数据集中正态分布的残差(residual)

1.数据准备

差异分析是两两数据集间的比较，一个是对照组样本(CK)，一个是处理组样本(Treat)。

2.数据分析

2.1 导入所需的包

## 导入R包
library(limma)
library(dplyr)

2.2 加载数据

df <- read.table("Counts_data.txt", header = T, sep = "\t", row.names = 1, check.names = F)
head(df)

2.3 样本信息注释

list <- c(rep("Treat", 66), rep("CK",148)) %>% factor(., levels = c("CK", "Treat"), ordered = F)
##--------------
> head(list)
  CK Treat
1  0     1
2  0     1
3  0     1
4  0     1
5  0     1
6  0     1
.............
list <- model.matrix(~factor(list)+0)  #把group设置成一个model matrix
colnames(list) <- c("CK", "Treat")
df.fit <- lmFit(df, list)  ## 数据与list进行匹配

2.4 差异分析

df.matrix <- makeContrasts(Treat - CK , levels = list)
fit <- contrasts.fit(df.fit, df.matrix)
fit <- eBayes(fit)
tempOutput <- topTable(fit,n = Inf, adjust = "fdr")
> head(tempOutput)
               logFC  AveExpr         t      P.Value    adj.P.Val        B
gene11796 -1.2620803 5.551402 -8.392278 5.886136e-15 7.083964e-11 23.25729
gene7364  -0.9825962 7.471963 -7.598215 8.629186e-13 5.192613e-09 18.51600
gene11454 -1.3204341 6.238318 -7.368543 3.472935e-12 1.301549e-08 17.19354
gene11689 -1.0769861 4.967290 -7.331950 4.325880e-12 1.301549e-08 16.98503
gene11227 -1.2035217 5.925234 -7.263016 6.531456e-12 1.572122e-08 16.59390
gene3259  -2.3695741 9.467290 -6.962632 3.829470e-11 5.760959e-08 14.91576

到这部分，差异分析就结束了。后面的部分就是提取差异结果，和绘制火山图。

参考：
https://zhuanlan.zhihu.com/p/437712423
https://www.jianshu.com/p/f009bea514af
https://zhuanlan.zhihu.com/p/388984969
https://zhuanlan.zhihu.com/p/653841949