2.3 全基因组鸟枪法拼装

无论采用 2.2 节所述的哪一类平台，测序仪的直接产物都是海量的短读段：一条条数百碱基的字符串，既不知道它们来自基因组的哪个位置，也不知道彼此的先后关系。全基因组鸟枪法(whole-genome shotgun, WGS)的策略是：索性不去追踪每条片段的来源——把整个基因组随机打碎、构建文库、大量测序，然后完全依靠计算，从读段之间的序列重叠把它们重新拼接回去。拼装 (assembly) 因此成为一个纯粹的组合学问题：给定数百万条读段，寻找能最好地解释它们的一组基因组路径。

这一路线在概念上早已被小基因组证实：1995 年，流感嗜血杆菌成为第一个以纯鸟枪法完成全基因组测序的自由生活生物（约 1.83 Mb、约 2.4 万条读段，Fleischmann et al., 1995）。1997 年，Weber 与 Myers 正式主张以 WGS 测定人类基因组，引发与公共计划分级克隆策略的著名论争（Weber & Myers, 1997）；2000 年果蝇基因组与 2001 年人类基因组序列的发表，最终证明两条路线在技术上均可行（Adams et al., 2000；Venter et al., 2001）。本节先建立鸟枪法的层级结构与计算框架（2.3.1），再给出其概率论基础（2.3.2），继而讨论重复序列的破坏作用与配对末端的解法（2.3.3）、拼装质量的度量与工具（2.3.4），最后讨论从支架到染色体的图谱锚定与两条策略的定量比较（2.3.5）。

2.3.1 鸟枪法思想与拼装层级

鸟枪法的思想可以概括为一句话：以冗余换连续。设想把一部只有孤本的书送入碎纸机：若只碎一份，碎片之间的衔接无从校验；若先把书复印多套再粉碎，则同一处文字会被多张碎片重复覆盖，相邻碎片只要共享一段足够长的相同文字，就可以对齐、粘连。基因组测序中的“复印”就是覆盖度冗余（2.3.2 节），“相同文字”就是读段末端的重叠 (overlap)，而“足够长”这一限定词至关重要——太短的重叠在随机背景下随处可见，会把毫不相干的读段错误地连在一起。

经典拼装程序因此普遍遵循重叠—布局—一致(overlap–layout–consensus, OLC)三阶段框架。第一阶段重叠检测：对所有读段做两两比对，寻找并评估候选重叠。直接比对的计算量随读段数平方增长，实践中先以共享 k-mer（k 个连续碱基的短词）预筛，仅对命中数超过阈值的读段对做精细比对；同时要求重叠长度、一致性比例都达标，以排除纯属巧合的短重叠与来自高拷贝重复的虚假重叠。第二阶段布局：把读段视为节点、合格的重叠视为边，构建重叠图(overlap graph)，在其中寻找一组协调的读段排列（铺叠路径，tiling path）。早期程序多采用贪心策略——自最长的重叠起不断延伸，遇到歧义即断开；更严格的实现把布局表述为在图上求解某种意义下的最优路径，在重复与噪声造成的分支处宁可结束当前重叠群，也不冒错拼之险。第三阶段一致序列：对铺叠路径上的每一列，按各读段碱基的质量值加权投票，得到一致性碱基与相应的一致性质量（见 2.3.4 节）。图 2.3-1 概括了这一流程与产物层级。

2.3.2 Lander–Waterman 覆盖度理论

随机打碎没有偏爱，也就没有保证：有的位置会被反复测到，有的位置可能始终无人问津。Lander 与 Waterman（1988）为此给出了经典模型。设基因组长度为 G，随机测得 N 条长为 L 的读段（L ≪ G），定义覆盖度(coverage / depth)为平均每碱基被测的次数 c = NL/G。由于读段起点近似构成泊松过程，任一基因组位置被 k 条读段覆盖的概率趋于泊松分布 P(k) = e^−cc^k/k!，于是该位置从未被覆盖的概率为：

式 (2.3-1) 的含义是：未覆盖概率随覆盖度指数衰减。这正是“以冗余换连续”的定量形式——覆盖度每增加 1 倍，漏测概率乘以 e⁻¹ ≈ 0.37。表 2.3-1 给出常用覆盖度下的数值。

P₀ 回答“有多少碱基漏掉”，工程上更直观的则是“留下多少个缺口”。未覆盖碱基的期望总数为 G·e^−c；若按缺口平均长度与读长 L 同一量级估算，缺口数目约为：

