日本生態学会誌59巻3号

生態学と集団遺伝学理論

 生態学とは「生物の分布と個体数を決定する相互作用 の科学的研究」と定義される（Begon et al. 2006）。そのため、 生態学において、個体の分布とその数はもっとも基礎的 で重要な情報である。さらに近年では、保全生態学の見 地からも、人間活動が生物の分布・個体数に及ぼす影響 を観測する試みがさかんに行われている。同様に、過去 に起きた個体数の増減や集団の分化・融合といった個体 群動態（demography）を知ることも、（進化）生態学にお いて、また保全生態学において今後の動態の振る舞いを 予測する上でも非常に重要である。  過去のイベントは、直接観測されていなかった場合、 人間が知る事は難しい。そのような状況下では、DNA 配 列に残された多型データ（SNP やマイクロサテライト） から過去の動態を推定する試みがなされてきた。しかし ながら、このような推定は想像するほど直感的で簡単な ものではなく、実は非常に危険な行為である。なぜなら、 DNAの多型パターンが生成される背景には、過去の個体 群動態だけでなく、その他の要因が複雑に絡み合ってい るためである（図 1A）。そのなかでも、偶然の要素が占 める割合は非常に大きい。この偶然の要素のおかげで、 ある遺伝子 A の DNA 多型パターンから推定した個体群 動態と、別の遺伝子 B の DNA 多型パターンから推定し た個体群動態が食い違うといったような「事件」が起こ るのである。このような失敗を減らすには、その推定方 法が開発された背景にある「理論」を十分に理解しなけ ればいけない。  ここでいう「理論」とは、集団遺伝学の理論である。 集団遺伝学の大きなテーマの 1 つは、DNA 多型パターン がどのような進化的要因によって形成されてきたかを理 解することである。したがって、DNA 多型データから過 去のイベントを推定したいと思う生態学者が、集団遺伝 学に興味を持つのは当然のなりゆきだろう。しかし、一 般的な生態学者は、集団遺伝学に対して微妙に距離を感 じているのではないだろうか。例えば、DNA 多型データ から「過去の動態の推定」をするという場合、生態学で は、分集団や近縁種の分岐や移住、個体数の増加や減少、 ボトルネックといった動態を、記述的に扱うことが多い。 一方、集団遺伝学では、モデルを仮定したパラメータ推 定というような定量的な解析が常とされる。これが可能 なのは、集団遺伝学が理論的に非常に厳格な分野であり、 その理論が多くの仮定のもとに構築されているからであ る。例えば、集団遺伝学理論の多くは、「個体数一定でラ ンダムに交配する（panmictic）集団」という仮定をして いる。そして、それを発展させた形で、個体数の増減や 移住を組み込むことも出来るが、そのような発展的な理 論はさらに複雑な仮定をもうけている。そして、このよ うな仮定を伴う理論のもとで、パラメータを推定するの である（図 1B）。  このような集団遺伝学に対する（それに精通していな い）生態学者の反応は、主に 2 通りに分けられるのでは ないだろうか。1 つめは、成り立たない仮定（例えば、 完全にランダム交配をしている集団はこの世に 1 つもな い）の伴う理論は、使えないから無視するという場合で ある。その状況では、モデルを用いて分岐年代などのパ ラメータ推定を行うことをせずに、例えば DNA 多型デー タからレフュージアを推論するなどといった定性的な議 1）e-mail: [email protected]

始めよう！エコゲノミクス（3）

集団内変異データが語る過去：解析手法と理論的背景（その 1）

山道 真人 *

,

_{**・印南 秀樹 *}

1） *総合研究大学院大学生命共生体進化学専攻 **日本学術振興会特別研究員 DC

論に留まることになる。もう 1 つは、仮定をしっかり理 解しないまま応用する場合である。例えば、集団遺伝学 のソフトウェアを、そのモデルの仮定にあてはまらない データに対して乱用するなどといった行為がそれにあた る。しかし、うまくつきあえば集団遺伝学は強力な解析 ツールを提供することができるため、興味深いデータを 持っている場合、上記のような例は非常に惜しい行為で ある。うっとうしく感じる仮定は所詮、仮定でしかない ので、データに合うようにうまく仮定を再設定してやれ ば、どのような場合でも対応することが可能になる。そ のことを理解していただくことを今回のメインテーマに したい。  過去の動態を推定するというテーマに対して、本連載 では今回を含めて複数回を当てたい。今回はその基礎で ある coalescent 理論の概要を紹介する。Coalescent 理論と は、遺伝子系統樹（genealogy）を現在から過去に遡って 考えることで、効率的なモデリングを可能にする理論で ある。さらに、シンプルなモデルを仮定して、coalescent 理論にもとづいたパラメータ推定法を解説する。複雑な モデルを用いて過去の動態を推定する手法の紹介は次回 に譲る事としたい。  今回のテーマである過去の動態の推定は、連載第 1 回 で挙げた「ゲノム情報がいかに生態学に貢献しうるか」 でリストアップした（1）∼（5）のポイントの内、（5） に相当する（山道・印南 2008）。これは、最後にリスト アップしたにもかかわらず、実はその他の（1）∼（4） にも大きく関わる。例えば、連載第 2 回でも取り上げた 様に、ゲノムワイド関連マッピング（2）を行う際には、 図 1．A．DNA 多型パターンを決定する要因。 まず、ゲノム全体に同様に影響を及ぼす要 因として、個体群動態・環境がある。次に、 個々の遺伝子に特異的に影響する要因とし ては、組み換え・選択・突然変異・遺伝的 浮動がある。それらのプロセスの影響を受 けて、現在見られる DNA の多型パターン が形成されてきた。多型パターンは、多型 サイト数 S や塩基多様度 p といった要約統 計量で表現される。B．要約統計量をもとに、 さまざまな仮定にもとづいた理論モデルに よって、モデルのパラメータを推定するこ とができる。

