macでインフォマティクス

macでインフォマティクス

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

ハイブリッドアセンブリにも対応したショートリードアセンブラ Unicycler

2019 1/11 追記

2019 4/12 dockerリンク追加

2019 5/27 help all追加

 

 現在、バクテリアゲノムはIlluminaシーケンシングプラットフォームによって支配されている。イルミナのリードは正確で、ベースあたりのコストが低く、全ゲノムシーケンシングの普及が可能である。しかし、多くのイルミナシーケンシングは、バクテリアゲノムの多くのリピートエレメントよりも短いフラグメント(500 bp以下)を使用している[ref,1]。これにより、ショートリードのアセンブリツール(アセンブラ)による完全なゲノムのresolutionを防ぎ、代わりに数十のアセンブリ(contigs)に断片化する。結果として、利用可能なバクテリアゲノムのほとんどは不完全であり、これは大規模な比較ゲノム研究を妨げる。

 Pacific Biosciences(PacBio)およびOxford Nanopore Technologies(ONT)シーケンシングプラットフォームは、10kbp以上のDNA断片をシーケンシングできるが、イルミナのプラットフォームよりも1塩基あたりのコストが高くなる。 PacBioおよびONTロングリードは、バクテリアゲノムアセンブリを完成させるにはしばしば十分だが、イルミナのリード(5-15%vs <1%)よりもずっと高いベースエラー率を持っている[ref.2,3]。したがって、ほとんどの研究者は、安価なIlluminaシーケンシングを使用して多くの分離株の断片化されたドラフトアセンブリを生成するか、高価なロングリード技術を使用して少ない分離株で完全なアセンブリを生成するかを選択する必要がある。ハイブリッドアセンブリは、代替手段としてショートリードとロングリードのコンビネーションを使用する。このアプローチでは、ショートリードは正確なコンティグを生成するために使用され、ロングリードはそれらを一緒にスキャフォールドするための情報を提供する。これは、完全なバクテリアゲノムへの最も費用効果の高い経路となり得る。

 近年のロングリード技術の進歩にもかかわらず、イルミナのリードりは、public health and research laboratoriesで広く使用されており[ref.4]、高精度で低コストであるためにしばらくの間続くだろう。さらに、イルミナのデータは既に数十万のバクテリア分離株で利用可能であり、これらの大部分はロングリードシーケンシングデータに置き換えられる可能性は低い。(一部略)

 ハイブリッドアセンブリは、ショートリード優先またはロングリード優先のいずれかの方法で達成することができる。ショートリード優先の方法では、スキャフォールディングツールにロングリードを使用してIlluminaコンティグを一緒に結合する。しかし、スキャフォールディングミスは一般的であり、シークエンスの構造上のミスマッチ(misassemblies)につながる[ref,5]。ロングリード優先アプローチは、訂正されていないロング・リードのアセンブリ、それに続くショートリードを用いたアセンブリエラー訂正を含む[ref.3]。あるいは、最初にショートリードを使用してロングリードのエラーを修正し、その後に修正されたロングリード[ref,6,7]を使ってアセンブリする。エラー訂正がアセンブリの前または後に実行されるかどうかにかかわらず、ロングリード優先アプローチは、ショートリード優先のアプローチよりもロングリードのデプスを必要とする。

