逐渐感觉到,一些人正在将生信分析带偏,他们着重强调图片的观赏性,而忽视了生信分心的本质不是为了文章凑图片,而是能利用大数据分析观测到wet lab实验不能观测的新视角,提供新发现以供实验验证。火山图作为一款常见的结果展示图片似乎正在被玩坏。网上各种络绎不绝的帖子教你如何调整颜色,调整点的大小变化。我不否认漂亮的图片可以为你的文章增色,然而前提必须是建立在科学和准确的基础上的。
火山图是一个很好的用来显示差异分析结果的方式,通常情况下Y轴代表的是显著性的p值,为了方便观察显著性的,通常用-log10(p.val)或者-log10(p.adj),X轴数值一般是差异倍数FC(fold change)或log2FC。 举个例子,下图是我根据自己的蛋白质组数据用ggplot做的图。
你想要的是不是都在?
该图展示了以下几种作图的操作:
1.FC变化与点的大小成正相关.FC越大,点越大
2.颜色变化与p值的大小而变化,p值越小,颜色越深(渐变色),且左右两侧使用不同颜色。
3.使用SAM test将T-test 结果进行显著性区分,而非人为划定(FC=1,-log10(p.val)<0.05,虚线所示)。
4.在图中对感兴趣基因进行标注
为什么做这个图呢?最近注意到有些人开始传授如何在火山图画这个双曲线。这个双曲线在R里其实是个相对复杂的操作,并且牵扯到SAM test SAM (Significance Analysis of Microarray)。这个分析最初的参考文献见这里:Proc Natl Acad Sci USA. 2001 Apr 24;98(9):5116-21. SAM test 其实是一个 一种改进的T-tets,其引了Fudge factor 用以校正在t-test里一些差异较小的基因的假阳性问题。据我所知的其在火山图中的应用最常见的地方其实是来源于Perseus 软件里的火山图,其使用的统计学原理也是来自我上面提到的这篇文献。
由于这个双曲线和反比例函数曲线近似,于是就出现了我觉得比较荒唐的地方,有人开始用反比例函数来画这个图。这个倒也不是原创,研究了一下发现确实有人这么在文章里干过,例如: Mol Cell Proteomics. 2015 Jan;14(1):120-35. 文章里使用了一个自创的函数y=c/(x-x0),(c= curvature,x0 = minimum fold change) ,参数的取值似乎也是猜出来的。在我看来这个逆向工程有点不讲武德,为了区分显著性的蛋白就自己琢磨了一个分界线,还搞得有点玄幻。其本质上来讲,似乎还不如传统的-Log10(p.val)加FC分解来的有道理,虽然都是猜。所以问题来了,既然都是猜,是不是就可以乱来?是不是搞个指数函数也可以呢?因为如果你用 geom_function(fun = ~ 0.5*exp(-abs(.x))) 也能画出类似的曲线,如下图
geom_function(fun = ~ 0.5*exp(-abs(.x)))
我觉得答案显然是你不能想怎样就怎样,就如我一开始所言,数据分析是要有科学依据的,你做的是科学不是艺术,准确度的优先级一定是大于美学的。如果实在不确定怎么划界,那么就老实的用别人的工具,或者传统的-Log10(p.val)加FC。当然如果你有充分的理由不这么做的除外。
最后附上作图部分的代码,水平有限,仅供参考。
ggplot() +
geom_point(dat_red, mapping = aes(x = FC, y = -log10(p.val), color = -log10(p.val)), size = dat_red$FC) +
scale_size(limits = c(0, 6)) +
scale_colour_gradient(low = "#e1bb4e", high = "#c12131") +
guides(color = "none") +
ggnewscale::new_scale_color() +
geom_point(dat_blue, mapping = aes(x = FC, y = -log10(p.val), color = -log10(p.val)), size = abs(dat_blue$FC)) +
scale_size(limits = c(0, 5)) +
scale_colour_gradient(low = "#a7d3d5", high = "#1c0d51") +
ggrepel::geom_text_repel(data = lab_gene, max.overlaps = 100, min.segment.length = 0.1, force = 0.5, size = 2, colour = "black", aes(x = FC, y = -log10(p.val), label = rownames(lab_gene))) +
geom_point(dat_gray, mapping = aes(x = FC, y = -log10(p.val)), size = abs(dat_gray$FC), color = "gray") +
geom_vline(xintercept = c(-1, 1), lty = 2, col = "gray", lwd = 0.5) +
geom_hline(yintercept = -log10(0.05), lty = 2, col = "gray", lwd = 0.5) +
geom_line(aes(x = x, y = y), color = "#dd008f", scurve, inherit.aes = FALSE, linetype = "dotted", size = 0.8) +
ylim(c(0, 6)) +
ggtitle("SAM Test Soft Threshold") +
xlab("log2(FC)") +
ylab("-log10(p.val)") +
theme_classic() +
guides(color = "none")