本文首发于预印本平台：arXiv
摘要
目前可用的众多基因组展现了在大小、组成和结构方面的多样性，随着近年来多个全球生物多样性基因组计划的启动，未来还将有更多基因组问世。然而，即使有了近年来的技术进步，基因组测序仍然面临技术（如样本体积小、样本污染或缺乏合适的测序平台）和生物学（如生殖系限制性DNA、倍性水平变化、性染色体或超大基因组）等方面的挑战。近年来，基于k-mer的技术已经成为解决这些挑战的流行方法。这些技术基于将待分析的序列（如原始读长或基因组）分割为一组长度为k的子序列，称为k-mer。尽管这种方法看似简单，基于k-mer的分析可以快速直观地评估复杂的测序数据集。本文提供了k-mer在生物多样性基因组学中的理论性质和实际应用的全面综述，可作为参考手册使用。
背景介绍
过去的四分之一世纪里，基因组学领域取得了长足的进展。曾经，测序和组装一个部分基因组都是一个耗资数十亿美元的成就。如今，越来越多的物种已经完成了从端粒到端粒（T2T）的全基因组组装，全球范围内捕捉生命基因组多样性的重大项目也在稳步推进。然而，截至目前，绝大多数已知物种尚未拥有基因组组装，甚至没有相关的测序数据。实际上，根据“Genomes on a Tree tracker”追踪，仅有1%的已知物种拥有相应的基因组组装。
然而，基因组组装是一个非常复杂的过程，需要借助无组装方法的引导。这些方法通常基于一个流行的生物信息学概念——“k-mers”。k-mers已被证明是理解大量且持续增长的测序数据的高效且强大的工具。
每位基因组学家在某种程度上都会接触到k-mers。尽管这些名称多年来发生了很大变化（见框1），但其基本概念保持不变。例如，k-mers被广泛应用于许多生物信息学比对工具的后端算法中，如BLAST、bowtie2和minimap2，以及几乎所有的基因组组装软件。最近，k-mers作为一种分析基因组数据的流行框架得到了广泛应用，尤其是在基因组概况分析中，这得益于像GenomeScope这样用户友好的工具的开发。因此，k-mers现在已成为基因组学中的一个基础组成部分，这也激发了我们编写此指南的动机。我们定义并解释了k-mers的基本属性，回顾了基因组学中已确立的基于k-mer的技术，并展示了k-mers在解决复杂基因组问题中的最新创新应用。我们还提供了在线材料，读者可以通过真实的生物学示例来实践这些知识，在线材料可通过以下网址获取：
https://github.com/KamilSJaron/oh-know/wiki
知识卡片1：k-mers的历史背景
这一概念最早可以追溯到传奇人物克劳德 香农（Claude Shannon）的工作，他使用“N-gram”发展了通讯理论，后来用于计算自然语言的熵。在以数学为导向的研究群体中，这一概念通常被称为“k-tuples”，包括一些缩略的形式，如“ktup”或其他前缀如“L-tuple”或“ℓ-tuple”。为了吸引更广泛的受众，一些作者选择使用“k-word”或直接称为“word”。不过，最终“k-mer”成为了最流行且最常见的表达方式。杂交测序中曾使用“11-mers”作为寡核苷酸模板，尽管他们并未使用带k的通用形式。尽管k-mers的概念也在BLAST的文献中以“w-mers”出现，但k-mer成为更常用的术语花费了将近十年时间。在1990年代末期，k-mer的使用逐渐增多，包括在1999年具有里程碑意义的全基因组比对工具MUMmer的工作中。在2000年，Liu和Singh创造了“k-mer词频分布”这个术语，并将其描述为序列的“特征”。几年前，Mullikin和Ning发表了“词频图”，这是首次记录的k-mer频谱。在接下来的几年中，使用k-mer的出版物数量相比其他术语显著增加。k-mer这个表达随着一系列基因组组装软件、reads比对工具和专门用于计算k-mers的软件（如Jellyfish）的发布，在2000年代逐渐成为描述这一概念的主要方式。更多细节可参考补充文本1。
k-mer 基础
在基因组学背景下，k-mers是指包含在序列中的长度为k的核苷酸子串（例如单个read、参考基因组或任何其他序列）。k-mers通常用于DNA，但这一概念同样可以应用于RNA或蛋白质序列。任何基因组序列都可以分解为若干连续的k-mers，k-mers的数量取决于序列的长度(L)和k-mer的长度(k)。例如，在以下序列中：AAGTCCAT（L=8），有七个长度为2的k-mers（2-mer）、六个3-mer、五个4-mer、四个5-mer、三个6-mer和两个7-mer（见图1）。一个长度为L的序列中的k-mers数量等于L - k + 1。这是一个普遍原则，适用于任何序列，无论序列的长度或组成如何。
图1：k-mer分解的简单示例
该示例序列被分解为一系列k-mer，k值范围从2到7。注意，read中的k-mers数量为L - k + 1。
进一步将测序读长或基因组进行分割的意义何在？从最基本的层面来看，它们提供了一种方便的方式将基因组序列分解为简单的“词”。与自然语言不同，基因组没有空格或其他标点符号来标记词的开始或结束。然而，k-mers为基因组强加了这种结构，并产生了令人惊讶的强大效果。其一大好处是提高信噪比。在读长中出现一个测序错误意味着该读长与测序模板不完全匹配，但只有部分k-mers会继承这一错误。一个长度为L的读长包含L - k + 1个k-mers，但只有k个k-mers包含错误的序列（例如，对于100bp的读长使用21-mers进行分析时，一个测序错误会影响读长中的80个k-mers中的21个）。剩下的k-mers仍然代表真实的基因组序列。第二个好处则完全是计算层面的——在分析正确的k-mers时，我们可以使用完全匹配来进行k-mers计数，这比任何不完美匹配的算法（例如使用哈希表或二进制搜索，而不是计算完整的比对）要快得多。最值得注意的是，k-mers通过使用de Bruijn图来实现基因组组装，这种方法用于EULER或spades等组装软件。de Bruijn图是由读长中的k-mers构建的，随后进行一系列复杂的图转换。最后，第三个好处是，k-mers可以用于快速评估各种生物学特征。例如，它们可以用于估算基因组属性，如基因组大小、基因组重复水平和杂合性。它们还可以在不进行比对的情况下评估基因组之间的相似性。然而，在进一步扩展之前，我们需要了解k-mers的属性，以及k的选择如何影响数据集。
k-mers的基本属性
k的选择影响k-mers的两个重要属性：k-mer集合的复杂性和k-mer覆盖率。这里，复杂性指给定k时可能的不同k-mer的数量。对于四碱基字母表而言，其复杂性为4^k。然而，实际上我们通常不考虑序列的方向，因此反向互补序列通常与正向序列视为相同的序列，例如CAT和ATG将被一起计算。这些折叠的k-mer对称为“规范k-mers”，通常通过选择k-mer及其反向互补序列中按字典序排列较小的作为该对的代表序列（例如，ATG也表示CAT）。因此，当k为奇数时，可能的规范k-mers数量为4^k/2。在本文的其余部分中，除非另有说明，k-mers指的是规范k-mers。
k-mer空间的复杂性（即所有可能的k-mers）随着k的增加呈指数增长（参见补充图1）。然而，对于较小的k值，该空间将不足以代表或区分较长的序列。对于k=3时，使用四个碱基字母表时，只有32种可能的k-mers；对于k=7，有8192种。对于这些较小的k值，所有可能的k-mers在基因组中会多次出现，这些k-mers很可能在任何物种中都能被观察到。然而，这些短k-mers（k ≤ 11）的相对频率，通常被称为基因组特征，仍然携带着系统发育信号。
对于大多数应用，我们需要选择足够大的k值，使得基因组中的大多数k-mers不会是重复的。基因组中对应于唯一位置的k-mers比例随着k的增大而增加。对于人类基因组，k最小的值为11时，会有一些k-mers是唯一对应特定基因组位置的，尽管该空间中唯一位置的数量很少（在3x10^8个可能位置中少于2.1x10^6个k-mers）。相反，我们通过计算给定k-mer在长度为G的基因组中预期出现的次数来估计k的最小长度。其估计值为(G/4)^k，这意味着我们需要选择k至少为log4(G)。对于3Gb的人类基因组，建议的最小长度为16；然而，实际上人类基因组中有许多16bp的重复序列，因此选择更长的k-mers更为有利。
当k足够大以包含大量在基因组中唯一位置的k-mers时，我们可以估计平均有多少读长包含基因组中单拷贝的k-mer，这就是k-mer覆盖率Ck的定义。经典的测序覆盖率Cg定义为基因组中每个碱基的平均reads深度，
Cg = N * L / G
其中N是测序的reads数量，L是平均读长长度，G是基因组总大小。尽管平均基因组覆盖率和平均k-mer覆盖率并不相同，但它们之间有一个简单的关系可以近似计算它们的相对值：
Ck = Cg(L - k + 1) / L
因此，对于给定的reads集合，较小的k值会提供较高的k-mer覆盖率，不同k值的k-mer覆盖率差异在短读长中会更加明显（见补充图2）。例如，使用k=100对于100bp的reads来说存在问题，因为一个read只会表示一个k-mer，但使用k=51则每个read可以产生50个k-mers。
最后，测序过程的准确性有限，且随测序平台不同而有所变化，测序的核苷酸中有一部分会代表测序错误。这些错误会从reads传播到k-mers，并且受影响的k-mers比例同样取决于所选的k值。假设错误均匀分布（即每个read中的每个测序位置出错的概率相同），则k-mer表示真实基因组序列的概率为：
因此，对于给定的读段集合，较小的k值会提供较高的k-mer覆盖率，不同k值的k-mer覆盖率差异在短读段中会更加明显（见补充图2）。例如，使用k=100对于100bp的读段来说存在问题，因为一个读段只会表示一个k-mer，但使用k=51则每个读段可以产生50个k-mers。
最后，测序过程的准确性有限，且随测序平台不同而有所变化，测序的核苷酸中有一部分会代表测序错误。这些错误会从读段传播到k-mers，并且受影响的k-mers比例同样取决于所选的k值。假设错误均匀分布（即每个读段中的每个测序位置出错的概率相同），则k-mer表示真实基因组序列的概率为：
Pr(error free k-mer}) = (1 - e)^k
其中，e是每个碱基的测序错误率。因此，选择较小的k值时，会有更大比例的k-mers表示真实的基因组序列（见补充图2），从而在后续分析中更有用。
总之，k值的选择是捕获基因组复杂性、测序覆盖率和错误率之间的权衡。因此，“合适的”k值取决于具体的应用、测序深度、测序平台以及所测序基因组的类型。例如，长读长测序技术允许使用更大的k值，因为对于较大的k值，k-mer覆盖率不会显著下降（见上面的公式），因此在某些应用中可以设定k > 1000。然而，测序错误率仍然对k值的实际选择构成限制。对于短读长数据集，通常设置k在21到31碱基之间的范围，因为它能从几乎任何生物的已测序基因组的独特位置生成大量k-mers，同时提供高k-mer覆盖率。此外，设定k ≤ 31在计算上比较长的k值更为高效，因为这些序列可以在64位以下表示，使得计算机可以更快速地比较它们（即，不必逐字符比较，而是将其视为64位整数进行单次计算指令的比较）。最后，对于某些分析，选择奇数k值也有优势，因为正向和反向的k-mers将具有不同的序列，例如，AT的反向互补仍然是AT，但ACT的反向互补是AGT。
使用k-mer频谱对测序文库进行分析
k-mers在基因组学中最常见的直接应用之一是通过k-mer覆盖分布来表征测序数据集，这通常被称为k-mer频谱或k-mer直方图。一个中度杂合的二倍体生物的典型k-mer频谱通常有四个明显的覆盖峰（图2）：第一个峰代表测序错误（低覆盖率，图中为粉红色）；第二个峰代表来自杂合位点的独特基因组序列（1n，黄色）。该峰值将围绕Ck覆盖率，有时称为1n或单倍体k-mer覆盖率；第三个峰代表基因组中所有的纯合位点，围绕2*Ck覆盖率（2n，蓝色）；第四个峰，通常较小，代表基因组重复区域，围绕4*Ck覆盖率（4n，橙色）。
图2：从测序reads到k-mer频谱
此图展示了一个基本的示例，说明如何将基因组reads转化为可以计数的k-mers，并通过k-mer频谱表示。A) 两个单倍型中的基因组表示，带有一个重复区域（橙色）和六个杂合位点（黄色）。 B) 测序reads及其与基因组的比对，包含了基因组序列，同时也有测序错误（粉红色）的叠加。C) 40个碱基长的测序reads的k-mer（8-mer）分解展示。D) k-mer频谱表示在测序的基因组中，多少不同的k-mers（y轴）具有特定的覆盖率（x轴）。在绝大多数实际基因组数据集中，k-mer频谱中的峰值会重叠，并且还会有少量k-mers表示其他倍性（例如，杂合重复的3n）。
k-mer频谱分析已经成为基因组测序的标准步骤，生物信息学界也认识到了对高效工具的需求。因此，k-mer计数工具的性能得到了显著提升，并且在k-mer工具包中开发了许多附加功能（常用k-mer计数工具的概览见补充表1）。
然而，计算k-mer频谱只是第一步。可以通过拟合各种模型来估算基因组特征，例如基因组大小、杂合性或重复性，这些统称为基因组调查。
拟合基因组分析模型
最早的基因组分析技术之一是用于估算k-mer覆盖率和基因组大小。基因组大小估算的原理与计算k-mer覆盖率（公式1和2）几乎相同。这是因为基因组大小（G）是k-mer总数除以k-mer覆盖率（Ck）和倍性水平（p）。
G = N (L - k + 1) /(Ck * p)
其中，N是数据集中的reads总数，L是read长度，与k-mer大小k一起给出数据集中的k-mer总数。该公式忽略了测序错误，这些错误会通过将测序错误引入的k-mers计入基因组的一部分，进而人为增加估算的基因组大小。包含错误模型的工具可以减去估算的错误k-mers，从而允许更精确的基因组估算。
k-mer覆盖率（Ck）的估算是通过拟合一个或多个分布到经验计算的k-mer频谱上完成的。最简单的分布，也是基因组测序等采样过程的自然选择，是泊松分布。泊松分布易于拟合，因为它由单一参数决定，能够同时确定均值和方差。然而，在实践中，观察到的方差通常大于泊松分布拟合出的方差（即观察到的分布是“过度离散”的），因此更合适的做法是拟合负二项分布。后者是泊松分布的推广，具有两个参数：一个均值参数和一个独立的方差参数，但形状相似。对于高覆盖率值，可以使用正态或高斯分布，但对于低覆盖率（低于10x覆盖）则不适用。这是因为高斯分布在均值周围是对称的，所以在低覆盖率情况下，它会估算出小于0的覆盖率值（负值），这是不合理的。另一种分布——特别针对某些位点比其他位点具有较低测序概率的事实（例如异染色质）——是偏态正态分布。findGSE的作者提出了偏态正态分布来解决这个问题。没有一种分布适合所有情况，但通常来说，越复杂的覆盖模型越适合高覆盖率数据集，而对于低覆盖率数据，最好使用参数最少的分布。
大多数模型拟合方法期望具有相对较高的测序覆盖率（>10x），以便可以在几乎没有假设的情况下拟合分布参数。RESPECT的作者采用了不同的方法，专为低覆盖率（表面测序）数据设计，专门针对0.5-2x测序覆盖率的数据集进行优化。该模型估算了基因组中分别以一次、两次、三次等次数出现的k-mers数量，同时使用简单的错误和覆盖率模型。为了使模型拟合成为可能，作者对一些参数进行了约束；例如，错误被建模为在reads中仅观察到一次的k-mers的比例，但基因组中k-mers的个别计数也受限于几百个公开基因组中经验观察到的比率。这些经验约束是通过二倍体基因组的单倍体表示生成的，因此我们预期该方法仅对单倍体物种或具有低杂合性的二倍体基因组提供可靠的估算。
k-mer频谱还可以揭示基因组的其他特性。通过分析1n与2n峰的k-mer数量，可以确定全基因组的杂合率：在低杂合率情况下，1n峰相对较小，因为大多数k-mers将是纯合的，但随着杂合k-mers的增加，1n峰会变高。实际上，杂合率只需要相对适中的水平就足以使1n峰与2n峰相匹配，因为每个杂合变异会导致2*k个杂合k-mers（假设这些变异适当地间隔开）。在这种情况下，当k=21时，仅1.19%的杂合率就足以使1n峰与2n峰一样高。
三种基因组分析方法也包括对全基因组杂合率的估算。首先，在整体杂合率较低（<<0.5%）的基因组中，杂合位点通常距离其他杂合位点超过k个核苷酸。因此，每个杂合位点将生成2*k个杂合k-mers，这已被提出作为一种简便的杂合率估算方法。然而，许多异交物种、杂交种和其他物种的杂合率明显较高。GenomeScope通过解决相邻变异的问题（即距离小于k bp的多个杂合变异）来应对这一挑战。假设杂合位点和重复区域是独立且均匀分布在基因组中的，杂合k-mers和纯合k-mers的比例可以通过每个核苷酸杂合率估算来计算。第三种方法，Tetmer，基于无限位点模型和溯祖理论来估算种群的遗传多样性。有关这三种方法的详细比较，请参阅补充文本2。
最受欢迎的工具是GenomeScope 2.0，它还支持多倍体基因组和更高级的模型拟合方法。我们将用它来展示拟合基因组模型的各种用途和问题，但我们也鼓励读者考虑其他基因组拟合工具，因为它们可能更适合解决特定问题（参见补充表2中的概述）。
k-mer频谱的常见特征
为了生成高质量的二倍体基因组参考序列，建议至少进行25-30×覆盖的长读长测序，或更普遍的是每个单倍体15×覆盖。即便是简单的k-mer频谱目视检查也有助于快速评估是否达到了这一覆盖率。这样的覆盖率应产生一个显示明显覆盖峰值的k-mer频谱，如欧洲槲寄生（Viscum album）所示（图3A）。如果测序数据的覆盖率不足，峰值将定义不清，因为纯合和杂合基因组峰将在覆盖图的左侧混合在一起。如果峰值仍然可见，则可能仍能拟合出有意义的基因组模型，例如小龙虾（Procambarus virginalis）的情况下（图3B；数据来自【51】）。
通常，关于我们测序的物种已有一些已知信息，特别是通过细胞遗传学技术评估的倍性或基因组大小。将先前的知识与k-mer频谱得出的估算值进行对比，往往有助于识别数据中的潜在问题。例如，在小龙虾的案例中，k-mer频谱的基因组大小估算值与流式细胞仪结果几乎完全一致，这表明模型收敛良好。非常低的测序覆盖率或较高的错误率会导致峰值混合；基因组k-mers与错误k-mers变得难以区分。这可以在细香葱（Allium schoenoprasum）的案例中看到（图3C），其中模型（黑线）与数据（蓝色直方图）拟合得不佳。在这种情况下，估算值仅是模型收敛不佳的产物。预测的基因组大小远低于我们对葱属物种的预期值，其他物种的基因组大小在8.4-13.4Gbp之间，而覆盖率也远高于我们对这种频谱形状的预期。覆盖问题通常可以通过增加测序量解决，而高错误率可能需要更换测序技术或文库制备方法。
图3：k-mer频谱示例
A. Viscum album：二倍体频谱，数据量足以观察到两个明显的峰值，并拟合出一个准确反映基因组特征的模型，尽管该基因组体积较大。B. Procambarus virginalis：低覆盖率样本的k-mer频谱，几乎不足以进行模型拟合。值得注意的是，我们使用了k = 17来增加k-mer覆盖率，使模型拟合成为可能。 C. Allium schoenoprasum：该数据集的测序覆盖率大约为1x，错误k-mers与基因组k-mers完全混合，因此模型未能收敛到有意义的估算结果。D. Hypsibius dujardini：一个严重受污染的缓步动物样本。
污染是另一个常见原因，导致我们在k-mer频谱中观察到峰值混合，特别是对于那些难以在无菌环境中培养，甚至完全无法培养的小物种。在这些情况下，除了目标物种之外，还会同时测序其他物种的基因组序列，导致基因组剖面受到污染。例如，在缓步动物Hypsibius dujardini中，真正的基因组峰被多个共生体基因组叠加（图3D）。如果污染仅来自少数几种物种，并且覆盖率很高，我们可能会看到对应于个体来源的单独峰值。然而，这些峰值通常间隔不均匀（补充图3）。与低覆盖率数据集类似，基因组模型报告的数值很可能是不准确的。
不同的测序数据集在测序深度上显示出不同的变化，但一般规则是覆盖率越高，峰值分离越好。有时，尽管覆盖率相对较高，我们仍会观察到峰值混合。这可能由各种生物学或技术原因引起。例如，许多硬骨鱼类的测序运行显示k-mer频谱中存在一个“桥梁”连接了错误k-mers和1n峰（补充图4）。这一模式被归因于许多鱼类谱系中串联重复比例较高，以及PacBio HiFi测序在低复杂度A/G富集区域中的覆盖率下降。另一个例子是Illumina短读长测序众所周知的GC含量偏差。这种偏差因化学试剂的不同而变化，并且在无PCR策略中有所减少，但最显著的还是由物种的生物学特性引起。
上述偏差在使用相似样本和测序技术时得到了很好的描述和重复。然而，测序运行也受到样本处理、使用的防腐剂或偶发的测序流动槽制造问题的影响。即使是这些问题，有时也能在k-mer频谱中显现出来。例如，基于k-mer的海角蜜蜂样本的基因组大小估算值（142Mbp，补充图5）远小于已发表的基因组大小（236Mbp）。这种情况仅出现在三个样本中的一个，来自同一项目（和同一群体）的另外两个样本显示出预期的基因组大小。该异常样本的频谱还显示出更高的峰值混合，整体表明这种覆盖率下降不是由被测物种的生物学特性引起的，而是技术问题导致的。
总的来说，缺乏明确的、均匀间隔的峰值表明没有单一来源具有足够的覆盖率来生成预期的频谱。牢记目标物种的生物学特性，许多模式和潜在问题是可以预见的。
拟合模型时的常见陷阱
基因组模型拟合的质量在很大程度上依赖于数据的质量和覆盖率，但也与基因组的生物学特性有关。基因组模拟最常见的问题是单倍体（1n）k-mer覆盖率与“错误的基因组峰”收敛。这通常发生在1n覆盖率峰不够明显或显著高于二倍体峰时。这可能是由于基因组的杂合性极低（即1n信号非常弱，见补充图6A），或者因为覆盖率非常低，导致1n峰与错误峰大部分重叠（图3C）。
当1n覆盖率未正确拟合时，估算值将不再具有任何关于基因组的生物学信息。因此，必须对拟合进行目视检查，并确保估算值与已知的生物学信息相符。例如，如果我们测序一个二倍体自交植物，但估算的杂合率超过5%，极有可能真实的1n覆盖率约为估算值的一半。自交生物的预期遗传多样性较低，可能会让基因组分析工具忽略第一个峰。在GenomeScope中，用户可以添加标志“-l
另一个常见问题是，许多k-mer计数工具会在用户指定的值处停止计数，典型的默认值是10,000。k-mer直方图中的最后一个值通常表示k-mer计数覆盖率最高的k-mers或更高覆盖率的k-mers的数量。这会导致基因组中最常见重复序列的频率分辨率降低。对于许多基因组来说，这影响不大，但对于极度重复的基因组来说，基因组大小的估算可能会大大偏低。当使用默认值时，大理石纹龙虾的估算基因组大小几乎是实际大小的一半（补充图7）。如果观察到基因组大小出乎意料地小，通常是由于k-mer频谱被截断所致。值得注意的是，FastK并不会显式计算每个高覆盖率k-mer的覆盖率，取而代之的是，它通过报告最后一个值来表示所有更高覆盖率的k-mers总覆盖率，确保了更准确的基因组大小估算。
k-mer频谱和组装基因组的联合解释
在组装之后，可以使用从k-mer模型估算的基因组质量来评估组装的准确性。一个简单但信息丰富的指标是将估算的基因组大小与组装大小进行比较。这种简单的比较非常有用，因为它可以快速指示组装错误，尤其是由于杂合性引起的部分重复。在达尔文生命之树项目中测序的冬青（Ilex aquifolium）就是基因组大小估算与组装大小不匹配的一个例子。根据PacBio读长的k-mer直方图估算，I. aquifolium的基因组大小为815.6 Mbp。然而，初步组装稍大（830.1 Mbp；见补充图8A），并且该基因组中有7.2%的重复BUSCO基因，这表明该组装可能存在一些未折叠区域，导致基因数量和（单倍体）基因组大小膨胀。这可以通过后期去修正，在这种情况下使用了purge_dups工具。
现在有几种工具可以使用预组装的k-mers评估基因组组装的准确性，我们将特别讨论Merqury，但这并不是唯一的工具。这些工具的基本原理是，将测序reads中的k-mers与最终组装基因组中的k-mers进行比较，以确定组装的完整性。这是基于这样一个假设：最终的组装应该包含测序reads中大多数k-mers，排除那些可能由于测序错误导致的低覆盖率k-mers，以及二倍体基因组中被折叠组装时的一半杂合等位基因（见补充图8B中Merqury绘制的冬青基因组示例）。此外，Merqury估算了一致性质量值评分，假设最终组装中的所有k-mers应该存在于reads集中。然后它估算特定碱基为错误的概率，这可以以常用的Phred评分来报告。这种方法对于没有参考基因组或参考基因组组装不好的物种特别有用。然而，它确实需要高准确率的reads（如Illumina或PacBio HiFi）作为参考点，并且理想情况下，这些参考reads应该与用于组装的reads是正交的(orthogonal)。
使用k-mers对测序文库进行比较
k-mers 还可用于识别两个测序文库之间的重大基因组差异。这种分析方法可以帮助解决多种研究问题，最显著的是通过不同性别生成的文库识别性染色体。然而，它也可以用于识别组织或个体特异性染色体，如B染色体。k-mer 染色体识别技术可以独立使用，也可以与现有的基于覆盖率差异或两种测序文库之间序列组成差异（如SNP差异）的方法结合使用。k-mer 技术相较于其他方法的一个优势是对高质量参考基因组的依赖性较低。因此，k-mer技术可能是适用于参考基因组破碎或不存在的非模式生物的更好方法。
在使用k-mers 识别特定染色体的分析中，基本方法是比较两个染色体组成不同的测序文库中的k-mer频率。生成一个二维k-mer频谱，比较两个文库中k-mer的频率，可以识别出在两个样本中频率不同的染色体所属的k-mers。例如，如果在具有XO性别决定系统的物种中识别性染色体，X染色体的k-mers在雄性样本中的频率将是雌性样本的一半，而常染色体的k-mers频率则相同。X染色体所属的k-mers随后可以被分离，并映射到组装体上以识别属于该染色体的scaffolds，或映射到reads上以识别属于X染色体的reads。
图4：多个测序文库的k-mer比较
此图展示了两种物种的二维k-mer频谱示例，分别为：A. XO性别决定系统；B. XY性别决定系统。
两个图中的热图显示了雄性样本（x轴）与雌性样本（y轴）中的k-mers频率。两者中常染色体k-mers与X染色体k-mers可以区分，因为X染色体k-mers（绿色框）在雄性样本中的频率为雌性样本的一半，而常染色体k-mers（橙色框）则在两者中的频率相同。类似地，Y染色体k-mers可以在图B中识别出来，这些k-mers在雌性中缺失，而在雄性中存在（粉色框）。然而，需要注意的是，覆盖雄性常染色体杂合位点的k-mers可能被误认为属于性染色体的k-mers。
使用k-mers识别特定染色体在感兴趣的区域与其他基因组区域差异较大时最为有效。例如，这种方法难以用于识别XY性染色体系统中同型（homomorphic）Y染色体，因为X和Y染色体之间的分化较低，少数k-mers仅限于Y染色体。因此，这种情况可能需要更宽松的识别阈值。样本的杂合性也是一个重要的考虑因素。例如，如果从自然群体中收集的雄性和雌性样本中识别性染色体，覆盖杂合SNP的k-mers会表现出性染色体特异性的k-mers特征，因为它们的频率也将是常染色体纯合k-mers的一半。通过测序更多个体，可以区分个体特异性的杂合k-mers和性别特异性的k-mers，部分克服这一问题。最后，重要的是考虑每个样本的测序深度。测序覆盖率需要足够高，才能区分测序样本中不同倍性水平染色体的k-mer频谱峰值，并区分代表测序错误的k-mers，类似于基因组分析。
k-mers的创新使用案例
到目前为止，我们已经讨论了k-mers的一般特性以及如何将已建立的k-mer方法应用于测序数据。然而，生物学中有许多问题很难用现有工具解决，这迫使我们设计针对我们拥有的问题和数据开发出创新的方法。这里，我们展示了k-mers在三个不同背景下的三个问题中的创新应用：生殖细胞DNA的研究、异源四倍体的亚基因组，以及物种分类。
提取生殖细胞限制性染色体
Hodson等人通过比较文库的方法，识别了黑翅蕈[xùn]蚊（Bradysia coprophila）中属于生殖细胞限制性染色体的序列，以探究这些特殊染色体的进化起源。在该物种中，雄性表现出一系列特殊的染色体淘汰事件，这些事件在不同组织和生命阶段产生不同的染色体核型（最近的综述见【68】）。具体而言，精子携带两个X染色体、两个生殖细胞限制性染色体和一套常染色体，而体细胞只有一个X染色体、没有生殖细胞限制性染色体（因此得名），以及两套常染色体。在这种情况下，主要由精子组成的生殖组织与体细胞组织相比，所有染色体类型的频率都不同。为了表征该物种的生殖细胞限制性DNA，大约95个未交配雄性的睾丸被合并成一个文库，并与同一个体的头部分别测序，以生成一个纯净的体细胞文库，与生殖细胞文库进行比较。
通过2D k-mer频谱比较两个文库表明，确实有大量k-mers仅出现在睾丸文库中，属于生殖细胞限制性染色体。然而，意外的是，睾丸中常染色体k-mers的覆盖率仍高于X染色体k-mers的覆盖率，尽管常染色体在精子中的频率较低（预期在精子中的X染色体:常染色体k-mers比率为2:1）。这一现象主要是由于生殖细胞文库被体细胞污染。粗略估计表明，大约77%的生殖细胞测序文库由体细胞组成，导致X染色体的覆盖率低于常染色体，但并没有达到预期的1/2覆盖率（如果测序文库仅为体细胞组织）。因此，在这种情况下，两个组织类型中的三种染色体类型频率不同，有助于确定每种组织类型在测序文库中的相对组成。
染色体组特异性k-mers与contigs匹配，结果几乎完美地将contigs分配到染色体组（生殖细胞限制性染色体、X染色体、常染色体）。在该分析中，分配的质量特别高，原因有两个：首先，生殖细胞限制染色体与X染色体和常染色体（体细胞染色体）的同源区域差异很大；其次，该蕈蚊谱系的杂合性极低，因为它是在1910年代由Charles W. Metz分离的。基于k-mers的染色体排序方法能够在没有高质量参考基因组（当时）的情况下，准确识别出非模式物种的组织特异性染色体，促进了下游对生殖细胞限制性染色体起源的分析。
异源四倍体的亚基因组
多倍体基因组由于其庞大的规模、基因组震荡引起的重复内容，以及组成相似的染色体（同源染色体，即由物种分化产生并通过异源多倍体化重新结合在同一基因组中的染色体）的存在，组装难度较大。尽管存在这些挑战，已有几种异源多倍体的染色体级别组装被发表，其中亚基因组（即多倍体化事件后重新结合的基因组）得到了重建。在这一背景下，k-mer方法被用于分离亚基因组，推动了我们对基因组进化的理解，包括亚基因组的进化，如基因保留偏差、基因功能、自然选择和同线性等。
亚基因组的分离通过将一个谱系的进化历史划分为三个阶段来实现（补充图10）：
连续相似性矩阵（SSM）方法的示意图概览
(A) 遗传原理：在从二倍体祖先基因组到异源多倍体基因组的进化过程中，长末端重复序列（LTR）的分化，突出显示三个关键阶段。(B) 计算流程：全基因组长末端重复序列（LTRs）的识别，相似性矩阵的构建，基于R的聚类和可视化推断亚基因组。
该图来自Deciphering octoploid strawberry evolution with serial LTR similarity matrices for subgenome partition. Haomin Lyu et al. bioRxiv 2024.07.31.606053; doi: https://doi.org/10.1101/2024.07.31.606053
现在，考虑每个阶段中的TE积累情况：
阶段1：在这个阶段积累的TE预计在两个亚基因组中均匀分布；
阶段2：在分化事件之后、多倍化之前积累的TE将在两个亚基因组上表现出差异，因为两个独立的谱系会积累各自独特的TE；
阶段3：多倍体化事件后积累的TE将在两个亚基因组之间大致均匀分布。
在这一背景下，k-mers作为亚基因组分离的有力工具，尤其是在染色体水平组装中。由于进化信号往往会相互覆盖，序列的碎片化提供了一种有效的方式来解开过去发生的三个不同阶段中的过程（补充图10）。亚基因组的分离涉及两个步骤：首先，需要识别同源染色体（homeologs）。这可以通过UCEs/COS（超保守元件、保守直系同源序列）或基于比对的方法及同线性分析来实现。当同源染色体被识别后，k-mers可以用于识别TE在三个阶段中的积累情况。具体来说，通过进行k-mer频谱分析，寻找两个特征：首先，选择那些出现频率很高的k-mers（即覆盖率>100x），这些k-mers应该来自基因组的重复区域，如TE；其次，选择在同源染色体成员之间表现差异的k-mers（即在一个同源染色体中比另一个频率高两倍以上的k-mers）。高频k-mers指向TE，而第二个特征则挑选出表现差异的TE，从而有效地聚焦于第二阶段——即祖先基因组分离并积累各自TE的时期。通过基于这种k-mer分离的层次聚类，染色体可以按亚基因组进行分组。
使用短k-mers进行物种分类
在某些类型的分析中，使用非常短的k-mers（k<10）可能会有益。在此讨论的例子中，使用了靶向扩增子测序来分析平均仅为160bp的单倍型。分析的目标是通过将查询序列与参考序列集进行比较来识别物种。因此，分析的是重建单倍型中的k-mers，而不是直接分析reads。此外，由于单倍型可以通过引物进行定向，分析使用了完整的k-mer集，而不是规范k-mers。
选择k值的权衡在于序列变异的容忍度与捕获的复杂性之间的平衡。由于分析是基于重建的单倍型而非reads进行的，因此k-mer覆盖率在这一权衡中不发挥作用。较大的k值对查询序列与参考序列之间的变异容忍度较低，而较小的k值则很有可能在多个位置偶然找到相同的k-mer。例如，在149bp的序列中，5个均匀分布的SNP会导致没有任何25-mers与参考序列匹配。相反，在相同长度的序列中，所有4-mers都是唯一的几率极小（<10^-22）。基于这些权衡，我们选择8-mers作为合理的长度。在平均目标长度为160bp的情况下，单倍型内所有8-mers唯一的概率为84%。
为了进行物种分类，我们计算查询单倍型到参考序列集中每个单倍型的k-mer距离。k-mer距离量化了查询序列与参考序列之间匹配的k-mers的比例。最近邻序列是与查询单倍型的k-mer距离最小的参考单倍型。物种标签通过识别查询样本的所有扩增子目标的最近邻居并聚合它们的贡献进行分配。
该扩增子序列集和物种分类方法是为整个按蚊属（Anopheles）蚊子物种的分类而开发的。k-mers提供了一种客观的方法来比较高度分化的序列，而多序列比对或与单一参考基因组的比对往往会引入偏向参考序列集中代表性更好的类群。此外，k-mers还提供了一种自然的方式来结合小的插入缺失变异（indels）和SNPs，从而在处理小于10kb的序列时显著增强区分物种的能力。
k-mers在基因组学中的未来
正如词汇对语言学的重要性，k-mers对于生物信息学同样至关重要。它们是表示生物序列最基本的方法之一，然而它们的用途不容小觑，因为它们是支撑基因组学、转录组学和宏基因组学众多应用的基础数据类型。在这些成功的基础上，多个研究趋势已经出现，以使它们在若干重要应用中变得更加高效和有效。
在基因组学中，除了分析日益复杂的生物样本之外，k-mers现在被用于推动多个基因组组装和分析应用。一个强大的技术是“三亲分箱”（trio binning）的兴起，其中从父母基因组的未组装reads中识别出的k-mers用作标记，以对其子代的未组装reads进行单倍型分型。这使得在样本中进行全基因组分型和单倍型的de novo组装成为可能，从而在二倍体基因组组装中带来了新的进展，使其在连续性和准确性上超越了以往的方法。另一个强大的技术是“单一唯一核苷酸k-mers”（SUNKs）的应用，以帮助组装重复基因组和重复序列。SUNKs可以在未组装的reads中识别出来，作为“锚定”序列，以高度可信的方式将reads定位在基因组组装中。例如，在人类基因组中，SUNKs通过识别嵌入复杂重复阵列中的局部唯一序列，解决了片段重复的难题。另一个巧妙的应用是使用k-mers作为变异标记来推动基因型与表型关联研究。这在标记和分析结构变异时尤其强大，因为结构变异是从原始reads中最难研究的一类变异。展望未来，我们预计将通过相关的k-mer技术，在更复杂的基因组和泛基因组的组装和分析中取得进一步的进展。
同样，k-mers 在各种序列分类应用中也发挥着重要作用。在宏基因组学中，开创性算法Kraken利用k-mers作为个体物种的标记，使得无需比对即可实现快速且稳健的单个read分类。在转录组学中，开创性算法Sailfish展示了无需比对，通过使用k-mers作为单个转录本的标记就可以实现准确的转录本定量。我们预计未来k-mers将被广泛应用于作为分类和定量多样化样本的标记。同时，我们还预计未来将有更多应用会将k-mers嵌入抽象语义表示中，以用于机器学习应用，类似于word2vec及相关方法在自然语言处理中的核心地位。
最后，一个最重要的技术进展是采样和抽取技术的兴起，例如minimizers和minhash，以减少k-mers的计算复杂度。这些方法的核心思想是，对于许多应用而言，不必详尽考虑序列中的每一个可能的k-mer。相反，通常只需集中关注一个小的代表性子集即可进行分析，通常该子集只代表所有可能k-mers的几个百分点或更少。由于子集远小于完整列表，比较这些子集的速度快得多，并且所需的内存也大大减少。例如，minimizers最初是为减少基因组组装的计算需求而开发的，现在用于加速流行的比对工具如minimap2，通过将算法集中在一个较小的潜在种子列表上。同样，minhash推动了开创性的Mash算法，通过比较小型代表性k-mers列表，快速估算基因组对之间的相似性，使得能够在几分钟内用笔记本电脑比较成千上万个基因组。目前该领域的研究前沿集中在开发更高级的方案，用于选择或组合k-mers，以确保代表性k-mers子集既无偏差又能应对测序错误和重复元件的影响（详见【86】的综述）。还有正在进行的工作是开发高效的软件库和工具包，以快速高效地处理这些应用中的k-mers，特别是减少为大型k-mers集建立索引所需的内存，以及在大型基因组集合中建立k-mers索引。
总体而言，k-mers 的持续研究和应用前景非常光明。每当研究人员考虑对基因组进行分析、比较基因组或进行其他基于比对的分析时，他们都应思考是否可以仅使用k-mers来实现。
Cite as:
Katharine M. Jenike, Lucía Campos-Domínguez, Marilou Boddé, José Cerca, Christina N. Hodson, Michael C. Schatz, Kamil S. Jaron.
Guide to k-mer approaches for genomics across the tree of life
https://doi.org/10.48550/arXiv.2404.01519