ここでは、バクテリア分離株ゲノム用の新しいハイブリッドアセンブリパイプラインであるUnicyclerを紹介する。 Unicyclerは最初にショートリードを正確かつ接続されたアセンブリグラフ、つまりコンティグと相互接続の両方を含むデータ構造にアセンブルする[ref,8]。その後、ロングリードを使用して、グラフ内の最適なパスを検索する。ショートリード優先のアプローチに従うことによって、Unicyclerは少量のロングリードも有効に使用するが、ロングリードのデプスで十分であれば、完成したアセンブリ(1レプリコンあたり1コンティグ)を生成することができる。アセンブリグラフ接続を使用して可能なスキャホールディング構成を制約することにより、Unicyclerは、ショートリードファーストアセンブラよりもミスアセンブリ率が低くなる。

 

 

 

 Unicyclerはショートリードのde novo asssembly、ロングリードのde novo asssembly、ショートリードとロングリードを両方使ったハイブリッドのde novo asssemblyに対応したアセンブルパイプライン。単離されたバクテリア向けに設計されており、真核生物のゲノムやメタゲノムには対応していない。

 ショートリードからアセンブルのグラフを構築するために、内部でspadesを使用している。この時、Unicyclerはspadesのオプチマイザーとして機能する。具体的には、幅広いk-merを試し、dead-endとcontig数が両方とも最も少なくなる条件でアセンブルを行う。そのため、ショートリードだけのアセンブルでも、spadesよりパフォーマンスが上がる可能性がある。ロングリードだけを使ったアセンブルでは、内部でminiasmとraconが動く。miniasmでアセンブルすると、contigのエラー率は元のロングリードと同じくらいになるが、raconでくりかしpolishすることでエラーを取り除いている。

  Unicylcerの出力は、fastaとgfaファイルの両方で出力される。gfaはアセンブル情報を全て記載できるフォーマットであり、そのため、gfaの情報だけでbandageに読み込んでグラフパスを可視化することができる。また、Unicyclerはcontig末端の配列を刈り込むため、spadesのように末端にk-merサイズの重複ができることがない。そのためより正確にアセンブルの状態を把握することができる。

 

インストール

追記

2017 11/16 brewで導入できるようになった。

brew install unicycler

#Anaconda環境ならcondaを使う、例えば2018年12月現在最新の0.4.7
conda install -c bioconda unicycler==0.4.7

 

