macでインフォマティクス

macでインフォマティクス

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

(メタゲノム向け)高効率なプロテインレベルのアセンブリツール PLASS

 

 メタゲノム研究の主な制限は、ショートリードの大部分(土壌で80% - 90%[1])を、遺伝子およびタンパク質配列の予測を可能にするのに十分な長さの連続した配列(contigs)にアセンブリすることができないことである。
 低存在量のゲノムはアセンブリが困難であり、アセンブリされていないリードは遺伝的多様性の大きな不均衡と恐らくはさらに大きな割合の生物学的新規性を含み、それはその後の分析でほとんど失われる。
 この損失を減らし、リファレンスゲノムへの依存を少なくするために、遺伝子中心のアプローチが開発されてきた。 1つの環境からの何百ものサンプルがプールされ、コンティグ内の遺伝子が予測され、95%の同一性で遺伝子カタログにまとめられられる[ref.2-5]。各サンプル中の遺伝子存在量は、リードをリファレンス遺伝子クラスターにマッピングすることによって見出される。このようにして、メタゲノムサンプルの機能的および分類学的構成、ならびにそれらの環境パラメータへの依存性を調べることができる。また、ゲノムベースの分析はビンニングによって可能になる。
 メタゲノムのショートリードの最先端のアセンブラ[ref.7 - 9]は、de-Bruijnグラフのパスとしてコンティグを見つける。このグラフは、リード内の各k-merワードと、リード内で連続して発生するk-mer間のエッジに対するノードを持っている。メタゲノムのデータでは、de-Bruijnアセンブラは、感度と選択性のトレードオフが制限されている。偽のエッジでグラフが爆発するのを避けるために、k-merは長くspecificでなければならない。しかし、集団内の多様性が高く、重複したリードがSNPによるミスマッチを含むことが多い場合、長いk-merは感度を欠く。どのk値が選択されようとも、k-merは種間で保存されているゲノム領域では十分に特異的であるには短すぎ、集団内の多様性が高い領域では十分に敏感であるには長すぎる。このジレンマはアセンブリを断片化し短くする。

 微生物集団のほとんどのSNPはコードされたタンパク質配列中のアミノ酸変化をもたらさない(補足図1)。したがって、ORFome [ref.10]およびSFA-SPA [ref.11]は、ヌクレオチドの代わりにタンパク質をアセンブリすることを提案した。しかし、それらは大規模なメタゲノムに対して実行するには遅すぎる。

 タンパク質ははるかに少なく短いリピートしか有さない。さらに、非常に類似したタンパク質配列間のキメラアセンブリ(例えば、97%以上の配列同一性)は、ゲノム中でどの遺伝子が一緒に存在するかについての誤った結論をもたらさないため、DNAベースのアセンブリでの問題を回避できる。そのため、タンパク質レベルのアセンブリでは、DNA配列ではキメラアセンブリのリスクを抑え、アセンブリできない領域のアセンブリも可能になることでカバー率が向上する。
 Plassはこれらの利点を活用してアセンブリプロセスを大幅に簡素化する。Plassはグラフフリーのgreedy iterative assembly strategyを使い、all-versus-all のオーバーラップを線形時間で計算する。これにより、単一サーバー上で巨大なリードセットをオーバーラップベースでアセンブリすることが可能になる。最も重要なことは、Plassはショートk-merのみに頼るのではなくフルアラインメントのオーバーラップを計算することによって、de-Bruijnアセンブラの特異性と感度のトレードオフの限界を克服し、複雑なメタゲノムサンプルにおいて数倍以上のタンパク質配列を回収する。
 Plassはランダムディスクアクセスを回避するためにタンパク質配列をメインメモリに格納する必要がある。そのため、入力のリードから変換されたアミノ酸ごとに1バイトのメインメモリが必要で、2〜3億の2×150 bpリードのアセンブリには〜500 GBのRAMが必要になる。これと比較すると、オーバーラップグラフのアセンブラのメモリ要件とランタイムは、処理データと超直線的に比例する。したがって、Plassは、オーバーラップグラフアセンブラの高いspecificity- sensitivityと、de-Bruijnグラフアセンブラの線形のランタイムスケーリング性能および線形のメモリスケーリング性能を組み合わせたものである。(以下略)

 Plass(Protein-Level ASSembler)は、ショートシークエンスリードを6フレームで翻訳し、タンパク質レベルでアセンブリするソフトウェアである。 Plassの主な目的は、複雑なメタゲノムデータセットアセンブリである。 Plassは土壌メタゲノムにてMegahitより10倍多いタンパク質残基をアセンブリする。 Plassは、C ++で実装され、LinuxおよびmacOSで利用可能なGPLライセンスオープンソースソフトウェアである。 ソフトウェアは複数のコア上で効率的に動作するように設計されている。 著者らは2つの冗長性フィルタリングされたリファレンスタンパク質カタログ、 640の土壌試料(SRC)からの20億の配列、および775の海洋真核生物メタトランスクリプトーム(MERC)からの292ミリオンの配列をアセンブリし、最大のタンパク質配列の無料コレクションを構築した。

 


インストール

ubuntu18.04でテストした。

Github

#bioconda 
conda install -c biocore plass

#linuxバイナリ
wget https://mmseqs.com/plass/plass-static_sse41.tar.gz; tar xvfz plass-static_sse41.tar.gz; export PATH=$(pwd)/plass/bin/:$PATH

> plass assemble

$ plass assemble

plass assemble:

Extends sequence to the left and right using ungapped alignments.

 

Please cite:

Steinegger, M. Mirdita, M., & Soding, J. Protein-level assembly increases protein sequence recovery from metagenomic samples manyfold. biorxiv, https://doi.org/10.1101/386110 (2018)

© Martin Steinegger <martin.steinegger@mpibpc.mpg.de> 

 

Usage: <i:fast(a|q)File[.gz]> | <i:fastqFile1_1[.gz] ... <i:fastqFileN_1[.gz] <i:fastqFile1_2[.gz] ... <i:fastqFileN_2[.gz]> <o:fastaFile> <tmpDir> [options]

 

align options               default   description [value range]

  -e                        0.000     Extend sequences if the E-value is below [0.0, inf]         

  --min-seq-id              0.900     Overlap sequence identity threshold [0.0, 1.0]              

 

profile options             default   description [value range]

  --num-iterations          12        Number of assembly iterations [1, inf]                      

 

misc options                default   description [value range]

  --min-length              45        minimum codon number in open reading frames                 

  --max-length              2147483647 maximum codon number in open reading frames                 

  --protein-filter-threshold 0.200     filter proteins lower than threshold [0.0,1.0]              

  --filter-proteins         1         filter proteins by a neural network [0,1]                   

 

common options              default   description [value range]

  --threads                 56        number of cores used for the computation (uses all cores by default)

  -v                        3         verbosity level: 0=nothing, 1: +errors, 2: +warnings, 3: +info

 

An extended list of options can be obtained by calling 'plass assemble -h'.

3 Database paths are required

> plass nuclassemble

$ plass nuclassemble

plass nuclassemble:

Extends sequence to the left and right using ungapped alignments.

 

Please cite:

Steinegger, M. Mirdita, M., & Soding, J. Protein-level assembly increases protein sequence recovery from metagenomic samples manyfold. biorxiv, https://doi.org/10.1101/386110 (2018)

© Martin Steinegger <martin.steinegger@mpibpc.mpg.de> 

 

Usage: <i:fast(a|q)File[.gz]> | <i:fastqFile1_1[.gz] ... <i:fastqFileN_1[.gz] <i:fastqFile1_2[.gz] ... <i:fastqFileN_2[.gz]> <o:fastaFile> <tmpDir> [options]

 

align options              default   description [value range]

  -e                       0.000     Extend sequences if the E-value is below [0.0, inf]         

  --min-seq-id             0.900     Overlap sequence identity threshold [0.0, 1.0]              

 

profile options            default   description [value range]

  --num-iterations         12        Number of assembly iterations [1, inf]                      

 

common options             default   description [value range]

  --threads                56        number of cores used for the computation (uses all cores by default)

  -v                       3         verbosity level: 0=nothing, 1: +errors, 2: +warnings, 3: +info

 

An extended list of options can be obtained by calling 'plass nuclassemble -h'.

2 Database paths are required

dockerイメージも用意されている。

 

テストラン

git clone https://github.com/soedinglab/plass.git
cd plass/examples/

 ペアエンドfastq。ラン時は作業ディレクトリを指定する。

plass assemble reads_1.fastq.gz reads_2.fastq.gz assembly.fa tmp

 シングルエンドfastq

 plass assemble reads_1.fastq.gz assembly.fa tmp

 

感想

HIseqやNovaseqクラスの複数のメタゲノムデータセットアセンブリには、準スパコンクラスのマシンとスケーリング性能が極めて高いde novoアセンブラの両方が必要なのかと考えてましたが(Extreme Scale De Novo Metagenome Assemblyの"C. Performance Results"を参照)、Preprintを読むと、PLASSは、メモリ効率の良いmegahitがsegmentation errorを起こすようなデータセットでも(プロテインの)アセンブリが可能なようです。アセンブリされたプロテインのトータルサイズも、プロテインの結果はMEGAHITの結果よりずっと長く、特に複雑性の高いデータセットで差が顕著になっています。プロテインデータベースを構築してスクリーニングする用途にも使えそうです。

引用

Protein-level assembly increases protein sequence recovery from
metagenomic samples manyfold
Martin Steinegger, Milot Mirdita, Johannes Söding

bioRxiv preprint first posted online Aug. 7, 2018

 

関連


 

 

viral quasispeciesのアセンブリを行う SAVAGE

2019 5/20 誤字修正、コメント、ヘルプ追加

 

 HIV、Zika、Ebolaなどのウイルスは、一般的にウイルス準種(viral quasispecies, wiki)と呼ばれる、遺伝的に関連しているが異なる変異株の集団として宿主に存在する。それぞれ独自のハプロタイプ配列によって特徴付けられているこれらの株は、高い突然変異率および組換え率の傾向がある(Duffy et al、2008; Domingo et al、2012)。感染サンプルに存在するウイルス準種の遺伝的多様性を捉えることを目的とした次世代シーケンシング(NGS)に基づく方法が、治療選択肢や他の臨床的に関連性のある決定を選択する際に臨床医を支援する。
 理想的には、ウイルス準種のアセンブリは、それらの存在率と共にすべてのウイルスハプロタイプを提示することによって感染の遺伝的多様性を特徴付ける。これには2つの大きな課題がある。
