JP2005092719A

JP2005092719A - ハプロタイプ推定方法、推定装置、プログラム

Info

Publication number: JP2005092719A
Application number: JP2003327943A
Authority: JP
Inventors: Toshio Furuta; 利夫古田; Masao Yanagisawa; 政生柳澤; Naoyuki Kamatani; 直之鎌谷
Original assignee: Waseda University; NEC Corp
Current assignee: Waseda University; NEC Corp
Priority date: 2003-09-19
Filing date: 2003-09-19
Publication date: 2005-04-07
Also published as: US20050089906A1

Abstract

【課題】実行時間を指数関数オーダから多項式オーダへと削減し、大きなデータへの対応を可能とすること。
【解決手段】ＥＭアルゴリズムとグラフ構造を組み合せることにより、取り得る全てのハプロタイプ情報を陰に保持し、ハプロタイプ推定を最大スコアを持つ完全グラフを探索する問題に変換する。
【選択図】図９

Description

本発明は、ゲノム多型マーカを用いた研究に広く利用することができる、ハプロタイプ推定方法に関する。ここで、多型マーカを用いた研究とは、相関解析やテーラーメイド医療等を指す。

周知のように、ヒトゲノムの約３０億塩基対の配列の決定が報告された。人の遺伝子の塩基の配列（ＤＮＡ塩基配列）において、個体間に１％以上存在する変異の事を遺伝子多型（genetic polymorphism）と呼ぶ。遺伝子多型は単に多型とも呼ばれる。すなわち、多型とは、同じ種の個人個人の遺伝的違いのうち、遺伝子による違いのことをいう。多型にはさまざまなレベルのものがある。

尚、周知のように、ＤＮＡとはデオキシリボ核酸（deoxyribonucleic acid）の略称で、塩基成分はほとんどがアデニン（Ａ）、グアニン（Ｇ）、シトシン（Ｃ）とチミン（Ｔ）の４種類である。

特に、１個の塩基について存在する多型を一塩基多型（Single Nucleotide Polymorphism,ＳＮＰ）と呼び、遺伝子多型の中でも特に頻度が高いことから特に注目されている。ＳＮＰは数百塩基対から１，０００塩基対に1か所程度の割合で存在していると推測されているので、ゲノムの中には、３００万〜１，０００万のＳＮＰがあると考えられている。人における個人差は、このＳＮＰの違いにより大部分が決定されることが推測されている。

ＳＮＰは、（１）ゲノム内に数百万コピー以上存在すること、（２）判定が非常に容易であり、また、結果を（０，１）信号化することができるため情報処理が容易であること、（３）高速・大量のＳＮＰタイピング技術が実現化されつつあり、タイピングのオートメーション化も視野に入ってきていることから、多型マーカとして利用しやすく、臨床応用のための機器開発が比較的容易であると考えられている。

従って、以下では、多型マーカとして一塩基多型を用いる場合を例に挙げて説明するが、マイクロサテライト等の他の多型マーカにも適用可能である。

ＳＮＰの遺伝子機能や遺伝子発現に与える影響を解明することにより、病気にかかりやすい体質をつきとめたり、個人の体質に合わせたよりよい治療法、薬剤の選択や医薬品の開発が可能になると考えられている。

全ＳＮＰより遺伝子上に近接して存在する複数のＳＮＰを、遺伝子の機能単位で一括抽出した組をハプロタイプ（haplotype）と呼ぶ。

Ｘ，Ｙを塩基とすると、ＳＮＰは１塩基部位に対してＸ／Ｘ，Ｘ／Ｙ，Ｙ／Ｙの３パターンを取る。ただし、Ｘ，ＹはＡ（アデニン），Ｔ（チミン），Ｇ（グアニン），Ｃ（シトシン）のいずれかの塩基を現す。

遺伝子を解析する方法では、主に機能的アプローチおよび遺伝統計的アプローチの２つのアプローチが採用される。前者は、遺伝子に含まれるＳＮＰを分子生物学的に分析することによりＳＮＰの機能を解析するアプローチである。後者は、ＳＮＰを臨床情報と組み合わせ、遺伝統計的に関連を抽出するアプローチである。遺伝統計学の分野において、遺伝子データベースおよび臨床データベースを対象として、両データベースに潜在する関連を統計的に抽出する方式が研究されている。遺伝統計学は、両データベースに潜在する関連の傾向を統計的に抽出することを可能としている。

遺伝統計学の基礎となる法則はメンデルの法則と連鎖の法則である。メンデルの法則を確率事象に関する法則として捉えることが重要である。メンデルは、人の認識できる生物の固体の属性としての「表現型」の奥に、より本質で、人には認識できない「遺伝子型」が存在すると考えた。そして、その「遺伝子型」も個体にとってのみ安定なのである。そして、世代を超えて安定なアレル（すなわち遺伝子）を発見した。メンデルの法則とは以上の「表現型」、「遺伝子型」、「アレル」の関係に関する法則と捉えることができる。そして、その３つの概念を定義できる場としての座位が定義される。メンデルの法則は、後述する「分離の法則」、「優劣の法則」、「独立の法則」に要約される。

個体レベルで観察される個体の性格を表す言葉に表現型（phenotype）と形質（trait）がある。形質とはカテゴリーのことを示しており、表現型はその形質において観察される個々のタイプを示す。個々の個体の形質の関係を解析するにはメンデルの考えたように座位、アレル、遺伝子型の概念を定義する必要がある。

１つの交配集団に対して、多数の座位（Locus）が存在すると考える。座位は種においてほぼ安定である。個体において安定な単位は遺伝子型（Genotype）である。１つの個体については、１つの座位に１つの遺伝子型が存在する。染色体による表現では座位は組換えが起きる最小単位の領域にはさまれた部分である。分子生物学的には座位は一塩基である。座位は遺伝子座とも呼ばれる。

遺伝子型は個体のレベルでは安定であるが、世代を超えると、もはや安定ではない。個体を超えて安定して伝えられるのは遺伝子型ではなくアレル（Allele）である。アレルは対立遺伝子とも呼ばれる。

メンデルの分離の法則は遺伝子型を構成する要素のうち、世代を超えて安定な単位に関する法則である。分離の法則は「遺伝子型は２つのアレルの組合わせであるが、世代交代の際に、遺伝子型を構成するアレルが１個だけ、等確率で次世代に伝えられる」ことを主張する（つまり確率論）。そのため、世代を超えてはアレルが最小の安定した単位である。もちろん、突然変異によりアレルは破壊されるが、メンデルの法則は突然変異の存在を想定していない。染色体のレベルではアレルは片方の染色体上の部位であり、分子生物学では片方の染色体上のＳＮＰ、またはＳＴＲＰ（short tandem repeat polymorphism）、ＶＮＴＲ（variable number of tandem repeat）などの型である。例えばＳＮＰではアレルはＴ（チミン）あるいはＣ（シトシン）などの塩基の一方である。

メンデルの優劣の法則は「遺伝子型と表現型の対応に関する法則であり、個体の表現型が遺伝子型の関数であることを主張し」、その関数に優性と劣性が（場合によっては共優性）が存在することを示す。染色体のレベルでは、遺伝子型は２つの相同染色体上の相同部分に存在する２つの相同部位の組合せである。分子生物学のレベルでは遺伝子型は多型のあるＳＮＰ、またはＳＴＲＰ、ＶＮＴＲなどの組合せである。

連鎖は複数の座位に関するメンデルの独立の法則の例外である。メンデルの独立の法則は複数の座位の間の関係に関する法則である。「異なる座位にあるアレルの、子への分配は互いに独立である」ことを主張する。染色体のレベルでは、これは２つの座位が異なった染色体上にあるときにのみ正しい。同じ染色体上の２つのアレルは物理的に結合しているので、つながったまま同じ配偶子に入り、子に伝わる。すなわち、独立の法則が成立しない。しかし、同じ染色体上に存在するアレルが完全に結合して世代交代を続けるとは限らない。減数分裂（miosis）のときに、染色体が交差して結合していたアレルの組合せ変化することがあるからである（組換え）。組換えの可能性は２つの座位間の距離が離れていればいるほど大きい。

座位の異なる遺伝子においては必ずしもメンデルの独立の法則が当てはまらないことを示したのはモルガンである。すなわち連鎖（Linkage）の発見である。しかし、異なった染色体上の座位では依然として独立の法則は正しく、しかも連鎖の概念が成立するためには独立の法則の概念が存在することが前提であることに注意が必要である。独立の法則はメンデルの誤りであるという解釈は間違っている。

ここで、独立という意味は２つの異なった事象が同時に起きる確率が、それぞれの事象が起きる確率の積で表されるということである。遺伝統計学では、独立という意味は、同じ染色体上にあるという意味ではない。すなわち、独立とは確率、統計的な意味である。

連鎖が定義されると、座位１、座位２の各アレルの組合わせが定義できる（haplotype，ハプロタイプ）。連鎖する複数の座位のアレルのうち、１つの配偶子に存在するものの組合せがハプロタイプである。ハプロタイプは減数分裂での組換えがない限り保存される。組換えが起きると新たなハプロタイプが形成され、それが配偶子により次世代に伝えられ、次の組換えが起きるまで変化しない。ハプロタイプの存在は、その前に親から伝えられた配偶子のハプロタイプと、各座位間で組換えが起きたか否かの事象の関数である。組換えが起こる事象は組換え割合により確率的に起こるので、ハプロタイプの遺伝は組換え割合を確率分布関数にもつ確率事象と考えられる。家系の形質と遺伝子型に関する観測データに最も適合する座位の順番と距離を計算するのが連鎖解析である。

