多倍体基因组Kmer Survey填坑指南：从峰型判读到参数设置

干过几年基因组组装的人，大概都体会过那种感觉：拿到一个自然界里普普通通的多倍体物种，按着二倍体的Survey流程跑一遍Kmer，出来的数字怎么看怎么不对劲——基因组大小要么翻倍，要么直接砍半，杂合度高得离谱，有的甚至直接给你报一个“Optimization failed”就罢工了。这时候如果经验不足，第一反应通常是怀疑测序数据有问题，转头去重新测序、重新建库，白白烧掉不少钱和时间。但实际上，问题大概率出在Kmer分析模型与多倍体基因组结构不匹配上。这正是复杂基因组Survey分析中最值得讲清楚、也最容易被忽略的事情。

这篇文章不打算从Kmer原理的ABC开始讲，而是直接围绕“多倍体+Survey”这个组合展开。我会把二倍体跟多倍体在Kmer分布上的本质差异、动手前要确认的信息、完整实操流程，以及我这些年踩过的坑一一说清楚。无论你是刚开始接触基因组Survey的研究生，还是已经跑过几次流程但始终被多倍体折磨的从业者，相信都能从里面找到可用的东西。

1. 先把底层逻辑说透：Survey到底在算什么

1.1 Kmer的“小卡片记账法”

搞懂多倍体为什么难，先得清楚Kmer Survey在干什么。测序得到的reads是一段一段的碱基序列，Kmer就是把每条reads切成固定长度为K的短片段。比如一条150bp的reads，以K=21为例，会得到150-21+1=130个21bp的小片段。这些Kmer记完数之后，我们画一张直方图：横轴是某个Kmer在数据里出现的次数（深度），纵轴是这个深度对应的Kmer种类数量。

这张直方图之所以能用来估算基因组大小，背后的逻辑很简单。假设一个基因组大小为G，均匀测序覆盖度为C，那么基因组上每个位置平均会被C条reads覆盖，对应到Kmer深度也大约是C（准确说是C乘以(L-K+1)/L这个边界修正系数）。整个基因组有多少种Kmer？理论上大约是G种（每个碱基位置产生一个Kmer，边界效应忽略不计）。总Kmer条数则是G×C。于是基因组大小就可以用“总Kmer条数除以主峰深度”估算出来。

举个例子，一个500Mb的基因组，50X测序，150bp读长，K=21，那么：

• 总Kmer条数 ≈ 500×10^6 × 50 × 130/150 ≈ 2.17×10^10
• 主峰深度 ≈ 50 × 130/150 ≈ 43.3
• 估算基因组大小 ≈ 2.17×10^10 / 43.3 ≈ 500Mb

如果基因组里有重复序列，部分Kmer的真实拷贝数>1，这些Kmer会出现在更高深度位置，形成一个尾巴；如果物种杂合，杂合位点附近的Kmer因为存在两种等位基因，每种出现的深度只有平均水平的一半，于是会在主峰左边一半深度处拱起一个小峰。二倍体Survey的核心就是解读这两个信号。

1.2 多倍体如何把Kmer分布彻底搞乱

多倍体的麻烦在于，基因组不再是“每个位点一份”，而是两套甚至更多套。这里要先区分两类多倍体，因为它们的Kmer行为完全不同。

异源多倍体（allopolyploid）是由不同物种杂交后染色体加倍形成的，比如小麦是六倍体，棉花有很多异源四倍体种。它的特点是不同亚基因组之间序列分化程度较高，每个亚基因组内部的Kmer模式接近于一个独立的二倍体。整体来看，异源多倍体的Kmer分布约等于多套二倍体分布的直接叠加，主峰位置不会偏移，但低频区会有额外的信号累积，导致二倍体模型拟合出的杂合度虚高。

同源多倍体（autopolyploid）就麻烦得多。它是同一物种基因组直接加倍，各套同源染色体之间序列相似度极高，但又不是完全一致。一个四倍体物种，某个位点有四份拷贝，如果这四份拷贝完全相同，那么这个位置的Kmer会被四份拷贝同时贡献reads，Kmer深度直接变成4λ。如果拷贝之间存在SNP差异，有的Kmer只在其中一份拷贝中存在（深度λ），有的在两份中存在（深度2λ），有的三份（3λ），有的四份（4λ）。最终画出来的Kmer分布是λ、2λ、3λ、4λ多个峰的叠加，形态像一串连在一起的“糖葫芦”，或者干脆糊成一片宽峰。

