グループテストにもとづく事後確率の近似計算アルゴリズム

(1)

グループテストにもとづく事後確率の近似計算アルゴリズム

2008 年度

上原啓明

(2)

概要

1940

年代に

Dorfman

は，多数の検体から少数の陽性反応を示す検体

を効率よく選び出すには，グループテストを用いることが有効であることを示した．グループテストでは，すべての検体を個々に検査するのではなく，複数の検体をまとめていくつかのグループを作り，各グループに対して検査を行う．検査結果に誤りがないと仮定すると，あるグループに対する検査で陽性反応が検出されない場合は，

1

回のテストでそのグループに含まれる検体がすべて陰性であると判断できる．一方，陽性反応が検出された場合にはそのグループに含まれる検体のうち少なくとも

1

個の検体が陽性

(

ポジティブ

)

であることがわかる．また，背景の同じ検体が複数のグループに属するようにグループを作り，各グループに対する検査結果から陽性である可能性が高い，つまり，事後確率の高い検体に対して

2

次検査で個々に陽性反応を示すかどうかを検査し，ポジティブである検体を識別する．すなわち，グループテストは実行するテストの回数を少なくすることにより実験にかかるコストを減らし，かつポジティブである検体をできるだけ多く確実に発見することを目的として行われる．特に陽性の検体の割合が非常に少ない場合に大幅にテストの回数を減らすことができる．しかし実際には，テストの結果に誤りが生じることも少なくない．

また，各検体が陽性である事後確率を正確に求めるには検体の数に対して計算量が指数関数的に増加する．したがって，グループテストの結果に偽ポジティブ・偽ネガティブの誤りが生じる場合に検査結果からポジティブである可能性の高い検体を効率良く識別するアルゴリズムの研究が重要である．さらに，複数のグループが成す集合族の組合せ構造によって，ポジティブの識別能力および選別

(

)

能力が異なり，そのための最適な組合せ構造の研究も必要である．

グループテストは情報科学や生命科学などの分野などにも広く応用され

(5)

ており，

DNA

の塩基配列

(A

，

T

，

G

，

C

の配列

)

を特定する

DNA

ライブラリスクリーニングにも応用されている

(Sham

ら

(2002))

．

DNA

ライブラリスクリーニングの分野では，各検体をクローンと呼び，グループをプールと呼ぶ．また，さまざまなプールの集まりをプーリングデザインと呼んでいる．

Knill

ら

(1996)

は

DNA

ライブラリスクリーニングで必要となる，陽性である事後確率を効率的に求める

1

つの方法として

Markov chain pool result decoder (MCPD)

と称すアルゴリズムを提案している．さらに，グループテストは

DNA

マッピングにも用いられている

(Du

ら

(2000))

．近年では，

DNA

チップやマイクロアレイによる試験にも応用されている

(Schliep

ら

(2003)

，

Manoli

ら

(2006))

．

本論文では，グループテストの結果にもとづいて高速かつ高精度に各検体がポジティブである事後確率を計算する

2

つのアルゴリズムを提案し，

その能力を

MCPD

と比較する．

第

1

章では，グループテストに関連する研究の歴史や動向について概説し，特に

DNA

ライブラリスクリーニングに関して概説し，グループテストとその背景にある確率モデルについて述べる．

第

2

章では，まず

Pearl (1988)

により提案されたベイジアンネットワークの確率伝播アルゴリズムにもとづき考案した

Bayesian network pool re-

sult decoder (BNPD)

アルゴリズムを提案する

(Uehara and Jimbo (2007))

．確率伝播法は確率変数間の依存関係を表すグラフ上でいくつかの頂点での実現値が与えられたもとでの事後確率を効率良く計算するアルゴリズムである．グループテストのグループと検体の結合関係はタナーグラフと呼ばれる

2

部グラフで表される．タナーグラフに閉路がない，つまり木である場合には，

BNPD

は正確な事後確率を計算することができるが，タナーグラフに閉路が存在しても短い閉路が多数存在しない場合にはアルゴリズムは収束し事後確率の近似値を得ることができる．しかし，閉路が存在する場合には

BNPD

により求めた周辺事後確率と真の値との間には偏りが存在する．さらに，短い閉路が多数存在する場合にはアルゴリズムが収束しない恐れがある．この収束性に関する欠点を補うために

concave-convex pro-

cedure (CCCP)

にもとづき考案した

CCCP pool result decoder (CCPD)

アルゴリズムを提案する

(Uehara and Jimbo (2008))

．

CCCP

では，統

(6)

計力学における

Bethe

自由エネルギーをある制約条件のもとでラグランジュの未定乗数法を用いて最小化する問題として定式化することにより，

事後確率の推定値を求めている．タナーグラフに短い閉路が多数存在する場合でも

CCPD

は収束が保証されているアルゴリズムである．しかし，

CCPD

は

BNPD

に比べて計算時間が長くなる．また，

BNPD

と同様に

CCPD

により求めた周辺事後確率には偏りが存在する．

BNPD

および

CCPD

の欠点である求めた周辺事後確率の偏りについては，その偏りに関する

1

つの指標を導き，その結果を用いてタナーグラフに長さ

4

の閉路が存在しない場合にはその

2

次項が

0

になり，偏りが小さくなることを理論的に示す．

そして，第

3

章では

BNPD

および

CCPD

と

MCPD

とをさまざまな検体の数に対するタナーグラフのもとでシミュレーションにより比較し，

BNPD

，

CCPD

の高速性とスクリーニング効率を実証する．タナーグラフに長さ

4

の閉路が存在しない，つまり，どの

2

つの検体も同時に高々

1