個体は２つのハプロタイプをもつが、この２つのハプロタイプの組合わせをディプロタイプ形（diplotype configuration）ということもある。

１つの家系にｎ個の配偶子の伝達があるとし、ｍ個の、連鎖した座位があるとする。この家系において起きた減数分裂における組換えの事象を表すｍ個の継承ベクトルを考える。１つの継承ベクトルは１つの座位に対応し、それはｎ個の因子をもつ列ベクトルである。それぞれの因子は、順番をつけた各減数分裂に対応する。継承ベクトルの因子は０か１であり、０はその減数分裂の起きた個体が、その父親から継承したアレルを配偶子により子に伝達したことを示す。１は母親から継承したアレルを配偶子により子に伝達したことを示す。すなわち、継承ベクトルで表した場合、本質はアレルではなく、親由来のアレルをどちらの子に伝達したかという事象にあると考える。

一般に、家系のすべての構成員の順位つき遺伝子型（順列遺伝子型）は、すべての創始者の順位つき遺伝子型と継承ベクトルによりすべて決定できる。順位つき遺伝子型とは、どちらのアレルが父親由来かという情報つき２つのアレルの組合わせ（すなわち順列）である。遺伝子型１/２と２/１は異なった順位つき遺伝子型ということになる。通常、遺伝子型とはアレルの順列ではなく、組合わせであることに注意する（組合わせ遺伝子型）。

連鎖不平衡（Linkage disequilibrium）の概念は遺伝統計学で最も大切な概念の１つである。連鎖と連鎖不平衡とはどう異なるか、など、正確に理解する必要がある。連鎖の概念は、メンデルの独立の法則の例外であることは前述した。そして、独立とは純粋に確率的な概念であることを述べた。真の連鎖不平衡は連鎖を前提とした概念である。

複数の連鎖する座位があると、それらの座位に遺伝子型が存在する。そして家系情報が十分に得られれば個人のディプロタイプ形（２つのハプロタイプの組合わせ）が確定できることは前述した。ハプロタイプとは、１つの染色体上に存在する（すなわち連鎖している）それぞれの座位に１つのアレルの組合わせである。家系情報の存在しない場合、個体の遺伝情報についての完全情報が個体のディプロタイプ形である。

連鎖不平衡とは、集団においてのハプロタイプの分布に関する法則である。ディプロタイプ形は２つのハプロタイプの組合わせなので、集団には個体数の２倍のハプロタイプが存在する。

連鎖不平衡とは、連鎖する２つ以上の遺伝子座で観察される現象である。簡単に言えば、連鎖不平衡とは、異なる遺伝子座間における対立遺伝子の分配が独立ではない現象、すなわち、ハプロタイプ頻度が各遺伝子座における対立遺伝子の頻度を掛け合わせて得られる頻度からずれている現象である。換言すれば、ハプロタイプ頻度がアレル頻度の積で表されるとき（すなわち独立のとき）、連鎖平衡という。連鎖平衡にない状態を連鎖不平衡という。

もし、集団が十分に大きく、すべての突然変異の起きた後に十分な時間が経ち、しかもアレル頻度に変化がなければ連鎖平衡が達成されているであろうと予想される。しかし、実際には集団は十分に大きくなく、十分な時間は経っておらず、しかもアレル頻度には変化がある。また、近い過去に２つの連鎖平衡が達成された集団が混合したなら、混合集団には強い連鎖不平衡があると予想される。一般に、近い遺伝的距離の中では世界のほとんどの集団で強い連鎖不平衡が存在する。距離が離れた多型であるほど連鎖不平衡の程度は低い。新しい突然変異に関係したハプロタイプほど強い連鎖不平衡があると予想され、集団の有効な大きさが大きいほど弱い連鎖不平衡があると予想される。

尚、遺伝的距離（Genetic(map) distance）は「２つの座位間で起きる交差の回数の期待値」と定義される。すなわち、遺伝的距離は組換えではなく、交差をもとに定義される距離の概念である。単位は、Ｍ（モルガン）。１Ｍは１回の減数分裂において１回の交差が期待できる距離として定義される。ここで、１Ｍの距離の間に必ず１回交差が起きるのではなく、確率的な概念であることに注意が必要である。また、交差（Crossover）とは減数分裂の過程で起こる相同染色体同士の部分的な乗り換え現象のことを指す。交差は顕微鏡で染色体の現象としてみることができる。これは動作を示す言葉である。１回の減数分裂において起こる交差の回数は１回とは限らない。

連鎖不平衡の存在は、遺伝統計学のための多くの手法を提供する。中でもcommon disease common variant common originという予想が正しければ、疾患の変異と周囲の変異に強い連鎖不平衡が存在し、連鎖不平衡を利用して疾患関連部位を探し出すことができると予想される。

連鎖不平衡の強い領域は、歴史的に組み換えが多い領域に挟まれた形で存在する。この領域はblockと呼ばれ、物理的距離で平均１０〜１５ｋｂ程度であるといわれている。しばしば１００ｋｂほどのblockもあることが分かっている。ただし、ｋｂとは千塩基対（kilobase pair）であり、ｐｋｂとは、ｐ千塩基対の距離を示す。尚、物理的距離（Physical(map) distance）は遺伝子の塩基配列を読んだときの２座位間の物理的な距離を表す。

相関解析等においては、通常等間隔に多型マーカを配置する。ここで、相関解析（association study）は通常、ケース−コントロール研究（case-control study,患者−対照研究）といわれているものを指す。これは、最も容易に行うことができる解析法である。相関解析も連鎖不平衡解析と同様に、家系データのない個体の集団を取り扱う。患者集団と一般集団のサンプルを独立に収集し、その間のアレル頻度の違い、あるいは特定のアレルを持つ個体の頻度の差を分析する。よって、１００ｋｂのblockには１０〜１５ｋｂのblockと比較してより多くのマーカを配置することになる。例えば、１００ｋｂ以上にわたる巨大blockには３０個以上のＳＮＰをタイピングしてしまうこともある。

しかし、block領域ではおよそ５種類以下のハプロタイプで累積頻度が９５％以上となるため、１block当りに必要なＳＮＰ数は２^２＝４、２^３＝８より最低３個である。せいぜい５種類程度のハプロタイプを識別するために、３０個のＳＮＰは全く必要なく、無駄である。

ＳＮＰs（single nucleotide polymorphisms）は、ヒトゲノム中の２５０−３５０ｂｐごとに分布する（例えば、非特許文献１参照）。マイナーアレルの頻度が０．１以上のＳＮＰはおよそ６００ｋｂごとに存在する（例えば、非特許文献２参照）。ＳＮＰsは、突然変異率が低く、高スループットなタイピングが可能であるため、相関解析においてマイクロサテライトマーカより有効である。また、ＳＮＰのハプロタイプ情報は、連鎖不平衡の検定や疾患遺伝子の染色体へのマッピングなどに対して非常に有効である（例えば、非特許文献３、非特許文献４、非特許文献５参照）。ＳＮＰのハプロタイプ情報は、単体ＳＮＰよりも頑健かつ強力に連鎖不平衡（ＬＤ）の位置情報を推定することができる（例えば、非特許文献６、非特許文献７参照）。

分子レベルでハプロタイプを推定する手法としては、単分子希釈（single-molecular dilution）（例えば、非特許文献８参照）、広範囲ＰＣＲ（long-range PCR）（非特許文献９参照）、等温円回転増幅（isothermal rolling-circle amplification）（例えば、非特許文献１０参照）、倍数染色体から半数染色体への変換（例えば、非特許文献１１参照）などがある。しかしながら、これらの手法は自動化が困難、高コスト、低スループットなど多くの問題がある。

そのため、分子レベルで推定を行う代わりに統計学的にハプロタイプを推定する手法が注目されている。家系データが与えられたときは、相またはハプロタイプを決定する手法としては、Linkage Package（例えば、非特許文献１２参照）やGENEHUNTER（例えば、非特許文献１３参照）のようなソフトウェアや、いくつかの規則に基づいて家系内で最小の組み換えハプロタイプを推定する手法（非特許文献１４参照）などがある。

一方、家系データが存在しない場合は、Clarkのアルゴリズム（例えば、非特許文献１５参照）やハーディ・ワインバーグ平衡（Hardy-Weinberg equilibrium: HWE）を仮定して遺伝子を数えることによる（例えば、非特許文献１６参照）ＥＭアルゴリズム（例えば、非特許文献１７参照）を用いた最尤推定手法（例えば、非特許文献１８、非特許文献１９、非特許文献２０、非特許文献２１、非特許文献２２、非特許文献２３参照）、偽Gibbs Sampler (pseudo Gibbs Sampler)による手法（例えば、非特許文献２４参照）、ＰＬ法によるもの（例えば、非特許文献２５参照）など様々なものがある。

ＥＭアルゴリズムによる推定は、最も精度がよいとされており現在主流となっている。

詳述すると、ＥＭアルゴリズムによる推定は、ハプロタイプ頻度の初期値に適当な値を設定し、E-stepでＨＷＥを仮定してハプロタイプの期待値を求め，M-stepでハプロタイプ頻度の更新を行うことを繰り返しハプロタイプ頻度の最尤推定値を得るという手法で、現在最も多く用いられている。

非特許文献２２において、北村らは集団のハプロタイプ頻度を用いてベイズ法によって個体ごとのハプロタイプを決定する手法を提案している。非特許文献２１において、伊藤らは複数の個体のＤＮＡを含んだＤＮＡ poolから得られた遺伝子型より集団のハプロタイプ頻度およびＤＮＡ pool中のハプロタイプ（haplotype copies）の組合せを推定する手法を提案している。

