曼哈顿图上已经出现了一个显著峰，候选基因也标出来了。结果到了投稿阶段，审稿人仍然追问：

真正与表型相关的是哪个变异？

它落在蛋白的什么位置？

为什么这个位置可能影响耐药？

问题往往不在于图画得不够漂亮，而在于证据链停在了“哪里有信号”，还没有继续回答“到底是什么变异在驱动信号”。

这也是全基因组曼哈顿图与候选基因局部曼哈顿图最核心的区别：前者负责发现，后者负责解释。

全基因组曼哈顿图回答——全基因组范围内，哪些区域出现了显著关联信号？

候选基因局部图与序列比对图回答——候选基因内部，信号集中在哪一段、对应哪些等位变异、效应方向是什么？

不过，先纠正一个常见误区：局部放大图并不是所有细菌GWAS论文都必须机械追加。当研究目标只是进行预测，或者全局分析已经直接定位到单一明确位点时，未必一定需要；但当你要比较同一基因内的多个信号、解释不同表型的位点差异，或者进一步提出耐药机制假设时，局部精细解析就非常重要。

下面这两张曼哈顿图来自CRyPTIC Consortium于2022年发表在PLOS Biology的一项大规模结核分枝杆菌GWAS研究。研究共纳入10,228株来自27个国家的结核分枝杆菌，并测定了13种抗菌药物的最低抑菌浓度（MIC）；经过不同药物对应的质量控制后，每项GWAS实际使用6,388～9,418株菌。

这里有一个非常关键的前提：图中的每个点，并不是传统意义上的一个SNP。

作者从组装序列中统计了11个氨基酸长度的寡肽和31bp的寡核苷酸，使用它们在不同菌株中的存在/缺失状态进行关联分析。研究共测试了约1,051万种寡肽和553万种寡核苷酸，主要采用GEMMA实现的线性混合模型（LMM），并利用遗传相关矩阵校正群体结构。

因此，后面看到的显著峰，严格来说代表的是某一类短序列存在/缺失模式与log₂MIC的关联。

一

第一张图：全局扫描，先回答“信号在哪里”

图1将13种药物分别绘制成13个曼哈顿图面板。它的任务是在全基因组范围内完成第一轮无偏扫描。

图 1  13 种药物 MIC 的全基因组寡肽 GWAS 曼哈顿图（原文 Figure 3）

图片来源：The CRyPTIC Consortium, PLOS Biology, 2022。doi:10.1371/journal.pbio.3001755.g003

这张图应该怎么读？

· 每个面板代表一种药物；横坐标为寡肽比对到H37Rv参考基因组后的物理位置（Mb）。

· 纵坐标为−log₁₀P，数值越高，说明该寡肽存在/缺失状态与MIC的统计关联越显著。

· 黑色虚线为该药物对应的Bonferroni 校正阈值。作者并不是简单按全部原始特征数校正，而是按独特的“phylo-pattern”数量计算，因此不同药物的阈值并不完全相同。

· 橙色表示β＞0，即该寡肽的存在与更高的log₂MIC相关；蓝色表示β＜0，即与更低的log₂MIC相关。颜色深浅还反映效应值大小。

· 图上方的基因名和字母标记对应每种药物最显著的候选基因或基因间区，具体排名需结合 Table 1阅读。

· 无法可靠比对到H37Rv的寡肽被放在每个面板最右侧，以浅灰色显示。

这张图真正说明了什么问题

1. 它证明全局关联分析能够重新识别已知耐药信号，同时发现尚未被耐药目录充分收录的候选区域。

例如，研究在多种药物中重新识别了rrs、embB、katG、gyrA、rpoB等经典耐药相关基因，说明分析框架能够捕捉已知生物学信号；与此同时，部分显著区域并不在作者采用的耐药目录中，为后续挖掘新机制提供了候选。

2. 它展示了不同药物的遗传信号分布并不相同。

有些药物的显著信号高度集中在少数经典靶基因附近，有些药物则在多个区域出现较分散的峰。原文中异烟肼、左氧氟沙星和莫西沙星显著基因或基因间区相对较少；而其他药物的候选区域更加复杂。

3. 它为后续局部解析确定优先级，但不能单独证明因果。

全局图能够告诉我们“哪一段基因组值得继续看”，却不能仅凭一个高峰就断言某个基因一定导致耐药。峰值可能受到连锁不平衡、谱系共现、交叉耐药、低频变异和重复区域等因素影响，还需要候选区域整理、具体变异解析以及后续验证。

二

优先看哪些基因或区域

全基因组曼哈顿图上同时出现大量字母和基因标记，如果只看图，很难判断哪些信号是已知耐药基因、哪些是新候选、哪些又落在重复区域。表1的作用，就是把全局信号整理成可阅读、可筛选的候选清单。

作者以每个基因或基因间区中最显著的寡肽为依据，对每种药物的候选区域进行排序，通常列出前20个；表中还标记了哪些基因已被既往耐药目录收录、哪些区域属于重复序列，并通过字母与图2中的峰一一对应。

表1. 13种抗结核药物GWAS中最显著的候选基因和基因间区域

表格来源：The CRyPTIC Consortium, PLOS Biology, 2022。doi:10.1371/journal.pbio.3001755.g003

三

第二张图：局部解析，

把“gyrB有信号”拆成具体变异

全基因组曼哈顿图帮助我们找到值得关注的候选基因，而候选基因局部图则进一步把这个基因放大，观察显著变异究竟集中在哪个位置、对应哪些氨基酸变化，以及这些变异与MIC升高还是降低有关。只看图1，我们最多能说gyrA/gyrB区域与两种氟喹诺酮类药物的MIC相关，却无法知道gyrB内部是不是同一个位置在起作用。