個のグループにしか含まれないという条件を満たすとき，

unique collinearity condition

を満たすという．

unique collinearity condition

を満たす場合に，

BNPD

のスクリーニング効率の良さに言及する．また，

unique collinearity condition

を満たさない場合には，

BNPD

が収束しないことが多いためにその代替アルゴリズムとしての

CCPD

と

MCPD

のスクリーニング効率を比較する．次に，検体数とグループ数，およびタナーグラフの辺の数を固定したとき，タナーグラフが正則かつ

unique collinearity condition

を満たす場合，ランダムにタナーグラフを生成する場合よりスクリーニング能力が高くなることを実証する．さらに，

unique collinearity condition

を満たさない場合に

CCPD

により求めた周辺事後確率の偏り補正の効果を検証する．

最後の第

4

章は，結論として，タナーグラフが

unique collinearity con-

dition

を満たすか否かによって，

BNPD

，

CCPD

，

MCPD

を使い分けることが有用であることを述べる．

(7)

第 1 _章

グループテスト

本章では，まずグループテストの歴史や研究の動向およびその確率モデルなどについて述べる．グループテストは多数の検体から少数の目的物を効率良く識別するためにさまざまな分野で用いられているが，本論文ではグループテストの応用例として，

DNA

ライブラリスクリーニングを挙げ，その選別

(

)

効率を高めるためのグループテストを中心に述べる．グループテストのためのアルゴリズムには組合せ論的なアルゴリズムと確率論的なアルゴリズムがある．ここでは主に本論文で扱う確率論的な事後確率計算アルゴリズムのための確率モデルについて述べる．

1.1 _{グループテストの歴史}

グループテストは

1940

年代に

Robert Dorfman [15]

により血液検査にかかるコストを低く抑える目的で考案された方法である．グループテストを用いることにより実験や検査に要する費用や時間，作業量を削減することができる．特に陽性の検体の割合が非常に少ない場合に大幅にコストを減らすことができる．例えば，グループテストは罹患率の推定

(Thompson [62]

，

Sobel

ら

[57]

，

Tu

ら

[63]

，

Brookmeyer [9])

や患者対照研究における病原体や毒素の汚染の評価

(Weinberg

ら

[67])

のために用いられてきた．また，

HIV

の検査の際に匿名性を守るためにも用いられてきた

(Gastwirth

ら

[25]

，

Gastwirth

ら

[26])

．遺伝学の分野では，グループテストは

I

型

(

インスリン依存型

)

糖尿病の研究で初めて用いられた

(Arnheim

ら

[2])

．その後，植物の連鎖研究

(Michelmore

ら

[45])

や近交

(8)

集団の劣性疾患の同型接合性マッピング

(Sheﬃeld

ら

[55]

，

Carmi

ら

[12]

，

Nystuen

ら

[49]

，

Scott

ら

[53])

や変異検出

(Amos

ら

[1])

に用いられている．近年ではヒトゲノム計画の成功に見られるように，グループテストは

DNA

ライブラリスクリーニングのような分野での新しい応用が見出されている．

1.1.1 グループテストと DNA ライブラリスクリーニング

本論文では，グループテストの例として

DNA

ライブラリスクリーニングを取り上げる．まず，

DNA

ライブラリスクリーニングの実験の概要とそこにグループテストがどのように用いられるのかを説明する．

DNA

は

A (Adenine)

，

T (Thymine)

，

G (Guanine)

，

C (Cytosine)

の

4

種類の塩基から構成されている．

DNA

の塩基列の断片の複製をクローンと呼ぶ．

DNA

ライブラリとはある特定の生物の

DNA

を短い塩基列の断片に切り，各塩基列ごとにクローンのプールを作って得られる単一種のクローンのプールの集合である．

例えば，

Knill

ら

[35]

がグループテストによるスクリーニングを行った

Cla I

ライブラリと呼ばれる

DNA

ライブラリはヒトの第

16

染色体から

生成され，

1298

種のクローンからなる

(McCormick

ら

[43])

．

DNA

の列の長さは

98 Mbp (Mega base pair)

で複数の

DNA

の列をショットガンなどの方法で短いランダムな長さの塩基列に細分し，それらの塩基列をバクテリアに組み込んで増殖させ，単一種のクローンのコロニーを作る．それぞれの種のクローンを

microtiter plate

の

1

つの

well

に保存する．このような過程を経て，

Cla I

ライブラリの

1298

種のクローンが得られているが各クローンの平均長は

195 kbp (kilo base pair)

程度である

(Morton [46])

．複数の

DNA

の列を細分してクローンが作られているため

DNA

の同じ部位が複数種のクローンに重複して存在する．その平均重複度はこの場合には

195 kbp × 1298 / 98 Mbp = 2.6

程度である．これらのクローンに含まれる遺伝子などの塩基列を特定するために

PCR

法などを用いて目的のクローンのみを増幅させ反応試験を行う．この際にプローブと呼ばれる目的の塩基列に蛍光物質を添付しておき，プローブと相補関

(9)

係にあるクローンと次々にハイブリダイズと分離を繰り返すことにより，

目的のクローンのみがコピーされる．これにより目的のクローンがプールに含まれている場合にはプールが強い蛍光を発する．この蛍光の強さを測定することにより，目的のクローンの有無を識別することができる．

遺伝子機能の研究が非常に重要な研究分野に発展したことにより，高品質な遺伝子のライブラリが必要とされている．

DNA

ライブラリのクローンに含まれる遺伝子などを識別し，そのタグを付けてタグ付きの

DNA

