macでインフォマティクス

macでインフォマティクス

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

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

2019 5/27 追記

 

メタゲノミクスは、環境サンプルから直接採取した遺伝物質を研究する。 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)

*NOTICE: Genes are not predicted by default. Use --gene-finding or --glimmer option to enable it.
 

リファレンスを指定しない場合、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

 

関連