ＥＭアルゴリズムはE-stepをＨＷＥに依存しているが、Fallinらはハーディ・ワインバーグ不平衡（Hardy-Weinberg disequilibrium: HWD）のもとでのＥＭアルゴリズムによるハプロタイプ推定の頻度を評価し、ＨＷＤによって推定精度が下がるが、ＨＷＤに伴うホモ接合度の増加によって推定精度が上がるため、両者の間にバランスが生じていると述べている（非特許文献２６参照）。

ＥＭアルゴリズムを用いた手法は、多くの種類のデータに対して精度の高い推定が期待できるため現在主流になっている。しなしながら、ＥＭアルゴリズムを用いた手法では、取り得る全てのハプロタイプ頻度情報を保持しなければならず、計算量の問題から多座位の遺伝子型データを推定することが困難である。

その問題を解決するため、最近では、非特許文献２４において、StephensらがＰＧＳ（pseudo Gibbs Sampler）を用いたPhase法を提案している。これはClarkのアルゴリズムの改良とcoalescence modelに基づいたシミュレーションを組み合せたGibbs Samplerで、ＥＭアルゴリズムを用いたものと同等以上の精度を持ち、より多くの遺伝子座に対応している。

Zhangらは、ＳＮＰのうち、ハプロタイプを代表するtag ＳＮＰを用いることによってマーカ数を削減できることを示唆している。

ゲノムワイドはハプロタイプを考えたとき、ヒトゲノムはそれぞれ独立したハプロタイプブロックに分割することができることが知られている（例えば、非特許文献６、非特許文献２８、非特許文献２９、非特許文献３０参照）。

Dalyらは、非特許文献６において、Crohn病に関係すると考えられるヒト染色体5p31の５００ｋｂの範囲、マイナーアレル頻度が５％以上のＳＮＰを１０３個解析し、その範囲が１１のブロックに分割できることを発見した。また、Johnsonらは、非特許文献３０において、９遺伝子の１３５ｋｂの範囲におけるＳＮＰを１２２解析することにより、３８４人のヨーロッパ人のハプロタイプを特徴づけるためには３４ＳＮＰsで十分であることを示した。Niuらは、非特許文献２５において、ＰＬ（Partition-Ligation）法と呼ばれる、遺伝子座をハプロタイプブロックごとに区切りそれぞれについてハプロタイプを解析手法によって多くの遺伝子座を扱うことに成功している。

尚、連鎖不平衡は、一般的な集団において平均３ｋｂ以上になることはなく、全ゲノムを用いた相関解析には５００，０００以上のＳＮＰがいると考えられている（例えば、非特許文献３１参照）。

また、ヒト染色体5p31領域のＩＢＤ（Identity by descent）に関するＳＮＰデータが知られている（例えば、非特許文献３２参照）。尚、ＩＢＤとは、２人の間で対立遺伝子が先祖の同じものを共有する（これを同祖であるという）状態をいう。

一方、ＥＭアルゴリズムによる遺伝統計学的手法とマルコフモデルにおける辺重み最大化の問題をあわせることで実行時間を大幅に削減する、ハプロタイプ推定アルゴリズムが提案されている（例えば、非特許文献３３参照）。この提案されている手法は、入力データから同遺伝子座以外の任意の２遺伝子座のアレル間の遷移確率をＥＭアルゴリズムによって求める第１の処理と、患者の遺伝子型データから、各遺伝子座のアレルをノードと見たてて、マルコフモデルを作成し、その患者が尤も取り得るハプロタイプを推定する第２の処理とからなる。換言すれば、個体のディプロタイプを、アレルをノードとするグラフ構造として表していることになる。

また、ハプロタイプブロックによる分割とマルコフモデルにおける辺重み最大化により多遺伝子座の遺伝子型データに対応可能なハプロタイプ推定手法も提案されている（例えば、非特許文献３４参照）。ハプロタイプブロックを分割する点を除いて、上記非特許文献３３と同様な手法である。
L. Beaudet et al., "Homogenous assays for single-nucleotide polymorphism typing using AlphaScreen," Genome Res., Vol.11, pp.600-608, 2001 D. G. Wang et al., "Large scale identification, mapping and genotyping of single-nucleotide polymorphisms in the human genome," Science, vol.280, pp.1077-1082, 1998 S. E. Hodge, M. Boehnke and M. A. Spance, "Loss of information due to ambigunous haplotyping of SNPs," Nat. Genet., Vol.21, pp.360-361, 1999 M. J. Rieder, S. L. Taylor, A. G. Clark and D. A. Nickerson, "Sequence variation in the human angiotensin converting enzyme," Nat. Genet., vol.22, pp.59-62, 1999 N. Risch and K. Merikangas, "The future of genetics studies of complex human diseases," Science, vol.273, pp.1516-1517, 1996 M. J. Daly et al., "High-resolution haplotype structure in the human genome," Nat. Genet., vol.29, pp.229-232, 2001 J. K. Pritchard, "Are rare variants responsible for susceptibility to complex diseases?," Ann. Hum. Genet., vol.69, pp.124-137, 2001 G. Ruano, K. K. Kidd and J. C. Stephens, "Haplotype of multiple polymorphisms resolved by enzymatic amplification of single DNA molecules," in Proc. Natl. Acad. Sci. USA, vol.87, pp.6296-6300, 1990 S. Michalatos-Beloin et al., "Molecular haplotyping of genetic markers 10kb apart by allele-specific long-range PCR," Nucleic Acids Res., vol.24, pp.4841-4843, 1996 P. M. Lizardi et al., "Mutation detection and single-molecule counting using isothermal rolling-circle amplification," Nat. Genet., vol.19, pp.225-232, 1998 J. A. Dougles et al., "Experimentally-derived haplotypes substantially increase the efficiency of linkage disequilibrium studies," Nat. Genet., vol.28, pp.361-364, 2001 G. M. Lathrop et al., "Multilocus linkage analysis in humans: detection of linkage and estimation of recombination," Am. J. Hum. Genet., vol.37, pp.482-498, 1985 L. Kruglyak, M. J. Daly, M. P. Reeve-Daly and E. S. Lander, "Parametric and nonparametric linkage analysis: a unified multipoint approach," Am. J. Hum. Genet., vol.58, pp.1347-1363, 1996 M. Qian and L. Beckmann, "Minimum-Recombinant haplotyping in pedigrees," Am. J. Hum. Genet., vol.70, pp.1434-1445, 2002 A. G. Clark, "Inference of haplotypes form PCR-amplified samples of diploid populations," Mol Biol Evol, vol.7, pp.111-122, 1990 J. Ott, "Counting methods (EM algorithm) in human pedigree analysis: linkage and segregation analysis," Ann. Hum. Genet., vol.40, pp.443-454, 1977 A. P. Dempster, N. M. Laird and D. B. Rubin, "Maximum likelihood from incomplete data via the EM algorithm," J. Roy. Statist. Soc. Ser. B, vol.39, pp.1-38, 1977 M. N. Chiano and D. G. Clayton, "Fine genetic mapping using haplotype analysis and the missing problem," Ann. Hum. Genet., vol.62, pp.55-60, 1998 L. Excoffer and M. Slatkin, "Maximum-likelihood estimation of molecular haplotype frequencies in a diploid polulation," Mol. Biol. Evol. vol. 12, pp.921-927, 1995 M. Hawley and K. Kidd, "HAPLO: a program using the EM algorithm to estimate the frequencies of multi-cite haplotypes," J. Hered., vol.86, pp.409-411, 1995 T. Ito et al, "Estimation of haplotype frequencies, linkage-disequilibrium measures, and combination of haplotype copies in each pool by use of pooled DNA data," Am, J. Hum, Genet., vol.72, pp.384-398, 2003 Y. Kitamura et al., "Determination of probability distribution of diplotype configuration (diplotype distribution) for each subject from genotypic data using the EM algorithm," Ann, Hum. Genet., vol.66, pp.183-193, 2002 J. C. Long, R. C. Williams and M. Urbanek, "An E-M algorithm and testing strategy for multiple locus haplotytes," Am. J. Hum. Genet., vol.56, pp.799-810, 1995 M. Stephens, N. J. Smith and P. Donnely, "A new statistical method for haplotype reconstruction from population data," Am. J. Hum. Genet., vol.68, pp.978-989, 2001 T. Niu, Z. S. Qin, X. Xu and J. S. Liu, "Baysian haplotype inference for multiple linked single nucleotide polymorphisms," Am. J. Hum. Genet., vol.70, pp.157-169, 2002 D. Fallin and N. J. Schork, "Accuracy of haplotype frequency estimation for biallelic loci, via the expectation-maximization algorithm for unphased diploid genotype data," Am. J. Hum. Genet., vol.67, pp.947-959, 2000 K. Zhang, P. Calabrase, M. Nordborg and F. Sun, "Haplotype block structure and its applications to association studies: power and study designs," Am. J. Hum. Genet., vol.71, pp.1386-1394, 2002 E. Dawson et al., "A first-generation linkage disequilibrium map of human chromosome 22," Nature, vol.418, pp.544-548, 2002 S. B. Gabriel et al, "The structure of haplotype blocks in the human genome," Science, vol.296, pp.2225-2229, 2002 G. C. L. Johnson et al., "Haplotype tagging for the identification of common disease genes," Nat. Genet., vol.29, pp.233-237, 2001 L. Kruglyak, "Prospects for whole-genome linkage disequilibrium mapping of complex disease genes," Nat. Genet., vol.22, pp.139-144, 1999 ＵＲＬ：http://www-genome.wi.mit.edu/human/IDB5/ 下里二郎、甲藤二郎、"多遺伝子座対応可能なハプロタイプ推定アルゴリズム"、ＦＩＴ（情報科学技術フォーラム）２００２、Ａ−１５、第２９−３０頁下里二郎、甲藤二郎、"ブロック分割を用いた多遺伝子座のハプロタイプ推定に関する一考察"、２００３年電子情報通信学会総合大会、Ｄ−１２−４１、第２０２頁

