macでインフォマティクス

macでインフォマティクス

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

メタゲノムデータの高速なtaxonomy assignmentを行う kraken

metagenome 高速なツール 生物種の推定 (taxonomic profiling) Genome Biology 2014 taxonomic assignment abundance estimation in metagenomics data docker

2018 10/6 タイトル修正 

2018 11/17 簡単なテスト追加

2019 4/12 dockerリンク追加

2019 10/15リンク追加

 

krakenは2014年に発表されたメタゲノムデータの分類手法。fastqまたはfastaの入力からk-merの配列に分解し、構築したデータベースにアライメントを行う。BLASTと同等の精度を保ちながら、megablastより最大909倍高速と主張されている。似たツールにメタゲノムデータからtaxonomy情報を得るMetaPhlAnがあるが、MetaPhlAnと比較してもkrakenは最大11倍高速とされる。

 

公式ページ

https://ccb.jhu.edu/software/kraken/

 マニュアル

http://ccb.jhu.edu/software/kraken/MANUAL.html

 

インストール

cent OSに導入した。

ランにはJellyfishのver.1が必要。公式サイトからダウンロードしてビルドする。

wget http://www.cbcb.umd.edu/software/jellyfish/jellyfish-1.1.11.tar.gz
tar zxvf jellyfish-1.1.11.tar.gz
cd jellyfish-1.1.11/
./configure
make 
sudo make install

bin/jellyfishにパスを通しておく。

 

kraken本体はbrewで導入できる。

brew install kraken

または

git clone  https://github.com/DerrickWood/kraken.git
cd kraken/
./install_kraken.sh $KRAKEN_DIR #インストールディレクトリを指定

Anacondaを使っているならcondaで導入する。

#v1.1
conda install -c bioconda -y kraken==1.1

#v2 (bioconda)
conda install -c bioconda -y kraken2

 

NCBI refseqのURLが変更されているので、wgetで配列をダウンロードできなくなっている (2017年夏テスト時)。スクリプトを別途ダウンロードして修正する。

wget https://github.com/DerrickWood/kraken/blob/master/scripts/download_genomic_library.sh

download_genomic_library.shを開き、 43行目を以下のように変更。

wget $FTP_SERVER/genomes/archive/old_refseq/Bacteria/all.fna.tar.gz

さらに59行目を以下のように変更。

wget $FTP_SERVER/genomes/archive/old_refseq/Plasmids/plasmids.all.fna.tar.gz

 

brewでインストールしたkaken-buildを開き、

vi /usr/local/bin/kraken-build #viで開く。

262行目を修正したdownload_genomic_library.shのパスに変更する。ここでは以下のように修正した。

exec "/Volumes/3TB3/kraken_database/download_genomic_library.sh", $type;

保存して閉じる。

 

追記

dockerイメージ

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

 

ラン

 データベースの準備（初回だけ）

　↓

ラン

　↓

集計

という三段階の流れで進める。

 

--データベースの準備--

データベースのダウンロード。ここではdatabaseというディレクトリ名にした。

kraken-build --standard --db database

 database/にRefseqの全データとtaxonomy情報がダウンロードされる。

 

データベースのビルド。時間がかかるので、コア数があるだけ指定する。

kraken-build --standard --threads 22 --db database

環境によってはかなりの時間がかかる。

容量を確認。

user]$ du -sh database/

161G database/

2017年８月実行時のdatabase容量は161GBだった。

準備ができたらランする。

 

--taxonomyラベリングを行う--

入力のデフォルトはfastaだが、--fastq-inputのフラグをつけることでfastqにも対応している。また--gzip-compressedか--bzip2-compressedをつけることで、gzipやbzip2の入力も可能である。

 fastaを指定してラン。

kraken --preload --threads 20 --fastq-input --db database input.fa > output
  • --db　Name for Kraken DB (default: none)
  • --threads　Number of threads (default: 1)
  • --fasta-input　Input is FASTA format
  • --preload　Loads DB into memory before classification

 --preloadをつけることで、データベースがRAMディスクに読み込まれ高速化する。databaseは自分がつけたデータベースディレクトリ名。

  

ペアードエンドfastqを指定してラン。

kraken --preload --threads 20 --fastq-input --paired --check-names --db database R1.fq R2.fq > output
  • --fastq-input　Input is FASTQ format 
  • --paired　The two filenames provided are paired-end reads
  • --check-names　Ensure each pair of reads have names that agree with each other; ignored if --paired is not specified
  • --gzip-compressed　Input is gzip compressed
  • --bzip2-compressed　Input is bzip2 compressed
  • --quick　Quick operation (use first hit or hits)

krakenはペアリードをマージして１つにしてからk-merベースの探索を行う（NNNでつなぐ）。探索時、ATGC以外があるk-mer配列は除外される。そのためペアリードを入力に使うことで、シングルリードよりおよそ3%精度が上がるとされる。

  

出力は以下のようになる。

bin]$ head output

C M02077:58:000000000-APPA6:1:1101:15651:1463 1148 603 0:1 1148:16 0:31 1148:83 0:140 A:31 0:16 A:31 1148:4 0:220

