2.3.Differential Expression with DEseq2 and edgeR

本节我们将使用DESeq2edgeR 完成差异表达(Differential Expression)分析。

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

1) Getting software & data

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

软件依赖

  • 在基因差异表达分析中,DESeq2edgeR是两个最常用的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值

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

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都提供了极为详尽的文档,推荐希望进一步了解的同学阅读

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

Last updated