前述したように、ＥＭアルゴリズムによる推定は、最も精度がよいとされており現在主流となっている。しかしながら、ＥＭアルゴリズムには、次に述べるような欠点がある。すなわち、集団内で存在する可能性のあるハプロタイプを全て保持し、繰り返し計算時においてその全てを考慮して計算しなければならないため、座位数ｎに対して計算量がＯ（２^ｎ）で増加する点が挙げられる。ｎ＜３０であれば、計算機上でも実行可能であるが、ｎ＞３０を超えると、既存のＥＭアルゴリズムによるハプロタイプ推定は計算機上で実行困難となる。すなわち、ＥＭアルゴリズムでは集団のとりうるハプロタイプを全て保持するため、遺伝子座が増加すると、実行時間、メモリ使用量ともに指数関数的に増大するという欠点がある。

前述したように、全ゲノムを用いた相関解析には５００，０００以上のＳＮＰがいると考えられているため、ＥＭアルゴリズムによるゲノムワイドなハプロタイプ推定は不可能である。

換言すれば、ゲノム解析に対する需要の向上、シーケンスなどの技術的な向上により、解析規模が拡大し、一度に処理するデータ量は拡大していくものと考えられる。座位数ｎが３０を超えるようなデータを解析するためにはＭＣＭＣ法やＥＭアルゴリズムのように集団内のハプロタイプそれぞれに対して何かしらの値を持つ手法では、扱うデータ量が膨大となり、計算機上で実行することは困難になる。

また、前述したように、ＰＬ法により、遺伝子座をハプロタイプブロックごとに区切りそれぞれについてハプロタイプを解析することが提案されているが、現在ハプロタイプブロックの明確な分割法は確立されておらず、適切なブロックに分割するためアルゴリズムの開発が必要である。

更に、非特許文献３３，３４で提案された手法では、ノード（頂点）としてアレルを採用しているので、個体ごとのディプロタイプ形を一元的に扱うことが困難になるという問題がある。

本発明の課題は、上記のような背景に鑑みて、座位数ｎが大きい遺伝子型データに対応したハプロタイプ推定方法を提供することにある。

本発明の他の課題は、個体ごとのディプロタイプ形を一元的に扱うことができる、ハプロタイプ推定方法を提供することにある。

本発明者らは、ＥＭアルゴリズムに代表されるハプロタイプ推定アルゴリズムの多くが抱える計算量の問題を改善する新たなアルゴリズム（ハプロタイプ推定方法）の提案を行った。提案手法では、ＥＭアルゴリズムとグラフ構造を組み合せることにより、取り得る全てのハプロタイプ情報を陰に保持し、ハプロタイプ推定を最大スコアを持つ完全グラフを探索する問題に変換した。

具体的に述べると、本発明の第１の態様によれば、コンピュータシステムを用いて、座位数がｎの遺伝子型データから各個体のディプロタイプ形を推定する方法であって、個体のディプロタイプ形を、順列遺伝子型を頂点に対応させたｎ部完全グラフ構造として表すステップと、ｎ部完全グラフ構造の各辺の重みを、遺伝子型データから最尤推定によって推定したディプロタイプ形の確率とするステップと、所定のスコアを用いて、ｎ部完全グラフ構造の頂点集合を選択するステップと、各個体のディプロタイプ形を、最大スコアを持つｎ次完全グラフ探索問題に置き換えて推定するステップとを含むディプロタイプ形の推定方法が得られる。

上記ディプロタイプ形の推定方法において、最尤推定としてＥＭアルゴリズムを用いることができる。また、遺伝子型データとして欠測を含むデータを扱う個体については、欠測を除いてディプロタイプ形を推定し、その結果を用いて補完を行うステップを更に含んでも良い。

本発明の第２の態様によれば、コンピュータシステムを用いて、座位数がｎの遺伝子型データからハプロタイプを推定する方法であって、個体のディプロタイプ形を、順列遺伝子型を頂点に対応させたｎ部完全グラフ構造として表すステップと、ｎ部完全グラフ構造の各辺の重みを、遺伝子型データから最尤推定によって推定したディプロタイプ形の確率とするステップと、所定のスコアを用いて、ｎ部完全グラフ構造の頂点集合を選択するステップと、各個体のディプロタイプ形を、最大スコアを持つｎ次完全グラフ探索問題に置き換えて推定するステップと、その推定により得られた各個体のディプロタイプ形より集団のハプロタイプ頻度を求めるステップとを含むハプロタイプ推定方法が得られる。

上記ハプロタイプ推定方法において、最尤推定としてＥＭアルゴリズムを用いることができる。また、遺伝子型データとして欠測を含むデータを扱う個体については、欠測を除いてディプロタイプ形を推定し、その結果を用いて補完を行うステップを更に含でも良い。

本発明の第３の態様によれば、座位数がｎの遺伝子型データから各個体のディプロタイプ形を推定する装置であって、個体のディプロタイプ形を、順列遺伝子型を頂点に対応させたｎ部完全グラフ構造として表す手段と、ｎ部完全グラフ構造の各辺の重みを、遺伝子型データから最尤推定によって推定したディプロタイプ形の確率とする手段と、所定のスコアを用いて、ｎ部完全グラフ構造の頂点集合を選択する手段と、各個体のディプロタイプ形を、最大スコアを持つｎ次完全グラフ探索問題に置き換えて推定する手段とを含むディプロタイプ形推定装置が得られる。

上記ディプロタイプ形推定装置において、最尤推定としてＥＭアルゴリズムを用いることができる。また、遺伝子型データとして欠測を含むデータを扱う個体については、欠測を除いて前記ディプロタイプ形を推定し、その結果を用いて補完を行う手段を更に含むことが好ましい。

本発明の第４の態様によれば、座位数がｎの遺伝子型データからハプロタイプを推定する装置であって、個体のディプロタイプ形を、順列遺伝子型を頂点に対応させたｎ部完全グラフ構造として表す手段と、ｎ部完全グラフ構造の各辺の重みを、遺伝子型データから最尤推定によって推定したディプロタイプ形の確率とする手段と、所定のスコアを用いて、ｎ部完全グラフ構造の頂点集合を選択する手段と、各個体のディプロタイプ形を、最大スコアを持つｎ次完全グラフ探索問題に置き換えて推定する手段と、この推定により得られた各個体のディプロタイプ形より集団のハプロタイプ頻度を求める手段とを含むハプロタイプ推定装置が得られる。

上記ハプロタイプ推定装置において、最尤推定としてＥＭアルゴリズムを用いることができる。また、遺伝子型データとして欠測を含むデータを扱う個体については、欠測を除いて前記ディプロタイプ形を推定し、その結果を用いて補完を行う手段を更に含むことが好ましい。

本発明の第５の態様によれば、コンピュータに、座位数がｎの遺伝子型データから各個体のディプロタイプ形を推定させるためのプログラムであって、コンピュータに、個体のディプロタイプ形を、順列遺伝子型を頂点に対応させたｎ部完全グラフ構造として表す機能と、ｎ部完全グラフ構造の各辺の重みを、遺伝子型データから最尤推定によって推定したディプロタイプ形の確率とする機能と、所定のスコアを用いて、ｎ部完全グラフ構造の頂点集合を選択する機能と、各個体のディプロタイプ形を、最大スコアを持つｎ次完全グラフ探索問題に置き換えて推定する機能とを実現させるためのプログラムが得られる。

上記プログラムにおいて、最尤推定としてＥＭアルゴリズムを用いることができる。また、コンピュータに、遺伝子型データとして欠測を含むデータを扱う個体については、欠測を除いてディプロタイプ形を推定し、その結果を用いて補完を行う機能を更に実現させることが望ましい。

本発明の第６の態様によれば、コンピュータに、座位数がｎの遺伝子型データからハプロタイプを推定させるためのプログラムであって、コンピュータに、個体のディプロタイプ形を、順列遺伝子型を頂点に対応させたｎ部完全グラフ構造として表す機能と、ｎ部完全グラフ構造の各辺の重みを、遺伝子型データから最尤推定によって推定したディプロタイプ形の確率とする機能と、所定のスコアを用いて、ｎ部完全グラフ構造の頂点集合を選択する機能と、各個体のディプロタイプ形を、最大スコアを持つｎ次完全グラフ探索問題に置き換えて推定する機能と、この推定により得られた各個体のディプロタイプ形より集団のハプロタイプ頻度を求める機能とを実現させるためのプログラムが得られる。

本発明によれば、ＥＭアルゴリズムの持つ推定精度の高さ、小さいデータでの実行時間の短さといった利点を生かしながら、実行時間を指数関数オーダから多項式オーダへと削減し、大きなデータへの対応が可能となる。また、頂点としてアレルではなく順列遺伝子座を採用しているので、個体ごとのディプロタイプ形を一元的に扱うことができる。

以下、本発明にかかるハプロタイプ推定方法について詳細に説明する。

