macでインフォマティクス

macでインフォマティクス

NGS関連のインフォマティクス情報についてまとめています。

アセンブリのエラーやギャップ(NNN)を検出し、ポリッシュしたFASTAを出力するPilon

2018 8/31 タイトルと紹介文修正

2018 11/5 タイトル修正

2019 1/11  追記

2019 3/3ラストにnanopore long read追記

2019 4/12 ラストにpacbio long read追記

 

 何百というバクテリアゲノムをシーケンスできる時代になったが、それに伴い大量のデータを効率よく分析できる堅牢でスケール変化に対応できる手法が求められている。Pilonは全自動で動作するバクテリア向けのリファレンスベースのアセンブルツールである。draftゲノムかfinishゲノムのリファレンスとそのbamを指定することで、リファレンスと構造的に異なる部位を識別し、解析株のゲノムをリファンレンスから再構築して出力する。この時、構造が異なった部位はバリアントして出力する。リファレンスベースでアセンブルを行なっているため、SNVやsmall indelのみならず、リードの長さ以上の構造変化も検出することができる。また、ドラフトゲノムのポリッシュにも使われる。ポリッシュの過程で、ミスマッチ、indel、SV、scaffoldsのNNN...などが修復される。全てのプロセスは自動化されているため、ユーザーはFASTAとbamを指定するだけで、解析データのゲノムを得ることが可能になっている。2014年に論文が発表された。

  

公式サイト

http://software.broadinstitute.org/software/pilon/

Pilonに関するツイート

 

ダウンロード

Github

リリース

f:id:kazumaxneo:20171001220159j:plain

Javaのコンパリル済みバイナリも配布されている。バイナリをダウンロードして、実行権をつける。

#1.23をダウンロード(2019 3/3時点最新)
wget https://github.com/broadinstitute/pilon/releases/download/v1.23/pilon-1.23.jar

chmod u+x pilon-1.23.jar

#biocondaのサポートがあるのでcondaで導入できますが、エラーが起きた時はjavaの権限でダウンロードした.jarの実行ファイルを直接叩いてください。

 > java -jar pilon-1.22.jar --help

$ java -jar pilon-1.22.jar --help

