RDP Classifier の原理 - ほくそ笑む

1. はじめに

RDP Classifier は、RDP(Ribosomal Database Project) により開発された、16S rRNA 配列から菌種を判別するためのツールです。

RDP Classifier は、BLAST に比べて数百倍の速度で菌の判別(簡易同定)を行うことが可能です。また、この論文によると、RDP Classifier は misclassify はありますが、BLAST と同程度の正確性で菌の判別を行うことができます。

RDP Classifier は、菌の判定にナイーブベイズ分類器を使用します。今回は、RDP Classifier がどのように菌の判定を行っているかを説明してみることにします。

2. ナイーブベイズ分類器

菌の判定を行いたい 16S rRNA 配列を $S$ 、genus を $G$ とすると、ベイズの定理より、
　 $P(G|S) = \frac{P(S|G)P(G)}{P(S)}$
が成り立ちます。

$P(G|S)$ は、 $S$ が与えられた時に $G$ であるという条件付き確率を表します。したがって、 $P(G|S)$ をすべての $G$ に対して求めたとき、最も大きい確率となる $G$ が、 $S$ に対する最も可能性の高い genus であると判断できます。

今、 $W=\{w_1, w_2, \cdots, w_d\}$ を 8 塩基からなるすべての配列の集合とします*1。 $W$ の要素数は $d = \left|\{A,T,G,C\}\right|^8=4^8 = 65536$ となります。 $S$ に対して、ベクトル $w_S = (w_{S1}, w_{S2}, \cdots, w_{Sd})$ を次のように定義します。
　 $\begin{eqnarray}w_{Si} = \left\{\begin{array}{ll}1, & w_i \rm{ is subsequence of } S\\0, & w_i \rm{ is not subsequence of } S\end{array}\right.\end{eqnarray}$
このとき、
　 $P(S|G) = P(w_S|G)$
と仮定します*2。ナイーブベイズの独立性仮定により、
　 $P(S|G) = P(w_S|G) = \prod_{i=1}^d P(w_{Si}|G)$
が成り立ちます。

ここで、学習データを用意します。学習データは、あらかじめ分かっている genus と 16S rRNA 配列との対応を記したものです。
学習データ中において、 $G$ に対応する配列数を $M_G$ 、その中で $w_i$ を部分文字列として含む配列数を $m_G(w_i)$ とおくと、
　 $P(w_{Si}=1|G) = \frac{m_G(w_i)}{M_G}$
となります。しかし、このままではゼロ頻度問題が起こるのでスムージングを行います。
　 $\begin{eqnarray} P(w_{Si}=1|G) &=& \frac{m_G(w_i) + P_i}{M_G + 1} \\ P_i &=& \frac{n(w_i) + 0.5}{N + 1}\end{eqnarray}$
ここで、 $N$ は学習データの配列数、 $n(w_i)$ は学習データ中で $w_i$ を部分文字列として含む配列数を表します。
ここで行われているスムージングは単純な加算スムージングではないようです。不勉強な私には分かりませんが、 $0 < P(w_i|G) < 1$ となることは分かります。

一方、
　 $P(w_{Si}=0|G) = 1 - P(w_{Si}=1|G)$
となります。

以上により、具体的に $P(S|G)$ を計算することができます。
すべての $G$ に対して $P(S|G)$ を計算したとき、最大となる $G$ が求めるべき genus です。

3. ブートストラップ法

上記により、16S rRNA 配列が与えられたとき、ナイーブベイズ分類器により genus を判定することができるようになりました。しかし、この判定結果がどのくらい信頼のできるものかはわかりません。

ナイーブベイズ分類器は、もし未知の菌が入力されたとしても、必ず既知の菌から判定結果を出してしまいます。したがって一定の信頼性を得られない判定結果は、判別不能として扱うべきです。このような要求に対して、RDP Classifier は、判定結果の信頼性をブートストラップ法で算出します。

今、16S rRNA 配列 $S$ に対して、集合 $V_S$ を、 $S$ の部分文字列となる $W$ の要素の集合とします( $V_S \subseteq W$ )。この集合 $V_S$ から復元的リサンプリングを行うことで、ブートストラップサンプル $V_S^*$ を作成します。すなわち、 $V_S$ からランダムに要素を取り出し、新たな集合 $V_S^*$ を作るのです。(このとき同じ要素を何度取り出してもよい)

このブートストラップサンプルの意味するところは、 $S$ の単語セットを抽出することに対応します。 $S$ は単なる配列であり、単語の抽出に使用できる文法が無いため、上記ナイーブベイズ分類器の説明では、すべての部分文字列を単語としていました。しかし、すべての部分文字列を単語とすると、配列上で単語同士に重複が生じるため、ナイーブベイズの独立性仮定を満たさなくなってしまいます。

したがって、重複を生じないように単語を選ぶ必要があります。ブートストラップサンプルは、配列から単語をランダムに選んできたものであり、抽出回数を、配列から重複なく得ることのできる単語数に制限すれば、独立性仮定を満たした単語セットとみなせます。配列 $S$ から重複なく得ることのできる単語数は、 $S$ の長さを $n$ としたとき、 $n/8$ となります。( $w_i$ の塩基数は 8 なので)

このようにして作成したブートストラップサンプル $V_S^*$ を使って、genus を判定します。上記で行った説明と同様に、ベクトル $w_S^*=(w_{S1}^*, w_{S2}^*, \cdots, w_{Sd}^*)$ を下記のように定義します。
　 $\begin{eqnarray} w_{Si}^* = \left\{ \begin{array}{ll} 1, & w_i \in V_S^* \\ 0, & w_i \notin V_S^* \end{array} \right. \end{eqnarray}$
この $w_S^*$ を $w_S$ の代わりに使用すれば、ナイーブベイズ分類器による genus の判定を行うことができます。

ただし、ブートストラップサンプルはランダムに抽出したものであり、 $S$ を正確に特徴づけるものではありません。そこで、ブートストラップサンプルを 100 回作成し、genus の判別を行います。100 回中何回その genus に判別されたかをカウントすれば、多いものほど $S$ に対する genus の判別結果として信用できることになります。
この 100 回中何回判別されたかの割合が、その判別結果の信頼度として使えるということです。

genus より上の階層での判別は、ブートストラップで 100 個の genus の判別結果を出したあと、それらの親階層を GenBank Taxonomy を使ってたどり、genus と同様に多数決で決めます。*3