2.3.Differential Expression with DEseq2 and edgeR

本节我们将使用DESeq2和edgeR 完成差异表达（Differential Expression）分析。

注意：本章需要读者具有R的编程基础。

1) Getting software & data

本章不需要 docker，可以直接在自己的计算机上安装以下的R 包，并下载相应文件。

软件依赖

• 在基因差异表达分析中，DESeq2和edgeR是两个最常用的R package。

• 这两个package使用的统计模型非常相似，都是用基于负二项分布的广义线性模型直接对counts类型的数据建模，来计算统计显著性。

• 这两个package都host在bioconductor上，可以通过BiocManager::install来安装:

#如果未安装BiocManager，可使用install.packages("BiocManager")安装
BiocManager::install("DESeq2")
BiocManager::install("edgeR")

数据

count_exon.txt(可以从 **** Files needed 中的Files/ 路径下的相应文件夹中下载)是拟南芥野生型(WT)与uvr8突变型(uvr8)在光照前后的基因表达矩阵

• 第1列为基因名称

• 第2-4列为WT光照前的表达值

• 第5-7列为WT光照后的表达值

• 第8-10列为_uvr8_光照前的表达值

• 第11-13列为_uvr8_光照后的表达值

• 在本例中，我们考虑野生型在光照前后基因表达量的变化，所以只用到前7列的数据。

2) DESeq2

2a) 准备输入数据

• 读取野生型拟南芥的基因表达矩阵:

#读取数据
raw.counts <- read.table("count_exon.txt", sep='\t', header = T,row.names = 1)
#我们这里只使用野生型数据进行分析
wt.raw.counts <- raw.counts[,c("CD1_1", "CD1_2", "CD1_3", "CD0_1", "CD0_2", "CD0_3")]
#过滤掉表达量过低的基因
wt.filtered.counts <- wt.raw.counts[rowMeans(wt.raw.counts) > 5, ]

Tips: 由于deseq2和edgeR都是用负二项分布直接对counts建模，我们输入矩阵需为原始的count矩阵，而不应该是CPM/FPKM/TPM一类的数值。

2b) 提供样本条件信息

这里我们用一个叫做colData的数据框来存储样本信息。这里每行对应一个样本。我们只考虑一种条件，所以数据框只有一列。

# "CD1_1", "CD1_2", "CD1_3" 三个样本为control
# "CD0_1", "CD0_2", "CD0_3 三个样本对应treatment
conditions <- factor(c(rep("Control", 3), rep("Treatment", 3)),levels = c("Control","Treatment"))
colData <- data.frame(row.names = colnames(wt.filtered.counts),conditions=conditions)

2c) 差异分析

library(DESeq2)
# 我们到这里才开始使用DESeq2 package
# library(DESeq2)会输出一长串的提示信息，如果不需要可使用suppressPackageStartupMessages(library(DESeq2))

dds <- DESeqDataSetFromMatrix(wt.filtered.counts, colData, design = ~conditions)
#进行差异分析
dds <- DESeq(dds)
#获取结果
res <- results(dds)

• 以上我们就获得了差异分析的结果。可以看出，在我们这个例子中，只有dds2 <- DESeq(dds)这一行代码真正在计算两组之间的差异，以及差异的显著性，其他代码都是在准备输入输出。

• 其实dds2 <- DESeq(dds)的内部实现是比较复杂的，它实际上顺序的调用了DESeq2 package中的三个函数:

  • estimateSizeFactors: 对library size的大小进行估计，作为模型的一个系数

  • estimateDispersons: 对负二项分布模型进行参数估计

  • nbinomWaldTest: 用Wald test来检验两个条件之间差异的显著性，即广义线性模型中对应的系数不为0的显著性

2d) 保存结果

# 如果你想保存所有结果，也可以不过滤
write.table(res,"wt.light.vs.dark.all.txt", sep='\t', row.names = T, quote = F)
# 你也可以筛选出有差异的基因
# 过滤标准: padj < 0.05, log2 fold change > 1
diff.table <- subset(res, padj < 0.05 & abs(log2FoldChange) > 1)
# 将结果保存至"wt.light.vs.dark.txt
write.table(diff.table,"wt.light.vs.dark.txt", sep='\t', row.names = T, quote = F)

• 输出的结果是一个表格，每个基因对应一行，每列含义如下:

  • baseMean: 平均表达量

  • log2FoldChange: treatment相对于control的log2 fold change

  • lfcSE: 估计出的log2FoldChange的标准差

  • stat: 假设检验用到的统计量

  • pvalue: p值

  • padj: 经过multiple test correction的p值

