基因组组装

引言

基因组组装（genome assembly）是把一堆无序的测序读段拼回完整基因组序列的过程。它是所有「没有参考基因组」研究的第一步：新物种测序、微生物完成图、肿瘤重排结构解析，都要先过这一关。与重测序（有参考可比对）不同，组装面对的是「重建未知序列」这个组合优化问题——没有标准答案可对照，错误往往要等到下游分析才暴露。

真正的难点是重复序列。人类基因组约 50% 由重复元件构成，转座子、片段重复（segmental duplication）、着丝粒卫星序列遍地都是。短读（150 bp）根本无法跨越比它长的重复，图结构会在重复处「坍塌」成一条路径，导致组装出的序列在重复区被压缩或错连。这是二代测序时代组装碎裂的根本原因，也是长读测序改变游戏规则的逻辑起点。

三代长读（PacBio HiFi 读长 15-25 kb、准确率 99.9%；ONT 超长读可达 100 kb 以上但原始错误率 5-10%）让读段能整段跨越重复，使「端粒到端粒」（T2T）的完整组装成为现实。但代价是成本、算力与新的算法权衡：HiFi 准但贵且需要高分子量 DNA，ONT 便宜读长长但需要纠错或高覆盖共识。

本文按「问题本质 → 图算法 → k-mer → contig → scaffold → 长读工具选型 → Hi-C 挂载 → 质量评估 → T2T 与泛基因组」的顺序展开。命令基于 hifiasm 0.19、Flye 2.9、QUAST 5.2 与 BUSCO 5.6。读完你应该能独立设计一个组装项目，并判断一份组装结果「能不能用、用在哪」。

目录

1. 组装问题的本质与三代测序格局
2. 从读段到图：OLC 与 de Bruijn 的取舍
3. k-mer 选择、纠错与基因组特征估计
4. contig 构建与单元型分辨
5. scaffold 与 gap 填补
6. 长读组装工具选型与参数
7. Hi-C 挂载到染色体水平
8. 组装质量评估：N50、QUAST、BUSCO
9. T2T 与泛基因组概览

1. 组装问题的本质与三代测序格局

组装的形式化定义是：给定读段集合与覆盖度，输出一组 contig（连续无 gap 的序列）与 scaffold（由 contig 加 gap 拼成）。理论上它是「最短公共超串」（SCS）问题的变体，属 NP-hard；实践中靠图论启发式求解，因此「最优解」不可验证，只能靠质量指标间接评估。

三个根本困难贯穿始终：

• 重复：读段太短就跨不过去，图在重复处出现分叉或坍塌；
• 杂合：二倍体样本的两套等位基因会分叉成两个「气泡」，处理不当会把基因组撑大一倍或把等位差异当成错误抹掉；
• 错误：测序错误产生假 k-mer，污染图结构，需要纠错或靠覆盖度过滤。

覆盖度是首要设计参数。经验值：HiFi 建议 20-40x（hifiasm 在 30x 左右已能给出高质量单倍型组装）；ONT 需要 30-60x，若要用超长读打通着丝粒则需额外 10-20x 的 100 kb+ 读段；短读仅用于 polishing 或混合纠错时，20-30x 即可。

常见物种的组装难度与推荐策略
  细菌（5 Mbp，低重复）      ：HiFi 30x 或 ONT 50x，单次运行即可完成图
  真菌（30-100 Mbp）         ：HiFi 40x + Hi-C，染色体级
  植物（500 Mbp-10 Gbp，高重复）：HiFi 40x + ONT 超长 + Hi-C，分染色体
  哺乳动物（3 Gbp，50% 重复） ：HiFi 30x + trio/Hi-C 相位，T2T 需 Verkko

覆盖度与成本的换算

覆盖度不是「越高越好」的单调函数，它存在明显的收益拐点。HiFi 从 20x 提到 30x，contig N50 通常有数倍提升；从 30x 提到 50x，提升幅度迅速收敛，但成本线性上升。这是因为组装的瓶颈已从「覆盖度不足导致断裂」变成「重复结构无法用当前读长解开」——此时该增加的是读长而非覆盖度。

