コア遺伝子のアミノ酸配列を使って系統解析を行う bcgTree

2020 4/1  関連ツールリンク追加

　DNAシーケンシングデータによる生物の進化的および分類学的関係の再現は、バクテリアにおいて長い歴史を持つ（Cavalier-Smith、1993; Woese and 33Fox、1977; Woese、1987）。バクテリアは形態学的に区別し分類するのが難しく、DNAバーコードと分子系統学が細菌の関係を決定しようとする研究者のための選択肢となっている。しかしながら、DNAシークエンスの使用による系統発生的関係を解決することは困難な作業であり得る。適切な遺伝マーカーの選択、分類を区別するのに十分な情報、比較を有効かつ決定的にするために十分な相同性を有するものを選ぶのが正しい再構成のために不可欠になる（Capella-Gutierrez et al、402014; Wu et al。しばしば、情報と保存の間のこのようなトレードオフを補うために、低レベル（strain/species/genus）および高レベル（family/class/order/phylum）の異なるマーカーが系統解析に使用される（Capella-Gutierrez et al 、2014; Wu et al。、2013）。
　バクテリアでは、現在、16S rDNAがほとんどの系統学的および生態学的研究のための選択肢として、比類のない普遍的に適用されるマーカーである。このマーカーはいくつかの可変領域および保存領域が存在し、情報の豊富な領域の増幅が可能でclosely relatedな分類群の区別を可能にする。 16S rDNAはすべての原核生物種にわたってよく保存されており、phyla間の比較が可能である。しかし、このアプローチには、ハイレベル分析とローレベル分析（Capella-Gutierrez et al、2014; Wu et al、2013）の両方に関して、その限界がある。 バクテリアの系統発生は依然としてその基礎の枝の部分で未解決であり、これらは単一のマーカーを用いて解決される可能性は低い。さらに、16S配列は、生態学的機能を有する特定の遺伝的に異なる株を区別するためのこのマーカーの可能性を排除する（Jaspers and Overmann、2004）、すなわち、大規模なゲノム再構成にもかかわらず株間で同一であり得る。さらに、リボソームrDNAは各ゲノムに複数コピー存在し、コピーが同一ゲノム内で可変なこともある（Větrovskýand Baldrian、2013）。両方の効果は、特に、機能的または生態学的な解析において、また相互主義的または病原性の宿主 - バクテリア関係と同様に、分類学的帰属の解釈を混乱させる可能性がある。
　16S rDNAベースの系統発生解析は、バクテリアの同一性と進化的関係を理解する上で非常に重要であり、1つのマーカー配列で系統樹を計算することには欠陥がある。ハイスループットシークエンシング技術における現在の進歩により、この単一のマーカーコンベンションよりも広範な分析が可能になる。バクテリアゲノムは、大多数の真核生物と比べると、通常130kbpから14Mbpの範囲にサイズが制限されている。bp当たりのコストダウンと小さなゲノムサイズにより、限られた財源しか有さないワーキンググループであっても、完全なバクテリアゲノムのシーケンシングが可能になっているKeller et al、2014; Metzker、2009）。これはまた、シーケンシングされ、depositされるコンプリートなバクテリアゲノムの増加を導く（Pruitt et al、2007; Uchiyama et al。、2013）。
　しかしながら、ゲノム全体から系統発生を導き出すことは、それ自体チャレンジングである。 例えば、大抵の場合、明確な類似性がない大きなゲノム領域が存在する。 また、相同性を有する領域については、下流の系統分析のために抽出する必要があり、これは広範なバイオインフォマティクスの専門知識を必要とするプロセスである。この課題に取り組む一つの方法は、あらかじめアライメントが計算されたタイトなゲノムクラスター（ATGC、Novichkov et al、2009）のデータベースを構築し、高分解能なmicro-evolutionary analysesを可能にすることである。しかし、利用可能なゲノムのうちのほんの一部しか含まれておらず、ユーザが提供もできないため、このアプローチは限定されている。この問題の1つの解決策は、関心のある生物の大多数に存在する保存された領域に集中することである（Ciccarelli et al、822006）。
　ここでは、全ゲノムデータのアミノ酸配列から必須の107個のシングルコピー遺伝子のセットを同定および抽出するツールであるbcgTreeを報告する。ここで使用される「essential core genes」の定義は、Dupont et al（2012）（pubmed）の研究に基づいており、生物学的ではなく統計的な議論に従う。本発明者らのソフトウェアは、コア遺伝子配列データを自動的にコンパイルし、これを用いて分割最尤解析を用いて系統樹を再構築する。 妥当性確認の目的のために、本発明者らは、Lactobacillalesで現在利用可能なほとんどのゲノムを含む、90のLactobacillus属の系統発生の解明にbcgTreeを適用した。さらに、評価目的のために、高レベルおよび低レベルの系統発生分析テストを対応する16S rDNAツリーと直接比較した。

bcgTreeパイプライン。論文より転載。

依存

• perl5、Java Runtime Environment (version8)

以下のperlモジュールとjavaフレームワークはレポジトリに含まれている。

• Getopt::Long
• Pod::Usage
• FindBin
• File::Path
• File::Spec
• SwingX 1.6.5 (LGPL)

必要な外部ツール

• hmmsearch (HMMER version 3.1b1) - Eddy et al, 2010. HMMER3: a new generation of sequence homology search software.
• muscle (v3.8.31) - Edgar et al, 2004. MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res. 32, 1792–1797.
• Gblocks (version 0.91b) - Castresana et al, 2000. Selection of conserved blocks from multiple alignments for their use in phylogenetic analysis. Mol. Biol. Evol. 17, 540–552.
• RAxML (version 8.2.4) - Stamatakis et al, 2014. RAxML version 8: A tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics 30, 1312–1313.

#Anaconda環境にて
conda install -c bioconda raxml==8.2.4 
conda install -c bioconda gblocks==0.91b 
conda install -c bioconda muscle==3.8.31 
conda install -c bioconda hmmer==3.1b2

 本体 github

git clone https://github.com/molbiodiv/bcgTree.git
# Now you can run bcgTree.pl
bcgTree/bin/bcgTree.pl --help
# Or start the java GUI
cd bcgTree/bcgTreeGUI
java -jar bcgTree.jar

 CUI版のヘルプ

> perl bcgTree/bin/bcgTree.pl --help

$ perl bcgTree.pl -help

Usage:

      $ bcgTree.pl [@configfile] --proteome bac1=bacterium1.pep.fa --proteome bac2=bacterium2.faa [options]

Options:

    [@configfile]            Optional path to a configfile with @ as prefix.

                             Config files consist of command line parameters

                             and arguments just as passed on the command

                             line. Space and comment lines are allowed (and

                             ignored). Spreading over multiple lines is

                             supported.

    --proteome <ORGANISM>=<FASTA> [--proteome <ORGANISM>=<FASTA> ..]

                             Multiple pairs of organism and proteomes as

                             fasta file paths

    [--outdir <STRING>]      output directory for the generated output files

                             (default: bcgTree)

    [--help]                 show help

    [--version]              show version number of bcgTree and exit

    [--check-external-programs]

                             Check if all of the required external programs

                             can be found and are executable, then exit.

                             Report table with program, status (ok or

                             !fail!) and path. If all external programs are

                             found exit code is 0 otherwise 1. Note that

                             this parameter does not check that the paths

                             belong to the actual programs, it only checks

                             that the given locations are executable files.

    [--hmmsearch-bin=<FILE>] Path to hmmsearch binary file. Default tries if

                             hmmsearch is in PATH;

    [--muscle-bin=<FILE>]    Path to muscle binary file. Default tries if

                             muscle is in PATH;

    [--gblocks-bin=<FILE>]   Path to the Gblocks binary file. Default tries

                             if Gblocks is in PATH;

    [--raxml-bin=<FILE>]     Path to the raxml binary file. Default tries if

                             raxmlHPC is in PATH;

    [--threads=<INT>]

        Number of threads to be used (currently only relevant for raxml).

        Default: 2 From the raxml man page: PTHREADS VERSION ONLY! Specify

        the number of threads you want to run. Make sure to set "-T" to at

        most the number of CPUs you have on your machine, otherwise, there

        will be a huge performance decrease!

    [--bootstraps=<INT>]

        Number of bootstraps to be used (passed to raxml). Default: 100

    [--min-proteomes=<INT>]

        Minimum number of proteomes in which a gene must occur in order to

        be kept. Default: 2 All genes with less hits are discarded prior to

        the alignment step. This option is ignored if --all-proteomes is

        set.

    [--all-proteomes]

        Sets --min-proteomes to the total number of proteomes supplied.

        Default: not set All genes that do not hit all of the proteomes are

        discarded prior to the alignment step. If set --min-proteomes is

        ignored.

    [--hmmfile=<PATH>]

        Path to HMM file to be used for hmmsearch. Default:

        <bcgTreeDir>/data/essential.hmm

    [--raxml-x-rapidBootstrapRandomNumberSeed=<INT>]

        Random number seed for raxml (passed through as -x option to raxml).

        Default: Random number in range 1..1000000 (see raxml command in log

        file to find out the actual value). Note: you can abbreviate options

        (as long as they stay unique) so --raxml-x=12345 is equivalent to

        --raxml-x-rapidBootstrapRandomNumberSeed=12345

    [--raxml-p-parsimonyRandomSeed=<INT>]

        Random number seed for raxml (passed through as -p option to raxml).

        Default: Random number in range 1..1000000 (see raxml command in log

        file to find out the actual value). Note: you can abbreviate options

        (as long as they stay unique) so --raxml-p=12345 is equivalent to

        --raxml-p-parsimonyRandomSeed=12345

    [--raxml-aa-substitiution-model "<MODEL>"]

        The aminoacid substitution model used for the partitions by RAxML.

        Valid options for RAxML 8.x are: DAYHOFF, DCMUT, JTT, MTREV, WAG,

        RTREV, CPREV, VT, BLOSUM62, MTMAM, LG, MTART, MTZOA, PMB, HIVB,

        HIVW, JTTDCMUT, FLU, STMTREV, DUMMY, DUMMY2, AUTO, LG4M, LG4X,

        PROT_FILE, GTR_UNLINKED, GTR bcgTree will not check whether the

        provided option is valid but rather pass it to RAxML literally.

        Default: AUTO

    [--raxml-args "<ARGS>"]

        Arbitrary options to pass through to RAxML. The ARGS part should be

        in quotes and is appended to the RAxML command as given.