前提条件として集団構造がわかっていることが必要にな る（山道・印南 2009）。同様に、自然選択の痕跡を多型 データから検出する際（3）にも、個体群動態の効果を考 慮しなければならない。このように、単に生態学的興味 だけでなく、高度な集団遺伝学的解析を生態学に持ち込 むにあたっても、過去の動態の推定は非常に大きな役割 を果たす。

Coalescent

理論とは

 Coalescent 理論は、1980 年代初頭に Kingman・Hudson・ 田嶋によって、ほぼ同時期に独立に提唱された集団遺伝 学の理論である（Kingman 1982；Hudson 1983；Tajima 1983）。この理論がユニークなのは、1．現在から時間を 遡って後ろ向きに考えるため、現在の状態が理論の初期 状態として与えられる。2．初期状態を、集団全体の構造 ではなく、そのサブセット（サンプル）として設定する ことも可能である、という点である。例えば、10,000 個 体からなる大集団においても、そこからサンプリングし てきた 20 個体の遺伝的変異の状態を初期状態とするこ とができるのである（図 2）。これは非常に便利で、とく にコンピューターシミュレーションを行う際には、計算 時間を大幅に短縮することができる。近年では CPU の計 算速度も向上してきているため、どのような複雑な動態 でも（理論上は）計算することが可能になってきている （Rosenberg and Nordborg 2002）。Coalescent 理論の弱点は、

基本的に中立進化を仮定することであり、実際、理論に よって解析的に自然選択を考慮することは難しい。しか し、シミュレーションによって、この問題も解消されつ つあり、選択を組み込んだソフトウェアなども開発され ている（Teshima and Innan 2009）。ただし、今回のテーマ である、中立マーカーを用いた過去の個体群動態の推定 には特に関係ないので、ここでは自然選択は考慮しない。  本稿ではまず始めに、coalescent 理論の最も基本的な部 分を説明する。とりあえず、古典的な「個体数一定でラ ンダムに交配する」集団という仮定を置いて話を進める 図 2．Coalescent 理論の基礎となる遺伝子系統樹。個体数が一定（N ＝ 10）で、交配がランダムに起こる半数体生物の 中立な遺伝子の系統樹の例。時間について前向きに考えると、1 世代目からスタートして、第 9 世代目で集団中の 全個体が 1 世代目の 1 個体の子孫となる。後ろ向きに考えると、MRCA は現在から遡って第 7 世代前で現れる。現 在の世代でサンプルされた 5 個体とそれらの祖先個体を黒で、それ以外の個体を白で表している（Rosenberg and Nordborg（2002）より改変）。

が、次回以降、徐々にその制限を緩めてゆく。また、理 論的な部分に関して少し数式が出てくるが、これはシミ ュレーションが理論の結果と合うということを示すため だけに用いるのであり、数式が理解できなくても今後い っさい問題ないので、安心して読み飛ばしてほしい。こ こで紹介する理論は、個体数一定でランダム交配という 仮定のときのみ成り立つものであるが、非常に有意義な 直感的理解を促すものなので、あえて割愛しない。  一方、個体数一定でランダム交配という仮定を緩める と、解析的理論は難解になる。この問題を解決するため、 ここでは、計算機によるシミュレーションを活用する方 法を中心に紹介する。もちろん、シミュレーションには 難しい数式は出てこない。そして、データ解析にも柔軟 に応用することができる。弱点は、理論的に厳密でなく なることだが、これも「うまく」やればそれほど問題は ない。このように、いかにシミュレーションベースの理 論がデータ解析に有用であるかを示したい。

Coalescent