这就是多倍体Survey的第一个认知门槛：分布变得极其复杂，经典“单主峰+半峰”的二倍体判读经验基本失效。套用二倍体模型，软件会自动把最大的那个峰当成单拷贝主峰，要是最大峰实际是4λ位点贡献的，算出来的基因组大小就会整体偏小或偏离真实单倍体大小。更麻烦的是，低频区（λ峰附近）的大量真实Kmer会被误当成测序错误或者高杂合信号，导致杂合度估计全面失真。

2. 做多倍体Survey之前，必须先确认的三件事

2.1 倍性不是猜的，也不是抄的

这句话听起来像废话，但我在实际项目里见过太多次因为倍性信息没核实就开跑，最后整个Survey白做的案例。很多物种在文献里写的倍性资料其实是几十年前的细胞学观察，有的物种不同地理种群之间存在倍性变异，还有的所谓“四倍体”其实是异源四倍体，只是形态上接近同源四倍体。这些都会导致你的Kmer分析策略完全不同。

动手之前，建议至少确认三件事。第一，物种的染色体基数和倍性，最好有流式细胞术或者染色体压片的数据支撑；第二，要弄清楚是同源还是异源多倍体，这直接决定后面用哪个模型去拟合；第三，如果物种存在多倍体复合群，要确认你手头的材料跟文献记载的倍性一致。如果条件允许，最靠谱的做法是从同一份DNA样品里取一小份做个低深度的Kmer预扫描，用smudgeplot这类工具先快速判断一下倍性和杂合模式，再决定正式分析参数。这一步的成本极低，但能省掉后续大量返工。

2.2 K值选择：多倍体里第一个坑

二倍体Survey里，K值选择相对宽容，17到25之间一般影响不大。但多倍体对K值非常敏感。

原因在于Kmer是“序列身份”的最小单位。K越大，单个Kmer包含的信息量越大，对序列差异越敏感——只要同源拷贝之间有一个SNP差异，跨越这个位点的Kmer就会无法比对到一起，被拆分成两条不同的Kmer。在二倍体中，这种拆分程度跟杂合度有关；在同源多倍体中，四条同源染色体之间每多一个SNP，就可能把Kmer拆成两到四份，低频峰大规模膨胀。如果把K选得太大（比如27、31），多倍体物种的低频Kmer会占到总量的30%以上，主峰被严重压缩，基因组大小估计自然就跑偏。

反过来，K也不能太小。K=15、17这类低K值虽然对序列差异不敏感，但随机出现在基因组上的概率增大，而且测序错误产生的错误Kmer比例大幅上升。多倍体基因组往往比较大，需要较高的内存来做Kmer计数，如果大量内存花在处理错误Kmer上，也是种浪费。

我自己的经验是，多倍体物种最稳妥的K值是21或23，具体可以两个都跑一遍对比。K=21对于大多150bp的二代测序数据来说，既能保持对同源拷贝差异的容忍度，又不会产生太多随机匹配和错误Kmer。如果物种基因组高度重复或者套数特别高（六倍体以上），可以试试K=19；如果倍性不高且同源拷贝之间SNP密度较低，K=25也可以接受。关键是要对比不同K值下估算的基因组大小是否稳定，稳定才是可信的。

2.3 数据量与质控：多倍体经不起脏数据

二倍体Survey要求的数据量通常30-50X就够了，但多倍体最好是50X以上，能到100X更稳。原因是多倍体的Kmer分布峰会因为拷贝数差异而拉宽，如果测序深度不够，λ、2λ这些小深度峰会被淹没在背景噪声里，根本分辨不出来。我做过一个同源四倍体的项目，25X数据量时峰会糊成一团，补测到约60X后才勉强能拆出几个峰。多倍体的Survey真的不建议省数据。

质控方面，有几个点要特别留意。接头残留和低质量碱基会产生大量低频错误Kmer，多倍体本身的低频信号就多，两者叠加之后很难区分，所以正式分析前一定要做严格的质量修剪。最容易被忽略的是PCR重复：在做Survey时，如果建库过程中PCR循环数较多，会引入深度偏倚，导致Kmer主峰展宽。有条件的话用PCR-free文库，没有的话可以在Kmer计数前用工具去重，不过要小心过度去重把真实的低覆盖区域也去掉。另外，叶绿体、线粒体或者共生微生物污染也会在Kmer分布里形成独立的高深度峰，虽然一般不影响主峰判读，但如果污染比例很高，就会干扰模型拟合。建议先用BlobTools之类工具看看组分情况，或者至少比对一下线粒体叶绿体序列评估污染比例。

3. 多倍体Kmer Survey的完整实操流程

3.1 Kmer计数：Jellyfish与KMC的取舍