C M02077:58:000000000-APPA6:1:1101:17046:1522 1148 601 0:1 1148:11 0:31 1148:61 0:31 1148:82 0:54 A:31 1148:11 0:5 A:31 1148:5 0:217

C M02077:58:000000000-APPA6:1:1101:15487:1536 1148 603 1148:14 0:72 1148:22 0:163 A:31 0:16 A:31 0:224

C M02077:58:000000000-APPA6:1:1101:17442:1545 1148 603 0:1 1148:113 0:31 1148:44 0:31 1148:16 0:35 A:31 1148:16 A:31 1148:36 0:188

C M02077:58:000000000-APPA6:1:1101:22458:1603 1148 601 0:16 1148:16 0:44 1148:14 0:181 A:31 0:1 1148:7 0:8 A:31 0:222

C M02077:58:000000000-APPA6:1:1101:19281:1604 1148 602 0:1 1148:64 0:31 1148:24 0:31 1148:28 0:31 1148:10 0:50 A:31 1148:16 A:31 1148:31 0:193

C M02077:58:000000000-APPA6:1:1101:23186:1616 1148 601 0:1 1148:168 0:31 1148:5 0:66 A:31 1148:16 A:31 0:5 1148:21 0:196

C M02077:58:000000000-APPA6:1:1101:18130:1621 1148 603 1148:213 0:58 A:31 1148:16 A:31 1148:81 0:31 1148:31 0:81

C M02077:58:000000000-APPA6:1:1101:9815:1622 1148 599 1148:208 0:62 A:31 1148:16 A:31 1148:46 0:175

C M02077:58:000000000-APPA6:1:1101:21784:1636 1148 602 0:1 1148:130 0:31 1148:83 0:26 A:31 1148:16 A:31 1148:33 0:190

１列目　C or Uは、分類されたか（C）、されなかった（U）かを表す。

２列目　シーケンスID（インプットのヘッダー名）

３列目　taxonomy ID（wiki）

４列目　リードサイズ

５列目　↓のマニュアルの説明参照。

 

  1. A space-delimited list indicating the LCA mapping of each k-mer in the sequence. For example, "562:13 561:4 A:31 0:1 562:3" would indicate that:
    • the first 13 k-mers mapped to taxonomy ID #562
    • the next 4 k-mers mapped to taxonomy ID #561
    • the next 31 k-mers contained an ambiguous nucleotide
    • the next k-mer was not in the database
    • the last 3 k-mers mapped to taxonomy ID #562

 

 --taxonomy nameに変換--

taxonomy IDからtaxonomy nameに変えるコマンドが用意されている。以下のように打つ。

kraken-translate --db database output > sequences.labels

出力

f:id:kazumaxneo:20170830114628j:plain

入力が純粋培養のバクテリアのデータなので、同じ生き物ばかり表示されている。

 

--mpa-formatをつけると、superkingdom, kingdom, phylum, class, order, family, genus, speciesだけの表記になる。

kraken-translate --mpa-format --db database output > sequences.labels

f:id:kazumaxneo:20170830115002j:plain

 

  --集計--

集計レポートを作る。

kraken-report --db database output > summary

 出力を確認。

head -20 summary 

  0.71 2036 2036 U 0 unclassified

 99.29 285920 3 - 1 root

 99.29 285916 1 - 131567   cellular organisms

 99.29 285915 65 D 2     Bacteria

 98.87 284692 36 - 1783272       Terrabacteria group

 98.55 283771 0 - 1798711         Cyanobacteria/Melainabacteria group

 98.55 283771 3 P 1117           Cyanobacteria

 98.54 283759 2 O 1890424             Synechococcales

 98.54 283753 0 F 1890428               Merismopediaceae

 98.54 283753 0 G 1142                 Synechocystis

 98.54 283753 283695 S 1148                   Synechocystis sp. PCC 6803

  0.02 57 57 - 1080228                     Synechocystis sp. PCC 6803 substr. GT-I

  0.00 1 1 - 1080229                     Synechocystis sp. PCC 6803 substr. PCC-N

  0.00 3 0 F 1890436               Pseudanabaenaceae

  0.00 3 0 G 1152                 Pseudanabaena

  0.00 3 3 S 82654                   Pseudanabaena sp. PCC 7367

  0.00 1 0 F 1890426               Synechococcaceae

  0.00 1 0 G 1129                 Synechococcus

  0.00 1 1 S 110662                   Synechococcus sp. CC9605

  0.00 7 1 - 1301283             Oscillatoriophycideae

１列目　そのcladeがカバーされたリードの割合（全体に対するパーセンテージ）

２列目　そのcladeがカバーされたリード数

３列目　そのtaxonに直接アサインされたリード数

４列目　rank code（↓の説明参照）

５列目　NCBI taxonomy ID

６列目　scientific name

  1. A rank code, indicating (U)nclassified, (D)omain, (K)ingdom, (P)hylum, (C)lass, (O)rder, (F)amily, (G)enus, or (S)pecies. All other ranks are simply '-'.

 

 

