軟件安裝
官方提供了編譯好的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有如下作用
- 對初步組裝進(jìn)行polish
- 尋找同一物種不同株系間的變異,包括結(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.gz和read_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)容,基本上選snps和indels就夠了 -
--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
參考資料