DNAメタバーコーディングに頻繁に用いられるQiime2では、ナイーブベイズ分類器を訓練して生物情報の割り当てを行うことができます。また、16S RNAなどメジャーな領域では既にトレーニングされた訓練器を用いることができます。ここでは導入としてナイーブベイズ分類器について簡単な説明をしたいと思います。
___________________
1)前提(ベイズの定理)
①ベイズの定理
ナイーブ分類器にはベイズの定理が用いられています。
\( P(A|B) = \displaystyle\frac{P(B|A) \times P(A)}{P(B)} \)
ベイズの定理は、あるデータがある場合に、それが各クラスに属する確率を計算するものです。例えば、あるメールがスパムメールであるかそうでないかを分けるような場合に用いられ、この場合はスパムメールである確率と、スパムメールではない確率が計算されます。
③事前確率と事後確率
ある病気の検査を考えてみましょう。この病気は全人口の1%が罹患しているとしましょう。しかし、特定の年齢層、たとえば60歳以上の人々の中では、病気の発生率が3%という新しい情報が得られました。この新しい情報をもとに計算した確率が事後確率です。事後確率の前の段階である何の情報もない状態(この場合1%)を事前確率といいます。
では実際にナイーブベイズ分類器の動作概要(アルゴリズム)を見ていきましょう。
2)事前準備
①事前確率の準備
塩基配列を元にタクソノミー情報を割り当てる場合を想定してみます。ある様々な生物の塩基配列情報とその生物のタクソノミー情報がセットになったデータベースがあるとします。この塩基配列とタクソノミー情報(どのクラスに割り当てるかセットになったすべてのデータの数に対して、ある生物のデータがいくつあるかその割合を訓練データで事前に計算しておきます。もちろんその割合は生物毎に存在します。これが事前確率となります。
②尤度の準備
訓練データで生物のクラス毎に特徴量が現れる確率を計算しておきます。Qiime2のq2-feature-classifierプラグインのデフォルト設定では連続する7つの塩基対(7-gram)をセットにして特徴量としています。これは--p-feat-ext--ngram-range オプションで確認することができます。例えば以下のようになります。
シーケンス: "ATCGAATC"
2-gram: "AT", "TC", "CG", "GA", "AA", "AT", "TC"
3-gram: "ATC", "TCG", "CGA", "GAA", "AAT", "ATC"
3)分類
ある生物の塩基配列情報があり、それが属するクラスを推定します。
基本式になるのはベイズの定理です。
\( P(X|A) = \displaystyle\frac{P(A|X) \times P(X)}{P(A)} \)
Aという塩基配列があった場合にXという生物に分類される確率は?ということが求めたいことです。これを各生物分類において求め、最も確率が大きいものをその生物種であるとして採用します。
①まず分母ですが、その塩基配列Aである確率です。これは無視できるのですが、一応置いておきます。
②次に分子のP(X)ですがこれは生物Xになる確率です。これは先ほど事前確率で求めた手順で求めます。
③さらに分子のP(A|X)を求めます。(「Aという塩基配列情報があった場合にXというクラスにそれが該当する可能性は?」と逆になっていることに注意してください。
④このような要領で各カテゴリとなる確率を求めていきます。
(具体例)
\( P(X|A) = \displaystyle\frac{P(A|X) \times P(X)}{P(A)} \)
\( P(Y|A) = \displaystyle\frac{P(A|Y) \times P(Y)}{P(A)} \)
\( P(Z|A) = \displaystyle\frac{P(A|Z) \times P(Z)}{P(A)} \)
これらの式を見て頂くと分母は共通しているので、大小比較をするなら分子だけに注目すれば良いことが分かります。
P(X)は訓練データからその生物分類の数が全体の生物分類数に占める割合となります。
次にもう一つの分子P(A|X)を求めます。少々やっかいです。P(A|X)は生物分類が既知である生物Xの配列がその特徴量となっている確率です。
生物Xの配列は、先ほどの7-gramで捉えた特徴量の集合として考えることができます。この特徴量が生物Xの配列情報Aに対してN個含まれているとします。
さらに配列情報Aのそれぞれの7-gramに名前をつけていなかったのでg1, g2 ,g3 という名前にします。そして
\[ P(A|X) = P(g1|X) \times P(g2|X) \times P(g3|X) \]
というように掛け算で求めていきます。例えばP(g1|X)は生物Xにg1という7-gramが含まれる確率です。
実際の計算では計算しやすいように両辺でlogをとります。logを取ると右辺が足し算になり、計算しやすくなります。
ここでg1、g2、g3...g3はそれぞれ独立に生じるとみなします。現実はg1が出たときはg2が出やすくなるなど各7-gramは独立ではなく正しくないのですが便宜上そのように仮定します。カテゴリ分類であればこれでうまくいきます。これが「ナイーブ(単純)ベイズ」と呼ばれる背景です。
P(g1|X)は
分母に訓練データに存在する生物Xのすべての配列の7-gramの数とし、分子にg1が出てきた回数とします。但し、g1が出てこない場合、掛け算で求めれられる
P(A|X)がゼロになってしまうので、こちらも便宜上
分母=(生物Xのすべての配列の7-gramの数+1)
分子=(g1が出てきた回数+1)
とします。
これらのソースコードはRやPythonで調べると出てきますので、試してみるとよいと思います。
4)メリットとデメリット
紹介したように比較的簡易に計算が可能で、高速に計算が可能です。その割に比較的良い結果が得られるとされています。一方で、現実にはありえませんが全ての特徴量が独立として仮定しています。さらに事前確率を仮定するのですが、データセットの偏りがあると出力結果も偏りがある方に働きます。DNAメタバーコーディングの結果は精査が必要と考えておいた方が良いと思われます。
【追記】
feature-classifier fit-classifier-naive-bayesには事前学習を行わないオプションを使うことができます。データセットに偏りがあり、その影響を回避したい場合はこのオプションを利用して学習モデルを作成することが考えられます。。
具体的には --p-classify--fit-prior オプションをFalse にします。この場合すべてのクラスが同じ確率で発生するという均等な仮定の下でモデルがトレーニングされます。デフォルトは[null]になっており、事前確率を自動的に推定します。
以下はFalseオプションを使った場合のコマンド例です。
qiime feature-classifier fit-classifier-naive-bayes \
--i-reference-reads NCBI_ITS2_Viridiplantae_fasta_file_2023_10_17.qza \
--i-reference-taxonomy NCBI_ITS2_Viridiplantae_taxonomic_lineages_2023_10_17.qza \
--o-classifier NCBI_ITS2_Viridiplantae_classifier_2023_10_17.qza \
--p-classify--fit-prior False