ライブラリが得られる．これらのライブラリは実験を行う研究施設などに保存されマイクロアレイの解析のような多くの遺伝子情報解析のために繰り返し利用される．したがって，誤りのないタグが付いた品質の良いライブラリ生成への要求は強く，そのための高い識別能力を持つ効率の良いスクリーニング手法の開発は非常に重要である．

スクリーニング実験ではクローンがある特定の

DNA

の塩基列を含んでいるか否かを判定する．クローンがある特定の

DNA

の配列を含んでいるときそのクローンはポジティブであるといい，そうでないときネガティブであるという．また，ポジティブクローンを単にポジティブと呼ぶこともある．

DNA

ライブラリスクリーニングの目的は大量のクローンの中からポジティブクローンを識別すること，つまり，クローンのライブラリがある特定の塩基列を含むか否かを識別することである．その際に時間とコストを削減することが必要であり，そのためにグループテストを用いることが有効である．複数種のクローンからなるグループをプールと呼び，このプールに対して試験を行う．さまざまなクローンの組合せからなるプールの集まりをプーリングデザインと呼ぶ．もし各プールに対する試験結果に誤りがないとすると，あるプールに対する試験結果がネガティブならば，

そのプールに含まれるすべてのクローンがネガティブであることが一回の検査で判定でき，個々のクローンに対して試験するより実験回数を減らすことができる．一方，もしプールの検査結果がポジティブならそのプールに少なくとも一つのポジティブクローンが含まれていることがわかる．同じ断片から複製されたクローンが複数のプールに属するようにプールを作り，各プールに対して検査を行う方法をグループテストと呼ぶ．グループテストの結果からポジティブクローンを特定するスクリーニングを行う．

(10)

上記のようなグループテストはポジティブクローンの割合が非常に小さい場合にはプールに対する試験結果がネガティブとなる確率が高く，実験回数を減らすことが可能となる．グループテストの効率の良さは

Barillot

ら

[5]

，

Berger

ら

[6], Bruno

ら

[10]

，

Du

ら

[16, 17, 18]

，および

Sham

ら

[54]

で述べられている．

グループテストの分類の仕方のひとつに

adaptive

と

non-adaptive

という分類がある．

adaptive

なグループテストでは，それまでのグループテストの観測値の情報を得てから次のプールを作成することができる．つまり，試験の結果の観測値によってプールを順次作成することができる．

したがって，試験とプールの作成という手順を繰り返すことによってポジティブクローンを識別する．一方，

non-adaptive

なグループテストでは，

実験を行う前にすべてのプールを決定しておく．つまり，試験の結果の観測値によってプールを順次作成することができない．このグループテストの方法ではすべての試験を同時に行い，プールに対する観測値を解析しポジティブクローンを識別する．この

2

つの方法のうちどちらを用いるかは状況による．

adaptive

なグループテストは

non-adaptive

なグループテストに比べてプールを作成するための情報が多いので，一般にテストの回数を少なくすることができる．一方，実験と解析が異なる人や組織によって行われるときは

adaptive

な方法を用いることが難しく，同時にすべてのテストを行う

non-adaptive

な方法を用いる．本論文では，

non-adaptive

な場合を考える．

通常，観測値には偽ポジティブや偽ネガティブのような誤りが起こることは避けられない．偽ポジティブとはネガティブなクローンあるいはプールがポジティブと判断される誤りであり，偽ネガティブとはポジティブなクローンあるいはプールがネガティブと判断される誤りである．このような試験結果の観測値の誤りが無視できない状況のもとで多量のクローンの中からポジティブクローンを識別するための効率の良いスクリーニングアルゴリズムの開発とプーリングデザインの計画が重要である．本論文では，プーリング実験において高い精度でポジティブクローンを識別するために次の

2

つの問題に焦点をあわせる．

(i)

効率的なポジティブ識別アルゴリズム

(11)

クローン

プール

c

₁

c

₂

c

₃

c

₄

c

₅

c

₆

· · ·

c

_n

G

₁

G

₂

· · ·

G

_m

図

1.1:

プーリングデザインおよび

(ii)

効率的なプーリングデザインの組合せ論的な性質．

なお，効率の良いプーリングデザインの組合せ論的な性質は識別アルゴリズムに依存する可能性があることに注意する必要がある．

1.1.2 DNA ライブラリスクリーニングのためのプーリングデザインと 2 段階グループテスト

C

を

n

個のクローンの集合，

G

^を

m

個のプールの族とし，クローンを

c ∈ C

でプールを

G ∈ G

で表すことにする．このとき，組

(C, G )

をプーリングデザインと呼ぶ．プーリングデザインは

m × n { 0, 1 } -

行列または図

1.1

のような

2

部グラフで表すことができる．

H = (h

ij

)

を

2

元

m × n

行列とする．

H

の各行は

m

個のプールに対応しており，

H

の各列は

n

個のクローンに対応している．

H

の各成分

h

ij

を

h

_ij

=

 



 

1, c

j

∈ G

i のとき，

0, c

_j

∈ / G

_i のとき．

(12)

c

₁

c

₂

c

₃

c

₄

c

₅

c

₆

c

₇

c

₈

c

₉

G

₁

G

₂

G

₃

G

₄

G

₅

G

₆

図

1.2:

式

(1.1)

の結合行列

H

に対するタナーグラフ

と定義する．このとき，

H

をプーリングデザイン

(C, G )

の結合行列またはテスト行列と呼ぶ．

E = { (c, G) | c ∈ C, G ∈ G

^かつ

c ∈ G }

^{とすると，}

E

は

2

部グラフ

(C, G ; E)

の辺の集合とみなすことができる．このような

2

部グラフをタナーグラフ

(Tanner graph)