反过来，覆盖度不足会以「静默」的方式伤害组装：图在低覆盖区断裂，组装器把低覆盖节点当作错误剪除，结果是一份看起来干净、实际丢失了部分基因组的组装。这类错误只有 BUSCO 完整性检查能暴露，N50 完全看不出来。

成本换算的经验数字：PacBio Revio 单张 SMRT cell 约产出 25-30 Gbp HiFi（人类约 8-10x），3 Gbp 基因组做 30x 约需 3-4 张；ONT PromethION 单张 flow cell 约 50-100 Gbp，但产出波动大。一个常被低估的成本项是高分子量 DNA 提取——HiFi 与超长 ONT 都要求 DNA 长度大于 50 kb，常规柱式提取的片段化 DNA 会让长读比例骤降，前期投入不足会直接毁掉项目。

2. 从读段到图：OLC 与 de Bruijn 的取舍

组装算法分两大流派，选错流派会浪费大量算力：

OLC（Overlap-Layout-Consensus） 把读段当作节点，两两计算重叠（overlap），再用重叠关系构图、找路径、做共识。它的优势是「保留读段级信息」，能处理较长读段和低覆盖；劣势是重叠计算复杂度接近 O(n²)，读段数上百万时爆炸——所以现代实现都用 minimizer 索引把重叠搜索降到近线性。

de Bruijn 图把 k-mer 当作节点（或把 (k-1)-mer 当作节点、k-mer 当作边），用「相邻 k-mer 共享 k-1 个碱基」连边。它的优势是内存与覆盖度无关（只与 k-mer 总数有关），适合短读高覆盖；劣势是丢失读段级信息，且错误 k-mer 会制造大量假分叉，必须先纠错。

现代长读组装器多数是 OLC 的变体（string graph）：hifiasm、Flye、Canu 都属于这一类。de Bruijn 仍活跃在短读组装（SPAdes 做细菌完成图）与混合纠错（用短读给长读纠错）场景。选型口诀：读段能跨越重复就用 OLC，只有短读才用 de Bruijn。

string graph 为什么能省内存

OLC 的朴素实现要存储所有读段对的重叠，空间复杂度接近 O(n²)。string graph 的改进在于「传递性削减」（transitive reduction）：若读段 A 与 B 重叠、B 与 C 重叠，且 A 与 C 的重叠被 B 完全蕴含，则 A-C 这条边冗余，可删。削减后边数从「接近读段数的平方」降到「接近读段数」，这是长读组装能在普通服务器上处理 3 Gbp 基因组的根本原因。

但传递性削减有一个隐含假设：重叠关系可传递。重复序列恰恰破坏这个假设——A 与 C 的重叠可能通过某个重复拷贝建立，而非真实相邻关系。所以 string graph 在重复区仍需「按读长与覆盖度区分重复拷贝」，这正是组装器复杂度最高的部分，也是不同工具结果差异的主要来源。

3. k-mer 选择、纠错与基因组特征估计

k-mer 是组装前最重要的「免费诊断」。取一段读段的所有长度为 k 的子串，统计每个 k-mer 出现的次数，得到 k-mer 谱（k-mer spectrum）：

from collections import Counter

def kmer_spectrum(seq, k=21):
    """统计单条序列的 k-mer 频次，用于估计基因组特征"""
    return Counter(seq[i:i+k] for i in range(len(seq) - k + 1))

# 真实数据应扫描全部读段，并对双端读段的 k-mer 计数除以 2（每个位点被两条读段覆盖）
counts = kmer_spectrum("ACGTACGTACGTACGTACGTACGT", k=5)
print(counts.most_common(3))

谱形能直接读出三个参数：

• 主峰位置 ≈ 测序覆盖度（如峰值在 30，说明约 30x）；
• 左侧小峰 ≈ 杂合度（杂合位点把读段分成两套，各占覆盖度一半，形成 0.5x 处的小峰）；
• 主峰后的长尾 ≈ 重复率（重复序列贡献了 2x、3x 的倍数峰）。

GenomeScope2 把这个谱拟合出基因组大小、杂合度、重复率三个数字，是项目启动前的必跑步骤：

谱形异常的诊断

