2019-02-01

種の存在量を推定する Bracken

2023/0710 merge追記

2024/08/22 追記

　メタゲノミクスは、急速にそして安価に膨大な量のDNA配列を生成する能力に部分的に起因して、急速に成長している研究分野である。ヒトゲノムが2001年に最初に発表されて以来（The International Human Genome Sequencing Consortium、2001; Venter et al、2001）、シーケンシング技術はおよそ100万倍高速かつ安価になり、個々の研究室ができるだけ多くの配列を生成することを可能にした。メタゲノミクス実験では、これは「ショットガン」シークエンシングによって微生物の複雑なmixtureをサンプリングすることを可能にし、それは単純にDNAを単離し、シークエンシングのためDNA調製し、そしてmixtureをできるだけディープにシークエンシングすることを含む。ショットガンシークエンシングは、広く使用されている16Sリボソーム RNAシークエンシングを含むターゲットシークエンシング法（Venter et al、2004）と比較して比較的バイアスがなく、真核生物を含むDNAベースのあらゆる種のゲノムDNAを捕獲するというさらなる利点を有する。偏りがないので、ショットガンシークエンシングは、各分類群に属するリード数を数えることによって、元のサンプル中の各分類群（種、属、門など）の存在量を推定するためにも使用できる。

　技術の進歩とともに、完成したドラフトゲノムの数も過去10年間で指数関数的に増加している。現在、公共のGenBank アーカイブには、数千のcompleteバクテリアゲノム、2万のドラフトバクテリアゲノム、および8万の完全または部分的なウイルスゲノムがある（Benson et al、2015）。（一部略）リードを微生物ゲノムのデータベースに迅速かつ正確にアラインメントすることができるいくつかの正確な方法が出現しているが、このステップだけではどれだけの種が存在するかを推定するのに十分ではない。同じサンプル内に closey relatedな種が存在すると、厄介な問題が発生する。多くのリードが、複数の種に等しくアライメントされるためで、これは頻繁に発生する。これを解決するには、別の存在量推定アルゴリズムが必要になる。本論文では、個々のリードを単純に分類し、メタゲノミクス実験で収集されたDNA配列から種、属、またはその他の分類学的カテゴリーの存在量を計算することを超えた新しい方法、Brackenについて説明する。

　Krakenメタゲノミクス分類器は、それが2014年に最初に発表されたとき、大きなメタゲノミクスシーケンスデータを処理することができる大きな速度向上を表し（Wood＆Salzberg、2014）、その時点で最も近い競合だったMegaBlastよりも900倍以上速く走った（Morgulis et al、2008）（kraken紹介）。 Krakenの成功と正確さは、長さkの短い、非常に効率的な短い配列のインデックス使用に依存している。これは特殊なデータベースに組み込まれる。 k値が適切に選択された場合、データベース中の長さkのほとんどの配列は単一の種に固有のものとなり、そして多くは特定の株またはゲノムに固有のものとなる。 kの値が大きいほど、各ゲノムのさらに多くがk-merによって一意的にカバーされているデータベースが得られる。しかし明らかに、kはシーケンシングリードの長さより長くてはいけない、そしてメタゲノミクスプロジェクトは現在75〜100bpほどの短いリードを生成する。また、長いk-merはエラーを含む可能性が高くなる。つまり、kが長すぎると、より多くのリードが分類されないままになる。対照的に、より小さいk-merは、最小一致長がより短いのでより高い感度をもたらすであろう。

　Krakenシステムは非常に速く正確であるが（Wood＆Salzberg、2014）、しかし、rawシーケンシングリードを分類するとき、多くのリードは2つ以上のゲノム間の同一領域に対応する。 Krakenは、この問題を解決するために、後述するように、シーケンスを共有するすべての種のlowest common ancestor（LCA）でシーケンスをラベル付けする。

微生物種と菌株のあいまいさ

　バクテリアゲノムのデータベースが大きくなるにつれて、ますます多くのゲノムが他のゲノムとそれらの配列の大部分を共有している。多くの場合、これらのゲノムはほぼ同一である。確かに、シーケンシングは研究者たちに、多くの異なった種と属がシーケンシング前に知られていたよりはるかに近いことを明らかにしている。結果として、継続的に多くの種名が変更されてきたが、他の多くの種は、しばしば歴史的な理由またはその他の理由で古い名前を保持している。

