加载数据加载差异分析的结果 rm(list = ls()) load(file='limma差异分析结果.Rda') library(ggplot2) class(DEG_limma_voom) volcano<-DEG_limma_voom 数据处理我们先把差异分析结果增加一列,以区分上调下调和无差异的基因。我们以adj.P.Val < 0.05,logFC>2或者logFC<-2作为有明显差异的基因 volcano$type[(volcano$adj.P.Val > 0.05|volcano$adj.P.Val=="NA")|(volcano$logFC < 2)& volcano$logFC > -2] <- "none significant" volcano$type[volcano$adj.P.Val <= 0.05 & volcano$logFC >= 2] <- "up-regulated" volcano$type[volcano$adj.P.Val <= 0.05 & volcano$logFC <= -2] <- "down-regulated" 纵坐标一般取-log10,绘图看一下 p = ggplot(volcano,aes(logFC,-1*log10(adj.P.Val),color=type)) p + geom_point() 给横坐标设定界限,自定义颜色,加上临界值的虚线 x_lim <- max(volcano$logFC,-volcano$logFC) gg=p + geom_point( size=1) + xlim(-x_lim,x_lim) +labs(x="log2(FC)",y="-log10(adjP)")+ scale_color_manual(values =c("#A52A2A","grey","#f8766d"))+ geom_hline(aes(yintercept=-1*log10(0.05)),colour="black", linetype="dashed") + geom_vline(xintercept=c(-2,2),colour="black", linetype="dashed") print(gg) 看起来还不错
如何把差异最明显的基因名标在火山图上呢?这里提供了几种方法可供选择 方法一volcano1<-volcano volcano1$gene_name<-row.names(volcano1)
library(dplyr) top <- volcano1 %>% group_by(type) %>% top_n(n = 10, wt = -1*log10(adj.P.Val)) %>% filter(type !='none significant') # 取上下调中显著差异表达基因前10个
gg + geom_text(data=top,aes(x=logFC, y=-1*log10(adj.P.Val),label=as.character(gene_name)),size=3) 基因名重叠的可以下载PFD格式,在AI里面修改 如果想把背景的灰色去掉们也有办法 library(ggthemes) gg + geom_text(data=top,aes(x=logFC, y=-1*log10(adj.P.Val),label=as.character(gene_name)),size=3)+theme_few() 方法二
这样基因名就不会重叠在一起,还有箭头指示 方法三library(ggrepel) up <- subset(volcano1, type == 'up-regulated') up <- up[order(up$adj.P.Val), ][1:10, ] down <- subset(volcano1, type == 'down-regulated') down <- down[order(down$adj.P.Val), ][1:10, ] p2 <- gg + theme(legend.position = 'right') + geom_label_repel(data = rbind(up, down), aes(label = gene_name), show.legend = FALSE) print(p1) 比方法二多了基因的框,看起来也不错。 最后再介绍一个画热图的在线网站:http://www./ImageGP/index.php/Home/Index/Volcanoplot.html,只要准备好数据输入进去就行了,不过要注意,数据的行名和列名有要求,不能包含特殊字符,比如空格,逗号等,最好就只用字母表示就行。比如我们把数据整理成下面格式 粘贴到网站中,出图 参数可以调节,图片可以下载。 好了,今天的分享就到这了。 另外,最近收集了一些很好的资源,想分享给大家,顺便能涨一些粉,主要有
R语言实战(中文完整版) R数据科学(中文完整版) ggplot2:数据分析与图形艺术 30分钟学会ggplot2
前期从https:///datapages/ (UCSC Xena)数据库下载的TCGA数据,传到了百度云上备份。 感兴趣的话,转发朋友圈或者100人以上的微信群,截图发到公众号,即可获取全部资源的百度云链接,链接7天有效,希望大家赶紧下载。你们的支持是我前进的动力,感谢。 |
|