と呼ぶ．この「タナーグラフ」という用語は符号理論の分野で用いられている．ここでは，次章で説明する

BNPD

アルゴリズムと「

low density parity check

」

(LDPC)

の復号アルゴリズムを対応させるために同じ用語を用いることにする．

例えば

c

₁

c

₂

c

₃

c

₄

c

₅

c

₆

c

₇

c

₈

c

₉

H =

G

₁



 

 

0 1 1 0 0 1 1 0 0 

 

 

G

₂

0 0 0 1 0 0 0 1 1

G

₃

0 0 0 0 1 0 1 0 0

G

₄

0 1 0 0 0 0 1 1 0

G

5

1 0 1 0 1 0 0 0 1

G

6

1 0 0 1 0 1 0 1 0

(1.1)

とする．この場合，プール

G

₁ は

4

個のクローン

c

₂

, c

₃

, c

₆

, c

₇ を含み，

G

₅ は

4

個のクローン

c

₁

, c

₃

, c

₅

, c

₉を含んでいる．式

(1.1)

の結合行列

H

に対応するタナーグラフは図

1.2

で表される．

通常，グループテストでは

2

段階グループテストあるいは多段階グループテストが行われることが多い．

2

段階グループテストでは図

1.3

のよう

(13)

グループテスト実験

c

1

c

2

c

3

· · · c

n

G

1

s

1

G

2

s

2

· · ·

G

m

s

m

プーリングデザイン

スクリーニングアルゴリズムアルゴリズム

個別テスト

観測値

s = (s

1

, . . . , s

m

)

を出力

ポジティブクローンの候補を出力

図

1.3:

グループテストの流れ

にまず各プール

G

_i に対して反応試験が行われ，その観測値

s

_i が得られる．得られた観測値

s

₁

, . . . , s

_m からスクリーニングアルゴリズムにより，

ポジティブの確率が高いクローンを抽出し，第

2

段階の個別テストに渡す．個別テストでは各クローンごとに反応試験を行い，最終的にポジティブクローンを識別する．本論文では，主にスクリーニングのための事後確率計算アルゴリズムに焦点をあてて効率の良い事後確率計算アルゴリズムを開発する．また，事後確率計算アルゴリズムに適したプーリングデザインにも注目し，効率の良いプーリングデザインの組合せ的性質を明らかにする．

1.1.3 グループテストの結果にもとづく決定論的なスクリーニング

組合せ論的グループテストでは

n

種のクローンのうちポジティブの数が事前に固定されているかある固定された正整数

d

以下であるという仮定がなされており，決定論的なモデルが用いられる．さらに観測値に誤りがある場合には，誤りの最大数が与えられている．組合せ論的グループテストの主

(14)

な問題は，ポジティブの最大数を仮定したもとでポジティブクローンとネガティブクローンを識別できるように効率的なプーリングデザインを構成することである．組合せ論的グループテストでは以下のような「

d-disjunct

」な行列がよく用いられる．

{ 0, 1 } -

行列を

H = (h

₁

, . . . , h

_n

)

とし，その添え字の集合を

I = { 1, . . . , n }

とする．このとき，任意の

D ⊂ I

，

| D | ≤ d

と任意の

j / ∈ D

に対して

∨

i∈D

h

_i

≯ h

_j となるとき，行列

H

は

d-disjunct

であるという．ただし，

2

つの

{ 0, 1 } -

ベクトル

h = (h

₁

, . . . , h

_m

)

^⊤，

h

^′

= (h

^′₁

, . . . , h

^′_m

)

^⊤ に対して，

h ∨ h

^′

= (h

₁

∨ h

^′₁

, . . . , h

_m

∨ h

^′_m

)

^⊤ であり，任意の

j

に対して

h

_j

> h

^′_j のとき

h > h

^′ と書く．

d-disjunct

な行列は

d

個までのポジティブクローンを識別できるだけでなくその識別アルゴリズムも単純である．すなわち，あるクローンがポジティブであるのはそのクローンを含むすべてのプールがポジティブであるときまたそのときに限るということである。したがって，試験結果の観測値に誤りがないときは

d-disjunct

な行列に対するポジティブクローンの識別アルゴリズムでは，

各クローン

c ∈ C

に対して

{ G ∈ G : G ∋ c }

のプールがすべてポジティブのとき

c

はポジティブであり，ひとつでもネガティブなプールがあれば

c

はネガティブである．

d-disjunct

な行列は

Kautz

ら

[33]

が提案した強さ

d

の「

superimposed code

」の結合行列と同値である．極値集合論の分野では「

d-cover-free fam- ily

」という名前で

Erd˝ os

ら

[21]

により研究された．

さらに，

d-disjunct

な行列

H

のどの

d

列の成分ごとの論理和ベクトルもほかのどの列の少なくとも

e + 1

個の

1

を

cover

しないとき，

H

は

d

^e

-disjunct

であるという．

d

^e

-disjunct

な行列は

Macula [40]

により導入された．

d-disjunct

な行列は

d

⁰

-disjunct

である．

d

^e

-disjunct

な行列は

e − 1

個の偽ポジティブや偽ネガティブの誤りを検出することができ，

⌊ e/2 ⌋

個の誤りを訂正することができる．

d

^e

-disjunct

な行列は強さ

d

，距離

e + 1

の

superimposed code

の結合行列である．

D’yachkov

ら

[19]

と

D’yachkov

ら

[20]

は

d

^e

-disjunct

な行列の構成法を提案している．

これまで組合せ論的グループテストは多くの研究者により研究されてきた．

Sobel

ら

[58]

は多くの工業分野への応用を紹介し，

Li [37]

は組合せ論的なグループテストを導入した．

Bush

ら

