hifiasm大概是目前為止支持PacBio HiFi數(shù)據(jù)組裝的所有軟件中表現(xiàn)最優(yōu)異的軟件了。它不但能輸出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
上面打了五角星★的三個文件是后續(xù)重點關注的文件。根據(jù)
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數(shù)據(jù)（測基因組的話大概率會有HiC數(shù)據(jù)的吧？）可以嘗試基于HiC的證據(jù)的分型，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
格式轉(zhuǎn)換
用gfatools或者自己寫個腳本把gfa格式的文件轉(zhuǎn)換成fasta就得到最終的組裝結果啦。
如果寫腳本轉(zhuǎn)換
參考自：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轉(zhuǎn)換
先用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
私貨時間
- 我曾經(jīng)嘗試過把多個PacBio HiFi的gz壓縮文件直接用cat給連接成一個文件然后去組裝，但是組裝結果和保持多個文件一起輸入給hifiasm的結果是不一樣的。cat到一起的結果的N50要略低于把三個文件分別輸給hifiasm的版本。我和老板冥思苦想了好一會兒沒想明白原因在哪兒，也就作罷了。所以建議如果有多個文件，就保持多個文件，不需要合并成一個。反正hifiasm是可以接受輸入多個文件的。
- 建議有HiC的情況下盡量加上HiC數(shù)據(jù)，結果會更好一些。
- PacBio HiFi的組裝結果不需要再經(jīng)過其他的軟件polish，不像傳統(tǒng)的三代組裝軟件那樣流程冗長而繁瑣。
Hifiasm does not need polishing tools like pilon or racon, either.
引用自：
（https://hifiasm.readthedocs.io/en/latest/index.html）
