日常瞎掰
火山图作为展示差异基因的首选，可以说是生信分析常见的图形了。常规的火山图会在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两列，且列名也要保持一致，有没有其他的列并没有影响。当然了，上面的火山图没有包含标记差异基因的功能，有需要的话可以参考常规火山图添加基因名注释的方法。哦了，今天就到这里了~~~
往期回顾
ChIP-seq数据质控
ChatGPT！见证AI的力量！
ChIPseeker绘图函数借用
R语言书籍免费领
可视化：网络图
