做基因组组装这几年我最大的体会是很多人不是不会跑工具而是不知道什么样的数据组合能拿到什么样的结果跑完也不知道怎么判断组装得好不好。有次拿到一套Q20级别的PacBio HiFi数据N50在17kb左右我照着常规流程直接丢给组装软件出来的contig N50才几兆和预期差了一个数量级。后来把同一条DNA的ONT超长读取加进去重新组装同样的HiFi数据长青序列N50直接到了几十兆还能看到朊粒的轮廓。这就是PacBio HiFi和ONT超长读取组合的价值HiFi贡献准确的基础单元ONT超长读取把单元串成一整条长链两者合起来才能搞定复杂基因组里那些重复区和结构变异区域。这篇保姆级教程我从拿到下机数据开始写覆盖数据质控、覆盖度估算、hifiasm参数调优、结果验证和踩坑经验适合第一次完整组装基因组但还没摸清整套流程的人。1. 先弄明白为什么HiFi和ONT超长读取要组合而不是二选一1.1 两种数据的技术定位完全不同PacBio HiFi读长的核心优势是准确率高单碱基准确率一般能到Q20以上甚至不少跑得好的批次能到Q30。它把同一酶分子上的CCS子读段循环编码校正相当于把同一个位置读了很多遍然后取一致序列所以既有长读长的跨度又有接近短读长的准确度。但HiFi有一个天然限制它的读长取决于测序过程中DNA聚合酶的存活时间实际下来通常就是15到25kb这个区间能上30kb的样本不多。ONT的特点是读长几乎没有上限。只要你把DNA提取做好高分子量DNA能保持完整一次测序读出100kb甚至200kb以上的读长并不稀奇。比如我之前一个植物样本用特殊的低剪切力提取方案之后ONT超长读长的N50做到了80kb最长的读长超过200kb。代价是单读准确率不如HiFi即便用了最新的Q20化学单读水平也就在97%到99%之间波动逐碱基错误里还有不少是清单式缺失。过去有人纠结到底是选HiFi还是ONT实际做过复杂基因组之后会更务实这不是二选一而是互补。HiFi负责提供质量高的“砖块”ONT超长读取负责把砖块按正确的顺序串起来。简单基因组可能单个HiFi就够了但遇到高度重复的着丝粒、rDNA阵列、大型结构变异区域短一点的读长进入重复区就不知道该往哪边延伸而那个一二百kb的ONT读数可以横跨整个重复区告诉它们“两个拷贝的顺序是这样连的”。1.2 从组装结果倒推哪一类数据解决哪一类问题把目标拆开看会更清楚连续性contiguity主要取决于最长读长能否跨越大的重复单元。ONT超长读取往往能覆盖完整的转座子簇、rDNA串所以对N50的提升极其明显。碱基准确性取决于单碱基层面的共识序列质量这主要是HiFi的贡献。HiFi准确度高出现在最终组装里的单核苷酸错误就少。单碱基分辨率的结构变异两种数据结合后可以在重复区域边界上获得准确的手性。单倍型hifiasm这种工具本身就能处理HiFi中的杂合位点信息谁按照序列组装出两套单倍型。如果你还想区分哪个来自父本、哪个来自母本可配合trio数据。有一个很常见的问题“我做了HiFi是不是ONT就完全不必要了”根据我做过的经验对于1Gb左右的基因组长、重复序列占比超过50%的物种只给HiFi数据组装结果会有大量类托盘状的图结构加上5到10倍的ONT超长读取后大部分可能因折叠起来图变得更线性。自己手上的项目是不是需要可以先看参考物种的重复序列占比和基因组大小再决定要不要投入ONT。某些真菌基因组只有几十Mb重复区域少那就完全没必要加超长读取。2. 下机数据先体检原始数据质量决定了能否组装出染色体级别很多人拿到数据就急着直接跑组装结果出来一堆问题再回头查数据质量。我的习惯是组装前先做一轮完整的数据体检花的时间不多但能省下后面几百个小时的排错时间。2.1 四个必做的初步检查第一步读长N50和读长分布。用seqkit stats或NanoPlot分别统计HiFi和ONT数据。HiFi我一般要求至少看到N50在15kb以上如果低于这个值可能是片段化或者测序酶活下降需要重新评估。ONT超长读取则重点看N50和最长读长N50低于20kb的话对组装的连续性帮助有限我的心理门槛是至少30kb以上且需要有相当比例超过50kb的读取。第二步碱基质量分布。HiFi需要看Q值分布是否集中在Q20以上如果大片低于Q20说明循环共有序列的深度不够或者是样本的碱基修饰干扰了酶活性。ONT超长读取直接看平均质量值最好在Q10以上同时要看是否有低质量尾段。第三步读段比对已知近缘物种。虽然你不用优质参考基因组但对比一下能很快发现数据有没有受污染、insert size是否异常。如果比对的覆盖度均匀性很差有可能是异常拷贝数或者样本混入了其他DNA。第四步k-mer分析。这一步容易被跳过但我强烈建议做。对HiFi数据直接跑一个小的GenomeScope2或KMC估计基因组大小、杂合度和重复度。这可以让你在组装前就知道这个基因组大概有多大、重复比例是多少、杂合度高不高甚至能提前发现样本是不是混了多材料。曾经有个项目用户坚持说测的是一个单株的个体样本但k-mer谱图上明显看到两个峰后来一查是采样时混了两棵不同的植株还好组装前发现不然整个组装流程都要白跑。2.2 覆盖度到底算多少才够用别被通量忽悠了我常用的一套定量标准是HiFi数据做到30x到60xONT超长读取做到5x到15x。看起来ONT超长读取得倍数很低为什么会有用关键在于它的角色不是提高共识准确率而是提供“骨架”和“跨越重复”的信息。理论上一个只覆盖一次的超长读取也能连接两个相距很远的区域所以不必把它堆到和HiFi一样的深度。例如一个1Gb的物种基因组30x HiFi大概对应30Gb的HiFi数据加上10x的ONT超长读取也就是10Gb左右。这么点ONT数据量也许听上去不多但在hifiasm里加入--ul后组装连续性提升会非常直观。覆盖度也不是越高越好HiFi堆到80x以上消耗的是真金白银的测序费用而且对组装质量的提升已经非常有限ONT超长要多的话对建库要求就变得很严格洗DNA片段化很严重的情况下配对产生许多短读反而浪费了机器时间。2.3 数据要不要质控与修剪我的处理原则HiFi数据一般不需要按碱基质量逐段修剪ccs流程输出的已经是共识序列保留原始分数即可。但有一个例外如果某条HiFi读段内部有巨大的质量低谷比如某个200bp窗口Q值跌到Q10以下这种读段在组装时有可能会造成后续local misassembly我会用seqkit按质量筛掉超低质量读段而不是去“修剪”。ONT超长读取的质控要看你是用哪个工具做完孔处理。如果用Dorado或Guppy做完碱基识别后直接输出的是带质量的FASTQ一般我会先用Dorado自带的--min-qscore选项或者chopper做一次低质量过滤我常用Q8作为下限。剪接trim这块需要谨慎ONT读段中的接头序列需要去掉但哪些部分是因为测序噪声哪些是真实序列很难剪得干净。我自己的经验是只要不明显的adaptor序列尽量不要主动把读段尾部截断宁可用软屏蔽的方式处理。超长读段的末端哪怕质量低只要没有严重错误配成一个个对齐时也不会让组装崩掉反而能提供延伸的线索。3. 组装跑了哪些工具hifiasm加ONT超长读取是我目前的默认流程3.1 为什么不推荐再用老一代组装工具硬拼HiFi如果你是第一次做这块建议直接上手hifiasm。它支持纯HiFi模式也支持用--ul接收ONT超长读取是目前处理HiFi最主流的工具之一。它的核心思路是从HiFi读段中的杂合位点来区分两条单倍型构建一个具有单倍型感知haplotype-aware的字符串图组装过程中能直接分型并输出两个单倍型版本以及一个primary版本。早期很多人拿HiFi去跑wtdbg2或者Canu的组装流程也不是不行但设置不少而且产量在重复序列和单倍型分离上不明显。hifiasm这个工具自带省心它不会要你不用把数据先纠错成某种中间格式一次跑完就能输出图组装结果也不需要人为搭建这样的三代-四代流程。所以后面我只讲hifiasm的用法。3.2 具体命令与参数怎么看假设你的目录结构是这样的project/ 00_tools/ 01_data/ sample.ccs.fastq.gz sample.ulont.fastq.gz 02_out/先写好输入文件的软链接然后执行如下命令cd 02_out ln -s ../01_data/sample.ccs.fastq.gz ./hifi.fastq.gz ln -s ../01_data/sample.ulont.fastq.gz ./ulont.fastq.gz hifiasm -o sample -t 64 \ --ul ulont.fastq.gz \ hifi.fastq.gz \ 2 hifiasm.log核心参数只有几个这里解释一下为什么这样设置-t 64线程数。这个按你的服务器CPU核数来设一般建议32到96之间。hifiasm对计算资源不大敏感但超长读段的比对阶段和后续图构建阶段多核能明显缩短时间。--ul ulont.fastq.gz这个参数让hifiasm把ONT超长读段作为数据比对时的“超长辅助”信息不需要你把它们拆成和HiFi一样的单位。它是hifiasm在组装过程中用来解开重复区域分支关系的关键。-o sample指定输出前缀运行完成后会在当前目录下生成sample.bp.hap1.p_ctg.gfa、sample.bp.hap2.p_ctg.gfa和sample.bp.p_ctg.gfa等文件。主集合一般是p_ctg你可以直接拿来做后续分析。如果你只想输出一个主组装不用管phased文件直接用sample.bp.p_ctg.gfa即可。现象想同时输出比对信息用于后续分析而加--write-ec或者不开沟的选项会占用较多内存跑之前先预估一下。纯命令行看到“⌛**时间提示**”是不是进度卡住了实际hifiasm的log可能比较朴素它不像某些工具那样一行行刷新漂亮进度条只要日志文件在正常增长CPU占用率不低就不用担心。3.3 超长读取在hifiasm中到底做了什么从图结构理解参数组装过程中hifiasm的输入读取先构建一个string graph。图中每个节点是一段连续的基因组序列节点之间的边表示这些序列在读取中能形成一个延伸路径。重复区域会在图里形成气泡状、菱形或发卡状的结构组装算法需要一个准则来决定“走哪条边”。HiFi读段只有15到25kb遇到两个相同的串联重复单元时读段往往落在同一个重复内部无法告诉组装算法重复单元结束之后应该接哪边。ONT超长读段如果长达80kb以上可能一次横跨几个重复单元和两侧的独特序列hifiasm就可以依据它来把正确的路径“串”起来。这也是为什么ONT超长读取不需要很多覆盖它不需要覆盖每个位置它只需要在重复区域上提供一条正确顺序的“长距离证据”。但有一点要特别说明--ul只是hifiasm的辅助功能不代表它会用ONT去替换HiFi做共识序列。所有最终输出出来的碱基仍是以HiFi为主体ONT的作用是告诉算法如何拆解图结构。所以你在评估数据时不要觉得“ONT质量低所以不能用”它参与的方式跟HiFi完全是两种逻辑。3.4 计算资源估算与运行时间对我常用这个规模进行估算1Gb基因组30x HiFi加10x ONT在64线程机器上占大约90到120GB内存运行时间通常10到20小时具体取决于读段的错误率和重复度。如果你给的线程数很少比如8线程内存使用不会明显减少因为hifiasm的内存与图规模相关而图规模又取决于基因组复杂度和读取覆盖。这里最可怕的是把内存给得不够导致跑到一半OOM白白浪费几十个小时。我个人的铁律是服务器内存至少按目标基因组大小的150倍来预备。比如1Gb基因组至少准备80GB可用内存更保险上到128GB。如果你有6Gb的复杂繁殖基因组那就需要准备512GB到1TB内存才能比较从容地跑。别看着hifiasm软件小巧图结构构建时内存吃得比你预想大得多。4. 组装完成不算完必须过这三道验证关卡组装出序列只是第一步没有验证的组装结果是不能拿去下结论的。我曾经见过有些初学同学看到hifiasm跑完就高高兴兴搞下游注释结果后来一查很多核心单拷贝基因直接缺失一个污染源让整个结果报废。验证过这些关卡才叫真正“组装完”。4.1 第一关连续性指标和读回比对先要把GFA转成FASTA。这一步我常用gfatools gfa2fa也可以直接用awk提取gfatools gfa2fa sample.bp.p_ctg.gfa sample.p_ctg.fa拿到FASTA后看两个指标contig N50对于1Gb的基因组如果N50在10Mb以下说明组装可能被打得很碎需要怀疑数据或者图结构有问题达到几十Mb是正常的染色体级别的甚至可以把N50做到单个染色体的长度。通过已样本读段回帖比对的覆盖比例用minimap2把原始HiFi数据比对回你的组装结果再统计比对率。正常情况比对率应该在98%以上一致性错误率尽量在0.2%以内。minimap2 -ax map-hifi sample.p_ctg.fa hifi.fastq.gz -t 32 | samtools sort - 8 -o hifi.bam samtools flagstat hifi.bam samtools depth hifi.bam | awk {sum$3;n} END{print sum/n}这个操作同时还能帮你检查有没有组装错误。如果某段区域上读段的覆盖度明显比其他区域低同时有很多硬裁剪hard clip那这个地方很可能有局部错误组装需要进一步看。4.2 第二关BUSCO完整性评估BUSCO是我最信任的“基因层面完整性”检查。它看的是这个门类里普遍高度保守的单拷贝基因在你组装里能否找到。busco -i sample.p_ctg.fa \ -l embryophyta_odb10 \ -o busco_plant \ -m genome \ -c 16这里的-l参数要根据物种选动物用metazoa_odb10哺乳动物用mammalia_odb10真菌用fungi_odb10植物用embryophyta_odb10。不要随便套用别的物种的库否则结果没有意义。对于一个新的高质量基因组BUSCO的complete比例应该在95%以上如果只有80%多大概率存在三种可能数据覆盖不足、基因组杂合导致单倍型相分离后某些基因被拆到不同单倍型上、或者组装里出现了长段缺失。还要注意看duplicated那一列。如果duplicated比例偏高可能是单倍型没有合并好hifiasm输出的两个单倍型被拼进了同一个主组装里。最重要的一个联想BUSCO结果是基因组完整性的必要不充分条件。BUSCO完整不代表组装处处正确因为BUSCO基因只占全基因组很小一部分而那些重复区是否组装正确BUSCO根本看不到。所以还得有第三关。4.3 第三关端粒/着丝粒信号和结构特征这一关我习惯叫“结构感检查”。一个接近完整染色体的组装会在一端或两端出现端粒重复信号一般植物用Arabidopsis型端粒重复序列TTTAGGG来找哺乳动物用TTAGGG来找。你可以在组装出的最长几条contig的序列里搜一下这些重复是否出现在末端。如果大多数长contig两端都摸不到端粒说明它们还是断在重复区内部离染色体级别还差得远。如果两条染色体的两个端分别出现在同一条contig的两端而且中间区域BUSCO完整、覆盖度均匀那这条contig很可能已经对应一整条染色体了。这个时候再用ONT超长读取把空隙填一填配合后续的Hi-C辅助搭建染色体基本就是一篇高质量基因组的骨架了。结构特征之外还有一个小细节我每次必看组装里有没有线粒体和叶绿体序列被错误拼入核基因组。用seqkit搜索已知线粒体基因如COX1或植物叶绿体基因rbcL如果发现核基因组contig里混入了这些细胞器序列也没太大关系可以在后续流程里用mitoHiFi或MitoFinder分离出来。但要记录清楚否则下游基因注释时会多出一堆假基因。5. 组装从零开始的路上这五个坑我踩过希望你一次绕开5.1 拼命追求更高HiFi覆盖度结果只是浪费测序预算我以前接过一个项目对方说已经测了60x的HiFi数据但组装连续性还是上不来于是想再加到100x。我先看了他们下机数据质量发现读长N50只有9kb明显不如其他样本而且他们没做ONT超长读取。后来补做了一次长片段文库提取DNA之前多加了一步低熔点琼脂糖包裹处理ONT N50达到70kb只用10x左右的数据就把原有组装提升到染色体级别。加HiFi覆盖度不是万能解绕开读长瓶颈才是。5.2 hifiasm跑到一半内存不足大多数情况不是服务器不行是数据结构选择不对我在64GB内存的服务器上跑过一个1Gb的基因组以为足够结果跑到一半OOM。后来反复对比才发现问题出在我同时开启了--write-paf和保留中间文件再加上没有限制缓冲区。hifiasm在做图构建时会将大量中间结果缓存在内存中如果你并不需要调试信息建议关闭多余输出选项只保留必要的--ul和输出路径。遇到OOM时直接开启交换文件不是一个好主意速度会慢到让你怀疑人生更实际的做法是增加物理内存或分两步跑先用纯HiFi模式出初步图再用--ul模式从同一数据继续这样内存峰值可控。5.3 超长读取数据量很大但里面掺了很多短读和接头导致组装时间里浪费严重有一次ONT下机数据看着量很大但洗孔basecalling后的平均读长只有6kb我一开始还没警觉。后来用NanoPlot看长度分布才发现50kb以上的真正超长读段占比不到1%。这种数据加入--ul带给hifiasm的帮助非常有限。排查下来罪魁祸首是DNA提取过程中涡旋过度大分子DNA被机械切断。所以建库前一定要做脉冲场电泳或者至少跑一个质检胶看高分子量DNA条带不要只看核酸浓度。一旦确认DNA已降解再怎么调整组装参数都补救不了。5.4 杂合度过高导致主组装里混入两个单倍型BUSCO duplicated异常你可以想象两个不同个体间的SNP很多hifiasm会把它分别建图到两个单倍型上。在输出p_ctg时如果算法没有完全处理好可能会出现某些区域来自单倍型1、另一些区域来自单倍型2的情况体现在BUSCO上就是duplicated比例升高但实际上不是真正的串联重复。这种情况我一般先看同一物种近缘参考的杂合度水平如果杂合度高于1%到2%在做组装前就明确把hifiasm的phasing输出留给后续分析。做下游基因注释时许多人会直接在主组装上注释这会导致等位基因被当成两个不同基因。最省事的做法是用hap1和hap2分别跑BUSCO看哪个更适合作为代表再把另一个作为辅助资源。5.5 不要盲目相信N50N50很高也有可能包含大量错误连接N50只是一个单值指标如果组装有错位连接N50可能还是很漂亮。最典型的情况是旁系同源区域之间发生了错误连接比如两个不同基因家族共享高相似序列ONT读段又没覆盖到足够的独特序列来区分组装算法就误把来自不同染色体的片段拼在一起。这种错误在检查时读回覆盖率和BUSCO都可能是正常的但你会发现某些contig远端序列在其他地方的比对中出现非常奇特的断点。这个时候我最常用的一招是把长时间测序读段重新比对回组装结果专门看那些能跨过contig边界的“连线”。如果一条ONT读段在contig断点处左右各覆盖了很长一段那说明这条边可能还有别的组装证据反过来如果几十条覆盖该断点的读段都只能锚定一侧另一侧完全没有映射那这个连接基本是错的。这里要耐心尤其是做染色体级参考基因组早期多投入一些手动审查时间后面文章投稿时就能少一堆硬伤。今天这篇整体算是一套很朴素的流程。好数据就是好数据坏数据再怎么调参也救不回来反过来数据质量扎实哪怕只用hifiasm一条命令也能得到相当可靠的组装结果。我每次拿到新项目的原始数据都会先花半天把质控k-mer和读长分布看清楚了再启动组装这个习惯帮我挡掉了太多后续的返工。你如果真的决定从零开始组装一个基因组我建议不必一上来追求完美染色体级别先把一条完整可用的pipeline跑通再用ONT超长读取和Hi-C逐步升级这样既稳又能真正理解每一层数据在组装里扮演的角色。