の基礎理論

 本来、coalescent 理論は 1 世代前の個体数が与えられれ ば成立する柔軟さを持っているが、ここでは簡単のため に、個体数が一定（N）で、離散世代であり、全個体が同 じ確率で子孫を残す 2 倍体の有限集団を考えてみよう（こ れは Wright-Fisher model と呼ばれる）。このような集団で は、各個体の残す子の数の期待値は 1 になるが、これは 確率的なモデルであるため、中には子を残せないもの・2 個体の子を残すもの・3 個体以上の子を残すものがあら われる。子の数は正の整数なので、期待値 1 のポアソン 分布をすることになる。すなわち、約 37％の確率で 0 個 体、約 37％の確率で 1 個体、約 18％の確率で 2 個体、約 6％の確率で 3 個体、約 2％の確率で 4 個体を、非常にま れに 5 個体以上を、子供として残すような確率分布にな る。従って、世代を重ねると徐々に最初の世代から存続 している系統の数は減少し、時間が充分経つと、最終的 には集団中の全ての個体が、最初の世代のある 1 個体の 子孫という状況になる（図 2）。以上で述べてきたような 考え方が時間的に前向き（forward）な考え方である。  この考え方にもとづくモデルでシミュレーションを行 うには 2 つの難点がある。1 つ目は、現在からは知り得な い第 1 世代の状態を仮定して出発しなくてはならないとい う点である。2 つ目は、現在まで子孫を残さず、サンプル されることがない系統についても考えなくてはならないた め、無駄が多いという点である。図 2 で言えば、現在まで 存続し、サンプリングされるのは太線の系譜のみだが、前 向きのモデルならば、不必要である細線の系譜も含めた全 てについて計算しなくてはならない（Nordborg 2001）。  そこで考え出されたのが時間的に後ろ向き（backward） な coalescent 理論である。先程のモデルを逆に現在の状 態から考えてみよう。現在の世代 t の状態から出発して、 1つ前の世代 t-1 を考えてみる。すると、集団中からサン プルしてきた 2 本の塩基配列が親世代の同じ個体の同じ 染色体から得られる（coalesce する）確率は、1 本目の 親となる染色体を選んだうえで、2 本目の親となる染色 体が 1 本目のそれと同じであるという条件付き確率なの で、1/(2N) になる。逆に、それらが別の染色体から得ら れる確率は 1-1/(2N) になる。このようにして時間を遡っ ていけば、どれだけ遠縁のサンプル同士であっても、い ずれはもっとも最近の共通祖先（Most Recent Common Ancestor; MRCA）において coalesce することになる。こ の coalescence までの待ち時間の期待値は 2N 世代となる （Box 1）。このような後ろ向きのモデルを用いれば、祖先 状態を仮定する必要がなく、かわりに、現在の集団が初 期条件になる。また、現在まで存続せず途中で消える系 統を考える必要がなくなるという利点がある（図 2）。  この理論は、前述のようにサンプルベースなので、サ ンプルからなるDNA多型データ解析に非常に有用である。 ここでは、まず、DNA 多型データ解析の基本中の基本で ある変異量に関する coalescent 理論を紹介する。多型デー タの変異量は、しばしば塩基多様度（nucleotide diversity; p） と多型サイト数（segregating site; S）という要約統計量で 表現される（第 1 回参照）。Coalescent 理論は、このモデ ルの仮定のもとでの期待値をたちどころに与える。  まず、同一集団からサンプリングしてきた独立な 2 本 の配列を考えよう。この配列の長さを Lbp（base pairs） とする。前述のように（Box 1 も参照）、2 本の配列が時 間を遡って MRCA に至るまでの時間の期待値は 2N 世代 である。突然変異率を m（世代あたり・サイトあたり） とすると、MRCA から現在までに 2 本の系統に蓄積され る突然変異の数の期待値は、 （1） となる。ここで、集団サイズと突然変異率の積に 4 をか けたもの（4Nm）をθ として定義する。集団遺伝学において、 この θ は最も重要なパラメータの一つである。Wright-Fisher modelにおいては、θ はその集団中の変異を決める 唯一のパラメータとなる。ここで、θ は集団サイズ N と

突然変異率 m の積であることに注意されたい。つまり、 個体数 10,000・突然変異率 10-8_{の時と、個体数 100,000・} 突然変異率 10-9_{の時、期待される変異量はまったく等し} くなるということである。  θ が充分小さければ、新しい突然変異は常に別のサイ トに起き、同じサイトに 3 つ以上の対立遺伝子は存在し ないと仮定できる。これを無限サイトモデル（infinite-sites model）と言う。このモデルは θ が 10-2_{以下のオーダーな} らば、非常によい近似を与える。この条件下では、2 本 の配列を考えたとき、多型サイト数と塩基多様度の期待 値は等しくて、両方とも θ L となる。しかし、3 本以上の 配列を比較する時、その期待値は異なる。3 本以上であ っても、塩基多様度は総当たりペアを比較した時の多型 サイト数の平均値なので、 （2） となる。一方、多型サイト数は、n 本の塩基配列間を比 較した時、その内の 1 本でも他の配列と異なる変異（多型） を持っている塩基対（サイト）の数なので、期待値は、 （3） となる。このような式になる理由を、3 本の塩基配列を 比較する状況で説明してみよう（Box 1 も参照）。3 本の 塩基配列の内、2 本が coalesce するまでの待ち時間の期待 値は 2N/3C2＝ 2N/3 である。その期間中は、存在する 3 本 の系統に突然変異が蓄積するので、変異数は 2N/3 × 3 × m × L ＝ 2Nm L となる。従って、多型サイトの期待数は、 4Nm L ＋ 2Nm L ＝ θ L(1+1/2) となる。これを一般化すると、 n 本の配列が (n-1) 本に coalesce するまでの時間の期待 値が 4N/[n(n-1)] であり、その間に蓄積する突然変異数が θ L/(n-1) ということになるため、無限サイトモデルにお ける多型サイト数の期待値は（3）式で表せる（Watterson 1975）。このように、coalescent 理論は、集団中の DNA レ ベルの変異量を非常に簡単に表現できる。