Kmer计数的工具现在主流就是Jellyfish和KMC两个，多说一句，它们都足够成熟，选哪个更多看数据规模和计算资源。

Jellyfish的特点是内存友好、IO开销适中，适合单机分析。它用一个自适应的哈希表来存Kmer计数，命令简单直接，判定自由度比较高。对于几个Gb的中等基因组，Jellyfish完全够用。它的一个潜在问题是，如果设置的内存上限（-s参数）过小，会出现Kmer被扔掉的情况，对后续直方图影响很大。说白了就是桶不够大，元素放不下溢出丢弃了。

KMC在构建Kmer计数时利用了磁盘临时存储加排序的思路，能用相对小的内存处理超大基因组，适合动辄几十Gb的高杂合多倍体基因组。它把输入序列先切好，利用临时文件分批排序归并，最终得到的直方图统计同样稳定。对于复杂的多倍体项目，我更倾向用KMC，因为大基因组下不用担心哈希表满了丢Kmer。

两条典型命令如下。Jellyfish：

jellyfish count -m 21 -s 10G -t 24 -C -o reads.jf <(zcat reads_1.fq.gz reads_2.fq.gz) jellyfish histo -h 1000000 reads.jf > reads.histo

KMC：

# 先准备一个文件列表，每行一个fastq路径（注意是传给KMC的中间参数格式） echo -e "reads_1.fq.gz\nreads_2.fq.gz" > fq_list.txt kmc -k21 -m10 -t24 -ci1 @fq_list.txt kmc_out workdir kmc_tools transform kmc_out histogram reads_k21.histo -cx1000000

这里有两个参数要专门解释一下。-C对于Jellyfish是把正负链的Kmer合并计数，双端测序默认建议加上，不然后续直方图深度会变成单链值，所有估算结果都得乘2，容易出错。KMC的-h后面是直方图最大深度，多倍体因为存在高拷贝峰，建议设得大一些，比如至少1,000,000，不然高深度尾巴被截断，看不到完整的分布形态。

3.2 直方图判读：怎么看懂多倍体的峰型

数据跑完之后，第一步不是急着丢进GenomeScope，而是先把直方图画出来肉眼看一下。我自己习惯用R画一个简单的频数分布图，横轴深度截到200左右就够了，纵轴用对数坐标更清楚。

看多倍体直方图，要看几个核心特征。

第一，看有没有一个明显的“主峰”以及它的位置。这个峰对应的深度是后面所有计算的地基。如果是二倍体，主峰代表单拷贝深度λ；如果是同源多倍体，直方图往往在λ、2λ、3λ、4λ位置依次出现峰或肩部，需要判断哪个是真正的单拷贝λ峰。判断依据是：λ峰通常是深度最低的那个峰，且它的位置约等于平均测序深度的一半以下。假如你的总测序深度是60X，二倍体的主峰应该在60附近；同源四倍体则可能在30、60、90、120都会看到信号，其中最低的30才是单拷贝深度。

第二，看低频区域。深度1-5之间如果有一个非常陡峭的下降，这是测序错误Kmer的特征信号，正常处理时可接受；如果低频区域在深度10-20还是明显鼓起来，说明要么杂合度很高，要么同源拷贝的分化导致大量低频Kmer，这是多倍体的典型表现。

第三，看主峰右侧是否有长长的尾巴。高深度尾巴说明存在高拷贝重复序列家族。多倍体物种经常伴随大量转座子扩张，这个尾巴会比二倍体更明显。

画完图，心里大致有谱之后，再进入模型拟合阶段。这里强烈建议把直方图和倍性假设结合起来看，不要盲目相信软件输出的数字。如果直方图形态跟软件报出的倍性参数明显矛盾，比如软件说这是二倍体高杂合，但直方图在4个位置都有等间距峰，那就要警惕了。

3.3 基因组大小与杂合度估算：GenomeScope与findGSE

当前用得最多的模型拟合工具是GenomeScope和findGSE，两个思路不太一样。GenomeScope默认用二倍体模型，通过负二项分布叠加拟合杂合Kmer和纯合Kmer的混合分布，输出的参数包括基因组大小、杂合度、重复序列比例等。它后来也加了polyploid模式的选项（在网页版和R脚本里可以通过参数指定倍性），但实际使用中，多倍体模式对同源四倍体的拟合仍然不稳定。

findGSE则是专门为处理复杂基因组设计的方法，核心思路是把Kmer频率分布看作多个泊松分布的混合，每个泊松分量对应Kmer在基因组中不同拷贝数的状态。它用期望最大化算法去估计每个拷贝数分量的参数，最终估算基因组大小。相比于GenomeScope，findGSE不需要预先指定精确的倍性，而是让数据说话，所以对同源多倍体更友好一些。