ＥＭアルゴリズムは少ない座位数であれば、高精度にハプロタイプ推定を高速に実行することが可能である。ＥＭアルゴリズムの利点を生かすために、ＥＭアルゴリズムを用いた既存のハプロタイプ推定手法をグラフ構造に組み込むことにより、時間計算量・空間計算量を削減し、同程度の推定精度を確保している。本発明により提案されたハプロタイプ推定方法を、ここでは“ldlight”と呼ぶことにする。

座位数ｎ＞３０の遺伝子型データに対応するためには、既存のアルゴリズムのように集団内で考えられている全てのハプロタイプを考慮する方法では、必要となるメモリ量が爆発的に増加するため実現不可能である。

この問題を解決するためには、個々のハプロタイプに関する情報を持たずに、全てのハプロタイプ情報を考慮することのできる新たなデータ構造が必要となる。

遺伝子型データの各座位間には、連鎖不平衡が形成されている。連鎖不平衡は、連続した座位間のみに生じるものではなく、離れた領域間に強い連鎖不平衡が生じる場合もある。このため、ハプロタイプを推定するためには、全ての座位間に生じる連鎖不平衡を考慮する必要がある。

集団では１つの座位においてアレルは複数存在するため、データ構造は異なる座位のアレル間の連鎖不平衡を考慮する必要がある。よって、図１に示されるように、各座位のアレルを頂点（ノード）とし、異なる座位の頂点間を辺で結んだデータ構造が集団のハプロタイプを考慮する上で有効であると考えられる。

座位数をｎとすると、集団が持つ新たなデータ構造はｎ部完全グラフとして扱うことができる。ｎ部完全グラフは、下記の数１３で表される。

ここで、ｍ_ｉとは座位ｉにおけるアレル数である。各座位の頂点には各座位で考えられるアレルを与え、頂点間を結ぶ各辺には辺が結ぶ頂点間の結びつきの強さを表す数値を辺の重みとして与える。結び付きの強さとしては、ＥＭアルゴリズムによって２座位間を結ぶ辺に対応する集団のハプロタイプ頻度を用いる。これは、２座位にまたがるアレル間の結びつきの強さはもともと集団内に存在するアレルの割合、および連鎖不平衡の強さを考慮する必要があるためであり、ＥＭアルゴリズムにより推定されたハプロタイプ頻度はアレル頻度、連鎖不平衡の強さの両方を考慮しているからである。

本発明よるハプロタイプ推定方法では、集団のハプロタイプ分布を直接求めるのではなく、集団内の各個体のハプロタイプをそれぞれ推定し、それらを集計することにより求める方法を採用している。よって個体のハプロタイプを推定するために集団の持つ特徴から構成された、数１３で表されるグラフ構造から個人の持つディプロタイプ形を推定するためには、個人の持つ遺伝子型に合わせた新たなグラフ構造を数１３から生成する必要がある。

各個人では１つの座位においてアレルは２つ存在するため、個体ごとのグラフ構造は集団のグラフ構造の座位ごとに２ずつアレルを抽出したものに相当する。そして個体ごとに相の特定された遺伝子型（順列遺伝子型）を考えると、順列遺伝子型はホモの場合１種類、ヘテロの場合２種類存在する。よって個体ごとに生成されるグラフ構造は、下記の数１４で表されるｎ部完全グラフ構造となる。

辺の重みとしては、集団のグラフ構造のようなハプロタイプ頻度ではなく、ＥＭアルゴリズムによって求めた集団における２座位間でのディプロタイプ頻度（ディプロタイプ形の確率）を与える。

ノード（頂点）としてアレルではなく順列遺伝子型を採用することによって、個体ごとにディプロタイプ形を一元的に扱うことができ、個体の持つ遺伝子型データから、ある１つのディプロタイプ形を選択すると、グラフ構造ではｎ次完全グラフが１つ選択されたことになる。

図２に４座位のデータ構造から１つの４次完全グラフ構造を選択した場合を示す。

座位ｉの頂点をｖ_ｉ＝（ｖ_ｉ１，ｖ_ｉ２）と表すとすると、座位ｉの頂点ｖ_ｉ１と座位ｊの頂点ｖ_ｊ１とを結ぶ辺には個人がｖ_ｉ１とｖ_ｊ１を持つ同時確率ｐ（ｖ_ｉ１，ｖ_ｊ１）が付加されていると考えることができる。個人が実際に持つディプロタイプ形を構成するｎ次完全グラフでは、どの辺も高い同時確率を持ち、低い同時確率を持つ辺がグラフ内に入ることはないと予想される。この予想より、頂点集合Ｓ＝{ｖ_１’，ｖ_２’，・・・，ｖ_ｎ’}を選択したときの評価値（スコア）として下記の数１５を提案する。

個人indが持つ実際のディプロタイプ形Ｄ（ind）は、予想から最大の評価値をとるはずであるから、個人のディプロタイプ形は最大の評価値を持つｎ次完全グラフを探索する問題に置き換えることが可能である。その問題は下記の数１６で表される。

グラフ構造から個人のディプロタイプ形を推定するために、頂点の全ての組み合わせ探索する方法では、座位数ｎに対して時間計算量がＯ（２^ｎ）となってしまい、実行困難に陥ってしまうことが予想される。よって、問題を辺の持つ値を対数値とし、ｎ部完全グラフから辺重みの和が最大となるｎ次完全グラフを求めるグラフ問題に変換し、高速にハプロタイプを推定する手法が必要である。

座位ｉの順列遺伝子型に相当する頂点をｖ_ｉ、座位ｉと座位ｊの頂点を結ぶ辺をｅ_ｉｊとし、辺集合Ｅは判定時に使用する辺集合Ｅ_ａと使用しない辺集合Ｅ_ｐに分けられているものとする。辺ｅ_ｉｊはディプロタイプ頻度を持つため、ｅ_ｉ１ｊ１とｅ_ｉ２ｊ２、ｅ_ｉ１ｊ２とｅ_ｉ２ｊ１はそれぞれ同じ値を持つこととなる。尚、頂点ｖ_ｉと辺ｅ_ｉｊは、それぞれ、下記の数１７、数１８で表される。

図３に、以下に個人の持つディプロタイプ形を推定するための、数１４で表されるｎ部完全グラフ構造から辺の重みの和が最大となるｎ次完全グラフ構造を抽出するアルゴリズムを示す。

最初に、全てのｅ_ｉｊ（０≦ｊ≦ｎ−１、ｉ＋１≦ｊ≦ｎ）において（ｅ_ｉ１ｊ１，ｅ_ｉ２ｊ２）と（ｅ_ｉ１ｊ２，ｅ_ｉ２ｊ１）の値を比較し、大きい方をＥ_ａ、小さい方をＥ_ｐに追加する（Step-1）。

次に、頂点集合と辺集合Ｅ_ａで構成されるグラフ構造内にｎ次完全グラフが存在するかを判定する（Step-2）。詳述すると、まず、座位０と第２頂点のみを除き、全ての頂点をactiveにする（Step-2.1）。座位ｉと座位ｉ＋１において以下の操作を行う（初期値ｉ＝０）（Step-2.2）。

座位ｉ＋１の頂点のうち、座位ｉの頂点とＥ_ａの辺で接続されていないものはpassiveとする（Step-2．2.1）。ｊ＞ｉ＋１なるｖ_ｊにおいて、接続関係にあるｖ_ｉ，ｖ_ｉ＋１の両方に接続している頂点のみをactiveとし、残りをpassiveとする。その座位、ある座位のｖ_ｊが２つともpassiveとなったら、Step-3へと進む（Step-2．2.2）。ｉ＝ｉ＋１として、Step-2．2.2を繰り返す（Step-2．2.3）。

次に、ｉ＝ｎとし、Step-2.1の操作を逆順に行う。ｉ＝１まで終了したならばStep-4へ進む（Step-2.3）。

引続いて、Ｅ_ｐの辺の中で最大値を持つ辺のペアをＥ_ａに移動し、Step-2の判定を繰り返す（Step-3）。

最後に、残っている各座位の頂点集合から個人のディプロタイプ形を推定する。もし、同一座位の２頂点が残っている場合は、判定の繰り返し時に先に選択されたものを優先する（Step-4）。

このｎ次完全グラフ判定アルゴリズムにより、ＥＭアルゴリズムによるハプロタイプ推定が座位数ｎのＳＮＰデータに対して計算複雑度、空間複雑度ともＯ（２^ｎ）であるのに対し、本発明によるハプロタイプ推定方法では計算複雑度をＯ（ｎ^４）に、空間複雑度をＯ（ｎ^２）に削減している。

以上のデータ構造、アルゴリズムでは個体中の欠測を考慮せずに説明を行った。しかしながら、実際に得られる大規模な遺伝子型データでは、すべての座位において遺伝子型が観測できることは稀であり、各個体において少なからず遺伝子型が観測できない部位が存在することになる。実際に得られるデータに対応するためには、欠測した部分を補完して推定する必要があると考えられる。

ＥＭアルゴリズムによる推定では、欠測した座位において考えらえる全ての順列遺伝子型を仮定しハプロタイプ候補を上げることにより対応することが可能である。これに対して、本発明によるハプロタイプ推定方法では全ての順列遺伝子型を考慮すると個人の持つグラフ構造が数１４ではなくなり、辺の重みが最大となる完全グラフ構造の探索が複雑になる問題点が存在する。