Coalescent

シミュレーション

 シミュレーションとは、コンピュータ上で実験的に検 証する手段である。シミュレーションを用いることで、 解析理論と同様の結果を理論なしで得ることができる。 つまりシミュレーションは、ある設定のもと、ある程度 の回数の試行を繰り返し行うことによって、理論と同様 の役目を果たすのである。例えば、θ ＝ 0.01・L ＝ 1,000 と決めて、coalescent シミュレーションを行い、塩基多様 度 p の理論値を得てみよう。1 回目は p ＝ 9.1 で、2 回目 は 11.9、3 回目は 5.6、4 回目は 4.3、5 回目は 13.0、とい うように繰り返し計算して、それらの値を平均すると、 徐々に理論値である 10 に近付く。だいたい 10,000 回か ら 100,000 回繰り返すと、平均値はほぼ理論値に等しく なる。このように、シミュレーションによって理論値を 得るには、時間がかかる。これがシミュレーションの最 大の弱点である。一方、シミュレーションは解析的理論 ではうまく取り扱えないようなさまざまな状況に柔軟に 対応することが可能で、非常に有用である。シミュレー ションの利点である柔軟さについては、次回に複雑な個 体群動態を考慮したモデルを用いて解説することにする。 今回は、シミュレーションがいかに理論の代役をこなす か、またそれがどこまで正確か、という点について、以 下の θ の推定を例にとって説明する（coalescent シミュレ ーションの手順については、Box 2 にまとめたので参照さ れたい）。

パラメータ推定：モーメント推定法

 上で紹介した coalescent 理論は、直接データ解析に使え る。最も直感的な応用は、「集団突然変異率」と呼ばれる パラメータ θ（＝ 4Nm）の推定である。これは、塩基多 様度（p）および多型サイト数（S）から推定可能である。 式（2）及び（3）からそれぞれ、 （4） （5） となる。ここで、データから推定された θ はハットを付け て区別している。下付き文字はそれぞれ、p にもとづく、 あるいは S にもとづくθ の推定値ということを示している。  例として、図 3A のようにサンプル数が n ＝ 20 で、配 列長が L ＝ 1,000 bp の時、多型サイト数が S ＝ 8、塩基 多様度がおよそ p ＝ 2.9 である仮想データを考える。p ＝ 2.9を使うと、（4）式から θ の推定値は ˆθp＝ 2.9 × 10-3と なり、S ＝ 8 を使うと、（5）式から θ の推定値はおよそ θˆS＝ 2.3 × 10-3となる。  この例の場合、推定法は異なっても同じデータをつか っているため、似たような推定値を与える（もちろん、 方法が違うので一致はしない）。このような、単に割り算

して求めた推定値をモーメント推定値という。  このモーメント推定値を用いて、さらに集団サイズの 推定を行うこともできる。突然変異率 m が既知であり、 10-8_{（世代あたり・サイトあたり）であったとしよう。θ ＝} 4Nm という関係から、p にもとづく N の推定値はおよそ 74,000個体となる。同様に、S にもとづく N の推定値は およそ 56,000 個体となる。  しかし、このような理論ベースの推定はあまりデータ 解析には向いていない。なぜなら、理論の仮定が合って いない場合、推定値自体にあまり意味がないためである。 本連載ではむしろ、そのような仮定をデータに合うよう に再設定することがテーマである。これが、シミュレー ションを使うことが望ましい理由のひとつである。

パラメータ推定：最尤法

解析的理論とシミュレーション

 あるモデルにもとづいてパラメータを推定する際、尤 度（likelihood）という考え方が重要になる。尤度とは、「尤 も（もっとも）らしさの度合い」という名前の通り、パ ラメータセット m を与えた条件でデータ D が得られる確 率であり、条件付き確率の形式で表せば P(D|m) となる。 この確率がもっとも大きくなる、すなわち最大尤度を持 つパラメータセットを求める推定方法が最尤法である。 最尤法は、複雑なモデルのもとで複数のパラメータを推 定する状況で力を発揮するため、広く用いられている。  まず、解析的に多型サイトの尤度を求めてみよう。配 列が coalesce するまでの時間は指数分布、coalesce するま 図 3．A．仮想的な DNA 塩基配列の多型データ（SNPs）。サンプルサイズは n ＝ 20 で、多型サイトが 8 つ見つかった例。20 個 のサンプル中、もっとも頻度が多い対立遺伝子をそのサイトのコンセンサス（Cons.）として最上列に示した。B．多型サ イト数を要約統計量とした θ の最尤推定の結果。横軸が θ × 10-3_{、縦軸が最大値を 0 とした相対的な対数尤度。理論的な尤} 度関数によって算出した尤度を実線で、シミュレーションによって算出した尤度を四角形の点で表している。C．塩基多様 度を要約統計量としたθ の最尤推定の結果。B と同様、横軸が θ × 10-3_{、縦軸が最大値を 0 とした相対的な対数尤度。シミ} ュレーションによって近似尤度を計算する際、許容するずれ δ が 1％のものを四角形、10％を青い菱形、100％を赤い三角 形の点で表している。

