本节概览: 1.DESeq2、 edgeR、limma的使用 2.三类差异分析软件的结果比较——相关性、韦恩图 3.选取差异基因绘制火山图和热图 一、DESeq2、 edgeR、limma的使用强烈建议查看官方说明书进行这三种差异分析的学习,链接在文章末尾给出。 注意,这三个包都需要输入counts进行分析,不能用tpm、fpkm等归一化后的数据。 承接上节 RNA-seq入门实战(三):在R里面整理表达量counts矩阵 和 RNA-seq入门实战(二):上游数据的比对计数——Hisat2+ featureCounts 与 Salmon 正式分析前先进行目录设置、实验组和对照组的指定: rm(list = ls()) options(stringsAsFactors = F)
setwd("C:/Users/Lenovo/Desktop/test") load(file = '1.counts.Rdata') dir.create("3.DEG") setwd("3.DEG")
##设定 实验组exp / 对照组ctr exp="primed" ctr="naive"
1. DESeq2DESeq2是目前最常用的差异分析R包。除了可以导入counts外,如果上游使用salmon,DESeq2官方还给出了直接导入tximport生成的txi对象的方法。counts与txi的获取见 RNA-seq入门实战(三):在R里面整理表达量counts矩阵 和 RNA-seq入门实战(二):上游数据的比对计数——Hisat2+ featureCounts 与 Salmon library(DESeq2) library("BiocParallel") #启用多核计算
##构建dds DESeqDataSet if(T){ dds <- DESeqDataSetFromMatrix(countData = counts, colData = gl, design = ~ group_list) } if(F){ #若上游为salmon dds <- DESeqDataSetFromTximport(txi, colData = gl, design = ~ group_list) }
dds$group_list <- relevel(dds$group_list, ref = ctr) #指定 control group
keep <- rowSums(counts(dds)) >= 1.5*ncol(counts) #Pre-filtering ,过滤低表达基因 dds <- dds[keep,] dds <- DESeq(dds,quiet = F) res <- results(dds,contrast=c("group_list", exp, ctr)) #指定提取为exp/ctr结果 resOrdered <- res[order(res$padj),] #order根据padj从小到大排序结果 tempDEG <- as.data.frame(resOrdered) DEG_DEseq2 <- na.omit(tempDEG)
2. edgeR使用EdgeR时注意选择合适的分析算法,官方建议bulk RNA-seq选择quasi-likelihood(QL) F-test tests,scRNA-seq 或是没有重复样品的数据选用 likelihood ratio test。 library(edgeR) #install.packages("statmod") library(statmod)
#分组矩阵design构建 group <- factor(group_list) group <- relevel(group,ctr) #将对照组的因子设置为1 design <- model.matrix(~0+group) rownames(design) <- colnames(counts) colnames(design) <- levels(group)
## 表达矩阵DGEList构建与过滤低表达基因 dge <- DGEList(counts=counts, group=group) keep.exprs <- filterByExpr(dge) #自动筛选过滤低表达基因 dge <- dge[keep.exprs,,keep.lib.sizes=FALSE] dge <- calcNormFactors(dge, method = 'TMM') #归一化因子用于 normalizes the library sizes dge <- estimateDisp(dge, design, robust=T)
## To perform quasi-likelihood(QL) F-test tests: bulk RNA-seq fit <- glmQLFit(dge, design, robust=T) #拟合模型 lt <- glmQLFTest(fit, contrast=c(-1,1)) #统计检验#注意比对顺序:实验-1 /对照1
## To perform likelihood ratio test:scRNA-seq and no replicates data # fit <- glmFit(dge, design, robust=T) # lt <- glmLRT(fit, contrast=c(-1,1)) #比对:顺序实验/对照,已设对照为1
tempDEG <- topTags(lt, n = Inf) #sort by PValue, n is the number of genes/tags to return tempDEG <- as.data.frame(tempDEG) DEG_edgeR <- na.omit(tempDEG)
3. limmalimma进行差异分析有两种方法 :limma-trend和voom,在样本测序深度相差不大时两种方法差距不大,而测序深度相差大时voom更有优势,因此我们一般都选择voom的方法进行差异分析。Limma-voom vs limma-trend (bioconductor.org) library(limma) library(edgeR)
#分组矩阵design构建 design <- model.matrix(~0+factor(group_list)) #构建分组矩阵 colnames(design) <- levels(factor(group_list)) rownames(design) <- colnames(counts)
## 表达矩阵DGEList构建与过滤低表达基因 dge <- DGEList(counts=counts) keep.exprs <- filterByExpr(dge,design=design) #过滤低表达基因 dge <- dge[keep.exprs,,keep.lib.sizes=FALSE] dge <- calcNormFactors(dge) #归一化基因表达分布,得到的归一化系数被用作文库大小的缩放系数 cont.matrix <- makeContrasts(contrasts=paste0(exp,'-',ctr), #比对顺序实验/对照 levels = design)
## DE分析 limma-trend(logCPM,有相似文库大小) or voom(文库大小差异大) # de <- cpm(dge, log=TRUE, prior.count=3) #如选择logCPM,则eBayes设trend=TRUE de <- voom(dge,design,plot=TRUE, normalize="quantile") fit1 <- lmFit(de, design) #线性拟合 fit2 <- contrasts.fit(fit1,cont.matrix) #统计检验 efit <- eBayes(fit2, trend=F) #Apply empirical Bayes smoothing to the standard errors
tempDEG <- topTable(efit, coef=paste0(exp,'-',ctr), n=Inf) #padj值从小到大排列 DEG_limma_voom <- na.omit(tempDEG)
4. 保存DEG数据#### 保存DEG数据 ####save(DEG_limma_voom,DEG_DEseq2,DEG_edgeR, file ='test_DEG_results.Rdata') #### 合并三种DEG文件交集至csv allg <- intersect(rownames(DEG_limma_voom),rownames(DEG_edgeR))#取交集 allg <- intersect(allg,rownames(DEG_DEseq2)) ALL_DEG <- cbind(DEG_limma_voom[allg,c(1,4,5)], DEG_edgeR[allg,c(1,4,5)], DEG_DEseq2[allg,c(2,5,6)]) colnames(ALL_DEG) colnames(ALL_DEG) <- c('limma_log2FC','limma_pvalue','limma_padj', 'edgeR_log2FC','edgeR_pvalue','edgeR_padj', 'DEseq2_log2FC','DEseq2_pvalue','DEseq2_padj') write.csv(ALL_DEG,file = 'test_DEG_results.csv')
二、3种差异分析结果比较由于本次样品两组间差异十分显著,差异基因很多,因此筛选阈值调整为:FoldChang=10,padj=0.001。一般情况下选择FoldChang=1.5~4,padj<=0.05即可,根据样本情况而定。 下面查看三种差异分析结果的相关性和差异基因的重叠情况。 #### 比较一下三种DEG方法结果 #### load(list.files(pattern = 'DEG_results.Rdata')) a <- read.csv(file = list.files(pattern = 'DEG_results.csv'),header = T,row.names = 1)
#查看相关性、一致性 colnames(a) print(cor(a[,c(1,4,7)]))
#查看显著差异基因重叠性,绘制韦恩图 #BiocManager::install("RBGL") #安装依赖包 #install.packages("Vennerable", repos="http://R-Forge.R-project.org") #安装Vennerable包 library(Vennerable)
log2FC_cutoff=log2(10); p_cutoff=0.001 #筛选显著差异基因比较 if(T){#根据Padj筛选 DEseq2_deg <- rownames(DEG_DEseq2[with(DEG_DEseq2,abs(log2FoldChange)>log2FC_cutoff & padj<p_cutoff),]) edgeR_deg <- rownames(DEG_edgeR[with(DEG_edgeR,abs(logFC)>log2FC_cutoff & FDR<p_cutoff),]) limma_deg <- rownames(DEG_limma_voom[with(DEG_limma_voom,abs(logFC)>log2FC_cutoff & adj.P.Val<p_cutoff),]) }
mylist <- list(DEseq2=DEseq2_deg, edgeR=edgeR_deg, limma=limma_deg) str(mylist) Vennplot <- Venn(mylist) Vennplot1 <- Vennplot[,c('DEseq2','edgeR')] Vennplot2 <- Vennplot[,c('DEseq2','limma')] Vennplot3 <- Vennplot[,c('limma','edgeR')]
pdf(file = paste0('3DEG_Vennplot_lg2FC',log2FC_cutoff,'.pdf')) plot(Vennplot, doWeights = F) plot(Vennplot, doWeights = T) #doWeights=T设置为按数量比例绘图 plot(Vennplot1, doWeights = T) plot(Vennplot2, doWeights = T) plot(Vennplot3, doWeights = T) dev.off()
从以下图中可以看出,三种方法所得结果相关性很高,都在0.99以上,显著差异基因的重叠也很高。(所以一般来说大家无需纠结使用哪种方法,都是认可的)
三、选取差异基因绘制火山图和热图以下示范选取DESeq2差异分析结果进行绘制, 筛选阈值设置为:FoldChang=10,padj=0.001 1. 火山图的绘制library(ggplot2) library(pheatmap) ##筛选条件设置 log2FC_cutoff = log2(10) padj_cutoff = 0.001 ##选取差异分析结果 need_DEG <- DEG_DEseq2[,c(2,6)] #选取log2FoldChange, padj信息 colnames(need_DEG) <- c('log2FoldChange','padj')
need_DEG$significance <- as.factor(ifelse(need_DEG$padj < padj_cutoff & abs(need_DEG$log2FoldChange) > log2FC_cutoff, ifelse(need_DEG$log2FoldChange > log2FC_cutoff ,'UP','DOWN'),'NOT'))
title <- paste0(' Up : ',nrow(need_DEG[need_DEG$significance =='UP',]) , '\n Down : ',nrow(need_DEG[need_DEG$significance =='DOWN',]), '\n FoldChange >= ',round(2^log2FC_cutoff,3))
g <- ggplot(data=need_DEG, aes(x=log2FoldChange, y=-log10(padj), color=significance)) + #点和背景 geom_point(alpha=0.4, size=1) + theme_classic()+ #无网格线 #坐标轴 xlab("log2 ( FoldChange )") + ylab("-log10 ( P.adjust )") + #标题文本 ggtitle( title ) + #分区颜色 scale_colour_manual(values = c('blue','grey','red'))+ #辅助线 geom_vline(xintercept = c(-log2FC_cutoff,log2FC_cutoff),lty=4,col="grey",lwd=0.8) + geom_hline(yintercept = -log10(padj_cutoff),lty=4,col="grey",lwd=0.8) + #图例标题间距等设置 theme(plot.title = element_text(hjust = 0.5), plot.margin=unit(c(2,2,2,2),'lines'), #上右下左 legend.title = element_blank(), #不显示图例标题 legend.position="right") #图例位置
ggsave(g,filename = 'volcano_padj.pdf',width =8,height =7.5)
2. 热图的绘制##选择要展示基因表达量的数据 # dat <- log2(edgeR::cpm(counts)+1) dat <- log2(tpm+1) # dat <- read.table("../2.check/Deseq2_rld.txt"); colnames(dat) <- rownames(gl) #R读取数据列名可能会出错,需要重新对应一下
gene_up <- rownames(need_DEG[with(need_DEG,log2FoldChange>log2FC_cutoff & padj<padj_cutoff),]) gene_down <- rownames(need_DEG[with(need_DEG,log2FoldChange< -log2FC_cutoff & padj<padj_cutoff),]) cg <- c(head(gene_up, 50), #取前50 padj上下调基因名 head(gene_down, 50)) cg <- na.omit(match(cg,rownames(dat)))
pheatmap::pheatmap(dat[cg,], show_colnames =T, show_rownames = F, #scale = "row", fontsize = 7 , cluster_cols = T, annotation_col=gl, filename = 'heatmap_top50up&down_DEG.pdf')
到此就完成了基因差异分析的基本过程,得到了不同分组间的差异基因相关信息,接下来就要对差异基因进行富集分析啦。 目前简单的差异分析流程,基本上转录组测序技术和芯片技术拿到的表达量矩阵后续分析大同小异,公众号推文在: 参考资料 三种R包的官方说明书: Analyzing RNA-seq data with DESeq2 (bioconductor.org) edgeR: differential analysis of sequence read count data User's Guide (bioconductor.org) limma usersguide.pdf (bioconductor.org) bioconductor的worlk flow: 使用limma、Glimma和edgeR,RNA-seq数据分析易如反掌 (bioconductor.org) 部分代码修改参考自: GitHub - jmzeng1314/GEO 【生信技能树】转录组测序数据分析_哔哩哔哩_bilibili 【生信技能树】GEO数据库挖掘_哔哩哔哩_bilibili
|