[11]

と

Hwang

ら

[29]

によ

(15)

り

non-adaptive

なグループテストに関連した結果が得られている．

Vakil

ら

[66]

はポジティブが

2

個の場合のグループテストのためのプーリングデザインを提案し，

Chang

ら

[13]

はポジティブが

2

個および

3

個の場合を扱った．

Wolf [68]

は通信へのグループテストの応用について述べた．

Hwang [27]

が初めてグループテストの組合せ論的スクリーニングアル

ゴリズムのひとつを提案し，

Bonis

ら

[7]

は

2

段階グループテストのアルゴリズムについて述べた．また，

Balding

ら

[4]

は誤り訂正可能なプーリングデザインを提案し，

Macula [41]

と

Ngo

ら

[48]

は

DNA

ライブラリスクリーニングのための組合せ論的グループテストについて述べ，

Bonis

ら

[8]

や

Hwang

ら

[28]

は阻害物質と呼ばれる物質がある場合のアルゴリズムやプーリングデザインを扱った．さらに，

Thierry-Mieg [61]

はプーリングデザインのための

shifted transversal design

の有効性に言及しており，

Wu

ら

[69]

はプーリングデザインの分子生物学の分野への応用を見出している．

M´ ezard

ら

[44]

は誤りが存在するときに組合せ論的グループテストの確率論的な解析を与えた．

組合せ論的グループテストでは主に効率の良いプーリングデザインの構成法に主眼が置かれている．

1.2 グループテストの結果にもとづく確率論的なスクリーニング

図

1.3

の手順において観測値からポジティブクローンを識別するには観測値

s

を用いてポジティブクローンを識別する良いスクリーニングアルゴリズムが必要不可欠であるが，偽ポジティブ・偽ネガティブの確率が高く，ポジティブの数が想定した数を上まわることが頻繁にあるような場合には組合せ論的グループテストに対するアルゴリズムだけではポジティブクローンをスクリーニングすることが困難である．

確率論的グループテストでは，ポジティブクローンが事前確率に従って確率的に現れるという仮定のもとで，確率論的なモデルが設定される．偽ポジティブや偽ネガティブのような誤りも確率的に発生すると仮定する．

確率論的アルゴリズムでは，これらの確率モデルに従ってポジティブク

(16)

ローンをスクリーニングする．確率論的な手法は

Bruno

ら

[10]

や

Knill

ら

[35]

により開発された．

Knill

ら

[35]

は

Markov Chain Monte Carlo

(MCMC)

と呼ばれるシミュレーションの方法を利用してポジティブであ

る事後確率を求めるアルゴリズムを提案した．この方法では，プーリングデザインに

d-disjunct

な行列を用いた際にポジティブの数がプーリングデザインで定められたポジティブの最大数

d

を超えていてもそれ以上のポジティブを高い確率で識別することができる可能性がある．

通常，組合せ論的グループテストはポジティブクローンの最大数や誤りの最大数が事前に与えられているという条件のもとで行われるので，組合せ論的方法はこれらの条件を満たしている限り効率が良い．しかし，多くの組合せ論的なグループテストは偽ポジティブや偽ネガティブの数が少ない場合にしか適用できない．また，組合せ論的グループテストはプールの観測値が

2

値の場合にしか適用できない．観測値が多値や連続値をとる場合に組合せ論的グループテストを適用するためには観測値をポジティブとネガティブの

2

値に変換する必要があり実験の情報が失われてしまう．一方，確率論的な方法はポジティブである確率や誤りの確率が事前に与えられていれば，多値や連続した値をとる場合にも適用することが可能である．さらに，組合せ論的グループテストはクローンがポジティブかネガティブかという

2

値の判定しかできないが，確率論的な方法はポジティブである事後確率の大小でポジティブであるか否かを判定する．したがって，誤り確率が大きいときに高い確率ですべてのポジティブを見つけたいときには組合せ論的なスクリーニングアルゴリズムは適切ではないかもしれない．しかし，組合せ論的なプーリングデザインと確率論的な事後確率計算アルゴリズムを組み合わせることにより，効率的な識別アルゴリズムを開発することが可能であることを本論文を通して主張したい．

1.2.1 プーリング実験の確率モデル

本節ではプーリング実験の確率モデルについて述べる．

(17)

X

c を

X

c

=

 



 

0,

クローン

c

がネガティブのとき，

1,

クローン

c

がポジティブのとき，

となる確率変数とする．通常，各クローンがポジティブである事前確率

P (X

_c

= 1)

は小さく，例えば，

P(X

_c

= 1) = 0.0001, 0.001

などの値をとる．

1.1

節で述べた

Cla I

ライブラリの場合には

P (X

_c

= 1) ≈ 2.6/1298 ≈ 0.002

である．さらに，

Z

_G を

Z

G

= ∨

c∈G

X

c

で定義される確率変数とする．ただし，

∨

は

G

に属するすべてのクローンに対する

X

c の論理和である．プール

G

がネガティブなクローンしか含まないとき，

Z

G

= 0

となり，

G

が少なくとも一つのポジティブクローンを含むとき，

Z

G

= 1

となる．

S

G をプール

G

に対する観測値を表す確率変数とする．スクリーニングの実験結果の観測値は蛍光の強さなどにより測定され，以下の

4

つの値で表されることがある．

(Knill

ら

[35]

による．

)

S

_G

=

 

 

 

 



0,

プール

G

がネガティブのとき，

1,

プール

G

が弱ポジティブのとき，

2,

プール

G

が中ポジティブのとき，

3,

プール

G

が強ポジティブのとき．

このとき，誤り確率は

P (S

_G

= s | Z

_G

= z)

