第二季 · 第 54 章

把破碎的生命之书重新拼起来

《算法的故事》第二季 · 第五篇：破解生命的密码——从桑格的DNA测序、全基因组鸟枪法，到德布鲁因图与欧拉路径组装算法

1977年，英国剑桥。一群科学家终于完成了一件长期令人难以想象的工作。他们测定了一种噬菌体的完整DNA基因组序列。这种噬菌体叫作：φX174。它是一种感染细菌的病毒。它的遗传物质是单链环状DNA。完整序列包含大约五千四百个核苷酸。对于今天的测序技术来说，这是一段很短的基因组。但在1977年，它代表着一项巨大的成就。因为科学家第一次能够将一个DNA生物体的完整基因组，以确定的字符顺序写下来。这项研究的主要人物之一，正是我们在第51章认识的：弗雷德里克·桑格（Frederick Sanger）。1977年2月24日，《Nature》发表：Nucleotide Sequence of Bacteriophage φX174 DNA《噬菌体φX174 DNA的核苷酸序列》。论文由桑格及八位合作者共同署名。他们报告了约5375个核苷酸的DNA序列，并揭示了噬菌体基因组中某些基因在不同阅读框里发生重叠的现象。 ^{[Nature][PubMed]}

桑格时期的DNA测序凝胶与早期序列记录。那个时代的研究者，需要把分散的实验读数整理成连续的核苷酸顺序。但是，完成基因组序列并不意味着实验人员拥有一台能够一次读出整条DNA的机器。他们获得的，是许多较短的测序结果。每一条结果都像从一本书中撕下的一段文字。其中一些片段彼此重叠。科学家必须根据这些重叠部分，重新判断它们原来的先后关系。例如，一段DNA原本是：

ATGCGAT

但我们只能得到：

ATGC
TGCG
GCGA
CGAT

而且，这些片段可能以完全打乱的顺序出现。问题是：如果不知道每段DNA原来位于什么位置，我们还能把它们拼回去吗？当基因组只有几千个字符时，研究人员还能够结合实验图谱、重叠分析和人工核查逐步完成任务。但随着测序技术发展，科学家希望读取的DNA越来越长。几万个碱基。几百万个碱基。甚至数十亿个碱基。每一条短读段都只提供局部信息。而组装程序必须从成千上万甚至更多的局部片段中，重建原来的顺序。到1995年，科学家利用全基因组随机测序与组装，获得了一种自由生活细菌的完整基因组。这种细菌叫作流感嗜血杆菌（Haemophilus influenzae Rd）。其基因组长度为：

\boxed{1\,830\,137\text{ bp}}

1995年7月28日，《Science》发表了这项成果。论文标题是：Whole-Genome Random Sequencing and Assembly of Haemophilus influenzae Rd《流感嗜血杆菌Rd的全基因组随机测序与组装》。 ^[PubMed]在这项工作中，科学家不必事先将整个染色体精确划分成一个个已经确定位置的大片段，再按照这些片段的位置依次测序。他们可以先随机获取大量DNA片段，再依靠计算与进一步实验恢复其对应关系。这便是：Whole-Genome Shotgun Sequencing。全基因组鸟枪法测序。但成功组装一个基因组，并没有消除所有困难。因为生命的DNA并不是一本没有重复词语的书。很多片段会重复出现。有的序列重复两次。有的重复成百上千次。两段DNA可能具有完全相同的短词，却属于基因组中不同的位置。这时，组装程序不只是要解决字符重叠。它还必须面对一种图论上的歧义：经过同一个重复片段之后，究竟应该接向哪一条后续路径？2001年8月14日，《Proceedings of the National Academy of Sciences》发表了一篇重要论文。An Eulerian Path Approach to DNA Fragment Assembly《一种基于欧拉路径的DNA片段组装方法》。作者是：帕维尔·佩夫兹纳（Pavel A. Pevzner）。唐海旭（Haixu Tang）。迈克尔·沃特曼（Michael S. Waterman）。他们提出了一套名为EULER的组装算法思想。其中重要的转变是：与其将每条长短不一的读段都直接当作图里的节点，不如将它们拆解成长度固定、具有前后缀关系的短词。然后，把组装问题转化为寻找经过相应图边的路径。这项研究与一种经典图论问题紧密相连。欧拉路径（Eulerian Path）。而用于表示短词重叠的图，则与荷兰数学家尼古拉斯·德布鲁因（Nicolaas de Bruijn）的名字联系在一起。2001年的论文第98卷第17期，第9748—9753页，明确介绍了这种图论组装路线。 ^{[PubMed][PMC]}

在这一章，我们将看到：一段DNA怎样被打碎成k-mer。一批短词怎样变成一张有向图。为什么节点代表前后缀，边却代表完整的k-mer。为什么组装问题有时变成寻找一条欧拉路径，而不是寻找经过每个节点一次的哈密顿路径。以及，当重复序列和测序错误出现时，为什么一张图既能帮助我们恢复信息，也会暴露我们根本无法仅凭现有数据解决的歧义。这个故事从生物化学实验开始。却最终回到了18世纪的一道数学难题。

第一幕：从1977年到1995年，科学家为什么愿意先打碎DNA，再尝试拼起来？

如果我们已经拥有一条完整DNA，为什么还要把它切成碎片？答案其实非常简单。在真实测序实验中，我们往往并不知道完整序列是什么。那是我们最终想获得的答案。但测序方法的单次读取能力有限。所以，实验会生成许多较短的序列读数。这些读数叫作：Reads（读段）。

假设原始DNA是一本完整的书

原始序列：

ATGCGAT

我们得到四条读段：

ATGC
TGCG
GCGA
CGAT

如果恰好知道它们原来的先后顺序，就很容易拼接：

ATGC
 TGCG
  GCGA
   CGAT

最后得到：

ATGCGAT

但这是一种极其理想的情况。实际测序时，读段可能按照随机顺序出现。我们可能先得到：

CGAT

随后：

ATGC

再得到：

GCGA

最后：

TGCG

计算机需要根据读段本身的序列关系，判断怎样排列。

为什么需要很多读段覆盖同一个位置？

因为测序数据可能出现错误。一个碱基可能被误读。某段DNA可能没有被成功测到。某些读段质量很差。因此，研究者通常希望同一个基因组位置，被多个相对独立的读段覆盖。这叫作：Coverage（覆盖深度）。如果基因组长度是：

G

我们一共测得：

N

条读段。每条平均长度为：

L

那么，平均测序深度常简化为：

\boxed{ C=\frac{NL}{G} }

例如：基因组长：

G=100\,000

测得：

N=10\,000

条读段。每条平均：

L=100

个碱基。于是：

C=\frac{10\,000\times100}{100\,000}

得到：

\boxed{C=10}

即平均约10×覆盖深度。但它只是平均值。不代表每个位置都刚好被覆盖10次。真实测序可能存在明显的覆盖偏差。

图54-1：多条读段覆盖同一段DNA

查看 SVG 配图原件

54-005-7e6ac26f6a9d.svg · 点击查看大图

实际组装时，完整DNA通常未知；图中的完整长条只是帮助理解读段覆盖关系。覆盖深度越高不一定越好，还要考虑重复区域、错误率、测序偏差和读段长度。

1995年的突破发生在哪里？

1995年发表的流感嗜血杆菌基因组研究，真正令人振奋的地方是：研究者利用随机获取的全基因组片段，结合计算与实验补缺，完成了一个约183万碱基的细菌基因组。这项工作证明：在合适的数据量、算法与实验验证条件下，不必先为整个基因组构建完整的逐段物理测序路线，也可以进行全基因组组装。但请注意：它并不是只靠一条神奇的字符串算法，就自动生成最终完整基因组。测序质量控制、重叠组装、缺口填补和实验验证都发挥了作用。

为什么把DNA打碎以后，问题反而更复杂？

因为原本连续的信息被转变成了许多局部观察。组装程序不能直接看到：读段A在读段B之前。读段C属于染色体的哪个位置。某个重复序列究竟出现了几次。它必须从重叠关系中推断顺序。而这种推断有时并不唯一。这就是组装问题的核心困难。

第二幕：最直观的办法，为什么会在重复序列面前失去方向？

假设我们手中有四条读段：

ATGC
TGCG
GCGA
CGAT