1、異なる系統の数は通常不明である。さらに、2つの異なる株はごくわずかな量の区別できる突然変異によって異なり得る。最後に述べるが決して軽んずべきでないが、豊富さの割合はシークエンシングエラー率と同じくらい低くなる可能性があり、これは低頻度で存在する真の突然変異の検出を妨げる。
2、遺伝的多様性と高突然変異率の両方、高品質のコンセンサスゲノム配列を表すリファレンスゲノムは、疾患発生時には時代遅れになる可能性がある。適切なゲノムの欠如は、多くのウイルス準種アセンブリにとって大きな障害となる。
 既存のすべてのアセンブリ方法では、最初のポイントまたは2番目のポイントのどちらにも対処できないことを理解することが重要である。ウイルス準種の構築に特化した最近のリファレンスガイドアプローチは、ウイルス準種の進化の根底にある原動力をモデル化する統計的枠組みを示唆した。以前のアプローチは主にハプロタイプの局所的再構築に焦点を当てていたが(Zagordi et al、2010、2011; Quince et al、2011; Huang et al、2012)、例えばハプロタイプのグローバル再構築を目的としたより高度なアプローチとして、Dirichletプロセス混合モデル(Prabhakaran et al、2014)、隠れマルコフモデル(Töpferet al、2013)、またはサンプリングスキーム(Prosperi and Salemi 2012)がある。オーバーラップグラフでパスを計算する(Astrovskaya et al、2011)、オーバーラップグラフで最大クリークを列挙する、またはコンフリクトグラフで最大独立集合を計算する(Mangul et al、2014)最近の組み合わせアプローチもある。 これらのアプローチはしっかりとポイント1を扱っているが、それらの方法のバックボーンとしてそれらの大多数は高品質のリファレンスシーケンスに依存している。これが今度はポイント2を扱えない理由である。変異パターンを逸脱すると、これらのアプローチは十分にうまく機能しない。

 一方、de novoアセンブリアプローチはリファレンスゲノムに依存しない。哺乳動物ゲノム構築のための多数のデノボアプローチが存在するが(例えば、比較評価についてはSalzberg et al、2011; Bradnam et al、2013; Gurevich et al、2013を参照)、これらの一般的方法はウイルス準種構築にはあまり適していない。主な違いは、ウイルスの突然変異率が真核生物よりも数桁高いため、1つのリードに複数のpolymorphicな部位が生じることである(Duffy et al、2008; Domingo et al、2012)。これにより、突然変異を別々のハプロタイプにphasingすることが可能になる。しかし、一般的なアセンブリ手法はこの性質を利用しない。一般的なアセンブラは、単一のコンセンサス配列を再構築することを目的としており、ひどい倍数体のゲノムを処理するようには設計されていない。これに関して、ウイルスゲノム構築を専門とする新規のアセンブラが既にあることに注意されたい(Yang et al、2012; Hunt et al、2015)。しかしながら、これらの特定のアプローチはまた、株特異的な配列よりもむしろコンセンサスゲノムをアセンブリすることを目的としており、その目的は個々の配列よりもむしろ新しいリファレンスを構築することである。著者らの知る限りでは、haplotype-resolvedのウイルス準種アセンブリの既存の唯一のアプローチはMLEHaploである(Malhotra et al、2016b)。結果として、ポイント2に対処していると、大部分の既存のデノボアセンブリ方法は満足できる程度までポイント1に対処することができない。
 考えられる主な問題は、上記の特殊なde novoウイルス準種アプローチを含むほぼすべてのNGSベースのゲノムアセンブラが、アセンブリパラダイムとしてde Bruijnグラフに依存していることである。それによって、リードはk-merに分解される。ここで、kは通常リード長よりかなり小さい。この概念を一般化したものとして、paired de Bruijn graphが導入された(Medvedev et al、2011)、これは後処理段階でメイトペアを分析する代わりにメイトペア情報をグラフ構造自体に組み入れ、結果としてより大きなcontigsを作り出す。上記のように、ウイルス種間のアセンブリでは、低頻度の変異とシーケンシングエラーを区別することが不可欠である。低頻度の変異は遺伝的に関連しており、したがって異なるリード内でco-occurするが、シーケンシングエラーは co-occurrenceのパターンを示さない。 Co-occurrenceパターンの検出は、リードを全長で調べることによって決定的にサポートされるが、この情報はde Bruijnグラフでは利用できない。一方、オーバーラップグラフは、完全長のリードを利用し、それらを小さい部分に分解しない。したがって、オーバーラップグラフのパラダイムはウイルス準種アセンブリ問題により適していると考えられる。
 オーバーラップグラフに基づくウイルス準種アセンブリのための唯一の既存の方法はHaploCliqueである(Töpferet al、2014)。このメソッドはリファレンスガイドだが、オーバーラップグラフを作成するためのアンカーポイントを提供するためだけにリファレンスを使用する。他の多くのアプローチとは異なり(Zagordi et al、2010、2011;Töpferet al、2013; Di Giallonardo et al、2014)、ハプロタイプ配列はリファレンスからではなくリードからアセンブリされる。インスピレーションを提供する一方で、HaploCliqueアルゴリズムは、比較的低いカバレッジ(1000倍以上)のデータセットにも既に過剰な計算リソースを必要とすることが証明されている。その理由は、最大クリークの列挙に基づいているためである。これは、ランタイムとメモリ空間の両方ともシーケンシングカバレッジの対数になる。したがって、アセンブリアルゴリズムのクリーク列挙のための新規でより効率的なアルゴリズムを提示する。

 リファレンスゲノム欠如問題を解決する2つの出口戦略がある。最初の戦略は、利用可能なde novoコンセンサスゲノムアセンブラのうちの1つ(その中で最も人気のあるツールはVICUNA [Yang et al、2012])を使用して、患者サンプル自体からコンセンサスゲノム配列を構築することである。このアドホックコンセンサスをリファレンスとして使用し、リファレンスガイドアプローチの1つに適用する。この戦略はMangulら (2014)によっても示唆されていて、本著者らはそれをここでさらに探求する。 2番目の戦略は、患者サンプルリードから直接オーバーラップグラフを作成することである。続いて、倍数性を意識したアセンブリアルゴリズムを使用してオーバーラップグラフから株特異的配列を抽出する。問題は、オーバーラップグラフを作成するにはすべてのリードのペアワイズ比較を行う必要があることである。このため、カバレッジのディープなデータセットで実行可能にするためには、洗練されたインデックス作成手法が必要になる。ここでは、リファレンスゲノム無しでオーバーラップグラフを作成するため、FMインデックスベースの手法を効率的に利用する方法を示す(Välimäkiet al、2012)。そのようにして、オーバーラップグラフに基づいたウイルス準種のde novoアセンブリの最初のアプローチを提供する。
要約すると、以下に関連する貢献をする。
1、ディープカバレッジシーケンシングデータセットからのオーバーラップグラフ作成
2、オーバーラップグラフアセンブリパラダイムを用いたウイルス準種アセンブリ
これらを組み合わせて、本著者らはSAVAGE (strain aware viral genome assembly)を提示する。これは、ディープカバレッジ(20,000倍以上)のシーケンスデータセットからウイルス準種のリファレンスなしのアセンブリを可能にする方法である。ここでは、オーバーラップグラフに基づく最初の本物のde novoウイルス準種アセンブリアプローチを提供するだけでなく、例えばVICUNA(Yang et al、2012)によって計算されるように、患者サンプルから生成されたアドホックコンセンサス配列を利用した最初の方法も提供する。

 最後に、ジカウイルスとC型肝炎ウイルスに感染した患者からの2つのウイルスサンプルについてSAVAGEを使ったアセンブリ結果を示す。(以下略)

 

 

f:id:kazumaxneo:20190519211711p:plain

Figure 1. An overview of the workflow and algorithms of SAVAGE.   論文より転載。

 

 インストール

依存

Python2

  • rust-overlaps
  • ncbi-blast
  • bwa mem
  • Kallisto

本体 Bitbucket 

#bioconda
conda install -c bioconda -y savage

#またはconda仮想環境を作って導入する
conda create -n savage -c bioconda -y savage python=2.7.15
conda activate savage

> savage -h

$ savage -h

usage: savage.py [-h] [-s INPUT_S] [-p1 INPUT_P1] [-p2 INPUT_P2]

                 [-m MIN_OVERLAP_LEN] [-t THREADS] --split SPLIT_NUM

                 [--revcomp] [-o OUTDIR] [--ref REFERENCE] [--no_stage_a]

                 [--no_stage_b] [--no_stage_c] [--no_overlaps]

                 [--no_preprocessing] [--no_assembly] [--count_strains]

                 [--ignore_subreads] [--merge_contigs MERGE_CONTIGS]

                 [--min_clique_size MIN_CLIQUE_SIZE]

                 [--overlap_len_stage_c OVERLAP_STAGE_C]

                 [--contig_len_stage_c CONTIG_LEN_STAGE_C] [--keep_branches]

                 [--sfo_mm SFO_MM] [--diploid]

                 [--diploid_contig_len DIPLOID_CONTIG_LEN]

                 [--diploid_overlap_len DIPLOID_OVERLAP_LEN]

                 [--average_read_len AVERAGE_READ_LEN] [--no_filtering]

                 [--max_tip_len MAX_TIP_LEN]

 

Program: SAVAGE - Strain Aware VirAl GEnome assembly

Version: 0.4.1

Release date: July 17, 2017

Contact: Jasmijn Baaijens - j.a.baaijens@cwi.nl

 

SAVAGE assembles individual (viral) haplotypes from NGS data. It expects as

input single- and/or paired-end Illumina sequencing reads. Please note that the

paired-end reads are expected to be in forward-forward format, as output by

PEAR.

 

Run savage -h for a complete description of required and optional arguments.

 

For the complete manual, please visit https://bitbucket.org/jbaaijens/savage

 

optional arguments:

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

 

basic arguments:

  -s INPUT_S            path to input fastq containing single-end reads

  -p1 INPUT_P1          path to input fastq containing paired-end reads (/1)

  -p2 INPUT_P2          path to input fastq containing paired-end reads (/2)

  -m MIN_OVERLAP_LEN, --min_overlap_len MIN_OVERLAP_LEN

                        minimum overlap length required between reads

  -t THREADS, --num_threads THREADS

                        allowed number of cores

  --split SPLIT_NUM     split the data set into patches s.t. 500 < coverage/split_num < 1000

  --revcomp             use this option when paired-end input reads are in forward-reverse orientation;

                        this option will take reverse complements of /2 reads (specified with -p2)

                        please see the SAVAGE manual for more information about input read orientations

  -o OUTDIR, --outdir OUTDIR

                        specify output directory

 

reference-guided mode:

  --ref REFERENCE       reference genome in fasta format

 

advanced arguments:

  --no_stage_a          skip Stage a (initial contig formation)

  --no_stage_b          skip Stage b (extending initial contigs)

  --no_stage_c          skip Stage c (merging maximized contigs into master strains)

  --no_overlaps         skip overlap computations (use existing overlaps file instead)

  --no_preprocessing    skip preprocessing procedure (i.e. creating data patches)

  --no_assembly         skip all assembly steps; only use this option when using --count_strains separate from assembly (e.g. on a denovo assembly)

  --count_strains       compute a lower bound on the number of strains in this sample; note: this requires a reference genome.

  --ignore_subreads     ignore subread info from previous stage

  --merge_contigs MERGE_CONTIGS

                        specify maximal distance between contigs for merging into master strains (stage c)

  --min_clique_size MIN_CLIQUE_SIZE

                        minimum clique size used during error correction

  --overlap_len_stage_c OVERLAP_STAGE_C

                        min_overlap_len used in stage c

  --contig_len_stage_c CONTIG_LEN_STAGE_C

                        minimum contig length required for stage c input contigs

  --keep_branches       disable merging along branches by removing them from the graph (stage b & c)

  --sfo_mm SFO_MM       input parameter -e=SFO_MM for sfo: maximal mismatch rate 1/SFO_MM

  --diploid             use this option for diploid genome assembly

  --diploid_contig_len DIPLOID_CONTIG_LEN

                        minimum contig length required for diploid step contigs

  --diploid_overlap_len DIPLOID_OVERLAP_LEN

                        min_overlap_len used in diploid assembly step

  --average_read_len AVERAGE_READ_LEN

                        average length of the input reads; will be computed from the input if not specified

  --no_filtering        disable kallisto-based filtering of contigs

  --max_tip_len MAX_TIP_LEN

                        maximum extension length for a sequence to be called a tip

