manhattan · Folio

什么时候用它

GWAS 或基因水平关联分析跑完，手上是几十万到几百万行检验结果。Manhattan plot 把它们压进一张图：横轴按染色体和物理位置排，纵轴是 -log10(P)——P 越小点越高，显著信号就像高楼一样从背景里立出来。

它回答的是"哪里值得跟进"，不回答"效应有多大"。

本质上它就是一张散点图（位置 vs -log10(P)），只是横轴按染色体分了段——所以这里归在散点那一族，而不是单开一档。

怎么读

先看 y 轴。 是 -log10(P)，不是效应量。点越高只说明 P 越小。
横轴不是连续坐标轴。 它按染色体分段，段内才是物理位置。相邻两点可能隔着几百万碱基。
看阈值线。 超过线的点是候选。用的是哪条线（Bonferroni / FDR / suggestive）必须在图注里说清楚，三者含义完全不同。
看峰不看点。 同一区域连续多个点一起升高，比一个孤立高点可信得多——后者常常是基因分型错误。
标签只是候选。 标注出来的位点还需要后续精细定位和功能解释。

常见陷阱

高点不等于效应大。 样本量足够时，一个很小的效应也能给出极小的 P。
阈值必须说明。 5e-8 是全基因组显著性的惯例（约等于 100 万次独立检验的 Bonferroni 校正），1e-6 只是提示性水平。
标签别标太多。 标注超过十来个就会互相遮挡，反而看不见信号峰。
P 值必须大于 0。 有 0 或缺失时 -log10(P) 会出 Inf 或 NA，整张图会歪。画之前先查。
染色体和位置列错了图就没意义。 而且不会报错——它照样画得出来，只是排序全乱。
CMplot() 默认会自己往磁盘写文件。 file.output 的默认值是 TRUE，不设就会在工作目录里落下一堆图。上面两段都关掉了它，是因为要把图画到当前设备上；但真要导出时反过来该用它——file = "pdf"、dpi = 300、file.name = "..." 交给 CMplot 自己处理，它会按图的内容挑合适的画布尺寸。手动 png() → CMplot() → dev.off() 也行，但尺寸得自己试。

配色

这条 recipe 使用站内定性色板 cancer_mosaic。深蓝与橙色交替区分相邻染色体，红色标出 Top hits，粉色表示 genome-wide 阈值与信号，亮蓝色表示 suggestive 阈值与信号。这个组合接近 Manhattan plot 常见的蓝橙背景与红粉蓝信号层次，矩形与圆形布局保持一致。

Manhattan plot 的颜色不是五个平级分类：背景染色体、重点位点和两级阈值承担不同语义。因此这里显式按名称分配色板中的五个位置，而不是把色板直接循环到所有元素上；图例和标签仍是判断信号级别的必要依据。

圆形图最外圈的标记密度是连续量，不再混用多种分类色。它使用 cancer_mosaic 深蓝色派生的浅蓝至深蓝梯度：颜色越深表示局部标记越密。梯度截掉了最接近白色的约 42%，让最低非零档从清楚可见的浅蓝开始，也避免把密度误读成额外的类别或信号等级。

配方

方法	绘图系统	什么时候选它
ACMplot(plot.type = "m")	base graphics	要标注位点、比较阈值、写进正文——绝大多数情况用这个
BCMplot(plot.type = "c")	base graphics	版面只有一个方块，或者只想给一眼全基因组的密度概览

数据准备

两种画法用的是同一份数据、同一套阈值，所以先备好。

library(CMplot)
library(dplyr)
library(biopalette)

gwas_data <- read.csv("content/tessera/data/csv/gwas_manhattan.csv")

threshold_line <- c(5e-8, 1e-6)              # 全基因组显著 / 提示性
manhattan_cols <- setNames(
  get_palette("cancer_mosaic", type = "qualitative")[c(1, 3, 9, 7, 13)],
  c("odd", "even", "genome_wide", "top", "suggestive")
)
threshold_cols <- unname(manhattan_cols[c("genome_wide", "suggestive")])
density_cols <- colorRampPalette(
  c("white", unname(manhattan_cols["odd"]))
)(171)[72:171]  # 以旧图约 4 档的浅蓝作为最低非零密度色

top_hits <- gwas_data |> arrange(P) |> slice_head(n = 8)   # 只标 Top 8

真实数据只要整理成 Marker / Chr / Pos / P 这四列，就能直接替换掉 gwas_data。

方法 A · 直线型