で表すことができる．つまり，

P (S

_G

= 1, 2, 3 | Z

_G

= 0)

が偽ポジティブで，

P(S

_G

= 0 | Z

_G

= 1)

が偽ネガティブである．

ここで，各

X

_c は互いに独立で観測値

S

_G は

Z

_G にのみ依存し，

S

_G はすべての

Z

_G が知られているという条件のもとで独立であると仮定する．

(18)

仮定

:

プール

G

1

, . . . , G

m

∈ G

に対して

P (S

_G₁

= s

_G₁

, . . . , S

_G_m

= s

_G_m

| Z

_G₁

= z

_G₁

, . . . , Z

_G_m

= z

_G_m

)

= P (S

_G₁

= s

_G₁

| Z

_G₁

= z

_G₁

) × · · · × P(S

_G_m

= s

_G_m

| Z

_G_m

= z

_G_m

), (1.2)

および

P (S

_G

= s

_G

| Z

_G

= z

_G

, X

c1

= x

c1

, . . . , X

cn

= x

cn

)

= P(S

_G

= s

_G

| Z

_G

= z

_G

)

(1.3)

が成り立つ．

つまり，

Z

G が既知のとき，式

(1.2)

はそれぞれの

S

G が

Z

G にのみ依存することを示しており，式

(1.3)

は

Z

G が既知のもとで

S

G が

X

c と独立であることを示している．

1.2.2 陽性反応の事後確率

プーリング実験の確率モデルの陽性反応の事後確率について述べる．

X = (X

c1

, . . . , X

cn

)

と

S = (S

G1

, . . . , S

Gm

)

を確率ベクトルとする．

S

が測定されたときの事後確率は

P(X = x | S = s) = P(X = x, S = s) P (S = s)

= KP (X = x, S = s)

= KP (X = x)P(S = s | X = x)

= K ∏

c∈C

P (X

c

= x

c

) ∏

G∈G

P (S

G

= s

G

| Z

G

= z

G

) (1.4)

と書くことができる．ただし，

K = P (S = s)

⁻¹は定数で，

z

G

= ∨

c∈G

x

c

である．

(19)

したがって，

X

c の条件付周辺事後確率は

P (X

_c

= x | S = s)

= ∑

x∈{0,1}ⁿ s.t.xc=x

P (X = x | S = s)

= K ∑

x∈{0,1}ⁿ s.t.xc=x

( ∏

c^′∈C

P (X

_c′

= x

_c′

) ∏

G∈G

P (S

_G

= s

_G

| Z

_G

= z

_G

) )

,

(1.5)

となる．ただし，

∑

x∈{0,1}ⁿs.t.xc=x は

x

c

= x

を満たすすべての

x = (x

c1

, . . . , x

cn

) ∈ { 0, 1 }

ⁿ についての和を意味している．さらに，

P (S

G

= s

G

| Z

G

= z

G

)

の値は事前の実験から経験的に知られていると仮定する．

D ⊂ C

をクローンの部分集合とし，集合

D

を「

D

に属するクローンがポジティブで

C \ D

に属するクローンがネガティブである事象」と同一視する．すべてのプールの観測値ベクトル

s

が与えられたもとでの

D

の事後確率は

P (D | S = s) = P (D)P (S = s | D) P (S = s)

となる．ただし，

S = (S

G1

, . . . , S

Gm

)

である．また，

D

の生起確率は

P (D) = ∏

c∈D

P(X

_c

= 1) ∏

c∈C\D

P (X

_c

= 0)

である．さらに，プール

G

の観測値

S

_Gはプール

G

がポジティブクローンを含むか否かにのみ依存すると仮定すれば

P (S = s | D) = ∏

G∈G

P (S

G

= s

G

| D) = ∏

G∈G

P (S

G

= s

G

| G ∩ D)

が成り立つ．ただし，

G ∩ D

は

G ∩ D ̸ = ∅

^のとき

Z

_G

= 1

，

G ∩ D = ∅

のとき

Z

_G

= 0

である．したがって，観測値が得られたときの事後確率は

P (D | S = s)

= KP (D)P(S = s | D)

= K ∏

c∈D

P (X

c

= 1) ∏

c∈C\D

P (X

c

= 0) ∏

G∈G

P (S

G

= s

G

| G ∩ D)

(20)

となる．

1.2.3 スクリーニング

式

(1.5)

の

P (X

c

= 1 | S = s)

は観測値

s = (s

G1

, . . . , s

Gm

)

が与えられたもとでクローン

c

がポジティブである事後確率である．観測値が

s

のときの事後確率

P (X

c

= 1 | S = s)

が大きいとき，クローン

c

はポジティブである確率が高いとみなすことができる．第

2

段階の実験では，事後確率が高い順に一定数のクローンを抽出し，ポジティブクローンを決定するためにそれらを個別にテストする．

しかし，式

(1.5)

から事後確率

P (X

_c

= x | S = s)

を計算するためには

2

ⁿ⁻¹ 回の和が必要である．したがって，クローンの数

n

が増加すると，計算量は指数関数的に増加する．そのため，この値を計算するための効率的なアルゴリズムが必要である．

1996

年に

Knill

ら

[35]

が

Markov

Chain Monte Carlo (MCMC)

法にもとづく確率論的な事後確率計算アルゴリズムを提案した．次章では，

Knill

ら

[35]

によるアルゴリズムについて概説し，この従来のアルゴリズムの欠点を補うために考案したアルゴリズムを提案する．

(21)

第 2 _章

BNPD _{アルゴリズムと} CCPD _アルゴリズム

本章では，まず従来の

Markov Chain Monte Carlo (MCMC)