しかしながら、もし入力座位数に対して欠測数が十分に少ない数であると仮定できるならば、欠測を除いた状態で推定した結果と欠測を含めて推定した結果が異なることは極めて稀であると推測される。よって本発明によるハプロタイプ推定方法では、欠測を除いた状態で一度推定を行い、得られた結果を用いて補完を行う方法を考える。

簡単のため座位数がｎであり、欠測数が１つであるとする。欠測を含まない座位ではldlightによる推定の結果から順列遺伝子型が決定しているとする。以上の仮定を図式化したものを図４に示す。

座位ｉを除いた座位では順列遺伝子型（グラフ中では頂点）が決定しているため、座位ｉに対して、下記の数１９のように、それぞれの同時確率（グラフ中では辺の重み）を得ることが可能である。

ldlightのアルゴリズムと同様に考えるならば、座位ｉとその他の各座位との同時確率は、推定結果が尤もらしならばどれも高くなり、１つでも極端に小さい値をとるならば推定結果が実際に個体が持つハプロタイプとは異なる可能性が高いと考えられる。以上から、欠測部分の順列遺伝子型を評価する評価値（スコア）score（Ｖ_ｉ）として、下記の数２０で表されるものを提案する。

入力データにＳＮＰを仮定すると、欠側座位において考えられる順列遺伝子型はアレル１，２に対して１１，１２，２１，２２の４種類存在することになる。この４種類の順列遺伝子型を示す頂点をそれぞれＶ_ｉ１，Ｖ_ｉ２，Ｖ_ｉ３，Ｖ_ｉ４とするならば、補完するべき頂点は、下記の数２１のように、上記数２０の最大値をとるものとなる。

入力データが、マルチアレルであればscoreの算出回数が増加することになるが、ＳＮＰデータの解析ならば１座位につき４回のscore算出のみで済むため、実行時間の面では大きく問題になることは無いと予想される。

次に、本発明によるハプロタイプ推定方法の有効性を示すために、本発明者らは、公開されている実データとシミュレーションデータの２種類のデータについて解析を行い、ＥＭアルゴリズムによる手法と精度および実行時間の比較を行った。特に実行時間に関しては、他の多遺伝子座のハプロタイプ推定に対応したソフトウェアと比較を行った。

本発明者らは、非特許文献３２に開示されているヒト染色体5q31領域のＩＢＤに関するＳＮＰデータを用いて実行時間の測定を行った。その測定結果を図５に示す。ＥＭアルゴリズムを用いたソフトウェアとして、非特許文献２２に記載されているLDSUPPORTを使用している。

図５から少ない座位数ではＥＭアルゴリズムの高速性が顕著に表れているが、遺伝子座数１０以上のデータでは、本発明に係るハプロタイプ推定方法（idlight）の実行時間が常に短い結果となった。１８座位のデータの処理時間の比較では約１３５０倍の高速性を実現している。

推定精度の比較実験として、本発明者らは、報告されているハプロタイプブロックのうち、最初の５ブロックに相当する部分のＳＮＰデータを使用して、そのうちの１００人分の解析を行った。集団をcase集団（罹患者集団）とcontrol集団（非罹患者集団）の２種類に分類し、そのそれぞれについて本発明に係るハプロタイプ推定方法とＥＭアルゴリズムによって推定したハプロタイプ頻度を比較した。その比較結果を下記の表１に示す。

各ブロックにおいて集団内に１０％以上存在するものをメジャーハプロタイプとして列挙し、その他のマイナーハプロタイプはothersとした。表１から、ブロックにおいては多少の差が生じているものの、ＥＭアルゴリズムによる推定と比較しても遜色の無い推定が実現できていると考えられる。

次に、本発明者らは、本発明に係るハプロタイプ推定方法と従来の方法について、シミュレーションデータに対する実行時間及び推定精度の比較を行った。実行時間を比較するために、本発明者らは、染色体5q31領域のハプロタイプブロックのうち、3-7ブロックに相当する58遺伝子座のハプロタイプを用い、その中からサンプリングを行うことによってデータ生成を行った。また、多遺伝子座のハプロタイプを推定するためのソフトウェアとして、非特許文献２４に開示されているPHASEおよび非特許文献２５に開示されているPL-EMを比較対象に加えた。個体数を１００として、遺伝子座数を５から５０まで変化させたときの実行時間を図６に示す。

図６から、シミュレーションデータの場合、１０遺伝子座以下ではオーバヘッドの少ないＥＭアルゴリズムによる推定の方が高速に実行できる。これは少ない遺伝子座ではとりうるハプロタイプ数が少なくなり、ＥＭアルゴリズムのE-stepが高速に実行されることによるものである。それに対してldlightは、座位間の重みを求めて個体ごとにグラフ構造を構成するため、少ない座位ではそのオーバーヘッドが大きく出る。

一方、１０遺伝子座より多い遺伝子座においては、本発明に係るハプロタイプ推定方法の効果が現れていることが分かる。これは遺伝子座が多くなると、とりうるハプロタイプが指数関数的に増大するためＥＭアルゴリズムは実行速度が指数関数的に増加するが、ldlightはグラフ構造を作った後の解析が非常に高速であるため、遺伝子座数１５以降ではその差が顕著に表れている。

また、idlightと同じく多遺伝子座を扱うことができるPHASE、PL-EMと比較すると、５０遺伝子座においておよそ１００倍近くの高速化が達成されていることが分かる。これは、従来の高速化手法はＥＭアルゴリズムで考慮されていた集団全体のハプロタイプに対して、それを減らすというアプローチであったため、本質的にはＥＭアルゴリズムと同様の問題を抱えているからである。それに対して、本発明に係るハプロタイプ推定方法では、新しいデータ構造を用いることによって、集団の取り得るハプロタイプを考慮せず、完全グラフが発見されたときに初めてハプロタイプとして構成されるため、根本的に異なるアプローチである。

本発明に係るハプロタイプ推定方法では１００遺伝子座以上のデータも高速に処理可能であり、１００個体１００座位の遺伝子型データからハプロタイプを推定した場合では、Celeron８００ＭＨｚにおいて約４秒、６．６ＭByteのメモリ使用量で処理が可能である。

次に精度に関する比較を行うために、本発明者らは、coalescent with recombinationにもとづくシミュレーションによって、擬似的にＳＮＰデータを生成して精度の比較を行った。個体数を１００に固定して、遺伝子座数を２から１６まで２ずつ変化させたときの平均エラー率を図７および下記の表２に示す。

平均エラー率は１０００個のデータセットを解析したときのエラー率の平均値である。

図７および表２より、idlightは、遺伝子座数８においてidsupportよりエラー率が高くなったが、その他の遺伝子座数では全てのidsupportよりも精度の高い結果が得られた。遺伝子座数２の場合では、idlightは全ての個体に対して正しいハプロタイプを推定することができた。

すなわち、本発明に係るハプロタイプ推定方法は、精度に関しては、５q３１領域の解析においてＥＭアルゴリズムとほぼ同様の推定結果を示している。coalescent-based modelにおいては、本発明に係るハプロタイプ推定方法は、遺伝子座数１４では０．２以下のエラー率で解析可能となっている。ＥＭアルゴリズムと比較してもこれは遜色のない結果で、本発明に係るハプロタイプ推定方法は、平均的にはＥＭアルゴリズムよりよい結果を示している。coalescent-based modelは集団遺伝学の理論に基づいてシミュレーションデータを生成しているため、実際の集団を解析した場合についてもある程度高い精度を示すことが期待できる。従って、本発明に係るハプロタイプ推定方法は、精度の面でもＥＭアルゴリズムの代替となり得ることが示された。

このように、本発明に係るハプロタイプ推定方法は、ＥＭアルゴリズムを用いたアルゴリズムとの実行時間、推定精度の比較により、ＥＭアルゴリズムの持つ推定精度と同等の精度を確保し、大幅な実行時間の削減を実現していることが確認された。

次に、図８を参照して、本発明によるハプロタイプ推定方法を実現するコンピュータシステムについて説明する。

図８に示されるように、図示のコンピュータシステムは、ＣＰＵ１１と、主記憶装置１２と、入力装置１３と、出力装置１４と、表示装置１５と、記録媒体１６とを備え、これらは中央伝送路１７によって互いに接続されている。

入力装置１３は、データを入力するための装置で、例えば、キーボードやマウスなどのポインティングデバイス等から成る。尚、入力するデータとしては、遺伝子型データがある。出力装置１４はＣＰＵ１１の実行により得られた結果を出力するための装置で、例えば、プリンタ等から成る。表示装置１５はＣＰＵ１１の実行に得られた結果を表示するための装置で、例えば、ディスプレイ等から成る。記録媒体１６はハードディスク等からなる。