での時間に起こる突然変異の数はポアソン分布に従って 確率的にばらつく。そこで、2 本の配列を比較した時（n ＝ 2）の多型サイト数の尤度を式で表してみる。集団サイズ を N、遺伝子系統樹の全長を T とすると、全長 T は期待 値 4N の指数分布に従う。また、突然変異の数（＝多型サ イト数）は期待値 m TL のポアソン分布をする。そこで尤 度関数は、coalescent までの時間と突然変異数の積を時間 で積分して、 （6） となる（Watterson 1975）。3 本以上の配列を比較する場合 （n ＝ 3, 4, 5, …）も、同様にして畳み込み積分（convolution） を繰り返し、尤度を求めることができる。  ここでもう一度、図 3A の仮想データをつかう。（5） 式から、多型サイト数が 8 のとき、θ のモーメント推定 値はおよそ 2.3 × 10-3_{となる。では、最尤法ではどのよ} うな推定値になるだろうか。（6）式のような畳み込みに よって得られた式を用いて n ＝ 20 の時の尤度 P(S ＝ 8|θ ) を計算し、その尤度曲線を図 3B に示す。ここで、x 軸に θ（× 103_{）の値を、y 軸に最大値を 0 とした相対的な対} 数尤度を取った。理論的な尤度関数に基づき計算された 対数尤度を実線で示している。モーメント推定値の 2.3 × 10-3_{と比較すると、実際に同じ値（θ ＝ 2.3 × 10}-3_）で最 大の尤度が得られている。ちなみに、最大対数尤度と対 数尤度の差の 2 倍は χ2_{乗分布に従うので、（最大対数尤} 度−対数尤度）がおよそ 1.9 以内の範囲（図 3B 横線）が おおよその 95％信頼区間になる。  次に、シミュレーションによって多型サイトの尤度を 求める。この際には、ある θ の値を与えて多数回（今回 は 200,000 回）coalescent シミュレーションを行って得ら れた SNP データから多型サイト数 S を計算し、S がちょ うど 8 になったシミュレーションの割合を尤度とする。 そして θ の値を 10-4_{ずつ変化させ、最も尤度が高くなる} θ を探索する。その結果が図 3B の四角い点によって示さ れている。最尤推定値は理論値と同様 2.3 × 10-3_{となった。} 図 3B において、シミュレーションによって算出した尤度 が理論的に求めた尤度関数ときれいに一致することがわ かるだろう。このように、尤度を求めるという点で、シ ミュレーションが完全に解析理論の代役を果たすことが 解っていただけたと思う。  では、塩基多様度 p を用いて θ を推定するときはどう なるだろうか。解析的に尤度関数を求めることは（不可 能ではないが）かなり難しい。そこで、シミュレーショ ンによって塩基多様度の尤度を求めるしかないが、シミ ュレーションによって求める際にも、多型サイト数のよ うに簡単ではない。なぜなら、先程と同様に、尤度をデ ータと同じ値を取る試行の割合とすると、大変なことに なるためである。図 3A の例では、正確な p の観察値は 2.947368421…であり、これとぴったり同じ p の値がシミ ュレーションで出てくることは滅多にない。S からの推 定がうまく行ったのは、多型サイト数が正の整数しか取 らないからである。  そこで、考えだされたのが近似最尤法である（Tavaré et al. 1997；Pritchard et al. 1999；Marjoram et al. 2003）。これ は、データとまったく同じ塩基多様度の値を取る試行だ けを「データとマッチした」としてカウントするのでは なく、値がそこそこ近い試行もカウントすることによっ て、シミュレーションの効率をあげる方法である。この ずれの絶対値を δ と表現する。δ ＝ 0 の場合が正確な最 尤法にあたる。そして、若干のずれを認める場合（δ ＞ 0）、 それから求められた尤度は近似的なものとなる。  問題は、どれくらいのズレなら許されるのか、という 点である。これは実は非常に難しい問題なのだが、ここ では、直感的にそれを理解してもらうために以下の実験 をおこなう。まず、かなり厳密な場合（δ が p の値の 1％） を試してみる（図 3C の四角い点）。これは真の尤度にか なり近いと考えられる。ここから δ を増やすと、どんど ん真の尤度から離れるのは想像に容易い。そこで次に、δ＝ 10％を試してみる（青い菱形の点）と、δ ＝ 1％の近似尤 度曲線とほぼ同じものが得られた。これは何を意味する かというと、δ ≤ 10％の範囲で、シミュレーションは非常 に良い近似尤度を与えるということである。それなら、δ＝ 1％より δ ＝ 10％を採用した方が、計算速度（シミュレ ーションの時間）を 10 倍短縮できる。では、もっと効率 を上げるために δ ＝ 100％にしてみよう。直感的にこれ はやり過ぎと思われるであろうが、もしうまくいけば儲 け物である（計算速度が 100 倍速くなるため）。しかし、 残念ながら 100％の場合（赤い三角形の点）の近似尤度 曲線は、δ ＝ 1％のそれとは大きくかけ離れるものになり、 使い物にならない。  このように、うまく δ を設定できれば、シミュレーシ ョンは解析理論の代役を見事に果たすので、非常に有用 である。たとえ解析的に尤度が求められない場合でも、 シミュレーションを使えば最尤法と実質的に同様の結果 が得られるということは特筆すべきポイントである。こ れを示すために、不本意ながらいくらかの難解な数式を 用いた。理解しようと挑戦した方には申し訳ないが、こ