实际操作时，我建议两个工具都跑一遍，并且用以下几种配置做交叉验证：

GenomeScope命令（R脚本版）：

Rscript genome_scope.R reads_k21.histo 21 150 output_genomescope

如果二倍体模型拟合失败，可以先尝试加上最小深度截断：

Rscript genome_scope.R reads_k21.histo 21 150 output_genomescope -l 5

这个-l 5的意思是忽略低于5x的Kmer，可以把低频错误和部分多倍体低频信号过滤掉，让模型更容易收敛。不过要小心，多倍体的单拷贝峰可能就在深度几到十几的位置，如果阈值设得太大，会把真实的低拷贝信号也滤掉，导致基因组大小低估。

findGSE命令：

Rscript findGSE.R reads_k21.histo 21

findGSE的输出包含不同拷贝数状态的分布，可以看它估计的“1-copy kmers”占比和基因组大小。如果它报告的多拷贝Kmer比例很高，那基本坐实了这是一个多倍体或者高重复基因组。

3.4 用smudgeplot验证倍性假设

如果说GenomeScope和findGSE是从“总量”角度估算，那smudgeplot就是从“关系”角度验证倍性。这个工具的思路很巧妙：它利用双端reads的插入片段信息，找出那些在reads两端成对出现的、序列高度相似的Kmer对，然后根据这些Kmer对的深度比值画一张二维图。每个点的横纵坐标是两个Kmer分别的覆盖度，点的聚集区域对应不同的倍性模式。

具体来说，同源四倍体里，如果一个Kmer存在于四份拷贝（AAAA型），另一个Kmer只存在于一份拷贝（A型），那么两个Kmer的深度比会很悬殊；如果两个Kmer各自存在于两份拷贝（AA和AA），深度比接近1:1。通过观察图中点的分布模式，可以很直观地判断样本是二倍体、同源四倍体，还是异源四倍体。

smudgeplot的使用流程大致是：先用Jellyfish或KMC产出Kmer计数，再从中提取支持配对关系的Kmer，最后绘图。这个工具的优势在于，它不需要知道真实的基因组大小，就可以给出倍性和杂合模式判断。对于多倍体项目中“倍性说不清”的情况，smudgeplot基本是救命稻草级别的验证手段。我在很多项目里都是靠它的图说服合作团队修正倍性认知的。

4. 避坑清单与问题排查实录

4.1 低频峰到底是杂合还是测序错误

多倍体Survey最常见的困惑就是低频区的一堆峰，分不清哪些是真实信号，哪些是测序错误。我的经验是三步走。

第一步，看深度。测序错误Kmer的深度极低，通常集中在1-3x，而且随深度增加急剧衰减。如果你看到深度10-30之间有一片鼓包，这不太可能是纯错误信号，更可能是杂合或同源多拷贝信号。

第二步，看K值。跑两个不同K值（比如21和25），如果某个低频峰在两个K值下都稳定存在，那它大概率是生物信号；如果只在低K值下出现，可能是随机匹配和错误Kmer的产物。

第三步，看材料。如果一个被文献描述为高纯合的材料在低频区有巨量信号，那就要怀疑是不是样品搞混了、存在严重近交衰退或者实际倍性跟认知不符。反过来，如果一个已知高杂合的物种低频信号反而很弱，那也可能是样品来源单一导致杂合度意外偏低。生物数据一定要结合生物学背景来解读，不能只盯着统计输出。

4.2 基因组大小算出来离谱，先查这四件事

如果估算的基因组大小跟流式细胞术结果或者预期值差得离谱，别慌，按顺序排查以下四件事。

第一，主峰识别是不是错了。多倍体的主峰如果选错（比如把2λ峰当成λ峰），大小会直接翻倍或减半。解决方法是回到直方图，用总测序深度反推λ值。比如你的数据平均深度是50X，不考虑边界修正，λ应该在50附近；如果在50附近的峰并不是最高峰，而是后面还有一个更大的峰，那说明最高峰可能是重复或多拷贝信号。

第二，Kmer计数有没有丢数。如果Jellyfish内存设置太小，大量稀有Kmer被丢弃，Kmer种类数被低估，算出来的基因组偏小。KMC模式下一般不会出现这种问题，但要注意临时文件磁盘空间是否够用。

