RNA-Seq(7):对差异分析结果绘制火山图进行可视化
2022-01-11 本文已影响0人
Z_bioinfo
1.既然知道有多少基因上下调表达了，但是结果并不直观，对他们做火山图可视化
1.将比较所产生的差异情况反映到火山图中，灰色为非显著性差异的基因，红色和绿色为显著性差异基因；X轴为log2 FoldChange的展示，Y轴方向为-log10 Padj的,只需要其中的logFC和padj就可以了。在绘图之前，我们需要对padj进行转换，将它的值变成-1 * log10，这样的话可以拉开差异表达基因之间的间距。
install.packages("ggpubr")
install.packages("ggthemes")
library(ggplot2)
library(ggpubr)
library(ggthemes)
>
#对 Padj进行log转换
res$logp <- -log10(res$padj)
#绘制基本热图
res = data.frame(res)将res转为表格形式，不然后面会报错
ggscatter(res,x= "log2FoldChange", y="logp") +theme_base()
image.png
2.对火山图进行美化，美化一下，首先需要区分其中哪些是显著差异表达基因。因此，我们需要对logFC和padj两列进行过滤。我们设置的过滤的条件为，padj小于0.05并且logFC大于2（4倍差异）为显著上调差异表达基因，padj小于0.05并且logFC小于-2（4倍差异）为显著下调差异表达基因。大家可以根据实验结果，适当调整logFC。
新加一列type，将padj小于0.05，logFC大于2（4倍差异）为显著上调差异表达基因
res$type <- ifelse(res$padj < 0.05,ifelse(abs(res$log2FoldChange) > 2, ifelse(res$log2FoldChange < -2,'down','up'),'noSig'),'noSig')
查看结果
> table(res$type)
down noSig up
103 14789 123
改变火山图的颜色（palette）和大小(size)，使用geom_hline和geom_vline分别添加横向和纵向的辅助线,并对坐标轴标注
ggscatter(res,x= "log2FoldChange", y="logp",color="type", palette=c("green", "gray", "red"),size=1, xlab="log2FoldChange", ylab="-log10(padj)") +theme_base()+geom_hline(yintercept=-log(0.05,10), linetype="dashed")+geom_vline(xintercept=c(-2,2), linetype="dashed")
或者
ggplot(res,aes(x = log2FoldChange,y = -log10(padj)))
+geom_point(aes(color = type),alpha = 0.5,size = 2) +
theme_bw(base_size = 16) + theme(aspect.ratio = 1, plot.title =
element_text(hjust = 0.5), axis.text = element_text(color = 'black'),
axis.title = element_text(color = 'black')) +scale_color_manual(name = '', values = c('up'='#DA1212','noSig'='grey','down'='#3E7C17'), label = c('up'='up (123)','noSig'='noSig (14789)','down'='down (103)')) + geom_hline(yintercept = -log10(0.05),lty = 'dashed',size = 0.8) +geom_vline(xintercept = c(-2,2),lty = 'dashed',size = 0.8) +
scale_x_continuous(breaks = c(-6,-4,-2,-1,0,1,2,4,6))
image.png
3.如果想在火山图上给差异表达的基因加标签，可以继续往下做
新加一列label
res$label = ""
对padj进行排序
res = res[order(res$padj), ]
选择表达水平高的前10个基因名字
up.gene = head(res$external_gene_name[which(res$group=="up")],10)
down.gene = head(res$external_gene_name[which(res$group=="down")],10)
将up.gene 和down.gene合并，添加到label中
res.top10.gene =c(as.character(up.gene),as.character(down.gene))
res$label[match(res.top10.gene,res$external_gene_name)]=res.top10.gene
> ggscatter(res,x= "log2FoldChange", y="logp",color="group", palette=c("#2f5688", "#BBBBBB", "#CC0000"),size=1, label=res$label, font.label=8, repel=T, xlab="log2FoldChange", ylab="-log10(p-value)") +theme_base()+geom_hline(yintercept=-log(0.05,10), linetype="dashed")+geom_vline(xintercept=c(-2,2), linetype="dashed")
image.png
image.png
