双曲线火山图一键拿捏

日常瞎掰

火山图作为展示差异基因的首选，可以说是生信分析常见的图形了。常规的火山图会在x、y轴方向上添加垂直参考线，以方便区分满足阈值的差异基因。常规的火山图这里就不多了，今天我们来说说双曲线火山图，也许该图没有那么高的出镜率，但其却有比较实用的价值。相对于常规火山图来说，双曲线火山图使用两条曲线作为阈值的参考线，如此更有利于筛选到更为真实的差异基因。那么，下面我们就来说说如何绘制双曲线火山图。

绘图

下面的代码里面定义了一个绘图函数，包括三个输入参数：数据框，pvalue阈值和foldchange阈值。其中，输入的数据框，必须包含log2foldchange和pvalue两列；pvalue_threshold为pvalue阈值，默认为0.05；foldchange_threshold为foldchange阈值，默认为1。代码如下：

library(ggplot2) volcano_plot <- function(df, pvalue_threshold = 0.05, foldchange_threshold = 1) { xmax <- max(abs(na.omit(df$log2foldchange))) + 0.2 xmin <- min(abs(na.omit(df$log2foldchange)), 0.0001) x <- seq(xmin, xmax, by = 0.0001) y <- 1/x + (-log10(pvalue_threshold)) curve_xy <- rbind(data.frame(xpos = x + foldchange_threshold, ypos = y), data.frame(xpos = -(x + foldchange_threshold), ypos = y)) df$curve_y <- ifelse(df$log2foldchange > 0, 1/(df$log2foldchange - foldchange_threshold) + (-log10(pvalue_threshold)), 1/(-df$log2foldchange - foldchange_threshold) + (-log10(pvalue_threshold))) df$curve_group <- ifelse(-log10(df$pvalue) > df$curve_y & df$log2foldchange > foldchange_threshold, 'up', ifelse(-log10(df$pvalue) > df$curve_y & df$log2foldchange < -foldchange_threshold, 'down', 'nosignif')) df$pvalue <- -log10(df$pvalue) p <- ggplot(df, aes(x = log2foldchange, y = pvalue, color = curve_group)) + geom_point(size = 1) + geom_line(data = curve_xy, aes(x = xpos, y = ypos), lty = 3, col = "black", lwd = 0.6) + scale_color_manual(values = c('up'='red', 'down'='blue', 'nosignif'='gray')) + xlim(-xmax, xmax) + ylim(0, 30) + labs(x = "log2(FoldChange)", y = "-log10(P-value)") + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.spacing.x = unit(0.05, 'cm'), plot.title = element_text(hjust = 0.5), legend.text = element_text(size = 8)) + guides(color = guide_legend(override.aes = list(size = 2), title = NULL)) return(p) } data <- read.table('desktop/sample_dge.txt',header=T,stringsAsFactors=F,sep='\t') head(data) gene pvalue log2foldchange 1 ENSG00000000003.15 TSPAN6 6.954955e-04 1.0305811 2 ENSG00000000005.6 TNMD 1.103522e-01 -2.1289526 3 ENSG00000000419.12 DPM1 7.168680e-02 0.5515042 4 ENSG00000000457.14 SCYL3 5.743836e-01 0.1453620 5 ENSG00000000460.17 C1orf112 1.173320e-06 2.1643651 6 ENSG00000000938.13 FGR 1.388476e-13 -4.0345022 p <- volcano_plot(data) p

结果如下：

结束语

双曲线火山图绘制的关键，就是根据反比例函数确定参考线的坐标位置以及差异基因的定义。上面的绘图代码为了方便起见所以采用了硬编码的方式，所以对输入的数据框格式要求稍微严格一些，需要含有log2foldchange和pvalue两列，且列名也要保持一致，有没有其他的列并没有影响。当然了，上面的火山图没有包含标记差异基因的功能，有需要的话可以参考常规火山图添加基因名注释的方法。哦了，今天就到这里了~~~