第三，测序数据有没有混入非目标DNA。细菌、真菌、叶绿体污染的白菜价，都会贡献额外Kmer，导致基因组大小高估。可以看高频尾巴是不是有独立于主峰的尖锐峰，这种往往是细胞器DNA的标志。

第四，K值影响是否一致。如果K=19、21、23三个K值估算出的基因组大小差异超过20%，那说明数据里低频信号太强，模型拟合不稳定，需要尝试findGSE或者调整GenomeScope的过滤参数，而不是纠结于某个K值的单次结果。

4.3 多工具结果打架怎么办

GenomeScope给出300Mb，findGSE给出600Mb，流式说500Mb——这种情况在多倍体项目里太常见了。不要试图找一个“正确”的工具，而是要理解每个工具在模型的哪里做了简化。

GenomeScope对同源多倍体的二倍体拟合，倾向于把多拷贝Kmer压进“重复序列”参数里，因此它报出的基因组大小接近单倍体基因组含量，但对总基因组大小（含多套染色体）会低估。findGSE如果把低拷贝Kmer都当成独立拷贝，报出的基因组大小会偏大，更接近流式测到的总DNA含量。流式细胞术给的是核DNA总量，不区分多拷贝同源的序列冗余。三者其实在测量不同层面的“基因组大小”，互相印证才能得到完整判断。

具体来说，如果同一个材料同时拿到流式、GenomeScope、findGSE三个结果，我建议这样解读：流式给的是上限参考；findGSE如果接近流式，说明多拷贝Kmer的处理比较合理，同源序列在测序层面的确大多独立存在；GenomeScope如果显著小于flow值，则说明有一个不小的比例被归入了重复或多拷贝分量中。对下游组装来说，组装团队更关心的是“总共有多少碱基要组装”，所以一般以接近流式或稍低的值作为预算；对注释和分析来说，则更关注单拷贝核心基因组的量。

4.4 与流式结果对不上的排查思路

如果Kmer结果和流式差异实在太远，比如相差一倍，大概率是倍性判断或者样品一致性问题。

最典型的情况：材料实际是同源四倍体，但四套同源染色体之间的序列相似度极高（在Kmer层面几乎完全一致）。这种情况下，流式测到的是四倍DNA量，Kmer计数却把多套同源拷贝合并当作重复处理，算出来接近单套基因组大小。做Survey的人如果只拿到流式结果，看到Kmer估算差一半，很容易误判为测序样品搞错。排查方法是看直方图是否存在等间距多峰，同时跑smudgeplot看倍性图谱。

还有一种情况是流式取样和测序取样不是同一个个体或同一个组织。多倍体物种常有混合倍性现象，不同部位的倍性都可能不同。如果流式材料和测序材料来源不一致，差异可能是真实的生物学差异，而不是技术问题。这种时候只能重新确认材料，没有别的捷径。

下面是一个我整理的速查表，供大家在实际项目中快速定位问题。

4.5 实操中的几个补充心得

最后分享几个SOP里不会写、但实际项目里很有用的经验。

第一，多倍体Survey不要只跑一个K值，至少跑19、21、23三个，把三个K值下估算的基因组大小列个表对比。稳定在10%差异范围内的结果才值得信任，偏差大就要回头找原因。这个习惯花不了多少计算时间，却能让你免于被单个K值的偶然误差带偏。

第二，直方图的纵轴一定用对数坐标。多倍体数据的动态范围极大，主峰可能高出低频区几个数量级，线性坐标下低频区域直接糊成一条线，什么细节都看不见。

第三，测序深度判断要结合Kmer深度和碱基深度两个概念。Kmer深度大约等于碱基深度乘以(L-K+1)/L，对150bp读长和K=21，这个系数是0.867。有的工具报告直接用了Kmer深度，有的用碱基深度，换算关系搞混会带来8%左右的偏差，在几百Mb的尺度上就是几十Mb的误差，足以干扰判断。

第四，如果项目最终目的是做染色体级别组装，Survey阶段一定要保留好Kmer计数文件和直方图，后续组装后的完整性评估（比如BUSCO对比）还需要用到这些信息，别做完就删。

我自己在一次同源六倍体项目中吃了大亏，当时用默认二倍体流程跑出来基因组大小是流式的三倍，杂合度报出8%，整个人懵了一整天。后来核对发现是六倍体的2x峰被当成了主峰，低频峰全被算成了杂合信号，整个Survey结论完全错误。从那次以后，我对多倍体材料一律先画直方图肉眼盯着看半小时，再决定用哪个模型拟合。仪器和软件再智能，最终对生物学理解的把关还是在人自己身上。这个习惯，建议所有做复杂基因组Survey的朋友都培养起来。