ラン

GUI版を起動する（違いはない。CUI版のラッパー、いわゆるガワ）。

cd bcgTree/bcgTreeGUI/
java -jar bcgTree.jar

Add proteomesから比較する生物種のproteomeのfastaファイル (amino-acids.fasta) を入力する。

バクテリアのprotein.fasta５つを登録したところ。

GIthubより。

SupplementaryにはbcgTreeのランニングタイムの図もあります。

引用

bcgTree: automatized phylogenetic tree building from bacterial core genomes
Ankenbrand MJ, Keller A

Genome. 2016 Oct;59(10):783-791. Epub 2016 May 11.

PhyloPhlAnは以前簡単に紹介しています。

関連

ラージゲノムにも対応したde novo assembly評価ツール QUAST-LG

2019 7/28 help追記、タイトル修正、コマンド例追記

2019 10/20 リンク追加

2020 1/11 インストール修正

　現代のDNAシーケンシング技術は染色体の全配列を読み取ることができない。代わりに、それらはゲノムの異なる部分からサンプリングされた多数のリードを生成する。低コストで高品質の第2世代シーケンシング（次世代シークエンシングまたはNGSとも呼ばれる）の出現により、科学者は以前にシーケンシングれていない数多くの生物を解読することができるようになった。しかしながら、ゲノムシーケンシングの市場は急速に発展しており、NGS技術の支配力も、第3世代（long-read sequencing（LRS）の絶え間ない改善のために、現在疑問視されている。 Pacific BiosciencesとOxford Nanopore TechnologiesのLRS技術によって作成されたリードの長さは、NGSによって生成された長さより数桁高く10万塩基を超えることができる。同時に、これらの技術は、NGS競合他社と比較してエラー率が大幅に高いリードを生成し、かなり高価なままである。従って、両方の種類のDNAシーケンシング技術は重大な欠点があるが、しかし、計算方法で克服することができる。

　ゲノムアセンブリソフトウェアは、リードをコンティグと呼ばれるより大きな配列に組み合わせる。ゲノムの長いリピートは、完全な染色体の再構成を妨げる。完全な染色体の再構成は、NGS技術によって生成されたショートリードのみを使用するとほとんどの場合不可能である。しかし、ショートリードアセンブリは長距離のメイトペアライブラリー（Chaisson et al、2009; Vasilinetc et al、2015）を用いて有意に改善することができる。典型的なインサートサイズが1キロベース（kb）未満のペアエンドライブラリとは対照的に、現在のメイトペアプロトコルは、より長いインサート（最大20kb）を生成する。アセンブリソフトウェアは、その情報を利用して、長いリピートとカバレッジギャップをまたぎ、コンティグをscaffoldsと呼ばれるより長いシーケンスに拡張する。ロングリードのアセンブリアルゴリズムは、リード長が複雑なリピートを解決するのに十分なので、通常、ペアは必要ない。 LRS技術の最も重要な欠点は、高いベース・エラー率であり、より高いシーケンスデプス（Koren et al。、2017）に向け投資するか、高品質のショートリードを使用してエラーを訂正する（Walker et al、2014）。　

　ゲノムアセンブリの課題に取り組むため、様々なアセンブリツールが異なるヒューリスティックアプローチを使用しており、その結果、生じるコンティグおよびScaffoldsのに大きな違いが生じる。それは、アセンブリを互いに比較するという問題を引き起こす。 Assemblathon1（Earl et al、2011）は、主要ゲノムアセンブリソフトウェアを評価する最初の試みの1つだった。主催者は、未知のシミュレートされたゲノムを使いシミュレーションのショートリードデータセットを提供した。このデータセットは、世界中のアセンブラ開発者とユーザーの17グループによってアセンブリされ、結果として得られた41のアセンブリが慎重に評価された。その後の研究であるAssemblathon 2（Bradnam et al、2013）は同様の競合モデルを使用したが、3つの脊椎動物種からの実際のシーケンシングデータを提供した。残念なことに、両方のアセンブリプロジェクトで使用された評価スクリプトは、ほとんどがそれらの特定のゲノムに集中しており、ルーティンの品質評価に簡単には適用できなかった。 GAGEの研究（Salzberg et al、2012）は、ゲノムアセンブリ評価のゴールデンスタンダードを確立し、それを実装するプログラムを作成した。 GAGEツールは、その後のバクテリアゲノムアセンブラ（Magoc et al、2013）のGAGE-B競争の主な評価ツールとなったQUAST（Gurevich et al。、2013）ソフトウェアでさらに拡張され改良された。メタゲノムの解釈の研究（Sczyrba et al、2017）の最近の批判的評価は、メタゲノム評価のためQUASTの拡張であるMetaQUAST（Mikheenko et al、2016b）を用いた。

　QUASTパッケージは幅広く使用されているが、普遍的なゲノムアセンブリ評価ツールとして扱うことはほとんどない。その理由は、QUAST設計の元々の目的がバクテリアや小型の真核生物アセンブリの評価であり、100メガよりも小さいゲノムに制限されているからである。この研究では、ラージゲノムの配列を解析できる大幅にアップグレードされたQUASTのバージョンであるQUAST-LGを提示する。著者らのツールは、通常のサーバー上で数時間以内に哺乳動物サイズのアセンブリを評価することができる。 QUAST-LGには、トランスポゾンの豊富さなど、真核生物ゲノムの特定の特徴を考慮に入れたスピードアップの改善と新しい品質メトリックの両方が含まれている。 REAPR（Hunt et al、2013）（紹介）、ALE（Clark et al、2013）およびLAP（Ghodsi et al、2013）のような純粋なリファレンスフリーアプローチとは対照的に、元のQUASTは既知のリファレンスゲノムに依存していた。 QUAST-LGは、パイプラインの一部にリファレンスフリーのツールを組み込むことで、新規種の評価を可能にする。しかし、リファレンスフリー解析はQUAST-LGの主要なユースケースではない。

　QUAST-LGの別の改良点は、ゲノムの長いリピートおよび低いカバレッジの領域のために、リファレンスゲノムが未処理のリードから完全に再構築され得ないという事実に基づくアセンブリ上限の概念である。 near-optimalな アセンブリ問題（Bresler et al、2013; Lam et al、2014）に関する以前の研究で、いくつかの理論的仮定の下で、ゲノム再構成を成功させるのに必要なショットガンシーケンシングデータのリード長とカバレッジ深さが推定されている。著者らは反対の問題にアプローチする：データセットが与えられた場合、QUAST-LGは、この特定のリードセットからアセンブリソフトウェアによって理論的に到達可能な完全性および連続性の上限を推定する。さらに、本方法は、データセットがNGSとLRSの両方の技術によって生成された複数のシークエンシングライブラリで構成されていることを考慮している。　

　QUAST-LGの機能性を実証するために、著者らはミディアムサイズおよびラージサイズのゲノムの最先端のゲノムアセンブリツールの能力を評価する。2つのライブラリ、すなわちイルミナペアエンドリードとハイクオリティなイルミナメイトペア、またはイルミナペアエンドリードとロングリードライブラリーのいずれかを含むWGSデータセットを使用する。パシフィックバイオサイエンスの１分子リアルタイム（PacBio SMRT）シークエンシングとOxford Nanoporeシークエンシングの両方の主要LRS技術は、これらのテストデータセットの中に存在する。結果の簡便な再現を可能にするために、この研究で使用されたすべてのデータセットおよびソフトウェアツールは自由に利用可能となっている。

QUAST-LG HP（論文内のアセンブリツール比較に使われたデータと、QUAST-LG resultsのリンクあり（S. cerevisiae、human、Fruit fly、C. elegans））

QUAST-LG – Center for Algorithmic Biotechnology

マニュアル

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

プレゼンテーションPDF 

http://cab.spbu.ru/files/quast/quast-lg/talk_and_poster/Gurevich_QUAST-LG_ISMB2018.pdf

Poster

http://cab.spbu.ru/files/quast/quast-lg/talk_and_poster/QUAST-LG_poster_ISMB2018.png

QUAST-LGに関するツイート

インストール

Linux and Mac OS X are supported.

依存

For the main pipeline:

• Python2 (2.5 or higher) or Python3 (3.3 or higher)
• Perl 5.6.0 or higher
• GCC 4.7 or higher
• basic UNIX tools (make, sh, sed, awk, ar)

For the optional submodules:

• Time::HiRes perl module for GeneMark-ES (needed when using --gene-finding --eukaryote)
• Java 1.8 or later for GRIDSS (needed for SV detection)
• R for GRIDSS (needed for SV detection)
• QUAST draws plots in two formats: HTML and PDF. If you need the PDF versions, make sure that you have installed Matplotlib.

本体　GIthub

バージョン５からQUASTのパッケージに QUAST-LGも含まれるようになった（conda link）。

#依存が多いので仮想環境に入れる
conda create -n quast -y
conda activate quast
conda install -c bioconda -y quast

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

> quast -h   # v5.0.2

# quast -h

WARNING: Python locale settings can't be changed

QUAST: Quality Assessment Tool for Genome Assemblies

Version: 5.0.2

Usage: python /root/.pyenv/versions/miniconda3-4.0.5/bin/quast [options] <files_with_contigs>

Options:

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

-r                <filename>      Reference genome file

