manhattan · Folio

Introduction

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

它回答的是"哪里值得跟进"，不回答"效应有多大"。本质上它就是一张散点图（位置 vs -log10(P)），只是横轴按染色体分了段——所以归在散点那一族，而不是单开一档。

读的时候按这个顺序：

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

Example Data

Manhattan plot 要的输入是：标记名 + 染色体 + 位置 + P 值，四列，一行一个标记。P 是原始 P 值，取不取 -log10 由 LOG10 参数决定，不用自己先转。

gwas_manhattan 是模拟的 GWAS 结果，12 条染色体 × 500 个标记。真实数据只要整理成 Marker / Chr / Pos / P 这四列，直接替换掉 gwas_data 就能跑。

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

# 公开地址，和数据集页上「下载 CSV」给的是同一个 —— 不写仓库相对路径：
# 那个目录不进仓库，读者 clone 下来也没有这个文件，这段代码就跑不了
gwas_data <- read.csv("https://assets.evanzhou.org/tessera/csv/gwas_manhattan.csv")

# 5e-8 是全基因组显著性的惯例（约等于 100 万次独立检验的 Bonferroni 校正），
# 1e-6 只是提示性水平。两条线含义差得很远，所以下面必须补图例
threshold_line <- c(5e-8, 1e-6)

# 只标 Top 8。标注超过十来个就会互相遮挡，反而看不见信号峰本身
top_hits <- gwas_data |> arrange(P) |> slice_head(n = 8)

Palettes

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

这五个颜色不是五个平级分类：背景染色体、重点位点和两级阈值承担的语义完全不同。所以这里按名称显式分配色板里的五个位置，而不是把色板循环到所有元素上——循环取色的结果是"第几个"决定颜色，而这里需要"是什么"决定颜色。

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

manhattan_cols <- setNames(
  get_palette("cancer_mosaic", type = "qualitative")[c(1, 3, 9, 7, 13)],
  c("odd", "even", "genome_wide", "top", "suggestive")
)

# CMplot 各参数收的是裸色值向量，名字得去掉
threshold_cols <- unname(manhattan_cols[c("genome_wide", "suggestive")])

# 白 → 深蓝的连续梯度，砍掉最白的一段：最低非零密度也要看得见
density_cols <- colorRampPalette(
  c("white", unname(manhattan_cols["odd"]))
)(171)[72:171]

Recipe

No.	Method	Input Data	Palettes
1	`CMplot`	`gwas_data`	`manhattan_cols`
2	`CMplot`	`gwas_data`	`manhattan_cols` / `density_cols`

1 · 直线型

有一条完整的横轴可以读坐标，所以标注位点、比较阈值线、写进正文配图都靠它。

#| 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
)

2 · 环形

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

#| 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

导出。 上面两段都写了 file.output = FALSE，因为要把图画到当前设备上；这个参数的默认值是 TRUE，不设就会在工作目录里落下一堆文件。但真要导出时反过来该用它——file = "pdf"、dpi = 300、file.name = "..." 交给 CMplot 自己处理，它会按图的内容挑合适的画布尺寸。手动 png() → CMplot() → dev.off() 也行，但尺寸得自己试。

Constraints

高点不等于效应大。 y 轴是 -log10(P)。样本量足够时，一个很小的效应也能给出极小的 P。
阈值线本身不带含义。 Bonferroni、FDR、suggestive 三者画出来都是一条虚线，图上分不出来，必须在图注或图例里说明用的是哪一条。
P 值必须大于 0。 有 0 或缺失时 -log10(P) 会出 Inf 或 NA，整张图的纵轴范围跟着歪。画之前先查一遍。
染色体和位置列错了不会报错。 图照样画得出来，只是横轴上的排序全乱——而这件事看图看不出来。
环形图读不出精确位置。 圆周上没有可对齐的刻度，别指望在上面比较两个区域的距离。

两个 recipe 怎么比

绝大多数情况 → recipe 1。要标注位点、要比较阈值线、要写进正文配图，直线型都更合适——它有一条完整的横轴可以读坐标。
recipe 2 只在两种场合更好：版面被限制成方形（比如放进多图拼版的一格），或者要传达的是"信号在全基因组上的分布密度"而不是"具体哪几个位点"。