Pilon | 利用二代测序数据优化三代测序数据组装结果

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

# 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

# 下载二代测序数据用于纠错
wget \
-O illumina.sra \
https://sra-pub-run-odp.s3.amazonaws.com/sra/SRR8482586/SRR8482586

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

# 输入参数
--genome : 输入基因组序列文件，即需要纠错的文件；
--frags : 输入Illumina小片段文库（FR方向）测序数据比对到参考基因上的BAM文件路径。该BAM文件是需要按coordinate排序，且具有.bai索引；
--jumps ： 输入Illumina大片段文库（RF方向）测序数据比对到参考基因上的BAM文件路径。该BAM文件是需要按coordinate排序，且具有.bai索引；
--unpaired ： 输入Illumina单端测序数据比对到参考基因上的BAM文件路径。该BAM文件是需要按coordinate排序，且具有.bai索引；
--bam ： 输入未知类型的Illumina测序数据比对到参考基因上的BAM文件路径。该BAM文件是需要coordinate排序，且具有.bai索引。程序能自动识别测序文库类型。
# 输出参数
--output : 设置输出文件前缀；
--outdir : 设置输出文件夹路径；
--changes ： 加入该参数会生成一个包含基因组序列改变信息的文本文件；
--vcf : 加入该参数后，会生成一个vcf格式文件；
# 运行参数
--fix : 该参数用于设定pilon对基因组进行指定类型的修正，该参数的值是用逗号分隔的多个单词。pilon可以对基因组进行多种类型的修正。
--threads : 设置线程数

# 对拼接结果建立索引（如何获得assembly.fasta详见Flye三代基因组推文）
bwa index assembly.fasta
# illumina与assembly.fasta进行比对，生成assembly_illumina.sam结果文件
bwa mem -t 12 assembly.fasta  illumina_1.fastq.gz illumina_2.fastq.gz > assembly_illumina.sam
# 将assembly_illumina.sam进行排序，生成assembly_illumina.sorted.bam 文件
samtools sort -@ 12 -O bam -o assembly_illumina.sorted.bam assembly_illumina.sam
# 对排序后的文件建立索引，生成assembly_illumina.sorted.bam.bai文件
samtools index assembly_illumina.sorted.bam
# 运行Pilon
pilon --genome assembly.fasta \
--fix all  \
--changes  \
--frags assembly_illumina.sorted.bam  \
--output pilon  \
--outdir pilon_result \
--vcf

java -jar  \
/opt/biosoft/GENOME/pilon/pilon-1.23.jar  \
--genome assembly.fasta  \
--fix all  \
--changes  \
--frags assembly_illumina.sorted.bam  \
--output pilon  \
--outdir pilon_result   \
--vcf

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/miniconda3/envs/NGS/bin/pilon

pilon.changes  # 该文件列出了纠错的位点
pilon.fasta  # 最终纠错后文件

# 有多少行代表有多少错误被纠正
wc -l  pilon.changes
# 统计纠错前后文件信息
seqkit stats pilon.fasta assembly.fasta