依存

  • Linux or macOS
  • Python 3.4 or later
  • C++ compiler with C++14 support:
    • GCC 4.9.1 or later
    • Clang 3.5 or later
    • ICC also works (though I don't know the minimum required version number)
  • setuptools (only required for installation of Unicycler)
  • For short-read or hybrid assembly:
    • SPAdes v3.6.2 or later (spades.py)
  • For long-read or hybrid assembly:
  • For polishing
  • For rotating circular contigs:
    • BLAST+ (makeblastdb and tblastn)

依存するものが多いが、spades、racon、pilon、samtools、bowtie2、blast+あたりはbrewで導入できる。raconはbrewで導入できるバージョンが古いので、gitから 直接インストールする。

git clone https://github.com/isovic/racon.git && cd racon && make modules && make tools && make -j 

 

python3.4環境を想定しているので、python2環境しかない人は、pyenvなどでpythonのバージョンを切り替えできるようにしておくと楽である。GCC4.9.1以降を使うため。、mac OSX標準のLLVMだとビルドできなない。一時的にGCCに切り替えてビルドした(参考にしたページ)。GCCのバージョンは

gcc --version

でチェックできる。

 

本体 Github

git clone https://github.com/rrwick/Unicycler.git
cd Unicycler
python3 setup.py install #GCCに切り替えてビルド

macOS10.11、10.12に加え、ubuntu14.04でもビルドできることを確認した(python3.4.0)。

 

テストランして、依存が全て揃っているか確認する。

f:id:kazumaxneo:20171026133944j:plain

ubuntuではjavaが古いと言われた。ここを参考にjava8にアップデートするとgoodになった。

>unicycler -h

$ unicycler -h

usage: unicycler [-h] [--help_all] [--version] [-1 SHORT1] [-2 SHORT2] [-s UNPAIRED] [-l LONG] -o OUT [--verbosity VERBOSITY] [--min_fasta_length MIN_FASTA_LENGTH] [--keep KEEP] [-t THREADS] [--mode {conservative,normal,bold}]

                 [--linear_seqs LINEAR_SEQS] [--vcf]

 

       __

       \ \___

        \ ___\

        //

   ____//      _    _         _                     _

 //_  //\\    | |  | |       |_|                   | |

//  \//  \\   | |  | | _ __   _   ___  _   _   ___ | |  ___  _ __

||  (O)  ||   | |  | || '_ \ | | / __|| | | | / __|| | / _ \| '__|

\\    \_ //   | |__| || | | || || (__ | |_| || (__ | ||  __/| |

 \\_____//     \____/ |_| |_||_| \___| \__, | \___||_| \___||_|

                                        __/ |

                                       |___/

 

Unicycler: an assembly pipeline for bacterial genomes

 

Help:

  -h, --help                           Show this help message and exit

  --help_all                           Show a help message with all program options

  --version                            Show Unicycler's version number

 

Input:

  -1 SHORT1, --short1 SHORT1           FASTQ file of first short reads in each pair (required)

  -2 SHORT2, --short2 SHORT2           FASTQ file of second short reads in each pair (required)

  -s UNPAIRED, --unpaired UNPAIRED     FASTQ file of unpaired short reads (optional)

  -l LONG, --long LONG                 FASTQ or FASTA file of long reads (optional)

 

Output:

  -o OUT, --out OUT                    Output directory (required)

  --verbosity VERBOSITY                Level of stdout and log file information (default: 1)

                                         0 = no stdout, 1 = basic progress indicators, 2 = extra info, 3 = debugging info

  --min_fasta_length MIN_FASTA_LENGTH  Exclude contigs from the FASTA file which are shorter than this length (default: 100)

  --keep KEEP                          Level of file retention (default: 1)

                                         0 = only keep final files: assembly (FASTA, GFA and log), 1 = also save graphs at main checkpoints, 2 = also keep SAM (enables fast rerun in different mode),

                                         3 = keep all temp files and save all graphs (for debugging)

  --vcf                                Produce a VCF by mapping the short reads to the final assembly (experimental, default: do not produce a vcf file)

 

Other:

  -t THREADS, --threads THREADS        Number of threads used (default: 8)

  --mode {conservative,normal,bold}    Bridging mode (default: normal)

                                         conservative = smaller contigs, lowest misassembly rate

                                         normal = moderate contig size and misassembly rate

                                         bold = longest contigs, higher misassembly rate

  --linear_seqs LINEAR_SEQS            The expected number of linear (i.e. non-circular) sequences in the underlying sequence (default: 0)

help all

> unicycler --help_all 

$ unicycler --help_all 

usage: unicycler [-h] [--help_all] [--version] [-1 SHORT1] [-2 SHORT2] [-s UNPAIRED] [-l LONG] -o OUT [--verbosity VERBOSITY] [--min_fasta_length MIN_FASTA_LENGTH] [--keep KEEP] [-t THREADS] [--mode {conservative,normal,bold}]

                 [--min_bridge_qual MIN_BRIDGE_QUAL] [--linear_seqs LINEAR_SEQS] [--min_anchor_seg_len MIN_ANCHOR_SEG_LEN] [--spades_path SPADES_PATH] [--no_correct] [--min_kmer_frac MIN_KMER_FRAC] [--max_kmer_frac MAX_KMER_FRAC] [--kmers KMERS]

                 [--kmer_count KMER_COUNT] [--depth_filter DEPTH_FILTER] [--spades_tmp_dir SPADES_TMP_DIR] [--no_miniasm] [--racon_path RACON_PATH] [--existing_long_read_assembly EXISTING_LONG_READ_ASSEMBLY] [--no_rotate]

                 [--start_genes START_GENES] [--start_gene_id START_GENE_ID] [--start_gene_cov START_GENE_COV] [--makeblastdb_path MAKEBLASTDB_PATH] [--tblastn_path TBLASTN_PATH] [--no_pilon] [--bowtie2_path BOWTIE2_PATH]

                 [--bowtie2_build_path BOWTIE2_BUILD_PATH] [--samtools_path SAMTOOLS_PATH] [--pilon_path PILON_PATH] [--java_path JAVA_PATH] [--min_polish_size MIN_POLISH_SIZE] [--vcf] [--bcftools_path BCFTOOLS_PATH]

                 [--min_component_size MIN_COMPONENT_SIZE] [--min_dead_end_size MIN_DEAD_END_SIZE] [--contamination CONTAMINATION] [--scores SCORES] [--low_score LOW_SCORE]

 

       __

       \ \___

        \ ___\

        //

   ____//      _    _         _                     _

 //_  //\\    | |  | |       |_|                   | |

//  \//  \\   | |  | | _ __   _   ___  _   _   ___ | |  ___  _ __

||  (O)  ||   | |  | || '_ \ | | / __|| | | | / __|| | / _ \| '__|

\\    \_ //   | |__| || | | || || (__ | |_| || (__ | ||  __/| |

 \\_____//     \____/ |_| |_||_| \___| \__, | \___||_| \___||_|

                                        __/ |

                                       |___/

 

Unicycler: an assembly pipeline for bacterial genomes

 

Help:

  -h, --help                            Show this help message and exit

  --help_all                            Show a help message with all program options

  --version                             Show Unicycler's version number

 

Input:

  -1 SHORT1, --short1 SHORT1            FASTQ file of first short reads in each pair (required)

  -2 SHORT2, --short2 SHORT2            FASTQ file of second short reads in each pair (required)

  -s UNPAIRED, --unpaired UNPAIRED      FASTQ file of unpaired short reads (optional)

  -l LONG, --long LONG                  FASTQ or FASTA file of long reads (optional)

 

Output:

  -o OUT, --out OUT                     Output directory (required)

  --verbosity VERBOSITY                 Level of stdout and log file information (default: 1)

                                          0 = no stdout, 1 = basic progress indicators, 2 = extra info, 3 = debugging info

  --min_fasta_length MIN_FASTA_LENGTH   Exclude contigs from the FASTA file which are shorter than this length (default: 100)

  --keep KEEP                           Level of file retention (default: 1)

                                          0 = only keep final files: assembly (FASTA, GFA and log), 1 = also save graphs at main checkpoints, 2 = also keep SAM (enables fast rerun in different mode),

                                          3 = keep all temp files and save all graphs (for debugging)

  --vcf                                 Produce a VCF by mapping the short reads to the final assembly (experimental, default: do not produce a vcf file)

 

Other:

  -t THREADS, --threads THREADS         Number of threads used (default: 8)

  --mode {conservative,normal,bold}     Bridging mode (default: normal)

                                          conservative = smaller contigs, lowest misassembly rate

                                          normal = moderate contig size and misassembly rate

                                          bold = longest contigs, higher misassembly rate

  --min_bridge_qual MIN_BRIDGE_QUAL     Do not apply bridges with a quality below this value

                                          conservative mode default: 25.0

                                          normal mode default: 10.0

                                          bold mode default: 1.0

  --linear_seqs LINEAR_SEQS             The expected number of linear (i.e. non-circular) sequences in the underlying sequence (default: 0)

  --min_anchor_seg_len MIN_ANCHOR_SEG_LEN

                                        If set, Unicycler will not use segments shorter than this as scaffolding anchors (default: automatic threshold)

 

SPAdes assembly:

  These options control the short-read SPAdes assembly at the beginning of the Unicycler pipeline.

 

  --spades_path SPADES_PATH             Path to the SPAdes executable (default: spades.py)

  --no_correct                          Skip SPAdes error correction step (default: conduct SPAdes error correction)

  --min_kmer_frac MIN_KMER_FRAC         Lowest k-mer size for SPAdes assembly, expressed as a fraction of the read length (default: 0.2)

  --max_kmer_frac MAX_KMER_FRAC         Highest k-mer size for SPAdes assembly, expressed as a fraction of the read length (default: 0.95)

  --kmers KMERS                         Exact k-mers to use for SPAdes assembly, comma-separated (example: 22,33,44, default: automatic)

  --kmer_count KMER_COUNT               Number of k-mer steps to use in SPAdes assembly (default: 10)

  --depth_filter DEPTH_FILTER           Filter out contigs lower than this fraction of the chromosomal depth, if doing so does not result in graph dead ends (default: 0.25)

  --spades_tmp_dir SPADES_TMP_DIR       Specify SPAdes temporary directory using the SPAdes --tmp-dir option (default: make a temporary directory in the output directory)

 

miniasm+Racon assembly:

  These options control the use of miniasm and Racon to produce long-read bridges.

 

  --no_miniasm                          Skip miniasm+Racon bridging (default: use miniasm and Racon to produce long-read bridges)

  --racon_path RACON_PATH               Path to the Racon executable (default: racon)

  --existing_long_read_assembly EXISTING_LONG_READ_ASSEMBLY

                                        A pre-prepared long read assembly for the sample in GFA format. If this option is used, Unicycler will skip the miniasm/Racon steps and instead use the given assembly (default: perform long read assembly

                                        using miniasm/Racon)

 

Assembly rotation:

  These options control the rotation of completed circular sequence near the end of the Unicycler pipeline.

 

  --no_rotate                           Do not rotate completed replicons to start at a standard gene (default: completed replicons are rotated)

  --start_genes START_GENES             FASTA file of genes for start point of rotated replicons (default: start_genes.fasta)

  --start_gene_id START_GENE_ID         The minimum required BLAST percent identity for a start gene search (default: 90.0)

  --start_gene_cov START_GENE_COV       The minimum required BLAST percent coverage for a start gene search (default: 95.0)

  --makeblastdb_path MAKEBLASTDB_PATH   Path to the makeblastdb executable (default: makeblastdb)

  --tblastn_path TBLASTN_PATH           Path to the tblastn executable (default: tblastn)

 

Pilon polishing:

  These options control the final assembly polish using Pilon at the end of the Unicycler pipeline.

 

  --no_pilon                            Do not use Pilon to polish the final assembly (default: Pilon is used)

  --bowtie2_path BOWTIE2_PATH           Path to the bowtie2 executable (default: bowtie2)

  --bowtie2_build_path BOWTIE2_BUILD_PATH

                                        Path to the bowtie2_build executable (default: bowtie2-build)

  --samtools_path SAMTOOLS_PATH         Path to the samtools executable (default: samtools)

  --pilon_path PILON_PATH               Path to a Pilon executable or the Pilon Java archive file (default: pilon)

  --java_path JAVA_PATH                 Path to the java executable (default: java)

  --min_polish_size MIN_POLISH_SIZE     Contigs shorter than this value (bp) will not be polished using Pilon (default: 10000)

 

VCF:

  These options control the production of the VCF of the final assembly.

 

  --bcftools_path BCFTOOLS_PATH         Path to the bcftools executable (default: bcftools)

 

Graph cleaning:

  These options control the removal of small leftover sequences after bridging is complete.

 

  --min_component_size MIN_COMPONENT_SIZE

                                        Graph components smaller than this size (bp) will be removed from the final graph (default: 1000)

  --min_dead_end_size MIN_DEAD_END_SIZE

                                        Graph dead ends smaller than this size (bp) will be removed from the final graph (default: 1000)

 

Long read alignment:

  These options control the alignment of long reads to the assembly graph.

 

  --contamination CONTAMINATION         FASTA file of known contamination in long reads

  --scores SCORES                       Comma-delimited string of alignment scores: match, mismatch, gap open, gap extend (default: 3,-6,-5,-2)

  --low_score LOW_SCORE                 Score threshold - alignments below this are considered poor (default: set threshold automatically)

 

追記

dockerイメージ

https://hub.docker.com/r/staphb/unicycler/

 

実行方法

ショートリードのアセンブリ

unicycler -1 R1.fastq.gz -2 R2.fastq.gz -o output_dir -t 12

 もちろん非圧縮のfastqも使用できる。

 

 ハイブリッドアセンブリ

unicycler -1 short_reads_1.fastq.gz -2 short_reads_2.fastq.gz -l long_reads.fastq.gz -o output_dir

 

 アセンブルの精度とcontiguityは‑‑modeオプションで設定できる。defaultではブリッジのクオリティの閾値 (cutoffを決める値)は10だが、--mode conservativeでは25になる。同時にショートリードをブリッジ情報として使われなくなる。--mode boldでは反対にcutoffの値が最低の1まで下がり、ショートリードもブリッジ情報として利用される。エラーのリスクも出てくるが、アセンブルをできる限り進めるならこのモードも検討すると良いかもしれない。

 

アセンブルFASTAとgfaファイルは必ず出力されるが、他の関連ファイルの数はkeepフラグの数値の大小で調整できる。--keep 0で関連ファイルはほぼ残らなくなる。数値をあげると関連ファイルも残るようになり、--keep 3で中間ファイルも含め全てのファイルが出力されるようになる。詳細はGitで確認してください。

 

追記

以下のまとめが大変参考になります。raconやnanopolish、pilonの単独、または併用polishよりunicyclerの方が良好な結果を出していますね。

What is a good genome assembly? – Albertsen Lab

ただし2年前の情報なので、現在改めて行ったらもう少しロングリード単独のアセンブリもよくなっているのではないでしょうか。

 

引用

Unicycler: Resolving bacterial genome assemblies from short and long sequencing reads

Ryan R. Wick , Louise M. Judd, Claire L. Gorrie, Kathryn E. Holt

PLoS Comput Biol. 2017 Jun 8;13(6):e1005595.