第一种直觉是：两条读段只要有足够长的前后缀重叠，就将它们连接起来。例如：ATGC的后三个字符是：

TGC

TGCG的前三个字符也是：

TGC

因此：

ATGC → TGCG

具有三字符重叠。同理：

TGCG → GCGA

以及：

GCGA → CGAT

于是可以形成一条路径：

ATGC → TGCG → GCGA → CGAT

这便是常见的：Overlap–Layout–Consensus（OLC）思想的简化版本。其中：Overlap：寻找读段之间的重叠。Layout：根据重叠关系排列读段。Consensus：根据多个读段确定最终的共识序列。

为什么这套方法看起来很合理？

因为它直接利用读段中真实存在的重叠信息。而且，不必先将全部读段统一切成固定长度短词。它也容易结合更复杂的重叠评分，以及测序错误模型。因此，OLC并不是一种落后的错误方法。在适合的测序数据和实现条件下，它长期具有重要价值。但是，它存在一项实际困难。如果有：

N

条读段。朴素地对每一对读段都检查重叠，候选对数量约为：

\frac{N(N-1)}2

数量级：

\boxed{O(N^2)}

实际组装器可以通过索引、种子查找和过滤等方法减少比较成本。所以不能说所有OLC实现都必然执行

N^2

次昂贵比对。但当读段数量非常庞大时，寻找和管理大量重叠关系，确实是一项工程挑战。

更麻烦的是重复序列

假设原始DNA包含两段完全相同的：

ATAT

而它们位于不同的位置。来自这两处的短读段，可能拥有完全相同的字符内容。对于只看到读段序列的程序来说，它们可能无法区分：某个ATAT究竟来自第一处。还是第二处。如果将全部相同内容简单视为同一个位置，就可能错误合并不同的基因组区域。如果将它们任意拆开，又可能错误增加重复次数。这便是：Repeat Problem（重复序列问题）。

2001年论文真正想改变什么？

Pevzner、Tang与Waterman在论文开头明确表示，经典片段组装长期使用OLC范式。他们提出：针对基因组随机测序和大量重复序列，可以使用另一种图表示。将读段转换成固定长度短词，再由短词的前后缀关系形成图。这样，原本基于读段之间连接关系的问题，就变成了在图中寻找符合边使用条件的路径。^[PMC]不过，这里需要提前纠正一个历史误解。2001年并不是科学家第一次想到用短词图解决DNA组装问题。1995年，Ramana M. Idury与Michael S. Waterman已经发表：A New Algorithm for DNA Sequence Assembly。他们将鸟枪法测序与固定长度短词集合的思想结合起来，研究了新的序列组装算法。^{[Sage Journals]}2001年的EULER，是在这些更早研究基础上发展并改进了欧拉路径式组装框架，特别关注重复问题与实际组装。历史上不能把这条研究路线全部归功于某一篇论文。

第三幕：1736年的七座桥，为什么会出现在21世纪的DNA组装论文里？

现在，我们突然把时间拨回1736年。距离DNA双螺旋结构被阐明，还有两个多世纪。那一年，数学家：莱昂哈德·欧拉（Leonhard Euler）研究了著名的：柯尼斯堡七桥问题。城市里的河流将土地分成几个部分。七座桥将不同土地连接起来。人们提出：能不能从某个地方出发，每座桥恰好经过一次，最后走完全部七座桥？欧拉发现，我们不需要记录每一座桥周围所有街道的具体形状。只要将：陆地表示成节点。桥表示成边。就可以将问题转化为图论问题。从此，一类与“经过边”有关的路径问题，拥有了清晰的数学形式。

什么是欧拉路径？

在一张图中，一条：经过每条边恰好一次的路径叫作欧拉路径。如果最终回到起点，就是：欧拉回路（Eulerian Circuit）。这里最重要的是：边恰好走一次。而不是：节点恰好访问一次。这与哈密顿路径（Hamiltonian Path）有根本区别。

路径类型	必须满足的访问条件
欧拉路径	每条边恰好经过一次
哈密顿路径	每个节点恰好访问一次

一个节点可以在欧拉路径中出现多次。例如：

A → B → C → A → D

其中A出现两次。这并不违反欧拉路径的定义。只要每条边都只使用一次即可。

有向图里的欧拉路径条件

对于我们后面使用的有向图，定义：

\deg^+(v)

表示节点v的出度。

\deg^-(v)

表示节点v的入度。如果存在一条覆盖全部边的有向欧拉路径，通常需要满足：起点：

\deg^+(s)=\deg^-(s)+1

终点：

\deg^-(t)=\deg^+(t)+1

其他节点：

\boxed{\deg^+(v)=\deg^-(v)}

如果所有节点的入度等于出度，则有可能存在欧拉回路。但要注意：度数条件是必要条件之一，还需要满足非孤立节点之间适当的连通性条件。只检查入度和出度，并不能保证能够经过所有边。

为什么DNA组装更适合研究“边”？

因为我们马上会构造一种图：每条边表示一个观察到的DNA短词。每个节点表示它的前后缀。如果我们希望使用全部短词的相应出现次数，就相当于希望：让路径经过每一条代表短词出现的边。因此，DNA重建问题与欧拉路径形成了直接联系。这不是因为DNA真的在数学图上移动。而是因为我们构造了一种表示，使DNA字符的连续重叠关系，对应了有向边之间的连接关系。这种表示方式，正是本章最重要的数学转化。

第四幕：1946年，德布鲁因的组合数学为什么成为DNA图的名字？

时间来到1946年。荷兰数学家：尼古拉斯·德布鲁因（Nicolaas Govert de Bruijn）发表了一篇题为：A Combinatorial Problem《一个组合数学问题》的论文。论文刊登于《荷兰皇家科学院科学院报》相关卷册，第49卷第7期，第758—764页。^{[埃因霍温大学研究门户]}这项研究与一种特殊的循环序列有关。

能不能构造一条尽可能短的序列，包含全部长度为k的短词？

假设字母表只有两种字符：

\{0,1\}

我们希望找到一个循环序列，让全部长度为：

k=2

的二进制短词，都能恰好出现一次。全部短词是：

一共有：

2^2=4

个。考虑循环序列：

如果允许从尾部重新接回开头，那么长度为2的连续循环短词分别是：

四种短词全部出现了。这类循环序列，后来被称为：de Bruijn Sequence（德布鲁因序列）。

为什么它与DNA有关？

DNA也由有限字母表构成。主要是：

\{A,C,G,T\}

如果我们想研究长度为：

k

的所有可能DNA短词，就会发现类似的组合结构。一个长度为k的短词：

ATG

可以拆成：前缀：

AT

以及后缀：

TG

于是可以画出一条有向边：

\boxed{AT\xrightarrow{ATG}TG}

它表示：如果当前已经到达前缀AT，那么沿着代表ATG的边，就可以转移到后缀TG。再看：

TGC

它对应：

\boxed{TG\xrightarrow{TGC}GC}

现在，这两条边可以连接：

AT\rightarrow TG\rightarrow GC

对应的DNA便是：

ATGC

这就是DNA组装中非常常见的一种：de Bruijn图（德布鲁因图）的边表示法。

要区分两种常见的图表示

在不同教材和论文中，de Bruijn图可能使用略有差别的表示方式。本章采用最适合教学的一种：

节点是长度为 $k-1$ 的字符串。
有向边代表一个长度为 $k$ 的k-mer。
边连接这个k-mer的前缀与后缀。

因此：

\boxed{ v_{\mathrm{from}}=x_1x_2\cdots x_{k-1} }

\boxed{ v_{\mathrm{to}}=x_2x_3\cdots x_k }

而整条边对应：

\boxed{x_1x_2\cdots x_k}

有些组装软件还会压缩无分叉路径，将多个连续k-mer合并为较长的unitig节点。那是进一步的图压缩表示，不等同于我们当前的最基本模型。

图54-2：一个3-mer怎样变成一条有向边对k=3，节点长度为2，边长度为3。沿边前进时，只需要把边对应短词的最后一个字符追加到已经构造的序列后面。

原始de Bruijn序列与基因组组装不是同一个问题