Pilon version 1.22 Wed Mar 15 16:38:30 2017 -0400

 

    Usage: pilon --genome genome.fasta [--frags frags.bam] [--jumps jumps.bam] [--unpaired unpaired.bam]

                 [...other options...]

           pilon --help for option details 

 

 

         INPUTS:

           --genome genome.fasta

              The input genome we are trying to improve, which must be the reference used

              for the bam alignments.  At least one of --frags or --jumps must also be given.

           --frags frags.bam

              A bam file consisting of fragment paired-end alignments, aligned to the --genome

              argument using bwa or bowtie2.  This argument may be specifed more than once.

           --jumps jumps.bam

              A bam file consisting of jump (mate pair) paired-end alignments, aligned to the

              --genome argument using bwa or bowtie2.  This argument may be specifed more than once.

           --unpaired unpaired.bam

              A bam file consisting of unpaired alignments, aligned to the --genome argument 

              using bwa or bowtie2.  This argument may be specifed more than once.

           --bam any.bam

              A bam file of unknown type; Pilon will scan it and attempt to classify it as one

              of the above bam types.

         OUTPUTS:

           --output prefix

              Prefix for output files

           --outdir directory

              Use this directory for all output files.

           --changes

              If specified, a file listing changes in the <output>.fasta will be generated.

           --vcf

              If specified, a vcf file will be generated

           --vcfqe

               If specified, the VCF will contain a QE (quality-weighted evidence) field rather

               than the default QP (quality-weighted percentage of evidence) field.

           --tracks

               This options will cause many track files (*.bed, *.wig) suitable for viewing in

               a genome browser to be written.

         CONTROL:

           --variant

              Sets up heuristics for variant calling, as opposed to assembly improvement;

              equivalent to "--vcf --fix all,breaks".

           --chunksize

              Input FASTA elements larger than this will be processed in smaller pieces not to

              exceed this size (default 10000000).

           --diploid

              Sample is from diploid organism; will eventually affect calling of heterozygous SNPs

           --fix fixlist

              A comma-separated list of categories of issues to try to fix:

                "snps": try to fix individual base errors;

                "indels": try to fix small indels;

                "gaps": try to fill gaps;

                "local": try to detect and fix local misassemblies;

                "all": all of the above (default);

                "bases": shorthand for "snps" and "indels" (for back compatibility);

                "none": none of the above; new fasta file will not be written.

              The following are experimental fix types:

                "amb": fix ambiguous bases in fasta output (to most likely alternative);

                "breaks": allow local reassembly to open new gaps (with "local");

                "circles": try to close circlar elements when used with long corrected reads;

                "novel": assemble novel sequence from unaligned non-jump reads.

           --dumpreads

              Dump reads for local re-assemblies.

           --duplicates

              Use reads marked as duplicates in the input BAMs (ignored by default).

           --iupac

              Output IUPAC ambiguous base codes in the output FASTA file when appropriate.

           --nonpf

              Use reads which failed sequencer quality filtering (ignored by default).

           --targets targetlist

              Only process the specified target(s).  Targets are comma-separated, and each target

              is a fasta element name optionally followed by a base range.

              Example: "scaffold00001,scaffold00002:10000-20000" would result in processing all of

              scaffold00001 and coordinates 10000-20000 of scaffold00002.

              If "targetlist" is the name of a file, each line will be treated as a target

              specification.

           --threads

              Degree of parallelism to use for certain processing (default 1). Experimental.

           --verbose

              More verbose output.

           --debug

              Debugging output (implies verbose).

           --version

              Print version string and exit.

         HEURISTICS:

           --defaultqual qual

              Assumes bases are of this quality if quals are no present in input BAMs (default 15).

           --flank nbases

              Controls how much of the well-aligned reads will be used; this many bases at each

              end of the good reads will be ignored (default 10).

           --gapmargin

              Closed gaps must be within this number of bases of true size to be closed (100000)

           --K

              Kmer size used by internal assembler (default 47).

           --mindepth depth

              Variants (snps and indels) will only be called if there is coverage of good pairs

              at this depth or more; if this value is >= 1, it is an absolute depth, if it is a

              fraction < 1, then minimum depth is computed by multiplying this value by the mean

              coverage for the region, with a minumum value of 5 (default 0.1: min depth to call 

              is 10% of mean coverage or 5, whichever is greater).

           --mingap

              Minimum size for unclosed gaps (default 10)

           --minmq

              Minimum alignment mapping quality for a read to count in pileups (default 0)

           --minqual

              Minimum base quality to consider for pileups (default 0)

           --nostrays

              Skip making a pass through the input BAM files to identify stray pairs, that is,

              those pairs in which both reads are aligned but not marked valid because they have

              inconsistent orientation or separation. Identifying stray pairs can help fill gaps

              and assemble larger insertions, especially of repeat content.  However, doing so

              sometimes consumes considerable memory.

 

ラン

リファレンス(近縁種でも可能)と、そのbamファイルが必要である。mate-pairやシングルエンドのbamにも対応している。

java -Xmx16G -jar pilon-1.23.jar --genome ref.fa --frags pair.bam --vcf --tracks --changes --verbose --threads 8 --outdir output_dir

INPUTS:

  • --genome genome.fasta The input genome we are trying to improve, which must be the reference used for the bam alignments. At least one of --frags or --jumps must also be given.
  • --frags frags.bam A bam file consisting of fragment paired-end alignments, aligned to the --genome argument using bwa or bowtie2. This argument may be specifed more than once.
  • --jumps A bam file consisting of jump (mate pair) paired-end alignments, aligned to the --genome argument using bwa or bowtie2. This argument may be specifed more than once.
  • --unpaired A bam file consisting of unpaired alignments, aligned to the --genome argument using bwa or bowtie2. This argument may be specifed more than once.
  • --bam  A bam file of unknown type; Pilon will scan it and attempt to classify it as one of the above bam types.

OUTPUTS:

  • --output Prefix for output files
  • --outdir Use this directory for all output files.
  • --changes If specified, a file listing changes in the <output>.fasta will be generated.
  • --vcf  If specified, a vcf file will be generated
  • --vcfqe If specified, the VCF will contain a QE (quality-weighted evidence) field rather than the default QP (quality-weighted percentage of evidence) field.
  • --tracks This options will cause many track files (*.bed, *.wig) suitable for viewing in a genome browser to be written.

 解析が終わると、指定したディレクトリにFASTAと、リファレンスと異なる部位を示したファイル(.change)が出力される。また、--taracksもつけていれば、異なる部位を記載したbedやGCカバレッジ、copy number variationなどのファイルも出力される。

変更箇所が記載されたpilon.cganges

f:id:kazumaxneo:20180831133820j:plain

 