于是作者将gyrB单独放大，并进一步把显著寡肽比对到参考蛋白序列上，形成图2。

图 2  gyrB 中与左氧氟沙星和莫西沙星 MIC 显著相关的寡肽及其序列比对（原文 Figure 4）

图片来源：The CRyPTIC Consortium, PLOS Biology, 2022。doi:10.1371/journal.pbio.3001755.g004

A/B：候选基因局部曼哈顿图

· 横坐标不再是整个基因组的Mb坐标，而是寡肽在GyrB蛋白中的位置，因此可以观察基因内部的信号分布。

· 纵坐标仍为−log₁₀P，黑色虚线仍表示全基因组Bonferroni显著性阈值。

· 黑色寡肽表示其序列与gyrB注释的正确阅读框一致；灰色寡肽表示来自其他翻译阅读框。

· 图右侧单独排列的灰色寡肽，是先被nucmer分配到该区域、但未通过BLAST再次确认精确比对的序列。

C/D：显著寡肽序列比对图

· 每一条横向短条代表一个与MIC相关的11 aa寡肽，底部为H37Rv的参考氨基酸序列。

· 彩色方块标记发生变化的氨基酸残基；相同颜色对应相同氨基酸。

· 背景深灰表示β＞0，与更高MIC相关；背景浅灰表示β＜0，与更低MIC 相关。

序列按显著性排列，越靠上关联越显著；虚线以上为达到全基因组显著的寡肽。

这张图真正说明了什么问题？

1. 左氧氟沙星（LEV）和莫西沙星（MXF）虽然都指向gyrB，但显著信号集中在不同位置。

两个局部图中都可以看到相邻的两组峰，但对于每一种药物，只有其中一组超过显著性阈值，而且两种药物对应的显著峰并不相同：左氧氟沙星的主要信号集中在约461 aa，莫西沙星的主要信号集中在约501 aa。

2. 局部序列比对把“基因水平关联”进一步拆解为“具体氨基酸等位变化”。

左氧氟沙星的显著寡肽主要捕捉457和461位点，其中标记D461N的寡肽与更高MIC相关，β＝2.46（log₂MIC尺度），在7300株菌中出现于15株；莫西沙星的显著寡肽主要捕捉499和501位点，其中标记E501D的寡肽与更高MIC相关，β＝1.86，在6388株菌中出现于23株。

与之相对，携带参考等位氨基酸的高频寡肽通常与较低MIC相关。这也解释了为什么同一位置附近会同时出现β＞0和β＜0的寡肽：它们往往分别代表突变等位序列与参考等位序列。

3. 两种药物的局部差异提示了药物特异性的变异效应。

GyrB的461和501位点均位于蛋白与氟喹诺酮结合界面附近，因此这些关联具有一定结构生物学合理性。作者据此认为，gyrB尤其是E501D值得纳入未来莫西沙星耐药目录的评估。

需要注意的是，局部图可以帮助我们缩小范围、提出机制假设，但不能仅凭一张图就证明某个突变一定导致耐药，后续仍需要结合功能注释、已有文献或实验验证。

四

两张图合起来，完整证据链是什么？

全基因组曼哈顿图把研究从“数百万个特征”缩小到“少数候选区域”；

局部图和序列比对再把“一个候选基因”缩小到“少数具体变异”；

真正的机制结论，还要靠后续验证完成最后一步。

结语：曼哈顿图不是终点，

而是候选机制的入口

很多细菌GWAS项目做到显著峰和候选基因列表就停止了，但对读者和审稿人来说，他们真正关心的通常是：这个峰由什么变异构成？不同变异的效应方向是否一致？信号能否落到具有生物学意义的蛋白位置？

这篇结核分枝杆菌研究给出了一个很清晰的写作范式：先用全基因组曼哈顿图完成无偏发现，再用候选清单确定优先级，最后通过局部关联图和序列比对把信号落到具体氨基酸。

所以，真正完整的细菌GWAS结果，不是画出一个峰就结束，而是要把峰逐步转化为可解释、可验证的候选机制。

如果你正处于：

· 已经积累了一批同一物种的菌株样本，也有耐药性、毒力、宿主来源或临床结局等表型数据，想进一步开展细菌GWAS；

· 手中已有细菌全基因组测序数据，却不清楚应该选择SNP-GWAS还是PAN-GWAS；

· 不确定现有样本量和表型数据是否适合做GWAS，也不知道连续型、二分类表型应该如何整理。

· 担心菌株亲缘关系和群体结构造成假阳性，不清楚应该如何建树、构建亲缘矩阵和选择关联模型；

· 已经完成细菌GWAS，却不知道如何解释曼哈顿图、Q-Q图和显著位点表；

那么，开展细菌GWAS的第一步，并不是直接运行某个软件，而是先判断样本、基因组数据和表型数据是否匹配，再根据研究目标选择合适的变异类型、群体结构校正方法和关联模型。

目前，对于需要个性化分析、复杂数据处理和结果深度解读的项目，唯那生物可以提供细菌 GWAS 一对一分析与项目服务，根据你的菌株数量、物种、表型类型和研究目标，制定具体分析方案。

与此同时，我们也正在将其中标准化程度较高的 GWAS 分析环节逐步整理为在线云分析流程。

参考文献

[1] The CRyPTIC Consortium. Genome-wide association studies of global Mycobacterium tuberculosis resistance to 13 antimicrobials in 10,228 genomes identify new resistance mechanisms. PLoS Biology. 2022;20(8):e3001755. doi:10.1371/journal.pbio.3001755

支持单位：