# 2.3.Differential Expression with DEseq2 and edgeR 本节我们将使用[**DESeq2**](https://bioconductor.org/packages/release/bioc/html/DESeq2.html)和[**edgeR**](https://bioconductor.org/packages/release/bioc/html/edgeR.html) 完成差异表达（Differential Expression）分析。 > 注意：本章需要读者具有R的编程基础。 ## 1) Getting software & data 本章不需要 docker，可以直接在自己的计算机上安装以下的R 包，并下载相应文件。 ### 软件依赖 * 在基因差异表达分析中，[DESeq2](https://bioconductor.org/packages/release/bioc/html/DESeq2.html)和[edgeR](https://bioconductor.org/packages/release/bioc/html/edgeR.html)是两个最常用的R package。 * 这两个package使用的统计模型非常相似，都是用基于负二项分布的广义线性模型直接对counts类型的数据建模，来计算统计显著性。 * 这两个package都host在bioconductor上，可以通过`BiocManager::install`来安装: ```bash #如果未安装BiocManager，可使用install.packages("BiocManager")安装 BiocManager::install("DESeq2") BiocManager::install("edgeR") ``` ### 数据 `count_exon.txt`(可以从 \*\*\*\* [**Files needed** ](https://courses.ncrnalab.org/files)中的**Files/** 路径下的相应文件夹中下载)是拟南芥野生型(WT)与uvr8突变型(*uvr8*)在光照前后的基因表达矩阵 * 第1列为基因名称 * 第2-4列为WT光照前的表达值 * 第5-7列为WT光照后的表达值 * 第8-10列为\_uvr8\_光照前的表达值 * 第11-13列为\_uvr8\_光照后的表达值 * 在本例中，我们考虑野生型在光照前后基因表达量的变化，所以只用到前7列的数据。 ## 2) DESeq2 ### 2a) 准备输入数据 * 读取野生型拟南芥的基因表达矩阵: ```r #读取数据 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`的数据框来存储样本信息。这里每行对应一个样本。我们只考虑一种条件，所以数据框只有一列。 ```r # "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) 差异分析 ```r 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) 保存结果 ```r # 如果你想保存所有结果，也可以不过滤 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没有区别，如果你之前已经读取了表达矩阵，就不用做这步操作了。 ```r #读取数据 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) 提供样本条件信息 ```r 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为例进行介绍。 ```r 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) 选取差异显著的基因 ```r # 这里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: * edgeR: > Here is a github page, [用edgeR寻找差异基因](https://yj-mo.github.io/2020/10/26/Team/#%E7%94%A8edger%E5%AF%BB%E6%89%BE%E5%B7%AE%E5%BC%82%E5%9F%BA%E5%9B%A0), shared by students of our course.