基因组学

hifiasm对HiFi PacBio进行组装

hifiasm是一个能有效利用PacBio HiFi测序技术，在分型组装图(pahsed assembly gprah)中可靠的表示单倍体信息的算法。

流程介绍

hifiasm的分析流程如下，主要分为3个阶段

第一阶段：通过所有序列的相互比对，对前在测序错误进行纠正。如果一个位置只存在两种碱基类型，且每个碱基类型至少有3条read支持，那么这个位置会被当作杂合位点，否则，视作测序错误，将被纠正。

第二阶段：根据序列之间的重叠关系，构建分型的字符串图(phased string graph)。其中调整朝向的序列作为顶点(vertex)，一致重叠作为边(edge)。字符串图中的气泡(bubble)则是杂合位点。

第三阶段：如果没有额外的信息，hifiasm会随机选择气泡的一边构建primary assembly，另一边则是alternate assembly. 该策略和HiCanu，Falcon-Unzip一样。对于杂合基因组而言，由于存在一个以上的纯合haplotype，因此primary assembly可能还会包含haplotigs。HiCanu依赖于第三方的purge_dups, 而hifiasm内部实现了purge_dups算法的变种，简化了流程。如果有额外的信息，那么hifiasm就可以正确的对haplotype进行分型。

hifiasm workflow

安装

hifiasm仅仅依赖 g++和zlib，以及git

1
2
3
4
5
6
7
# 依赖g++和zlib
git clone https://github.com/chhylp123/hifiasm
cd hifiasm && make

# bioconda
conda install -c bioconda hifiasm

通过源码编译的方式安装，需要将hifiasm移动到你的软件目录下，或者将hifiasm的路径加入到环境变量PATH中。

如果是 trio-binning模式，需要额外安装yak

1
2
3
4
5
6
#source code
git clone https://github.com/lh3/yak
cd yak && make
# bioncda
conda install -c bioconda yak

案例展示

hifiasm的使用非常简洁明了，根据已有的数据分为，仅HiFi数据模式，有双亲二代测序的Trio-binning模式和有Hi-C数据的Hi-C模式。

仅有HiFi数据

最基本的用法，会得到两个部分分型的组装

1
2
3
wget https://github.com/chhylp123/hifiasm/releases/download/v0.7/chr11-2M.fa.gz
hifiasm -o test -t 32 chr11-2M.fa.gz 2> test.log

其中 -o定义输出文件的文件名前缀， -t是线程数

运行结束后生成的一系列文件中，我们只需要关注如下几项 (prefix表示前缀）

  • prefix.bp.r_utg.gfa: haplotype-resolved raw unitig graph，记录所有的单倍型信息

  • prefix.bp.p_utg.gfa: 在raw unitig graph基础上过滤小的bubble，

  • prefix.bp.p_ctg.gfa: 主要contig的assembly graph

  • prefix.bp.hap1.p_ctg.gfa: haplotype1的部分分型的contig graph

  • prefix.bp.hap2.p_ctg.gfa: haplotype2的部分分型的contig graph

如果并不需要部分分型的组装，而只想要primary和alternate的组装结果，可以在之前的命令的基础上，加上 --primary参数。

1
hifiasm --primary -o test -t 32 chr11-2M.fa.gz 2> test.log2

由于hifiasm运行时会将步骤中纠错和相互比对的结果保存成 bin 文件，因此重新这一次运行速度会很快

primay模式下输出的文件和之前的类似，唯一的不同在于没有 bp

  • ``prefix.r_utg.gfa: haplotype-resolved raw unitig graph

  • ``prefix.p_utg.gfa: haplotype-resolved processed unitig graph without small bubbles.

  • ``prefix.p_ctg.gfa: assembly graph of primary contigs.

  • ``prefix.a_ctg.gfa: assembly graph of alternate contigs.

我们关心的，可能就是 主要的contig，通过awk进行提取

1
awk '/^S/{print ">"$2;print $3}' test.p_ctg.gfa > test.p_ctg.fa

Trio-binning模式

如果测了双亲，则可以使用trio-binning方法进行更加可靠的分型。分为两步，先用yak统计k-mers， 然后用hifiasm进行组装

案例代码如下

1
2
3
yak count -k31 -b37 -t16 -o pat.yak paternal.fq.gz
yak count -k31 -b37 -t16 -o mat.yak maternal.fq.gz
hifiasm -o NA12878.asm -t 32 -1 pat.yak -2 mat.yak NA12878.fq.gz

输出的文章和之前类似，主要关心其中文件名带 dip 的输出gfa文件

  • prefix.dip.hap1.p_ctg.gfa: 完成分型的父源单倍体 contig图.

  • prefix.dip.hap2.p_ctg.gfa: 完全分型的母源单倍体contig图.

整合Hi-C数据

由于Hi-C数据能够提供远距信息，因此也能用于单倍体分型。只需要加上两个参数， h1接受Hi-C的read1, h2 接受Hi-C的read2

1
2
hifiasm -o NA12878.asm -t32 --h1 read1.fq.gz --h2 read2.fq.gz HiFi-reads.fq.gz

在该模式下，每个contig要么是来自于父亲，要么是来自于母亲。hifiasm会将同一来源的contig放在同一个组装中。需要注意的是，hifiasm未必能够处理好着丝粒附近的区域，另外hifiasm中Hi-C也不会用于进行scaffold。

输出结果中，我们重点关注其中名字带hic的文件

  • prefix.hic.p_ctg.gfa: 主要contig的组装图

  • prefix.hic.hap1.p_ctg.gfa: 完全分型的haplotype1的contig图

  • prefix.hic.hap2.p_ctg.gfa: 完全分型的haplotype2的contig图

  • prefix.hic.a_ctg.gfa : 如果设置了 --primary参数，还会输出该次要contig的组装图

日志和参数调整

绝大部分的时候，我们只需要使用默认参数即可得到相对比较好的结果。但是当默认参数无法达到自己的目的，那我们就需要检查日志信息，阅读相关参数从而优化结果。

日志信息主要分为三项

  • k-mer图: 纯合样本只有一个peak，杂合样本则会有2个peak。

  • 纯合峰的覆盖度: [M::purge_dups] homozygous read coverage threshold: X , 一般会由hifiasm自动推断。

  • 杂合/纯合碱基数目（Hi-C模式）: 在Hi-C模式下，如果纯合的碱基数超过杂合碱基数，那么hifiasm就不容易找对纯合read的所在峰。

对于日志信息，我们最主要关注的就是k-mer图，从而判断hifiasm是否能够正确的找到纯合峰，杂合峰的所在位置。如果hifiasm没有找对纯合峰所在的位置，就需要我们根据k-mer图手动指定 --hom-cov。

对于一个组装结果，最直接的评估标准就是基因组大小是否符合预期，分型的两套基因组是否相差不大，序列是否足够长，是否存在错误组装的情况。

如果基因组大小不符合预期，一般都是hifiasm找错了纯合峰的位置，我们需要手动指定 --hom-cov；如果分型的两套基因组差别过大，则通过降低 -s 调整。如果序列不够长，片段化明显，则可以尝试增加 -D 和 -N, 虽然会增加运行时间，但是会提高重复区域的分辨率。如果后续的Hi-C，或者BioNano发现hifiasm组装结果有比较多错误组装，则可以适当降低 --purge-max, -s和 -O。或者设置 -u 关闭post-join 步骤，hifiasm通过该步骤提高组装的连续性。

参考资料

https://github.com/chhylp123/hifiasm

https://hifiasm.readthedocs.io/en/latest/index.html

https://www.nature.com/articles/s41592-020-01056-5