出力解説

https://github.com/broadinstitute/pilon/wiki/Output-File-Descriptions

 

追記

 de novo assemblyで得たcontigを使えば、polishも行える。

まずショートリードをcontigにマッピングする。

#マッピング
bwa mem -t 20 contig.fasta illumina_read1.fq illumina_read2.fq |samtools view -@ 20 -bS - > aln.bam

#ソート
samtools sort -@ 20 aln.bam > aln_sorted.bam

#index
samtools index aln_sorted.bam

#ポリッシュ
java -Xmx16G -jar pilon-1.23.jar --genome racon_polished3.fasta --frags aln_sorted.bam --vcf --changes --verbose --threads 20 --outdir polish_dir

 

 

 

tracks情報をigvで確認する。

cd directory/
igv -g ../ref.fa pair.bam,,pilonBadCoverage.wig,pilonChanges.wig,pilonClippedAlignments.wig,pilonCopyNumber.wig,pilonCoverage.wig,pilonDeltaCoverage.wig,pilonDipCoverage.wig,pilonGC.wig,pilonPctBad.wig,pilonPhysicalCoverage.wig,pilonPilon.bed,pilonUnconfirmed.wig,pilonWeightedMq.wig,pilonWeightedQual.wig &

 リファレンスは出力ではなく、参照したレファンレンスのfastaの方を指定する。

 

.changesファイルの部位をIGVの窓にコピペすれば移動が速い(IGVフォーマットに対応している)。

SNV部位

f:id:kazumaxneo:20171001224403j:plain

 

large Insertion部位

f:id:kazumaxneo:20171001224207j:plain

 

 large deletion部位

f:id:kazumaxneo:20171001224318j:plain

 

copy number variation部位

f:id:kazumaxneo:20171001224507j:plain

ただしこの部位は検出できていない。

 

 

ゲノム構造の違いも確認してみる。

prokkaでアノテーションをつける。

prokka pilon.fasta
prokka ref.fa     #originalもgbkを作る。

gbkファイルが出力されたら、これを使いBryan Weeさんのラッパーツールでblast解析(リンク

bwast.py ref.gbk pilon.gbk -a
  •  -a Run ACT after performing BLAST

ACTが自動起動する。

f:id:kazumaxneo:20171001233004j:plain

今回のデータは欠損が真ん中付近に落ちている。

f:id:kazumaxneo:20171001233252j:plain

 

ACTはコード領域のprotein情報を見ているので、ゲノムを直接比較するならMauveを使う。

修正 ACTはblastnの結果でも比較できます。ラッパーツールを使うなら

bwast.py ref.fa pilonA.fa pilonB.fa -a

という感じでやってみてください(リンク)。

 

他にも k-merの値、領域、差分、欠損のcoverage閾値など多様なオプションがある。詳細は公式サイトで確認してください。

https://github.com/broadinstitute/pilon/wiki/Requirements-&-Usage

 

主観ですがエンドユーザにとって良いツールだと思います。 ただし、この手法で認識されなかった変異はアセンブルされたFASTAに再現されないことは注意してください。複雑な構造変化の検出は難しいでしょう。

 

追記

他のツールと比較されています。とても参考になります。

What is a good genome assembly? – Albertsen Lab

2019 4/12追記

v1.23にはlong readを使ったpolishオプションが実装されています。

ショートリードのbamとロングリードのbamを両方指定する。#vcfは遅くなるので消す。

nanopre

java -Xmx24G -jar pilon-1.23.jar --genome assembly.fasta --frags short_resd.bam --nanopore long_read.bam --tracks --changes --verbose --threads 20 --outdir pilon_output_dir

pacbio

java -Xmx24G -jar pilon-1.23.jar --genome assembly.fasta --frags short_resd.bam --pacbio long_read.bam --tracks --changes --verbose --threads 20 --outdir pilon_output_dir

 

テストした限り、ショトーリードもロングリードも両方あった方がpolishできる部位はかなり増えます。ロングリード情報があった方がNNNの解消にも有利なはずです。

引用

Pilon: An Integrated Tool for Comprehensive Microbial Variant Detection and Genome Assembly Improvement

Bruce J. Walker , Thomas Abeel , Terrance Shea, Margaret Priest, Amr Abouelliel, Sharadha Sakthikumar, Christina A. Cuomo, Qiandong Zeng, Jennifer Wortman, Sarah K. Young, Ashlee M. Earl

PLoS One. 2014 Nov 19;9(11):e112963.