前言

三代測序錯誤率比較高，一般組裝后需要進(jìn)行糾錯來提高準(zhǔn)確度。本次介紹使用Pilon通過引入二代測序數(shù)據(jù)來對三代基因組進(jìn)行糾錯。

Pilon官網(wǎng)

https://github.com/broadinstitute/pilon/wiki

Pilon軟件安裝

#conda 安裝pilon
conda install -y pilon

#編譯安裝
wget https://github.com/broadinstitute/pilon/releases/download/v1.24/pilon-1.24.jar
chomd 755 pilon-1.24.jar

Pilon示例數(shù)據(jù)下載

#下載二代測序數(shù)據(jù)用于糾錯
wget \
-O illumina.sra \
https://sra-pub-run-odp.s3.amazonaws.com/sra/SRR8482586/SRR8482586

本期需要糾錯的基因組選擇上期推文中Flye組裝的nanopore數(shù)據(jù)進(jìn)行演示，即下文assembly.fasta

Pilon示例數(shù)據(jù)處理

fastq-dump --split-files --gzip illumina.sra

fastq-dump會將sra格式轉(zhuǎn)化成fastq格式，同時--gzip參數(shù)會對fastq進(jìn)行壓縮，示例illumina.sra最終會被轉(zhuǎn)化為illumina_1.fastq.gz 和 illumina_2.fastq.gz

Pilon常用參數(shù)

--genome ： 設(shè)置需要糾錯的基因組
--fix : 參數(shù)可選snps、indels、gaps、local、all等（默認(rèn)all）
--changes ： 列出糾錯位點(diǎn)
--frags ： 輸入paired-end比對文件（不同測序數(shù)據(jù)該選項(xiàng)不同，具體查看該軟件幫助文檔；若不知道，可直接使用--bam
--output ： 輸入結(jié)果前綴
--outdir ： 輸出文件
--vcf : 生成vcf格式文件

Pilon使用案例

示例使用的是conda安裝的Pilon

#對拼接結(jié)果建立索引（如何獲得assembly.fasta詳見Flye三代基因組推文）
bwa index assembly.fasta
#illumina與assembly.fasta進(jìn)行比對，生成assembly_illumina.sam結(jié)果文件
bwa mem -t 12 assembly.fasta  illumina_1.fastq.gz illumina_2.fastq.gz > assembly_illumina.sam
#將assembly_illumina.sam進(jìn)行排序，生成assembly_illumina.sorted.bam 
samtools sort -@ 12 -O bam -o assembly_illumina.sorted.bam assembly_illumina.sam
#運(yùn)行Pilon
pilon --genome assembly.fasta --fix all --changes --frags assembly_illumina.sorted.bam --output pilon --outdir pilon_result  --vcf

可能會遇到下面的報(bào)錯信息，這是由于軟件設(shè)定的內(nèi)存不足造成的

Exception in thread "main" java.lang.OutOfMemoryError: Java heap space
    at org.broadinstitute.pilon.BaseSum.<init>(BaseSum.scala:24)
    at org.broadinstitute.pilon.PileUp.<init>(PileUp.scala:27)
    at org.broadinstitute.pilon.PileUpRegion.$anonfun$new$1(PileUpRegion.scala:30)
    at org.broadinstitute.pilon.PileUpRegion$$Lambda$52/0x0000000100178840.apply$mcVI$sp(Unknown Source)
    at scala.collection.immutable.Range.foreach$mVc$sp(Range.scala:190)
    at org.broadinstitute.pilon.PileUpRegion.<init>(PileUpRegion.scala:30)
    at org.broadinstitute.pilon.GenomeRegion.initializePileUps(GenomeRegion.scala:150)
    at org.broadinstitute.pilon.GenomeFile.$anonfun$processRegions$4(GenomeFile.scala:104)
    at org.broadinstitute.pilon.GenomeFile.$anonfun$processRegions$4$adapted(GenomeFile.scala:102)
    at org.broadinstitute.pilon.GenomeFile$$Lambda$51/0x0000000100169840.apply(Unknown Source)
    at scala.collection.immutable.List.foreach(List.scala:333)
    at org.broadinstitute.pilon.GenomeFile.processRegions(GenomeFile.scala:102)
    at org.broadinstitute.pilon.Pilon$.main(Pilon.scala:111)
    at org.broadinstitute.pilon.Pilon.main(Pilon.scala)

解決辦法如下：

#查詢pilon路徑
which pilon
#修改pilon配置
vim /home/xiaoli/miniconda3/envs/NGS/bin/pilon

修改下圖紅色框，將 -Xmsg和-Xmx對應(yīng)的數(shù)值調(diào)大，再次運(yùn)行即可成功。

pilon debug.png

Pilon主要結(jié)果文件

pilon.changes  #該文件列出了糾錯的位點(diǎn)
pilon.fasta  #最終糾錯后文件

查看Pilon糾錯效果

#有多少行代表有多少錯誤被糾正
wc -l  pilon.changes
#統(tǒng)計(jì)糾錯前后文件信息
seqkit stats pilon.fasta assembly.fasta

PS.糾錯可以進(jìn)行多次，即：將第一次糾錯結(jié)果作為第二次需要糾錯的文件再次糾錯