　例えば、Mycobacterium bovis種はMycobacterium tuberculosisと99.95％以上同一であり（Garnier et al、2003）、ヒト結核の多くの症例はM. tuberculosisではなくM. bovis（これもウシに感染する）によって引き起こされる（Grange、 2001）。それらの高い配列同一性は、それらが単一の種の２つの株として考慮されるべきであることを示すが、それらは異なる種名を保持する。妥協案として、分類学者たちは、カテゴリーMycobacterium tuberculosis complex（Brosch et al、2002）を作成して、現在5つの異なる種の100以上の株を含む分類群コレクションを表している。このカテゴリーは、現在の微生物分類学において種のレベルより上で属のレベルより下に位置しているが、それは種として最もよく説明できる。

　他の例は数多くあり、まだ増え続けている。炭疽菌の3つの種、Bacillus anthracis、Bacillus cereus、およびBacillus thuringiensisは、99％をはるかに超えて同一であり、すべて単一の種として命名される必要がある（Helgason et al、2000）。シーケンシングによりそれらがほぼ同一の配列であることが明らかにされたにもかかわらず。妥協案として、分類学者はこれらの3つの種と少なくとも5つの他の種を、種と属の間のカテゴリーとなるBacillus cereus groupを作成した（Liu et al、2015 link）。場合によっては、同じ種で呼ばれるべき2つの生物の属名が異なることさえある。例えば、Escherichia coliとShigella flexneriは異なる属に分類されているが、我々はそれらが同じ種を表すことをシーケンシングから知っている（Lan＆Reeves、2002 link）。

　バクテリア系統分類学の"変異性"を認識していないと、メタゲノム分類器の性能について誤った結論を導き得る。例えば、ある最近の研究（Peabody et al、2015）は11種のモックコミュニティを作成した。そのうちの1つはAnabaena variabilis ATCC 29413で、この種名は変更されていて、Nostoc属の種（Thiel et al、2014 link）と同義であることはわからなかった。Anabaenaがデータベースから削除されたとき、Krakenは正しくNostocとしてリードを識別した。しかしPeabody et al、2015は、これらすべてのリードが誤分類されていると誤って見なした。

分類 vs 存在量推定

　Krakenは、メタゲノム解析サンプルのすべてのリードにtaxonomyラベルをアサインしようとする。このデータベースには、ユーザーが選択した種をすべて含めることができる。現在の一連のcompleteバクテリアゲノムおよびアーキアゲノムの中で、それらの配列の大部分が異なる株、種、または属の他のゲノムと同一である何百もの種が見いだされる。これらの種の共通領域から生じるリードは、Krakenの分類アルゴリズムを用いて分析された場合に同じスコアとなるため、KrakenはLCAのみを正しく報告する（Wood＆Salzberg、2014）。その結果、ゲノムの多様性が低い集団のクレードについては、Krakenは固有の領域からのリードについて種レベルのアサインのみを報告し、真の総存在量は、種および属（またはそれ以上）のレベルのアサインを考慮することによってのみ可能になる。これは、いくつかの種では、大多数のリードがより高いレベルの分類で分類されるかもしれないことを意味する。このようにKrakenは種レベルを超えて多くのリードを「取り残された」ままにする。つまり、種に直接分類されるリード数は実際に存在する数よりはるかに少ないかもしれない。

　したがって、より高い分類レベルのリードを無視すると一部の種を著しく過小評価するため、Krakenのrawリードのアサインで種レベルまたは株レベルの存在量推定値に直接変換できるという仮定（Schaeffer et al。、2015）には欠陥がある。そしてKrakenのアサイン自体が間違っていたという誤った印象を作成する。

　それにもかかわらず、メタゲノミクス分析はしばしば特定のサンプル中の種の存在量を推定することを含む。リードを種に明確にアサインすることはできないが、具体的にはサンプル内のリード数または割合を推定することによって、各種の存在量を推定したい。メタゲノミクスサンプル中の種の存在量を推定するためのソフトウェアツールがすでにいくつか開発されている[MetaPhlAn, ConStrains, GAAS, GASiC, TAEC, GRAMMy] (Angly et al., 2009; Lindner & Renard, 2012; Luo et al., 2015; Segata et al., 2012; Sohn et al., 2014; Xia et al., 2011)。しかし、これらのツールは、Krakenのk-merアプローチほど正確で効率的とは限らない、リードレベルの分類に異なる戦略を採用している（Lindgreen、Adair＆Gardner、2016 link）。あいまいなリード分類問題に対処し、存在量推定値を直接提供するようにKrakenを再設計するのではなく、ここで説明する新しい種レベルの存在量推定方法を別のプログラムとして実装した。これにより、既存のKrakenユーザーとの後方互換性が維持され、Krakenによって既に処理されたデータセットに対してより正確な種の存在量の推定値を生成することができる。 Krakenが種の識別に失敗した場合（例えば、その種がKrakenデータベースから欠落していた場合）、Brackenもその種を識別しない。