还有一项历史区别必须保留。1946年的de Bruijn研究，属于组合数学。它不是为了DNA测序而提出的。DNA组装算法使用的是与这一组合结构相关的图表示。而且，生物组装中的图，通常并不包含所有可能的k-mer。它只包含从测序数据中观察到或推断出的k-mer。因此，真实组装图未必是规则、完整的经典de Bruijn图。它可能存在：缺失的边。重复的边。错误产生的分叉。以及多个相互断开的连通分量。但德布鲁因的组合数学，为这些图的结构提供了一条重要的思想来源。

第五幕：把ATGCGAT切成五条边，一张图突然出现了

现在，我们回到本章开篇的小DNA。

\boxed{\text{ATGCGAT}}

设：

k=3

我们要做的第一件事，是提取其中所有长度为3的连续短词。原始序列长度：

L=7

线性序列中k-mer的出现次数为：

\boxed{L-k+1=5}

因此，得到：

ATG
TGC
GCG
CGA
GAT

每个k-mer对应一条有向边。

k-mer	前缀节点	后缀节点
ATG	AT	TG
TGC	TG	GC
GCG	GC	CG
CGA	CG	GA
GAT	GA	AT

于是，图中的边是：

AT\rightarrow TG

TG\rightarrow GC

GC\rightarrow CG

CG\rightarrow GA

GA\rightarrow AT

图54-3：ATGCGAT的3-mer图蓝色节点表示长度为2的字符串；每条有向边对应一个3-mer。节点AT同时是本例路径的起点和终点，因此图形成了一个有向回路。

一件有趣的事：这张图居然形成了回路

我们从：

AT

出发。沿：

ATG

到达：

TG

再沿：

TGC

到：

GC

如此继续。最后：

GA\xrightarrow{GAT}AT

回到了起点。因此，这张图存在一条欧拉回路。但这里要格外小心。图存在欧拉回路，并不意味着原始DNA一定是一个环状基因组。我们的输入明明是线性字符串：

ATGCGAT

之所以形成回路，是因为前两个字符和最后两个字符恰好都是：

AT

它们在k−1字符的节点表示里被合并了。这正好展示了一项重要的组装现象：k-mer图的路径结构，不一定能够唯一恢复原始DNA的生物学边界。如果知道原始起点是AT，并按上面边的顺序展开，就能恢复：

ATGCGAT

但是，如果不知道起点，则可以从回路中的不同位置开始读取。图本身不会告诉我们：原始线性DNA到底从哪里开始。这并不是算法实现错误。而是可观测信息不足造成的歧义。

换一条DNA，就能看见普通的欧拉路径

考虑：

GATTACA

k仍然等于3。得到：

GAT
ATT
TTA
TAC
ACA

对应：

GA → AT → TT → TA → AC → CA

这一次，起点：

GA

出度比入度多1。终点：

CA

入度比出度多1。其他节点入度等于出度。这便是一条典型的有向欧拉路径，而不是欧拉回路。

为什么知道起点以后，可以恢复完整序列？

假设我们从：

GA

开始。第一条边是：

GAT

它的最后一个字符是：

于是：

GA + T = GAT

第二条边是：

ATT

只需添加最后一个：

得到：

GATT

继续：

TTA → A
TAC → C
ACA → A

最终：

\boxed{\text{GATTACA}}

因此，重建序列的规则可以写成：

\boxed{ \text{完整序列} = \text{起始}(k-1)\text{-mer} + \text{每条路径边的最后一个字符} }

在没有额外闭合处理的线性展开表示里，如果路径使用了：

E

条边，那么输出字符串长度为：

\boxed{E+k-1}

这条公式将图的边数与重建DNA的长度连接了起来。

第六幕：重复序列第一次让我们看清——存在欧拉路径，不等于答案唯一

现在，真正困难的问题来了。假设我们已经把全部k-mer转换成一张图。也确认这张图存在欧拉路径。是不是就能宣布：完整DNA已经恢复成功？不一定。因为一张图可能存在不止一条欧拉路径。不同欧拉路径，可能对应不同的DNA序列。而这些DNA序列，却能够产生完全相同的k-mer集合，甚至相同的k-mer出现次数。

一组只有八个字符的反例

考虑第一条DNA：

AAAACAAC

以及第二条：

AAACAAAC

两者显然不同。但如果：

k=3

把它们分别拆成全部3-mer，会发生什么？第一条：

AAA
AAA
AAC
ACA
CAA
AAC

第二条：

AAA
AAC
ACA
CAA
AAA
AAC

你会发现：不仅短词的种类相同。就连每种短词出现的次数也完全相同。

3-mer	第一条出现次数	第二条出现次数
AAA	2	2
AAC	2	2
ACA	1	1
CAA	1	1

它们拥有相同的k-mer多重集合。甚至两条序列的开头两个字符都是：

AA

结尾两个字符也都是：

AC

所以，即使已知这些边界信息，仍然无法仅凭这组六条3-mer，确定哪一条是原始DNA。

图54-4：相同的短词集合，可以对应不同的完整序列

候选组装结果A

AAAACAAC

候选组装结果B

AAACAAAC

两者均可产生相同的3-mer多重集合与相同的首尾2-mer

这项歧义是算法不够聪明吗？

不是。假设我们只有上述k-mer数据。并且没有其他实验信息。那么，两条不同的原始DNA都与数据完全相容。这种情况下，没有任何算法能够从这些数据中，凭空获得原本不存在的区分信息。算法可以：列出多个可能结果。指出图中的分叉。利用额外假设给不同候选评分。或者要求更多实验数据。但它不能证明其中某一个候选必然是真实原序列。这叫作：不可辨识性（Non-identifiability）。对于实际基因组组装，重复区域造成的歧义，是最重要的困难之一。

那么，2001年的EULER究竟解决了什么？

这正是理解原始论文时必须谨慎的一点。Pevzner、Tang与Waterman并不是只提出：“找到任意一条欧拉路径即可。”他们进一步研究了：Eulerian Superpath Problem（欧拉超路径问题）。普通欧拉路径只要求经过图中的全部边。但真实测序读段具有比独立k-mer更多的信息。例如，一条长度为100的读段，本身已经告诉我们：其中连续的几十个k-mer，必须以特定顺序出现。如果我们将读段切碎成独立k-mer，却完全忘记这些k-mer原本属于同一条读段，就可能丢失解决重复区域的关键信息。因此，EULER还将：原始读段在图中的路径约束保留下来。再寻找一条既经过相关边、又满足读段路径约束的完整路径。原始论文将此表述为欧拉超路径问题，而非仅仅寻找一条任意欧拉路径。^{[布朗大学计算机科学系]}

一个重复区为什么需要更长的读段？

想象一段DNA结构：

左侧独特区域 → 重复区域 → 右侧独特区域

如果一条读段只落在重复区域内部，它可能无法区分这个重复序列来自哪个位置。但如果另一条更长的读段：从左侧独特区域开始。穿过整个重复区域。最终进入右侧独特区域。那么，这条读段就能够为重复区两侧的连接关系提供额外信息。它告诉算法：这个入口应该与那个出口相连。这并不能保证所有重复问题都能解决。如果重复区域比所有可用读段和片段跨度都长，歧义仍然可能存在。但它指出了正确方向：面对不可辨识性，真正有价值的办法不是凭空猜测，而是引入能够区分候选的新观测。

第七幕：互动实验室——亲手把DNA切成短词，再沿欧拉路径拼回去

现在，我们把本章最核心的图论过程，做成一个可逐步操作的实验。你可以选择不同的DNA。改变k-mer长度。尝试删去一个短词，模拟局部数据缺失。然后观察：图中出现哪些节点？每条边对应什么k-mer？是否存在覆盖全部边的欧拉路径？以及最终能够拼出怎样的序列？

互动实验：De Bruijn图与欧拉路径

图54-5：DNA组装图论实验台把已知的教学序列切成k-mer，打乱短词顺序，再由图寻找欧拉路径。原始DNA仅用于检查重建是否正确。

选择DNA样本

k-mer长度

3

删除一条k-mer，模拟局部数据缺失

不同节点

5

k-mer边

5

欧拉路径

存在

教学用原始DNA

ATGCGAT

打乱顺序后的k-mer多重集合

ATG

CGA

GAT

GCG

TGC

De Bruijn有向图节点为(k−1)-mer，连线为k-mer边。重复k-mer在图中按次数汇总显示。橙色边表示已经用于当前回溯顺序，线条并不表示道路的真实长度或DNA的三维结构。