以两个尺度感受式 (2.3-2)：5 Mb 的细菌基因组，L = 800 bp、c = 8 时，G·e^−c/L ≈ 2.1，缺口只有个位数量级——高覆盖之下细菌基因组近乎封闭；而 3.1 Gb 的人类基因组，L = 700 bp、c = 5 时，缺口数约 3×10⁴ 量级（计入因子 c 约 1.5×10⁵），这就是人类草图需要数以万计 contig 的原因。还要注意缺口数量随 c 指数下降的同时，单个缺口的期望长度 L/c 也随 c 缩短，故提高覆盖度同时减少缺口的数目与长度，是 finishing（2.4 节）成本的第一决定因素。

下面的交互演示把式 (2.3-1) 变成可以亲手操作的实验：一条 1 000 bp 的环形迷你基因组（首尾相接以消除末端效应），随机撒上 150 bp 的读段。选择目标覆盖度或重新采样，观察覆盖深度的随机涨落、直方图与泊松理论曲线的吻合程度，以及实测零覆盖比例与 e^−c 的接近程度。

2.3.3 重复序列的破坏作用与配对末端

重叠合并有一条不可回避的前提：重叠必须唯一。若基因组里两段序列完全相同，来自两个位点的读段就无从分辨，重叠图在重复内部成为稠密的团块。当重复长度明显超过读长时，没有任何一条读段能“看穿”重复的两端，拼装器只能把两个拷贝坍缩(collapse)成一个——重复两侧的唯一序列因此出现多种同样讲得通的连接方式，任选其一都是潜在的错误拼装；更隐蔽时，坍缩直接把两段本应分开的序列压成一段，拷贝数被低估，基因组的该部分在拼装中缩短。重复序列因此被称为鸟枪法的“阿喀琉斯之踵”，也是式 (2.3-1) 无法刻画的限制：无论覆盖度多高，长于读段的重复都无法仅凭重叠穿越。

这一问题在人类基因组中尤为尖锐：重复序列约占基因组的一半，其中散在重复（转座子来源，如长约 6 kb 的 LINE-1 与约 300 bp 的 Alu 等 SINE）合计约 45% 量级，卫星 DNA 等串联重复另占数个百分点（Lander et al., 2001；重复家族的详细分类见 2.5 节）。两件事随之成为拼装成败的关键：其一是把重叠检测的阈值设得足够严格，宁可在重复边界断开、留下大量短 contig，也不冒险错拼——这是“重叠群数偏多”的自觉选择；其二是以配对末端(paired-end)测序为重复提供外部的方向与距离信息。

配对末端的机制值得稍作展开：构建基因组文库时，每个克隆的插入片段长度在一定区间内分布（如 2–3 kb），从插入片段两端各测一条读段，这对读段之间的实际距离虽有个体涨落，但服从已知的文库分布，且两条读段方向相向。于是每一对读段都携带一条独立于序列重叠的拓扑约束：它们必然位于同一染色体上，间距约等于文库均值，方向相对。落在重复两侧唯一区的一对读段，就能把两个 contig 的相对位置与朝向钉住——重复内部依旧读不出，但连接关系已被唯一确定。实践中同时构建多档插入长度的文库（数 kb 的质粒文库、约 10 kb 乃至 40–150 kb 的 BAC 末端序列），短的保证覆盖密度，长的跨越大段重复与缺失区域，构成分级的“跨越能力”。

配对末端还提供一道一致性检验。拼装完成后，逐对检查配对读段：间距与方向均符合文库分布的称为“正常对”；间距显著小于预期（被压缩）的一对往往指示该处发生坍缩——重复的两个拷贝被压近，两侧唯一序列随之被拉近；间距显著大于预期（被拉伸）则提示中间可能缺失一段序列或存在真实的长插入；方向异常则可能来自倒位或错误的连接。把全部配对关系在支架坐标上作图，异常簇出现的位置就是需要人工复核的候选位点，这是 Arachne 一类拼装器内置的复核步骤，也是 finishing 阶段（2.4 节）排定优先次序的依据（Batzoglou et al., 2002）。

2.3.4 拼装质量的度量与工具

给定一批读段，不同程序（或同一程序的不同参数）会给出不同的拼装，需要一个不依赖真值的指标刻画其连续性。行业通用标准是 N50。