れらの式はシミュレーションを使う限りあまり意味がな い、というか忘れてもらって結構で、問題なく数式をフォ ローできた方は、本格的な集団遺伝学への転向を勧める。  なお今回の例で、δ ＝ 10％程度が、正確性を失わない 範囲でもっとも計算効率を上げる δ であると結論したが、 これは常に当てはまる法則ではないことを肝に銘じてほ しい。完全にケースバイケースである。これが、シミュ レーションによる近似最尤法の最も批判される点である。 そのため、実際に近似尤度を用いて推定する際には、こ こで紹介したような簡単な事前実験で適当な δ を決める プロセスを行うことを強く推奨する。この点さえ気をつ ければ、シミュレーションは万能で、（理論家としては残 念なことではあるが）エレガントな解析的理論の価値も かすんでしまう程である。

より正しく推定するために必要な努力

 推定値というものは常に分散を持つ。これは当たり前 ではあるが、意外と軽く扱われている事実である。分散 がある以上、推定値はあやふやなものである（すなわち、 誤差がある）。直感的な理解としては、例えば推定値が 10でその 95％信頼区間が 5 から 15 であった場合、真の 値は 20 回に 1 回は信頼区間の外にあると考えてもらえば 良い。これを考えれば、いわゆる「推定値」が、いかに 怪しげなものであるかは想像に難くない。しかしながら、 自分が得た推定値をあたかも正しいと仮定してその後の 議論を続ける行為は、生物学では頻繁にみられる。これは、 あまり褒められたものではない。重要な事は、その誤差 の存在を認識し、それがどのように結論に影響するかを 理解することである。さらには、より正確に小さい誤差 で推定値を得るように心がけることである。  ここでは、多型サイト数及び塩基多様度から θ を推定 する場合を例にとって、この問題を考えたい。図 3A の データから得られた S を用いた θ の推定では、最尤推定 値は 2.3 × 10-3_{であるが、その 95％信頼区間は広くおよ} そ 1 ∼ 6 × 10-3_{になる（図 3B）。一般的に p を用いた θ} の推定値はもっと大きな信頼区間を持つ（図 3C）。この ように、図 3A の例程度のデータ量では満足に信頼でき る θ の推定値を与えないのである。では、どうしたらい いだろうか。もちろんデータ量を増やせばいいのだが、 やみくもに増やすのは賢くない。ここでは、効果的な増 やし方を理論的に考えてみよう。  まず、データを増やせと言われれば、サンプル数（個 体数）を増やすか、シークエンスする領域の長さを増や すことを考えるだろう。しかし、これはあまりコストパ フォーマンスが高くない。図 4A では、サンプル数と推 定値の分散の理論的関係をプロットした。サンプル数を 増やせば分散は減っていくが、（特に塩基多様度の）分散 はなかなかゼロに近くならない。同様に、シークエンス する領域の長さを増やした場合も、効果は低いことがわ かるだろう（図 4B）。もっと効率の良い方法はないもの だろうか。  最も効果的な方法は、ひとつの遺伝子座のデータを使 うのではなく、複数（出来るだけ多く）の遺伝子座のデ ータを使うものである。突然変異率はゲノム中でわりと 一定であるため、θ もどの遺伝子座でもほぼ一定である と言える。統計学的に、独立な標本の数を増やすと、分 散は標本数の逆数で減っていく（図 4C）。従って、もし 10遺伝子座を使えば分散は 1/10 になり、100 遺伝子座を 使えば 1/100 になる。このように、遺伝子座数を増やす と、分散はどんどん減っていく。図 5 では、遺伝子座数 が 1・10・100 のときの θ の最尤推定値の頻度分布をプ ロットしている。真の θ は 0.01 と設定して仮想データを 作り、それぞれに対して最尤法で θ を推定した。遺伝子 座数が 1 の場合の分布は広く、正しく θ ＝ 0.01 を推定す ることは稀であるが、遺伝子座数が 100 の場合の分布は、 ほとんどの試行において非常に良い推定値を得られると 期待できる。100 遺伝子座というデータ量は、理想的な 推定値を与える基準と考えてよい。もちろん、遺伝子座 数を増やすことは、新たなプライマーを開発するなどと いった手間がかかり、常に最も効率的であるとは言えな い。しかし、次回に紹介するような、より踏み込んだ個 体群動態の推定を目指すならば、ある程度の遺伝子座数 は絶対に必要になる。逆に言えば、1 遺伝子座を用いた 個体群動態の推定は、非常に信頼度の低いものといえる （詳しくは次回に解説する予定である）。