欧拉路径逐步展开

已走过的边0 / 5

当前拼接结果

AT

从起始(k−1)-mer开始，尚未添加任何边。

本实验使用确定性Hierholzer算法寻找一条合法欧拉遍历。为了演示，数据来自一条已知字符串，并假设每个k-mer出现次数都恰好已知、没有测序冗余。实际读段的重复观测并不等于该k-mer在基因组中的真实拷贝数。若存在多条欧拉路径，这里只显示其中一条。

实验一：首尾重合并不意味着生物学上的环状DNA

选择：ATGCGAT以及：

k=3

观察图中的首尾节点。它们都是：

AT

因此，图形成回路。程序可以找到一条欧拉回路。但是，原始教学序列依然是一个线性字符串。这提醒我们：k-mer图本身可能丢失序列的绝对起止位置。

实验二：改变k的大小

将：

k=3

改成：

k=5

你会发现：k-mer数量减少。节点变得更长。短词重叠需要匹配更多字符。这通常可以帮助减少由较短重复结构引起的歧义。但真实测序中的k也不能任意增大。因为如果读段长度有限，或者测序错误较多，过大的k可能导致许多真实重叠无法建立连接。

实验三：删除一条边

勾选：“删除一条k-mer”。观察图发生了什么。可能出现：路径被截断。某些节点失去平衡。或者图被分成无法相互连接的部分。但要注意，删除一条边不一定总会让欧拉路径消失。有些图删除一条边之后，仍然存在一条覆盖剩余全部边的欧拉路径。只是，这条新路径对应的信息已经不完整。因此：图上存在合法组装路径，不等于测序数据已经覆盖真实基因组的全部区域。

实验四：重复结构为什么造成歧义？

选择：“存在歧义”。尝试改变k-mer长度。你会看到一些节点拥有多条候选边。这些分叉意味着：算法可能需要额外的读段顺序信息，才能确定正确的连接方式。欧拉路径保证的是一种图论上的边覆盖。但它不能无条件保证拼出的DNA就是原始的唯一真实序列。

第八幕：如果有上亿条k-mer，怎样在线性时间里找到欧拉路径？

现在，我们已经有了一张有向图。但还没有真正解释：如何高效找到一条经过所有边的欧拉路径？如果使用最直接的暴力方法，可能在遇到分叉时尝试：第一条边。第二条边。第三条边。然后不断回溯，尝试各种可能的完整边顺序。这显然可能十分昂贵。幸运的是，欧拉路径属于一类能够高效求解的经典图论问题。

1873年：Hierholzer算法

19世纪，数学家：卡尔·希尔霍尔策（Carl Hierholzer）相关的欧拉回路构造方法得到发表。今天，寻找欧拉路径时，一种非常常用的方法就是：Hierholzer Algorithm。它的核心思想是：沿着尚未使用的边前进。如果走到没有可用边的节点，就将当前节点加入最终路径。再向前回退。通过栈管理尚未完成的路线，最终反向构造整个欧拉路径。

为什么不能简单地“一路贪心向前走”？

假设一张图中某个节点有两条可选边。第一条通向某个局部回路。第二条通向真正的终点。如果直接选择第二条，并且只记录当前行走结果，就可能提前走到终点，却留下第一条回路尚未访问。因此，真正的欧拉遍历算法，需要能够在局部路线结束后，将其他回路插入最终结果。Hierholzer的栈与逆序输出方式，恰好完成了这项工作。

一个适合实现的算法结构

我们保存：stack：当前正在探索的节点路径。path：已经确认的逆序结果。如果栈顶节点还有尚未使用的出边：沿着一条边前进。将新节点压栈。如果栈顶节点已经没有未使用出边：弹出节点。加入最终结果。反复执行，直到栈为空。最后将结果反转。

为什么必须给每条边一个独立编号？

这是组装任务里一个容易犯错的地方。假设：

AAA

出现两次。那么在多重图里，它代表两条具有相同起点和终点、但具有不同身份的边。如果我们只使用：

visited["AAA"] = True

就可能把两次出现错误地当成一条边。因此，需要使用边的唯一编号：

edge 0: AAA
edge 1: AAA

欧拉路径必须按问题给定的边重数使用这些边。它们虽然字符完全相同，但在多重图中仍然是不同的边。

算法复杂度

如果使用合适的邻接表，并确保每条边只取出一次，那么Hierholzer算法的时间复杂度为：

\boxed{O(V+E)}

其中：

V

是节点数量。

E

是边数量。空间复杂度通常也是：

\boxed{O(V+E)}

这意味着：寻找一条合法欧拉路径，本身可以非常高效。但我们必须保持清醒。在真实基因组组装中，困难往往不是单纯找到任意一条欧拉路径。更难的是：测序错误怎样修正？重复边的真实出现次数怎样估计？多条欧拉路径哪一条更符合原始读段？哪些分叉属于真实生物变异？以及有没有足够的信息区分不同候选？这些都超出了基础Hierholzer算法的职责。

第九幕：2001年的EULER，为什么不是“找到任意一条欧拉路径就大功告成”？

2001年的EULER论文，是本章真正的历史主角。现在，我们已经有能力更准确地理解它。作者在论文开篇使用了一个很生动的拼图类比。传统组装方法像是：拿到一堆拼图碎片。逐对寻找碎片之间的重叠。再设法把它们拼起来。而他们反过来做了一件看似不合理的事。将原本已经破碎的片段，再切成更小且长度统一的片段。为什么这样做反而可能更好？因为一旦将序列组装转变成固定长度短词组成的图，就可以利用欧拉路径算法的高效结构。但论文的研究内容远远不只是这一项数学转化。

第一项：重复区域会造成入口与出口配对歧义

一条重复序列可能被多个不同区域连接。在图上，它可能形成：两个入口。两个出口。我们需要判断：第一个入口究竟应该连接第一个出口？还是第二个出口？如果没有足够长的读段跨越重复区域，就可能无法判断。因此，EULER使用原始读段在图中的路径信息，约束可能的出口选择。这形成：Eulerian Superpath Problem。原始论文将其作为核心的广义图论问题。^{[布朗大学计算机科学系]}

第二项：测序错误会制造虚假的边

如果真实DNA是：

AATGCA

某条错误读段却显示：

AATACA

那么，使用：

k=3

时，真实短词为：

AAT ATG TGC GCA

错误短词则可能包含：

AAT ATA TAC ACA

它们共享：

AAT

并且在图的后面重新汇合。于是，形成一条局部的替代路径。这种图结构经常被形象地称为：Bubble（气泡）。

图54-6：一次测序替换错误可能产生图中的气泡

查看 SVG 配图原件

54-013-11ce29b492f3.svg · 点击查看大图

为便于阅读，将若干连续节点压缩成两条路径。真实测序错误可能生成气泡，也可能生成短分叉死端；真实遗传变异也可能形成类似气泡，因此不能见到气泡就一律删除。

如果大量高质量读段支持上方路径，而只有极少数低质量读段支持下方路径，那么，程序可能有理由认为下方路径属于测序错误。但如果两条路径都有充分支持，也可能存在真实的生物学差异。因此，错误修正不能简单地等同于：永远选择出现次数最多的边。还需要考虑数据来源、覆盖分布、测序质量及研究目标。

第三项：原始论文的实际数据实验

2001年论文不仅讨论理论，还评估了EULER在大型测序项目上的表现。作者分析了真实测序错误对短词集合的影响。并将其与利用原始读段信息进行的纠错联系起来。论文中一项误差修正分析报告：原始测序数据平均每条读段约有4.8处错误。在所讨论的处理结果中，平均错误数降低到约0.11处。这属于论文采用的特定测试数据与纠错流程结果。不能当成所有测序技术或所有组装器的一般保证。^{[布朗大学计算机科学系]}

为什么作者把“重复问题”看得那么重要？

因为重复序列会造成真正的信息歧义。即使数据完全没有测序错误，图中也可能存在多条与短词集合一致的路径。这意味着：纠正错误，不等于解决重复。二者是不同的问题。错误修正解决的是：哪些碱基可能被测错？重复解析解决的是：多个相同或相似片段，究竟应该连接到哪些不同的基因组位置？EULER通过图结构与读段路径约束，将这两项问题放进同一套算法框架。但不能把原论文关于特定测试中解决重复问题的成果，夸大成：它在数学上保证能够唯一恢复任意复杂、任意重复的真实基因组。如果数据缺乏足够的跨重复信息，唯一重建本来就可能不可能。