# 1. 生成 k-mer 直方图（k=21 是常用起点）
jellyfish count -m 21 -s 100M -t 32 -C -o mer_counts.jf reads.fq.gz
jellyfish histo mer_counts.jf > reads.histo

# 2. 用 GenomeScope2 拟合（R 脚本）
Rscript genomescope2.R reads.histo 21 150 output_dir

k 的选择是经典权衡：k 太小（< 17）会让重复区坍塌，图被压缩；k 太大（> 主峰读长的 2/3）会让覆盖不足的位点断裂，且内存与错误率上升。经验起点是 k ≈ 读长的 1/2 到 2/3，短读 150 bp 常用 k = 21-31，长读纠错常用 k = 19-25。实践中要「试几个 k 看 N50 与总长是否稳定」——若 k 变化导致总长剧烈波动，说明重复结构复杂，需要更长读段而非调 k。

纠错策略取决于平台：HiFi 读段本身准确率已 > 99.9%，不需要纠错，直接组装；ONT 原始读段需要纠错（Canu 的 -correct、NextDenovo 的 correct），或走「不纠错但靠高覆盖共识」的路线（Flye 的 --nano-raw）；短读纠错长读（如 Ratatosk）在混合数据下效果最好但流程复杂。

4. contig 构建与单元型分辨

从图到 contig 的核心是「找路径并解分叉」。图上的分叉有两类：

• 气泡（bubble）：两条平行路径，通常来自杂合位点或测序错误；
• 尖端（tip）：从主路径伸出又断掉的短枝，几乎都是错误。

处理策略是「按覆盖度与长度过滤 + 解气泡」。低覆盖节点（覆盖度远低于主峰）判为错误并剪除；长度相近、覆盖度各半的气泡判为杂合，此时面临关键抉择：是保留两套单倍型（phased assembly）还是合并成一套（collapsed/primary）。

hifiasm 的做法值得学习：它默认输出两套单倍型（--primary 之外还有 .hap1/.hap2），若提供 trio 数据（父母样本）或用 Hi-C 数据，能通过「读段来源相位」把气泡正确分配到两条单倍型上，得到真正的二倍体组装：

# HiFi 单样本组装（不提供 trio/Hi-C 时输出 primary + 未定相 alt）
hifiasm -o asm -t 32 hifi.fq.gz

# 提供 Hi-C 数据做相位（推荐，能显著提高单倍型连续性）
hifiasm -o asm -t 32 --h1 hic_R1.fq.gz --h2 hic_R2.fq.gz hifi.fq.gz

# 提供 trio 数据做相位（最准，需要父母样本）
hifiasm -o asm -t 32 -1 pat.yak -2 mat.yak hifi.fq.gz

-1/-2 需要先用 yak 从父母短读构建 k-mer 计数库（yak count -k31 -b37 -o pat.yak pat_R1.fq.gz pat_R2.fq.gz）。trio 相位是当前最可靠的单元型分辨手段，成本上要多测两个样本，但对杂合度高的物种（如很多植物、鱼类）几乎是必需品——否则两套单倍型会在组装里混成一团。

单元型组装还是合并组装

这是组装项目最早要做的决策之一，且很难事后反悔：

实践中 hifiasm 的默认输出（primary + 两个 hap）覆盖了大多数需求：做参考比对时用 primary，做等位特异分析时用 hap1/hap2。但要注意 primary 不是「真正的单倍型」，它是「把两套等位合并后的一条嵌合序列」——在杂合位点处，primary 的碱基选择是任意的，这对下游的等位特异分析是致命的，必须改用 hap 输出。

5. scaffold 与 gap 填补

contig 之上是 scaffold：用「能连接两个 contig 的证据」把它们串起来，中间用 N 占位。证据来源有三种：

1. 配对读段 / mate pair 文库：同一条长插入文库的两端比对到不同 contig，说明它们相邻；插入长度分布给出 gap 长度估计；
2. 长读：一条长读横跨两个 contig，直接给出顺序与 gap 大小（最可靠）；
3. Hi-C 接触：用于染色体级挂载（见第 7 节）。

# LR_Gapcloser：用长读填补 scaffold 内部的 N gap
LR_Gapcloser.sh -i scaffolds.fa -l ont.fq.gz -s fill -t 32 -o filled.fa

