之前那我由于需要對(duì)PacBio的組裝結(jié)果進(jìn)行polish，于是寫了「三代組裝」Pacbio組裝后如何用自身數(shù)據(jù)進(jìn)行polish。最近發(fā)現(xiàn)自己又有了需求，于是重新回顧了我之前寫的這篇文章，但是在實(shí)踐的時(shí)候，卻發(fā)現(xiàn)新版本pb-assembly(0.0.8)沒(méi)有之前的pbalign和variantCaller這兩個(gè)命令，也就是意味著我之前的文章隨著版本更新而失效了。

此時(shí)我就面臨著兩個(gè)選擇，一個(gè)用conda安裝之前的版本，逃避可恥但有用。另一個(gè)就是探索新版本下的代碼，勇敢面對(duì)慘烈的人生。由于時(shí)間不是特別的緊張，所以我選擇了第二條路。

通過(guò)我一波檢索，我發(fā)現(xiàn)第二條路是一條沒(méi)有人走過(guò)的路。https://github.com/PacificBiosciences/pb-assembly只講了falcon組裝后，使用falcon-unzip進(jìn)行polish，而沒(méi)有介紹使用其他軟件(如Canu, Mecat2) 應(yīng)該怎么做。https://github.com/PacificBiosciences/GenomicConsensus提到的quiver命令并不存在于我們pb-assembly中。

我雖然把pb-assembly下的命令都看了個(gè)遍，但是也沒(méi)有想到哪些命令可以組合在一起用來(lái)polish。但是我突然想到，既然falcon-unzip里面會(huì)有polish，那么的它運(yùn)行日志里面一定會(huì)有它如何進(jìn)行polish蛛絲馬跡。經(jīng)過(guò)一波的分析和追蹤，終于被我找到了相關(guān)的代碼。修改之后如下

set -e
REF=$1
BAM=$2
THREADS=80

# align
pbmm2 align --sort  $REF $BAM aln.bam
pbindex aln.bam
# arrow
gcpp --algorithm=arrow -x 5 -X 120 -q 0 -j $THREADS -r $REF  aln.bam -o cns.fa,cns.fq,cns.vcf

新的腳本不再用BLASR進(jìn)行比對(duì)，而是擁抱了minimap2，提高了比對(duì)速度，并且避免了之前BLASR可能會(huì)出現(xiàn)的內(nèi)存泄露問(wèn)題。使用gcpp替代了之前variantCaller，但是參數(shù)幾乎沒(méi)有變動(dòng)，你可以認(rèn)為就是改了一個(gè)名字而已。

另外這是一個(gè)簡(jiǎn)化的腳本，如果是多個(gè)bam文件，則需要分別比對(duì)然后合并。原來(lái)的falcon流程會(huì)將參考基因組按照染色體拆分然后并行處理，但是我們的腳本設(shè)置了更高的線程數(shù)，最后結(jié)果速度也差不多。