第十幕：2008年，Velvet如何处理短读段带来的错误与分叉？

2001年之后，基因组测序技术继续发展。越来越多测序设备能够在一次运行中产生大量短读段。在某些早期高通量测序平台上，一条读段可能只有几十个碱基。例如：25个碱基。36个碱基。50个碱基。相比此前常见的较长测序读段，这些短读段数量可能更多，但每条读段所包含的连续上下文更少。于是，一项新的组装困难出现了。同样的基因组，现在可能被切成数量巨大、长度很短的片段。这使de Bruijn图重新获得了重要地位。

2008年：Velvet

2008年3月18日，《Genome Research》在线发表了一项研究：Velvet: Algorithms for de Novo Short Read Assembly Using de Bruijn Graphs《Velvet：利用德布鲁因图进行从头短读段组装的算法》。作者是：丹尼尔·泽尔比诺（Daniel R. Zerbino）。以及：尤恩·伯尼（Ewan Birney）。论文后来刊登于该刊第18卷第5期，第821—829页。^[PMC]

他们面对的是：大量短读段。大量重叠k-mer。测序错误制造的图分叉。以及基因组中的真实重复区域。因此，Velvet不仅构造de Bruijn图。还特别重视：图简化与纠错。

第一项操作：Tip Removal（剪除短死端）

假设一条本来连续的真实路径，在某个位置因为测序错误，出现了一条很短的分叉。这条错误分叉可能只延伸几个k-mer。然后突然终止。这就像一棵树上多出来的一根短枝。在组装图里，这种结构常被称为：Tip（短死端）。如果它具有：很短的长度。很低的读段支持。以及其他符合错误模型的特征。那么，它可能是测序错误产生的。组装算法可以考虑将其剪除。但必须谨慎。因为某些真实基因组区域也可能拥有低覆盖，甚至形成短分支。因此，实际纠错不能简单地规定：“凡是短分支，一律删除。”

第二项操作：Bubble Removal（气泡处理）

另一类错误是：两条路径从同一个节点分开。经过不同的中间序列。最后又重新汇合。这便是刚才看到的气泡结构。Velvet提出的：Tour Bus方法，会寻找某些相似的替代路径，并结合覆盖支持与序列关系进行合并和纠错。原始论文描述，这项方法使用一种受图搜索思想启发的策略，比较候选分支的长度与支持情况。^[PMC]需要注意：气泡也可能来自真实的生物变异。例如，在含有多个等位基因的样本中，两条路径可能对应真实的序列差异。因此，把所有气泡都合并成一条路径，未必符合研究目标。特别是当目标涉及：变异检测。多倍体组装。或者区分不同单倍型。时，真实分支信息可能非常重要。

第三项操作：Unitig Compression（无分叉路径压缩）

考虑：

AA → AT → TG → GC → CA

如果中间每个节点：只有一条入边。也只有一条出边。且没有其他分叉。那么，这条路径可以压缩成一个较长的连续序列单元。这种较长、没有内部选择歧义的图路径，通常称为：Unitig。压缩后的图可以减少节点数量，使后续分析更加高效。但这里的Unitig并不自动等于完整基因。也不自动等于整条染色体。它只是图结构中能够在当前信息下明确连接的一段序列。

k为什么不能随意选择？

Velvet论文非常明确地讨论了k-mer长度的权衡。在一个具体实现中，对于25碱基的读段，研究者可能选择接近21的k。因为较大的k能够减少一些不必要的图连接。但如果k过大，短读段可能提供不了足够多的重叠。而测序错误还会破坏包含错误位置的多个k-mer。因此：较小的k，更容易建立连接，但可能增加重复歧义。较大的k，可以减少部分重复歧义，但可能增加图断裂和对错误的敏感性。这项权衡，与我们在第52章BLAST里选择种子长度时遇到的问题非常相似。短词长度是算法参数，同时也是信息保留方式的选择。

第十一幕：2012年的SPAdes，为什么不能只看一张k-mer图？

到2012年，基因组组装研究面临更多复杂场景。其中一项是：单细胞测序（Single-cell Sequencing）。如果我们从大量细胞中获取DNA，某些测序偏差可能通过样本量和重复测量有所缓解。但是，如果研究对象是一个细胞中非常有限的DNA，数据可能存在更加明显的不均匀覆盖。某些区域被大量测到。某些区域却几乎没有读段。这对组装算法提出了新的要求。

SPAdes的出现

2012年，Anton Bankevich、Sergey Nurk、Dmitry Antipov等研究者发表：SPAdes: A New Genome Assembly Algorithm and Its Applications to Single-Cell Sequencing《SPAdes：一种新的基因组组装算法及其在单细胞测序中的应用》。论文发表于《Journal of Computational Biology》第19卷第5期，第455—477页。其中也包括Pavel Pevzner等共同作者。^[PubMed]SPAdes针对复杂覆盖情况，发展了多k-mer图处理、错误校正及其他组装策略。它不是简单地执行：建立一张图。找到任意欧拉路径。然后结束。真实组装通常需要多个阶段共同配合。

为什么多个k可能比单独一个k更有价值？

假设较小的k能够连接一些低覆盖区域。但较小的k又造成了太多重复分叉。于是，我们希望利用更大的k，帮助区分某些重复结构。可是，如果只使用大k，某些低覆盖区域可能直接断裂。因此，多k策略尝试利用不同k提供的互补信息。这个思想与图像处理中使用不同尺度识别结构，具有某种方法论上的相似性。不过，基因组组装中的多k不是简单地将不同k的图随意相加。不同尺度的图必须通过专门的算法组织、校正和连接。

2016年：长读段又让研究者重新考虑重叠图

随着单分子测序技术发展，科学家开始获得长度远超早期短读段的测序结果。例如：数千碱基。甚至更长的连续读段。长读段可能跨越很多短读段无法跨越的重复区域。这就重新改变了组装算法的权衡。2016年，Heng Li发表：Minimap and Miniasm: Fast Mapping and de Novo Assembly for Noisy Long Sequences。介绍了面向长读段的快速映射与组装方法。Miniasm使用读段之间的重叠关系构造组装图，并特别强调运行速度。^[PMC]这项发展说明：De Bruijn图并没有永久替代所有重叠组装方法。对于不同测序技术，算法设计仍然需要考虑：读段长度。测序错误类型。覆盖深度。重复结构。内存限制。以及对最终碱基准确度的要求。即使在今天，基因组组装也不是一项已经由单一算法彻底解决的问题。

第十二幕：Python实战——用Hierholzer算法把DNA短词重新拼接起来

现在，我们真正实现一套可以运行的微型DNA组装器。它不依赖Biopython。也不需要额外的图论库。只使用Python标准库。但它会完整实现本章最重要的计算过程：提取k-mer。构造有向多重图。计算每个节点的入度和出度。利用Hierholzer算法寻找欧拉路径。最后将路径重新转换成DNA。

第一部分：提取全部k-mer


from collections import Counter, defaultdict


def kmers_from_sequence(sequence, k):
    """
    从一条已知DNA中提取全部长度为k的短词。

    保留每个k-mer的出现次数。
    """

    sequence = sequence.upper()

    if not sequence:
        raise ValueError("DNA不能为空")

    if set(sequence) - set("ACGT"):
        raise ValueError("DNA仅允许A、C、G、T")

    if not 2 <= k <= len(sequence):
        raise ValueError("k必须介于2和序列长度之间")

    return [
        sequence[i:i + k]
        for i in range(len(sequence) - k + 1)
    ]


dna = "ATGCGAT"

kmers = kmers_from_sequence(
    dna, k=3
)

print(kmers)

assert kmers == [
    "ATG",
    "TGC",
    "GCG",
    "CGA",
    "GAT",
]

注意，这里的函数接收一条已经知道完整顺序的DNA。它是为了模拟：“如果真实序列被转换成短词，我们能不能再把它恢复？”而设计的。真实组装中，通常只拥有读段，不拥有完整参考序列。因此，后面的重建测试属于可验证的教学实验，而不是真实未知基因组的直接组装。