(savage) kazu@kazu:~/taniguchi_datadir/20190518$ 

 

python freq_est.py -h

$ python freq_est.py -h

usage: freq_est.py [-h] -c CONTIGS [-o OUT_FILE] [-m MIN_LEN]

                   [--select_ids SELECT_IDS] [--kallisto]

                   [--kallisto_path KALLISTO_PATH] -l FRAGMENTSIZE -d STDDEV

                   [-f FORWARD] [-r REVERSE] [-s SUBREADS_FILE]

                   [-k LEN_CORRECTION]

 

Estimate relative frequencies for contigs produced by SAVAGE. The frequency

estimation procedure has two modes: 1. Compute frequency estimates for given

contigs by running Kallisto, then parse Kallisto output (abundance.tsv) and

compute frequency estimates from TPM counts. 2. Quick estimates based on the

relative number of reads contributing to a contig during SAVAGE assembly.

Option 1 obtains most accurate results, but requires Kallisto to be installed.

Note: For option 1, please make sure that the input fasta contains all

assembled sequences; contigs of sufficient length will be selected after

estimating TPM counts (Kallisto). For option 2, please make sure to run SAVAGE

with the option --use_subreads for better frequency estimates of stage b/c

contigs.

 

optional arguments:

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

 

general arguments:

  -c CONTIGS, --contigs CONTIGS

                        fastq or fasta file containing contigs

  -o OUT_FILE, --out OUT_FILE

                        write output to file; if not specified, output is

                        written to stdout.

  -m MIN_LEN, --min_len MIN_LEN

                        only consider contigs at least this length

  --select_ids SELECT_IDS

                        select target contigs for frequency estimation by

                        providing a comma-separated list of contig identifiers

 

kallisto-mode arguments:

  --kallisto            use Kallisto for improved frequency estimation

  --kallisto_path KALLISTO_PATH

                        path to kallisto executable; needed only if Kallisto

                        path is not included in PATH

  -l FRAGMENTSIZE, --fragmentsize FRAGMENTSIZE

                        Estimated average fragment size

  -d STDDEV, --stddev STDDEV

                        Estimated standard deviation of fragment size

  -f FORWARD, --forward FORWARD

                        original forward input reads (before assembly) in

                        fastq format

  -r REVERSE, --reverse REVERSE

                        original reverse input reads (before assembly) in

                        fastq format

 

quick-mode arguments:

  -s SUBREADS_FILE, --subreads SUBREADS_FILE

                        subreads file produced by savage, corresponding to the

                        fastq file provided

  -k LEN_CORRECTION, --correction LEN_CORRECTION

                        (optional) correction term for computing effective

                        contig length: use c = fragment_size - 2*(read_len-

                        min_overlap_len)

(savage) kazu@kazu:/media/kazu/1/virus/savage$

 

 freq_est.pyはcondaによるインストールでパスは通っていない。 

 

テストラン

git clone https://bitbucket.org/jbaaijens/savage.git
cd savage/example/

リードのstatistics