• 我们一般会根据log2FoldChange和padj两列信息对基因进行筛选，并进行下游的通路富集等分析

3) edgeR

3a) 准备输入数据

• 读取表达矩阵。这一步和使用deseq2没有区别，如果你之前已经读取了表达矩阵，就不用做这步操作了。

#读取数据
raw.counts <- read.table("count_exon.txt", sep='\t', header = T,row.names = 1)
#我们这里只使用野生型数据进行分析
wt.raw.counts <- raw.counts[,c("CD1_1", "CD1_2", "CD1_3", "CD0_1", "CD0_2", "CD0_3")]
#过滤掉表达量过低的基因
wt.filtered.counts <- wt.raw.counts[rowMeans(wt.raw.counts) > 5, ]

3b) 提供样本条件信息

conditions <- factor(c(rep("Control", 3), rep("Treatment", 3)),levels = c("Control","Treatment"))
#获取design矩阵
design <- model.matrix(~conditions)

3c) 差异分析

edgeR对差异分析提供了三种检验方法:

• exact test

• likelihood ratio test

• quasi-likelihood glm

其中exact test针对两组样本的情形，另外两种检验适用于更复杂的实验设计。我们这里likelihood ratio test为例进行介绍。

library(edgeR) # 至此我们开始用到edgeR package
y <- DGEList(counts = wt.filtered.counts) # 定义edgeR用于存储基因表达信息的DGEList对象

# TMM标准化 (TMM 实际上是edgeR的默认参数)
y <- calcNormFactors(y, method="TMM")

# 估计dispersion
y <- estimateDisp(y,design = design)
# edgeR内部进行了以下三步调用，有兴趣的同学可以查阅文档，看一看它们分别在做什么事情
# y <- estimateGLMCommonDisp(y,design = design)
# y <- estimateGLMTrendedDisp(y,design = design)
# y <- estimateGLMTagwiseDisp(y,design = design)

# 拟合广义线性模型
fit <- glmFit(y, design = design)

# 似然比检验
# coef = 2指的是对design矩阵的第二列（即是否照光）对应的系数进行检验
lrt <- glmLRT(fit,coef=2) 

• 还有一种做法是用design <- model.matrix(~0+conditions)定义design matrix，再用lrt <- glmLRT(fit,contrast=c(-1,1))进行检验，结果应当是相同的，大家如果希望进一步了解请自行参考edgeR的文档。

3d) 选取差异显著的基因

# 这里tag就是基因的意思，topTags意思是返回变化最top的基因
# 默认返回10个基因，按p值排序
# 这里我们用n = nrow(y)要求它返回所有基因的结果
diff.table <- topTags(lrt, n = nrow(y))$table

# 保存差异分析结果
write.table(diff.table, file = 'edger.wt.light.vs.dark.txt', sep = "\t", quote = F, row.names = T, col.names = T)

# 当然你也可以只挑选显著变化的基因
diff.table.filtered <- diff.table[abs(diff.table$logFC) > 1 & diff.table$FDR < 0.05,]

topTag返回的结果中包括的差异表格每行对应一个基因，由以下几列组成：

• logFC: 相对control的log2 fold change

• logCPM: 表达量

• LR: 似然比统计量

• PValue: p值

• FDR: 经过multiple test correction的p值

4) Homework

1. 什么是Multiple test correction? 并解释 q value(很多时候也叫FDR) 和 p value 的差别。

2. 请结合上课时所讲的知识阐述DESeq2和edgeR中如何对数据进行 normalization，列出并解释具体的公式 。

3. 利用我们以上介绍的方法和数据，分别使用DESeq2和edgeR找出uvr8突变型（uvr8）在光照前后的差异基因，保存为文本文件

4. 对于uvr8突变型的差异基因，定义|log2FC|>1，FDR<0.05的基因为差异表达基因。比较两个软件得到的差异基因有多少是重合的，有多少是不同的，用venn图的形式展示

5. 对于edgeR找出的FDR<0.05的基因，选出log2FoldChange最大的10个基因和最小的10个基因。计算原始表达量矩阵的log10CPM值并对每个基因进行Z-score处理，使用刚才筛选出来的20个基因绘制热图（heatmap）作为最后结果输出。

5) References

我们介绍的两个R package都提供了极为详尽的文档,推荐希望进一步了解的同学阅读

• deseq2: https://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html
• edgeR: https://www.bioconductor.org/packages/release/bioc/vignettes/edgeR/inst/doc/edgeRUsersGuide.pdf

Here is a github page, 用edgeR寻找差异基因, shared by students of our course.