第二部分：建立有向多重图


def build_debruijn_graph(kmers):
    """
    使用(k-1)-mer作为节点，
    每一个k-mer出现都产生一条独立有向边。

    返回：
        adjacency
        indegree
        outdegree
        edges
    """

    if not kmers:
        raise ValueError("至少需要一个k-mer")

    k = len(kmers[0])

    if any(len(word) != k for word in kmers):
        raise ValueError("全部k-mer长度必须相同")

    adjacency = defaultdict(list)

    indegree = Counter()
    outdegree = Counter()

    edges = []

    for edge_id, word in enumerate(kmers):
        source = word[:-1]
        target = word[1:]

        edges.append({
            "id": edge_id,
            "word": word,
            "source": source,
            "target": target,
        })

        adjacency[source].append(edge_id)

        outdegree[source] += 1
        indegree[target] += 1

        # 终点即使没有出边，也需要保留
        adjacency[target]

    return (
        adjacency,
        indegree,
        outdegree,
        edges,
    )


kmers = [
    "ATG",
    "TGC",
    "GCG",
    "CGA",
    "GAT",
]

graph = build_debruijn_graph(kmers)

adjacency, indegree, outdegree, edges = graph

for edge in edges:
    print(
        edge["source"],
        "->",
        edge["target"],
        "标签:",
        edge["word"],
    )

预期输出：

AT -> TG 标签: ATG
TG -> GC 标签: TGC
GC -> CG 标签: GCG
CG -> GA 标签: CGA
GA -> AT 标签: GAT

这里：每个k-mer出现都分配一个独立的：

edge_id

这非常重要。例如：

AAA
AAA

虽然两个短词完全相同，但在我们的数学多重图里，它们仍然可以是不同的边。

第三部分：寻找有向欧拉路径


def find_eulerian_edge_path(
    adjacency,
    indegree,
    outdegree,
    edges,
):
    """
    Hierholzer算法，返回边ID的欧拉遍历顺序。

    验证：
        1. 节点入度与出度条件
        2. 最后确实使用全部边

    总体时间复杂度 O(V+E)。
    """

    if not edges:
        return []

    nodes = set(adjacency)

    start_candidates = []
    end_candidates = []

    for node in nodes:
        diff = outdegree[node] - indegree[node]

        if diff == 1:
            start_candidates.append(node)

        elif diff == -1:
            end_candidates.append(node)

        elif diff != 0:
            return None

    valid_path = (
        len(start_candidates) == 1
        and len(end_candidates) == 1
    )

    valid_circuit = (
        len(start_candidates) == 0
        and len(end_candidates) == 0
    )

    if not (valid_path or valid_circuit):
        return None

    if valid_path:
        start = start_candidates[0]

    else:
        # 对欧拉回路，选择一个有出边的节点。
        # 这个起点选择不代表已确定DNA的
        # 原始生物学边界。
        start = edges[0]["source"]

    # 拷贝邻接表，避免破坏调用者的数据
    remaining = {
        node: edge_ids[:]
        for node, edge_ids in adjacency.items()
    }

    # 栈中记录：当前节点、抵达它的边编号
    stack = [(start, None)]

    reversed_path = []

    while stack:
        node, incoming_edge = stack[-1]

        if remaining[node]:
            edge_id = remaining[node].pop()

            next_node = edges[edge_id]["target"]

            stack.append((
                next_node,
                edge_id,
            ))

        else:
            node, incoming_edge = stack.pop()

            if incoming_edge is not None:
                reversed_path.append(incoming_edge)

    # 如果仍有边无法从起点覆盖，
    # 即使度数条件正确，也不能视为完整组装。
    if len(reversed_path) != len(edges):
        return None

    reversed_path.reverse()

    return reversed_path

这里还有一项不太显眼、但非常重要的正确性检查。即使每个节点满足相应的入度和出度条件，也可能存在彼此断开的图分量。例如：一个小回路。加上另一个完全不相连的小回路。两个回路内部都满足：入度等于出度。但不存在一条连续路径，能够从一个回路走到另一个回路，并使用全部边。因此，我们最后必须检查：

len(reversed_path) == len(edges)

如果无法使用全部边，就不能宣称成功得到完整欧拉遍历。

第四部分：从边路径恢复DNA


def assemble_from_kmers(kmers):
    """
    假设每一个输入k-mer出现次数，
    都对应真实序列中的一次出现。

    返回：
        一条与全部输入k-mer多重集合一致的字符串。

    如果没有覆盖全部边的欧拉路径，返回None。
    """

    if not kmers:
        return None

    graph = build_debruijn_graph(kmers)

    adjacency, indegree, outdegree, edges = graph

    path = find_eulerian_edge_path(
        adjacency,
        indegree,
        outdegree,
        edges,
    )

    if path is None:
        return None

    if not path:
        return None

    first_edge = edges[path[0]]

    result = first_edge["source"]

    for edge_id in path:
        result += edges[edge_id]["word"][-1]

    return result


dna = "GATTACA"

kmers = kmers_from_sequence(dna, k=3)

# 模拟打乱短词顺序
shuffled = sorted(kmers, reverse=True)

reconstructed = assemble_from_kmers(shuffled)

print("原始DNA:", dna)
print("打乱k-mer:", shuffled)
print("组装结果:", reconstructed)

assert reconstructed == dna

预期输出：

原始DNA: GATTACA
组装结果: GATTACA

这说明：对于这一组短词，即使打乱原来的出现顺序，算法仍然可以利用前后缀关系恢复一条与它们一致的完整序列。

第五部分：回到ATGCGAT


dna = "ATGCGAT"

kmers = kmers_from_sequence(dna, 3)

assembled = assemble_from_kmers(kmers)

print("原始:", dna)
print("组装:", assembled)

assert Counter(
    kmers_from_sequence(assembled, 3)
) == Counter(kmers)

这里的断言故意没有仅仅检查：

assembled == dna

因为本例的k-mer图形成欧拉回路。如果程序从不同节点开始，可能得到一个循环位移后的字符串表示。它们不一定与输入的线性字符串逐字相同。但仍然可能使用相同的k-mer边集合。这再次说明：组装图中的循环，并不自动指定DNA的绝对起点。如果额外已知序列起点，便可以利用这项信息选择相应的表示。

第六部分：验证同一组k-mer可能产生两种不同DNA


a = "AAAACAAC"
b = "AAACAAAC"

kmers_a = kmers_from_sequence(a, 3)
kmers_b = kmers_from_sequence(b, 3)

assert a != b

assert Counter(kmers_a) == Counter(kmers_b)

assert a[:2] == b[:2]
assert a[-2:] == b[-2:]

print("序列A:", a)
print("序列B:", b)

print("共同k-mer多重集合:")

for word, count in sorted(Counter(kmers_a).items()):
    print(word, count)

print("两条不同DNA，拥有相同3-mer多重集合")

这个测试非常有价值。因为它证明：即使一个欧拉组装程序运行完全正确。也不代表仅凭这份数据一定可以唯一恢复原始序列。有时候，多条不同的DNA都能够满足全部已知的图边约束。我们需要更多信息。例如：更长的读段。已知的读段路径。成对末端测序提供的距离与方向信息。或者其他独立实验约束。

第七部分：完整性验证应该怎样做？

对于一个合法的数学组装结果，我们至少可以检查：输入k-mer的出现次数，是否与输出序列拆解得到的k-mer多重集合一致？


def verify_assembly(kmers, assembled):
    if assembled is None:
        return False

    k = len(kmers[0])

    reconstructed_kmers = kmers_from_sequence(
        assembled, k
    )

    return (
        Counter(kmers)
        == Counter(reconstructed_kmers)
    )


for dna in [
    "GATTACA",
    "ATGCGAT",
    "ATATAT",
    "AAAACAAC",
]:
    k = 3
    reads = kmers_from_sequence(dna, k)

    result = assemble_from_kmers(reads)

    assert verify_assembly(reads, result)

    print(dna, "=>", result)

print("全部基础组装验证通过")

这套验证可以检查：输出序列是否与输入短词多重集合一致。但它不能证明：真实测序数据没有错误。基因组没有缺失区域。重复次数一定正确。以及输出序列就是唯一真实染色体。因为这些问题，需要超出当前数学模型的额外证据。