# hifiasm 内部已用长读做 gap 填补，通常无需额外步骤

gap 填补（gap filling）是 scaffold 的收尾：用读段填掉 N 区，把 scaffold 还原成完整序列。评估时要注意报告N 的数量与总 gap 长度——一个 N50 很高但 gap 很多的组装，实际可用性可能不如 gap 少而略碎的版本，因为 gap 区的基因无法分析。

一个判断「gap 该不该填」的实用标准是看它落在哪里：若 gap 落在基因密集区（可从近缘物种的注释推断），填补价值高；若落在已知的重复区或着丝粒，填补成功率低且容易引入错误序列，不如保留 N 并如实报告。把 gap 的位置、长度、是否可填补写进组装报告，比一味追求「零 gap」更有工程价值。

一个容易踩的坑是「scaffold 的错误连接」：mate pair 的假连接会把两个不相邻的 contig 硬串起来，形成「染色体级但序列错误」的产物。这类错误在 N50 指标上表现为「好看」，却会让下游的共线性分析、SV 检测全面失真。因此现代流程更倾向「以长读/Hi-C 为主证据，谨慎使用 mate pair」。

6. 长读组装工具选型与参数

选型逻辑很直接：有 HiFi 就用 hifiasm；只有 ONT 且追求速度用 NextDenovo/Flye；目标 T2T 且有三代混合数据用 Verkko。

资源画像（以 3 Gbp 哺乳动物基因组为参照）：

这类任务对内存带宽与并行文件系统吞吐敏感，通常放在 HPC 集群上跑，参见 HPC 与并行计算 。工程上有一个常被忽略的点：组装过程会产生大量中间文件（string graph 的 dump、纠错的读段副本），峰值磁盘占用可能是输入数据的 5-10 倍。把工作目录放在本地 NVMe 而非共享并行文件系统上，通常能缩短 20%-40% 的总时间。

一个实用参数细节：hifiasm 的 -t 是线程数，--n-hap 指定倍性（默认 2，单倍体/细菌设 1），-l 用于调整主次组装的分离阈值（值越大越倾向把等位差异归为杂合）。对于高度杂合（杂合度 > 1.5%）的样本，-l 3 到 -l 5 能明显改善单元型分离。

7. Hi-C 挂载到染色体水平

Hi-C 通过甲醛交联把空间上邻近的 DNA 固定，再测序接触对。它的关键信号是：同一条染色体内部的接触频率远高于染色体之间，且接触频率随基因组距离衰减。这两个性质分别用于「分组」（哪些 contig 属于同一染色体）和「定向」（contig 的方向与顺序）：

# YaHS：快速 Hi-C 挂载（Juicer 的轻量替代）
yahs contigs.fa hic_R1.fq.gz hic_R2.fq.gz -o yahs_out

# 或用 Juicer + 3D-DNA 流程（老牌，生态完整）
juicer.sh -g genome_id -s MboI -z ref.fa -y site_file -p chrom.sizes \
  -D /path/to/juicer ...

挂载质量靠接触矩阵可视化判断：好的挂载在矩阵上呈现清晰的「方块」（每条染色体一个方块），若出现横跨两个方块的条纹，说明存在错误连接（misjoin）。工具如 Juicebox 支持手工修正（在矩阵上拖拽调整 contig 顺序与方向），是发表级组装的常规步骤。

Hi-C 挂载的局限：它给出的是「相对顺序」，不保证绝对正确——小 contig 的定向容易出错，近端着丝粒区域的接触信号弱。因此 Hi-C 挂载后应回到「长读是否支持这个连接」做交叉验证，尤其对小于 1 Mb 的 contig。

Hi-C 数据的质量同样需要先评估：看「有效接触对比例」（通常 30%-60%）、「同一染色体接触占比」（应显著高于染色体间），以及接触矩阵是否随距离单调衰减。一份建库失败的 Hi-C 数据（交联不充分、酶切不完全）会给出随机噪声矩阵，用它挂载只会把 contig 排成错误的顺序，且错误很难通过统计指标发现——必须肉眼看矩阵。

8. 组装质量评估：N50、QUAST、BUSCO