上述した本発明によるハプロタイプ推定方法を実施するためのプログラム(以下、本発明のプログラムともいう。）と入力データとは、記録媒体１６に保存されている。ＣＰＵ１１は、入力装置１３からの信号を受け、本発明のプログラムを主記憶装置１２に呼び出し、実行する。

ＣＰＵ１１は、入力データを主記憶装置１２に呼び出し、本発明のプログラムを実行する。得られた結果であるハプロタイプ推定結果を表示装置１５に表示し、記録媒体１６に保存する。得られた結果は、出力装置１４から出力することもできる。

図９は、本発明のプログラムを説明するためのフローチャートである。

最初に、ＣＰＵ１１は、集団に関するグラフを作成する（ステップＳ１）。すなわち、個体のディプロタイプ形を、順列遺伝子型を頂点（ノード）に対応させたｎ部完全グラフ構造で表す。ここで、ノード（頂点）には順列遺伝子形が対応する。エッジ（辺）にはＥＭアルゴリズムによって求めたディプロタイプ頻度（ディプロタイプ形の確率）を付加する。

次に、ＣＰＵ１１は、集団に関するグラフより、未探索の個体に関するグラフを抽出する（ステップＳ２）。そして、ＣＰＵ１１は、個体に関するグラフより枝を抜き出して部分グラフを作成する（ステップＳ３）。このとき、ＣＰＵ１１は、２遺伝子座間でスコアの高い方の枝を抜き出す。ＣＰＵ１１は、部分グラフが完全グラフを含むか否かを判断する（ステップＳ４）。部分グラフが完全グラフを含まない場合（ステップＳ４のｎｏ）、ＣＰＵ１１は、選ばれていない枝のうち、最もスコアの高いものを部分グラフに加え（ステップＳ５）、ステップＳ４に戻る。

一方、部分グラフが完全グラフを含むなら（ステップＳ４のｙｅｓ）、ＣＰＵ１１は、完全グラフに対応するディプロタイプ形を得る（ステップＳ６）。引続いて、ＣＰＵ１１は、欠測部分の推定を行う（ステップＳ７）。ＣＰＵ１１は、未探索の個体があるか否かを判断する（ステップＳ８）。未探索の個体がある場合（ステップＳ８のｙｅｓ）、ＣＰＵ１１はステップＳ２に戻る。

未探索の個体がない場合（ステップＳ８のｎｏ）、ＣＰＵ１１は全ての個体のディプロタイプ形より各ハプロタイプの出現回数をカウントし、集団のハプロタイプ頻度を求める（ステップＳ９）。

詳述すると、集団のハプロタイプ頻度は、集団内の各個体のディプロタイプ形をそれぞれ推定し、ハプロタイプ毎に集計して割合を求めることによって得られる。例えば、集団の個体数がＮの場合を考える。推定された各個体のディプロタイプ形の中に、あるハプロタイプがａ本存在すれば、そのハプロタイプ頻度はａ／（２×Ｎ）となる。

図９に示されたプログラムを図８に示されたコンピュータシステムが実行することで、本発明によるハプロタイプ推定装置が実現される。すなわち、本発明によるハプロタイプ推定装置は、座位数がｎの遺伝子型データからハプロタイプを推定する装置である。ハプロタイプ推定装置は、個体のディプロタイプ形を、順列遺伝子型を頂点に対応させたｎ部完全グラフ構造として表す手段と、ｎ部完全グラフ構造の各辺の重みを、遺伝子型データから最尤推定によって推定したディプロタイプ形の確率とする手段と、所定のスコアを用いて、ｎ部完全グラフ構造の頂点集合を選択する手段と、各個体のディプロタイプ形を、最大スコアを持つｎ次完全グラフ探索問題に置き換えて推定する手段と、この推定により得られた各個体のディプロタイプ形より集団のハプロタイプ頻度を求める手段とを含む。

このハプロタイプ推定装置において、最尤推定としてＥＭアルゴリズムを用いることができる。また、座位ｉ（ｉは１からｎ−１までの整数）の頂点をｖ_ｉ＝{ｖ_ｉ１，ｖ_ｉ２}と表し、座位ｉの頂点ｖ_ｉ１と座位ｊ（ｊはｉ＋１からｎまでの整数）の頂点ｖ_ｊ１を結ぶ辺に個体が頂点ｖ_ｉ１と頂点ｖ_ｊ１を持つ同時確率ｐ（ｖ_ｉ１，ｖ_ｊ１）が付加されていると考えた場合に、スコアscore（Ｓ）として上記数１５を用いることができる。

また、本発明によるハプロタイプ推定装置は、遺伝子型データとして欠測を含むデータを扱う個体については、欠測を除いてディプロタイプ形を推定し、その結果を用いて補完を行う手段を更に含んでも良い。その場合、欠測した座位ｉに対してのそれぞれの同時確率をＰ（Ｖ_１，Ｖ_ｉ），Ｐ（Ｖ_２，Ｖ_ｉ），…，Ｐ（Ｖ_ｎ，Ｖ_ｉ）としたとき、欠側部分の順列遺伝子型のスコアscore（Ｖ_ｉ）として上記数２０を用いることができる。

以上、本発明について好ましい実施の形態によって説明してきたが、本発明は上述した実施の形態に限定しないのは勿論である。例えば、上述した実施の形態では、多型マーカがＳＮＰである場合を例に挙げて説明したが、アレルが複数あるマイクロサテライト等の他の多型マーカに対しても当業者であれば容易に適用可能であることは明らかである。また、最尤推定としてＥＭアルゴリズムを用いた場合についてのみ説明しているが、他の最尤推定を用いても良いのは勿論である。更に、個体が取りうるディプロタイプ形の確率を推定することもできる。

グラフのデータ構造を示す図である。４座位のデータ構造から１つの４次完全グラフ構造を選択した場合の例を示す図である。個人の持つディプロタイプ形を推定するためのｎ部完全グラフ構造から辺の重みの和が最大となるｎ次完全グラフ構造を抽出するアルゴリズムを示す図である。欠測した部分を補完して推定する例を説明するための図である。ヒト染色体５q３１領域のＩＢＤに関するＳＮＰデータを用いて実行時間の測定を行った結果を示す図である。本発明に係るハプロタイプ推定方法と従来方法について、シミュレーションデータに対する実行時間を比較するために、個体数を１００として、遺伝子座数を５から５０まで変化させたときの実行時間を示す図である。本発明に係るハプロタイプ推定方法と従来方法について、シミュレーションデータに対する推定精度を比較するために、個体数を１００に固定して、遺伝子座数を２から１６まで２ずつ変化させたときの平均エラー率を示す図である。本発明によるハプロタイプ推定方法を実現するコンピュータシステムのハードウェア構成を示すブロック図である。本発明によるハプロタイプ推定方法を図８に示すコンピュータシステムで実施するためのプログラムを説明するためのフローチャートである。

符号の説明

１１ＣＰＵ
１２主記憶装置
１３入力装置
１４出力装置
１５表示装置
１６記録媒体
１７中央伝送路

Claims

コンピュータシステムを用いて、座位数がｎの遺伝子型データから各個体のディプロタイプ形を推定する方法であって、
個体のディプロタイプ形を、順列遺伝子型を頂点に対応させたｎ部完全グラフ構造として表すステップと、
前記ｎ部完全グラフ構造の各辺の重みを、前記遺伝子型データから最尤推定によって推定したディプロタイプ形の確率とするステップと、
所定のスコアを用いて、前記ｎ部完全グラフ構造の頂点集合を選択するステップと、
各個体のディプロタイプ形を、最大スコアを持つｎ次完全グラフ探索問題に置き換えて推定するステップと
を含むディプロタイプ形の推定方法。
前記最尤推定としてＥＭアルゴリズムを用いることを特徴とする請求項１に記載のディプロタイプ形の推定方法。
座位ｉ（ｉは１からｎ−１までの整数）の頂点をｖ_ｉ＝{ｖ_ｉ１，ｖ_ｉ２}と表し、前記座位ｉの頂点ｖ_ｉ１と座位ｊ（ｊはｉ＋１からｎまでの整数）の頂点ｖ_ｊ１を結ぶ辺に個体が前記頂点ｖ_ｉ１と前記頂点ｖ_ｊ１を持つ同時確率ｐ（ｖ_ｉ１，ｖ_ｊ１）が付加されていると考えた場合に、前記スコアscore（Ｓ）として下記の数１

を用いることを特徴とする請求項１又は２に記載のディプロタイプ形の推定方法。
前記遺伝子型データとして欠測を含むデータを扱う個体については、前記欠測を除いて前記ディプロタイプ形を推定し、その結果を用いて補完を行うステップを更に含む、請求項１に記載のディプロタイプ形の推定方法。
欠測した座位ｉに対してのそれぞれの同時確率をＰ（Ｖ_１，Ｖ_ｉ），Ｐ（Ｖ_２，Ｖ_ｉ），…，Ｐ（Ｖ_ｎ，Ｖ_ｉ）としたとき、前記欠側部分の順列遺伝子型のスコアscore（Ｖ_ｉ）として下記の数２

を用いることを特徴とする請求項４に記載のディプロタイプ形の推定方法。
コンピュータシステムを用いて、座位数がｎの遺伝子型データからハプロタイプを推定する方法であって、
個体のディプロタイプ形を、順列遺伝子型を頂点に対応させたｎ部完全グラフ構造として表すステップと、
前記ｎ部完全グラフ構造の各辺の重みを、前記遺伝子型データから最尤推定によって推定したディプロタイプ形の確率とするステップと、
所定のスコアを用いて、前記ｎ部完全グラフ構造の頂点集合を選択するステップと、
各個体のディプロタイプ形を、最大スコアを持つｎ次完全グラフ探索問題に置き換えて推定するステップと、
該推定により得られた各個体のディプロタイプ形より集団のハプロタイプ頻度を求めるステップと
を含むハプロタイプ推定方法。
前記最尤推定としてＥＭアルゴリズムを用いることを特徴とする請求項６に記載のハプロタイプ推定方法。
座位ｉ（ｉは１からｎ−１までの整数）の頂点をｖ_ｉ＝{ｖ_ｉ１，ｖ_ｉ２}と表し、前記座位ｉの頂点ｖ_ｉ１と座位ｊ（ｊはｉ＋１からｎまでの整数）の頂点ｖ_ｊ１を結ぶ辺に個体が前記頂点ｖ_ｉ１と前記頂点ｖ_ｊ１を持つ同時確率ｐ（ｖ_ｉ１，ｖ_ｊ１）が付加されていると考えた場合に、前記スコアscore（Ｓ）として下記の数３

を用いることを特徴とする請求項６又は７に記載のハプロタイプ推定方法。
前記遺伝子型データとして欠測を含むデータを扱う個体については、前記欠測を除いて前記ディプロタイプ形を推定し、その結果を用いて補完を行うステップを更に含む、請求項６に記載のハプロタイプ推定方法。
欠測した座位ｉに対してのそれぞれの同時確率をＰ（Ｖ_１，Ｖ_ｉ），Ｐ（Ｖ_２，Ｖ_ｉ），…，Ｐ（Ｖ_ｎ，Ｖ_ｉ）としたとき、前記欠側部分の順列遺伝子型のスコアscore（Ｖ_ｉ）として下記の数４

を用いることを特徴とする請求項９に記載のハプロタイプ推定方法。
座位数がｎの遺伝子型データから各個体のディプロタイプ形を推定する装置であって、
個体のディプロタイプ形を、順列遺伝子型を頂点に対応させたｎ部完全グラフ構造として表す手段と、
前記ｎ部完全グラフ構造の各辺の重みを、前記遺伝子型データから最尤推定によって推定したディプロタイプ形の確率とする手段と、
所定のスコアを用いて、前記ｎ部完全グラフ構造の頂点集合を選択する手段と、
各個体のディプロタイプ形を、最大スコアを持つｎ次完全グラフ探索問題に置き換えて推定する手段と
を含むディプロタイプ形推定装置。
前記最尤推定としてＥＭアルゴリズムを用いることを特徴とする請求項１１に記載のディプロタイプ形推定装置。
座位ｉ（ｉは１からｎ−１までの整数）の頂点をｖ_ｉ＝{ｖ_ｉ１，ｖ_ｉ２}と表し、前記座位ｉの頂点ｖ_ｉ１と座位ｊ（ｊはｉ＋１からｎまでの整数）の頂点ｖ_ｊ１を結ぶ辺に個体が前記頂点ｖ_ｉ１と前記頂点ｖ_ｊ１を持つ同時確率ｐ（ｖ_ｉ１，ｖ_ｊ１）が付加されていると考えた場合に、前記スコアscore（Ｓ）として下記の数５

を用いることを特徴とする請求項１１又は１２に記載のディプロタイプ形推定装置。
前記遺伝子型データとして欠測を含むデータを扱う個体については、前記欠測を除いて前記ディプロタイプ形を推定し、その結果を用いて補完を行う手段を更に含む、請求項１１に記載のディプロタイプ形推定装置。
欠測した座位ｉに対してのそれぞれの同時確率をＰ（Ｖ_１，Ｖ_ｉ），Ｐ（Ｖ_２，Ｖ_ｉ），…，Ｐ（Ｖ_ｎ，Ｖ_ｉ）としたとき、前記欠側部分の順列遺伝子型のスコアscore（Ｖ_ｉ）として下記の数６

を用いることを特徴とする請求項１４に記載のディプロタイプ形推定装置。
座位数がｎの遺伝子型データからハプロタイプを推定する装置であって、
個体のディプロタイプ形を、順列遺伝子型を頂点に対応させたｎ部完全グラフ構造として表す手段と、
前記ｎ部完全グラフ構造の各辺の重みを、前記遺伝子型データから最尤推定によって推定したディプロタイプ形の確率とする手段と、
所定のスコアを用いて、前記ｎ部完全グラフ構造の頂点集合を選択する手段と、
各個体のディプロタイプ形を、最大スコアを持つｎ次完全グラフ探索問題に置き換えて推定する手段と、
該推定により得られた各個体のディプロタイプ形より集団のハプロタイプ頻度を求める手段と
を含むハプロタイプ推定装置。
前記最尤推定としてＥＭアルゴリズムを用いることを特徴とする請求項１６に記載のハプロタイプ推定装置。
座位ｉ（ｉは１からｎ−１までの整数）の頂点をｖ_ｉ＝{ｖ_ｉ１，ｖ_ｉ２}と表し、前記座位ｉの頂点ｖ_ｉ１と座位ｊ（ｊはｉ＋１からｎまでの整数）の頂点ｖ_ｊ１を結ぶ辺に個体が前記頂点ｖ_ｉ１と前記頂点ｖ_ｊ１を持つ同時確率ｐ（ｖ_ｉ１，ｖ_ｊ１）が付加されていると考えた場合に、前記スコアscore（Ｓ）として下記の数７

を用いることを特徴とする請求項１６又は１７に記載のハプロタイプ推定装置。
前記遺伝子型データとして欠測を含むデータを扱う個体については、前記欠測を除いて前記ディプロタイプ形を推定し、その結果を用いて補完を行う手段を更に含む、請求項１６に記載のハプロタイプ推定装置。
欠測した座位ｉに対してのそれぞれの同時確率をＰ（Ｖ_１，Ｖ_ｉ），Ｐ（Ｖ_２，Ｖ_ｉ），…，Ｐ（Ｖ_ｎ，Ｖ_ｉ）としたとき、前記欠側部分の順列遺伝子型のスコアscore（Ｖ_ｉ）として下記の数８

を用いることを特徴とする請求項１９に記載のハプロタイプ推定装置。
コンピュータに、座位数がｎの遺伝子型データから各個体のディプロタイプ形を推定させるためのプログラムであって、前記コンピュータに、
個体のディプロタイプ形を、順列遺伝子型を頂点に対応させたｎ部完全グラフ構造として表す機能と、
前記ｎ部完全グラフ構造の各辺の重みを、前記遺伝子型データから最尤推定によって推定したディプロタイプ形の確率とする機能と、
所定のスコアを用いて、前記ｎ部完全グラフ構造の頂点集合を選択する機能と、
各個体のディプロタイプ形を、最大スコアを持つｎ次完全グラフ探索問題に置き換えて推定する機能と
を実現させるためのプログラム。
前記最尤推定としてＥＭアルゴリズムを用いることを特徴とする請求項２１に記載のプログラム。
座位ｉ（ｉは１からｎ−１までの整数）の頂点をｖ_ｉ＝{ｖ_ｉ１，ｖ_ｉ２}と表し、前記座位ｉの頂点ｖ_ｉ１と座位ｊ（ｊはｉ＋１からｎまでの整数）の頂点ｖ_ｊ１を結ぶ辺に個体が前記頂点ｖ_ｉ１と前記頂点ｖ_ｊ１を持つ同時確率ｐ（ｖ_ｉ１，ｖ_ｊ１）が付加されていると考えた場合に、前記スコアscore（Ｓ）として下記の数９

を用いることを特徴とする請求項２１又は２２に記載のプログラム。
前記コンピュータに、前記遺伝子型データとして欠測を含むデータを扱う個体については、前記欠測を除いて前記ディプロタイプ形を推定し、その結果を用いて補完を行う機能を更に実現させる、請求項２１に記載のプログラム。
欠測した座位ｉに対してのそれぞれの同時確率をＰ（Ｖ_１，Ｖ_ｉ），Ｐ（Ｖ_２，Ｖ_ｉ），…，Ｐ（Ｖ_ｎ，Ｖ_ｉ）としたとき、前記欠側部分の順列遺伝子型のスコアscore（Ｖ_ｉ）として下記の数１０

を用いることを特徴とする請求項２４に記載のプログラム。
コンピュータに、座位数がｎの遺伝子型データからハプロタイプを推定させるためのプログラムであって、前記コンピュータに、
個体のディプロタイプ形を、順列遺伝子型を頂点に対応させたｎ部完全グラフ構造として表す機能と、
前記ｎ部完全グラフ構造の各辺の重みを、前記遺伝子型データから最尤推定によって推定したディプロタイプ形の確率とする機能と、
所定のスコアを用いて、前記ｎ部完全グラフ構造の頂点集合を選択する機能と、
各個体のディプロタイプ形を、最大スコアを持つｎ次完全グラフ探索問題に置き換えて推定する機能と、
該推定により得られた各個体のディプロタイプ形より集団のハプロタイプ頻度を求める機能と
を実現させるためのプログラム。
前記最尤推定としてＥＭアルゴリズムを用いることを特徴とする請求項２６に記載のプログラム。
座位ｉ（ｉは１からｎ−１までの整数）の頂点をｖ_ｉ＝{ｖ_ｉ１，ｖ_ｉ２}と表し、前記座位ｉの頂点ｖ_ｉ１と座位ｊ（ｊはｉ＋１からｎまでの整数）の頂点ｖ_ｊ１を結ぶ辺に個体が前記頂点ｖ_ｉ１と前記頂点ｖ_ｊ１を持つ同時確率ｐ（ｖ_ｉ１，ｖ_ｊ１）が付加されていると考えた場合に、前記スコアscore（Ｓ）として下記の数１１

を用いることを特徴とする請求項２６又は２７に記載のプログラム。
前記コンピュータに、前記遺伝子型データとして欠測を含むデータを扱う個体については、前記欠測を除いて前記ディプロタイプ形を推定し、その結果を用いて補完を行う機能を更に実現させる、請求項２６に記載のプログラム。
欠測した座位ｉに対してのそれぞれの同時確率をＰ（Ｖ_１，Ｖ_ｉ），Ｐ（Ｖ_２，Ｖ_ｉ），…，Ｐ（Ｖ_ｎ，Ｖ_ｉ）としたとき、前記欠側部分の順列遺伝子型のスコアscore（Ｖ_ｉ）として下記の数１２

を用いることを特徴とする請求項２９に記載のプログラム。