干过几年基因组组装的人大概都体会过那种感觉拿到一个自然界里普普通通的多倍体物种按着二倍体的Survey流程跑一遍Kmer出来的数字怎么看怎么不对劲——基因组大小要么翻倍要么直接砍半杂合度高得离谱有的甚至直接给你报一个“Optimization failed”就罢工了。这时候如果经验不足第一反应通常是怀疑测序数据有问题转头去重新测序、重新建库白白烧掉不少钱和时间。但实际上问题大概率出在Kmer分析模型与多倍体基因组结构不匹配上。这正是复杂基因组Survey分析中最值得讲清楚、也最容易被忽略的事情。这篇文章不打算从Kmer原理的ABC开始讲而是直接围绕“多倍体Survey”这个组合展开。我会把二倍体跟多倍体在Kmer分布上的本质差异、动手前要确认的信息、完整实操流程以及我这些年踩过的坑一一说清楚。无论你是刚开始接触基因组Survey的研究生还是已经跑过几次流程但始终被多倍体折磨的从业者相信都能从里面找到可用的东西。1. 先把底层逻辑说透Survey到底在算什么1.1 Kmer的“小卡片记账法”搞懂多倍体为什么难先得清楚Kmer Survey在干什么。测序得到的reads是一段一段的碱基序列Kmer就是把每条reads切成固定长度为K的短片段。比如一条150bp的reads以K21为例会得到150-211130个21bp的小片段。这些Kmer记完数之后我们画一张直方图横轴是某个Kmer在数据里出现的次数深度纵轴是这个深度对应的Kmer种类数量。这张直方图之所以能用来估算基因组大小背后的逻辑很简单。假设一个基因组大小为G均匀测序覆盖度为C那么基因组上每个位置平均会被C条reads覆盖对应到Kmer深度也大约是C准确说是C乘以(L-K1)/L这个边界修正系数。整个基因组有多少种Kmer理论上大约是G种每个碱基位置产生一个Kmer边界效应忽略不计。总Kmer条数则是G×C。于是基因组大小就可以用“总Kmer条数除以主峰深度”估算出来。举个例子一个500Mb的基因组50X测序150bp读长K21那么总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会被四份拷贝同时贡献readsKmer深度直接变成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也不能太小。K15、17这类低K值虽然对序列差异不敏感但随机出现在基因组上的概率增大而且测序错误产生的错误Kmer比例大幅上升。多倍体基因组往往比较大需要较高的内存来做Kmer计数如果大量内存花在处理错误Kmer上也是种浪费。我自己的经验是多倍体物种最稳妥的K值是21或23具体可以两个都跑一遍对比。K21对于大多150bp的二代测序数据来说既能保持对同源拷贝差异的容忍度又不会产生太多随机匹配和错误Kmer。如果物种基因组高度重复或者套数特别高六倍体以上可以试试K19如果倍性不高且同源拷贝之间SNP密度较低K25也可以接受。关键是要对比不同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。两条典型命令如下。Jellyfishjellyfish 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.histoKMC# 先准备一个文件列表每行一个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在基因组中不同拷贝数的状态。它用期望最大化算法去估计每个拷贝数分量的参数最终估算基因组大小。相比于GenomeScopefindGSE不需要预先指定精确的倍性而是让数据说话所以对同源多倍体更友好一些。实际操作时我建议两个工具都跑一遍并且用以下几种配置做交叉验证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 21findGSE的输出包含不同拷贝数状态的分布可以看它估计的“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值影响是否一致。如果K19、21、23三个K值估算出的基因组大小差异超过20%那说明数据里低频信号太强模型拟合不稳定需要尝试findGSE或者调整GenomeScope的过滤参数而不是纠结于某个K值的单次结果。4.3 多工具结果打架怎么办GenomeScope给出300MbfindGSE给出600Mb流式说500Mb——这种情况在多倍体项目里太常见了。不要试图找一个“正确”的工具而是要理解每个工具在模型的哪里做了简化。GenomeScope对同源多倍体的二倍体拟合倾向于把多拷贝Kmer压进“重复序列”参数里因此它报出的基因组大小接近单倍体基因组含量但对总基因组大小含多套染色体会低估。findGSE如果把低拷贝Kmer都当成独立拷贝报出的基因组大小会偏大更接近流式测到的总DNA含量。流式细胞术给的是核DNA总量不区分多拷贝同源的序列冗余。三者其实在测量不同层面的“基因组大小”互相印证才能得到完整判断。具体来说如果同一个材料同时拿到流式、GenomeScope、findGSE三个结果我建议这样解读流式给的是上限参考findGSE如果接近流式说明多拷贝Kmer的处理比较合理同源序列在测序层面的确大多独立存在GenomeScope如果显著小于flow值则说明有一个不小的比例被归入了重复或多拷贝分量中。对下游组装来说组装团队更关心的是“总共有多少碱基要组装”所以一般以接近流式或稍低的值作为预算对注释和分析来说则更关注单拷贝核心基因组的量。4.4 与流式结果对不上的排查思路如果Kmer结果和流式差异实在太远比如相差一倍大概率是倍性判断或者样品一致性问题。最典型的情况材料实际是同源四倍体但四套同源染色体之间的序列相似度极高在Kmer层面几乎完全一致。这种情况下流式测到的是四倍DNA量Kmer计数却把多套同源拷贝合并当作重复处理算出来接近单套基因组大小。做Survey的人如果只拿到流式结果看到Kmer估算差一半很容易误判为测序样品搞错。排查方法是看直方图是否存在等间距多峰同时跑smudgeplot看倍性图谱。还有一种情况是流式取样和测序取样不是同一个个体或同一个组织。多倍体物种常有混合倍性现象不同部位的倍性都可能不同。如果流式材料和测序材料来源不一致差异可能是真实的生物学差异而不是技术问题。这种时候只能重新确认材料没有别的捷径。下面是一个我整理的速查表供大家在实际项目中快速定位问题。现象可能原因排查手段估算大小约为预期的2倍主峰误判为半深度峰异源多倍体被当单倍体算检查峰位用平均深度反推λ跑smudgeplot确认倍性估算大小约为预期的一半同源多倍体序列高度一致Kmer重复分量过大看直方图低频峰形态与流式结果交叉验证杂合度报告特别高5%同源多倍体低频Kmer被当杂合信号样品污染换findGSE模型用smudgeplot确认杂合模式GenomeScope无法收敛数据低频信号太复杂初始参数离真实值太远加-l截断低频换低K值把直方图裁到合理深度范围直方图呈现多个等间距峰同源多倍体的1x/2x/3x/4x峰用findGSE拟合不要强行套二倍体模型直方图有尖锐高深度峰细胞器DNA或污染物种富集检查是否需去除细胞器reads评估污染来源4.5 实操中的几个补充心得最后分享几个SOP里不会写、但实际项目里很有用的经验。第一多倍体Survey不要只跑一个K值至少跑19、21、23三个把三个K值下估算的基因组大小列个表对比。稳定在10%差异范围内的结果才值得信任偏差大就要回头找原因。这个习惯花不了多少计算时间却能让你免于被单个K值的偶然误差带偏。第二直方图的纵轴一定用对数坐标。多倍体数据的动态范围极大主峰可能高出低频区几个数量级线性坐标下低频区域直接糊成一条线什么细节都看不见。第三测序深度判断要结合Kmer深度和碱基深度两个概念。Kmer深度大约等于碱基深度乘以(L-K1)/L对150bp读长和K21这个系数是0.867。有的工具报告直接用了Kmer深度有的用碱基深度换算关系搞混会带来8%左右的偏差在几百Mb的尺度上就是几十Mb的误差足以干扰判断。第四如果项目最终目的是做染色体级别组装Survey阶段一定要保留好Kmer计数文件和直方图后续组装后的完整性评估比如BUSCO对比还需要用到这些信息别做完就删。我自己在一次同源六倍体项目中吃了大亏当时用默认二倍体流程跑出来基因组大小是流式的三倍杂合度报出8%整个人懵了一整天。后来核对发现是六倍体的2x峰被当成了主峰低频峰全被算成了杂合信号整个Survey结论完全错误。从那次以后我对多倍体材料一律先画直方图肉眼盯着看半小时再决定用哪个模型拟合。仪器和软件再智能最终对生物学理解的把关还是在人自己身上。这个习惯建议所有做复杂基因组Survey的朋友都培养起来。