理想的なデータについて

 非モデル生物の過去の個体群動態を推定する際には、 しばしばミトコンドリアや葉緑体といったオルガネラゲ ノムの塩基配列が用いられる。しかし、このような推定 には問題がある。ミトコンドリアは組み換えのない 1 つ の遺伝子座に相当するので、いかに長い配列を読んでも、 ミトコンドリアだけで推定を行うと、図 5 のように信頼 区間が広く、正確な推定値が得られないという結果にな る。ミトコンドリアのみに基づいた個体群動態の推定が いかに意味のない危険な行為であるかがおわかりいただ

けると思う（この点に関しては、もう一度次回で触れる）。 一方、同様に非モデル生物でよく用いられる核ゲノム上 のマイクロサテライトでは、複数のマーカーが用いられ ることが多い。複数のマイクロサテライトは、複数の遺 伝子座に相当するので、ミトコンドリア単体よりもずっ と良い推定ができる。  マイクロサテライトの問題は、そこで起こっている突 然変異のメカニズムが、点突然変異と違って複雑である ことである。マイクロサテライトでは、反復配列が DNA 複製時にその反復数を変えるかたちで変化する。その際、 反復数がどのように増減するのか（1 個ずつか、または 一気に 5 個か）といったメカニズムを、実際の現象に合 うようにモデル化できるなら、今回紹介したようなシミ ュレーションを用いた近似最尤法は簡単に応用できる （Box 2）。ちなみに集団遺伝学では、一回の突然変異で増 減する反復配列の数を幾何分布で近似するモデルがよく 使われる。そして、S や p でなく、マイクロサテライト のヘテロ接合頻度などを要約統計量に用い、シミュレー ションをとおして θ やさらには個体群動態の推定を行う ことが出来る。  結論としては、ミトコンドリアのみにもとづいて推定 を行うことは問題があり、マイクロサテライトがベター だが、理想的には核の塩基配列多型（の複数の領域）を 用いるのが望ましい、ということになる。

まとめと次回予告

 今回の連載のポイントは、coalescent シミュレーション が十分理論の代用になり、解析理論なしで最尤法と事実 上同じアプローチがとれるということである。さらにシ ミュレーションには、柔軟性という理論にはない利点が ある。これを生かして、θ だけでなく複雑な個体群動態 を推定することが可能である。また、今回あえて無視した、 シークエンスする領域内での組換えも簡単に組み込むこ とが出来る。このような、より複雑な coalescent シミュレ ーションとその応用は、次回以降で扱う予定である。 図 4．サンプル数・配列長・遺伝子座数を増やした時の、多型サイト数・塩基多様度から推定した θ の分散の変化。黒点 が塩基多様度、灰色の点が多型サイト数にもとづく値。A．サンプル数 n を n ＝ 2, 3, …, 100 と増やした場合。固定 条件：L ＝ 1,000 bp で θ ＝ 0.01、1 遺伝子座。B．配列長 L を 100 bp から 10,000 bp まで増やした場合。固定条件：n ＝ 20で θ ＝ 0.01、1 遺伝子座。C．遺伝子座数 l を 1 から 10 まで増やした場合。固定条件：n ＝ 20、L ＝ 1,000 bp で θ ＝ 0.01。図 4A ∼ C で、n ＝ 20、L ＝ 1,000 bp、l ＝ 1 という共通条件で縦線を引き、その際の分散の値に横線を引いた。 図 5．複数遺伝子座を用いた θ の最尤推定値の頻度分布。横軸が θ の最尤推定値で、縦軸が最尤推定値の頻度。遺伝子数が 1 の場合の最尤推定値の分布を点線で、遺伝子数が 10 の場合 を破線で、100 の場合を実線で表した。複数の遺伝子座を l 個用いて最尤推定する場合、全体の尤度は、各遺伝子座のデ ータから計算された尤度の積となる。つまり、遺伝子座 i の データを Diとすると尤度 P は となる。 遺伝子数が増えると、ばらつきが小さくなり、最尤推定値が 真の θ の値に近付く。サンプル数は 20、配列長は 1,000 bp、 真の θ の値は 0.01 で、シミュレーションの繰り返し回数は 10,000回。

