# 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. --- # Agent Instructions: Querying This Documentation If you need additional information that is not directly available in this page, you can query the documentation dynamically by asking a question. Perform an HTTP GET request on the current page URL with the `ask` query parameter: ``` GET https://book.ncrnalab.org/teaching/part-iii.-ngs-data-analyses/2.rna-seq/2.3.differential_expression_with_deseq2-edger.md?ask= ``` The question should be specific, self-contained, and written in natural language. The response will contain a direct answer to the question and relevant excerpts and sources from the documentation. Use this mechanism when the answer is not explicitly present in the current page, you need clarification or additional context, or you want to retrieve related documentation sections.