-g  --features [type:]<filename>  File with genomic feature coordinates in the reference (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 GeneMarkS (prokaryotes, default) or GeneMark-ES (eukaryotes, use --eukaryote)

    --mgm                             Use MetaGeneMark for gene prediction (instead of the default finder, see above)

    --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)

    --est-ref-size <int>              Estimated reference size (for computing NGx metrics without a reference)

    --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]

    --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 QUAST on the data from the test_data folder, output to quast_test_output

    --test-sv                         Run QUAST with structural variants detection on the data from the test_data folder,

                                      output to quast_test_output

-h  --help                            Print full usage message

-v  --version                         Print version

Online QUAST manual is available at http://quast.sf.net/manual

root@47647a654823:/data# 

ラン

 アセンブリした配列３つを評価する。

quast contig1.fa contig2.fa contig3.fa \
 -R reference.fa.gz \
 -G genes.txt \
 -O operons.txt \
 -1 reads1.fastq.gz -2 reads2.fastq.gz \
 -o quast_test_output \
 -t 8 --large

論文中のデータとしてQUAST-LG  HPに登録されている、cell line HG001 (NA12878)の解析結果を見てみる（ reportへのダイレクトリンク ）。

下の方のExtended reportをクリックするとより詳細なmetricsが表示される。

Icarusを使い、contigのリファレンスへのアライメント結果を確認する。上のView in Icarus contig browserをクリックする。

表示する染色体を選択する。

1番染色体を表示する。

拡大縮小したり、左右にスクロールできる。

ブロックをクリックすれば、アライメントされた領域、ミスアセンブリのポジションなどを確認できる。

上のメニューからは、表示するミスアセンブリの種類を限定できる。ローカルミスアセンブリを表示。

公式ページにも記載されている評価項目を簡単にまとめておく。

リファレンスの有無に関係なく出力されるmetrics