组装没有「正确答案」可比，所以质量评估是判断可用性的唯一依据。指标分三层：

连续性指标：N50 是「把 contig 按长度降序累加，累计长度达到总长 50% 时的那个 contig 长度」。它比「平均长度」更能反映长片段情况，但有严重误导性——错误连接会人为抬高 N50，一个把多条染色体误串成一条的组装 N50 会非常漂亮。所以 N50 必须与 L50（达到 50% 所需的 contig 条数）、NG50（以基因组大小为分母的版本，需已知基因组大小）联合看。

完整性指标：BUSCO 用「单拷贝直系同源基因集」（如哺乳动物用 mammalia_odb10 的 9225 个基因）扫描组装，报告 Complete（完整）/ Fragmented（片段）/ Missing（缺失）比例。它回答的是「组装丢了多少保守基因」，是跨物种可比的硬指标：

quast.py -t 32 -o quast_out -m 1000 --large asm.fa          # 连续性 + 统计报告
busco -i asm.fa -l mammalia_odb10 -m genome -c 32 -o busco_out   # 完整性
merqury.sh reads.meryl asm.fa merqury_out                    # k-mer 一致性 → QV

碱基级准确度：k-mer 一致性（Merqury）用「组装中出现的 k-mer 是否都能在读段中找到」估计 QV（质量值）；LAQ（Long-read Alignment Quality）用长读比对到组装上统计错配与 indel，是当前公认最贴近「真实碱基错误率」的指标。HiFi 组装的 QV 通常 > 50（错误率 < 1/100000），ONT 组装经 polishing 后约 QV 40-50。

需要注意 BUSCO 的适用边界：它依赖「直系同源基因集」与「单拷贝」假设，因此对多倍体（基因集里的单拷贝基因实际有多份）、近缘杂交物种、细胞器基因组都会给出偏低或误导的分数。跨物种比较时要用同一套 BUSCO 数据集与版本，否则数字不可比。另一个细节是 BUSCO 报告里的 Duplicated 比例——若高达 20% 以上，通常意味着「两套单倍型被当成重复」或组装里混入了污染，是一个很有价值的诊断信号。

指标之间的冲突

这些指标并非总是一致，冲突时往往正是问题所在：

• N50 高但 BUSCO 低：错误连接把不该合并的区域串起来，或把重复区压缩了。查 Hi-C 接触矩阵与近缘物种共线性。
• BUSCO 高但 QV 低：序列结构对了，但碱基有大量小错误（通常是 ONT 组装未充分 polishing）。补一轮 polishing 即可。
• 总长明显小于预期：重复区被坍塌。用「读段回贴覆盖度」检查——若某些区域覆盖度是 2x 而非 1x，说明那里本应是两条重复序列被压成了一条。
• 总长明显大于预期：杂合被拆成两套（未合并）或污染混入。先查 k-mer 谱的杂合峰。

一个务实的做法是「把组装当假设检验」：用组装结果回贴读段，检查覆盖度是否均匀（理想是 1x，性染色体与细胞器除外）。覆盖度异常的区域，就是组装可疑的区域。

9. T2T 与泛基因组概览

T2T-CHM13（2022）是组装史上的里程碑：它首次给出人类基因组的完整序列，补全了 GRCh38 遗留的约 8% 缺口——着丝粒卫星阵列、rDNA 簇、片段重复区。实现路径是「HiFi + ONT 超长读 + trio 相位 + Verkko 组装 + 手工修正」，把人类基因组的 N 从数万个降到只剩几个。

T2T 的价值不只是「完整」：着丝粒、rDNA 这些区域此前从未被系统研究，现在成为新的发现空间（着丝粒变异与癌症、衰老的关联）。对下游分析的影响也很大——用 T2T 作参考，重复区的比对与 SV 检测都更准确，代价是坐标与 GRCh38 不通用（需 liftOver 且重复区不可靠）。

泛基因组（pangenome） 是对「单一线性参考」的反思。人类参考基因组本质上是「少数个体的嵌合体」，无法代表全球人群的结构变异。HPRC（Human Pangenome Reference Consortium）用 47 个样本构建了图基因组（pangenome graph），用节点与边的图结构表示「哪些序列是共享的、哪些是特定单倍型独有的」：