MetaPhlAnのようなタブ仕分けのファイルを出力させるコマンドも用意されている。

kraken-mpa-report --db database output > summary_tab

bin]$ head -10 summary_tab 

d__Viruses 1

d__Bacteria 285915

d__Bacteria|p__Proteobacteria 419

d__Bacteria|p__Fusobacteria 5

d__Bacteria|p__Deferribacteres 2

d__Bacteria|p__Spirochaetes 7

d__Bacteria|p__Synergistetes 1

d__Bacteria|p__Bacteroidetes 724

d__Bacteria|p__Firmicutes 625

d__Bacteria|p__Deinococcus-Thermus 1

 

 

手持ちのfastaを追加して、カスタムデータベースを作ることもできる。作り方は公式マニュアルを参照。カスタムデータベース作成の流れが丁寧に説明されてます。

 

 

 

テスト

SRAの湖水からのメタゲノムデータでテストしてみた。15 GB x 2のデータをthread40でテストしたところ３０分ほどかかった。kraken-reportの出力の先頭を貼っておく。

bin]$ head -30 summary2

 74.03 29617719 29617719 U 0 unclassified

 25.97 10391739 6719 - 1 root

 25.91 10364486 1126 - 131567   cellular organisms

 25.90 10361130 51629 D 2     Bacteria

 25.37 10150686 3934 - 1783272       Terrabacteria group

 25.22 10090850 0 - 1798711         Cyanobacteria/Melainabacteria group

 25.22 10090850 129569 P 1117           Cyanobacteria

 22.78 9114469 8754 - 1301283             Oscillatoriophycideae

 22.59 9036657 21480 O 1150               Oscillatoriales

 21.70 8681472 0 F 1892254                 Oscillatoriaceae

 21.70 8681472 934 G 1158                   Oscillatoria

 21.59 8639339 0 S 482564                     Oscillatoria nigro-viridis

 21.59 8639339 8639339 - 179408                       Oscillatoria nigro-viridis PCC 7112

  0.10 41199 0 S 118323                     Oscillatoria acuminata

  0.10 41199 41199 - 56110                       Oscillatoria acuminata PCC 6304

  0.39 155174 0 F 1892249                 Cyanothecaceae

  0.39 155174 23122 G 43988                   Cyanothece

  0.15 59279 59279 S 497965                     Cyanothece sp. PCC 7822

  0.12 46609 46609 S 65393                     Cyanothece sp. PCC 7424

  0.03 12335 12335 S 43989                     Cyanothece sp. ATCC 51142

  0.01 5589 5589 S 395962                     Cyanothece sp. PCC 8802

  0.01 4347 4347 S 395961                     Cyanothece sp. PCC 7425

  0.01 3893 3893 S 41431                     Cyanothece sp. PCC 8801

  0.29 114192 578 F 1892252                 Microcoleaceae

  0.13 53893 0 G 35823                   Arthrospira

  0.13 53893 0 S 118562                     Arthrospira platensis

  0.13 53893 53893 - 696747                       Arthrospira platensis NIES-39

  0.09 36770 0 G 44471                   Microcoleus

  0.09 36770 36770 S 1173027                     Microcoleus sp. PCC 7113

  0.06 22951 0 G 1205                   Trichodesmium

unclassifiedが74%あったが、淡水性のシアノバクテリアなどが様々hitしている。

 

 追記

こちらも参考にしてください。

追記

Visualization tool

 

追記

簡単なテスト

in silico mock metagenome cummunityのデータを作成し、テストしてみた。以下の比率で10のゲノムを混ぜ込む。

f:id:kazumaxneo:20181120110800j:plain

Kronaで可視化

f:id:kazumaxneo:20181120110927j:plain

2%ほどしか混ぜ込んでいないRickettsia canadensisも検出できている（*1）。どの種の定量結果も、妥当な数値が出ている。

 

こちらも参考になります。

追記

Building a Kraken database with new FTP structure and no GI numbers

http://www.opiniomics.org/building-a-kraken-database-with-new-ftp-structure-and-no-gi-numbers/

 

引用

Kraken: ultrafast metagenomic sequence classification using exact alignments

Derrick E WoodEmail author and Steven L Salzberg

Genome Biology 2014 15:R46

 

ダウンロード時のエラーについて。

https://github.com/DerrickWood/kraken/issues/40

 

*1

リード数は50万x2。ファイルサイズは313MBx2と、メタゲノムにしては非常に少ない。

 

関連

krakenの出力を分析

krakenを利用

マルチサンプル


2023/03/04追記

 "Kraken2とMetaPhlAn 3を使用して、ヒト由来または環境由来のメタゲノム内のリードを分類したところ、分類されたリードの割合と同定された種の数の両方に大きな相違があることがわかりました。次に、様々な模擬サンプルや模擬サンプルを用いて、どのツールがメタゲノムサンプルの実際の組成に最も近い分類を与えるかを調べ、ツール、パラメータ、データベースの選択の組み合わせが、与えられた分類に与える影響を検討しました。"

https://pubmed.ncbi.nlm.nih.gov/36867161/