随着测序技术的发展及新的组装工具的不断开发应用,基因组denovo测序及组装进入了Genomic2.0时代,我认为Genomic2.0时代的标志有两点:1. 三代长读长测序及Hi-C测序技术在基因组denovo测序上的用;2.组装方法上,Canu和Hifisam等工具不断被开发应用出来,有的工具极大的降低了算力要求,有的工具能够将基因组组装到单体型水平,也就是将同源或非同源的两套多套染色体分别组装出来,因此,最近几年,不仅很多物种的基因组被公布,而早些年间即使被公布了的基因组,也都利用新的测序及组装策略进行了更新。今天我先学习Hifiasm工具。
一.Hifiasm工具简介:
Hifiasm是哈佛大学李恒团队提出的一种全新的单倍体基因组组装算法, 2021年2月份发表在Nature Methods上[ref1]。它可以多线程运行,对计算资源消耗教少,组装快,结果准确性和连续性较高。Hifiasm (Hi-C) 针对PacBio HiFi (High-Fidelity) 长读长测序技术和Hi-C (High-Throughput Chromatin Confirmation Capture) 测序技术进行了全新的设计。该算法结合了HiFi数据中精确的局部单倍型信息和Hi-C数据中的长距离互作用信息以达到全局定相 (phasing),从而获得不依赖亲本信息的染色体级别的单倍型组装结果。为了进一步提高组装质量,作者充分利用了组装图中的结构信息,以及其前期研究中的Graph-binning等策略。
二.算法简介
Hifiasm组装主要分为三步。
Step1: 测序错误碱基纠错
尽管Hifi reads准确性已经很高了,但仍然会有部分测序(<1%)错误,Hifiasm会先通过所有序列的相互比对(all vs all),对测序错误进行纠正。在比对中,基于reads间的overlap关系,如果同一个位置的reads出现两种碱基类型,且每个碱基类型至少有3条reads支持,那么这个位置会被当作杂合位点,即一个SNP被保留,否则,视作测序错误,将被纠正(默认三轮纠错)。值得注意的是,Hifiasm只使用相同单倍型的数据进行纠错,从而避免过度校正,保留来自不同单倍型的杂合变异信息。在这一步,Hifiasm可以对杂合SNP进行定相(phasing)。
Step2: 构建分型字符串图(phased string graph)
根据序列之间的重叠关系,构建分型字符串图string-graph。Hifiasm以reads作为顶点,一致的overlap重叠区域作为边,保留全部的气泡(bubble)即保留了所有的杂合位点(图1),因而可以保留下来基因组上全部的单倍型信息,以便后续对于单倍型的处理。
图1. Hifiasm组装算法示意图
Step3: 单倍体分型组装
如果没有额外的信息,Hifiasm在输出序列时会任意选择气泡的一侧构建初级组装,删除多余的单倍体,输出结果类似Falcon unzip和HiCanu的主要组装结果(primary contigs)。优于HiCanu需要依赖第三方工具去除dups序列,Hifiasm内部实现了去除dups的算法优化,简化了流程。如果有来自父母本的测序数据,Hifiasam可以利用亲本特有的Kmer在图上识别出了父母本的序列,从而得到来自父母本的单倍体基因组序列。
在基于父母本特有Kmer时,区别于TrioCanu软件的trio-binning策略,先将三代reads区分为来自父本、母本以及部分无法区分的reads后对区分后的reads分别组装获得了子代的两套单倍体序列,Hifiasm使用了graph-binning的策略对此进行了改进。它不预先划分reads,而是在string-graph中对reads进行标记。因此在一个较长的bubble中,即使只有一小部分reads被正确标记,hifiasm也可以正确地将其定相。通过这种方式,可以避免因为reads划分错误而引入的错误位点和组装断裂,从而获得更完整和更准确的单倍体组装结果[ref2]。
三.软件使用
1.软件及测试数据下载
Github链接:https://github.com/chhylp123/hifiasm;
下载后make编译;
下载测试数据:
wget https://github.com/chhylp123/hifiasm/releases/download/v0.7/chr11-2M.fa.gz
2.运行程序;
hifiasm使用时根据已有的数据分为三种模式: 2.1.只有HiFi数据(基本)模式; 2.2.有Hi-C数据的Hi-C模式;2.3.有双亲二代测序的Trio-binning模式。
2.1# Run on test data,基本模式,
./hifiasm -o test -t4 -f0 chr11-2M.fa.gz 2> test.log
awk '/^S/{print ">"$2;print $3}' test.bp.p_ctg.gfa > test.p_ctg.fa # get primary contigs in FASTA
参数解释:-o 输出文件前缀, -f0 小数据使用,-t 线程数
awk提取主要的contig,这句话意思是对S开头行处理,提取序列名称$2和序列$3,获得超长的contig序列;
可选参数--primary: 不组装分型,只有primary和alternate的组装结果
运行完成后需要关注的结果 (prefix表示前缀):
test.bp.hap1.p_ctg.gfa: haplotype1的部分分型的contig graph;
test.bp.hap2.p_ctg.gfa: haplotype2的部分分型的contig graph;
test.bp.p_ctg.gfa (Primary assembly contig graph):主要contig的assembly graph, 对于低杂合度物种来说,优先选择该文件;对于高杂合度物种,该结果代表其中一个单倍型;
test.bp.p_utg.gfa(Haplotype-resolved processed unitig graph without small bubbles): 无小气泡的单倍型解析, 在raw unitig graph基础上过滤小的bubble,去掉由于体细胞突变和数据背景噪音引起的small bubbles(这个并不是真正的单体型信息),对于高度杂合基因组物种优先选择这个结果;
test.bp.r_utg.gfa(haplotype-resolved raw unitig graph in GFA format): 保留了所有的单倍型信息,包括体细胞突变和重复测序错误;
*.bin文件:运行时的纠错和相互比对的结果;
其它结果:有的网友还提到了一个结果,我这次没有生成:
prefix.a_ctg.gfa(Alternate assembly contig graph):组装出来的另一套单体型基因组结果。
对于2.2.有Hi-C数据的Hi-C模式;2.3.有双亲二代测序的Trio-binning模式,过段时间我再跑。
四.日志信息及参数调整
通常使用默认参数就可以,要根据日志信息判断是否需要进行参数调整,最主要的日志信息是Kmer图,从而判断hifiasm是否能够正确的找到纯合峰,杂合峰的所在位置。如果hifiasm没有找对纯合峰所在的位置,会导致基因组大小不符合预期,
对于杂合率高的样本,一个常见的问题是分型的结果两套基因组差别较大,需要为-s设置更小的值(默认值:0.55)。
还有其它参数引用ref3,xuzhougeng的分享:
如果序列不够长,片段化明显,则可以尝试增加 -D 和 -N, 虽然会增加运行时间,但是会提高重复区域的分辨率。如果后续的Hi-C,或者BioNano发现hifiasm组装结果有比较多错误组装,则可以适当降低 --purge-max, -s和 -O。或者设置 -u 关闭post-join 步骤,hifiasm通过该步骤提高组装的连续性。
五.参考:
Ref1:Haplotype-resolved de novo assembly using phased assembly graphs with hifiasm
Ref2: https://zhuanlan.zhihu.com/p/283131167
Ref3:https://www.jianshu.com/p/6d79690dce5d?ivk_sa=1025883j
Ref4: https://hifiasm.readthedocs.io/en/latest/trio-assembly.html
本文使用 文章同步助手 同步