火山图是散点图的一种,它将统计测试中的统计显著性量度(如 p value)和变化幅度 FC 相结合,从而能够帮助快速直观地识别那些变化幅度较大且具有统计学意义的数据点(基因等)。常应用于转录组、基因组、蛋白质组、代谢组等统计数据。
1)什么是 fold change?
翻译成中文是差异倍数,简单来说就是基因在一组样品中的表达值的均值除以其在另一组样品中的表达值的均值。所以火山图只适合展示两组样品之间的比较。
2)为什么要做 Log2 转换?
两个数相除获得的结果(fold change)要么大于 1,要么小于 1,要么等于
1。那么对应于基因差异呢?简单说,大于 1
表示上调(可以描述为上调多少倍),小于 1
表示下调(可以描述为下调为原来的多少分之多少)。大于 1
可以到多大呢?多大都有可能。小于 1 可以到多小呢?最小到 0。用原始的 fold
change
描述上调方便,描述下调不方便。绘制到图中时,上调占的空间多,下调占的空间少,展示起来不方便。所以一般会做
Log2 转换。
3)什么是 P-value?
统计检验获得的是否统计差异显著的一个衡量值,约定成俗的 P-value < 0.05
为统计检验显著的常规标准。
本文讨论如何绘制火山图以及如何对其进行解读。
必须包含 3 列数据(建议有行名):名称列、FC 列、PValue
列。Marker
列为非必须。如果下方参数中选择
Marker
为标记方式,则需要有此列数据。
下方代码为你提供的原始逻辑,已整理进可执行的 R Markdown 代码块(未改动算法与参数)。
## ========= 依赖 =========
need_pkgs <- c("ggplot2","ggrepel","dplyr")
to_install <- setdiff(need_pkgs, rownames(installed.packages()))
if (length(to_install) > 0) install.packages(to_install, dependencies = TRUE)
library(ggplot2); library(ggrepel); library(dplyr)
dat_f <- read.delim(
"~/Documents/小红书/火山图/volcano_filtered.tsv",
sep = "\t", header = TRUE, check.names = FALSE,
stringsAsFactors = FALSE, quote = "", comment.char = ""
)
head(dat_f)
## Name FC PValue Marker P_safe log2FC negLog10P
## 1 ID1 0.138995 0.007489 0 0.007489 -2.8468951 2.1255762
## 2 ID2 1.192851 0.801349 0 0.801349 0.2544138 0.0961783
## 3 ID3 0.777325 0.001722 0 0.001722 -0.3634102 2.7639669
## 4 ID4 2.416923 0.039086 1 0.039086 1.2731715 1.4079788
## 5 ID5 0.839945 0.746802 0 0.746802 -0.2516332 0.1267945
## 6 ID6 1.804295 0.001731 1 0.001731 0.8514352 2.7617029
FC_cut <- 1.5
P_cut <- 0.05
method <- "IQR" # 可选: "IQR" 或 "quantile"
q_tail <- 0.995 # method="quantile" 时,保留 99.5%(去掉最极端 0.5%)
## ========= 显著性分组与Top10标注 =========
dat_f <- dat_f %>%
mutate(
sig = case_when(
negLog10P >= -log10(P_cut) & log2FC >= log2(FC_cut) ~ "Up",
negLog10P >= -log10(P_cut) & log2FC <= -log2(FC_cut) ~ "Down",
TRUE ~ "NotSig"
),
sig = factor(sig, levels = c("NotSig","Up","Down"))
)
label_df <- dat_f %>%
dplyr::filter(sig %in% c("Up","Down")) %>%
mutate(score = negLog10P * abs(log2FC)) %>%
arrange(desc(score)) %>%
slice_head(n = 10)
dat_f <- dat_f %>% mutate(label = ifelse(Name %in% label_df$Name, Name, ""))
## ========= 配色与作图 =========
cols <- c("NotSig"="#BDBDBD","Up"="#E64B35","Down"="#4DBBD5")
set.seed(1)
p <- ggplot(dat_f, aes(x = log2FC, y = negLog10P)) +
geom_point(aes(color = sig), alpha = 0.85, size = 1.6) +
geom_hline(yintercept = -log10(P_cut), linetype = 2, linewidth = 0.4) +
geom_vline(xintercept = c(-log2(FC_cut), log2(FC_cut)), linetype = 2, linewidth = 0.4) +
ggrepel::geom_text_repel(
data = subset(dat_f, label != ""),
aes(label = label),
box.padding = unit(0.5,"lines"),
point.padding = unit(0.15,"lines"),
max.overlaps = Inf,
min.segment.length = 0,
segment.color = "black",
size = 3, show.legend = FALSE
) +
labs(title="Volcano Plot (filtered)", x="log2(FC)", y=expression(-log[10](P))) +
scale_color_manual(values = cols, name = NULL) +
coord_cartesian(clip = "off") +
theme_classic(base_size = 12) +
theme(
panel.background = element_rect(fill = "transparent", colour = NA),
plot.background = element_rect(fill = "transparent", colour = NA),
legend.background = element_rect(fill = "transparent", colour = NA),
legend.key = element_rect(fill = "transparent", colour = NA),
plot.margin = margin(5.5, 60, 5.5, 5.5, unit = "pt")
)
print(p)
需要数据分析联系xhs:飞高高啦