法を用いた事後確率計算アルゴリズム

Markov chain pool result decoder (MCPD)

について概説し，その上で新しい事後確率計算アルゴリズムである

Bayesian network pool result decoder (BNPD)

，および

CCCP pool result decoder

(CCPD)

を提案する．さらに，

BNPD

あるいは

CCPD

で求まる周辺事

後確率の偏りに関する

1

つの指標を導く．

2.1 事後確率計算アルゴリズム MCPD

Markov Chain Monte Carlo (MCMC)

法にはさまざまな手法があるが，

本節ではギブスサンプリングという基本的な手法による周辺事後確率の計算について述べる

(Gamerman [23]

，

Geman [24])

．本手法を用いた事後確率計算アルゴリズムは

Knill

ら

[35]

により導入され，

Markov chain pool result decoder (MCPD)

と呼ばれている．

マルコフ連鎖の定常確率を利用してモンテカルロシミュレーションを行う方法を

MCMC

法という．まず，事後確率

P(X = x | S = s)

を定常確率として持つマルコフ連鎖を次のようにして構成する．

x = (x

₁

, . . . , x

_n

) ∈ { 0, 1 }

ⁿ ^とし

p(x) = P (X = x | S = s)

とおく．

(22)

{0, 1}

ⁿ 上のマルコフ連鎖を考え，

x

^′ から

x

への推移確率を

p(x | x

^′

) =

∏

n j=1

p(x

j

| x

1

, . . . , x

j−1

, x

^′_j+1

, . . . , x

^′_n

)

とする．このマルコフ連鎖は既約であり，非周期的である．実際

p(x

j

| x

1

, . . . , x

j−1

, x

^′_j+1

, . . . , x

^′_n

)

= P (X = (x

₁

, . . . , x

_j₋₁

, x

_j

, x

^′_j+1

, . . . , x

^′_n

) | S = s)

∑

x^′_j

P (X = (x

1

, . . . , x

j−1

, x

^′_j

, x

^′_j+1

, . . . , x

^′_n

) | S = s)

であるので，

P (X

j

= x

j

) > 0

，

j = 1, . . . , n

に注意すれば

P (X = x | S = s) > 0

となることより，

p(x | x

^′

) > 0

．また

p(x | x) > 0

も明らかである．さらに

∑

x^′

p(x | x

^′

)p(x

^′

)

= ∑

(x^′₁,...,x^′_n)

∏

n j=1

p(x

_j

| x

₁

, . . . , x

_j₋₁

, x

^′_j+1

, . . . , x

^′_n

)p(x

^′₁

, . . . , x

^′_n

)

= ∑

(x^′₂,...,x^′_n)

∏

n j=2

p(x

j

| x

1

, . . . , x

j−1

, x

^′_j+1

, . . . , x

^′_n

)p(x

1

, x

^′₂

, . . . , x

^′_n

)

= ∑

(x^′₃,...,x^′_n)

∏

n j=3

p(x

_j

| x

₁

, . . . , x

_j₋₁

, x

^′_j+1

, . . . , x

^′_n

)p(x

₁

, x

₂

, x

^′₃

, . . . , x

^′_n

)

· · ·

= ∑

x^′_n

p(x

_n

| x

₁

, , . . . , x

_n₋₁

)p(x

₁

, . . . , x

_n₋₁

, x

^′_n

)

= p(x

₁

, . . . , x

_n

)

= p(x).

が成り立つことより

p(x)

はこのマルコフ連鎖の定常確率になり，このマルコフ連鎖が既約かつ非周期的であることから

p(x)

はこのマルコフ連鎖の極限分布になる．

このことを利用した計算法がギブスサンプリングであり，アルゴリズムは次のようになる．

(23)

1.

初期値として

{ 0, 1 } -

ベクトル

x

⁽⁰⁾

= (x

⁽⁰⁾₁

, x

⁽⁰⁾₂

, . . . , x

⁽⁰⁾_n

)

を発生させる．

2. x

^(τ)

= (x

^(τ)₁

, x

^(τ)₂

, . . . , x

^(τ)n

)

が与えられたもとで

(a) x

^(τ+1)₁ を

p(x

₁

| x

^(τ)₂

, . . . , x

^(τ)n

)

から発生させる．

(b) x

^(τ+1)₂ を

p(x

₂

| x

^(τ+1)₁

, x

^(τ)₃

, . . . , x

^(τ)n

)

(c) x

^(τ+1)₃ を

p(x

₃

| x

^(τ+1)₁

, x

^(τ+1)₂

, x

^(τ)₄

, . . . , x

^(τ)n

)

同様に

x

^(τ+1)₄

, . . . , x

^(τ+1)n を順次発生させていく．

という手順により，

x

^(τ+1)

= (x

^(τ+1)₁

, x

^(τ+1)₂

, . . . , x

^(τ+1)_n

)

を得る．

この操作を

τ = 1, 2, 3, . . .

と繰り返して

x

⁽⁰⁾

, x

⁽¹⁾

, x

⁽²⁾

, . . .

が得られる．

この発生のさせ方から

x

⁽⁰⁾

, x

⁽¹⁾

, x

⁽²⁾

, . . .

はマルコフ連鎖に従う乱数であるとみなせる．したがって，

lim

t→∞

(1/t) ∑

_t

τ=0

I

_(x(τ)=x)

= p(x)

が成り立つ．ただし，

I

_(x=x′) は

x = x

^′ のとき

1

，

x ̸= x

^′ のとき

0

である．

X = (X

1

, . . . , X

n

)

にポジティブであるクローンの集合

D = {c

j

∈ C : X

j

= 1}

を対応させれば，マルコフ連鎖