Sanger 时代的拼装器谱系大致如下：Phrap（Green 实验室）以质量值加权的方式完成重叠检测与一致性投票，是 2.1 节所述 Phred–Phrap–Consed 流水线的中枢；Celera Assembler（Myers 等）为果蝇与人类的全基因组鸟枪拼装而设计，首次在大基因组尺度上系统使用配对末端约束复核布局（Myers et al., 2000）；Arachne（Batzoglou et al., 2002）以“重叠群合并—支架构建—配对一致性复核”三段式流程成为当时 WGS 拼装的代表性学术实现；Newbler 则是 454 平台的官方拼装器，针对焦磷酸测序的均聚物插入/缺失错误做了专门处理。进入短读长 NGS 世代后，主流框架由重叠图转向 de Bruijn 图（把读段拆成 k-mer、以图论方式重组），但其目标函数——在重复与误差干扰下寻找最可能的真实序列——与本节讨论并无二致（原书出版于 2009 年，正处于这一过渡期）。

工具之外，覆盖冗余本身即是纠错机制。若单条读段在某位的判读错误率约为 p，则深度 d 的一列中出现多数错误读段的概率随 d 迅速下降（量级为 p^d/2）；按质量值加权投票后，一致序列在该位的错误率可远低于任何单条读段，其置信度可折算成一致性质量值(consensus quality)——深度越大、各读段质量越高，一致性质量值越高。Phrap 即按列统计各碱基的加权似然给出该值。需要警惕的是：覆盖冗余只能消除随机的判读错误，无法消除系统性的偏差——同一克隆被反复测序时，其位点特异的错误会以“一致”的面目出现；同理，二倍体基因组中的杂合位点会呈现两条读段群体的稳定分歧，这是真实的多态而非错误，处理方式留待第 3 章讨论。

2.3.5 从支架到染色体：图谱锚定与策略之争

即便支架构建成功，拼装产物仍是一堆长短不一、方向不定的序列块——它们在染色体上的位置与次序尚属未知。把支架落到染色体坐标上，靠的是第 1 章 1.3 节所述的遗传图与物理图：以遗传图（STS、SSR 等标记，单位为厘摩尔根 cM）提供粗框架；以物理图（重叠 BAC 克隆群、STS 含量图谱、FISH 直接观察大片段在染色体上的位置）提供局部顺序与间距。支架内的序列可用 PCR 或杂交比对到标记上，每个被锚定的支架同时获得染色体归属、坐标区间与方向；未被锚定的支架则以 N50 等指标单独统计。人类基因组计划正是先以数万个 STS 标记建立“骨干”，再让序列沿骨干落位——图谱在此的意义与 1.3 节所述完全一致：为海量序列碎片提供独立于拼装计算的外部坐标。还应指出，染色体级的“完成”永远是相对的：着丝粒与近端粒的大段高度重复异染色质至今仍以缺口形式保留在人类参考基因组中。

由此可以定量地对全基因组鸟枪法与分级鸟枪法（先构建物理图、逐个 BAC 分别鸟枪测序再依图拼接，亦称 map-based 或 clone-by-clone 策略，见 1.4 节的建库与克隆框架）作一比较。成本与时限上，WGS 无需先行构建高质量的 BAC 重叠群图谱，文库制备与测序可高度并行：公共计划分级路线自 1990 年启动、2003 年宣布完成，总成本约 30 亿美元量级；Celera 的 WGS 路线在约一年内产出约 5 倍覆盖的自产读段并结合公共计划数据完成拼装，测序成本约数亿美元量级（量级估计）。代价则集中在重复处理：分级法把每个 BAC 内部的拼装限制在约 100–200 kb 的窗口内，长重复几乎总能被克隆边界截断，错拼风险小、缺口可定位到具体克隆逐一歼灭；WGS 则让拼装器直接面对全基因组所有重复的混淆，必须依赖多档配对文库与严格的复核规则。两条路线在人类基因组上的竞逐，最终以两篇论文同刊发表、公共分级草图与 Celera WGS 拼装互为印证告终（IHGSC, 2001；Venter et al., 2001）。

本节至此完成了从“读段”到“支架”的全部计算环节：覆盖度理论告诉我们需要测多深（式 2.3-1、2.3-2），重复序列与配对末端决定能拼得多准，N50 与配对一致性检验描述拼得多好，图谱锚定则把结果送回染色体坐标。剩下的缺口——无论物理缺口还是顺序缺口——如何逐一封闭，是 2.4 节 finishing 的主题；而把支架上的一长串 ATGC 变成基因与功能的目录，则从 2.5 节的基因组注释开始。