# 构建图基因组（minigraph-cactus 流程，简化示意）
minigraph -cxggs -t 32 ref.fa sample1.fa sample2.fa > graph.gfa
vg convert -g graph.gfa -x > graph.xg      # 转成 vg 格式供比对
vg giraffe -Z graph.gfa -f reads.fq.gz > aln.gam   # 把读段比对到图上

图基因组的取舍：它更准确地表示结构变异，但工具生态仍在成熟中。比对（vg giraffe、GraphAligner）、变异检测（vg call）、可视化都有一套新工具链，学习成本高；下游的注释工具对图格式支持有限。当前实践是「线性参考做常规分析，图基因组做 SV 富集区域的精细分析」，两者并存而非替代。

对组装项目而言，泛基因组带来两个具体影响。其一，「参考」不再是唯一真理：同一物种不同个体的序列差异（尤其结构变异）会被图结构显式表达，比对时读段可以走到「自己的那条路径」上，减少参考偏倚（reference bias）导致的假变异。其二，组装目标从「一条序列」变成「一个图」：HPRC 的做法是先各自组装 47 个样本的单元型，再用 minigraph-cactus 对齐合并成图，组装质量直接决定图的质量。

权衡取舍

常见坑清单

1. 只看 N50 判断好坏：错误连接会抬高 N50；必须同时看 L50、BUSCO 完整性与 QV。
2. k-mer 估计被污染样本误导：样本带细菌污染会让基因组大小虚高；先看谱形是否有异常峰。
3. 忽略杂合度直接组装：高杂合物种会得到「总长为基因组 1.5 倍」的假组装；用 GenomeScope 先估杂合度。
4. 用短读组装高重复基因组：重复区必然坍塌；改用长读或混合策略。
5. polishing 用错参考：用近缘物种读段 polishing 会引入偏差；用同物种数据。
6. Hi-C 挂载后不做验证：小 contig 定向易错；用长读支持度交叉检查。
7. 把 scaffold 的 N 当成序列：N 区无法做基因分析；报告时必须单列 gap 统计。
8. 基因组大小预估错误导致覆盖度不足：大基因组按小基因组量测序，组装必然碎裂；先跑 k-mer 估计。
9. 忽视细胞器基因组：叶绿体/线粒体序列常被误拼进核基因组；组装前先分离或事后过滤。
10. 不同工具结果直接混合：hifiasm 与 Flye 的 contig 命名与坐标系不同；混合前统一命名与去冗余。

小结

基因组组装是「无参考条件下的序列重建」，它的全部难点都源于重复、杂合与错误这三件事。图算法（OLC 与 de Bruijn）决定了工具的基本能力边界，k-mer 谱决定了你对基因组的先验认知，相位与 Hi-C 决定了能否把两套单倍型和染色体结构还原出来，而质量评估决定了这份组装能支撑到什么级别的结论。

工程上最该建立的认知是「指标会骗人」。N50 可以被错误连接抬高，总长可以被污染撑大，唯一的办法是多指标交叉：连续性看 N50/L50，完整性看 BUSCO，准确度看 QV 与 LAQ，结构正确性看 Hi-C 接触矩阵与共线性。任何一个指标单独看都可能得出错误结论。

下一步可以看 测序原理与数据产出 理解长读与短读的技术差异，或 变异检测与 VCF 处理 了解组装好的参考如何用于下游变异分析。如果你的组装「总长是预期的两倍」，先查杂合度与相位设置；如果「N50 高但 BUSCO 低」，多半是错误连接把不该合并的区域串在了一起。

常见问题

Q：细菌基因组用短读能组装完成图吗？
可以。细菌基因组小（通常 2-7 Mbp）、重复少，SPAdes 用 100x 短读往往能给出接近完成图的结果。只有含大量 rRNA 操纵子或插入序列的菌株才需要长读。

Q：组装好的序列可以直接当参考用于重测序分析吗？
可以，但要先做评估：BUSCO 完整性、QV、以及「与近缘参考的共线性」。低质量组装当参考会系统性地引入假变异——这是自建参考最常见的隐性错误来源。参见 SAM/BAM 与 samtools 实战 了解比对后的验证手段。