1. コンティグ数
2. 最大コンティグ長（bp）
3. N50、N75（アセンブリの連続性に関するmetrics。N50ならアセンブリの少なくとも50％をカバーする最短コンティグの長さで定義される）
4. L50、L75（こちらはコンティグの数。L50なら、N50を満たす最小contig数として定義される）
5. 予測された遺伝子数（GeneMark.hmm (for prokaryotes、or GlimmerHMM (for eukaryotes), or MetaGeneMark (for metagenomes)
6. Nの数

リファレンスがある時のみ出力されるmetrics

1.  ミスアセンブリの数 (inversions, relocations, translocations, interspecies translocations (metaQUAST only) ).
2. アライメントされたcontigの数とそのトータルサイズ
3. アライメントされなかったcontigの数とそのトータルサイズ
4. ミスマッチとindelの数
5. Genome fraction %（リファレンスのカバー率）
6. Duplication率（the total number of aligned bases in the assembly divided by the total number of those in the reference.）
7. 全長 / 一部カバーされている遺伝子の割合（ユーザーが提供したリファレンスの遺伝子リストを使う）
8. NG50（N50に似ているが、NG50はリファレンスの少なくとも50％をカバーする最短コンティグの長さ）
9. NGA50, a reference-aware version of N50 metric. It is calculated using aligned blocks instead of contigs. Such blocks are obtained after removing unaligned regions, and then splitting contigs at misassembly breakpoints. Thus, NGA50 is the length of a block, such that all the blocks of at least the same length together cover at least 50% of the reference.
10. Complete BUSCO (%)
11. K-mer-based completeness

論文内では、最新のmummer4（紹介）やminimap2（紹介）によるcontigのアライメント、KMC3（紹介）によるk-mer評価を元にしたアセンブリ完全性の評価、BUSCO（紹介）によるハウスキーピング遺伝子のhit率に基づいたアセンブリ評価、などの詳細が書かれています。またSupplementary documentsには、アセンブラのパフォーマンス、ラージゲノム評価時のメモリ使用量などもまとめられています。

追記

eukaryoteのラージゲノム（fungiは"--fungus"）、遺伝子予測、 rRNA予測、buscoによる評価、も行う。

quast contig1.fa contig2.fa contig3.fa contig4.fa contig5.fa contig6.fa \
 -R reference.fa.gz \
 --labels miniasm,redbeans,canu,flye,ra,shasta \
 --eukaryote --large\
 --gene-finding \
 --rna-finding \
 --conserved-genes-finding \
 -1 reads1.fastq.gz -2 reads2.fastq.gz \
 --nanopore ONT.fq.gz --pacbio pacbio.fq.gz\
 -t 50 \
 -o quast_output

引用

Versatile genome assembly evaluation with QUAST-LG
Mikheenko A, Prjibelski A, Saveliev V, Antipov D, Gurevich A

Bioinformatics. 2018 Jul 1;34(13):i142-i150

 
Icarus: visualizer for de novo assembly evaluation.
Mikheenko A, Valin G, Prjibelski A, Saveliev V, Gurevich A

Bioinformatics. 2016 Nov 1;32(21):3321-3323.

関連

MetaBAT

2019 8/28 追記

2019 9/30 metabat2紹介リンク追加

　ハイスループットのメタゲノムショットガンシークエンシングは、環境から採取された微生物群集を直接研究するための強力なツールであり、それによって培養から解放され、また培養から生じる可能性のあるバイアスを回避する。ショートメタゲノムショットガンリードのアセンブリでは、ショートリードアセンブラによって大きなゲノムフラグメント（コンティグ）を組み立てるが、完全長ゲノムの作成にはよく失敗する（Pevzner＆Tang、2001; Pevzner、Tang＆Waterman、2001）。メタゲノムコンティグのメタゲノムビニングによるドラフトゲノム予測は、完全長ゲノムの代替となる（Mande、Mohammed＆Ghosh、2012; Mavromatis et al、2007）。断片にもかかわらず、これらはしばしば個々の種（または異なる株のコンセンサス配列を表す「population genomes」（Imelfort et al、2014））のほぼ完全な遺伝子セットを持ち、ドラフトゲノムを近似している。

　2つのメタゲノムビニング手法が開発されている（Mande、Mohammed＆Ghosh、2012）。supervised（教師あり）のビニングアプローチは、既知のゲノムをリファレンスとして使用し、ビニングのため配列相同性または配列組成類似性に頼っている（Krause et al、2008; Wu＆Eisen、2008）。このアプローチは、多くの微生物が既知のゲノムとclosely relatedな種を持たない環境サンプルではうまく機能しない。対照的に、unsupervised（教師なし）のアプローチは、区別可能な配列組成（Teeling et al、2004b; Yang et al、2010）または種（またはゲノム断片）のco-abundance(Cotillard et al., 2013; Le Chatelier et al., 2013; Nielsen et al., 2014; Qin et al., 2012; Wu & Ye, 2011)、またはその両方(Albertsen et al., 2013（紹介）; Alneberg et al., 2014; Imelfort et al., 2014; Sharon et al., 2013; Wrighton et al., 2012; Wu et al., 2014) をビニングに使っている。近年の研究では、利用可能なサンプルが多数ある場合には、種の co-abundanceのフィーチャが複雑なコミュニティをデコンボリューションするのに非常に有効であることが示されている (Albertsen et al., 2013; Alneberg et al., 2014; Cotillard et al., 2013; Imelfort et al., 2014; Karlsson et al., 2013; Le Chatelier et al., 2013; Nielsen et al., 2014; Sharon et al., 2013) 。最近の少しだけ、CONCOCT（Alneberg et al、2014）とGroopM（Imelfort et al、2014）のような完全自動ビニング方法も報告されている。

　上記のツールの多くは、大きなメタゲノミックデータセットには適していない。この研究では、数千のサンプルから数百万のコンティグをビニングすることができる効率的で完全自動化されたソフトウェアツールMetaBAT（Metagenome Binning with Abundance and Tetra-Nucleotide Frequency）を開発した。テトラヌクレオチド頻度（TNF）とコンティグの存在量確率を組み合わせるための新しい統計的枠組みを用いることにより、MetaBATは高品質のゲノムビンを産生することを示した。

　ビニングの前提条件として、ユーザーは、各サンプルのリードをアセンブルされたメタゲノムに個別にマッピングすることによってBAMファイルを作成する必要がある（論文　図1のステップ1-3 link）。 MetaBATはアセンブリファイル（fastaフォーマット、必須）とソートされたbamファイル（サンプルごとに1つ、任意）を入力として受け入れる。 メタゲノムアセンブリ中の各コンティグペアについて、MetaBATは、テトラヌクレオチド頻度（TNF）および存在量（すなわち、平均ベースカバレッジ）に基づいて確率的距離を計算し、次いで、2つの距離を1つの合成距離に統合する。 全てのペアワイズ距離はマトリックスを形成し、その後、modifyされたk-medoidクラスタリングアルゴリズムに供給され、コンティグを反復的かつ完全にゲノムビンにビニングする（論文　図1）。 

MetaBAT overview。論文図１より転載

インストール

mac os10.13でテストした。

依存

• boost >= 1.55.0
• python >= 2.7
• scons >= 2.1.0
• g++ >= 4.9
• zlib >= 1.2.4
• binutils >= 2.2.2

本体　BitBucket

#Anaconda環境ならcondaで導入可能。version2が入る（paper）。
conda install -c ursky metabat2

#dockerでもランできる。
docker pull metabat/metabat
docker run metabat/metabat:latest runMetaBat.sh

> runMetaBat.sh 

 $ runMetaBat.sh 

/Users/user/.pyenv/versions/anaconda2-4.2.0/bin/runMetaBat.sh <select metabat options> assembly.fa sample1.bam [ sample2.bam ...]

You can specify any metabat options EXCEPT:

  -i --inFile

  -o --outFile

  -a --abdFile

For full metabat options: metabat2 -h

——

> metabat2

 $ metabat2

MetaBAT: Metagenome Binning based on Abundance and Tetranucleotide frequency (version 2.12.1; Aug 31 2017 21:02:54)

by Don Kang (ddkang@lbl.gov), Feng Li, Jeff Froula, Rob Egan, and Zhong Wang (zhongwang@lbl.gov) 

Allowed options:

  -h [ --help ]                     produce help message

  -i [ --inFile ] arg               Contigs in (gzipped) fasta file format [Mandatory]

  -o [ --outFile ] arg              Base file name and path for each bin. The default output is fasta format.

                                    Use -l option to output only contig names [Mandatory].

  -a [ --abdFile ] arg              A file having mean and variance of base coverage depth (tab delimited; 

                                    the first column should be contig names, and the first row will be 

                                    considered as the header and be skipped) [Optional].

  -m [ --minContig ] arg (=2500)    Minimum size of a contig for binning (should be >=1500).

  --maxP arg (=95)                  Percentage of 'good' contigs considered for binning decided by connection

                                    among contigs. The greater, the more sensitive.

  --minS arg (=60)                  Minimum score of a edge for binning (should be between 1 and 99). The 

                                    greater, the more specific.

  --maxEdges arg (=200)             Maximum number of edges per node. The greater, the more sensitive.

  --pTNF arg (=0)                   TNF probability cutoff for building TNF graph. Use it to skip the 

                                    preparation step. (0: auto).

  --noAdd                           Turning off additional binning for lost or small contigs.

  --cvExt                           When a coverage file without variance (from third party tools) is used 

                                    instead of abdFile from jgi_summarize_bam_contig_depths.

  -x [ --minCV ] arg (=1)           Minimum mean coverage of a contig in each library for binning.

  --minCVSum arg (=1)               Minimum total effective mean coverage of a contig (sum of depth over 

                                    minCV) for binning.

  -s [ --minClsSize ] arg (=200000) Minimum size of a bin as the output.

  -t [ --numThreads ] arg (=0)      Number of threads to use (0: use all cores).

  -l [ --onlyLabel ]                Output only sequence labels as a list in a column without sequences.

  --saveCls                         Save cluster memberships as a matrix format

  --unbinned                        Generate [outFile].unbinned.fa file for unbinned contigs

  --noBinOut                        No bin output. Usually combined with --saveCls to check only contig 

                                    memberships

  --seed arg (=0)                   For exact reproducibility. (0: use random seed)

  -d [ --debug ]                    Debug output

  -v [ --verbose ]                  Verbose output

[Error!] There was no --inFile specified

[Error!] There was no --outFile specified

> metabat1

$ metabat1

MetaBAT: Metagenome Binning based on Abundance and Tetranucleotide frequency (version 0.32.5; Aug 31 2017 21:02:53)

by Don Kang (ddkang@lbl.gov), Jeff Froula, Rob Egan, and Zhong Wang (zhongwang@lbl.gov) 

Allowed options:

  -h [ --help ]                     produce help message

  -i [ --inFile ] arg               Contigs in (gzipped) fasta file format [Mandatory]

  -o [ --outFile ] arg              Base file name for each bin. The default output is fasta format. Use -l 

                                    option to output only contig names [Mandatory]

  -a [ --abdFile ] arg              A file having mean and variance of base coverage depth (tab delimited; 

                                    the first column should be contig names, and the first row will be 

                                    considered as the header and be skipped) [Optional]

  --cvExt                           When a coverage file without variance (from third party tools) is used 

                                    instead of abdFile from jgi_summarize_bam_contig_depths

  -p [ --pairFile ] arg             A file having paired reads mapping information. Use it to increase 

                                    sensitivity. (tab delimited; should have 3 columns of contig index 

                                    (ordered by), its mate contig index, and supporting mean read coverage. 

                                    The first row will be considered as the header and be skipped) [Optional]

  --p1 arg (=0)                     Probability cutoff for bin seeding. It mainly controls the number of 

                                    potential bins and their specificity. The higher, the more (specific) 

                                    bins would be. (Percentage; Should be between 0 and 100)

  --p2 arg (=0)                     Probability cutoff for secondary neighbors. It supports p1 and better be 

                                    close to p1. (Percentage; Should be between 0 and 100)

  --minProb arg (=0)                Minimum probability for binning consideration. It controls sensitivity. 

                                    Usually it should be >= 75. (Percentage; Should be between 0 and 100)

  --minBinned arg (=0)              Minimum proportion of already binned neighbors for one's membership 

                                    inference. It contorls specificity. Usually it would be <= 50 

                                    (Percentage; Should be between 0 and 100)

  --verysensitive                   For greater sensitivity, especially in a simple community. It is the 

                                    shortcut for --p1 90 --p2 85 --pB 20 --minProb 75 --minBinned 20 

                                    --minCorr 90

  --sensitive                       For better sensitivity [default]. It is the shortcut for --p1 90 --p2 90 

                                    --pB 20 --minProb 80 --minBinned 40 --minCorr 92

  --specific                        For better specificity. Different from --sensitive when using correlation

                                    binning or ensemble binning. It is the shortcut for --p1 90 --p2 90 --pB 

                                    30 --minProb 80 --minBinned 40 --minCorr 96

  --veryspecific                    For greater specificity. No correlation binning for short contig 

                                    recruiting. It is the shortcut for --p1 90 --p2 90 --pB 40 --minProb 80 

                                    --minBinned 40

  --superspecific                   For the best specificity. It is the shortcut for --p1 95 --p2 90 --pB 50 

                                    --minProb 80 --minBinned 20

  --minCorr arg (=0)                Minimum pearson correlation coefficient for binning missed contigs to 

                                    increase sensitivity (Helpful when there are many samples). Should be 

                                    very high (>=90) to reduce contamination. (Percentage; Should be between 

                                    0 and 100; 0 disables)

  --minSamples arg (=10)            Minimum number of sample sizes for considering correlation based 

                                    recruiting

  -x [ --minCV ] arg (=1)           Minimum mean coverage of a contig to consider for abundance distance 

                                    calculation in each library

  --minCVSum arg (=2)               Minimum total mean coverage of a contig (sum of all libraries) to 

                                    consider for abundance distance calculation

  -s [ --minClsSize ] arg (=200000) Minimum size of a bin to be considered as the output

  -m [ --minContig ] arg (=2500)    Minimum size of a contig to be considered for binning (should be >=1500; 

                                    ideally >=2500). If # of samples >= minSamples, small contigs (>=1000) 

                                    will be given a chance to be recruited to existing bins by default.

  --minContigByCorr arg (=1000)     Minimum size of a contig to be considered for recruiting by pearson 

                                    correlation coefficients (activated only if # of samples >= minSamples; 

                                    disabled when minContigByCorr > minContig)

  -t [ --numThreads ] arg (=0)      Number of threads to use (0: use all cores)

  --minShared arg (=50)             Percentage cutoff for merging fuzzy contigs

  --fuzzy                           Binning with fuzziness which assigns multiple memberships of a contig to 

                                    bins (activated only with --pairFile at the moment)

  -l [ --onlyLabel ]                Output only sequence labels as a list in a column without sequences

  -S [ --sumLowCV ]                 If set, then every sample that falls below the minCV will be used in an 

                                    aggregate sample

  -V [ --maxVarRatio ] arg (=0)     Ignore any contigs where variance / mean exceeds this ratio (0 disables)

  --saveTNF arg                     File to save (or load if exists) TNF matrix for each contig in input

  --saveDistance arg                File to save (or load if exists) distance graph at lowest probability 

                                    cutoff

  --saveCls                         Save cluster memberships as a matrix format

  --unbinned                        Generate [outFile].unbinned.fa file for unbinned contigs

  --noBinOut                        No bin output. Usually combined with --saveCls to check only contig 

                                    memberships

  -B [ --B ] arg (=20)              Number of bootstrapping for ensemble binning (Recommended to be >=20)

  --pB arg (=50)                    Proportion of shared membership in bootstrapping. Major control for 

                                    sensitivity/specificity. The higher, the specific. (Percentage; Should be

                                    between 0 and 100)

  --seed arg (=0)                   For reproducibility in ensemble binning, though it might produce slightly

                                    different results. (0: use random seed)

  --keep                            Keep the intermediate files for later usage

  -d [ --debug ]                    Debug output

  -v [ --verbose ]                  Verbose output

[Error!] There was no --inFile specified

[Error!] There was no --outFile specified

——

使い方

アセンブリしたfastaファイルと、そのfastaにリードをマッピングして作ったbamを指定する。

runMetaBat.sh assembly.fasta sample1.bam [sample2.bam ...]

assembly.fasta .metabat-binsディレクトリができ、その中にビニングされたfastaが出力される。また、カレントにcontigそれぞれのリードデプスを記したテキストファイルassembly.fasta.depth.txtが出力される。

デプスファイルを指定し、より詳細なパラメータを指定してランすることもできる。"-a"で先ほどのリードデプス出力ファイルを指定する。

metabat2 -i input.fasta -a depth.txt -o output -t 12 -m 1500 -v --unbinned

• -i    Contigs in (gzipped) fasta file format [Mandatory]
• -o   Base file name and path for each bin. The default output is fasta format.
  Use -l option to output only contig names [Mandatory]
• -t     (=0) Number of threads to use (0: use all cores).
• -m   (=2500) Minimum size of a contig for binning (should be >=1500).
• -v    Verbose output
•  --unbinned   Generate [outFile].unbinned.fa file for unbinned contigs
• --saveCls       Save cluster memberships as a matrix format
• --noBinOut    No bin output. Usually combined with --saveCls to check only contig memberships

binningされた各fastaファイルが出力される。"--unbinned"をつけると、binningされなかったcontigのfastaファイルも別途出力される。

既知25種のゲノムからなるmock metagenome communityのbamファイルダウンロードリンク。

http://portal.nersc.gov/dna/RD/Metagenome_RD/MetaBAT/Software/Mockup/

引用

MetaBAT, an efficient tool for accurately reconstructing single genomes from complex microbial communities
Kang DD, Froula J, Egan R, Wang Z

PeerJ. 2015 Aug 27;3:e1165.

参考ページ

Binning and Visualization

https://github.com/HRGV/Marmics_Metagenomics/wiki/Module-6-Binning-and-Visualization

関連

メタゲノムのカバレッジやGCを可視化する簡単なwebツール gbtlite

gbtliteは、kbseahさんが作ったメタゲノムのカバレッジやGCのplotを描画してグラフ出力できるwebツール 。

I’ve written up a simple browser-based visualization for rendering coverage-GC% plots, called gbtlite.

https://kbseah.wordpress.com/2016/12/14/visualize-metagenomes-in-a-web-browser/

Github

使い方

１、BBtoolsのbbmapを使い、De novoアセンブリして得たFASTA配列に、アセンブリに使ったリードをマッピングする（BBtools紹介）。

bbmap.sh ref=assembly.fa nodisk in1=pair1.fq in2=pair2.fq covstats=sample1.stats

出力（sample1.stats）

２、https://kbseah.github.io/gbtlite/にアクセスする。

３、ファイルを選択をクリックし、fasta配列をアップロードする。Draw graphsをクリック。横軸がGC、縦軸がlogカバレッジのグラフが表示される。

グラフは、Trackpadの上下で拡大縮小したり、ドラッグで移動したりできる。上のメニューからは縦軸のスケールやプロットのサイズなどを変更できる。

感想

プロットが1000を超えてくるとブラウザがハングアップすることがあるそうなので、真のメタゲノムデータに使うには厳しそうですが、簡易チェックには役に立ちそうです。

gbtools（紹介）

https://kbseah.wordpress.com/2015/12/07/new-paper-visualizing-metagenome-binning-with-r/

引用

https://github.com/kbseah/gbtlite

バクテリアのシーケンシングデータ分析ツール GenomePeek

　シーケンシングコストが低下するにつれて、バクテリアゲノムの配列が増加している。現在、NCBI（Benson et al、2009; Sayers et al、2009）、SEEDデータベース（Overbeek、Disz＆Stevens、2004）には約15,000種類の原核生物ゲノムがあり、約75,000種類のアセンブリされていないデータがSequence Read Archiveに保存されている。 NCBIには約35,000のメタゲノムがあり、MG-RASTからは約90,000のメタゲノムが入手可能である（Meyer et al、2008）。完全なゲノムシーケンシングは、単一の原核生物種についての詳細な知識を提供するが、メタゲノムシーケンシングは、微生物環境の概要を私たちに提供する（Dinsdale et al、2008）。主要な目標の1つは、存在する種の分類学的起源を同定することである（Belda-Ferre et al、2012; Mande、Mohammed＆Ghosh、2012; Carr、Shen-OrrおよびBorenstein、 2013; Silva et al、2014）。

　メタゲノムに存在する種を同定するための2つの典型的なアプローチがある。最も一般的な方法は、既知の分類学的系統のデータベースに対して相同性検索を使用する事である（Altschul et al、1997; Meyer et al、2008; Segata et al、2012）。対照的に、ensembleアプローチは、タンパク質ドメイン頻度またはk-mer構成（Meinicke、Achhauer＆Lingner、2011; Silva et al、2014）のような、全リードのsignatureデータを使用する。相同性に基づく方法は、原核生物ゲノムの非常に多様で変異性の性質のために、一般に塩基ではなくタンパク質レベルのアライメントを使用する。このアプローチの問題は、メタゲノムシーケンシングのリードが、タンパク質のオープンリーディングフレームと比較して比較的短い傾向があることである。タンパク質をコードする原核生物遺伝子の平均長は約750bpである（Brocchieri＆Karlin、2005）。現在のシーケンシング技術は、30〜700 bp（Buermans＆den Dunnen、2014）の平均長のリードを生成する。シーケンシングデータにはサンプリングの確率的性質のためにオープンリーディングフレームの断片が含まれる。得られた予測タンパク質は、部分タンパク質と非タンパク質コード領域の部分的誤翻訳の混合物である。タンパク質断片が短ければ短いほど、原核生物間の共有された相同性のため（Wommack、Bhavsar＆Ravel、2008）、その断片に対するトップヒットを同定することが難しくなる。原核生物の分類学的同定に使用される一番の遺伝子は、全原核生物に普遍的に存在し、高度に保存されている16S rRNAスモールサブユニット遺伝子である（Woese＆Fox、1977; Lane et al、1985）。 16Sに基づく分析の問題の1つは、ゲノム中の16S遺伝子のコピー数がゲノムあたり1〜15コピーまで異なるため、定量に用いる前にゲノムあたりのコピー数に基づいて存在量を正規化しなければならない（Angly et al、2014）。そのため、代替案として、recA、rpoB、groEL、sodA、gyrB、nifD、fusA、およびdnaJを用いることが提案されている（Holmes、Nevin＆Lovley、2004;Adékambi ＆Drancourt、2004; Ghebremedhin et al、2008; Weng et al、2009）。

　メタゲノムサンプルとは対照的に、ゲノムデータの分類は、通常はコンティグへのアセンブリ後に実行されるが、シーケンシングのために培養された生物の先験的な知識に基づいており、影響を受ける。コンティグへのリードのアセンブリにより、完全なORFが同定される。これらの完全な遺伝子を用いて、シーケンシングされた株の分類を同定することができる。 RAST（Aziz et al、2008; Overbeek et al、2013）のようなツールは、BLASTP検索に基づく累積類似性に基づいて類似生物のリストを提供する。最初の培養が純粋でない場合、ゲノム中に存在する種の同定、および実際にはアセンブリおよび下流のアノテーションが妨げられる。不完全な培養は微生物学的技術が乏しいことから生じることがあるが、環境生物の研究では、常に複数の生物を含むいくつかの分離株が見つかる。著者らは、これらの生物が密接な相互関係を形成し、それゆえそれらの共培養が継続していると考えている（M Doane＆EA Dinsdale、2014、未発表データ）。

　メタゲノムデータを分析するための現在の方法論に関連する制限および問題を克服し、ゲノムシーケンシングデータをプレスクリーニングするため、著者らはWebベースのツールGenomePeek（https://edwards.sdsu.edu/GenomePeekから入手可能）を開発した。 GenomePeekは、原核生物の区別に有用な高度に保存された遺伝子セットと相同なシーケンシングデータを全て見つけることによって、存在する原核生物種を迅速に同定する。これらのリードをコンティグに組み立て、完全なORFを使用して、組み立てられた遺伝子の系統発生を決める。 GenomePeekは現在、4つの原核生物遺伝子、すなわち16S、recA、rpoB、およびgroEL、および4つの真核生物遺伝子（18S、RAD51、HSP60、RPB2）を分析する。ゲノムシーケンシングデータをGenomePeekで分析する場合、培養の同一性および純度が測定され、メタゲノムシーケンシングデータが分析される場合、その環境での分類学的分布が測定される。

FAQ

http://edwards.sdsu.edu/GenomePeek/faq.php

使い方

GenomePeekにアクセスする。

http://edwards.sdsu.edu/GenomePeek/index.php

fastqをアップロードする。

Add Files => Start upload。uploadが終わったらSubmit Sequencesをクリック。

 （対応拡張子: fastq, fq, fasta, fa, fna, qual, q. Files may be gzipped）

アセンブリされ、taxonomy classificationが実行される。さらに16S, recA, rpoB, and groELの定量が行われる。

1GB以上のシーケンシングデータはサブサンプリングして分析する。

引用

GenomePeek-an online tool for prokaryotic genome and metagenome analysis.
McNair K, Edwards RA

PeerJ. 2015 Jun 16;3:e1025.

環状ゲノムのシーケンシングデータのマッピングを改善する CircularMapper

　Graphフォーマットを使えば環状のリファレンスゲノムを正しく表現できるが、プレーンのFASTA形式には環状のリファレンスゲノムを正しく表現する方法が整備されていない。そのため、環状ゲノムのシーケンシングデータを"線状の"リファレンスゲノムFASTAにマッピングすると、FASTAの末端付近に相当する領域を読んだリードのアライメントが難しくなる（*1）。この環状ゲノムを考慮しないマッピングは、特に古いマッパーで起こりやすい。

　末端領域のカバレッジが落ち込むと、アライメント後の分析ステップで複数の問題が生じ得る。１つ目に、リードがアライメントされないため末端付近に変異（SNV、indel、CNV）が起きていても検出されなくなる（false negative）。２つ目に、ペアエンドリードが5'末端と3'末端にジャンプしてマッピングされると、実際は連続した領域にも関わらず、大きな構造変化が起きているとSV callerが勘違いする（false positiveの問題）。３つ目に、アライメントされなかったリードがゲノムのセカンドベストな部位にアライメントされ、そこで偽のバリアントコールを引き起こす（false positiveの問題）。１つ目の問題は、一般に5'末端と3'末端に遺伝子のコード領域がまたがった状態でリファレンスが登録されることはないため、たとえ変異の見逃しがあっても、大きな問題にはなりにくい（目的による）。２つ目の問題も、例えば、"1-7191318  Translocation"というようなゲノム全域に渡るSVコールは目視チェックで除けるため問題にならない。３つ目はより大きな問題かもしれない。オルガネラゲノムはコピー数が変動し、クロモソームよりコピー数が非常に多い場合もある。アライメントされなかったリードがリファレンスの別の領域にアライメントされると、そこで高度にミスアライメントを引き起こす（FPコールのホットスポットになる）可能性がある（*3）。このようなアライメントの誤りに由来するリードは、オルガネラゲノムのカバレッジの深さから、バリアントのフィルタリングを行っても信頼性の高いコールとして残る可能性がある（strand biasが強かったり、カバレッジが異常ならフィルタリングできるかもしれない）。この問題が根深いと表現したのは、例えば第三者が解析したデータにだけそのようなFPコールがあると、それが特異的なコールとして抽出され、candidatesのリストに残ってしまうことがありえるためである。

  CircularMapperは環状ゲノムのシーケンシングデータのリファレンス末端のアライメントを改善できるツール。延長したリファレンスにマッピングし、その後リアライメントして補正することで環状ゲノムのマッピングを改善する。Documentsを読む限り、本ツールはbwa alnで検証されている。

マニュアル

https://circularmapper.readthedocs.io/en/latest/contents/userguide.html

インストール

mac os10.13 java8でテストした。

依存

• Java Runtime Environment 8 or later

本体　Github

#Anaconda環境ならcondaで導入できる。それ以外はコンパイル済み実行jarファイルをダウンロードするか、ビルドする。
conda install -c bioconda circularmapper 

> circulargenerator

$ circulargenerator 

usage: CircularGenerator

 -e,--elongation <ELONGATION>   the elongation factor [INT]

 -h,--help                      show this help page

 -i,--input <INPUT>             the input FastA File

 -s,--seq <SEQ>                 the names of the sequences that should to

                                be elongated

Missing required options: i, e, s

> realignsamfile 

$ realignsamfile 

usage: RealignSAMFile

 -e,--elongation <ELONGATION>                          the elongation

                                                       factor [INT]

 -f,--filterCircularReads <FILTER>                     filter the reads

                                                       that don't map to a

                                                       circular identifier

                                                       (true/false),

                                                       default false

 -h,--help                                             show this help page

 -i,--input <INPUT>                                    the input SAM/BAM

                                                       File

 -r,--reference <REFERENCE>                            the unmodified

                                                       reference genome

 -x,--filterNonCircularSequences <FILTERSEQUENCEIDS>   filter the sequence

                                                       identifiers that

                                                       are not circular

                                                       identifiers

                                                       (true/false),

                                                       default false

Missing required options: i, r, e

ラン

まず500-bp modifyしたリファレンスゲノムを作る。ターゲットはリファレンスのミトコンドリアゲノムとする（FASTAのヘッダーで指定する *2）。

circularGenerator -e 500 -i ref.fa -s "chrMT"

 modifyされたfasta (ref_500.fa) ができる。

 modifyされたfastaにマッピングする。

#indexをつける。
bwa index ref_500.fa

#bwa alnでmapping
bwa aln -t 8 ref_500.fa merged.fastq -n 0.04 -l 32 > output.sai
bwa samse -r "specify_read_group_shere" ref_500.fa output.sai merged.fastq > output_circmapper.sam

samファイルができる。 

リアライメントを実行してオリジナルsamにマッピングした状態を再現する。オリジナルfastaを指定する。 

realignSAMFile -e 500 -i output_circmapper.sam -r ref.fa -f true -x false

 bamができる。 

bamをソートする。

samtools sort -@ 8 output_circmapper_realigned.BAM -o output_circmapper_realigned.sorted.bam
samtools index output_circmapper_realigned.sorted.bam

IGVで末端付近を表示した。上がCircularMapperでリアライメントを実行したbam、下がbwa aln=> samtoolsで特別な作業なしに作ったbam。

引用

GitHub - apeltzer/CircularMapper: A method to improve mappings on circular genomes, using the BWA mapper

参考

*1

結果は同じように見えるが、線状ゲノムの末端がシーケンスができなくてマッピングされない事とは全く話が異なる。

*2

headerが不明ならsamtools faidxでヘッダー情報を抜き出すか、grep -n ">" 等を実行する。

*3

NGSのマッパーはシーケンスエラーも考慮に入れた高感度なマッパーなので、アルゴリズムにもよるが、少々ミスアライメントがあってもアンマップにはなりにくい。そのため、ベストマッピング部位がリファレンスから意図的に取り除かれると、セカンドベストな部位にアライメントされることがままある。Exomeキャプチャのシーケンシングデータでも全ゲノムをリファレンスにしてマッピングを行うのは、キャプチャ時に少量のexome外領域が入り込み、そのような領域由来のリードがexomeに誤ってマッピングされfalse positiveな情報が発生するのを防ぐためである（他の理由も考えられる）。

review article要約 genome assembly reconciliation toolsの比較

　真核生物ゲノムの大部分は、それらを組み立てるというアルゴリズム上の課題のために未完成である。 様々なアセンブリやスキャフォールディングツールが利用できるが、特定のゲノムサイズや複雑さにどのツールやパラメータを使用するかは必ずしも明らかではない。 したがって、異なるアセンブラとパラメータを使用して複数のアセンブリを作成し、公開用に最適なアセンブリを選択するのが一般的な方法である。 より魅力的なアプローチは、より高品質なコンセンサスアセンブリを作成する目的で、複数のアセンブリをマージすることである。

本文要約

Introduction

　現在市場に出回っているシーケンシングマシンの莫大なスループットにもかかわらず、アセンブリ問題は、大きなゲノムのリピート、不均一なシーケンシングカバレッジ、および（不均一な）シークエンシングエラーおよびキメラリードの存在のために、非常に困難なままである。 Pacific Biosciences [ref.1]やOxford Nanopore [ref.2]のような第3世代シーケンシング技術は、1塩基あたりのコストが非常に高く、シーケンシングエラー率ははるかに高い。

• 　相当数のde novoゲノムアセンブラが利用可能である。最も適切なアセンブラの選択は、アセンブリされるゲノムのサイズおよび複雑さ（リピート内容、倍数性など）、リード生成に使われるシーケンシングマシンのタイプ（例えば、Sanger、454、Illumina、PacBio 、Nanoporeなど）、ペアエンドまたはロングインサートのメイトペアなどの利用可能性によって変わってくる。
• 各アセンブラは、ゲノムのリピート、不均一なカバレッジ、シーケンシングエラー、キメラリードに対処するために、わずかに異なるヒューリスティックを実装している。最終的なアセンブリはめったにコンプリートにはならず、典型的な出力は、contigsと呼ばれる順番や向きの揃っていない連続した領域のセットになる。
• 使用するアセンブラの選択は簡単な作業ではない。Genome Assembly Gold-Standard Evaluation（GAGE）[ref.3]やAssemblathon[ref.4]などのアセンブリコンペティションが開催され、共通データセットを使い複数のアセンブラを評価している。そのような比較評価は一般的なガイドラインを提供することができるが、特定のゲノムおよび特定のデータセットに対して最良のアセンブリを生成するために使用するアセンブラおよびパラメータ設定を決定する体系的な方法はない。
• 結果として、いくつかの異なるアセンブラおよび/またはパラメータ（例えば、de Bruijnグラフのk-merサイズ）から複数のゲノムアセンブリを生成し、アセンブリ統計に基づいて最良のアセンブリを推測しようとする試みが一般的である。実際、最良のアセンブリの概念は明確に定義されていない。アセンブリエラーのない、ゲノム全体をカバーする完璧なアセンブリを得ることはできないので、コンティグおよびscaffoldの長さを最大化すること（より多くのミスアセンブリを導入する犠牲を払い）が重要かどうかを決定しなければならない。
• ドラフトアセンブリの品質評価は、統計的測定およびリファレンスゲノムへのアライメント（利用可能な場合）によって行われる。 N50は、アセンブリの連続性を評価するために広く使用されるメトリックであり、これはコンティグがアセンブリの少なくとも50％をカバーする最短コンティグの長さによって定義される。NG50は、メトリックがアセンブリサイズではなくゲノムサイズに関連する点を除きN50に似ている（wiki）。
• アセンブリの正しさは、ミスマッチ、indels、misjoinsなどのアセンブリのエラーを検出することによって測定される。
• Misjoinsは、ゲノム内ではるかに離れている座位が不適切にアセンブリに参加している、ミスアセンブリの中で最も望ましくないタイプのエラーである[ref.5]。Misjoinsには、逆位、リアレンジメント、転座がある。
• Assembly reconciliationアルゴリズムはコンプリートゲノムを目指し、さらに一歩を踏み出す。Assembly reconciliationツールは、複数のドラフトアセンブリの中から最適なアセンブリを推測しようとするのではなく、魅力的な選択肢を提供する。これらのツールは、2つ以上のドラフトアセンブリをマージして、より高品質のコンセンサスアセンブリを生成することを約束するものである。Assembly reconciliationアルゴリズムの主な目的は、アセンブリの連続性を向上させると同時に、アセンブリエラーの発生を回避することである。
• 本論文では、異なるクオリティ属性を持ついくつかのコンセンサスアセンブリの品質を測定することにより、Assembly reconciliationツールの最初の包括的評価を実施した。

Assembly reconciliation tools

• Assembly reconciliationの概念は、Ziminらによって最初に導入された（pubmed）。その論文では、著者らはもはやメンテナンスされていない（2007年に最後に更新された）RECONCILIATORと呼ばれるAssembly reconciliationツールを導入していた。この論文では、もはや維持されていないAssembly reconciliationツールは評価から除外した。また、GAM_NGSに取って代わり、GAMも除外した。
• eRGA ref.[7]、MAIA [ref.8]、RAGOUT [ref.9]、Minimus2 [ref.10]などの他のツールも、これらのツールが異なる問題に対処しているため、比較評価には含まなかった。Reference-guided assembly（eRGA、RAGOUT、およびMAIA）とハイブリッドアセンブリ（Minimus2）は、アセンブリの調整の問題に関連しているが、全く同じではない。前者はターゲットに非常に近縁なリファレンスを使用してゲノムの保存された領域をアセンブルし、デノボアセンブリの非保存部分の複雑さを軽減する。後者のハイブリッドアセンブリでは、ユーザは異なるシーケンシング技術（例えば、イルミナショートリードとPacBioロングリード）からのリードをアセンブリする。
• closely related な複数のリファレンスゲノムが利用可能である場合、MAIAはde novoアセンブリをマージする能力も有する。
• この論文では、CISA、GAA、GAM_NGS、GARM、Metassembler、MIX、ZORROの7つのassembly reconciliationツールをベンチマークした。論文表1には、7つのassembly reconciliationツールの主な目標と機能をまとめている。これらのアルゴリズムのいくつかは、compression–expansion （CE）統計を使用してアセンブリのcompression（誤った欠失による）またはアセンブリのexpansion（誤った挿入による）を検出する[ref.6]。
• CE統計値を得るために、ペアエンドリードは評価されるアセンブリにマッピングされる。 CE統計は、実際にマップされたペアの距離と予想されるインサートサイズを比較することによって計算される。

各ツールの特徴

• CISAの目的は、バクテリアのゲノムアセンブリをreconcileすることである[ref.12]。各入力のドラフトアセンブリのコンティグが与えられると、CISAは代表コンティグ（すなわち最も長いコンティグ）を選択し、代表的なコンティグを延長しようとする。またクエリのコンティグを代表拡張コンティグにアライメントしてミスアセンブリを検出する。複数ポジションにアライメントされるコンティグは誤ったアセンブリと考えられ、別の代表的なコンティグが選択される。この過程でアライメントされなかった部分のコンティグは分割される。最後に、得られたコンティグが繰り返しマージされる。[ref.13,14]の論文で３つの異なるアセンブラの出力をマージするために用いられた。
• GAAのユーザは、ターゲットアセンブリの品質が高くなると予想されるターゲットアセンブリとクエリアセンブリを指定しなければならない。 GAAの目的は、クエリアセンブリを使用してターゲットアセンブリのギャップを埋めることである。少なくとも2つのターゲットコンティグにアンカーされていないクエリコンティグは無視される。GAAPは、[ref.17]の論文でSOAPdenovoアセンブリをNewblerアセンブリとマージするために用いられ、[ref.18]の論文ではNewblerアセンブリをPCAPアセンブリとマージするために用いられた。[ref,20]の論文ではABySSアセンブリをCLCアセンブリとマージするために用いられた。
• GAM_NGSの入力は、各リードのアライメントとアセンブリである[ref.22]。GAM_NGSは最初に、リードの同じセットを共有する両方の入力アセンブリ（ブロック）の最大部分を識別する。 GAM_NGSは、各頂点がコンティグに対応し、（i）それらが異なるアセンブリに属し、（ii）それらが少なくとも1つのブロックを共有する場合、エッジは2つのコンティグを連結する重み付き無向グラフを構築する。このグラフから、GAM_NGSは、両方の入力アセンブリに対するブロックの一貫した順序付けと方向付けを計算する。次に、別の有向グラフ（アセンブリグラフと呼ばれる）を構築し、各頂点はブロックを表し、各エッジは同じコンティグに属する場合に2つのブロックを接続する。アセンブリグラフの競合を解決した後、GAM_NGSは、少なくとも1つのブロックを共有する2つのコンティグ間のセミグローバル整合を計算する。 2つのコンティグの同一性が少なくとも95％である場合、GAM_NGSは最高のCE統計値を持つアセンブリを選択することによってアセンブリをマージする。 GAM_NGSはVelvet-SC、SPAdes、IDBA-UDで作成された3つのアセンブリをマージするために[ref.23]で使用された。 [ref.27]では、NewblerアセンブリをSOAPdenovoアセンブリとマージするために使用された。
• GARM [ref,28]はアセンブリを非対称に操作するが、ユーザーはより良いアセンブリーを事前に知る必要はない。このツールは、さまざまなアセンブリ統計に基づいてリファレンスアセンブリを決定する。 GARMは、（nucmer [ref.29]を使用して）オーバーラップを検出するためにアセンブリを互いにアライメントさせ、（ii）互いの（ほぼ）完全に包含されるあいまいなオーバーラップおよびコンティグを除去し、（iii）レイアウトおよびコンセンサスのスコアを生成し、 （iv）コンティグをマージし、そして（v）オリジナルのscaffoldsの順序と向きにマージされたコンティグをオーダーする。GARMは[ref.30]でIDb-UDアセンブリとNewblerアセンブリをマージするために使用された。
• 2つの入力アセンブリのCE統計は、Metassembler [ref.31]（紹介）でも使用されている。まず、Metassemblerはnucmer [ref.29]を使用して2つの入力アセンブリをアライメントする。これらの境界線はブレークポイントと呼ばれる。ブレークポイント間の各領域に対して、CE統計に基づいて2つのアセンブリの1つが選択される。 Metassemblerを使用すると、ユーザーは3つ以上のアセンブリを入力できるが、progressive pairwise の方式でマージされる。 [ref.32]では、Metassemblerを使用して、ALLPATHS-LGアセンブリをIllumina Moleculo [ref.33]合成ロングリードに基づくアセンブリとマージした。
• MIX [ref.5]（紹介）は、入力グラフ内のコンティグ間のprefix–suffix 重複を表現するために様々なウエイトを付けられた拡張グラフと呼ばれる有向グラフを使用する。 MIXはコンティグをマージするために、拡張グラフ上の重複しない最大の独立した最長パスセットを決定する。いずれの経路にも含まれていないコンティグは重複について検査される。含まれているかほぼ含まれているコンティグは削除され、残りはアセンブリに追加される。 MIXはエラー訂正を実行せず、むしろ連続性を高めることに焦点を合わせている。 MIXは[ref.34]と[ref.35]で2つのアセンブリをマージし、[ref.36]で3つのアセンブリをマージするために使用された。
• ZORRO [ref.37]は、k-mer統計を用いて反復領域を同定し、マスキングすることから始まる。繰り返し領域がマスクされると、Minimus [ref.10]を使用して2つのアセンブリ間のオーバーラップが検出される。それから繰り返し領域のマスクを解除し、重複するコンティグをマージする。最後に、ZORROはBambus [ref.38]を使用してペアエンドを使いコンティグを順序付けと方向付けを行う。 ZORROは[ref.39]と[ref.40]で2つのアセンブリをマージするために使用された。

Assembly reconciliationツールが大規模で複雑な（真核生物の）ゲノムに特に有益であることを期待する一方、Assembly reconciliationツールは多くのバクテリアゲノムのアセンブリプロジェクトで活用されている。

結果

　上記の7つのアセンブリ調整ツールの比較評価を行うために、著者らはGAGE用に作成された公的に入手可能なアセンブリを使用した[ref.3]。また構造変化を含むSaccharomyces cerevisiae S288c [37]の合成アセンブリを作成した。GAGEアセンブリを選択した動機は、このデータセットがAssembly reconciliationツールで最も一般的に使用されていたことである。 GAM_NGSの著者は実験結果にこのデータセットを使用し、CISAはS. aureusとR. sphaeroidesのアセンブリでテストされ、MIXはS. aureusとR. sphaeroidesのアセンブリを含むGAGE_B [ref.41]のデータを使用している。他のAssembly reconciliationツールはAssemblathonデータセット[ref.4]を使用している。たとえば、MetassemblerはAssemblathon 1とAssemblathon 2の両方のデータセットを使用している。

• すべてのAssembly reconciliationツールはデフォルトのパラメータで実行され、Quast [ref.42]は広範なアセンブリ統計の収集に使用された（詳細は論文のメソッド参照）。すべての統計に関する完全なレポートは、追加ファイル1: 表S1-S19に報告されている。本文では連続性/正確性のトレードオフの結果のグラフのみをまとめる。
• 入力アセンブリおよび出力アセンブリは、x座標が連続性（NGA50）を表し、y座標がミスアセンンブリの数である散布図上の点として表される。図1はプロットの解釈方法を示している。Assembly reconciliationツールは、入力点をプロットの右下隅に向かって「移動」させる、すなわち、連続性を高め、アセンブリ誤差の数を減らすことが期待される（図1に例 link）。

• すべてのアセンブリのペアで各アセンブリのreconciliationツールを実行すると数百のアセンブリが生成され、一般的な結論を導き出すことは困難になる。代わりに、6つの異なる基準に基づいて入力アセンブリのペアを選択し、選択したペアの結果を比較することにした。プレゼンテーションを合理化するため、正常に実行されなかったツールについてはコメントしない。これらのツールの実行に関連するその他の制限事項は、「論文の追加ファイル1：注6」に記載されている。最後に、対応するアセンブリに加えて、raw シーケンシングリードを利用できるツールもある（「論文の追加ファイル1：注7」を参照）。

（結果はより簡潔に、影響の大きい結果についてだけまとめる）

連続性の高い配列と精度の高い配列のマージ（GAGEより）

　最初の一連の実験の目的は、contiguity（以後、連続性）と精度の間のトレードオフを探すことだった。具体的には、最初の入力アセンブリが連続性に優れ、２番目のアセンブリが精度に優れる場合、1番目の入力アセンブリの連続性と2番目の入力アセンブリの正確さでマージするreconciliationツールの能力をテストしたかった。マージされる2つの入力アセンブリは、片方は高いN50値（ただし、誤ったアセンブリエラーを持つ可能性が高い）と誤ったアセンブリが少ない（おそらくN50が低い）からなるように選択された。

• いずれのassembly reconciliationツールもABySSアセンブリに比べてアセンブリエラーを改善できなかった。
• CISAはミスアセンブリの数が最も少なかった。
• 入力アセンブリがcontigではなくscaffoldsで構成されている場合（図2の下段）、すべてのassembly reconciliationツールは連続性をわずかに改善した（5％未満）。
• 予想通り、マスターアセンブリに依存するツールは、入力をランク付けしなかったツールよりも誤ったアセンブリ数が少なかった。

入力アセンブリの並べ替え（GAGEより）

入力の順番による影響。（省略）

ハイクオリティな入力（GAGEより）

　3番目の実験では、2つの高品質アセンブリを結合するassembly reconciliationツールの機能をテストした。 2つのハイクオリティアセンブリ（すなわち、少数のコンティグ / scaffoldsおよび高いN50と、よりミスアセンブリ数が少ない、の２つを選択した。ALLPATHS-LGを最初の入力とし、MSR-CA、SOAPdenovo、CABOGのいずれかを2番目の入力としてマージした。

Table S7より（link）

S. aureus

• S. aureusではcontigを入力とするとGAM_NGSがベストで、ミスアセンブリなく連続性を66％改善させた。２番目に最良のアセンブリは107％の連続性増加を伴うMetassemblerによるものだったが、ALLPATHS-LGと比較してミスアセンブリ数はわずかに増加した。 MIXは多数のミスアセンブリを生成したが（MSR-CAより高い）、連続性を4％向上させた。 CISAは連続性を11％改善したが、ALLPATHS-LGよりも多くのエラーを生成した。 ZORROは連続性を30％低下させた。
• Scaffoldsを入力とした場合、そもそもALLPATHS-LGはミスアセンブリがなく、NGA50もMSR-CAより高かった。asymmetricなツールとsymmetricなツールがあるが、傾向として、asymmetricなツールはよりミスアセンブリが少なく、N50を減少させた。 ZORROは非対称であるが、連続性は90％以上減少した。
•  GARMは、MSR-CAに近い多数のミスアセンブリを生じたが、連続性を最も高めた（108％）。 MIXはミスアセンブリを導入しなかったが、ゲノム配列の25％しかカバーしなかった。
• GAM_NGSとMetassemblerは連続性を66.5％改善し、ミスアセンブリを導入しなかった。これらはアセンブリの驚異的な改善が見られた2つのまれな例である。

R. sphaeroides

• MetassemblerだけがNGA50を大幅に増やした。他のすべてのツールは、連続性を減少させた。正確さに関しては、ZORROとCISA（Scaffoldsを入力として使用）はミスアセンブリ数を減らしたが、連続性もそれぞれ99％と60％低下させた。
• GARMは、連続性を38％改善し、CISAは2％未満向上させた。ミスアセンブリ数を減らした唯一のツールはMIXだったが、そのアセンブリでもゲノムの約半分しかカバーしていなかった。
• いずれのツールも、連続性とミスアセンブリの両方は改善しなかった。

Hg_chr14

• GAAはNGA50を76％改善したが、入力のミスアセンブリ合計に等しい数のミスアセンブリを生成した。 GAM_NGSは隣接性を改善し（NGA50は28％の増加）、ミスアセンブリ数をわずかに減らした。
• Scaffoldsをを入力として使用すると、GAM_NGSとMetassemblerはALLPATHS-LGと同様の品質統計を維持した。 GARMはNGA50を9％減少させた。また、ミスアッセンブリの数も増加した。

高度に断片化された入力（GAGEより）

　この一連の実験の目的は、2つの高度に断片化された入力アセンブリが提供されている場合、アセンブリ調整ツールのパフォーマンスを評価することだった。入力アセンブリは、200bpsより短いコンティグ、多いコンティグ数と低いN50を有するように選択された。（結果省略）

De Bruijn対 string graph アセンブリ（GAGEより）

　異なるアセンブリ方法を使用して生成されたアセンブリをマージする効果をテストした。ALLPATHS-LG（de Bruijn graphに基づく）によって生成されたアセンブリを、SGAによって生成されたアセンブリ（string graphに基づく）とマージした結果を示す。全体として、GAM_NGS、Metassembler、およびMIXはALLPATHS-LGと同様のアセンブリ統計を維持した。

• Hg_chr14のコンティグを入力の場合、GAM_NGSはNGA50を2％増加させた。GAM_NGSとMetassemblerはALLPATHS-LGに近いミスアセンブリと連続性を維持した。 GARMは、misassembliesの数をALLPATHS-LGの455から496に増やし、NGA50を9％減少させた。

複数入力（GAGEより）

　3つ以上のアセンブリをマージするツールの能力をテストした。アセンブリ調整ツールが2つ以上のアセンブリを入力として許可しなかった場合、それらを反復的にマージした。たとえば、3つのアセンブリをマージする時は、最初に2つのアセンブリをマージし、その結果を3番目のアセンブリとマージした。 Metassemblerは、同様の戦略を使用している。ユーザーが複数のアセンブリを提供する場合、ツールは反復的にペアワイズ調整を実行する。 [ref.51]で提案されたアセンブリ品質メトリックであるフィーチャレスポンスカーブ（FRカーブ）に基づいて、入力アセンブリの順番を決めた。 FRカーブは、τの特徴を含まないコンティグと対応するゲノムカバレッジとの間の依存性を表す。 x軸はτを表し、y軸はゲノムカバレッジを表す。カーブが急であればあるほど、アセンブリは良好である。 [ref,22]のFRカーブを使用して、GAGEアセンブリのマージ順序を決定し、最高品質のアセンブリから始めた。

• S. aureusとR. sphaeroidesについては、CISAは一般に、マージしたアセンブリ数が増加するにつれて連続性を改善した。繰り返しによりエラー数は変動した。GAM_NGSの連続性も繰り返しマージで改善されたが、ミスアセンブリ数は減少しなかった。
• Metassemblerの連続性は、S. aureusについては繰り返しマージで改善されたが、ミスアセンブリ数も増加した。
• MIXはほとんどの繰り返しでミス・アセンブリはで少数を維持したが、NGA50は比較的貧弱だった。

合成アセンブリ

　この一連の実験では、特定の構造変異（詳細については方法を参照）を埋め込んだSaccharomyces cerevisiaeの合成アセンブリのアセンブリ調整ツールをテストした。 Decipher [ref.54]は、gradientsとして表示されるシンテシープロット生成に使用した。参照とクエリが一致しない場合、gradientsのカラーは中断される（Fig.6）。灰色の領域は、参照と一致しないブロックを示す。各実験では、（1）酵母ゲノムの4番染色体と15番染色体の2つのマージ、（2）1のバージョンにRSVSimによって生成された1つの構造変異、すなわち欠失、逆位 、転座を加えたもの。様々なサイズ（50,100,200および500kbp）の欠失および逆位を第4染色体に導入し、第4染色体から第15染色体へ様々なサイズ（再び50,100,200および500kbp）の転座を生成した。

• 図6（一番上の行）は、CISAが欠失を解決したが、第15染色体を出力しなかったことを示している。GARMは染色体4を出力しなかった。GAM_NGS、Metassembler、およびMIXは、欠陥のある入力アセンブリを生成した。ZORROは、欠失の位置でアセンブリを壊し、3つの個々のコンティグを生成し、欠失した配列を省略した。
• 図6（中段）は、CISAのみが逆位を解決したが、第15染色体は出力しなかったことを示している。GAAは逆位を修正しなかった 。再びGAM_NGS、Metassembler、MIXは欠陥のある入力アセンブリを生成した。 ZORROは、4番染色体の3つのコンティグと、15番染色体を表すコンティグを追加して逆位を起こした。
• 図6（一番下の行）に示すように、転座ツールの動作は転座のサイズに依存していた。 50,100、および200kbpの転座について、CISA、GAA、およびGAM_NGSは、第4染色体の正しいバージョンを産生した。 GARMはマージされたアセンブリを生成した。 MetassemblerとMIXの出力は、欠陥のある入力アセンブリのようだった。 ZORROはアセンブリを構造変化の点で分割する。 200および500 kbpの場合、350 GBを超えるRAMを割り当てた後、ZORROを停止した。いずれのツールも500kbpの転座を修正することはできなかった。

討論と結論

• 　デノボアセンブリの実際的な課題を考えると、アセンブリの調整のアイデアは非常に魅力的である。さまざまなアセンブリツールおよび/またはパラメータを使用して、同じデータセット上に複数のアセンブリを生成し、次にassembly reconciliation ツールを使用してすべてのアセンブリをマージし、高品質のコンセンサスアセンブリを得ることができる。
• 統合されたアセンブリの品質は、少なくとも入力における最適なアセンブリの品質と同じである必要がある。事実、両方の入力アセンブリが良好な品質のアセンブリ統計値を有する場合（例えば、一方がより連続し、他方がよりエラーが少ない場合）、コンセンサスアセンブリは両方の入力から良い面を継承することが期待される。
• 現実には、両方の入力アセンブリよりも一貫して優れている（または少なくとも少なくとも優れている）アセンブリを作成することは非常に難しい。この原稿で報告された一連の実験では、著者らが評価したツールのどれもがこの目標を一貫して達成することができなかった。出力が両方の入力より優れているケースは非常に少なかった。
• これらのアセンブリツールが一般的なアセンブリの調整問題を解決できなかったが、各ツールはアルゴリズムの進歩をもたらす可能性のあるいくつかの長所を示した。
• 例えば、CISAは一般的に、ほとんどの構造的なバリエーションを修正し、入力アセンブリの重複を無視できた（ただしマージしたアセンブリの数が増えるにつれてduplicationの割合が増加した）。
• GAAとGARMは、しばしば連続性を改善した（しかし、しばしばアセンブリエラーをもたらしduplicationの割合を増加した）。
• GAM_NGSは、通常、参照の質に非常に近いコンセンサスアセンブリを作成したが（それほど良くはない）、転座を解決することができた。
• MIXは連続性を適度に改善したが（しかしアセンブリエラー数は、最もエラーの多い入力に近いかそれ以上だった。場合によってはゲノムカバー率が低下した）
• Metassemblerは、しばしばミスあアセンブリが非常に少なく、両方の入力のアセンブリよりミスアセンブリが少ない場合もあった（しかし、N50は増加しなかった）。
• 最後に、ZORROは一般的に高いゲノムカバー率を維持した（しかし、それは連続性を有意に増加させなかった）。

引用
A comparative evaluation of genome assembly reconciliation tools
Hind Alhakami,corresponding author Hamid Mirebrahim, Stefano Lonardi

Genome Biol. 2017; 18: 93.

参考HP

http://compgenomics2016.biology.gatech.edu/index.php/Genome_Assembly_Group