日常瞎掰
火山图作为展示差异基因的首选，可以说是生信分析常见的图形了。常规的火山图会在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语言书籍免费领
可视化：网络图