「三代組裝」使用Pilon對基因組進(jìn)行polish

軟件安裝

官方提供了編譯好的jar包,方便使用

wget https://github.com/broadinstitute/pilon/releases/download/v1.23/pilon-1.23.jar
java -Xmx16G -jar pilon-1.23.jar

如果要順利運行程序,要求JAVA > 1.7, 以及根據(jù)基因組大小而定的內(nèi)存,一般而言是1M大小的基因?qū)?yīng)1GB的內(nèi)存。

總覽

Pilon有如下作用

  1. 對初步組裝進(jìn)行polish
  2. 尋找同一物種不同株系間的變異,包括結(jié)構(gòu)變異檢測

他以FASTA和BAM文件作為輸入,根據(jù)比對結(jié)果對輸入的參考基因組進(jìn)行提高,包括

  • 單堿基差異
  • 小的插入缺失(indels)
  • 較大的插入缺失或者block替換時間
  • 填充參考序列中的N
  • 找到局部的錯誤組裝

最后它輸出polish后的FASTA文件, 以及包含變異信息的VCF文件(可選)

分析流程

推薦使用PCR-free建庫測序得到的Illumina paired-end數(shù)據(jù),這樣子避免了PCR-duplication,有效數(shù)據(jù)更多,也不需要在分析過程中標(biāo)記重復(fù)。

下面步驟,假設(shè)你的組裝文件為draft.fa, 質(zhì)量控制后的illumina雙端測序數(shù)據(jù)分別為read_1.fq.gzread_2.fq.gz

第一步:比對

bwa index -p index/draft draft.fa
bwa mem -t 20 index/draft read_1.fq.gz read_2.fq.gz | samtools sort -@ 10 -O bam -o align.bam
samtools index -@ 10 align.bam

第二步:標(biāo)記重復(fù)(非PCR-free建庫)

sambamba markdup -t 10 align.bam align_markdup.bam

第三步:過濾高質(zhì)量比對的read

samtools view -@ 10 -q 30 align_markdup.bam > align_filter.bam
samtools index -@ 10 align_filter.bam

第三步:使用Pilon進(jìn)行polish

MEMORY= #根據(jù)基因組大小而定
java -Xmx${MEMORY}G -jar pilon-1.23.jar --genome draft.fa --frags align_filer.bam \
    --fix snps,indels \
    --output pilon_polished --vcf &> pilon.log

一般Pilon迭代個2到3次就夠了,所謂事不過三,過猶不及。

關(guān)于Pilon的一些參數(shù)說明:

  • --frags表示輸入的是1kb以內(nèi)的paired-end文庫,--jumps表示 大于1k以上的mate pair文庫, --bam則是讓軟件自己猜測
  • -vcf: 輸出一個vcf文件,包含每個堿基的信息
  • --fix: Pilon將會處理的內(nèi)容,基本上選snpsindels就夠了
  • --variant: 啟發(fā)式的變異檢測,等價于--vcf --fix all,breaks, 如果是polish不要使用該選項
  • minmq: 用于Pilon堆疊的read最低比對質(zhì)量,默認(rèn)是0。

閱讀日志輸出

這個日志文件是標(biāo)準(zhǔn)輸出而不是標(biāo)準(zhǔn)錯誤輸出,不過保險起見用&>

最開始,Pilon會輸出他的版本信息(如下示例),以及將會對基因組做的調(diào)整,

Pilon version 1.14 Sat Oct 31 14:30:00 2015 -0400
Genome: genome.fasta
Fixing snps, indels

其中Fixing后面的含義為:

  • "snps": 單堿基差異
  • "indels":小的indel的差異
  • "amb": 替換原有的N
  • "gaps": 填充基因組的gap
  • "local": 檢測和修改錯誤組裝
  • "all": 上述所有
  • "none": 不是上述的任何一種

接著Pilon會分染色體對BAM文件進(jìn)行處理,根據(jù)BAM文件進(jìn)行堆疊(pileup), 這個時候它會輸出有效reads的深度,這里的有效reads包括未成對的read和正確成對的read。

Processing ctg1:1-5414473
frags align_mkdup.bam: coverage 19
Total Reads: 808985, Coverage: 19, minDepth: 5

從Pilon v1.14開始,Pilon還會輸出基因組得到確認(rèn)的堿基比例。

Confirmed 5403864 of 5414473 bases (99.80%)

后續(xù)是Pilon將會對原參考基因組做的一些調(diào)整的總體情況,如下表示糾正2個snp, 2個小的插入,4個缺失。

Corrected 2 snps; 0 ambiguous bases; corrected 2 small insertions totaling 12 bases, 4 small deletions totaling 6 bases

最后聲明當(dāng)前部分處理結(jié)束

Finished processing ctg1:1-5414473

如果,在--fix中選了gaps, 那么輸出的內(nèi)容還有如下內(nèi)容。其中82048 -0 +276解釋為在坐標(biāo)82428處移除0個堿基,插入276個堿基。

# Attempting to fill gaps
fix gap: scaffold00001:82428-93547 82428 -0 +276 ClosedGap

參考資料

最后編輯于
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成,瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明:文章內(nèi)容(如有圖片或視頻亦包括在內(nèi))由作者上傳并發(fā)布,文章內(nèi)容僅代表作者本人觀點,簡書系信息發(fā)布平臺,僅提供信息存儲服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容