　著者らの新しい方法であるBracken（Bayesian Reestimation of Abundance after Classification with KrakEN）は、taxonomic tree内のリードを確率的に再分布させることによってメタゲノミクスサンプル中の種の存在量を推定する。種レベルより上のノードにアサインされたリードは、種ノードの下に分散されるが、株レベルでアサインされたリードは、それらの親の種上に再分散される。例えば、論文図1では、Mycobacteriaceae科とMycobacterium属にアサインされたリードをM. marinumとM. aviumに分配し、各M. marinum株にアサインされたリードはM. marinum種に再アサインされる。以下に示すように、Brackenは同じアルゴリズムを使用して他の分類レベル（例えば属や門）で存在量を簡単に再評価することができる。

インストール

Bracken is a companion program to Kraken 1.0 or Kraken 2.0.

本体　GIthub

git clone https://github.com/jenniferlu717/Bracken.git
cd Bracken/
sh install_bracken.sh

#追記　conda
mamba install -c bioconda bracken -y

> ./bracken -h

$ ./bracken -h

./bracken: illegal option -- h

Usage: bracken -d MY_DB -i INPUT -o OUTPUT -r READ_LEN -l LEVEL -t THRESHOLD

MY_DB location of Kraken database

INPUT Kraken REPORT file to use for abundance estimation

OUTPUT file name for Bracken default output

READ_LEN read length to get all classifications for (default: 100)

LEVEL level to estimate abundance at [options: D,P,C,O,F,G,S] (default: S)

THRESHOLD number of reads required PRIOR to abundance estimation to perform reestimation (default: 0)

> ./bracken-build -h

$ ./bracken-build -h

./bracken-build: illegal option -- h

Usage: bracken_build -k KMER_LEN -l READ_LEN -d MY_DB -x K_INSTALLATION -t THREADS

KMER_LEN kmer length used to build the kraken database (default: 35)

THREADS the number of threads to use when running kraken classification and the bracken scripts

READ_LEN read length to get all classifications for (default: 100)

MY_DB location of Kraken database

K_INSTALLATION location of the installed kraken/kraken-build scripts (default assumes scripts can be run from the user path)

**Note that this script will try to use kraken2 as default. If kraken2 is not installed, kraken will be used instead

データベース

kraken2データベースを作ったら、これを拡張してbrackenデータベースも構築する。40スレッド、k-mer長35（kraken2と同じ）、read length 100。

bracken-build -d KRAKEN2 -t 40 -k 35 -l 100

kraken2 buildの--cleanで不要なファイルを削除してしまうと実行できないので注意する。

実行方法

krakenレポートファイルを指定する。

bracken -d KRAKEN_DB -i kraken_report -r 100 -l S -o bracken

-i location of the built Kraken 1.0 or Kraken 2.0 database
-r read length
-t Default = 10. For species classification, any species with <= 10 (or otherwise specified) reads will not receive any additional reads from higher taxonomy levels when distributing reads for abundance estimation. If another classification level is specified, thresholding will occur at that level.
-l Default = 'S'. This specifies that abundance estimation will calculate estimated reads for each species. Other possible options are K (kingdom level), P (phylum), C (class), O (order), F (family), and G (genus).

追記

複数サンプルの結果をマージしてTSV出力（*1）

combine_bracken_outputs.py --files *_bracken -o output.tsv

condaでBrackenを導入していれば、combine_bracken_outputs.pyのパスは通っている。

引用

Bracken: estimating species abundance in metagenomics data

Jennifer Lu, Florian P. Breitwieser, Peter Thielen, Steven L. Salzberg

January 2, 2017 PeerJ

https://peerj.com/articles/cs-104/

参考

https://www.microbe.net/2017/04/27/why-use-bracken-instead-of-kraken/

macでインフォマティクス

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

種の存在量を推定する Bracken