一个非常重要的工程警告

真实测序数据中，同一个k-mer可能因为覆盖深度较高，被读取几十次。但是，这并不代表原始基因组里一定包含几十份相同k-mer。因此，不能简单地将每一条读段里出现的全部k-mer，都当作需要在最终欧拉路径中各走一次的真实基因组边。真实组装器需要利用覆盖深度、测序错误模型、重复结构和其他信息，推断图中的真实路径及重复次数。本章的Python实现适合于：已知短词多重集合代表原序列中的实际出现次数。且没有测序错误、缺失或未知方向问题。的教学场景。这是一个完整、可运行的图论算法实验，但不是可以直接用来组装真实未知基因组的工业级软件。

第十三幕：拼出来的DNA，怎样判断究竟好不好？

假设我们的程序已经成功输出一段DNA：

GATTACA

而且它的k-mer多重集合，与输入完全一致。是不是就说明组装任务已经彻底完成？对于本章的七字符教学例子，我们可以直接检查原始序列，判断是否相同。但真正的基因组组装通常没有这样的便利。因为：未知的真实基因组，本来就是我们希望得到的结果。因此，组装质量评估必须利用一组不同的指标。

Contig：能够连续确定的一段序列

组装器可能无法恢复整条染色体。它可能只能得到一些相对可靠的连续片段。这些片段通常称为：Contig（重叠群或连续组装片段）。例如，一次组装得到了：

Contig 1: 40,000 bp
Contig 2: 30,000 bp
Contig 3: 20,000 bp
Contig 4: 10,000 bp

它们合计：

100\,000\text{ bp}

但程序还不知道这些片段应该如何连接成完整染色体。这并不一定代表组装失败。它可能意味着：某些重复结构尚无法明确解析。或者测序数据不足以连接这些区域。

Scaffold：利用额外信息连接Contig

如果我们拥有成对末端读段或者其他长距离连接信息，就可能进一步判断：Contig 1与Contig 2应该相邻。中间可能存在一个尚未确定的区域。于是，可以将多个Contig连接成更长的：Scaffold（支架序列）。其中，尚未确定的碱基位置，常使用：

表示。所以，Scaffold很长，并不意味着其中每一个碱基都已经被真实测定并准确恢复。

N50：为什么组装报告里经常出现这个数字？

假设组装生成四条Contig，长度分别是：

40,\ 30,\ 20,\ 10

单位暂时记为kb。总长度：

100\text{ kb}

我们从最长的Contig开始累计。第一条：

40

还没有达到总组装长度的一半：

50

加入第二条：

40+30=70

已经达到50。因此：

\boxed{N50=30\text{ kb}}

也就是：使累计长度首次达到总组装长度50%时，当前那条Contig的长度。

图54-7：四条Contig怎样计算N50

查看图表数据

contig	Contig长度（kb）
Contig 1	40
Contig 2	30
Contig 3	20
Contig 4	10

组装总长度100 kb

累计达到50%时的Contig长度

30 kb

N50越高，就一定越好吗？

不一定。假设一个组装器错误地将两段毫无关系的DNA拼成了一条很长的序列。它可能获得一个非常漂亮的N50。但这种结果存在：Misassembly（错误组装）。因此，仅靠N50不能判断组装正确性。真正的质量评价，还需要考虑：组装总长度是否合理。是否存在大量缺口。测序读段能否一致地映射回组装结果。重复区域是否被错误合并或展开。有没有错误连接不同染色体区域。以及预期存在的基因和保守序列是否得到恢复。所以，组装评价同样体现出一项重要原则：不能把单个容易计算的数字，直接当作整个复杂任务的最终质量。

这与我们前两章的统计评价有什么联系？

第52章BLAST告诉我们：一个高分片段，不一定具有生物学显著性。第53章HMM告诉我们：模型对训练数据的似然提高，不一定意味着对未知数据的预测更准确。而本章组装问题告诉我们：一条Contig很长，不一定说明序列组装正确。这些例子都体现了一项共同思想。算法优化的数值目标，与科研真正关心的正确性，必须明确区分。

第十四幕：生命序列的组装，为什么最终是一场信息与算法共同完成的工作？

现在，我们已经经历了从1977年到2010年代的一条算法发展道路。我们看到：实验技术能够读取DNA片段。计算机能够寻找它们的重叠。图算法能够利用k-mer建立有向结构。欧拉路径能够在某些模型中高效恢复合法的字符串。而真实组装还需要处理测序错误、重复区域和不完整覆盖。这让我们可以重新思考一个问题：组装成功，究竟主要依靠测序技术，还是算法？答案不是二选一。

如果测序数据缺少关键信息，再优秀的算法也可能无法唯一恢复

回忆我们构造的两条序列：

AAAACAAC
AAACAAAC

当：

k=3

时，它们具有相同的k-mer多重集合。如果数据只有这些短词，那么算法无法仅靠它们证明哪一条才是真实来源。但假如我们得到一条更长的读段，能够直接覆盖它们之间不同的连接关系，就可能区分两种候选。因此，增加有用的测序跨度，可以帮助减少组装歧义。

如果算法设计很差，丰富的数据也可能得不到正确利用

另一方面，即使已经得到很多高质量读段，如果算法：错误地合并重复。忽略读段质量。随意删除真实分支。或者不正确处理图中的边重数。也可能产生错误组装。因此，测序数据与算法是相互依赖的。实验提供证据。算法组织证据。统计方法评价证据。最后，还需要独立检验。

一项从18世纪一直延续到现代的思想

1736年，欧拉研究的是七座桥。1946年，德布鲁因研究的是循环短词的组合结构。1977年，桑格团队完成了一个DNA基因组的完整测序。1995年，科学家将全基因组随机测序应用于完整细菌基因组。同年，Idury与Waterman发表了新的DNA序列组装算法。2001年，Pevzner、Tang与Waterman进一步发展欧拉路径式的组装框架。2008年，Velvet将de Bruijn图应用于大量很短的读段，并研究图纠错与重复解析。2012年，SPAdes继续拓展复杂测序条件下的图组装方法。2016年，长读段组装研究又展示了重叠图及其他方法的重要价值。这不是一条简单的“新算法不断淘汰旧算法”的直线。而是科学家不断根据新数据的特点，选择新的问题表示和计算方法。有时，我们需要从读段之间的重叠关系出发。有时，需要将读段切成更小的k-mer。有时，需要构造图。有时，需要压缩图。有时，还需要使用长距离信息，把局部无法确定的片段重新连接起来。最重要的创新，往往发生在重新表示问题的那一刻。

尾声：拼图真正的秘密，是把碎片变成一张图

1977年，剑桥。桑格和同事们完成了φX174噬菌体的完整DNA基因组序列。五千多个核苷酸，组成了一个连续的生命信息记录。但是，他们并不是直接获得一条完整的DNA字符带。实验提供了局部片段。研究者需要将这些证据组织起来，确定它们的相互关系。到了1995年，流感嗜血杆菌基因组的完整测序与组装，展示了全基因组随机测序方法的实际能力。曾经需要对小型病毒基因组进行的序列重建工作，已经扩展到数百万碱基规模。而数据规模的增长，又不断暴露新的困难。读段可能发生重叠。也可能缺失。可能出现错误。还可能来自多个完全相同的重复区域。如果仅仅把所有读段排成一列，然后依次尝试拼接，搜索空间可能十分复杂。于是，研究者开始寻找另一种表示方法。将DNA拆成固定长度的短词。把短词的前后缀看作节点。把每个短词看作一条有向边。让原本复杂的字符串排列问题，转化成图中的路径问题。这时，一项来自18世纪的数学思想突然重新出现。欧拉路径。如果我们的图满足适当条件，就能够使用高效算法，找出一条经过每条边恰好一次的路径。而沿着这条路径，我们可以逐步恢复一段DNA序列。但是，科学家的工作并没有因此结束。因为有些图拥有多条不同的合法路径。有些图因为缺失读段而断裂。有些图因为测序错误产生虚假分支。还有些图，即使算法完全正确，也缺少判断真实DNA连接顺序所需的信息。这提醒我们：数学上的合法路径，不一定就是生命真实的历史。

本章最值得记住的三个公式

第一个，DNA序列中k-mer的数量：

\boxed{E=L-k+1}

其中：

