hifiasm大概是目前为止支持PacBio HiFi数据组装的所有软件中表现最优异的软件了。它不但能输出primary assembly result
,还能分别组装出单倍体(haplotype)
基因组的组装结果。
今天的分享就从hifiasm的组装开始,其他的组装软件之后也会介绍到。
可用资源
- 文章地址:https://www.nature.com/articles/s41592-020-01056-5
- GitHub地址:https://github.com/chhylp123/hifiasm
- 官方教程:https://hifiasm.readthedocs.io/en/latest/index.html
- 洲更学长的笔记:https://xuzhougeng.top/archives/assemble-hifi-pacbio-with-hifiasm
软件安装
hifiasm的安装也非常简单,conda一条命令就可以了
conda install hifiasm
运行hifiasm
纯PacBio HiFi组装
把命令写到一个名为run_hifiasm.sh
的脚本里
vim run_hifiasm.sh
只要填入最基本的
hifiasm \
-o test.hifisam \
-t 24 \
/path/to/PacBio/test1.hifi_reads.fastq.gz \
/path/to/PacBio/test2.hifi_reads.fastq.gz \
/path/to/PacBio/test3.hifi_reads.fastq.gz
用nohup挂在后台运行就可以啦。
nohup bash run_hifiasm.sh &
输出结果如下:
test.hifisam.ovlp.source.bin
test.hifisam.ovlp.reverse.bin
test.hifisam.ec.bin
test.hifisam.bp.r_utg.noseq.gfa
test.hifisam.bp.r_utg.lowQ.bed
test.hifisam.bp.r_utg.gfa
test.hifisam.bp.p_utg.noseq.gfa
test.hifisam.bp.p_utg.lowQ.bed
test.hifisam.bp.p_utg.gfa
test.hifisam.bp.p_ctg.noseq.gfa
test.hifisam.bp.p_ctg.lowQ.bed
★ test.hifisam.bp.p_ctg.gfa
test.hifisam.bp.hap2.p_ctg.noseq.gfa
test.hifisam.bp.hap2.p_ctg.lowQ.bed
★ test.hifisam.bp.hap2.p_ctg.gfa
test.hifisam.bp.hap1.p_ctg.noseq.gfa
test.hifisam.bp.hap1.p_ctg.lowQ.bed
★ test.hifisam.bp.hap1.p_ctg.gfa
上面打了五角星★
的三个文件是后续重点关注的文件。根据
https://hifiasm.readthedocs.io/en/latest/interpreting-output.html
的描述:
Hifiasm generates the following assembly graphs only with HiFi reads in default:
prefix.bp.p_ctg.gfa
: assembly graph of primary contigs.
prefix.bp.hap1.p_ctg.gfa
: partially phased contig graph of haplotype1.
prefix.bp.hap2.p_ctg.gfa
: partially phased contig graph of haplotype2.
HiC data integrated Integrated
如果你有HiC数据(测基因组的话大概率会有HiC数据的吧?)可以尝试基于HiC的证据的分型,primary的结果不会有什么变化,但是haplotype的结果会有提高。
hifiasm \
-o test.hifisam_hic \
-t 24 \
--h1 /path/to/HiC/1_R1.fastq.gz --h2 /path/to/HiC/1_R2.fastq.gz \
/path/to/PacBio/test1.hifi_reads.fastq.gz \
/path/to/PacBio/test2.hifi_reads.fastq.gz \
/path/to/PacBio/test3.hifi_reads.fastq.gz
格式转换
用gfatools或者自己写个脚本把gfa格式的文件转换成fasta就得到最终的组装结果啦。
如果写脚本转换
参考自:https://github.com/linsalrob/EdwardsLab/blob/master/bin/gfa2fasta.sh
(假设命名为gfa2fasta.sh
)
#!/usr/bin/env bash
if [ -z $1 ]; then
echo "Usage: $0 <gfa file> > <fasta file>"
exit $E_BADARGS
fi
awk -v id=$1 '/^S/{print ">"id"_"$2"\n"$3}' $1
运行:
bash gfa2fasta.sh test.hifisam.bp.hap1.p_ctg.gfa > test.hifisam.bp.hap1.p_ctg.fasta
用gfatools转换
先用conda安装gfatools
conda install gfatools
接下来也写一个如上的脚本方便使用(假设命名为gfa2fa.sh
)
#!/usr/bin/env bash
if [ -z $1 ]; then
echo "$0 <gfa file>"
echo "It will generate the fasta file like ${1%.*}.fasta"
exit $E_BADARGS
fi
gfatools gfa2fa ${1} > ${1%.*}.fasta
使用起来就很简单啦:
bash gfa2fa.sh test.hifisam.bp.hap1.p_ctg.gfa
私货时间
- 我曾经尝试过把多个PacBio HiFi的gz压缩文件直接用cat给连接成一个文件然后去组装,但是组装结果和保持多个文件一起输入给hifiasm的结果是不一样的。cat到一起的结果的N50要略低于把三个文件分别输给hifiasm的版本。我和老板冥思苦想了好一会儿没想明白原因在哪儿,也就作罢了。所以建议如果有多个文件,就保持多个文件,不需要合并成一个。反正hifiasm是可以接受输入多个文件的。
- 建议有HiC的情况下尽量加上HiC数据,结果会更好一些。
- PacBio HiFi的组装结果不需要再经过其他的软件polish,不像传统的三代组装软件那样流程冗长而繁琐。
Hifiasm does not need polishing tools like pilon or racon, either.
引用自:
(https://hifiasm.readthedocs.io/en/latest/index.html)