#| fig: linear
#| fig-width: 12
#| fig-height: 6
par(mar = c(5, 5, 4, 2))   # 下、左边距留大：给轴标题腾地方

CMplot(
  gwas_data,
  type          = "p",         # 画点；"l" 是连线
  plot.type     = "m",         # m = manhattan（直线型）
  col           = unname(manhattan_cols[c("odd", "even")]),
  LOG10         = TRUE,        # 自动取 -log10(P)；已经转过就设 FALSE
  cex           = 0.7,         # 背景点大小
  band          = 0.5,         # 染色体之间的留白宽度
  threshold     = threshold_line,
  threshold.lty = c(2, 2),     # 线型，2 = 虚线
  threshold.col = threshold_cols,
  amplify       = TRUE,        # 超阈值的点放大重画一遍 —— 信号才跳得出来
  signal.col    = unname(manhattan_cols[c("genome_wide", "suggestive")]),
  signal.cex    = 1.2,
  signal.pch    = 19,          # 实心圆
  highlight     = top_hits$Marker,      # 要额外高亮的标记名
  highlight.col = unname(manhattan_cols["top"]),
  highlight.text = top_hits$Marker,     # 高亮点旁边写什么
  highlight.text.cex = 0.7,
  main          = "Simulated GWAS Manhattan Plot",
  main.cex      = 1.6,
  ylim          = c(0, max(-log10(gwas_data$P)) * 1.2),  # 顶上留 20%，不然标签被裁
  ylab          = "",          # 留空，下面用 mtext 自己画（要带下标排版）
  chr.border    = FALSE,
  file.output   = FALSE,       # 关键：不让 CMplot 自己往磁盘写文件
  verbose       = FALSE,
  box           = TRUE
)

# CMplot 的 ylab 不支持 expression，所以轴标题自己补
mtext("Chromosome Position", side = 1, line = 2.5, cex = 1.2, font = 2)
mtext(expression(-log[10](P)),  side = 2, line = 2,   cex = 1.2, font = 2)

# 阈值线画出来是两条无标识的虚线，必须自己补图例 —— 否则读者根本分不清
# 哪条是全基因组显著、哪条是提示性水平，而这两者的含义差得很远。
# 用 legend() 而不是 text() + par("usr") 手算坐标：后者换个 ylim 就跑位。
legend(
  "topleft",
  legend = c("Genome-wide (5e-8)", "Suggestive (1e-6)"),
  col    = threshold_cols,
  lty    = 2,
  bty    = "n",    # 不画外框，免得盖住背景点
  cex    = 0.9
)

方法 B · 环形

染色体沿圆周排列，在有限版面里塞下整个基因组。外圈那条彩带是标记密度。

#| fig: circular
#| fig-width: 8
#| fig-height: 8
par(mar = c(1, 1, 2, 1))   # 环形图四边都不需要轴，边距压到最小

CMplot(
  gwas_data,
  type          = "p",
  plot.type     = "c",         # c = circular
  LOG10         = TRUE,
  cex           = 0.45,        # 环形上点更密，要比直线型小
  chr.labels    = paste0("Chr", sort(unique(gwas_data$Chr))),
  col           = unname(manhattan_cols[c("odd", "even")]),
  r             = 0.1,         # 内圈半径 —— 太大中间空洞，太小点会挤在一起
  cir.band      = 0.8,         # 染色体之间的角度间隙
  cir.axis      = TRUE,
  cir.axis.col  = unname(manhattan_cols["odd"]),
  cir.chr.h     = 0.8,         # 最外层染色体带的厚度
  cir.axis.grid = TRUE,
  threshold     = threshold_line,
  threshold.lty = c(1, 2),
  threshold.col = threshold_cols,
  amplify       = TRUE,
  signal.col    = unname(manhattan_cols[c("genome_wide", "suggestive")]),
  signal.line   = 1,           # 从圆心往信号点拉一条线，帮助定位
  chr.den.col   = density_cols,
  file.output   = FALSE,
  verbose       = FALSE
)

manhattan — circular — manhattan-circular

两种方法怎么选

绝大多数情况用 A。 要标注位点、要比较阈值线、要写进正文配图，直线型都更合适——它有一条完整的横轴可以读坐标。
B 只在两种场合更好：版面被限制成方形（比如放在多图拼版的一格里），或者你想传达的是"信号在全基因组上的分布密度"而不是"具体哪几个位点"。
环形图读不出精确位置。 圆周上没有可对齐的刻度，别指望在上面比较两个区域的距离。