Box 1. 集 団 中 か ら サ ン プ リ ン グ し て き た 配 列 の coalescent 理論  1 世代前に遡った時、2 本の配列が coalesce する確 率は 1/(2N) である。したがって、T-1 世代前に遡った 時点でまだ coalesce しておらず、T 世代前で coalesce する確率は、   （7） となる。式の左は幾何分布と呼ばれるもので、ある 確率で成功する試行を繰り返して、初めて成功する までの待ち時間である。例えば、6 面体のさいころを 振って 1 が出るまでの待ち時間（試行回数＝ x）は、 で、試行回数 x の平均（期待値）は 6 回になる。 また、パチスロで 1/200 のボーナスを引くまでの試行 回数は で、期待値は 200 回となる。これ らの例と同様にして、確率 1/(2N) で coalesce するま での待ち時間の期待値は 2N 世代である。式の右は指 数関数で、N が 100 以上あれば非常によく幾何分布 を近似するため、集団遺伝学ではよく用いられる。  次に、集団中からサンプルしてきた 3 本の配列の coalescentプロセスを考えよう。3 本が 1 つ前の世代 で 2 本に coalesce する確率は、3 つの中から 2 つを選 んでくる組み合わせを考えて、   （8） となり、1 つ前の世代で 3 本のままである確率は 1-3/(2N) となる。そこで、T 世代前で 3 本の配列が 2 本へ coalesce する確率は、   （9） となる。この分布の平均値は 2N/3 となる。   同 様 に し て、n 本 の 配 列 が T 世 代 で n-1 本 へ coalesceする確率は、 （10） となり、平均は 2N/nC2である（Tajima 1983）。このよ うにして、配列が何本の時でも確率にしたがって遺 伝子系統樹を描くことができる。 Box 2. Coalescent シミュレーションによる多型デー タの生成  4 本の配列を、個体数 N の集団からサンプリング してきた状況を考えよう（図 6）。まず、4 本の内か らランダムに 2 本を選び、その 2 つが coalesce する までの時間 t4（世代数）を決定する。Coalescent 時間 は、平均 2N/4C2の指数分布をするので、指数乱数を 発生させて時間を得る。同様にして、3 本が 2 本に、 2本が 1 本に coalesce する過程で、どの個体がどれ くらいの時間（t3・t2）で coalesce するかをランダム に決定する。それぞれの値は、平均 2N/3C2と 2N/2C2 の指数分布をする（図 6-1）。そして、この遺伝子系 統樹のうえに、突然変異をランダムにばらまく（図 6-2∼ 4）。突然変異の個数は、遺伝子系統樹の全長 T （T ＝ 4t4+3t3+2t2）と突然変異率 m（世代あたり・サイ トあたり）、配列長 L に依存する。シミュレーション する配列の長さが L bp であれば、突然変異の個数は 平均 Tm L のポアソン分布からの乱数になる。突然変 異は、枝の長さに応じてランダムに遺伝子系統樹の 上に配置する（図 6-4）。突然変異が起こった DNA 塩 基の場所を L bp の領域上にランダムに決めれば、そ れが直接塩基配列の多型パターンに反映される（図 6-5）。同様にして、マイクロサテライトの多型もシ ミュレーションできる（図 6-5’）。この場合、例え ば祖先タイプの反復数を 20 と仮定し、それぞれの突 然変異が幾何分布に従って反復数を増減させるプロ セスをシミュレートすることになる。このようなシ ミュレーションは非常に簡単なプログラムで書くこ とが出来るが、web 上で公開されている ms（Hudson 2002）のようなソフトを使うことも出来る（http:// home.uchicago.edu/~rhudson1/source/mksamples.html を 参照）。

引 用 文 献

Begon M, Townsend CR, Harper JL (2006) Ecology: from individuals to ecosystem, 4th ed. Blackwell Publishing, Oxford

Hudson RR (1983) Properties of a neutral allele model with intragenic recombination. Theor Pop Biol 23:183-201 Hudson RR (2002) Generating samples under a Wright-Fisher

neutral model of genetic variation. Bioinformatics 18:337-338

Appl Probab 19:27-43

Marjoram P, Molitor J, Plagnol V, Tavaré S (2003) Markov chain Monte Carlo without likelihoods. Proc Natl Acad Sci USA 100:15324-15328

Nordborg M (2001) Coalescent theory. In: Balding DJ, Bishop M, Cannings C (eds) Handbook of statistical genetics. John Wiley & Sons, Ltd., West Sussex, pp 179-212

Pritchard JK, Seielstad MT, Perez-Lezaun A, Feldman MW (1999) Population growth of human Y chromosomes: a study of Y chromosome microsatellites. Mol Biol Evol 16:1791-1798

Rosenberg NA, Nordborg M (2002) Genealogical trees, coalescent theory and the analysis of genetic polymorphisms. Nat Rev Genet 3:380-390

Tavaré S, Balding DJ, Griffiths RC, Donnelly P (1997) Inferring

coalescence times from DNA sequence data. Genetics 145:505-518

Tajima F (1983) Evolutionary relationship of DNA sequences in finite populations. Genetics 105:437-460

Teshima KM, Innan H (2009) mbs: modifying Hudson's ms software to generate samples of DNA sequences with a biallelic site under selection. BMC Bioinformatics 10:166 Watterson GA (1975) On the number of segregating sites in

genetical models without recombination. Theor Popul Biol 7:256-276 山道真人・印南秀樹 (2008) 始めよう！エコゲノミクス (1) 局所適応と形質の分化. 日本生態学会誌 58:241-247 山道真人・印南秀樹 (2009) 始めよう！エコゲノミクス (2) ゲノムワイド関連マッピング. 日本生態学会誌 59:105-113 図 6．Coalescent シミュレーションから多型パターンを生成する方法。説明は Box 2 を参照。