单倍体组装工具Hifiasm简介及基本运行命令(一)
随着测序技术的发展及新的组装工具的不断开发应用，基因组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
本文使用 文章同步助手 同步