L

是原始线性DNA的长度。

k

是短词长度。该公式要求序列长度至少为k，并且每个连续位置都提取一次短词。第二个，k-mer到图边的转化：

\boxed{ x_1x_2\cdots x_k : (x_1\cdots x_{k-1}) \rightarrow (x_2\cdots x_k) }

这个规则将字符串重叠，转化成有向边之间的连接关系。第三个，欧拉路径遍历的时间复杂度：

\boxed{O(V+E)}

这说明：当问题已经转化成一张合适的图，而且只要求找到一条合法欧拉路径时，可以使用非常高效的算法。但它并不意味着，真实基因组组装的全部困难都能够在：

O(V+E)

时间内无条件解决。因为实际组装还包括错误修正、重复解析和大量不完整证据处理。

从第51章到第54章，我们究竟学会了什么？

《破解生命的密码》· 前四章

第51章

序列比对在两条已有序列之间，寻找最佳的全局或局部对应关系。

第52章

BLAST在大型序列数据库中，用短词种子与启发式延伸快速寻找值得研究的候选。

第53章

隐马尔可夫模型从观察到的DNA字符，推断隐藏状态与概率结构，并估计模型参数。

第54章

基因组组装从不完整的DNA读段出发，使用重叠关系、de Bruijn图和欧拉路径，尝试重建原始序列。

这四章，形成了一条完整的算法发展路线。第一步：比较已有的生命序列。第二步：在海量序列数据库中寻找相似片段。第三步：推断序列背后无法直接观察的统计状态。第四步：从不完整的分子片段中，尝试恢复原始连续信息。它们分别使用了：动态规划。启发式搜索。概率模型。以及图算法。但无论是哪一种算法，都无法脱离真实数据的质量与研究者对问题的正确理解。生命并没有主动为我们整理好全部资料。它留下的，可能只是散落的字符。一些重叠。一些缺失。以及许多重复。计算机科学家的任务，是找到一种合适的数学语言，将这些零散证据组织起来。有时候，解决一道复杂问题的最好方法，并不是直接处理原来的碎片。而是先改变碎片的表示方式，让隐藏在其中的结构显现出来。

历史档案：基因组组装算法的发展

年份	真实事件	意义
1736年	欧拉研究柯尼斯堡七桥问题	图论与欧拉路径研究的重要起点
1946年	德布鲁因发表《A Combinatorial Problem》	发展循环短词组合结构
1977年2月24日	桑格等发表φX174噬菌体DNA序列	完整DNA基因组测序的早期里程碑
1995年夏	Idury与Waterman发表DNA组装新算法	将短词方法与片段组装联系起来
1995年7月28日	流感嗜血杆菌全基因组随机测序与组装成果发表	完成一个自由生活细菌的全基因组序列
2001年8月14日	Pevzner、Tang、Waterman发表欧拉路径组装研究	发展EULER与欧拉超路径组装框架
2008年3月18日	Zerbino与Birney发表Velvet研究	将de Bruijn图用于高通量短读段组装与图纠错
2012年	Bankevich等发表SPAdes	拓展复杂覆盖和单细胞测序条件下的图组装
2016年	Heng Li发表Minimap与Miniasm	推动长读段快速映射与组装

主要参考文献与历史依据

1977年 · 第一项完整DNA基因组测序成果Sanger等 — Nucleotide Sequence of Bacteriophage φX174 DNA《Nature》第265卷，第687—695页。用于核查φX174完整DNA测序成果、论文日期及研究团队。

1946年 · 德布鲁因的组合数学论文N. G. de Bruijn — A Combinatorial Problem荷兰皇家科学院相关论文，第49卷第7期，第758—764页。帮助区分1946年的组合数学工作与后来的DNA组装应用。

1995年 · Idury–Waterman序列组装算法A New Algorithm for DNA Sequence Assembly《Journal of Computational Biology》第2卷第2期，第291—306页。用于核查2001年EULER之前已有的重要短词组装研究。

1995年 · 完整细菌基因组的全基因组随机测序Fleischmann等 — Whole-Genome Random Sequencing and Assembly of Haemophilus influenzae Rd《Science》第269卷第5223期，第496—512页。核查1,830,137 bp的完整基因组成果及其全基因组随机测序背景。

2001年 · 本章核心一手论文Pevzner、Tang、Waterman — An Eulerian Path Approach to DNA Fragment Assembly《PNAS》第98卷第17期，第9748—9753页。重点阅读论文关于重复区域、Eulerian Superpath、读段路径约束与测序错误修正的论述。

2008年 · Velvet短读段组装Zerbino & Birney — Velvet: Algorithms for de Novo Short Read Assembly Using de Bruijn Graphs《Genome Research》第18卷，第821—829页。介绍de Bruijn图、错误分叉处理、Unitig压缩和重复区域解析。

2012年 · SPAdesBankevich等 — SPAdes: A New Genome Assembly Algorithm and Its Applications to Single-Cell Sequencing《Journal of Computational Biology》第19卷第5期，第455—477页。进一步说明复杂覆盖、图构造和纠错为什么需要多阶段算法。

2016年 · 长读段组装Heng Li — Minimap and Miniasm: Fast Mapping and de Novo Assembly for Noisy Long Sequences《Bioinformatics》第32卷第14期，第2103—2110页。用于核查长读段如何推动重叠图组装方法的发展。

**史实与教学重构说明：**1977年的φX174研究、1995年的流感嗜血杆菌组装、1995年Idury–Waterman工作、2001年EULER论文及后续Velvet、SPAdes等算法，均依据对应的论文资料。开篇实验室叙述属于对历史研究情境的文学化组织，没有虚构具体人物之间的对话或未被记录的实验细节。本章使用的ATGCGAT、GATTACA、AAAACAAC等DNA，全部为教学构造。互动实验与Python程序假设已知完整k-mer出现次数，不模拟真实测序覆盖、反向互补链、质量分数、插入删除错误或复杂读段路径约束。它们能够验证图论上的欧拉组装，但不能代替真实基因组组装软件。尤其需要明确：2001年EULER的贡献不能缩减为“使用Hierholzer算法遍历k-mer图”，更不能宣称任何重复序列都能被欧拉路径唯一解决。 原始研究还包含错误修正、读段路径保留和欧拉超路径等重要内容。

下一章预告 · 第55章 · 第五篇收官

《生命的家谱：算法怎样重建进化树》

1987年，美国得克萨斯州休斯敦。分子进化研究者正在面对一项越来越重要的问题。科学家已经能够比较不同物种的DNA与蛋白质序列。他们知道哪些位置相同。哪些位置发生过替换。也能估计不同序列之间的进化距离。但是，仅仅知道两种生物之间有多相似，并不意味着我们已经知道它们在生命演化历史中的分支关系。如果有十种生物。究竟哪两种应该首先连接？哪个共同祖先节点应该放在哪里？怎样避免把偶然相似误认为近缘关系？两位研究者提出了一种新的方法。他们是：斋藤成也（Naruya Saitou）。根井正利（Masatoshi Nei）。1987年7月，他们在《Molecular Biology and Evolution》发表：The Neighbor-Joining Method: A New Method for Reconstructing Phylogenetic Trees《邻接法：一种重建系统发育树的新方法》。论文第4卷第4期，第406—425页。它利用物种之间的进化距离，逐步寻找适合合并的邻居，构造进化树。 ^[PubMed]而在此之前，1981年，Joseph Felsenstein已经发表了：Evolutionary Trees from DNA Sequences: A Maximum Likelihood Approach《利用DNA序列构造进化树：最大似然方法》。这项研究展示了如何使用概率模型，对不同进化树与DNA序列之间的关系进行统计计算。 ^[PubMed]下一章，我们将追踪系统发育树算法的真实历史，并用一组人工DNA序列，亲手完成：从序列比对到距离矩阵。从距离矩阵到Neighbor-Joining。从候选进化树到最大似然评分。最后，讨论为什么一棵根据有限DNA数据计算出的树，是对生命演化历史的统计推断，而不是直接拍摄到的祖先关系照片。这一章，我们试着恢复一条被打碎的DNA；下一章，我们将尝试从多条DNA之间的差异，重建它们可能拥有的共同进化历史。

——《算法的故事》第二季 · 第五篇《破解生命的密码》 · 第54章完