X

⁽⁰⁾

, X

⁽¹⁾

, X

⁽²⁾

, . . .

をクローン集合の列

D

⁽⁰⁾

, D

⁽¹⁾

, D

⁽²⁾

, . . .

に対応させることができる．

ギブスサンプリングのアルゴリズムにおいて

c

_j

∈ D

^(τ) であるとする．

すなわち

X

_j^(τ)

= 1

であると仮定する．このとき

c

_j

∈ / D

^(τ+1) となる確率は

p(0 | x

^(τ+1)₁

, . . . , x

^(τ+1)_j₋₁

, x

^(τ)_j+1

, . . . , x

^(τ_n⁾

)

= P(X = (x

^(τ+1)₁

, . . . , x

^(τ+1)_j₋₁

, 0, x

^(τ)_j+1

, . . . , x

^(τ)n

) | S = s)

∑

xj

P(X = (x

^(τ+1)₁

, . . . , x

^(τ+1)_j−1

, x

j

, x

^(τ)_j+1

, . . . , x

^(τ)n

) | S = s)

= P(D \ { c

_j

} | S = s)

P (D \ {c

j

} | S = s) + P (D | S = s)

となる．同様に

c

_j

∈ / D

^(τ) であるとすると，

c

_j

∈ D

^(τ+1) となる確率は

p(1 | x

^(τ+1)₁

, . . . , x

^(τ+1)_j₋₁

, x

^(τ)_j+1

, . . . , x

^(τ)_n

)

= P(D ∪ { c

_j

} | S = s)

P(D ∪ { c

j

} | S = s) + P (D | S = s)

(24)

となる．ここで

D △ { c } =

 



 

D \ { c } , c ∈ D

のとき

, D ∪ { c } , c / ∈ D

のとき

,

の記号を導入すれば，

c ∈ D

^(τ) のとき

c / ∈ D

^(τ+1) となる確率と

c / ∈ D

^(τ)

のとき

c ∈ D

^(τ+1) となる確率はまとめて次のようになる．

P (D △ { c } | S = s)

P (D △ { c } | S = s) + P (D | S = s)

= 1

1 + P (D | S = s)/P (D △ { c } | S = s)

(2.1)

ただし，

A △ B

は集合

A

と

B

の対称差，つまり，

A △ B = (A \ B) ∪ (B \ A)

である．

MCPD

アルゴリズムでは

D

⁽⁰⁾

, D

⁽¹⁾

, D

⁽²⁾

, . . .

を以下のように発生させる．すべてのクローン

c ∈ C

の対して

D

^(τ) に

c

を追加するか削除するかを式

(2.1)

の確率でランダムに決めることにより

D

^(τ+1) を生成す

る．

D

^(τ) から

D

^(τ+1) を生成することを

t

回繰り返したとき，クローン

c

を含む

D

^(τ) の数を

t

で割ることにより

c

がポジティブである事後確率

P (X

c

= 1 | S = s)

を推定する．したがって，

Knill

ら

[35]

による

MCPD

アルゴリズムは以下のようになる．

MCPD

アルゴリズム

ステップ

1 (

初期化

): D

⁽⁰⁾

= ∅

と初期化する．

t

を反復回数，

t

0 をウォームアップの反復回数とし，

τ = 1

とする．

ステップ

2 (D

^(τ) の初期化

): D

^(τ)

= D

^(τ⁻¹⁾ とする．

ステップ

3 (D

^(τ) の更新

):

各クローン

c ∈ C

に対して，乱数を発生させることにより式

(2.1)

の確率で

D

^(τ)

= D

^(τ)

△ { c }

^とする．

ステップ

4 (

反復

): τ < t + t

0 なら

τ = τ + 1

としステップ

2

に戻る．

τ = t + t

0 ならステップ

5

へ進む．

ステップ

5 (

周辺事後確率の計算

):

各

c ∈ C

に対して，

c

がポジティブ

(25)

である周辺事後確率の推定値として

Q

_c

(1) = # { D

^(τ)

∋ c | τ = t

₀

+ 1, . . . , t

₀

+ t } t

を出力する．

MCPD

によって求めた事後確率の推定値

Q

_c

(1)

が高い順に一定数のクローンを第

2

段階の個別テストにかける．

MCPD

はシミュレーションにより周辺事後確率を計算するアルゴリズムのため，繰り返し乱数を生成する必要があり，小さい確率を精度良く推定するためには反復回数

t

を多くする必要があるので実行時間が長くなる．これらの問題点を解決するために以下でベイジアンネットワーク上の確率伝播アルゴリズムにもとづいて考案したアルゴリズムを提案する．

2.2 ベイジアンネットワークと確率伝播法

本節では，ベイジアンネットワークの定義およびその上での確率伝播法について述べる．ベイジアンネットワークは確率変数間の依存関係を表すグラフィカルモデルである

(Jensen [32]

，

Lauritzen

ら

[36]

，

Neapolitan [47]

，および

Pearl [51])

．

2.2.1 ベイジアンネットワーク

ベイジアンネットワークは以下のように定義される．離散確率変数の集合

U = { U

₁

, U

₂

, . . . , U

_n

}

とこれらの変数上で定義される同時分布

P( · )

からなる確率モデルを考える．

同時分布が

P (U = u) =

∏

n i=1

P (U

i

= u

i

| V

i

= v

i

)

と分解されているとする．ただし，

U = (U

1

, U

2

, . . . , U

n

)

で各

V

i は

U

の部分ベクトルである．

U

^{の各確率変数}

U

i をグラフの頂点とし

V

i の各

グループテストにもとづく 事後確率の近似計算アルゴリズム