macでインフォマティクス

macでインフォマティクス

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

taxonomy ID、学名、系統情報など相互変換できるツール taxonkit

 

NCBI taxonomy databaseは、公共のシーケンスデータベースに含まれるすべての生物の分類(以後taxonomy)および命名法をまとめたものである(NCBI、2018)。taxonomyデータベースの一般的な操作には、分類名からのtaxonomy ID(taxid)問い合わせ、または分類群からの系統問い合わせなどがある。NCBI taxaonomyデータベース(https://www.ncbi.nlm.nih.gov/taxonomy)には、視覚化とデータクエリのための包括的なドキュメントとツールが用意されている。
 NCBIの公式Entrezプログラミングユーティリティ(E-utilities)(Sayers、2010)は、taxonomyデータベースから問い合わせるためのWebサービスを提供している。これは、バッチ処理などのプログラミング的な方法で使用できる。いくつかのツールキットもクエリ機能を提供している。 BioPythonのEntrezモジュール(Cock、et al、2009)は、データベースに対するさまざまな問い合わせのために、オンラインのEntrezサービスをpython関数にラップする。これらのオンラインクエリツールは最新のデータへのアクセスを提供するが、その間、処理速度はインターネット接続によって制限される。スタンドアロンツールはクエリを高速化するためにtaxonomy dump file(ftp://ftp.ncbi.nih.gov/pub/taxonomy/)をローカルデータベースにパースする。 ETE Toolkit(Huerta-Cepas et al、2016)は、NCBI taxonomy dump fileをSQLiteデータベースにして、pythonモジュールを介してクエリ機能を提供する。 Taxadb(Gourlé、2006)は、taxonomyデータをSQLiteMySQLPostgreSQLなどのリレーショナルデータベースに変換する。 Ncbitax2lin(Xue、2016)は、NCBI taxonomy dump fileを系統に変換するPythonコマンドラインスクリプトである。それでも、taxonomyデータベースは頻繁に更新されるため、ローカルデータベースの再構築が必要になる。リレーショナルデータベースからのSQLクエリは、バッチクエリでは特に効率的ではない。taxonomy dump fileのファイルサイズがわずかであることと、ストレージデバイスの読み取り速度が徐々に増加していることを考慮すると、一般的なtaxonomyのデータ操作には、インターネット接続の制限がないダンプファイルの直接解析が速い。この記事では、NCBI taxonomyデータを迅速に操作するためのコマンドライン・ツールキットTaxonKitについて説明する。

 taxonomyダンプファイルを効率的に解析するために、ツールキットを実装するためのコンパイル言語Go(https://golang.org/)を選択した。これは、WindowsLinux、およびMac OS Xを含む最も一般的なオペレーティングシステムのクロスコンパイルをサポートしている。 単一の実行可能バイナリファイルは、展開と実行が簡単である。 Taxonkitは「コマンドサブコマンド」の構造を採用している。サブコマンドは完全に独立した機能を提供し、標準ストリームまたはローカルファイルからのプレーンまたはgzip圧縮の入出力をサポートしている。 したがって、Taxonkitをコマンドラインパイプに簡単に組み合わせて、複雑な操作を実行できる。
 現時点では、names.dmpおよびnodes.dmpファイル内の命名法および階層情報のみが使用されている。 Taxonkitのすべてのサブコマンドは、必要なデータをランダムアクセスメモリ(RAM)にロードし直接解析する。デーモンプロセスは使用しない。

 

usage

https://bioinf.shenwei.me/taxonkit/usage/

tutorial

Tutorial - TaxonKit - NCBI Taxonomy Toolkit

 

taxonkitに関するツイート


インストール

mac os10.14向け64bitバイナリをダウンロードし、テストした。

TaxonKit is implemented in Go programming language

本体

Donwloadよりmaclinuxwindows向けバイナリをダウンロードできる。

sudo cp taxonkit /usr/local/bin/

#Anacondaを使っているならcondaでも導入可能
conda install -c bioconda taxonkit

#他にもbrew、goでのインストールがサポートされている。上のdownloadリンク参照

> taxonkit

$ taxonkit 

TaxonKit - A cross-platform and Efficient NCBI Taxonomy Toolkit

 

Version: 0.3.0

 

Author: Wei Shen <shenwei356@gmail.com>

 

Source code: https://github.com/shenwei356/taxonkit

Documents  : http://bioinf.shenwei.me/taxonkit

 

Dataset:

 

    Please download and decompress "taxdump.tar.gz":

    ftp://ftp.ncbi.nih.gov/pub/taxonomy/taxdump.tar.gz

 

    and copy "names.dmp" and "nodes.dmp" to data directory:

    "/Users/kazuma/.taxonkit"

 

    or some other directory, and later you can refer to using flag --data-dir,

    or environment variable TAXONKIT_DB

 

Usage:

  taxonkit [command]

 

Available Commands:

  genautocomplete generate shell autocompletion script

  help            Help about any command

  lineage         query lineage of given taxids

  list            list taxon tree of given taxids

  name2taxid      query taxid by taxon scientific name

  reformat        reformat lineage

  version         print version information and check for update

 

Flags:

      --data-dir string   directory containing nodes.dmp and names.dmp (default "/Users/kazuma/.taxonkit")

  -h, --help              help for taxonkit

      --line-buffered     use line buffering on output, i.e., immediately writing to stdin/file for every line of output

  -o, --out-file string   out file ("-" for stdout, suffix .gz for gzipped out) (default "-")

  -j, --threads int       number of CPUs. (default value: 1 for single-CPU PC, 2 for others) (default 2)

      --verbose           print verbose information

 

Use "taxonkit [command] --help" for more information about a command.

 

データベースの準備

使用前にNCBIの系統分類データベース(学名、階級など)を準備する。

wget ftp://ftp.ncbi.nih.gov/pub/taxonomy/taxdump.tar.gz
tar -zxvf taxdump.tar.gz

mkdir $HOME/.taxonkit
cp names.dmp nodes.dmp $HOME/.taxonkit/

 最新のデータが使いたいなら、定期的にダウンロードし直す。

 

実行方法

1、taxonkit list   taxonomy IDの照会

$ taxonkit list -h

list taxon tree of given taxids

 

Usage:

  taxonkit list [flags]

 

Flags:

  -h, --help            help for list

      --ids string      taxid(s), multiple values should be separated by comma

      --indent string   indent (default "  ")

      --json            output in JSON format. you can save the result in file with suffix ".json" and open with modern text editor

      --show-name       output scientific name

      --show-rank       output rank

 

Global Flags:

      --data-dir string   directory containing nodes.dmp and names.dmp (default "/Users/kazuma/.taxonkit")

      --line-buffered     use line buffering on output, i.e., immediately writing to stdin/file for every line of output

  -o, --out-file string   out file ("-" for stdout, suffix .gz for gzipped out) (default "-")

  -j, --threads int       number of CPUs. (default value: 1 for single-CPU PC, 2 for others) (default 2)

      --verbose           print verbose information

NCBI taxonomy ID照会サービス(link)と同じように使える。-

 

指定のtaxon IDリストを表示

taxonkit list -j 8 --ids 9605,239934

$ taxonkit list --ids 9605,239934

9605

  9606

    63221

    741158

  1425170

 

239934

  239935

    349741

  512293

    512294

    1131822

    1262691

    1263034

  1131336

  1574264

  1574265

  1638783

  1679444

  1755639

  1872421

  1896967

  1929996

  1945963

  1945964

  1945965

  1945966

  2478952

分類学の階級、 学名も表示する。

taxonkit list -j 8 --show-rank --show-name --ids 9605
  • --show-name   output scientific name
  • --show-rank   output rank
  • -j     number of CPUs. (default value: 1 for single-CPU PC, 2 for others) (default 2)

 

$ taxonkit list --show-rank --show-name --ids 9605,239934

9605 [genus] Homo

  9606 [species] Homo sapiens

    63221 [subspecies] Homo sapiens neanderthalensis

    741158 [subspecies] Homo sapiens subsp. 'Denisova'

  1425170 [species] Homo heidelbergensis

 

"--json"をつければ、JSON形式で出力できる。目的によってはかなり役立つオプション。

  

 

2、lineage    -  系統情報の照会

$ taxonkit lineage -h

query lineage of given taxids

 

Usage:

  taxonkit lineage [flags]

 

Flags:

  -d, --delimiter string      field delimiter in lineage (default ";")

  -h, --help                  help for lineage

  -t, --show-lineage-taxids   show lineage consisting of taxids

  -r, --show-rank             show rank of taxids

  -i, --taxid-field int       field index of taxid. data should be tab-separated (default 1)

 

Global Flags:

      --data-dir string   directory containing nodes.dmp and names.dmp (default "/Users/kazuma/.taxonkit")

      --line-buffered     use line buffering on output, i.e., immediately writing to stdin/file for every line of output

  -o, --out-file string   out file ("-" for stdout, suffix .gz for gzipped out) (default "-")

  -j, --threads int       number of CPUs. (default value: 1 for single-CPU PC, 2 for others) (default 2)

      --verbose           print verbose information

例えばファイル中のtaxonomy IDから系統を照会。

$ cat taxids.txt
9606
349741
239935

$ taxonkit lineage taxids.txt

$ taxonkit lineage taxids.txt

9606 cellular organisms;Eukaryota;Opisthokonta;Metazoa;Eumetazoa;Bilateria;Deuterostomia;Chordata;Craniata;Vertebrata;Gnathostomata;Teleostomi;Euteleostomi;Sarcopterygii;Dipnotetrapodomorpha;Tetrapoda;Amniota;Mammalia;Theria;Eutheria;Boreoeutheria;Euarchontoglires;Primates;Haplorrhini;Simiiformes;Catarrhini;Hominoidea;Hominidae;Homininae;Homo;Homo sapiens

349741 cellular organisms;Bacteria;PVC group;Verrucomicrobia;Verrucomicrobiae;Verrucomicrobiales;Akkermansiaceae;Akkermansia;Akkermansia muciniphila;Akkermansia muciniphila ATCC BAA-835

ランクを表示。

taxonkit lineage -j 8 -r taxids.txt | cut -f 1,3
  • -r    show rank of taxids

$ taxonkit lineage -r taxids.txt | cut -f 1,3

9606 species

349741 no rank

239935 species

系統と、照会元のtaxonomic IDも出力する。

taxonkit lineage -j 8 -t taxids.txt
  • -t    show lineage consisting of taxids

9606 cellular organisms;Eukaryota;Opisthokonta;Metazoa;Eumetazoa;Bilateria;Deuterostomia;Chordata;Craniata;Vertebrata;Gnathostomata;Teleostomi;Euteleostomi;Sarcopterygii;Dipnotetrapodomorpha;Tetrapoda;Amniota;Mammalia;Theria;Eutheria;Boreoeutheria;Euarchontoglires;Primates;Haplorrhini;Simiiformes;Catarrhini;Hominoidea;Hominidae;Homininae;Homo;Homo sapiens 131567;2759;33154;33208;6072;33213;33511;7711;89593;7742;7776;117570;117571;8287;1338369;32523;32524;40674;32525;9347;1437010;314146;9443;376913;314293;9526;314295;9604;207598;9605;9606

349741 cellular organisms;Bacteria;PVC group;Verrucomicrobia;Verrucomicrobiae;Verrucomicrobiales;Akkermansiaceae;Akkermansia;Akkermansia muciniphila;Akkermansia muciniphila ATCC BAA-835 131567;2;1783257;74201;203494;48461;1647988;239934;239935;349741

239935 cellular organisms;Bacteria;PVC group;Verrucomicrobia;Verrucomicrobiae;Verrucomicrobiales;Akkermansiaceae;Akkermansia;Akkermansia muciniphila 131567;2;1783257;74201;203494;48461;1647988;239934;239935

 

 

3、reformat    2の"lineage"コマンドの系統出力ファイルをリフォーマット

$ taxonkit reformat -h

reformat lineage

 

Output format can be formated by flag --format, available placeholders:

 

    {k}: superkingdom

    {p}: phylum

    {c}: class

    {o}: order

    {f}: family

    {g}: genus

    {s}: species

    {S}: subspecies

 

Output format can contains some escape charactors like " ".

 

This command appends reformated lineage to the input line.

The corresponding taxids of reformated lineage can be provided as another

column by flag "-t/--show-lineage-taxids".

 

Usage:

  taxonkit reformat [flags]

 

Flags:

  -d, --delimiter string               field delimiter in input lineage (default ";")

  -F, --fill-miss-rank                 fill missing rank with original lineage information (experimental)

  -f, --format string                  output format, placeholders of rank are needed (default "{k};{p};{c};{o};{f};{g};{s}")

  -h, --help                           help for reformat

  -i, --lineage-field int              field index of lineage. data should be tab-separated (default 2)

  -r, --miss-rank-repl string          replacement string for missing rank, if given "", "unclassified xxx xxx" will used, where "unclassified " is settable by flag -p/--miss-rank-repl-prefix

  -p, --miss-rank-repl-prefix string   prefix for estimated taxon level (default "unclassified ")

  -R, --miss-taxid-repl string         replacement string for missing taxid

  -t, --show-lineage-taxids            show corresponding taxids of reformated lineage

 

Global Flags:

      --data-dir string   directory containing nodes.dmp and names.dmp (default "/Users/kazuma/.taxonkit")

      --line-buffered     use line buffering on output, i.e., immediately writing to stdin/file for every line of output

  -o, --out-file string   out file ("-" for stdout, suffix .gz for gzipped out) (default "-")

  -j, --threads int       number of CPUs. (default value: 1 for single-CPU PC, 2 for others) (default 2)

      --verbose           print verbose information

 

2の出力フォーマットを再変換するために使う。

#まずlineageコマンドを実行、idから系統情報を得る
taxonkit lineage taxids.txt > lineage.txt

#情報が多くて分かりにくい。reformatコマンド出力で簡潔にする
taxonkit reformat -j 8 lineage.txt | cut -f 3

$ taxonkit reformat lineage.txt | cut -f 3

Eukaryota;Chordata;Mammalia;Primates;Hominidae;Homo;Homo sapiens

Bacteria;Verrucomicrobia;Verrucomicrobiae;Verrucomicrobiales;Akkermansiaceae;Akkermansia;Akkermansia muciniphila

Bacteria;Verrucomicrobia;Verrucomicrobiae;Verrucomicrobiales;Akkermansiaceae;Akkermansia;Akkermansia muciniphila

genusとspeciesのみ表示する

taxonkit reformat -j 8 -f "{g};{s}" lineage.txt | cut -f 3
  • -f     string output format, placeholders of rank are needed (default "{k};{p};{c};{o};{f};{g};{s}") 

$ taxonkit reformat -j 8 -f "{g};{s}" lineage.txt | cut -f 3

Homo;Homo sapiens

Akkermansia;Akkermansia muciniphila

Akkermansia;Akkermansia muciniphila

 

 

 

4、 name2taxid     - 学名からtaxonomic idを照合

$ taxonkit name2taxid -h

query taxid by taxon scientific name

 

Usage:

  taxonkit name2taxid [flags]

 

Flags:

  -h, --help             help for name2taxid

  -i, --name-field int   field index of name. data should be tab-separated (default 1)

      --show-rank        show rank

 

Global Flags:

      --data-dir string   directory containing nodes.dmp and names.dmp (default "/Users/kazuma/.taxonkit")

      --line-buffered     use line buffering on output, i.e., immediately writing to stdin/file for every line of output

  -o, --out-file string   out file ("-" for stdout, suffix .gz for gzipped out) (default "-")

  -j, --threads int       number of CPUs. (default value: 1 for single-CPU PC, 2 for others) (default 2)

      --verbose           print verbose information

学名のファイルを準備し、実行

$ cat names.txt
Homo sapiens
Akkermansia muciniphila ATCC BAA-835
Akkermansia muciniphila
Mouse Intracisternal A-particle
Wei Shen
uncultured murine large bowel bacterium BAC 54B
Croceibacter phage P2559Y

$ cat names.txt | taxonkit name2taxid

$ cat names.txt | taxonkit name2taxid

Homo sapiens 9606

Akkermansia muciniphila ATCC BAA-835 349741

Akkermansia muciniphila 239935

Mouse Intracisternal A-particle 11932

Wei Shen

uncultured murine large bowel bacterium BAC 54B 314101

Croceibacter phage P2559Y 1327037

taxonomic idが右端につく。

 

分類階級も表示 

cat names.txt | taxonkit name2taxid --show-rank

$ cat names.txt | taxonkit name2taxid --show-rank

Homo sapiens 9606 species

Akkermansia muciniphila ATCC BAA-835 349741 no rank

Akkermansia muciniphila 239935 species

Mouse Intracisternal A-particle 11932 species

Wei Shen

uncultured murine large bowel bacterium BAC 54B 314101 species

Croceibacter phage P2559Y 1327037 species

 

学名 => taxonomic id => 分類系統ファイル

cat names.txt | taxonkit name2taxid | taxonkit lineage --taxid-field 2

$ cat names.txt | taxonkit name2taxid | taxonkit lineage --taxid-field 2

Homo sapiens 9606 cellular organisms;Eukaryota;Opisthokonta;Metazoa;Eumetazoa;Bilateria;Deuterostomia;Chordata;Craniata;Vertebrata;Gnathostomata;Teleostomi;Euteleostomi;Sarcopterygii;Dipnotetrapodomorpha;Tetrapoda;Amniota;Mammalia;Theria;Eutheria;Boreoeutheria;Euarchontoglires;Primates;Haplorrhini;Simiiformes;Catarrhini;Hominoidea;Hominidae;Homininae;Homo;Homo sapiens

Akkermansia muciniphila ATCC BAA-835 349741 cellular organisms;Bacteria;PVC group;Verrucomicrobia;Verrucomicrobiae;Verrucomicrobiales;Akkermansiaceae;Akkermansia;Akkermansia muciniphila;Akkermansia muciniphila ATCC BAA-835

Akkermansia muciniphila 239935 cellular organisms;Bacteria;PVC group;Verrucomicrobia;Verrucomicrobiae;Verrucomicrobiales;Akkermansiaceae;Akkermansia;Akkermansia muciniphila

Mouse Intracisternal A-particle 11932 Viruses;Ortervirales;Retroviridae;unclassified Retroviridae;Intracisternal A-particles;Mouse Intracisternal A-particle

Wei Shen

uncultured murine large bowel bacterium BAC 54B 314101 cellular organisms;Bacteria;environmental samples;uncultured murine large bowel bacterium BAC 54B

Croceibacter phage P2559Y 1327037 Viruses;dsDNA viruses, no RNA stage;Caudovirales;Siphoviridae;unclassified Siphoviridae;Croceibacter phage P2559Y

 

cat names.txt | taxonkit name2taxid | taxonkit lineage --taxid-field 2 --show-lineage-taxids

$ cat names.txt | taxonkit name2taxid | taxonkit lineage --taxid-field 2 --show-lineage-taxids

Homo sapiens 9606 cellular organisms;Eukaryota;Opisthokonta;Metazoa;Eumetazoa;Bilateria;Deuterostomia;Chordata;Craniata;Vertebrata;Gnathostomata;Teleostomi;Euteleostomi;Sarcopterygii;Dipnotetrapodomorpha;Tetrapoda;Amniota;Mammalia;Theria;Eutheria;Boreoeutheria;Euarchontoglires;Primates;Haplorrhini;Simiiformes;Catarrhini;Hominoidea;Hominidae;Homininae;Homo;Homo sapiens 131567;2759;33154;33208;6072;33213;33511;7711;89593;7742;7776;117570;117571;8287;1338369;32523;32524;40674;32525;9347;1437010;314146;9443;376913;314293;9526;314295;9604;207598;9605;9606

Akkermansia muciniphila ATCC BAA-835 349741 cellular organisms;Bacteria;PVC group;Verrucomicrobia;Verrucomicrobiae;Verrucomicrobiales;Akkermansiaceae;Akkermansia;Akkermansia muciniphila;Akkermansia muciniphila ATCC BAA-835 131567;2;1783257;74201;203494;48461;1647988;239934;239935;349741

Akkermansia muciniphila 239935 cellular organisms;Bacteria;PVC group;Verrucomicrobia;Verrucomicrobiae;Verrucomicrobiales;Akkermansiaceae;Akkermansia;Akkermansia muciniphila 131567;2;1783257;74201;203494;48461;1647988;239934;239935

Mouse Intracisternal A-particle 11932 Viruses;Ortervirales;Retroviridae;unclassified Retroviridae;Intracisternal A-particles;Mouse Intracisternal A-particle 10239;2169561;11632;35276;11749;11932

Wei Shen

uncultured murine large bowel bacterium BAC 54B 314101 cellular organisms;Bacteria;environmental samples;uncultured murine large bowel bacterium BAC 54B 131567;2;48479;314101

Croceibacter phage P2559Y 1327037 Viruses;dsDNA viruses, no RNA stage;Caudovirales;Siphoviridae;unclassified Siphoviridae;Croceibacter phage P2559Y 10239;35237;28883;10699;196894;1327037

 

 

この他、シェルを自動補完してくれる"genautocomplete"がある。

引用

TaxonKit: a cross-platform and efficient NCBI taxonomy toolkit

Wei Shen, Jie Xiong

bioRxiv preprint first posted online Jan. 8, 2019;

 

Krona

 

 メタゲノム研究の重要な成果は、分類群または機能群の存在量の推定である。これらのグループへのアサインおける固有の不確実性は、それらの階層的コンテキストとそれらの予測信頼度の両方を考慮することを重要にしている。しかし、メタゲノムデータを視覚化するための現在のツールは、定量的な階層関係を省略または歪ませ、二次的な変数を表示するための機能を欠いている。

 ここでKrona、メタゲノム分類の複雑な階層内の相対存在量と信頼性の直感的な探査を可能にする新しい視覚化ツールを紹介する。 Kronaは、パラメトリックなカラーリングとインタラクティブ極座標ズーミングを備え、さまざまな放射状のスペースフィルディスプレイを組み合わせている。 HTML 5とJavaScriptの実装により、インストールされたソフトウェアやプラグインを必要とせずに、あらゆるモダンWebブラウザで閲覧できる完全にインタラクティブなチャートが可能になる。このWebベースのアーキテクチャでは、各チャートを独立したドキュメントにすることもできるため、電子メールで共有したり標準のWebサーバーに投稿したりすることが容易になる。 論文ではKronaの有用性を説明するために、さまざまなメタゲノムデータセットへのその適用と一般的なメタゲノム解析ツールとの互換性について説明する。Kronaレンダリングコードと変換ツールはどちらもBSDオープンソースライセンスの下で無料で提供され、http://krona.sourceforge.netから入手できる。

 

 

Kronaに関するツイート


インストール 

docker上のubuntu16.04でビルドし、動作確認した。

本体 Github

#Anacondaを使っているならcondaで導入できる
conda install -c bioconda -y krona

> ktImportTaxonomy

$ ktImportTaxonomy

                                                                                                                                                                                               ______________________________________________________________________________________________________________/ KronaTools 2.7 - ktImportTaxonomy \___

 

Creates a Krona chart based on taxonomy IDs and, optionally, magnitudes and scores. Taxonomy IDs corresponding to a rank of "no rank" in the database will be assigned to their parents to make the hierarchy less cluttered (e.g.

"Cellular organisms" will be assigned to "root").

                                                                                                                                                                                                                           __________________________________________________________________________________________________________________________________________/ Usage \___

 

ktImportTaxonomy \

   [options] \

   taxonomy_1[:magnitudes_1][,name_1] \

   [taxonomy_2[:magnitudes_2][,name_2]] \

   ...

 

   taxonomy    Tab-delimited file with taxonomy IDs and (optionally) query IDs, magnitudes and scores. By default, query IDs, taxonomy IDs and scores will be taken from columns 1, 2 and 3, respectively (see -q, -t, -s, and -m).

               Lines beginning with "#" will be ignored. By default, separate datasets will be created for each input (see [-c]).

 

   magnitudes  Optional file listing query IDs with magnitudes, separated by tabs. This can be used to account for read length or contig depth to obtain a more accurate representation of abundance. By default, query sequences

               without specified magnitudes will be assigned a magnitude of 1. Magnitude files for assemblies in ACE format can be created with ktGetContigMagnitudes.

 

   name        A name to show in the list of datasets in the Krona chart (if multiple input files are present and [-c] is not specified). By default, the basename of the file will be used.

                                                                                                                                                                                                                         ________________________________________________________________________________________________________________________________________/ Options \___

 

   [-o <string>]    Output file name. [Default: 'taxonomy.krona.html']

 

   [-n <string>]    Name of the highest level. [Default: 'Root']

 

   [-i]             Include a wedge for queries with no hits.

 

   [-c]             Combine data from each file, rather than creating separate datasets within the chart.

 

   [-q <integer>]   Column of input files to use as query ID. Required if magnitude files are specified. [Default: '1']

 

   [-t <integer>]   Column of input files to use as taxonomy ID. [Default: '2']

 

   [-s <integer>]   Column of input files to use as score. [Default: '3']

 

   [-m <integer>]   Column of input files to use as magnitude. If magnitude files are specified, their magnitudes will override those in this column.

 

   [-d <integer>]   Maximum depth of wedges to include in the chart.

 

   [-k]             Allow assignments to taxa with ranks labeled "no rank" (instead of moving up to parent).

 

   [-x <integer>]   Hue (0-360) for "bad" scores. [Default: '0']

 

   [-y <integer>]   Hue (0-360) for "good" scores. [Default: '120']

 

   [-u <string>]    URL of Krona resources to use instead of bundling them with the chart (e.g. "http://krona.sourceforge.net"). Reduces size of charts and allows updates, though charts will not work without access to this URL.

 

   [-qp <string>]   Url to send query IDs to (instead of listing them) for each wedge. The query IDs will be sent as a comma separated list in the POST variable "queries", with the current dataset index (from 0) in the POST

                    variable "dataset". The url can include additional variables encoded via GET.

 

   [-tax <string>]  Path to directory containing a taxonomy database to use. [Default: '/Users/kazuma/.pyenv/versions/miniconda3-4.0.5/opt/krona/taxonomy']

> ktImportText

$ ktImportText

                                                                                                                                                                                      

_____________________________________________________________________________________________________/ KronaTools 2.7 - ktImportText \___

 

Creates a Krona chart from text files listing quantities and lineages.

                                                                                                                                                                                                              

_____________________________________________________________________________________________________________________________/ Usage \___

 

ktImportText \

   [options] \

   text_1[,name_1] \

   [text_2[,name_2]] \

   ...

 

   text  Tab-delimited text file. Each line should be a number followed by a list of wedges to contribute to (starting from the highest level). If no wedges are listed (and just a quantity is given), it will

         contribute to the top level. If the same lineage is listed more than once, the values will be added. Quantities can be omitted if -q is specified. Lines beginning with "#" will be ignored. By default,

         separate datasets will be created for each input (see [-c]).

 

   name  A name to show in the list of datasets in the Krona chart (if multiple input files are present and [-c] is not specified). By default, the basename of the file will be used.

                                                                                                                                                                                                            

__________________________________________________________________________________________________________________________/ Options \___

 

   [-o <string>]  Output file name. [Default: 'text.krona.html']

 

   [-n <string>]  Name of the highest level. [Default: 'all']

 

   [-q]           Files do not have a field for quantity.

 

   [-c]           Combine data from each file, rather than creating separate datasets within the chart.

 

   [-u <string>]  URL of Krona resources to use instead of bundling them with the chart (e.g. "http://krona.sourceforge.net"). Reduces size of charts and allows updates, though charts will not work without access to

                  this URL.

 

 

データベース作成

condaパッケージマネージャを使ってインストールすると、最後に以下のようにメッセージが出る。

If you would like the taxonomic data stored elsewhere, simply replace

this directory with a symlink.  For example:

 

rm -rf /home/kazu/.pyenv/versions/miniconda3-4.0.5/opt/krona/taxonomy

mkdir /path/on/big/disk/taxonomy

ln -s /path/on/big/disk/taxonomy /home/kazu/.pyenv/versions/miniconda3-4.0.5/opt/krona/taxonomy

ktUpdateTaxonomy.sh

taxonomyデータベースをダウンロードする前に、以下のようにリンクを修正。ここではpyenvでpythonのバージョン管理しており、miniconda3-4.0.5に入れたkronaのデータベースリンクを修正している。

rm -rf ~/.pyenv/versions/miniconda3-4.0.5/opt/krona/taxonomy
mkdir -p ~/krona/taxonomy
ln -s ~/krona/taxonomy/ ~/.pyenv/versions/miniconda3-4.0.5/opt/krona/

 データベースダウンロード

ktUpdateTaxonomy.sh ~/krona/taxonomy

> head -n 30 ~/krona/taxonomy/taxonomy.tab

$ head -n 30 ~/krona/taxonomy/taxonomy.tab 

1 0 1 no rank root

2 2 131567 superkingdom Bacteria

6 7 335928 genus Azorhizobium

7 8 6 species Azorhizobium caulinodans

9 8 32199 species Buchnera aphidicola

10 7 1706371 genus Cellvibrio

11 9 1707 species Cellulomonas gilvus

13 7 203488 genus Dictyoglomus

14 8 13 species Dictyoglomus thermophilum

16 7 32011 genus Methylophilus

17 8 16 species Methylophilus methylotrophus

18 8 213421 genus Pelobacter

19 9 18 species Pelobacter carbinolicus

20 7 76892 genus Phenylobacterium

21 8 20 species Phenylobacterium immobile

22 7 267890 genus Shewanella

23 8 22 species Shewanella colwelliana

24 8 22 species Shewanella putrefaciens

25 8 22 species Shewanella hanedai

27 5 49928 species halophilic eubacterium NRCC 41227

28 5 49928 species halophilic eubacterium

29 6 28221 order Myxococcales

31 8 80811 family Myxococcaceae

32 9 31 genus Myxococcus

33 10 32 species Myxococcus fulvus

34 10 32 species Myxococcus xanthus

35 10 32 species Myxococcus macrosporus

38 10 47 species Archangium disciforme

39 8 80811 family Archangiaceae

40 9 39 genus Stigmatella

 

 

実行方法 

1、ktImportTaxonomy   分類学の階級に従って分類(link

フォーマットはクエリのID、NCBI taxonomy ID (taxIDs)、scoreの順番。異なる場合、下のランのようにフラグを立てて指定する。

f:id:kazumaxneo:20190120120821j:plain

kraken report出力を可視化する場合、NCBI taxonomy IDはカラム3、クエリのIDはカラム2になる。

f:id:kazumaxneo:20190120121022j:plain

よって、次のようにカラム指定する(*1)。

ktImportTaxonomy -q 2 -t 3 results.kraken -o output.plot.html
  • -q    Column of input files to use as query ID. Required if magnitude files are specified. [Default: '1']
  • -t    Column of input files to use as taxonomy ID. [Default: '2']

 output.plot.htmlが出力される。

f:id:kazumaxneo:20190120114722j:plain

 

2、ktImportText    テキストファイル入力から可視化(link

分類学の階級以外のデータの可視化に利用できる。例えば下のような食物の分類。

f:id:kazumaxneo:20190120115550j:plain

フォーマットは、左端のカラムが割合。2カラム目以降は使用している階級の上位から記載していく。例えば上のようにカラム2で脂質、カラム3で飽和脂肪酸、または不飽和脂肪酸

ランするには、用意したテキストファイルを指定する。

ktImportText results.kraken -o output.plot.html

f:id:kazumaxneo:20190120120340j:plain



 

引用

Interactive metagenomic visualization in a Web browser
Brian D Ondov, Nicholas H Bergman, Adam M Phillippy

BMC Bioinformatics. 2011; 12: 385.

 

参考

*1

またはcut -f 2-3でカラム2-3だけ取り出し、それをktImportTaxonomyの入力に使う。これでオプションのフラグ無しでランできる。

 

Question: Transforming kraken outputs to krona

Transforming kraken outputs to krona

 

multi-omicデータセットを可視化するCircosのラッパーツール MS-Helios

 

 マイクロアレイ、次世代シークエンシング、および質量分析(MS)などの革新的なハイスループット技術は、生物学的システムに対する我々の理解を大いに進歩させた。これらのすぐに利用可能で費用対効果に優れた包括的なデータ取得方法により、システム生物学はシングルオミクスからマルチオミクスデータ分析へと移行しつつある[論文より ref.1]。しかし、何千ものマルチオミクス分子プロファイルを統合して視覚化することは、システム生物学に新たな挑戦をもたらす。今日まで、ほとんどの多項式分析方法は、クラスタリング、相関[ref.2]、または次元削減方法、例えば視覚化前にデータを変換する主成分分析[ref.3]に依存している。

 Holisticで統合されたビューを提供するために、著者らは最近、プロテオーム中心の方法で円形のプロットでサンプルの特徴を視覚化する、 circular proteome maps (CPMs)を導入した[ref.4 link]。circular plotでは、ヒストグラム、散布図、折れ線プロットなどのよく知られているプロットタイプに依存して、直感的で審美的な方法で高次元データと特徴の関係を視覚化できる[ref.5、6]。さらに、データトラックは、複数のオミクスレベルにわたって特定の機能をコンテキスト化するための手段を提供する。circular plotを作成するための最も標準的なソフトウェアはCircosである。これはコマンドラインベースのPerlプログラムで、学習曲線がsteepである[ref.5]。円形のプロットの構築プロセスと視覚化を容易にするために、複数のRパッケージとツールが利用可能です[ref.6、7、8、9、10]。これらのツールは、ゲノムデータ用に構築されているか、他のデータソースをゲノムベースにマッピングしている。それらのいずれも、非ゲノムベースでのマルチオミクスデータ統合または視覚化を考慮していない。

 非ゲノムベースでcircular plotの作成を容易にするために、MS-Heliosと呼ばれるCircosラッパーを開発した。 MS-Heliosは、迅速なプロトタイピング、データ探索、そして簡単に高品質でpublish-ready な図を作成することを可能にするコマンドラインツールである。

 

 

 

インストール 

ubuntu16.04のminiconda3-4.0.5で動作確認した。

依存

#Anacondaを使っているならcodaで導入可能
conda install -c bioconda -y circos

本体 SourceForge

マニュアルPDFも本体に付属

f:id:kazumaxneo:20190119122844p:plain

> java -jar MS-Helios.jar -h 

$ java -jar MS-Helios.jar -h

usage: java -jar MS-Helios.jar -i <file> [OPTION] [<ARGUMENT>]

MS-Helios: A Circos wrapper to visualize multi-omic datasets

 -a,--add <plot>                define data track plot type; plot is

                                histogram, scatter, line, or highlight

 -c,--colors <scheme>           group colors, e.g., blues, brbg, accent,

                                see Circos brewer colors for more

                                [default: B/W]

 -d,--delimiter <character>     field delimiter character [default: ,]

 -e,--grid                      adds grid line to separate groups over

                                data tracks

 -f,--filter <method> <value>   apply filter; method is percentile v or

                                top-hat v

 -g,--group                     group features by sample occurrence

 -h,--help                      print this message

 -i,--input <file>              use given input file

 -l,--load <path>               load helios file [default:

                                MS-Helios.helios]

 -m,--missing <value>           missing value [default: 0]

 -n,--normalize <method>        use normalization; method is scale,

                                z_score, divide_min, divide_max,

                                divide_mean, divide_sd, or divide_sum

 -o,--output <path>             use given output path [default: jar path]

 -p,--partition                 plot samples separate

 -s,--sort <order>              sort data tracks; order is asc or dsc

 -t,--transform <method>        use transformation; method is shannon,

                                log2, or log10

 -v,--version                   print the version information

Please report bugs to: https://sourceforge.net/p/ms-helios/tickets

 

入力データフォーマット(マニュアルより)

f:id:kazumaxneo:20190119132024p:plain

 

テストラン

1、circosのランに必要なconfigファイル作成

cd MS-Helios-1/
java -jar MS-Helios.jar -i data/LFQ_protein.csv -g
  • -i       use given input file
  • -g      group features by sample occurrence

 

2、circosで可視化

circos –conf MS-Helios.config

 

引用

MS-Helios: a Circos wrapper to visualize multi-omic datasets

Harald Marx, Joshua J. Coon

BMC Bioinformatics201920:21

 

関連


pblat: マルチスレッドに対応したblat

 

 Blat [論文より ref.1 link]は、DNA、RNAおよびタンパク質配列をリファレンスゲノムにマッピングするように設計された配列アラインメントツールである。これは一般に、リファレンスゲノム内の配列の検索、closely relatedな種のゲノムからの相同配列の発見、mRNA配列からのエキソン - イントロンジャンクションの同定、および遺伝子構造の決定、ならびにゲノムおよびトランスクリプトーム配列の構築およびアノテーションを助けるために用いられる[ref.2]。 BWA [ref.3]やBowtie [ref.4]のような多くの高速シーケンスアライナーは、ハイスループットシーケンシングによって生成されたショートシーケンスリードマッピングするために開発されたが、ロングリードやギャップが豊富なシーケンスや不連続からのスプライスシーケンスのマッピングはできない。対照的に、blatは高感度と高精度を両立しており、そのようなアプリケーションに理想的なツールである[ref.67]。

 しかし、ハイスループットシークエンシングプロジェクトによって生成されるシークエンスの量が増えるにつれて、blatは大規模解析や定期的に更新されるアノテーションに必要な速度要件を満たすことができなくっている。例えば、脊椎動物の全トランスクリプトーム配列をリファレンスゲノムにマッピングするために使用される場合、blatを終了するまでに数日かかるだろう。これは、blatアルゴリズムがシングルスレッドであるため、最新のマルチコアプロセッサを最大限に活用できないためである。 GNU parallel [ref.8]ツールを使用して、1台以上のコンピューターを使用してblatの複数のインスタンスを並列に実行することはできる。しかしながら、各blatのプロセスは、全リファレンスゲノムのコピーをロードし、そしてゲノムのインデックスを構築してメモリに格納することになり、これは、複数のブラットプロセスが同時に実行される場合、従来のコンピュータの利用可能な物理メモリを超える可能性がある。

 これらの制限を克服するために、pblat(parallel blat)を提示する。これは、マルチスレッドおよびクラスタコンピューティングのサポートを実装することによって、blatのアライメントを高速化するように機能する。pblatでは、すべてのスレッドがリファレンスゲノム全体とインデックスの同じメモリコピーを共有する。そのため、pblatはblatとほぼ同じ量のメモリを使用する。実行時間は使用されるスレッドの数とともに減少し、pblatの出力結果はblatのそれと同じになる。 pblatのクラスタバージョン(pblat-cluster)はMPI(Message Passing Interface)をサポートするコンピュータクラスタ上で動作する。これにより、blatの実行時間を数日から数分に短縮することができる。

 

f:id:kazumaxneo:20190119101205p:plain

Fig.1Performance evaluation of pblat.  論文より転載。

 

インストール 

mac os10.14でビルドし、動作確認した。

pblat can run on Linux and Mac OS.

本体 Github

git clone https://github.com/icebert/pblat.git
cd  pblat/
make

> ./pblat

$ ./pblat 

pblat - BLAT with parallel supports v. 35 fast sequence search command line tool

 

usage:

   pblat database query [-ooc=11.ooc] output.psl

where:

   database and query are each a .fa file

   -ooc=11.ooc tells the program to load over-occurring 11-mers from

               and external file.  This will increase the speed

               by a factor of 40 in many cases, but is not required

   output.psl is where to put the output.

 

options:

   -t=type     Database type.  Type is one of:

                 dna - DNA sequence

                 prot - protein sequence

                 dnax - DNA sequence translated in six frames to protein

               The default is dna

   -q=type     Query type.  Type is one of:

                 dna - DNA sequence

                 rna - RNA sequence

                 prot - protein sequence

                 dnax - DNA sequence translated in six frames to protein

                 rnax - DNA sequence translated in three frames to protein

               The default is dna

   -prot       Synonymous with -t=prot -q=prot

   -ooc=N.ooc  Use overused tile file N.ooc.  N should correspond to 

               the tileSize

   -threads=N  number of threads to run

   -tileSize=N sets the size of match that triggers an alignment.  

               Usually between 8 and 12

               Default is 11 for DNA and 5 for protein.

   -stepSize=N spacing between tiles. Default is tileSize.

   -oneOff=N   If set to 1 this allows one mismatch in tile and still

               triggers an alignments.  Default is 0.

   -minMatch=N sets the number of tile matches.  Usually set from 2 to 4

               Default is 2 for nucleotide, 1 for protein.

   -minScore=N sets minimum score.  This is the matches minus the 

               mismatches minus some sort of gap penalty.  Default is 30

   -minIdentity=N Sets minimum sequence identity (in percent).  Default is

               90 for nucleotide searches, 25 for protein or translated

               protein searches.

   -maxGap=N   sets the size of maximum gap between tiles in a clump.  Usually

               set from 0 to 3.  Default is 2. Only relevent for minMatch > 1.

   -noHead     suppress .psl header (so it's just a tab-separated file)

   -makeOoc=N.ooc Make overused tile file. Target needs to be complete genome.

   -repMatch=N sets the number of repetitions of a tile allowed before

               it is marked as overused.  Typically this is 256 for tileSize

               12, 1024 for tile size 11, 4096 for tile size 10.

               Default is 1024.  Typically only comes into play with makeOoc.

               Also affected by stepSize. When stepSize is halved repMatch is

               doubled to compensate.

   -mask=type  Mask out repeats.  Alignments won't be started in masked region

               but may extend through it in nucleotide searches.  Masked areas

               are ignored entirely in protein or translated searches. Types are

                 lower - mask out lower cased sequence

                 upper - mask out upper cased sequence

                 out   - mask according to database.out RepeatMasker .out file

                 file.out - mask database according to RepeatMasker file.out

   -qMask=type Mask out repeats in query sequence.  Similar to -mask above but

               for query rather than target sequence.

   -repeats=type Type is same as mask types above.  Repeat bases will not be

               masked in any way, but matches in repeat areas will be reported

               separately from matches in other areas in the psl output.

   -minRepDivergence=NN - minimum percent divergence of repeats to allow 

               them to be unmasked.  Default is 15.  Only relevant for 

               masking using RepeatMasker .out files.

   -dots=N     Output dot every N sequences to show program's progress

   -trimT      Trim leading poly-T

   -noTrimA    Don't trim trailing poly-A

   -trimHardA  Remove poly-A tail from qSize as well as alignments in 

               psl output

   -fastMap    Run for fast DNA/DNA remapping - not allowing introns, 

               requiring high %ID. Query sizes must not exceed 5000.

   -out=type   Controls output file format.  Type is one of:

                   psl - Default.  Tab separated format, no sequence

                   pslx - Tab separated format with sequence

                   axt - blastz-associated axt format

                   maf - multiz-associated maf format

                   sim4 - similar to sim4 format

                   wublast - similar to wublast format

                   blast - similar to NCBI blast format

                   blast8- NCBI blast tabular format

                   blast9 - NCBI blast tabular format with comments

   -fine       For high quality mRNAs look harder for small initial and

               terminal exons.  Not recommended for ESTs

   -maxIntron=N  Sets maximum intron size. Default is 750000

   -extendThroughN - Allows extension of alignment through large blocks of N's

 

実行方法

リファレンスFASTAと入力のFASTAファイルを指定する(fastqには対応していない)。

pblat -threads=12 ref.fa input.fasta output
  • -threads=<N>   number of threads to run

> head output

f:id:kazumaxneo:20190119104004p:plain

 

samに変換するには以下を参考にしてください。

Question: Conversion Of Blat Output To Sam/Bam

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

引用

pblat: a multithread blat algorithm speeding up aligning sequences to genomes

Meng Wan, Lei Kong

BMC Bioinformatics 2019 20:28

 

 

 

ロングリードのシミュレーションやロングリードのエラーコレクションツールの評価を行う ELECTOR

 

 Pacific Biosciences(PB)とOxford Nanopore Technologies(ONT)のロングリードは、高いエラーレートと複雑なエラープロファイルにもかかわらず、さまざまなアプリケーションに急速に採用されてきている[論文 ref.1]。これらのリードは、エラー率が高く(テクノロジやライブラリによる9%から最大30%)、イルミナのリードをはるかに上回る。さらに、大部分のエラーが置換であるIlluminaとは対照的に、ロングリードのエラーには挿入および欠失(indels)エラーが主要に含まれる(ONTリードはより欠失されやすいのに対して、PBリードはより多く挿入を含む)。この問題には、新規かつ具体的なアルゴリズムの開発が必要である。この範囲で、これらの主にを直接対象とした数十のエラー訂正方法がこの5年間で出現した。ハイブリッド訂正器と呼ばれる第1の範囲のエラー訂正ツールは、ロングリードシーケンスを強化するためにショートリードの低いエラー率に頼って、エラー訂正を実行する。自己訂正と呼ばれる2番目のグループのメソッドは、唯一シーケンスに含まれる情報を使ってロングリードを修正することを目的としている(レビューについては[ref.2]を参照)。どちらのパラダイムにも、非常に多様なアルゴリズムソリューションが含まれているため、適切なベンチマークなしでは、訂正結果を(スループット、クオリティ、およびパフォーマンスに関して)グローバルに比較することは困難である。
 アセンブリや構造変異バリアントコールなど、ロングレンジ情報から利益を得るためにロングリードを採用するプロジェクトが増えている[ref.1]。ロングリードのエラー率を考えると、多くのアプリケーションの最初のステップはエラー訂正である。ただし、この段階は時間のボトルネックになる可能性がある[ref.1]。
加えて、エラー訂正品質は下流のプロセスにかなりの影響を及ぼす。したがって、どの訂正器が特定の実験計画(例えば、適用範囲、リードタイプ、またはゲノム)に最も適しているかを事前に知っておくことは興味深い。したがって、正確で信頼できる統計を用いてエラー訂正ツールを評価することを可能にする方法を開発することは、極めて重要である。
 そのような評価方法は、再現可能で包括的なベンチマークを実行することを可能にすべきであり、したがって、どのエラー訂正方法が所与のケースに最も適しているかを効率的に識別することを可能にするはずである。遭遇する可能性があるさまざまなシナリオを再現するために、それらはさまざまな複雑さ(細菌から真核生物まで)のデータセットで使用可能でなければならない。それらはまた速くて軽量であるべきであり、そしてこれらが評価する実際の訂正方法よりもコンピュータリソースおよび時間を消費するべきではない。この点は、訂正器の評価者が新しい訂正方法の開発の観点から立つ場合に特に重要である。それらは、最先端の訂正器と正確かつ迅速な比較を提供するのに役立つ。ユーザーだけでなく開発者にとっても、訂正器評価者は、その潜在的な落とし穴を特定するために、訂正器の振る舞い(すなわち、修正された塩基数、導入されたエラー、スループット)を正確に記述する必要がある。

 

対応したエラーコレクションツール

  • Proovread
  • LoRDEC
  • Nanocorr
  • NaS
  • CoLoRMap
  • HG-CoLoR
  • HALC
  • PBDAGCon
  • Canu
  • LoRMA
  • daccord
  • MECAT

対応したシミュレータ

  • SimLord

  • NanoSim

 

インストール 

docker上のubuntu16.04でビルドし、動作確認した。

依存

  • gcc and C++11
  • Python3 and Biopython
  • R
#PDFレポートを作成するなら必要
sudo apt update && sudo apt install texlive-latex-base

本体 Github

git clone --recursive https://github.com/kamimrcht/ELECTOR
cd ELECTOR/
./install.sh

> python3 elector.py 

$ python3 elector.py 

usage: elector.py [-h] [-threads [THREADS]] [-corrected [CORRECTED]] [-split]

                  [-uncorrected [UNCORRECTED]] [-perfect [PERFECT]]

                  [-reference [REFERENCE]] [-simulator [SIMULATOR]]

                  [-corrector [SOFT]] [-dazzDb [DAZZDB]]

                  [-output [OUTPUTDIRPATH]] [-remap] [-assemble]

                  [-minsize [MINSIZE]]

 

optional arguments:

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

  -threads [THREADS]    Number of threads

  -corrected [CORRECTED]

                        Fasta file with corrected reads (each read sequence on

                        one line)

  -split                Corrected reads are split

  -uncorrected [UNCORRECTED]

                        Prefix of the reads simulation files

  -perfect [PERFECT]    Fasta file with reference read sequences (each read

                        sequence on one line)

  -reference [REFERENCE]

                        Fasta file with reference genome sequences (each

                        sequence on one line)

  -simulator [SIMULATOR]

                        Tool used for the simulation of the long reads (either

                        nanosim or simlord)

  -corrector [SOFT]     Corrector used (lowercase, in this list: lorma, mecat,

                        pbdagcon, daccord). If no corrector name is provided,

                        make sure the read's headers are correctly formatted

                        (i.e. they correspond to those of uncorrected and

                        reference files)

  -dazzDb [DAZZDB]      Reads database used for the correction, if the reads

                        were corrected with Daccord or PBDagCon

  -output [OUTPUTDIRPATH]

                        Name for output directory

  -remap                Perform remapping of the corrected reads to the

                        reference

  -assemble             Perform assembly of the corrected reads

  -minsize [MINSIZE]    Do not assess reads/fragments chose length is <=

                        MINSIZE % of the original read

 

テストラン

#1

python3 elector.py -reference example/example_reference.fasta \
-corrected example/Simlord/correctedReads.fasta \
-uncorrected example/Simlord/simulatedReads \
-simulator simlord

summary.pdf

f:id:kazumaxneo:20190116163529j:plain

#2

python3 elector.py -perfect example/Simlord/simulatedReads_reference.fasta -corrected example/Simlord/correctedReads.fasta -uncorrected example/Simlord/simulatedReads.fasta

f:id:kazumaxneo:20190116164158j:plain

 

シミュレーションコマンドもテストして追記します。

 

引用

ELECTOR: Evaluator for long reads correction methods

bioRxiv preprint first posted online Jan. 7, 2019

Camille Marchet, Pierre Morisse, Lolita Lecompte, Antoine Limasset, Arnaud Lefebvre, Thierry Lecroq, Pierre Peterlongo

 

関連

 

 

 

Pacbioのロングリードのエラーコレクションツール pbdagcon

 

 イルミナなどの第2世代シーケンシング(2GS)プラットフォームは、ゲノムシークエンシングコストを劇的に削減しながら、スループットを飛躍的に向上させた(Shendure and Ji 2008)。 2GSプラットフォームの比較的低コストで大規模なスループットは、数千もの種のゲノムのシーケンシング、mおよび新たなアセンブリのための道を切り開いた(Alkan et al、2011)。

(一部略)ショートリードでは、リピートDNA配列は各リード長よりもはるかに長いことが多いため、de novoアセンブリは本質的に困難な計算上の問題である(Ukkonen 1992)。例えば、リピートDNA配列は誤ってアセンブリされ断片化された領域の数を増加させることがあるので、ショートリードのde novoアセンブリは配列情報の20%まで失われる可能性があると推定されている(Schatz et al。2010 link; Alkan et al。2011 link; Ukkonen 1992)。デノボアセンブリプロセスにおけるリピート DNA問題を軽減する1つの方法は、非常に長いインサート(>2kbp)を有する第2セットのメイトペアライブラリを組み込むことであった(Li et al、2010; Chaison et al、2009;Simpsonet al、2009; Alkan et al、2011; Butler et al、2008)。メイトペアライブラリーはリピートを解消し(Treangen and Salzberg 2012; Wetzel et al。2011)、scaffoldingを改善することができるが(van Heesch et al。2013)、ペアエンドのコンタミネーションとインサートサイズの誤推定もミスアセンブリを招く(Phillippy et al、 2008 link; Sahlin et al、2016 link)。

 最近では、Pacific Biosciences(PacBio)のSMRTシークエンスや、Oxford NanoporeのMinIONシーケンス(現在、はるかに長い最大54 kbpのリード長(論文執筆時点))などの第3世代(3GS)1分子シークエンシングテクノロジーは、2GS de novo assemblyのいくつかの欠点を克服することができる(Berlin et al、2015)。ロングリードシーケンシング技術は82.1%(Chin et al、2011)から84.6%の精度(Rasko et al、2011)の範囲の高いエラー率でリードを生成するが、シーケンシングエラーはロングリード全体でランダムな位置で起こる(Chin et al。2013)。2GSのショートリードデータ(Koren et al、2012)または過剰な3GS readを使用するセルフエラーコレクションでエラーは補正できる(Chin et al、2013)。

 本論文では、進化論的観点から特によく研究されているショウジョウバエDrosophila serrataのゲノムを新規にアセンブリするために、PacBioロングリードシークエンシングを使用する。 D. serrataD. monlanumサブグループのメンバーであり、これはD. melanogasterサブグループ〜40MYAから分かれており(Tamura et al、2004)、推定98種からなる(Brake andBächli、2008)。現在のところ、この種が豊富なサブグループからは、1つのドラフトゲノムアセンブリD. kikkawai)しか入手できない(Chen et al

、2014)。 D. serrataは、パプアニューギニアからオーストラリア南東部までの幅広い地理的分布を持ち、種の境界の進化などの進化の問題(Blows and Hoffman 1993; Hallas et al。2002; Magiafoglou et al、 2002)および気候適応(Frentiu and Chenoweth 2010; Latimer et al。2011; Kellermann et al。2009)に対処するための強力なモデルとして浮上している。

(一段落略)

 進化的研究のためのモデルとしてのD. serrataの重要性にもかかわらず、そのゲノムについての我々の不十分な理解は依然として重大な制限である。染色体の連鎖地図および物理地図が利用可能であり(Stocker et al、2012)、発現配列タグ(EST)ライブラリーが開発されている(Frentiu et al、2009)が、この種はドラフトゲノムを欠いている。ここでは、もっぱらPacBio SMRT技術を使用してD. serrataゲノムのシーケンシングおよびアセンブリを報告する。経験的なRNA seqデータによって裏付けられたin silco遺伝子予測因子に基づくゲノムの最初のアノテーションも提供する。我々(著者ら)のde novoゲノム配列とそのアノテーション情報は、この種の進行中の集団ゲノムと形質マッピング研究のためのリソースを提供するだけでなく、ショウジョウバエ科のゲノム進化のより広範な研究を促進する。

 

pbdagconに関するツイート

 

インストール

ubuntu14.04、miniconda3-4.0.5環境でテストした。

依存

  • blasr (optional)

The code now depends on C++11 features, in particular std::thread, std::move. GCC 4.8.1 or higher is known to work.

本体 Github

#Anacondaを使っているならcondaで導入可
conda install -c bioconda -y pbdagcon

#blasrもないなら入れておく
conda install -c bioconda -y blasr

> pbdagcon

$ pbdagcon --help

 

USAGE: 

 

   pbdagcon  [-v] [-a] [-t <uint>] [-m <uint>] [-c <uint>] [-j <int>] [--]

             [--version] [-h] <either file path or stdin>

 

 

Where: 

 

   -v,  --verbose

     Turns on verbose logging

 

   -a,  --align

     Align sequences before adding to consensus

 

   -t <uint>,  --trim <uint>

     Trim alignments on either size

 

   -m <uint>,  --min-length <uint>

     Minimum length for correction

 

   -c <uint>,  --min-coverage <uint>

     Minimum coverage for correction

 

   -j <int>,  --threads <int>

     Number of consensus threads

 

   --,  --ignore_rest

     Ignores the rest of the labeled arguments following this flag.

 

   --version

     Displays version information and exits.

 

   -h,  --help

     Displays usage information and exits.

 

   <either file path or stdin>

     (required)  Input data

 

 

   PBI consensus module

 

 

 

実行方法

2、blasrを使い、リファレンスまたはアセンブリして得たdcontig配列にロングリードをマッピングする。

blasr queries.fasta target.fasta --bestn 1 -m 5 --out mapped.m5
pbdagcon mapped.m5 > consensus.fasta

 

 コンセンサス配列を出力するには、blasrによるアライメント結果を使用する

pbdagcon mapped.m5 > consensus.fasta

 

 

引用
Single-Molecule Sequencing of the Drosophila serrata Genome
Scott L. Allen,* Emily K. Delaney,† Artyom Kopp,† and Stephen F. Chenoweth

G3 (Bethesda). 2017 Mar; 7(3): 781–788. Published online 2017 Jan 30

 

ONTのロングリードのハイブリッドエラーコレクションツール HG-CoLoR

 

 最近のPacific Biosciences やOxford Nanoporeのようなロングリードシーケンシング技術は、ショートリード技術で許容されるより大きくて複雑なゲノムのアセンブリ問題を解決する。しかし、これらのロングリードは非常にノイジーで、Pacific Biosciencesでは約10〜15%、Oxford Nanoporeでは最大30%のエラー率に達する。エラー訂正の問題は、ロングリードを自己訂正するか、ショートリードで補完するハイブリッドの手法を使って取り組まれてきたが、ほとんどのメソッドはPacific Biosciencesのデータにのみ焦点を当て、Oxford Nanoporeのリードには適用されない。さらに、Oxford Nanoporeの最近のケミストリーでは、エラー率を15%以下に下げることが約束されているが、それは実際にはまだ高く(論文執筆時点)、このようなノイズの多いロングリードを修正することは依然として課題である。
 HG-CoLoRは、ショートリードとロングリードのアライメント、可変オーダーのBruijn graphのトラバースに基づくシード・アンド・エクステンデッド・アプローチに重点を置いたハイブリッド・エラー訂正方法である。著者らの実験は、HG-CoLoRが44%という高いエラー率のOxford Nanoporeのロングリードを効率的に修正することを示している。他の最先端のロングリードエラー訂正方法と比較し、HG-CoLoRが実行時間と結果の品質の間で最良のトレードオフを提供し、真核生物に効率的に拡張できる唯一の方法であることを示す。

 

HG-CoLoRに関するツイート


インストール

ubuntu16.04に導入した。

依存

本体 Github

git clone https://github.com/pierre-morisse/HG-CoLoR 
git submodule init
git submodule update
cd KMC/ && make -j
cd ../PgSA/ && make build CONF=pgsalib
cd .. && make

#Anaconda環境ならcondaで導入可能(linux only
conda install -c bioconda hg-color

./HG-CoLoR

$ ./HG-CoLoR

Usage: ./HG-CoLoR [options] --longreads LR.fasta --shortreads SR.fastq --out resultPrefix -K maxK

 

Note: HG-CoLoR default parameters are adapted for a 50x coverage set of short reads with a 1% error rate.

Please modify the parameters, in particular the --solid and --bestn ones, as indicated below

if using a set of short reads with a much higher coverage and/or a highly different error rate.

 

Input:

LR.fasta:                   fasta file of long reads, one sequence per line.

SR.fastq:                   fastq file of short reads.

                            Warning: only one file must be provided.

                            If using paired reads, please concatenate them into one single file.

resultPrefix:               Prefix of the fasta files where to output the corrected, trim and split long reads.

maxK:                       Maximum K-mer size of the variable-order de Bruijn graph.

 

Options:

--minorder INT, -k INT:     Minimum k-mer size of the variable-order de Bruijn graph (default: K/2).

--solid INT, -S INT:        Minimum number of occurrences to consider a short read k-mer as solid, after correction (default: 1).

                            This parameter should be carefully raised accordingly to the short reads coverage and accuracy,

                            and to the chosen maximum order of the graph.

                              It should only be increased when using high coverage of short reads, or a small maximum order.

--seedsoverlap INT, -o INT: Minimum overlap length to allow the merging of two overlapping seeds (default: K - 1).

--branches INT, -b INT:     Maximum number of branches exploration (default: 1,250).

                            Raising this parameter will result in less split corrected long reads.

                            However, it will also increase the runtime, and may create chimeric linkings between the seeds.

--seedskips INT, -s INT:    Maximum number of seed skips (default: 3).

--mismatches INT, -m INT:   Allowed mismatches when attempting to link two seeds together (default: 3).

--bestn INT, -n INT:        Top alignments to be reported by BLASR (default: 50).

                            This parameter should be raised accordingly to the short reads coverage.

                            Its default value is adapted for a 50x coverage of short reads.

                              It should be decreased with higher coverage, and increased with lower coverage.

--nproc INT, -j INT:        Number of processes to run in parallel (default: number of cores).

--tmpdir STRING, -t STRING: Path where to store the directory containing temporary files (default: working directory).

--kmcmem INT, -r INT:       Maximum amount of RAM for KMC, in GB (default: 12).

--help, -h:                 Print this help message.

 

 

実行方法

 ショートリードとるんグリードを指定してランする。ペアエンドのショートリードはあらかじめインターレースにして、1ファイルで指定する必要がある。

HG-CoLoR --longreads long_reads.fasta --shortreads shortreads.fq --out resultPrefix -K <maxK>
  • -K    maxK: Maximum K-mer size of the variable-order de Bruijn graph.

  •  --shortreads   fastq file of short reads. Warning: only one file must be provided. If using paired reads, please concatenate them into one single file. 

     

 

引用
Hybrid correction of highly noisy long reads using a variable-order de Bruijn graph
Morisse P, Lecroq T, Lefebvre A

Bioinformatics. 2018 Dec 15;34(24):4213-4222