Comment on page
2.3.Differential Expression with DEseq2 and edgeR
注意：本章需要读者具有R的编程基础。
本章不需要 docker，可以直接在自己的计算机上安装以下的R 包，并下载相应文件。
#如果未安装BiocManager，可使用install.packages("BiocManager")安装
BiocManager::install("DESeq2")
BiocManager::install("edgeR")
- 第1列为基因名称
- 第2-4列为WT光照前的表达值
- 第5-7列为WT光照后的表达值
- 第8-10列为_uvr8_光照前的表达值
- 第11-13列为_uvr8_光照后的表达值
- 在本例中，我们考虑野生型在光照前后基因表达量的变化，所以只用到前7列的数据。
- 读取野生型拟南芥的基因表达矩阵:
#读取数据
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一类的数值。
这里我们用一个叫做
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)
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的显著性
# 如果你想保存所有结果，也可以不过滤
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两列信息对基因进行筛选，并进行下游的通路富集等分析
- 读取表达矩阵。这一步和使用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, ]
conditions <- factor(c(rep("Control", 3), rep("Treatment", 3)),levels = c("Control","Treatment"))
#获取design矩阵
design <- model.matrix(~conditions)
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的文档。
# 这里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值
- 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并作热图（heatmap）
我们介绍的两个R package都提供了极为详尽的文档,推荐希望进一步了解的同学阅读