> seqkit stats input_fas/*

file                     format  type  num_seqs  sum_len  min_len  avg_len  max_len

input_fas/paired1.fastq  FASTQ   DNA        200   49,970      245    249.9      250

input_fas/paired2.fastq  FASTQ   DNA        200   49,962      246    249.8      250

input_fas/singles.fastq  FASTQ   DNA      2,000  896,174      380    448.1      490

 

 

1、merge

PEARでペアエンドをマージする。テストデータではすでに用意されているため実行する必要はない。

pear -f pair1.fq -r pair2.fq -o merged

merged.fastqとunmpergedのfastqが出力される。

 

2-1、De novo

ペアエンドfastqとマージしたfastqを指定してDe novoアセンブリを実行。ペアエンドfastはPEARでmegeし、マージされなかったペアエンドfastq(forward- forwad)を指定する。error correctionが500-1000xのカバレッジに最適化されているため、それ以上ディープなデータの場合、データを 500-1000xカバレッジになるよう分割する"--split"フラグを立てる。

savage -s input_fas/singles.fastq -p1 input_fas/paired1.fastq -p2 input_fas/paired2.fastq -m 200 --split 1 -o output -t 12
  • -s    path to input fastq containing single-end reads  
  • -p1     path to input fastq containing paired-end reads (/1)  
  • -p2     path to input fastq containing paired-end reads (/2)  
  • -o       specify output directory
  • -m      minimum overlap length required between reads  
  • -t     allowed number of cores  
  • --split     split the data set into patches s.t. 500 < coverage/split_num < 1000
  • --merge_contigs    specify maximal distance between contigs for merging into master strains (stage c). By default this is set to 0.

出力

ls -lth out/

$ ls -lth out/

合計 112K

-rw-rw-r-- 1 kazu  19K  5月 19 21:57 contigs_stage_c.fasta

drwxrwxr-x 2 kazu 4.0K  5月 19 21:57 stage_c/

-rw-rw-r-- 1 kazu  19K  5月 19 21:57 contigs_stage_b.fasta

drwxrwxr-x 2 kazu 4.0K  5月 19 21:57 stage_b/

-rw-rw-r-- 1 kazu  60K  5月 19 21:57 contigs_stage_a.fasta

drwxrwxr-x 4 kazu 4.0K  5月 19 21:57 stage_a/

stage aがerror correction、stage bがオーバーラップグラフに基づくアセンブリ、stage cがマージしてhaplotypeを再構成するステージになる。stage cのマージは、"--merge_contigs"のパラメータでマージされるかどうか決定される。 例えばdefaultの0から0.01(=1%) まで増やせば、1%以下の違いしかない配列は1つのhaplotypeにマージされる。

 

2-2、Reference guide 

またはリファレンスを指定してReference guidedアセンブリ実行。

savage --ref /path/to/refeence.fasta --split patch_num --min_overlap_len M --s singles.fastq --p1 paired1.fastq --p2 paired2.fastq -t 12 -o output

リファレンスはfull pathで指定しないとエラーになる。

 

3、Counting virus strains

python freq_est.py --contigs output/contigs_stage_c.fasta --subreads output/stage_c/subreads.txt --min_len 1000 \
-l <Estimated average fragment size> \
-d <Estimated standard deviation of fragment size>

 

 

引用
De novo assembly of viral quasispecies using overlap graphs.

Baaijens JA, Aabidine AZE, Rivals E, Schönhuth A

Genome Res. 2017 May;27(5):835-848

 

関連


メタゲノムアセンブリを評価する MetaQUAST

 

メタゲノミクスは、環境サンプルから直接採取した遺伝物質を研究する。 NGS技術は、クローニングなしに少量の生物からDNAを抽出しショートリードシーケンシングすることを可能にする。しかし、そのような実験で作成されたデータは膨大でノイズが多く、そして豊富さと相同性が大きく異なる何千もの種からのDNA断片を含む傾向がある。これらの課題は、メタゲノムアセンブリの新たな計算上の問題をもたらし、それに続いて標準的なベンチマーク手順を必要とする(Boisvert et al、2012; Peng et al、2012; Haider et al、2014)。

 既存のほとんどのアセンブリ評価手法は、メタゲノムと連携するようには設計されていない。しかしながら、アセンブリに関してリードの尤度をカウントする方法(Clark e t al、2013; Ghodsi et al、2013)、または保存された普遍的なシングルコピー遺伝子の含量を決定する方法(Parks et al、2015; Simao et al、2005)が存在する。残念ながら、どれもclosely relatedなリファレンスゲノムへのコンティグアラインメントを使用していない。この論文では、QUASTのメタゲノム改良版MetaQUASTを紹介する(Gurevich et al、2013)。 QUASTは、与えられたclosely relatedなリファレンスゲノムとのアラインメントに基づいてエラーを検出し、さらにN50のようなコンティグ統計およびユーザー提供のリファレンス配列がなくても構成種の概観を与える遺伝子量をプロットする。メタゲノムアセンブリに対処するために、MetaQUASTはいくつか新しい機能を追加している:(i)無制限数のリファレンスゲノムを使用する能力、(ii)自動化された種含有量検出、(iii)キメラコンティグの検出。

 既知の種含有量を持つよく研究されたメタゲノムデータセット(Qin et al、2010)またはmockリードセット(Boisvert et al、2012; Namiki et al、2012)がある。それらは、リファレンスゲノムへのアラインメントに基づいてアセンブリ方法を評価するためにMetaQUASTと共に使用することができる。複数リファレンスパイプラインは、4つの主要なステップで構成されている(論文補足図S1)。

1、すべてのリファレンスゲノムは1つのファイルに結合される(結合リファレンス)。 QUASTは、結合リファレンスversus 入力アセンブリで評価される。 1つだけではなくすべてのあいまいなアラインメントも報告するように強制する。closely relatedな種を含むメタゲノムデータセットのために、すべてのあいまいなアラインメントは不可欠である。
2、生成されたアラインメントに基づいて、全てのコンティグをグループに分割する。複数ゲノムにマッピングされたコンティグは、すべてのマッチンググループに入る。アラインメントされていないコンティグは1つのextraグループに入る。
3、次に、各入力リファレンスに対してQUASTが別々に起動する。未アラインメントグループは、入力リファレンスなしで処理される。
4、最後に、すべてのQUAST実行の結果が要約レポートと視覚化してまとめられる。ユーザーは、各分析の詳細なフルQUAST出力と、データセット全体の結果の概観の両方を見ることができる。
スタンダードな品質評価metrics(N50、genome fractionなど)に加え、2つの測定基準を追加した。

 I nterspecies translocationsの数:隣接配列が異なるリファレンスにアラインメントする。

  誤ってアセンブリされた可能性があるコンティグの数:アラインメントサイズによらず未知のゲノムとの種間転座を含むコンティグ数。
GeneMarkSを使用する通常のQUASTとは対照的に、MetaQUASTはメタゲノム用に開発されたMetaGeneMarkを遺伝子予測に使用する(Zhu et al、2010)。

デノボ評価

ほとんどの実験的メタゲノミクス研究は、リファレンス情報が利用できないde novoアセンブリで動作する。 MetaQUASTが入力リファレンス配列または種リストなしで実行されると、それは種の内容を識別しそして自動的にリファレンス配列を引き出すことを試みる。このアルゴリズムは、研究者が主に微生物群集に興味を持っているという仮定の下で機能するので、検索は細菌と古細菌に限定されることに注意する。

 ワークフロー(論文補足図S2参照)は、BLASTn(Camacho et  al、2009)を適用して、SILVAデータベースから16S rRNA配列にコンティグをアラインメントさせることから開始する(Quast et al、2012)。ほとんど全ての微生物種に存在する16Sサブユニットは高度に保存された配列であるが、生物を分類群に分類するにも役立つ超可変領域も含む。検出された各種について、最高のスコアを持つ1つの株がアセンブリに残る。

 上位50の生物がNCBIに対して照会され、それぞれの種について最も断片化されていない配列がダウンロードされる。生物間のrRNAオペロンのコピー数の違い、および16S遺伝子のゲノム内の不均一性に関する既知の問題のために、ダウンロードされたゲノム配列のいくつかは評価中のアセンブリで表現できない可能性がある。 MetaQUASTは、10%未満のコンティグカバレッジラクションですべてのアセンブリに対してゲノムを除去することによって、偽陽性をフィルタリングしようとする。全ての配列が非常に低いゲノムfractionを有する特別な場合には、リストはフィルタリングされない。

 結果として、本発明者らは、おそらくアセンブリ配列によって表される一組のゲノムを得る。これらのシーケンスを使って(セクション2.1のように)MetaQUASTを起動し、通常のリファレンスベースの解析の場合と同じ出力ファイルを生成する。

私たち(著者ら)のアプローチは正確さと時間/メモリ消費の間の妥協である。より正確な結果を得るためには、MGTAXA(Williamson et al、2012)、または正確なリードアライメントに基づく方法、例えばKraken(Wood and Salzberg、2014)またはCLARK(Ounit et al、2015)を使う。非常に正確な結果を得るには、全NCBI-nrデータベースに対するBLASTx(Altschul et al、1990)検索によって得ることができる。得られた種名リストは、プレーンテキスト形式でMetaQUASTに送ることができ、NCBIデータベースから指定された配列をダウンロードしてそれらをリファンレンスに基づく評価に使用することができる。(以下略)

 

 

manual

http://quast.bioinf.spbau.ru/manual.html

 

インストール

依存

QUAST can be run on Linux (64-bit and 32-bit with slightly limited functionality) or macOS (OS X).

  • Python2 (2.5 or higher) or Python3 (3.3 or higher)
  • GCC 4.7 or higher
  • Perl 5.6.0 or higher
  • GNU make and ar
  • zlib development files

In addition, QUAST submodules require:

  • Time::HiRes perl module for GeneMark-ES (gene prediction in eukaryotes)
  • Java 1.8 or later for GRIDSS (SV detection)
  • R for GRIDSS (SV detection)

Github

#bioconda
#仮想環境を作って導入する。ここでは5.0.2を選択
conda create -c bioconda -n quast python=2.7.14 quast=5.0.2
source activate quast

metaquast.py -h 

$ metaquast.py -h

MetaQUAST: Quality Assessment Tool for Metagenome Assemblies

Version: 5.0.0

 

Usage: python /Users/kazuma/.pyenv/versions/miniconda3-4.3.14/bin/metaquast.py [options] <files_with_contigs>

 

Options:

-o  --output-dir  <dirname>       Directory to store all result files [default: quast_results/results_<datetime>]

-r   <filename,filename,...>      Comma-separated list of reference genomes or directory with reference genomes

--references-list <filename>      Text file with list of reference genome names for downloading from NCBI

-g  --features [type:]<filename>  File with genomic feature coordinates in the references (GFF, BED, NCBI or TXT)

                                  Optional 'type' can be specified for extracting only a specific feature type from GFF

-m  --min-contig  <int>           Lower threshold for contig length [default: 500]

-t  --threads     <int>           Maximum number of threads [default: 25% of CPUs]

 

Advanced options:

-s  --split-scaffolds                 Split assemblies by continuous fragments of N's and add such "contigs" to the comparison

-l  --labels "label, label, ..."      Names of assemblies to use in reports, comma-separated. If contain spaces, use quotes

-L                                    Take assembly names from their parent directory names

-e  --eukaryote                       Genome is eukaryotic (primarily affects gene prediction)

    --fungus                          Genome is fungal (primarily affects gene prediction)

    --large                           Use optimal parameters for evaluation of large genomes

                                      In particular, imposes '-e -m 3000 -i 500 -x 7000' (can be overridden manually)

-k  --k-mer-stats                     Compute k-mer-based quality metrics (recommended for large genomes)

                                      This may significantly increase memory and time consumption on large genomes

    --k-mer-size                      Size of k used in --k-mer-stats [default: 101]

    --circos                          Draw Circos plot

-f  --gene-finding                    Predict genes using MetaGeneMark

    --glimmer                         Use GlimmerHMM for gene prediction (instead of the default finder, see above)

    --gene-thresholds <int,int,...>   Comma-separated list of threshold lengths of genes to search with Gene Finding module

                                      [default: 0,300,1500,3000]

    --rna-finding                     Predict ribosomal RNA genes using Barrnap

-b  --conserved-genes-finding         Count conserved orthologs using BUSCO (only on Linux)

    --operons  <filename>             File with operon coordinates in the reference (GFF, BED, NCBI or TXT)

    --max-ref-number <int>            Maximum number of references (per each assembly) to download after looking in SILVA database.

                                      Set 0 for not looking in SILVA at all [default: 50]

    --blast-db <filename>             Custom BLAST database (.nsq file). By default, MetaQUAST searches references in SILVA database

    --use-input-ref-order             Use provided order of references in MetaQUAST summary plots (default order: by the best average value)

    --contig-thresholds <int,int,...> Comma-separated list of contig length thresholds [default: 0,1000,5000,10000,25000,50000]

-u  --use-all-alignments              Compute genome fraction, # genes, # operons in QUAST v1.* style.

                                      By default, QUAST filters Minimap's alignments to keep only best ones

-i  --min-alignment <int>             The minimum alignment length [default: 65]

    --min-identity <float>            The minimum alignment identity (80.0, 100.0) [default: 95.0]

-a  --ambiguity-usage <none|one|all>  Use none, one, or all alignments of a contig when all of them

                                      are almost equally good (see --ambiguity-score) [default: one]

    --ambiguity-score <float>         Score S for defining equally good alignments of a single contig. All alignments are sorted 

                                      by decreasing LEN * IDY% value. All alignments with LEN * IDY% < S * best(LEN * IDY%) are 

                                      discarded. S should be between 0.8 and 1.0 [default: 0.99]

    --unique-mapping                  Disable --ambiguity-usage=all for the combined reference run,

                                      i.e. use user-specified or default ('one') value of --ambiguity-usage

    --strict-NA                       Break contigs in any misassembly event when compute NAx and NGAx.

                                      By default, QUAST breaks contigs only by extensive misassemblies (not local ones)

-x  --extensive-mis-size  <int>       Lower threshold for extensive misassembly size. All relocations with inconsistency

                                      less than extensive-mis-size are counted as local misassemblies [default: 1000]

    --scaffold-gap-max-size  <int>    Max allowed scaffold gap length difference. All relocations with inconsistency

                                      less than scaffold-gap-size are counted as scaffold gap misassemblies [default: 10000]

    --unaligned-part-size  <int>      Lower threshold for detecting partially unaligned contigs. Such contig should have

                                      at least one unaligned fragment >= the threshold [default: 500]

    --skip-unaligned-mis-contigs      Do not distinguish contigs with >= 50% unaligned bases as a separate group

                                      By default, QUAST does not count misassemblies in them

    --fragmented                      Reference genome may be fragmented into small pieces (e.g. scaffolded reference) 

    --fragmented-max-indent  <int>    Mark translocation as fake if both alignments are located no further than N bases 

                                      from the ends of the reference fragments [default: 85]

                                      Requires --fragmented option

    --upper-bound-assembly            Simulate upper bound assembly based on the reference genome and reads

    --upper-bound-min-con  <int>      Minimal number of 'connecting reads' needed for joining upper bound contigs into a scaffold

                                      [default: 2 for mate-pairs and 1 for long reads]

    --est-insert-size  <int>          Use provided insert size in upper bound assembly simulation [default: auto detect from reads or 255]

    --plots-format  <str>             Save plots in specified format [default: pdf].

                                      Supported formats: emf, eps, pdf, png, ps, raw, rgba, svg, svgz

    --memory-efficient                Run everything using one thread, separately per each assembly.

                                      This may significantly reduce memory consumption on large genomes

    --space-efficient                 Create only reports and plots files. Aux files including .stdout, .stderr, .coords will not be created.

                                      This may significantly reduce space consumption on large genomes. Icarus viewers also will not be built

-1  --pe1     <filename>              File with forward paired-end reads (in FASTQ format, may be gzipped)

-2  --pe2     <filename>              File with reverse paired-end reads (in FASTQ format, may be gzipped)

    --pe12    <filename>              File with interlaced forward and reverse paired-end reads. (in FASTQ format, may be gzipped)

    --mp1     <filename>              File with forward mate-pair reads (in FASTQ format, may be gzipped)

    --mp2     <filename>              File with reverse mate-pair reads (in FASTQ format, may be gzipped)

    --mp12    <filename>              File with interlaced forward and reverse mate-pair reads (in FASTQ format, may be gzipped)

    --single  <filename>              File with unpaired short reads (in FASTQ format, may be gzipped)

    --pacbio     <filename>           File with PacBio reads (in FASTQ format, may be gzipped)

    --nanopore   <filename>           File with Oxford Nanopore reads (in FASTQ format, may be gzipped)

    --ref-sam <filename>              SAM alignment file obtained by aligning reads to reference genome file

    --ref-bam <filename>              BAM alignment file obtained by aligning reads to reference genome file

    --sam     <filename,filename,...> Comma-separated list of SAM alignment files obtained by aligning reads to assemblies

                                      (use the same order as for files with contigs)

    --bam     <filename,filename,...> Comma-separated list of BAM alignment files obtained by aligning reads to assemblies

                                      (use the same order as for files with contigs)

                                      Reads (or SAM/BAM file) are used for structural variation detection and

                                      coverage histogram building in Icarus

    --sv-bedpe  <filename>            File with structural variations (in BEDPE format)

 

Speedup options:

    --no-check                        Do not check and correct input fasta files. Use at your own risk (see manual)

    --no-plots                        Do not draw plots

    --no-html                         Do not build html reports and Icarus viewers

    --no-icarus                       Do not build Icarus viewers

    --no-snps                         Do not report SNPs (may significantly reduce memory consumption on large genomes)

    --no-gc                           Do not compute GC% and GC-distribution

    --no-sv                           Do not run structural variation detection (make sense only if reads are specified)

    --no-gzip                         Do not compress large output files

    --no-read-stats                   Do not align reads to assemblies

                                      Reads will be aligned to reference and used for coverage analysis,

                                      upper bound assembly simulation, and structural variation detection.

                                      Use this option if you do not need read statistics for assemblies.

    --fast                            A combination of all speedup options except --no-check

 

Other:

    --silent                          Do not print detailed information about each step to stdout (log file is not affected)

    --test                            Run MetaQUAST on the data from the test_data folder, output to quast_test_output

    --test-no-ref                     Run MetaQUAST without references on the data from the test_data folder, output to quast_test_output.

                                      MetaQUAST will download SILVA 16S rRNA database (~170 Mb) for searching reference genomes

                                      Internet connection is required

-h  --help                            Print full usage message

-v  --version                         Print version

 

2019年5月現在、最新版はv5.0.2。ロングリードを使うオプションも存在する。

 

データベースの準備

#インストール後 #To be able to use those, please run 
quast-download-gridss
quast-download-silva
quast-download-busco > quast

 

実行方法

アセンブリした配列とリファレンスを指定する。複数ある場合、リファレンスはカンマで繋ぐ。

metaquast.py contigs_1.fasta contigs_2.fasta contigs_3.fasta \
-r reference_1,reference_2,reference_3 \
-l asembly1,assembly2,assembly3 \
 -1 reads1.fq.gz -2 reads2.fq.gz \
-o output -t 40 --gene-finding
  • -o     Directory to store all result files [default: quast_results/results_<datetime>]
  • -r      Comma-separated list of reference genomes or directory with reference genomes
  •  -l     Names of assemblies to use in reports, comma-separated. If contain spaces, use quotes
  • --gene-finding   Predict genes using MetaGeneMark
  • --glimmer           Use GlimmerHMM for gene prediction (instead of the default finder, see above)

 

リファレンスを指定しない場合、16S rRNAの分析結果をもとにリファレンスがダウンロードされ、比較が行われる。

ペアエンドfastqも指定(マッピング済みならbamも使える)。

metaquast.py contigs_1 contigs_2 contigs_3 \
 -1 reads1.fq.gz -2 reads2.fq.gz \
-o output -t 40 --gene-finding

 

 引用
MetaQUAST: evaluation of metagenome assemblies
Mikheenko A, Saveliev V, Gurevich A

Bioinformatics. 2016 Apr 1;32(7):1088-90

 

関連


 

 

ターゲットアンプリコンシーケンシングのプライマーを除く pTrimmer

 

 ゲノムの変異検出は、臨床がん研究においてますます一般的になっている[ref.1]。多重アンプリコンに基づくディープシーケンシングは、特定の疾患関連遺伝子の突然変異検出のための主要なアプローチの1つである[ref.1、2、3]。がん関連遺伝子の変異を同定するために、多くのアルゴリズムが開発されている[ref.2、4]。マルチプレックスアンプリコンに基づく次世代シークエンシング(NGS)からのリードは通常2つの部分、すなわち遺伝子特異的プライマー(すなわち標的領域を増幅するために使用される配列)および目的の領域を含む[ref.5]。隣接するアンプリコン配列は通常、目的の遺伝子全体を網羅するように互いに重複するように設計されている[ref.6]。しかしながら、無効なオリゴヌクレオチド合成はプライマーにエラーを導入することになり、これは、トリミングされていないNGSリードをリファレンスゲノムに直接マッピングすることで信頼できない突然変異情報をもたらすことを意味する。具体的には、プライマーの合成エラーは、対応するゲノム部位における変異対立遺伝子頻度(VAF)の値を増加させ、したがって突然変異コーラーソフトウェアにより「真の」ヌクレオチド変異確率を増加させる変異情報を含むと誤って見なされる。

 Cutadapt [ref.7]やAlienTrimmer [ref.8]などのいくつかのツールが現在利用可能であるか、またはプライマートリミングと互換性がある。これらのツールは、アダプタシーケンス削除で効率的に機能する。しかし、マルチプレックスアンプリコンシークエンシングには通常、数百または数千のプライマー配列がある[ref.2]。これらのツールはマルチプレックスプライマー配列のトリミングにおいて非効率的になる。 bamベースのプライマー除去ツールであるBAMClipper(紹介)は、ユーザーが関心のある領域の端付近での挿入と欠失の検出に焦点を当てている[ref.5]。ただし、Perlの実装とソフトクリッピングベースのプライマーマッチングアルゴリズムにより、大規模なデータセットの処理がやや低くなる。Trimmomatic [ref.9]のような他のいくつかのツールは、5 '末端および3'末端の両方からのプライマー配列の代わりに3 '末端からの特定の配列のみをトリミングする。 cutPrimers[ref.10](紹介)は、マルチプルプライマー配列を削除するために特別に開発されたが、正規表現マッチングとPython実装のアルゴリズムは、感度、特異性、パフォーマンスが低いという結果になる。

 本研究では、マルチプレックスアンプリコン配列データからプライマー配列をトリミングする新しいツールpTrimmerを開発した。感度と特異度を高めるために、本著者らはk-mersアルゴリズム[ref.11]とNeedleman-Wunschアルゴリズム[ref.12]の両方を採用した。そしてC実装により高いパフォーマンスを保証する。このツールは、WindowsシステムとLinuxシステムの両方で使用できる。他の3つのツールを使用したベンチマーク分析により、pTrimmerはプライマートリミングにおいて非常に時間効率がよく、正確であることがわかった。

 

 

f:id:kazumaxneo:20190515201546p:plain

Workflow of pTrimmer. 論文より転載 

 


インストール

ubuntu16.04でテストした。

The program could run on a standard dual core laptops with 8 GB of RAM on both windows(win7 or win10) and linux(centos or ubuntu).

ビルド依存

  • zlib-1.2.7
  • the gcc compiler should be available on your server or laptop.

本体 Github

git clone https://github.com/DMU-lilab/pTrimmer.git
cd pTrimmer/
make 

> ./pTrimmer-1.3.1 -h

# ./pTrimmer-1.3.1 -h

[Err::ParseOpt::82]  Please give the [requied] parmeters!

 

Usage: pTrimmer [options]

Version: 1.3.1

 

Options:

       -h|--help        print help infomation

       -l|--keep        keep the complete reads if can't locate primer

                        sequence [default: discard the reads]

       -s|--seqtype     [required] the sequencing type [single|pair]

       -a|--ampfile     [required] input amplicon file [.txt]

       -f|--read1       [required] read1(forward) for fastq file [.fq|.gz]

       -r|--read2       [optional] read2(reverse) for paired-end seqtype [.fq|.gz]

       -o|--outdir      [required] output directory for trimed fastq file [dir]

       -q|--minqual     [optional] the minimum average quality to keep after triming [20]

       -k|--kmer        [optional] the kmer lenght for indexing [8]

       -m|--mismatch    [optional] the maxmum mismatch for primer seq [3]

 

テストラン

fastqとprimer配列を記載したファイル(例: pTrimmer/Example/data_amplicon.txt)を指定する。

cd pTrimmer/Example/
mkdir output
pTrimmer-1.3.1 -f data_R1.fq.gz -r data_R2.fq.gz \
-a data_amplicon.txt -s pair -o out

> ls -al output/

# ls -al output/

total 332

drwxr-xr-x  5 root root    170 May 13 18:11 .

drwxr-xr-x 36 root root   1224 May 13 18:11 ..

-rw-r--r--  1 root root 218742 May 13 18:11 Summary.ampcount

-rw-r--r--  1 root root  58400 May 13 18:11 data_trim_R1.fq

-rw-r--r--  1 root root  54878 May 13 18:11 data_trim_R2.fq

 

引用

pTrimmer: An efficient tool to trim primers of multiplex deep sequencing data

Xiaolong Zhang†, Yanyan Shao†, Jichao Tian†, Yuwei Liao, Peiying Li, Yu Zhang, Jun Chen, Zhiguang Li
†Contributed equally
BMC Bioinformatics 2019 20:236

 

RNA seqデータの正規化を行いアセンブリ負荷を軽減する ORNA

2019 5/17 誤字修正

 

 シーケンサスループットの増加および価格の低下に伴い、高カバレッジシーケンシングデータセットの生成は日常的になっている。これは、ゲノムおよびトランスクリプトームのデノボアセンブリのためのいくつかの異なるアプローチの開発を促進した(Miller et al、2010; Moreton et al、2016)。しかし、大きなゲノムやトランスクリプトームのアセンブリは、リソース集約的な作業である。

 データセットのサイズが大きいため、ゲノムおよびメタゲノムシーケンシングに適用される1つまたはいくつかのデータセットについて、de novoアセンブリのデータ構造をよりスペース効率的にすることに集中している(Chikhi et al、2016; Pell et al, 2012)。圧縮ゲノミクスと呼ばれる別のアプローチは、高速化計算のためにデータセットの圧縮表現を見つけることを扱い、これはリードアラインメントおよびSNPコールに首尾よく適用されている(Berger et al、2013; Loh et al、2012)。

 ここでは、de novo transcriptome assemblyのためのデータ削減アプローチが不均一なRNAシーケンスデータセットのパフォーマンスに与える影響を調べる。これは、de Bruijnグラフ(DBG)に依存する現在のアセンブリ方法が大量のメインメモリを消費するため、重要な問題である(Grabherr et al、2011; Robertson et al、2010; Schulz et al、2012)。しかし、それは情報理論の観点からも興味深い問題である。データのどの部分が実際にアセンブラによって使用されているのか?

 シーケンスエラーを取り除き、データサイズを減らすための簡単な方法は、リードサフィックスと低品質のプレフィックスを削除することである。これは一般にアセンブリ性能の低下につながり(MacManes、2014; Mbandi et al、2014)、現在のデータセットの高い冗長性には対処していない。他のアプローチは、RNA-seqデータセットにおけるシーケンシングエラーの効率的な訂正を可能にし、それはデータがエラー訂正されなければならないので実行時間は増加するが、一般にアセンブリ性能の改善につながる (Le et al., 2013; Song and Florea, 2015)。 

 冗長性を除去するための直接的なアプローチはCD-HIT(Fu et al、2012)のようなアルゴリズムを使用して配列類似性に従ってリードをクラスター化し、そしてクラスター中の非常に類似したリードを除去することである。最近アセンブリ転写物をクラスター化するの改良がなされたが(Srivastava et al、2016)、アセンブリ前に何億ものリードをクラスター化することは依然として困難である。

 特にアセンブリと組み合わせて使用​​される別のアプローチは、khmerパッケージに実装された「デジタル正規化」(Diginorm; Brown et al、2012)である(Crusoe et al、2015)。 Diginormは、最小カウントスケッチデータ構造を使用して、リードデータセットをストリーミングしてk-merの存在量を推定する。ユーザが選択した存在量のしきい値tを使用して、中央値のk-merカバレッジがtを超えると、リードが削除される。これに似たアイデアがTrinityのin silico正規化(TIS)で、これはTrinityアセンブラー・パッケージの一部である(Haas et al、2013)。データセット内のリードごとに、TISはリード内のk-merのカバレッジ中央値を計算する。カバレッジ中央値が目的のカバレッジよりも小さい場合、リードは常に保持される。そうでなければ、それは所望のカバレッジとmedianカバレッジの比に等しい確率で保たれる。さらに、リードの平均k-merカバレッジに対するk-merカバレッジのSDの比率がカットオフより高い場合、リードは削除される。最近開発されたNeatFreqアルゴリズム(McCorrison et al、2014)は、中央値のk-mer頻度に基づいてリードをビンにクラスタ化する。

 k-merカバレッジベースの正規化の利点は3 foldである。(i)冗長性の高いリードが削除され、アセンブリに必要なメモリおよびランタイム要件が軽減される。(ii)このプロセスにより誤ったリードが削除される。(iii) 正規化は高速で、アセンブラが消費するメモリのごく一部しか消費しない。アセンブリのパフォーマンスに大きな影響を与えることなく、多くの場合データの大部分を削除できることが示されているため、これによりアセンブリ問題の計算の複雑さが本質的に軽減される(Brown et al、2012; Haas et al、2013)。しかしながら、以前のアルゴリズムは、有用なk-merを含むデータの重要な部分を保存することに関していかなる確実性も与えない。少量ではあるが重要なk-merを含むリードは削除される可能性がある。これにより、DBG内の接続が失われる可能性があるため、断片化されたアセンブリが生成される可能性がある。これはRNA-seqやメタゲノミクスのように、カバレッジが不均一なシーケンシングデータセットにとって特に問題となる。

 ここで、DBGバックボーン(重み付けされていないノードおよびエッジ)を失うことなくリードを削減し、相対ノード存在量は元のDBGと比較して削減データセットで保持されるという考えに基づいてOptimized Read Normalization Algorithm (ORNA)を提案する。 各リードがm  k-mersからなるn個のリードセットが与えられると、リード正規化は、リードに関するset multi-cover(SMC)最適化問題として表現される。 この研究では、実際にうまくいくことが証明されている𝒪(nm log(nm))時間の発見的アルゴリズムが示唆されている。 正規化データとECデータを分析すると、実行時間とメモリ消費量を大幅に節約しながら、より優れたアセンブリを作成できることがわかる。 ソフトウェアはMITライセンスの下でhttps://github.com/SchulzLab/ORNAで無料で入手できる。

 

 

ORNAはDiginormのスピリットで開発されたリード正規化ソフトウェアである。 ORNAは計算コストが安く、オリジナルのデータセットからすべてのkmerを確実に保持する。ユーザーがカバレッジの高いデータセットを使用しているが、データの冗長性を排除するために、de novoアセンブリを実行するのに十分な計算能力(特にメモリだけでなく時間も限られている)を持っていない場合に使用できる。多くのシーケンスデータセットをマージするのにも使用できる。データセットに大きく依存するため、作成されるアセンブリに大きな影響を与える可能性があることに注意する必要がある。

 

インストール

ubuntu16.0.4でテストした(docker使用、ホストmacos10.12)。

依存

  • Linux operating system with gcc version >=4.7
  • All the analysis for the manuscript was performed on Debain 8 operating system
  • cmake

本体 Github

git clone https://github.com/SchulzLab/ORNA
cd ORNA/
bash install.sh
cd build/bin/

./ORNA -h

# ./ORNA -h

ERROR: Unknown parameter '-h'

 

[ORNA options]

       -cs       (1 arg) :    Collect stat  [default '0']

       -type     (1 arg) :    Type for the output file (fasta/fastq)  [default 'fasta']

       -binsize  (1 arg) :    Bin Size  [default '1000']

       -ksorting (1 arg) :    Kmer based Sorting  [default '0']

       -sorting  (1 arg) :    Quality Sorting  [default '0']

       -pair2    (1 arg) :    Second read of the pair  [default 'ORNAERROR']

       -pair1    (1 arg) :    First read of the pair  [default 'ORNAERROR']

       -base     (1 arg) :    Base for the logarithmic function  [default '1.7']

       -output   (1 arg) :    Prefix of the output File  [default 'Normalized']

       -input    (1 arg) :    Input File  [default 'ORNAERROR']

       -kmer     (1 arg) :    kmer required  [default '21']

       -nb-cores (1 arg) :    number of cores  [default '0']

       -verbose  (1 arg) :    verbosity level  [default '1']

       -version  (0 arg) :    version

       -help     (0 arg) :    help

 

 

実行方法

 ORNAは元のデータセットのすべてのkmerを保持しようとするため、誤ったkmerも保持される。 よって、前もってデータのエラー訂正しておく方がより多くの不要なリードを削減できる。 RNA-seqまたは他の不均一なデータの場合は、ORNAとうまく機能することが著者らによって確認されているSEECERアルゴリズム(github)を使用することが推奨されている。

 

fastqまたはfastaを指定する。

シングルエンドfastq

ORNA -input single.fq -output output -base 1.7 -kmer 21 -nb-cores 16 -type fastq
  • -input     Input File  [default 'ORNAERROR']
  • -base      Base for the logarithmic function  [default '1.7']
  • -kmer     kmer required  [default '21']
  • -output      Prefix of the output File  [default 'Normalized']
  • -type    Type for the output file (fasta/fastq) [default 'fasta']

 

ペアエンド fastq

ORNA -input -pair1 pair1.fq -pair2 pair2.fq -output output -base 1.7 -kmer 21 -nb-cores 16 -type fastq
  • -pair1     First read of the pair  [default 'ORNAERROR'] 
  • -pair2     Second read of the pair  [default 'ORNAERROR']

output_1.fqとoutput_2.fqが出力される。

 

引用

In silico read normalization using set multi-cover optimization
Dilip A Durai, Marcel H Schulz
Bioinformatics, Volume 34, Issue 19, 01 October 2018, Pages 3273–3280, 

 

関連

BBnorm

KATにはk-merとGCでフィルタリングする機能があります。


 

 

 

metaSPAdes

 

 メタゲノムシーケンシングは、細菌集団の分析ならびに新規な生物および遺伝子の発見のための選択技術として浮上している(Tyson et al, 2004、Venter et al, 2004、Yooseph et al, 2007、Arumugam et al, 2011)。初期のメタゲノミクス研究の1つにおいて、Venterら(2004)は複雑なSargasso Seaの微生物群集をアセンブリしようとした。しかし研究で述べられたように失敗した。メタゲノミクス研究のスペクトルの反対側で、Tysonら (2004)が少数の種から成る簡単な微生物群集のアセンブリに成功した。

 これらの画期的な研究(Tyson et al, 2004; Venter et al, 2004)は、従来のアセンブリツール、すなわちCelera(Myers et al、2000)およびJAZZ(Aparicio et al、2002)にわずかな修正を加えて使用した。それらが発表されて以来、多くの特殊なメタゲノムアセンブラが開発されてきた(Koren et al、2011; Laserson et al、2011; Peng et al、2011、2012; Boisvert et al、2012; Namiki et al、2012; Haider et al、2014; Li et al、2016)。しかしながら、バイオインフォマティシャンは、単純な微生物群集と複雑な微生物群集をアセンブリすることの間のギャップを埋めるのに依然として苦労している(総説については、Gevers et al、2012を参照 link)。一方、多くの研究者が複雑なメタゲノムからカバレッジ深度、配列構成、メイトペア情報、および他の基準に基づいて隣接配列をビンに分割してデノボアセンブリを補完する事で豊富な集団のゲノムを単離することに成功した (Hess et al. 2011; Dupont et al. 2012; Iverson et al. 2012; CL Dupont, D Kaul, A Bankevich, DB Rusch, RA Richter, J Zhang, J Stuzka, V Montel, A Young, AE Allen, in prep.) 。しかし、メタゲノムアセンブリの高い断片化は、ビンニングの精度と特定のビンに起因するゲノムの連続性の両方に悪影響を及ぼすため、このアプローチはしばしば困難に直面する。したがって、より良いアセンブラの開発は、メタゲノミクスにおける重要な目標であり続ける。

 シングルセル(Kashtan et al, 2014)およびTruSeq Synthetic Long Read(TSLR)(Sharon et al, 2015)などの最近の技術の応用は、様々な微生物群落内の関連菌株の膨大なmicrodiversityを明らかにした。株はゲノム配列の大部分を共有しているが、それらはしばしば突然変異、mobile elementの挿入、ゲノムリアレンジメント、または水平方向の遺伝子導入から生じる有意な変異を有する。例えば、シングルセルシーケンシングは、wildのプロクロロコッカス(地球上で最も豊富な光合成細菌)集団が何百もの異なるsubpopulations(5%未満の位置で異なる)の「連合」として見ることができることを明らかにした(Kashtan et al link、 2014; Biller et al、2015 link)。さらに、ほぼすべての分析されたシングルセルは、同じsubpopulation由来の他の細胞には見られない少なくとも1つの遺伝子カセットを保有していた。 TSLRを使用することによって、Sharonら(2015)は堆積物サンプルの最も豊富な種が関連した数十の株によって表されることを示した。さらに、研究者らは、このmicrodiversityが、ショートリードライブラリーからの再構築の貧弱さが原因であると主張した。しかし、microdiversityは、以下で論じる多くのメタゲノムアセンブリの課題のうちの1つにすぎない。

 第一に、微生物試料中の様々な種の広く異なる存在量レベルは、異なるゲノムにわたって非常に不均一なcoverageをもたらす。さらに、典型的なメタゲノムデータセットにおける大部分の種の網羅率は、典型的な培養サンプルのシーケンシングプロジェクトと比べてはるかに低い。結果として、coverageが高くそしてかなり均一なゲノムの標準的なアセンブリ技術は、メタゲノムアセンブリで断片化とエラーを起こしやすい。

 第二に、微生物群集内の様々な種はしばしば高度に保存されたゲノム領域を共有している。アセンブリを複雑にし、コンティグを断片化することに加えて、そのような「interspecies repeats」は、ほとんどの種の低いcoverageと共に、ゲノム間アセンブリエラーを引き起こす可能性がある。

 第三に、微生物試料中の多くの細菌種は株混合物、すなわち多様な存在量を有する複数の関連株によって表される(Biller et al, 2014、Kashtan et al, 2014、Rosen et al, 2015、Sharon et al, 2015)。メタゲノミクスの分野以外の様々な研究が、高度に多型の真核生物ゲノム内で2つのhaplomes をアセンブリするという同様の課題に広く取り組んできたが (Dehal et al. 2002; Vinson et al. 2005; Donmez and Brudno 2011; Kajitani et al. 2014; Safonova et al. 2015)、多くのclosely relatedな細菌株の集まりは、独特の計算上の難題を伴う多少異なる問題である。いくつかの研究は複雑な株変異体の同定に向けた最初のステップを記載しているが(Koren et al、2011; Peng et al、2011; Nijkamp et al、2013)、一般的なメタゲノムアセンブリツール(Boisvert et al、2012; Peng et al、2012; Li et al、2016)には、依然として、高レベルのmicrodiversityを有する株混合物をアセンブリするための初歩的な手順しか含まれていない。

 メタゲノミクスの分野以外の用途ではあるが、上記の各課題はすでにSPAdesアセンブリツールキットの開発の過程で解決されていることに注意する。 SPAdesは当初、シングルセルアセンブリの重要な課題の1つである、不均一なカバレッジのデータセットアセンブリするために開発された(Bankevich et al、2012 link; Nurk et al、2013)。 SPAdesのexSPAnder repeat resolutionモジュール(Prjibelski et al, 2014、; Vasilinetc et al, 2015; Antipov et al, 2016)は、様々な技術でシーケンシングされた複数のライブラリーを組み合わせることでゲノムリピートを正確に解決するべく開発された。最後に、dipSPAdes(Safonova et al。2015)は、高度に多型の2倍体ゲノム内で2ハプローム混合物を組み立てるという課題に取り組むために開発された。

 最近開発されたこれらのSPAdesツールは、困難なアセンブリ問題に対処しているが、メタゲノムアセンブリは、データセットサイズに関して、おそらく他のほとんどのDNAシーケンスプロジェクトを邪魔するさらに困難な問題である。SPAdesはメタゲノミクス応用のために設計されていなかったという事実にもかかわらず、様々なグループがSPAdesをメタゲノミクス研究に首尾よく適用した (McLean et al. 2013; Nurk et al. 2013; Coates et al. 2014; Cotten et al. 2014; Bertin et al. 2015; García-López et al. 2015; Kleigrewe et al. 2015; Kleiner et al. 2015; Miller et al. 2016; Tsai et al. 2016; Xie et al. 2016)。しかし、実際にSPAdesはシアノバクテリアフィラメントのような複雑さの低いメタゲノムのアセンブリ(Coates et al。2014 link)やランダムに選択された少数の細菌細胞のMDA増幅混合物(Nurk et al、2013)には有効だが、複雑な細菌群集ではその性能は低下する。

 本著者らの新しいソフトウェアmetaSPAdesは、新しいアルゴリズムのアイデアとSPAdesツールキットの実績あるソリューションを組み合わせて、メタゲノムアセンブリのさまざまな課題に対処する。以下では、metaSPAdesで使用されているアルゴリズム的アプローチを説明し、最先端のメタゲノム・アセンブラIDBA-UD(Peng et al、2012)、Ray-Meta(Boisvert et al、2012)、およびMEGAHIT(Li et al、2015)と比較する。

metaSPAdesパイプラインの概要
metaSPAdesはまずSPAdesを使用してすべてのリードのde Bruijnグラフを作成し、それをさまざまなグラフ単純化手順を使用してアセンブリグラフに変換し、メタゲノム内の長いゲノムフラグメントに対応するアセンブリグラフ内のパスを再構築する(Bankevich et al、2012; Nurk et al、2013)。metaSPAdesは広範囲のカバレッジデプスにわたって機能し、アセンブリの精度と連続性の間のトレードオフを維持しようとする。Microdiversityの課題に応えて、metaSPAdesは、株混合のコンセンサスバックボーンを再構築することに焦点を当てているため、まれな株に対応する株固有の特徴を無視する。

ベンチマークの課題
ゲノムアセンブラは通常、さまざまな測定基準を使用して既知のリファレンスゲノムを持つ分離株でベンチマークされている(Salzberg et al、2012; Gurevich et al、2013)。メタゲノムアセンブラベンチマークは、さらに複雑な複雑な微生物群集で、またリファレンスメタゲノムを利用することができないため、より難しい仕事である。

 この問題に取り組むための1つのアプローチは、メタゲノム中のいくつかのゲノムのclosely relatedなゲノムを同定することに依拠している(Koren et al、2011; Treangen et al、2013)。しかし、この方法は以下の2つの理由、(1)メタゲノムにおけるclosely relatedなリファレンスゲノムはごく一部の種についてしか利用できず、(2)メタゲノムにおける同定されたリファレンスとそれらの対応物との間の相違はしばしばアセンブリエラーとして誤解される (see “Analysis of the HMP Dataset” in the Supplemental Material)、によって制限される。メタゲノミクスアセンブラベンチマークするためのもう1つのアプローチは、既知のコミュニティメンバーの合成データセットを使用することである。そのようなデータセットは、既知の細菌のゲノムの混合物のシーケンシング(Turnbaugh et al、2007; Shakya et al、2013)、単離したサンプルのシーケンシングデータを混合する(Mavromatis et al、2007)、またはリファレンス配列からシミュレートする (Richter et al. 2008; Mende et al. 2012) ことによって得ることができる。しかし、合成データセットはさまざまなベンチマークの取り組みに役立つことが証明されているが、通常、実際のメタゲノムほど複雑ではない(Koren et al、2011; Peng et al、2012)。

3つの人気のメタゲノムアセンブラIDBA-UD(Peng et al。2012)、Ray-Meta(Boisvert et al。2012)、およびMEGAHIT(Li et al。2015)を使い、多様な合成データセットとリアルデータセットを使い、metaQUAST(Mikheenko et al,

2016)によってmetaSPAdesのベンチマークを行なった。データセットは、補足資料の「データの前処理」の説明に従って前処理されている。

 

 HP

http://cab.spbu.ru/software/meta-spades/

 


 

インストール

依存

SPAdes requires a 64-bit Linux system or Mac OS and Python (supported versions are Python2: 2.4–2.7, and Python3: 3.2 and higher) to be pre-installed on it. To obtain SPAdes you can either download binaries or download source code and compile it yourself.

本体 Github

#3.13.1 linux
wget https://github.com/ablab/spades/releases/download/v3.13.1/SPAdes-3.13.1-Linux.tar.gz
tar -zxf SPAdes-3.13.1-Linux.tar.gz
cd SPAdes-3.13.1-Linux/

#3.13.1 darwin
wget https://github.com/ablab/spades/releases/download/v3.13.1/SPAdes-3.13.1-Darwin.tar.gz


#bioconda
conda install -c bioconda -y sppades

./metaspades.py 

# ./metaspades.py 

/root/.pyenv/libexec/pyenv: line 44: cd: spades.py: Not a directory

SPAdes genome assembler v3.13.1 [metaSPAdes mode]

 

Usage: ./metaspades.py [options] -o <output_dir>

 

Basic options:

-o <output_dir> directory to store all the resulting files (required)

--iontorrent this flag is required for IonTorrent data

--test runs SPAdes on toy dataset

-h/--help prints this usage message

-v/--version prints version

 

Input data:

--12 <filename> file with interlaced forward and reverse paired-end reads

-1 <filename> file with forward paired-end reads

-2 <filename> file with reverse paired-end reads

-s <filename> file with unpaired reads

--merged <filename> file with merged forward and reverse paired-end reads

--pe<#>-12 <filename> file with interlaced reads for paired-end library number <#> (<#> = 1)

--pe<#>-1 <filename> file with forward reads for paired-end library number <#> (<#> = 1)

--pe<#>-2 <filename> file with reverse reads for paired-end library number <#> (<#> = 1)

--pe<#>-s <filename> file with unpaired reads for paired-end library number <#> (<#> = 1)

--pe<#>-m <filename> file with merged reads for paired-end library number <#> (<#> = 1)

--pe<#>-<or> orientation of reads for paired-end library number <#> (<#> = 1; <or> = fr, rf, ff)

--s<#> <filename> file with unpaired reads for single reads library number <#> (<#> = 1)

--pacbio <filename> file with PacBio reads

--nanopore <filename> file with Nanopore reads

--tslr <filename> file with TSLR-contigs

 

Pipeline options:

--only-error-correction runs only read error correction (without assembling)

--only-assembler runs only assembling (without read error correction)

--continue continue run from the last available check-point

--restart-from <cp> restart run with updated options and from the specified check-point ('ec', 'as', 'k<int>', 'mc', 'last')

--disable-gzip-output forces error correction not to compress the corrected reads

--disable-rr disables repeat resolution stage of assembling

 

Advanced options:

--dataset <filename> file with dataset description in YAML format

-t/--threads <int> number of threads

[default: 16]

-m/--memory <int> RAM limit for SPAdes in Gb (terminates if exceeded)

[default: 250]

--tmp-dir <dirname> directory for temporary files

[default: <output_dir>/tmp]

-k <int,int,...> comma-separated list of k-mer sizes (must be odd and

less than 128) [default: 'auto']

--phred-offset <33 or 64> PHRED quality offset in the input reads (33 or 64)

[default: auto-detect]

 

 

実行方法

 ペアエンドfastqを指定する(singleのfastqはサポートされていない)。

metaspades.py -1 pair_R1.fq -2 pair_R2.fq -o metaspades \
-t 24 -k auto -m 250

#spadesの"-meta"modeでもランできる。
spades.py -1 pair_R1.fq -2 pair_R2.fq -o metaspades \
-t 24 -k auto -m 250 --meta
  • -1    file with forward paired-end reads
  • -2    file with reverse paired-end reads
  • -m   RAM limit for SPAdes in Gb (terminates if exceeded)
    [default: 250]
  • -k    comma-separated list of k-mer sizes (must be odd and less than 128) [default: 'auto']
  • -t    number of threads
    [default: 16]
  • --meta    this flag is required for metagenomic sample data

SPAdesとは異なり、--carefulや--cov-cutoffフラグは提供されていない。

metaspadesはアダプター配列の残存に敏感なので、適切に前処理されたデータセットを使用する。

 

 

ロングリードを使いmetaspadesでハイブリッドアセンブリを行うことはできますが、現在試験中でパフォーマンスは保障されていないことに注意してください。

 引用

metaSPAdes: a new versatile metagenomic assembler

Nurk S, Meleshko D, Korobeynikov A, Pevzner PA

Genome Res. 2017 May;27(5):824-834

 

関連

 

 

複雑さの高いデータセットでは特にメモリーを要求されます。一番単純にはCPU指定数を減らせばメモリ使用量は減りますが、使用スレッドを減らしすぎると時間がかかり実効性が低くなります。複雑さの高いサンプルでは、実行時間と利用可能なメモリのバランスを考え、さらにdigital normalizationやGCsplitによる負荷の軽減、error correctionの工夫など使い分ける必要が出てきます。

 

 

 

https://www.biostars.org/p/268685/

リファレンスフリーで低メモリかつ高速にSNVとsmall indelを予測する DiscoSnp ++

 

 次世代シーケンス(NGS)データは生命メカニズムへの前例のないアクセスを提供する。特に、これらのデータは染色体、個体または種間の遺伝的差異を評価することを可能にする。そのような多型は、農学、環境または医学における多数の用途を有する生物学の多くの局面における基本的な情報源を表す。
 NGS技術によるシーケンシングの民主化に伴い、SNPまたはインデルなどの遺伝的差異を決定することは今や日常的な作業となっている。そのような多型を予測するために設計された多数のソフトウェアが存在する。ほとんどの場合、これらの方法はGATK [ref.6]やSamTools [ref.14]の場合のようにシーケンスリードマッピングすることによって、あるいはDISCOVAR [ref.22](紹介)やFERMI [ref.12](紹介)のように部分アセンブリマッピングすることによるリファレンスゲノムの使用に基づいている。基本的に、それらは最初にシーケンスをリファレンスゲノムにマッピングし、そして第二段階でそれらはリファレンス遺伝子座を走査して各遺伝子座についてリファレンス配列とマッピングされた配列との間の差異を分析する。
 これらの方法は広く受け入れられ広く使用されている。しかしながら、それらは2つの重大な欠点を提示する。最初にそれらは非常に敏感である。リファレンスゲノムのrepetitive regionは、十分な信頼度でマッピングすることが困難である。これらのリピート領域から検出された多型は、マッピングされたリードの定量が誤っており、さらにリピートの出現間の差異が誤って多型として解釈される可能性があるため、誤っている可能性がある。第二に、それらは高品質のリファレンスゲノムを必要とするという事実に苦しんでいる。この明白で強い条件は、十分に研究された種への適用を制限する。(一部略)
 うまくアセンブリされたリファレンス配列を必要とせずに、シーケンシングリードから直接SNPおよびindelを検出するリファレンスフリーの方法に対する重要かつ絶え間ない要求が依然としてある。 [ref.23]のように、代替的な方法は最初にリードをアセンブルしてリファレンスにその配列をマッピングする。しかしながら、そのような方法は、アセンブリおよびマッピングの両方の問題を累積する。ここでは、このような方法をハイブリッド戦略と呼ぶ。

 いくつかの方法[ref.19、8、11、18、15、21]が多型の新規検出のために提案された。これらすべての方法は、de Bruijnグラフ、すなわち頂点の集合がリードに含まれる長さk(k-mers)のワードの集合に対応し、2つのkの間に辺がある有向グラフの使用に基づいている。それらがk - 1ヌクレオチド上で完全に重複するならばこのデータ構造では、多相は「バブル」と呼ばれる認識可能なパターンを生成する。これらのツールは、その起源を解読するためにそのようなバブルを検出し分析する(シーケンシングエラー、不正確なリピートによる多型、実際のSNPまたはindel)。
 ここでは、DiscoSnp ++を提示する。これは、スキル、計算上のニーズ、および結果のクオリティの点で、他のリファレンスフリーの方法よりも優れた手法である。その主な特徴は、1、予測される変異型の範囲: 分離されているかどうかにかかわらずSNP、small indel。2、数十億のリードを6 GB以下のRAMメモリで処理できる非常に少ないメモリ使用量、3、その高い実行速度、4、優れたprecisionと recall、5、予測された各バリアントに割り当てられたfaithful score、6、出力ファイルはfastaおよびVCFフォーマットで、簡単にダウンストリーム分析に使用可能、7、1つだけも含み、任意の数のリードセットに適用できる、である。
 DiscoSnp ++はDiscoSnpの全く新しいバージョンである[ref.21]。このツールはGATBライブラリ[ref.7](Github)を使用して最初から再実装されたため、実行時間が大幅に短縮され、メモリ使用量が削減された。シーケンシングエラーをよりよく排除し、新しい種類のバリアントを検出するために、アルゴリズムモデルが再検討された。DiscoSnp ++は、必要に応じてその予測をリファレンスゲノムにマッピングした後に、一般的に使用されているVCFフォーマットで予測結果を出力する。この最後の点は、リファレンスのないアプローチに対して直感に反するように見えるかもしれないが、de novo予測したバリアントの位置を決めるためにはリファレンスゲノムが特に役立つ。これは、リファレンスがあまりにもひどくアセンブリされているか、またはシーケンシングされた種から遠すぎる場合には特に当てはまる。さらに、優れたリファレンスゲノムの場合でさえ、リファレンスフリーアプローチを用いたバリアント予測および遺伝子型決定は、リファレンス対立遺伝子によって決してバイアスされない。

 

DiscoSnpに関するツイート

 

インストール

ubuntu18.04でテストした。

依存

本体 Github

git clone --recursive https://github.com/GATB/DiscoSnp.git
cd DiscoSnp
sh INSTALL

#bioconda(mac osにも対応)
conda install -c bioconda -y discosnp

> ./run_discoSnp++.sh

$ ./run_discoSnp++.sh 

You must provide at least one read set (-r|--fof)

 ************

 *** HELP ***

 ************

run_discoSnp++.sh, a pipelining kissnp2 and kissreads for calling SNPs and small indels from NGS reads without the need of a reference genome

Version 2.3.X

Usage: ./run_discoSnp++.sh -r read_file_of_files [OPTIONS]

MANDATORY

-r|--fof <file name of a file of file(s)>

The input read files indicated in a file of file(s)

Example: -r bank.fof with bank.fof containing the two lines 

data_sample/reads_sequence1.fasta

data_sample/reads_sequence2.fasta.gz

 

OPTIONS

-k | --k_size value <int value>

Set the length of used kmers. Must fit the compiled value.

Default=31

-c | --min_coverage value <int value>

Set the minimal coverage per read set: Used by kissnp2 (don't use kmers with lower coverage) and kissreads (read coherency threshold).

This coverage can be automatically detected per read set (in this case use "auto" or specified per read set, see the documentation.

Default=3

-C | --max_coverage value <int value in 0, 1 or 2>

Set the maximal coverage for each read set: Used by kissnp2 (don't use kmers with higher coverage).

Default=2^31-1

-b | --branching value. 

0: forbid variants for which any of the two paths is branching (high precision, lowers the recall in complex genomes).

Default value

1: (smart branching) forbid SNPs for which the two paths are branching (e.g. the two paths can be created either with a 'A' or a 'C' at the same position

2: No limitation on branching (lowers the precision, high recall)

-s | --symmetrical value <int value>

In -b 2 mode only: maximal number of symmetrical crossroads traversed while trying to close a bubble. Default: no limit

-g | --graph <file name>

Reuse a previously created graph (.h5 file) with same prefix and same k and c parameters.

-X Stop discoSnp++ right after variant calling - the output is only a fasta file with no coverage information.

-D | --deletion_max_size <int>

discoSnp++ will search for deletions of size from 1 to D included. Default=100

-a | --ambiguity_max_size <int>

Maximal size of ambiguity of INDELs. INDELS whose ambiguity is higher than this value are not output  [default '20']

-P | --max_snp_per_bubble <int>

discoSnp++ will search up to P SNPs in a unique bubble. Default=3

--fof_mapping <file name of a file of file(s)>

If this option is used this fof is used when mapping back reads on the predicted variants instead of the original fof file provided by -r|--fof option

-p | --prefix <string>

All out files will start with this prefix. Default="discoRes"

-l | --no_low_complexity

Remove low complexity bubbles

-T | --contigs

Extend found polymorphisms with contigs (default: extend with unitigs)

-d | --max_substitutions <int>

Set the number of authorized substitutions used while mapping reads on found SNPs (kissreads). Default=1

-n | --no_genotype

Do not compute the genotypes

-u | --max_threads <int>

Max number of used threads. 0 means all threads

default 0

 

REFERENCE GENOME AND/OR VCF CREATION OPTIONS

-G | --reference_genome <file name>

Reference genome file (fasta, fastq, gzipped or nor). In absence of this file the VCF created by VCF_creator won't contain mapping related results.

-R

Use the reference file also in the variant calling, not only for mapping results

-B | --bwa_path <directory name>

bwa path. e.g. /home/me/my_programs/bwa-0.7.12/ (note that bwa must be pre-compiled)

Optional unless option -G used and bwa is not in the binary path.

-e Map variant predictions on reference genome with their unitig or contig extensions.

Useless unless mapping on reference genome is required (option -G). 

 

-w Wraith mode: only show all discoSnp++ commands without running them

-v <0 or 1>

Verbose 0 (avoids progress output) or 1 (enables progress output) -- default=1.

-h | --help

Prints this message and exist

 

Any further question: read the readme file or contact us via the Biostar forum: https://www.biostars.org/t/discosnp/

run_discoSnpRad.sh -help

$ run_discoSnpRad.sh -help

 

 

 

run_discoSnpRad.sh, pipelining kissnp2 and kissreads and clustering per locus for calling SNPs and small indels from RAD-seq data without the need of a reference genome

Version 2.3.X

Usage: ./run_discoSnpRad.sh -r read_file_of_files [OPTIONS]

MANDATORY:

-r read_file_of_files

    Example: -r bank.fof with bank.fof containing the two lines 

data_sample/reads_sequence1.fasta

data_sample/reads_sequence2.fasta.gz

-S **absolute** path to short_read_connector directory, containing the short_read_connector.sh'' file. 

-Note1: short read connector must be compiled.

-Note2: if this option is missing, discoSnpRad will still however provide a fasta file containing SNPs and INDELS, that won't be clustered by locus

DISCOSNPRAD OPTIONS:

-g: reuse a previously created graph (.h5 file) with same prefix and same k and c parameters.

-b value. 

1: (smart branching) forbid SNPs for which the two paths are branching (e.g. the two paths can be created either with a 'A' or a 'C' at the same position Default value

2: No limitation on branching (lowers the precision, high recall)

-s value. In b2 mode only: maximal number of symmetrical croasroads traversed while trying to close a bubble. Default: no limit

-D value. discoSnpRad will search for deletions of size from 1 to D included. Default=100

-a value. Maximal size of ambiguity of INDELs. INDELS whose ambiguity is higher than this value are not output  [default '20']

-P value. discoSnpRad will search up to P SNPs in a unique bubble. Default=5

-p prefix. All out files will start with this prefix. Default="discoRad"

-l: remove low complexity bubbles

-k value. Set the length of used kmers. Must fit the compiled value. Default=31

-c value. Set the minimal coverage per read set: Used by kissnp2 (don't use kmers with lower coverage) and kissreads (read coherency threshold). This coverage can be automatically detected per read set or specified per read set, see the documentation. Default=auto

-C value. Set the maximal coverage for each read set: Used by kissnp2 (don't use kmers with higher coverage). Default=2^31-1

-d value. Set the number of authorized substitutions used while mapping reads on found SNPs (kissreads). Default=1

-u: max number of used threads

-v: verbose 0 (avoids progress output) or 1 (enables progress output) -- default=1.

 

-h: Prints this message and exist

Any further question: read the readme file or contact us via the Biostar forum: https://www.biostars.org/t/discosnp/

 

 

テストラン

例えば1セットのリードを調べ、不完全SNVとsmall indelを検出する。fastq/fastaのファイル名をテキストで指定する。

git clone https://github.com/GATB/DiscoSnp.git
cd DiscoSnp/test/
run_discoSnp++.sh -r fof.txt -T
  • -r    read_file_of_files. Example: -r bank.fof with bank.fof containing the two lines
    data_sample/reads_sequence1.fasta
    data_sample/reads_sequence2.fasta.gz
  • -T    extend found polymorphisms with contigs

 

 

リードに加えリファレンスを指定する。さらにこのリファレンスを使ってバリアントコールを行う。bwaが必要。

run_discoSnp++.sh -r test/fof.txt -T -G test/reference_genome.fa -R
  • -G    reference genome file (fasta, fastq, gzipped or nor). In absence of this file the VCF created by VCF_creator won't contain mapping related results.
  • -R    use the reference file also in the variant calling, not only for mapping results

リファレンスを与えた場合でも、与えなかった場合でも、このデータセットでコールされるのは不完全なSNV三箇所になる。

f:id:kazumaxneo:20190514140437j:plain

 

リファレンスを与えると、リファレンスにmappingしてリファレンスでのポジションを報告してくれる。

引用

DiscoSnp ++ : de novo detection of small variants from raw unassembled read set(s)

Peterlongo, P., Riou, C., Drezen, E., Lemaitre, C.

bioRxiv preprint first posted online Oct. 27, 2017

 

Reference-free detection of isolated SNPs

Uricaru R., Rizk G., Lacroix V., Quillery E., Plantard O., Chikhi R., Lemaitre C., Peterlongo P. (2014).

Nucleic Acids Research 43(2):e11