limma,edgeR,DESeq2三大包基本是做转录组差异分析的金标准,大多数转录组的文章都是用这三个R包进行差异分析。
edgeR差异分析速度快,得到的基因数目比较多,假阳性高(实际不差异结果差异)。
DESeq2差异分析速度慢,得到的基因数目比较少,假阴性高(实际差异结果不差异)。
需要注意的是制作分组信息的因子向量时,因子水平的前后顺序,在R的很多模型中,默认将因子向量的第一个水平看作对照组。
logFC值:
limma使用对数化的表达矩阵作为输入,所以将零假设指定在平均log值(对数几何平均数)。
edgeR和DESeq2使用原始的count矩阵作为输入,所以将原假设指定在count的平均值上(算术平均数)。
logFC一般是在1.2到2之间,筛选到差异基因的数目在500-1000左右为宜•
根据logFC的统计指标确定阈值的方法是计算logFC绝对值的平均数与2倍标准差之和。(正态分布)
DEG表示差异表达基因,differential gene express
0. 准备
判断式安装R包:如果该R包存在,可以顺带加载该R包,不需要再次library。
if(!require(stringr))install.packages("stringr") # 常规R包 if(!require(ggplotify))install.packages("ggplotify") # 热图转换为ggplot的类型 if(!require(patchwork))install.packages("patchwork") # 拼图R包 if(!require(cowplot))install.packages("cowplot") # 拼图R包 if(!require(DESeq2))BiocManager::install("DESeq2") # 差异分析R包 if(!require(edgeR))BiocManager::install("edgeR") # 差异分析R包 if(!require(limma))BiocManager::install("limma") # 差异分析R包 rm(list = ls()) load("TCGA-CHOLgdc.Rdata") table(group_list)
差异分析输入数据:表达矩阵exp,分组信息group_list。
1. DESeq2
差异分析
library(DESeq2) library(stringr) colData <- data.frame(row.names = colnames(exp), condition = group_list) # 列出每个样品是肿瘤样品还是正常样品 dds <- DESeqDataSetFromMatrix(countData = exp, colData = colData, design = ~ condition) %>% DESeq() # 将数据框转为DESeq2的数据集类型,然后用DESeq函数做差异分析
提取结果,两两比较
res <- results(dds, contrast = c("condition",rev(levels(group_list)))) # 按设置的比较水平提取dds的数据(从DESeq分析中提取结果表,给出样本的基本均值,logFC及其标准误,检验统计量,p值和矫正后的p值)。rev表示调换顺序,levels(group_list)是因子水平 DEG <- res[order(res$pvalue),] %>% as.data.frame() # 将res按照P值从小到大排序,并转换回数据框。order函数获取向量每个元素从小到大排序后的第几位 head(DEG)
添加change列标记基因上调下调
logFC_cutoff <- with(DEG,mean(abs(log2FoldChange)) + 2*sd(abs(log2FoldChange))) # 设置logFC的阈值,可以计算得出,也可以设置固定值。abs函数计算log2FoldChange的绝对值。均数+2倍标准差可以包含log2FoldChange的95%置信区间 k1 = (DEG$pvalue < 0.05)&(DEG$log2FoldChange < -logFC_cutoff) # 判断下调基因 k2 = (DEG$pvalue < 0.05)&(DEG$log2FoldChange > logFC_cutoff) # 判断上调基因 DEG$change = ifelse(k1,"DOWN",ifelse(k2,"UP","NOT")) # 标记上、下调基因 table(DEG$change) # 查看上下调基因的数量 head(DEG) DESeq2_DEG <- DEG # 保存DESeq2差异分析的结果
2. edgeR
差异分析
library(edgeR) library(stringr) # 准备count表达矩阵和分组信息,构建edgeR的DGEList对象 dge <- DGEList(counts = exp, group = group_list) # 过滤low count数据,保留在至少在两个样本中cpm值(Counts per million)大约1的基因。 # keep <- rowSums(cpm(dge) > 1 ) >= 2 # table(keep) # dge <- dge[keep, , keep.lib.sizes = FALSE] # 由于在上一步整理数据的过程中已经过滤了一遍数据了,所以这里不再过滤。 dge$samples$lib.size <- colSums(dge$counts) # 标准化基因表达分布,以TMM标准化为例。 dge <- calcNormFactors(dge) # calcNormFactors函数并不会标准化数据,只是计算标准化因子 # 假设数据符合正态分布,构建线性模型。 design <- model.matrix( ~ 0 + group_list) # 0代表x线性模型的截距为0,需要修改 rownames(design) <- colnames(dge) colnames(design) <- levels(group_list) # 计算线性模型的参数 dge <- estimateGLMCommonDisp(dge,design) %>% estimateGLMTrendedDisp(design) %>% estimateGLMTagwiseDisp(design) fit <- glmFit(dge, design) %>% # 拟合线性模型 glmLRT(contrast = c(-1,1)) # 进行差异分析。(1,-1)意味着前比后
提取结果
DEG = topTags(fit, n=nrow(exp)) %>% as.data.frame() # 提取差异分析结果 # 添加change列标记基因上调下调 logFC_cutoff <- with(DEG,mean(abs(logFC)) + 2*sd(abs(logFC))) # 设置logFC的阈值 # logFC_cutoff <- 2 k1 = (DEG$PValue < 0.05)&(DEG$logFC < -logFC_cutoff) # 判断下调基因 k2 = (DEG$PValue < 0.05)&(DEG$logFC > logFC_cutoff) # 判断上调基因 DEG$change = ifelse(k1, "DOWN", ifelse(k2, "UP", "NOT")) # 标记上、下调基因 head(DEG) table(DEG$change) edgeR_DEG <- DEG
3. limma
差异分析
library(limma) library(stringr) # 创建设计矩阵和对比。假设数据符合正态分布,构建线性模型 design <- model.matrix( ~ 0 + group_list) colnames(design) = levels(group_list) rownames(design) = colnames(exp) dge <- DGEList(counts=exp) %>% # 将count表达矩阵转换为edgeR的DGEList对象 calcNormFactors() # 进行标准化和表达量计算。calcNormFactors函数并不会标准化数据,只是计算标准化因子。 # 标准化、线性模型拟合 fit <- voom(dge, design, normalize = "quantile") %>% # voom是limma中的一种标准化方法 lmFit(design) # 线性模型拟合,出图 # 设置需要进行对比的分组 comp = paste(rev(levels(group_list)), collapse = "-") cont.matrix <- makeContrasts(contrasts = comp, levels = design) # 比较每个基因 fit2 <- contrasts.fit(fit, cont.matrix) %>% # 根据对比模型进行差值计算 eBayes() # 贝叶斯检验 # 使用plotSA绘制了log2残差标准差与log-CPM均值的关系。平均log2残差标准差由水平蓝线标出 plotSA(fit2, main = "Final model: Mean-variance trend")
提取结果
DEG = topTable(fit2, coef = comp, n=Inf) %>% # 提取差异分析结果 na.omit() # 去除差异分析结果中包含NA值的行 length(which(DEG$adj.P.Val < 0.05)) # p值<0.05的基因有多少个 # 添加change列标记基因上调下调 logFC_cutoff <- with(DEG,mean(abs(logFC)) + 2*sd(abs(logFC))) # 设置logFC的阈值 k1 = (DEG$P.Value < 0.05)&(DEG$logFC < -logFC_cutoff) # 判断下调基因 k2 = (DEG$P.Value < 0.05)&(DEG$logFC > logFC_cutoff) # 判断上调基因 DEG$change = ifelse(k1,"DOWN",ifelse(k2,"UP","NOT")) # 标记上、下调基因 table(DEG$change) head(DEG) limma_voom_DEG <- DEG
4. 数据检查
整合三大R包差异分析的数据:
tj = data.frame(deseq2 = as.integer(table(DESeq2_DEG$change)), edgeR = as.integer(table(edgeR_DEG$change)), limma_voom = as.integer(table(limma_voom_DEG$change)), row.names = c("down","not","up") );tj
验证数据:
验证logFC是否由tumor/normal计算而来:
boxplot(as.integer(exp[rownames(DESeq2_DEG)[1],]) ~ group_list) # 取exp的第一行基因表达量,以group_list为分组依据,做箱线图。修改“DESeq2_DEG”即可分别分析三种方法。 DESeq2_DEG$log2FoldChange[1] # 得到DESeq2_DEG的第一个logFC值。如果箱线图中normal大tumor小,则logFC应为负值,反之为正值。 # 如果检查发现logFC的正负值反了,则修改为负值即可:DESeq2_DEG$log2FoldChange = -DESeq2_DEG$log2FoldChange
保存数据
save(DESeq2_DEG, edgeR_DEG, limma_voom_DEG, group_list, tj, file = paste0(cancer_type,"DEG.Rdata"))
5.三大R包差异分析可视化(热图、火山图)
rm(list = ls()) setwd("D:/WorkingDirectory/01.mRNA_Expdata") load("TCGA-CHOL_01.mRNA_Expdata_count.Rdata") # 加载表达数据整理结果 setwd("D:/WorkingDirectory/03.mRNA_Diff") load("TCGA-CHOL_DEG.Rdata") # 加载差异表达基因结果 if(!require(tinyarray))devtools::install_local("tinyarray-master.zip", upgrade = F) library(ggplot2) library(tinyarray) exp[1:4,1:4] # 原始count表达矩阵数据 dat = log2(exp+1) pca.plot = draw_pca(dat, group_list); pca.plot # 作PCA图 save(pca.plot, file = paste0(cancer_type, "_PCA.Plot.Rdata")) # 保存PCA图数据以便后续拼图 # 提取差异基因 cg1 = rownames(DESeq2_DEG)[DESeq2_DEG$change != "NOT"] cg2 = rownames(edgeR_DEG)[edgeR_DEG$change != "NOT"] cg3 = rownames(limma_voom_DEG)[limma_voom_DEG$change != "NOT"] # 画热图 h1 = draw_heatmap(dat[cg1,], group_list, scale_before = T) h2 = draw_heatmap(dat[cg2,], group_list, scale_before = T) h3 = draw_heatmap(dat[cg3,], group_list, scale_before = T) # 构建“均数+2倍标准差”函数 m2d = function(x){ mean(abs(x))+2*sd(abs(x)) } # 画火山图 v1 = draw_volcano(DESeq2_DEG, pkg = 1, logFC_cutoff = m2d(DESeq2_DEG$log2FoldChange)) v2 = draw_volcano(edgeR_DEG, pkg = 2, logFC_cutoff = m2d(edgeR_DEG$logFC)) v3 = draw_volcano(limma_voom_DEG, pkg = 3, logFC_cutoff = m2d(limma_voom_DEG$logFC)) # 拼图 library(patchwork) (h1 + h2 + h3) / (v1 + v2 + v3) + plot_layout(guides = "collect") # & theme(legend.position = "none") # +表示横向,/表示纵向,guides表示将图例统一放在右侧。 # &theme(legend.position = "none")表示不显示图例。 ggsave(paste0(cancer_type, "_Heat.Vo.png"),width = 15, height = 10)
6.三大R包差异基因对比(热图、韦恩图、PCA图)
rm(list = ls()) setwd("D:/WorkingDirectory/01.mRNA_Expdata") load("TCGA-CHOL_01.mRNA_Expdata_count.Rdata") setwd("D:/WorkingDirectory/03.mRNA_Diff") load("TCGA-CHOL_DEG.Rdata") load("TCGA-CHOL_PCA.Plot.Rdata") if(!require(tinyarray))devtools::install_local("tinyarray-master.zip", upgrade = F) library(ggplot2) library(tinyarray) # 构建提取上调基因的函数: UP = function(df){ rownames(df)[df$change == "UP"] } # 构建提取下调基因的函数: DOWN = function(df){ rownames(df)[df$change == "DOWN"] } # 取三大R包上调基因交集: up = intersect(intersect(UP(DESeq2_DEG), UP(edgeR_DEG)), UP(limma_voom_DEG)) # intersect函数表示取两个向量的交集。 # 取三大R包下调基因交集: down = intersect(intersect(DOWN(DESeq2_DEG), DOWN(edgeR_DEG)), DOWN(limma_voom_DEG)) # 取三大R包差异基因交集画热图: dat = log2(exp+1) hp = draw_heatmap(dat[c(up, down), ], group_list, scale_before = T) #上调、下调基因分别画维恩图: up_genes = list(Deseq2 = UP(DESeq2_DEG), edgeR = UP(edgeR_DEG), limma = UP(limma_voom_DEG)) down_genes = list(Deseq2 = DOWN(DESeq2_DEG), edgeR = DOWN(edgeR_DEG), limma = DOWN(limma_voom_DEG)) up.plot <- draw_venn(up_genes, "UPgene") down.plot <- draw_venn(down_genes, "DOWNgene") #维恩图拼图 library(patchwork) pca.plot + hp + up.plot + down.plot + plot_layout(guides = "collect") ggsave(paste0(cancer_type, "_Heat.Ve.PCA.png"), width = 15, height = 10)