JP2004192018A - Ｄｎａプールによるハプロタイプ頻度推定方法 - Google Patents
Ｄｎａプールによるハプロタイプ頻度推定方法 Download PDFInfo
- Publication number
- JP2004192018A JP2004192018A JP2002303395A JP2002303395A JP2004192018A JP 2004192018 A JP2004192018 A JP 2004192018A JP 2002303395 A JP2002303395 A JP 2002303395A JP 2002303395 A JP2002303395 A JP 2002303395A JP 2004192018 A JP2004192018 A JP 2004192018A
- Authority
- JP
- Japan
- Prior art keywords
- haplotype
- genotype
- frequency
- pool
- algorithm
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
Images
Landscapes
- Information Retrieval, Db Structures And Fs Structures Therefor (AREA)
- Complex Calculations (AREA)
Abstract
【課題】ハプロタイプ推定において、処理速度等のコンピュータにおけるスペックに依存せずに情報処理する。
【解決手段】集団に含まれる複数の検体に関する遺伝子型情報を集積してある遺伝子型プール情報を入力値として、期待値最大化アルゴリズムにより上記集団におけるハプロタイプ頻度を推定する。
【選択図】 図５
【解決手段】集団に含まれる複数の検体に関する遺伝子型情報を集積してある遺伝子型プール情報を入力値として、期待値最大化アルゴリズムにより上記集団におけるハプロタイプ頻度を推定する。
【選択図】 図５
Description
【０００１】
【発明の属する技術分野】
本発明は、最大尤度推定アルゴリズムとして知られている期待値最大化アルゴリズム（以下、ＥＭアルゴリズムと呼ぶ）を用いたハプロタイプ頻度推定方法、ハプロタイプ頻度推定アルゴリズム及びハプロタイプ頻度推定プログラムに関する。
【０００２】
【従来の技術】
遺伝統計学的手法とは、不明の疾患遺伝子を見出す方法の一つであり、個体間の遺伝子多型の違いと、各個体の形質データとから統計的な計算によって、形質に関連する遺伝子座を探す方法である。すなわち、遺伝統計学的手法によれば、遺伝的多型と形質（表現型、症状等）の情報のみから統計学を用いて原因（遺伝子座）を解明することができる。
【０００３】
遺伝統計学的手法としては、主としてパラメトリック連鎖解析とノンパラメトリック連鎖解析とに分類される。「連鎖解析」とは、メンデルの第３法則（独立の法則）の例外現象である連鎖を利用して表現型に関連する遺伝子座の染色体上に、疾患原因遺伝子座の存在領域を絞り込んでいく（位置をマッピングする）遺伝統計学的手法である。「パラメトリック」とは、解析に先立ち、遺伝子型に応じて固有の発症確率（＝浸透率）を仮定する（遺伝形式を仮定する）ことを意味する。遺伝子型における浸透率を仮定するということは、すなわち疾患原因遺伝子座位での遺伝子型によって発症確率が記述できるということである。浸透率が仮定できるときには、表現型（罹患の有無）からその疾患の原因遺伝子座での遺伝子型を簡単に関連づけることができる。「ノンパラメトリック」とは、遺伝形式を仮定せずに解析を行うことを意味する。
【０００４】
ところで、ハプロタイプとは、複数の遺伝子座における各対立遺伝子の組合せとして定義される。表現型を遺伝子領域にマッピングする研究など、ハプロタイプ推定は、多くの遺伝子研究において重要な意味を持っている。通常、集団からサンプリングされた複数個人の遺伝子型データに基づいて集団のハプロタイプ推定及び個人のディプロタイプ推定が行われている。
【０００５】
【発明が解決しようとする課題】
ところで、複数個人の遺伝子型データを用いた場合、例えば、最大尤度推定を行うためのＥＭアルゴリズムを使用して集団のハプロタイプ推定及び個人のディプロタイプ推定を行うことができる。しかしながら、この場合、複数個人の遺伝子型データを用いているため、入力するためのデータ量が膨大であり、現状使用されている高スペックのコンピュータでさえも処理することが困難であった。
【０００６】
そこで、本発明は、上述したような実状に鑑み、処理速度等のコンピュータにおけるスペックに依存せずに情報処理が可能なＥＭアルゴリズムを用いたハプロタイプ推定方法、ハプロタイプ頻度推定アルゴリズム及びハプロタイプ頻度推定プログラムを提供することを目的とする。
【０００７】
【課題を解決するための手段】
上述した目的を達成するため、本発明者が鋭意検討した結果、複数個人の遺伝子型データの代わりに、複数個人の遺伝子型データをプールしたデータを用いることにより、集団のハプロタイプ頻度推定に要する処理能力を低く抑えられ、また、複数個人の遺伝子型データをプールしたデータを用いても集団のハプロタイプ頻度推定を正確に行えることを見出し、本発明を完成するに至った。
【０００８】
すなわち、本発明は、以下を包含する。
（１）集団に含まれる複数の検体に関する遺伝子型情報を集積してある遺伝子型プール情報を入力値として、期待値最大化アルゴリズムにより上記集団におけるハプロタイプ頻度を推定する方法。
【０００９】
（２）上記複数の検体に由来するＤＮＡを含む複数のＤＮＡプールを調製し、当該複数のＤＮＡプールに含まれるＤＮＡにおける所定の遺伝子座に関するタイピングを行う工程を含むことを特徴とする（１）記載のハプロタイプ頻度を推定する方法。
【００１０】
（３）２以上６以下の検体に関する遺伝子型情報を集積することして遺伝子型プール情報を作成することを特徴とする（１）記載のハプロタイプ頻度を推定する方法。
【００１１】
（４）前記期待値最大化アルゴリズムはＥＭアルゴリズムであり、前記遺伝子型プール情報に最尤推定を行うことによりハプロタイプ頻度を推定することを特徴とする（１）記載のハプロタイプ頻度を推定する方法。
【００１２】
（５）集団に含まれる複数の検体に関する遺伝子型情報を集積してある遺伝子型プール情報を入力値として、期待値最大化アルゴリズムにより上記集団におけるハプロタイプ頻度を推定するアルゴリズム。
【００１３】
（６）上記複数の検体に由来するＤＮＡを含む複数のＤＮＡプールを調製し、当該複数のＤＮＡプールに含まれるＤＮＡにおける所定の遺伝子座に関するタイピングを行う工程を含むことを特徴とする（５）記載のハプロタイプ頻度を推定するアルゴリズム。
【００１４】
（７）２以上６以下の検体に関する遺伝子型情報を集積することして遺伝子型プール情報を作成することを特徴とする（５）記載のハプロタイプ頻度を推定するアルゴリズム。
【００１５】
（８）前記期待値最大化アルゴリズムはＥＭアルゴリズムであり、前記遺伝子型プール情報に最尤推定を行うことによりハプロタイプ頻度を推定することを特徴とする（５）記載のハプロタイプ頻度を推定するアルゴリズム。
【００１６】
（９）集団に含まれる複数の検体に関する遺伝子型情報を集積してある遺伝子型プール情報を入力値として、期待値最大化アルゴリズムにより上記集団におけるハプロタイプ頻度を推定するプログラム。
【００１７】
（１０）上記複数の検体に由来するＤＮＡを含む複数のＤＮＡプールを調製し、当該複数のＤＮＡプールに含まれるＤＮＡにおける所定の遺伝子座に関するタイピングを行う工程を含むことを特徴とする（９）記載のハプロタイプ頻度を推定するプログラム。
【００１８】
（１１）２以上６以下の検体に関する遺伝子型情報を集積することして遺伝子型プール情報を作成することを特徴とする（９）記載のハプロタイプ頻度を推定するプログラム。
【００１９】
（１２）前記期待値最大化アルゴリズムはＥＭアルゴリズムであり、前記遺伝子型プール情報に最尤推定を行うことによりハプロタイプ頻度を推定することを特徴とする（９）記載のハプロタイプ頻度を推定するプログラム。
【００２０】
【発明の実施の形態】
以下、本発明に係るＥＭアルゴリズムを用いたハプロタイプ頻度推定方法（以下、単に「ハプロタイプ頻度推定方法」と呼ぶ）を詳細に説明する。
【００２１】
ハプロタイプ頻度推定方法は、集団における複数の遺伝子型情報を集積（プール）した遺伝子型プール情報から、該集団におけるハプロタイプ頻度を算出する工程を含んでいる。
【００２２】
「遺伝子型」とは、所定の遺伝子座における一対の対立遺伝子の組み合わせ（常染色体における遺伝子座の場合）である。一対の対立遺伝子は、父性由来の一方の対立遺伝子と母性由来の他方の対立遺伝子とからなる。遺伝子型は、例えば所定の遺伝子座に多型部位が存在する場合、該多型部位に相当する一対の対立遺伝子についての遺伝子の個人についての特徴を示す。
【００２３】
「遺伝子型プール情報」とは、個人に関する遺伝子型情報を、複数の個人について集積した情報を意味する。複数の個人とは、特に限定されないが、２〜６人とすることが好ましく、２〜４人とすることがより好ましい。また、遺伝子型プール情報は、複数の遺伝子型情報を、複数の個人について集積した情報であってもよい。
【００２４】
遺伝子型プール情報となる遺伝子の特徴としては、例えば１塩基多型、マイクロサテライト及び挿入／欠失等の遺伝子変異を挙げることができる。遺伝子型に関する情報には、前記遺伝子変異のうち、１つの種類に限定されず、複数の種類の遺伝子変異が含まれていてもよい。
【００２５】
遺伝子型情報に複数の遺伝子座が含まれる場合、これら遺伝子座は連鎖していることが好ましい。従って、本発明においては、予め連鎖していることが知られている遺伝子座について遺伝子型プール情報を取得することが望ましいが、連鎖の有無が不明な場合もあり、必ずしもこの限りではない。
【００２６】
また、「ハプロタイプ」とは、複数の遺伝子座における各対立遺伝子の組合せを意味する。ハプロタイプは、複数の１塩基置換（ＳＮＰ； ｓｉｎｇｌｅ ｎｕｃｌｅｏｔｉｄｅ ｐｏｌｙｍｏｒｐｈｉｓｍ ）、マイクロサテライト、挿入／欠失等の多型部位が近接して存在し、これらが連鎖する場合に形成されるものである。また、常染色体上の遺伝子部位について、個人は、それぞれ父性或いは母性に由来する一対のハプロタイプを保有する。このような個人が保有する一対のハプロタイプの組み合わせを「ディプロタイプ」という。
【００２７】
「ハプロタイプ頻度」とは、複数の遺伝子座における各対立遺伝子について所定の組合せ（ハプロタイプ）が出現する頻度である。ハプロタイプは、複数の遺伝子座が狭い遺伝子領域に含まれる場合には世代交代によっても変化することは少ない。これは、複数の多型部位が存在する遺伝子領域内に交差が起きることによりハプロタイプは変化するが、当該遺伝子領域が狭い場合、交差が起きる確率は極めて低い（１世代あたり、１００，０００ｋｂに１回に過ぎない）ためである。そのため特定の集団内、例えば日本人では、それぞれのハプロタイプ頻度が決まっている場合が多い。ハプロタイプ頻度は、集団ごとに異なる場合が多い。
【００２８】
「集団」とは、遺伝子型情報を提供する個人の集まりを意味する。集団の人数は、算出されるハプロタイプ頻度の信頼性を高めるために、なるべく多くの人数から構成されることが望ましい。集団を構成する人数としては、２５名以上であることが好ましいが、特に制限されない。また、集団は均一（遺伝的背景が近接）であることが好ましい。例えば、ハプロタイプ頻度の推定を行う個人が日本人であるならば、遺伝子型情報を取得すべき「集団」は、日本人から構成される集団であることが好ましい。集団が均一か否かについては、それぞれの遺伝子座においてハーディー・ワインバーグ平衡が達成されているか否かを検定することにより調べることができる。
【００２９】
ハプロタイプ推定方法では、ＥＭアルゴリズムを使用して遺伝子型プール情報から、以下のようにしてハプロタイプ頻度を推定することができる。
【００３０】
以下の説明においては、先ず、Ｎ個のＤＮＡプールを作製し、各ＤＮＡプールにはＭ個の異なる個体からのサンプルを含むものとする。このとき、ＤＮＡをプールするためのサンプル選択は、ランダムに行われるものとする。またサンプル選択は、１個体から１度とする。
【００３１】
次いで、連鎖するＬ個の遺伝子座について定量的ＤＮＡタイピングを行い、各ＤＮＡプールに含まれるサンプルについて遺伝子型を決定する。したがって、ＤＮＡプール毎に決定された遺伝子型が遺伝子型プール情報となる。ここで、Ｌ個の遺伝子座は、二対立遺伝子であってもよいし、多対立遺伝子（ｍｕｌｔｉａｌｌｅｌｉｃ）であってもよい。各ＤＮＡプールについて定量的ＤＮＡタイピングを行うことによって、各ＤＮＡプールにおける各遺伝子座についての対立遺伝子コピー数が正確に決定される。
【００３２】
ここで、「対立遺伝子コピー」又は「ハプロタイプコピー」とは、所定の個体又は所定のプールにおける所定の遺伝子座の対立遺伝子又はハプロタイプを意味する。個体が所定の遺伝子座においてホモ接合体である場合、その個体は当該遺伝子座において１個の対立遺伝子（ただし２個の対立遺伝子コピー）を有すると判断される。すなわち、単一のプールにおいては、所定の遺伝子座で２Ｍ個の対立遺伝子コピー（２Ｍ個の対立遺伝子ではない）が存在することとなる。また、Ｍ＝１（単一の個体プール）である場合、単一のプールにおける定量的ＤＮＡタイピングは、各個体に対する一般的なＤＮＡタイピングと等しくなる。
【００３３】
以上のように、Ｎ個のＤＮＡプールについて遺伝子型プール情報を決定した後、ＥＭアルゴリズムを使用してＮ個の遺伝子型プール情報からハプロタイプ推定を行う。図５は、本発明において用いるＥＭアルゴリズムの各ステップを概略的に示すフローチャートである。以下、図５を参照しながら、ＥＭアルゴリズムによってＮ個の遺伝子型プール情報からハプロタイプ推定を行う方法を詳細に説明する。
【００３４】
第 １ ステップ：ハプロタイプ頻度への実数値の割り当て
Ａｉを第ｉの遺伝子座における対立遺伝子数とする。Ｌ個の遺伝子座に対して可能性があるハプロタイプ数は、
【００３５】
【数１】
である。まず本発明者らは、推定の第１ステップとしてハプロタイプ頻度に実数値を割り当てる。ｐｉを集団における第ｉのハプロタイプ頻度とする（ここでは、ｉ ＝ １，２，．．，Ｕにおいて、ｐｉ≧０）。自然と、
【数２】
となる。
【００３６】
第 ２ ステップ：ハプロタイプの組合せ
ＤＮＡプールは、Ｍ個の個体からのサンプルを含む。従って、２Ｍ個のハプロタイプコピーがプール中に存在するはずである。２Ｍ個のハプロタイプコピーがＵ個のハプロタイプ全体から選択される（反復的サンプリングを許す）場合、少なくとも１つのハプロタイプコピーの組合せが、プールにおいてＬ個の遺伝子座全てで認められたプールされた遺伝子型データに一致するはずである。Ｃｊｍをハプロタイプコピーの第ｍの組合せとする。これは、Ｌ個の遺伝子座について第ｊのプールに対して認められた遺伝子型データと一致する（ここでは、ｍ ＝ １，２，．．，Ｑｊ）。Ｑｊは、第ｊのプールに対して認められた遺伝子型データに一致するハプロタイプコピーの組合せ数を表す（ステップＳ１０１）。
【００３７】
ステップ ３ ：尤度計算
Ｈａｒｄｙ−Ｗｅｉｎｂｅｒｇ平衡であると仮定すると、Ｃｊｍの事前確率（ｐｒｉｏｒ ｐｒｏｂａｂｉｌｉｔｙ）は、
【数３】
［式中、ＲｊｍｉはＣｊｍ内の第ｉのハプロタイプコピー数を表し、ＴｊｍはＣｊｍ内の異なるハプロタイプ数を表す］である。任意のｊおよびｍについて、
【数４】
であることに注意されるべきである。ハプロタイプ頻度が与えられた第ｊのプールに対するデータの尤度は、
【００３８】
【数５】
として計算される。
【００３９】
Ｎ個のプール全てに対する全体の尤度は、
【数６】
であるはずである。なぜなら、Ｈａｒｄｙ−Ｗｅｉｎｂｅｒｇ平衡であると仮定すると、異なる個体におけるハプロタイプコピーの組合せの事象は独立しているはずであるからである（ステップＳ１０２）。
【００４０】
ステップ ４ ：期待値
第ｊのプールに対するＣｊｍの事後確率（ｐｏｓｔｅｒｉｏｒ ｐｒｏｂａｂｉｌｉｔｙ）は、以下のようにベイズ定理（Ｂａｙｅｓ’ ｔｈｅｏｒｅｍ）によって計算される。
【数７】
従って、全プールにおける第ｉのハプロタイプのコピー数の期待値は、
【数８】
である。
【００４１】
ステップ ５ ：最大化
最大化は、全てのｉについてｐｉをＥｉ／（２ＭＮ）で置き換えることによって行われた（ステップＳ１０３）。
【００４２】
ステップ ６ ：反復
ステップ２〜５をＬａｌｌが収束するまで繰り返す（ステップＳ１０４）。Ｌｍａｘは、Ｌａｌｌが収束したときの値を表す。
【００４３】
反復の最後のステップ後のｐｉを、ｐｉの最尤推定値
【数９】
とする。
【００４４】
以上のように、ＥＭアルゴリズムを使用することによって遺伝子型プール情報からハプロタイプ頻度を推定することができる。すなわち、本アルゴリズムによれば、Ｍ個のサンプルを含むＮ個のＤＮＡプールから、Ｌ個の遺伝子座に関するハプロタイプ頻度を推定することができる。
【００４５】
また、以下に説明するように、いわゆるブートストラップ法を実行し、推定されたハプロタイプ頻度に関する標準偏差を求めることができる。
【００４６】
最尤推定値が与えられた Ｃ ｊｍ の事後確率の計算
第ｊのプールに対するＣｊｍ （ハプロタイプの集団頻度がｉ ＝ １，２，．．，Ｕに対して、
【数１０】
であると与えられている）の事後確率は、得られた推定値をステップ２〜４に当てはめることによって得られる。従って、等式（４）によって得られたＢｊｍによって、第ｊのプールに対するＣｊｍ （ハプロタイプの集団頻度がｉ ＝ １，２，．．，Ｕに対して、
【００４７】
【数１１】
であると与えられている）の事後確率が得られる。
【００４８】
連鎖不平衡でないという仮定の下での尤度
連鎖不平衡でないという仮定の下でのデータの尤度計算は、以下のように行われた。
ｑｉｋを集団における第ｉの遺伝子座での第ｋの対立遺伝子の頻度とする。Ｖｉｊｋを第ｊのプールにおける第ｉの遺伝子座での第ｋの対立遺伝子のコピー数とする。またＷｉｊを第ｊのプールにおける第ｉの遺伝子座での異なる対立遺伝子数とする。任意のｉおよびｊに対して、
【数１２】
であることに注意されるべきである。連鎖不平衡でないという仮定の下での第ｉの遺伝子座での第ｊのプールデータの尤度は、
【００４９】
【数１３】
である。連鎖不平衡でない下では、異なる遺伝子座での対立遺伝子は独立しているので、全ての遺伝子座におけるデータの尤度は、
【数１４】
であるはずである。また全てのプールにおける全ての遺伝子座でのデータの尤度は、
【数１５】
であるはずである。
【００５０】
ｌｏｄ スコア
ｌｏｄスコアは、以下のようにして計算した。
【数１６】
【００５１】
連鎖不平衡でないという帰無仮説を排除するためのＰ値は、尤度比を等式（８）に組み込む、−ｌｎ（尤度比）とすることで計算された。なおこの統計値はχ２分布に漸近的に従うと仮定する。自由度は、
【数１７】
であるはずである。
【００５２】
プール法によるハプロタイプ頻度の変動および連鎖不平衡の尺度
プールされた遺伝子型データから推定されたハプロタイプ頻度は、サンプルの様々な組合せによって変動する。そのような変動を調べるために、本発明者らは、様々な個体からのＤＮＡサンプルの様々な組合せを作製し、ハプロタイプ頻度を推定した。従って、合計ＭＮ個の個体が存在する場合、Ｍ個の異なる個体からのサンプルが各プール中に存在し、Ｎ個のプールが作製されるはずである。サンプルの異なる組合せは、
【数１８】
個存在する。この数は非常に大きすぎて、全ての事例を調べることができない。従って、本発明者らはモンテカルロ法を用いて、全ての組合せが同じ確率であると仮定した際のＮ個のプールの組合せをサンプリングした。各サンプルからハプロタイプ頻度およびペアワイズ連鎖不平衡の尺度ＤおよびＤ’を、下記に説明するように推定した。ランダムに選択した１０００個の異なるサンプルの推定値から、尺度および標準偏差を計算した。
【００５３】
標準誤差を推定するための非パラメトリックブートストラップ法
非パラメトリックブートストラップ法を用いて、第ｉのハプロタイプ頻度の標準誤差、
【数１９】
を経験的に推定した。
【００５４】
最初のＤＮＡプールは、Ｎ個のプールから構成された。各プールはＭ個の個体からのＤＮＡを含んでいた。二重サンプリングを許すことで、最初のＮ個のプールから新しいセットのＮ個のプールを選び出すことによって、ブートストラップサンプルを構築した。次いで、新しいセットのプールにおけるデータを、ハプロタイプ頻度推定のアルゴリズム（すなわち、ｉ ＝ １，２，．．，Ｕに対するｐｉ）に当てはめた。
【００５５】
【数２０】
を、第ｂのブートストラップサンプルからの第ｉのハプロタイプ頻度推定値とする。ブートストラップサンプリングをＢ回繰り返し、推定値の平均値を以下のようにして計算した。
【００５６】
【数２１】
次いで、ｐｉに対する平均値の経験的標準誤差を、以下のようにして計算した。
【００５７】
【数２２】
通常、ブートストラップサンプリングを１０，０００回（Ｂ ＝ １０，０００）繰り返して、各ｐｉに対する経験的
【数２３】
を計算した。
【００５８】
Ｄ 、 Ｄ’ およびρ ２ の推定
２つの二対立遺伝子座に対する連鎖不平衡の尺度、Ｄ、Ｄ’およびρ２を、以下のように２つの異なる方法で推定した。ｆｉｊを、第１遺伝子座および第２遺伝座のそれぞれにおける第ｉの対立遺伝子および第ｊの対立遺伝子を含むハプロタイプ頻度とする。Ｄ、Ｄ’およびρ２は、以下のようにｉ ＝ １，２およびｊ ＝ １，２に対する推定値、
【数２４】
から計算された。
【数２５】
【００５９】
第１の方法では、利用可能な全ての遺伝子座についてのデータを用いたハプロタイプ頻度の最尤推定値を用いて、
【数２６】
を計算した。しかしながら、第２の方法では、関係する２つの遺伝子座のみの遺伝子型データを用いて、
【数２７】
を推定した。
【００６０】
上述したアルゴリズムは、仮想的なデータではなく、検体から得た真のデータを用いて集団におけるハプロタイプ頻度を推定することができる。ここで、上述したアルゴリズムにおいては、上述する実施例において実証したように、ハプロタイプの相対頻度が０．１より大であり、かつＤＮＡプールに含まれる検体数が４以上である場合、特に、ハプロタイプの頻度の推定はかなり正確であることが示唆された。
【００６１】
また、本アルゴリズムは、コンピュータープログラムにおいて実行することができる。本アルゴリズムを実行できるコンピュータープログラムは、汎用コンピュータが読取り可能な記録媒体に記録することができる。記録媒体としては、可搬型或いは固定型のいずれであっても良く、例えば、ＭＯ、ＣＤ−ＲＯＭ、ＣＤ−ＲＡＭ、ＤＶＤ−ＲＯＭ及びＤＶＤ−ＲＡＭ等の光学式記録媒体、フレキシブルディスク及びハードディスク等の磁気記録媒体及び半導体メモリ等を挙げることができる。
【００６２】
また、本アルゴリズムを実行できるコンピュータープログラムは、ネットワークに接続されたコンピュータの記録装置に格納しておき、ネットワークを介して他のコンピュータに転送することもできる。本アルゴリズムを実行するコンピュータープログラムを提供する提供媒体としては、様々な形式のコンピュータに読み出し可能な媒体として頒布可能であって、特定のタイプの媒体に限定されるものではない。
【００６３】
一方、本アルゴリズムは、コンピュータとして実現することができる。この場合コンピュータは、本アルゴリズムを実行できるコンピュータープログラムを格納したハードディスク装置等の補助記憶装置と、サンプルの遺伝子型に関する情報を入力できる入力装置と、ハプロタイプ頻度を算出する演算装置と、算出したハプロタイプ頻度を表示する表示装置とを備える。
【００６４】
【実施例】
以下、実施例を用いて本発明を更に詳細に説明するが、本発明の技術的範囲はこれらの実施例に限定されるものではない。
【００６５】
〔実施例１〕
ＳＡＡ 遺伝子におけるハプロタイプ頻度の推定
ＳＡＡ遺伝子については、６つのＳＮＰ座位についての１５６人の被験者から得たハプロタイプデータが開示されている（ＭｏｒｉｇｕｃｈｉらＡｒｔｈｒｉｔｉｓ Ｒｈｅｕｍ ４４：１２６６−１２７２ （２００１））。当該論文において、全被験者についてディプロタイプ構成が決定されており、これらのデータが実際のデータを反映しているものと判断された。
【００６６】
そこで当該論文に開示されているデータを適宜用いて、上述したアルゴリズムを用いたハプロタイプ頻度の推定を行った。ハプロタイプ頻度の推定を行うに際して、当該論文に開示されているハプロタイプデータを混合するとともに、相データを除いた。
【００６７】
多重連鎖座位における相未知の遺伝子型データを用いて、上述したアルゴリズムを実装した「ｌｄｐｏｏｌｅｄ」プログラムによりハプロタイプの頻度推定を行った。
【００６８】
表１には、必要なＣＰＵ時間、収束までの反復試行数、ロッドスコア、χ２値、及びＰ値を示す。それぞれのＰ値は、６つの座位の全てが互いに独立であるという帰無仮説を棄却することの危険率（ｒｉｓｋ）を表す。
【００６９】
【表１】
【００７０】
なお、表１において、Ｌｏｄ値及びχ２値の計算方法は、上述した方法に準じた。自由度は、上記式（９）により算出した場合、５７となるはずである。各プールに含まれる被験者の数が増加するにつれて、ロッドスコア及びχ２値は低下したが、Ｐ値は増大した（表１）。これらの結果は、おそらくは、プーリングのために情報内容が減少することによって生じるものである。
【００７１】
表２には、１人、２人及び４人の被験者から得た遺伝子型データを含む各プールについて、異なる頻度推定プロトコールを用いて推定したハプロタイプ頻度の結果を示す。各プールが２人以上の被験者から得たＤＮＡを含む場合には、各被験者についてのオリジナルデータからサンプリングされる被験者の様々な組み合わせを利用して推定した相対頻度の、平均及び標準偏差も示した。
【００７２】
【表２】
【００７３】
なお、表２における「ハプロタイプ」には、少なくとも１種の推定方法において相対頻度が＞０．００１となるハプロタイプが含まれる。表２における「相対頻度」は１回の推定を行って得られた推定頻度である。表２における「平均」及び「標準偏差」は種々の無作為抽出法を用いてプールを作製して得られた推定頻度の平均及び標準誤差である。
【００７４】
被験者を一人しか含まれないプール（１個体／プール）についての頻度推定の結果は、ＤＮＡがプールされなかった既報論文中で示された結果（ＫｉｔａｍｕｒａらＡｎｎ． Ｈｕｍ． Ｇｅｎｅｔ． ６６：１８３−１９３（２００２））と同じであった。表２から判るように、プール中の被験者数の変化とともに推定したハプロタイプ頻度（相対頻度）が変化したが、プールに含まれる被験者の数が４人までの範囲では、主たるハプロタイプであるＡＣＴＧＣＣ、ＡＣＣＧＴＣ、ＡＧＣＧＣＴの相対頻度は依然として適正であった。
【００７５】
また、表２から判るように、上記の主たる３つのハプロタイプにおいては、プールを作成するために様々なランダムサンプリングを行うことによって推定した相対頻度に対する標準偏差が平均の１０％未満であった。しかしながら、マイナーなハプロタイプ（相対頻度が０．１未満）については、頻度推定は正確ではなかった。例えば、ハプロタイプＡＣＴＧＴＣの頻度は、４人の被験者を含むプールについての頻度推定によっては０．０と推定されたが、被験者を１人含むプールについての頻度推定によっては０．０１３であった。後者の頻度推定が正確であれば、このハプロタイプは、１５６人（すなわちハプロタイプコピーが３１２）の間で４回しか出現しないこととなる。
【００７６】
また、上述した方法に準じて、表２に示した相対頻度から算出したＤ値及びＤ’値を表３に示す。
【００７７】
【表３】
【００７８】
なお、表３においてＤ値及びＤ’値は１回の推定を行って得られた値であり、括弧内の値は無作為抽出法によって算出される「平均±標準偏差」を示す。この場合、ｆｉｊ値、すなわち２つの遺伝子座についてのハプロタイプの推定頻度は、上述した方法に記載したように、６つの遺伝子座についてのハプロタイプの推定頻度から計算した。ｆｉｊを２つの遺伝子座についての遺伝子型データから算出した場合、Ｄ値及びＤ’値が非常に類似している場合もあったが、その２つの方法が全く異なる値をもたらした場合もあった（データは示していない）。また、表３において、各プールが２人以上の被験者を含む場合は、異なる組み合わせの被験者を用いて算出した値の平均及び標準偏差を示した。この値は、異なる頻度推定プロトコールの間で異なっていたが、プールに含まれる被験者の数が最大４人までであって、かつ｜Ｄ｜値が０．１を上回る限りは、依然としてかなり一致した値を示した。しかしながら、いくつかの場合においては、プールを作成するために様々なランダムサンプリングを行うことによって算出したＤ及びＤ’の標準偏差が、平均の約５０％であり、異なる頻度推定プロトコールにより算出される値は互いに大きく異なっていた。｜Ｄ｜値は０．１を上回るが｜Ｄ’｜値は０．１を上回らない場合は、Ｄ’値がそれほど異なっていなかった点は興味深い。
【００７９】
一方、各ＤＮＡプールに含まれるハプロタイプコピーの組み合わせについて頻度推定を行った結果を表４に示す。各プロトコールについて、最初の１２人の被験者に対応するデータ部分だけを示した。
【００８０】
【表４】
【００８１】
なお、表４において「確率」は、プールに対するハプロタイプの各組合せの事後確率を意味する。この事後確率（表４における「確率」）によって推定ハプロタイプ頻度が示される。
【００８２】
表４から判るように、大部分のプールにおいて、最も確率の高いハプロタイプコピーの組み合わせについての事後確率が１又はほぼ１であったことを示している。ハプロタイプコピーの頻度推定を行った組み合わせについてその内容を比較すると、異なる頻度推定プロトコールの間で、多くの場合、その内容は互いに一致していた。例えば、表４から判るように、２人の被験者を含むプールを用いるプロトコールの場合のプール番号１の内容は、被験者を１人含むプールを用いた場合のプール番号１及びプール番号２の内容を組み合わせたものと同となっていた。なお、あるプロトコールによって頻度推定を行ったプールの内容は、別のプロトコールによって頻度推定を行ったプールの内容と一致しない場合もあった。
【００８３】
一方、上述したように推定したハプロタイプ頻度には誤差があるため、上述した方法に準じてブートストラップ法による標準誤差を算出した。上述したハプロタイプ頻度の推定の際と同様に、１５６人の被験者から得たＳＡＡ遺伝子についての遺伝子型データを用いて、本発明者らは上述の通りに１人、２人、又は４人の被験者を含むプールを作成した。その後、上述した方法に準じてブートストラップ法に適用した。結果を図１に示す。
【００８４】
なお、図１には、１人、２人又は４人の被験者を含むプールを用いて推定したハプロタイプ頻度の平均及び標準誤差を示す。ブートストラップサンプリングは、各ＤＮＡプールを用いたハプロタイプ頻度の推定について、それぞれ１０，０００回繰り返した。図１から判るように、各ＤＮＡプールに含まれる被験者の数に拘わらず、ハプロタイプＡＣＴＧＣＣ、ＡＣＣＧＴＣ及びＡＧＣＧＣＴについての推定頻度は、かなり安定している。また、図１に示すように、これらのハプロタイプに関する標準誤差を示すバーの長さは、平均を示すバーと比較してかなり短い。さらに、図１から判るように、各ＤＮＡプールについて推定したハプロタイプ頻度の平均は、同一のハプロタイプについてはほぼ同じ値を示した。しかしながら、マイナーなハプロタイプ（ｐｉ＜ ０．１）については、各ＤＮＡプールについて推定したハプロタイプ頻度の平均が有意に異なっており、前記誤差を示すバーは平均と比較するとかなり長かった。
【００８５】
なお、表１に示したように、各計算を行うのに必要な時間及びメモリーを記録した。１ＧＨｚのＰｅｎｔｉｕｍ ＩＩＩ（ＣＰＵ）及び１．５ＧＢメモリーを搭載したコンピュータを用いる場合、プールに含める被験者の数は、遺伝子座の数が６であれば最大で６人であった。遺伝子座の数が１３であれば、プールに含む被験者は最大でたった２人となる可能性がある。
【００８６】
一方で、複数のＤＮＡサンプルをプールしない場合（換言すれば１個体の被験者を含むＤＮＡプールの場合）には、２５個の遺伝子座について計算することが可能であった。これは、上述したアルゴリズムを実装したコンピュータープログラムが、各ＤＮＡプールに含まれるハプロタイプコピーの可能な組み合わせを用いるものであって、このステップが大量のメモリーを消費するためである。そのような組み合わせの数は、所定の遺伝子座に存在する対立遺伝子の数の検定力関数に従って増加し、かつ、ＤＮＡプールに含まれる被験者数の階乗として増加する。従って、コンピュータを用いた計算に要する時間及びメモリーがこれらの要因に依存するのは、上記のような組み合わせのための余裕が必要となるためである。
【００８７】
〔実施例２〕
実施例２では、メチレンテトラヒドロ葉酸レダクターゼ（ＭＴＨＦＲ）遺伝子における連鎖する２つの遺伝子座についてハプロタイプ頻度を推定した。ＭＴＨＦＲ遺伝子は、葉酸代謝に関係するＭＴＨＦＲ酵素をコードしている。ＭＴＨＦＲ遺伝子におけるこれら遺伝子座に関する情報は、計８０人の被験者に関する発表済データ（ＵｒａｎｏらＰｈａｒｍａｃｏｇｅｎｅｔｉｃｓ １２：１８３−１９０ （２００２））を用いた。これらの発表済データを用いてＤＮＡプールを作成し、プールされたデータからパラメーターを算出した。結果を表５に示す。
【００８８】
【表５】
【００８９】
なお、表５における「相対頻度」は１回の推定を行って得られた推定頻度である。表５における「平均」及び「標準偏差」は種々の無作為抽出法を用いてプールを作製して得られた推定頻度の平均及び標準誤差である。
【００９０】
表５に示したように、ＤＮＡプールに含まれる被験者の人数に拘わらず、推定されたハプロタイプ頻度が安定していた。言い換えると、推定されたハプロタイプ頻度は、４人の被験者を含むＤＮＡプールを使用した場合でさえかなり正確だったことを示唆している。これは、これらの遺伝子座に関する連鎖不平衡が非常に強いためであると考えられる。また、表５からは、相対頻度が０．１を上回るハプロタイプ同士については、頻度推定がかなり正確であったことも示している。
【００９１】
図２は、ＭＴＨＦＲ遺伝子のハプロタイプデータに対するブートストラップ法の結果を示す。これらのデータは、本アルゴリズムによれば、ハプロタイプの頻度がかなり高い場合にはＭＴＨＦＲ遺伝子のハプロタイプ頻度をかなり正確に推定できることを示唆している。
【００９２】
ＭＴＨＦＲ遺伝子について推定したハプロタイプ頻度から、実施例１と同様に計算したＤ値及びＤ’値を表６に示す。また、２人及び４人の被験者を含むプールを用いるプロトコールによって得たデータに関して、様々なランダムサンプリングを行ったものに対する平均及び標準偏差も示す。
【００９３】
【表６】
これらのデータは、推定値の変動が大きいことも示している。これは、この場合の｜Ｄ｜値が全て０．１未満であったためであると考えられる。
【００９４】
〔実施例３〕
実施例３では、Ｎ−アセチルトランスフェラーゼ２（ＮＡＴ２）遺伝子における連鎖する７つの遺伝子座についてハプロタイプ頻度を推定した。ＮＡＴ２遺伝子は、Ｎ−アセチル残基の転移に関与するＮＡＴ２酵素をコードしている。ＮＡＴ２遺伝子におけるこれら遺伝子座に関する情報は、計１１６人の被験者に関する発表済データ（ＴａｎａｋａらＪ． Ｒｈｅｕｍａｔｏｌ （Ｉｎ ｐｒｅｓｓ））を用いた。これらの発表済データを用いてＤＮＡプールを作成し、プールされたデータからパラメーターを算出した。結果を表７に示す。
【００９５】
【表７】
【００９６】
なお、表７における「ハプロタイプ」には、少なくとも１種の推定方法において相対頻度が＞０．００１となるハプロタイプが含まれる。表７における「相対頻度」は１回の推定を行って得られた推定頻度である。表７における「平均」及び「標準偏差」は種々の無作為抽出法を用いてプールを作製して得られた推定頻度の平均及び標準誤差である。
【００９７】
表７に示したように、ＤＮＡプールに含まれる被験者の人数に拘わらず、推定されたハプロタイプ頻度が安定していた。言い換えると、推定されたハプロタイプ頻度は、４人の被験者を含むＤＮＡプールを使用した場合でさえかなり正確だったことを示唆している。
【００９８】
図３は、ＮＡＴ２遺伝子のハプロタイプデータに対するブートストラップ法の結果を示す。これらのデータは、本アルゴリズムによれば、ハプロタイプの頻度がかなり高い場合にはＭＴＨＦＲ遺伝子のハプロタイプ頻度をかなり正確に推定できることを示唆している。
【００９９】
ＮＡＴ２遺伝子について推定したハプロタイプ頻度から、実施例１と同様に計算したＤ値及びＤ’値を表８に示す。
【０１００】
【表８】
これらのデータは、｜Ｄ｜値が０．１を上回る場合には推定値の変動が小さいことも示している。
【０１０１】
〔実施例４〕
実施例４では、スムーセリン（ｓｍｏｏｔｈｅｌｉｎ）遺伝子における連鎖する３６個の遺伝子座についてハプロタイプ頻度を推定した。ｓｍｏｏｔｈｅｌｉｎ遺伝子は、第２２染色体に存在し、平滑筋細胞骨格を構成するｓｍｏｏｔｈｅｌｉｎをコードする遺伝子である。ｓｍｏｏｔｈｅｌｉｎ遺伝子におけるこれら遺伝子座に関する情報は、計１０２人の日本人の被験者に関するデータを用いた。
【０１０２】
本例では、遺伝子座の数が３６箇所と多かったため、全遺伝子座を３つの領域に任意に分けてハプロタイプ頻度の推定を行った。すなわち、ｓｍｏｏｔｈｅｌｉｎ遺伝子における５’末端から３’末端に向かって、１３個の遺伝子座を含む第１の領域、１３個の遺伝子座を含む第２の領域、及び１０個の遺伝子座を含む第３の領域となるように全遺伝子座を分けた。
【０１０３】
これら各領域について、前記データを用いてＤＮＡプールを作成し、プールされたデータからハプロタイプの相対頻度を算出した。結果を表９に示す。
【０１０４】
【表９】
【０１０５】
なお、表９において「ハプロタイプ」は、１個体を含むＤＮＡプールを用いた相対頻度が≧０．００１となるハプロタイプが含まれる。また、表９の「ハプロタイプ」欄における「Ｄ」は欠失を意味し、「Ｉ」は挿入を意味する。
【０１０６】
表９に示すように、１個体を含むＤＮＡプール及び２個体を含むＤＮＡプールのいずれを用いた場合でも、推定したハプロタイプの頻度（表９における相対頻度）は、そのハプロタイプ頻度が０．１を上回る限りにおいては、非常に近い値を示した。しかしながらマイナーなハプロタイプ（相対頻度が０．１未満）については、１個体を含むＤＮＡプールを用いた場合と２個体を含むＤＮＡプールを用いた場合とで類似しているとは言えなかった。
【０１０７】
次に、本例では、連鎖不平衡の強さを評価するため、上記方法に準じてペアワイズ方式で連鎖不平衡測定値ρ２を算出した。３６個の遺伝子座に関し（３６×３５）／２＝６３０ペアについて、算出したρ２を図４に示す。なお、本例では、ｆｉｊ値は２つの遺伝子座のみの遺伝子型データを用いて算出した。図４から判るように、多少の差異は観察されたものの、１個体を含むＤＮＡプールを用いて算出されたρ２値と２個体を含むＤＮＡプールを用いて算出されたρ２値とは、互いに非常に類似していた。様々なペアの遺伝子座について算出したρ２値の平均及び標準偏差は、１個体を含むＤＮＡプールを用いた場合はそれぞれ０．１１４及び０．２２４であり、２個体を含むＤＮＡプール用いた場合はそれぞれ０．１１８及び０．２２７であった。１個体を含むＤＮＡプールを用いて算出されたρ２値と２個体を含むＤＮＡプールを用いて算出されたρ２値との差の絶対値の平均及び標準偏差は、それぞれ０．０１１及び０．０２２であった。また、１個体を含むＤＮＡプールを用いて算出されたρ２値と２個体を含むＤＮＡプールを用いて算出されたρ２値との差の絶対値は、平均の約９．５％であった。
【０１０８】
【発明の効果】
以上、詳細に説明したように、本発明に係るハプロタイプ頻度推定方法によれば、複数のプールした塩基配列データを入力値として集団のハプロタイプ頻度を推定することができる。換言すると、本発明に係るハプロタイプ頻度推定方法は、プールすることなしに複数の塩基配列データをそれぞれ入力値として用いた場合と同様に、ハプロタイプ頻度を推定することができる。
【図面の簡単な説明】
【図１】ＳＡＡ遺伝子遺伝子についての遺伝子型データに関して１人、２人又は４人の被験者を含むプールを用いて、ブートストラップ法を適用して算出した、推定ハプロタイプ頻度の平均及び標準誤差を示す特性図である。
【図２】ＭＴＨＦＲ遺伝子についての遺伝子型データに関して１人、２人又は４人の被験者を含むプールを用いて、ブートストラップ法を適用して算出した、推定ハプロタイプ頻度の平均及び標準誤差を示す特性図である。
【図３】ＮＡＴ２遺伝子遺伝子についての遺伝子型データに関して１人、２人又は４人の被験者を含むプールを用いて、ブートストラップ法を適用して算出した、推定ハプロタイプ頻度の平均及び標準誤差を示す特性図である。
【図４】ｓｍｏｏｔｈｅｌｉｎ遺伝子についての遺伝子型データに関して、ペアワイズ方式で連鎖不平衡測定値ρ２を算出した結果を示す特性図である。
【図５】ＥＭアルゴリズムの各ステップを概略的に示すフローチャートである。
【発明の属する技術分野】
本発明は、最大尤度推定アルゴリズムとして知られている期待値最大化アルゴリズム（以下、ＥＭアルゴリズムと呼ぶ）を用いたハプロタイプ頻度推定方法、ハプロタイプ頻度推定アルゴリズム及びハプロタイプ頻度推定プログラムに関する。
【０００２】
【従来の技術】
遺伝統計学的手法とは、不明の疾患遺伝子を見出す方法の一つであり、個体間の遺伝子多型の違いと、各個体の形質データとから統計的な計算によって、形質に関連する遺伝子座を探す方法である。すなわち、遺伝統計学的手法によれば、遺伝的多型と形質（表現型、症状等）の情報のみから統計学を用いて原因（遺伝子座）を解明することができる。
【０００３】
遺伝統計学的手法としては、主としてパラメトリック連鎖解析とノンパラメトリック連鎖解析とに分類される。「連鎖解析」とは、メンデルの第３法則（独立の法則）の例外現象である連鎖を利用して表現型に関連する遺伝子座の染色体上に、疾患原因遺伝子座の存在領域を絞り込んでいく（位置をマッピングする）遺伝統計学的手法である。「パラメトリック」とは、解析に先立ち、遺伝子型に応じて固有の発症確率（＝浸透率）を仮定する（遺伝形式を仮定する）ことを意味する。遺伝子型における浸透率を仮定するということは、すなわち疾患原因遺伝子座位での遺伝子型によって発症確率が記述できるということである。浸透率が仮定できるときには、表現型（罹患の有無）からその疾患の原因遺伝子座での遺伝子型を簡単に関連づけることができる。「ノンパラメトリック」とは、遺伝形式を仮定せずに解析を行うことを意味する。
【０００４】
ところで、ハプロタイプとは、複数の遺伝子座における各対立遺伝子の組合せとして定義される。表現型を遺伝子領域にマッピングする研究など、ハプロタイプ推定は、多くの遺伝子研究において重要な意味を持っている。通常、集団からサンプリングされた複数個人の遺伝子型データに基づいて集団のハプロタイプ推定及び個人のディプロタイプ推定が行われている。
【０００５】
【発明が解決しようとする課題】
ところで、複数個人の遺伝子型データを用いた場合、例えば、最大尤度推定を行うためのＥＭアルゴリズムを使用して集団のハプロタイプ推定及び個人のディプロタイプ推定を行うことができる。しかしながら、この場合、複数個人の遺伝子型データを用いているため、入力するためのデータ量が膨大であり、現状使用されている高スペックのコンピュータでさえも処理することが困難であった。
【０００６】
そこで、本発明は、上述したような実状に鑑み、処理速度等のコンピュータにおけるスペックに依存せずに情報処理が可能なＥＭアルゴリズムを用いたハプロタイプ推定方法、ハプロタイプ頻度推定アルゴリズム及びハプロタイプ頻度推定プログラムを提供することを目的とする。
【０００７】
【課題を解決するための手段】
上述した目的を達成するため、本発明者が鋭意検討した結果、複数個人の遺伝子型データの代わりに、複数個人の遺伝子型データをプールしたデータを用いることにより、集団のハプロタイプ頻度推定に要する処理能力を低く抑えられ、また、複数個人の遺伝子型データをプールしたデータを用いても集団のハプロタイプ頻度推定を正確に行えることを見出し、本発明を完成するに至った。
【０００８】
すなわち、本発明は、以下を包含する。
（１）集団に含まれる複数の検体に関する遺伝子型情報を集積してある遺伝子型プール情報を入力値として、期待値最大化アルゴリズムにより上記集団におけるハプロタイプ頻度を推定する方法。
【０００９】
（２）上記複数の検体に由来するＤＮＡを含む複数のＤＮＡプールを調製し、当該複数のＤＮＡプールに含まれるＤＮＡにおける所定の遺伝子座に関するタイピングを行う工程を含むことを特徴とする（１）記載のハプロタイプ頻度を推定する方法。
【００１０】
（３）２以上６以下の検体に関する遺伝子型情報を集積することして遺伝子型プール情報を作成することを特徴とする（１）記載のハプロタイプ頻度を推定する方法。
【００１１】
（４）前記期待値最大化アルゴリズムはＥＭアルゴリズムであり、前記遺伝子型プール情報に最尤推定を行うことによりハプロタイプ頻度を推定することを特徴とする（１）記載のハプロタイプ頻度を推定する方法。
【００１２】
（５）集団に含まれる複数の検体に関する遺伝子型情報を集積してある遺伝子型プール情報を入力値として、期待値最大化アルゴリズムにより上記集団におけるハプロタイプ頻度を推定するアルゴリズム。
【００１３】
（６）上記複数の検体に由来するＤＮＡを含む複数のＤＮＡプールを調製し、当該複数のＤＮＡプールに含まれるＤＮＡにおける所定の遺伝子座に関するタイピングを行う工程を含むことを特徴とする（５）記載のハプロタイプ頻度を推定するアルゴリズム。
【００１４】
（７）２以上６以下の検体に関する遺伝子型情報を集積することして遺伝子型プール情報を作成することを特徴とする（５）記載のハプロタイプ頻度を推定するアルゴリズム。
【００１５】
（８）前記期待値最大化アルゴリズムはＥＭアルゴリズムであり、前記遺伝子型プール情報に最尤推定を行うことによりハプロタイプ頻度を推定することを特徴とする（５）記載のハプロタイプ頻度を推定するアルゴリズム。
【００１６】
（９）集団に含まれる複数の検体に関する遺伝子型情報を集積してある遺伝子型プール情報を入力値として、期待値最大化アルゴリズムにより上記集団におけるハプロタイプ頻度を推定するプログラム。
【００１７】
（１０）上記複数の検体に由来するＤＮＡを含む複数のＤＮＡプールを調製し、当該複数のＤＮＡプールに含まれるＤＮＡにおける所定の遺伝子座に関するタイピングを行う工程を含むことを特徴とする（９）記載のハプロタイプ頻度を推定するプログラム。
【００１８】
（１１）２以上６以下の検体に関する遺伝子型情報を集積することして遺伝子型プール情報を作成することを特徴とする（９）記載のハプロタイプ頻度を推定するプログラム。
【００１９】
（１２）前記期待値最大化アルゴリズムはＥＭアルゴリズムであり、前記遺伝子型プール情報に最尤推定を行うことによりハプロタイプ頻度を推定することを特徴とする（９）記載のハプロタイプ頻度を推定するプログラム。
【００２０】
【発明の実施の形態】
以下、本発明に係るＥＭアルゴリズムを用いたハプロタイプ頻度推定方法（以下、単に「ハプロタイプ頻度推定方法」と呼ぶ）を詳細に説明する。
【００２１】
ハプロタイプ頻度推定方法は、集団における複数の遺伝子型情報を集積（プール）した遺伝子型プール情報から、該集団におけるハプロタイプ頻度を算出する工程を含んでいる。
【００２２】
「遺伝子型」とは、所定の遺伝子座における一対の対立遺伝子の組み合わせ（常染色体における遺伝子座の場合）である。一対の対立遺伝子は、父性由来の一方の対立遺伝子と母性由来の他方の対立遺伝子とからなる。遺伝子型は、例えば所定の遺伝子座に多型部位が存在する場合、該多型部位に相当する一対の対立遺伝子についての遺伝子の個人についての特徴を示す。
【００２３】
「遺伝子型プール情報」とは、個人に関する遺伝子型情報を、複数の個人について集積した情報を意味する。複数の個人とは、特に限定されないが、２〜６人とすることが好ましく、２〜４人とすることがより好ましい。また、遺伝子型プール情報は、複数の遺伝子型情報を、複数の個人について集積した情報であってもよい。
【００２４】
遺伝子型プール情報となる遺伝子の特徴としては、例えば１塩基多型、マイクロサテライト及び挿入／欠失等の遺伝子変異を挙げることができる。遺伝子型に関する情報には、前記遺伝子変異のうち、１つの種類に限定されず、複数の種類の遺伝子変異が含まれていてもよい。
【００２５】
遺伝子型情報に複数の遺伝子座が含まれる場合、これら遺伝子座は連鎖していることが好ましい。従って、本発明においては、予め連鎖していることが知られている遺伝子座について遺伝子型プール情報を取得することが望ましいが、連鎖の有無が不明な場合もあり、必ずしもこの限りではない。
【００２６】
また、「ハプロタイプ」とは、複数の遺伝子座における各対立遺伝子の組合せを意味する。ハプロタイプは、複数の１塩基置換（ＳＮＰ； ｓｉｎｇｌｅ ｎｕｃｌｅｏｔｉｄｅ ｐｏｌｙｍｏｒｐｈｉｓｍ ）、マイクロサテライト、挿入／欠失等の多型部位が近接して存在し、これらが連鎖する場合に形成されるものである。また、常染色体上の遺伝子部位について、個人は、それぞれ父性或いは母性に由来する一対のハプロタイプを保有する。このような個人が保有する一対のハプロタイプの組み合わせを「ディプロタイプ」という。
【００２７】
「ハプロタイプ頻度」とは、複数の遺伝子座における各対立遺伝子について所定の組合せ（ハプロタイプ）が出現する頻度である。ハプロタイプは、複数の遺伝子座が狭い遺伝子領域に含まれる場合には世代交代によっても変化することは少ない。これは、複数の多型部位が存在する遺伝子領域内に交差が起きることによりハプロタイプは変化するが、当該遺伝子領域が狭い場合、交差が起きる確率は極めて低い（１世代あたり、１００，０００ｋｂに１回に過ぎない）ためである。そのため特定の集団内、例えば日本人では、それぞれのハプロタイプ頻度が決まっている場合が多い。ハプロタイプ頻度は、集団ごとに異なる場合が多い。
【００２８】
「集団」とは、遺伝子型情報を提供する個人の集まりを意味する。集団の人数は、算出されるハプロタイプ頻度の信頼性を高めるために、なるべく多くの人数から構成されることが望ましい。集団を構成する人数としては、２５名以上であることが好ましいが、特に制限されない。また、集団は均一（遺伝的背景が近接）であることが好ましい。例えば、ハプロタイプ頻度の推定を行う個人が日本人であるならば、遺伝子型情報を取得すべき「集団」は、日本人から構成される集団であることが好ましい。集団が均一か否かについては、それぞれの遺伝子座においてハーディー・ワインバーグ平衡が達成されているか否かを検定することにより調べることができる。
【００２９】
ハプロタイプ推定方法では、ＥＭアルゴリズムを使用して遺伝子型プール情報から、以下のようにしてハプロタイプ頻度を推定することができる。
【００３０】
以下の説明においては、先ず、Ｎ個のＤＮＡプールを作製し、各ＤＮＡプールにはＭ個の異なる個体からのサンプルを含むものとする。このとき、ＤＮＡをプールするためのサンプル選択は、ランダムに行われるものとする。またサンプル選択は、１個体から１度とする。
【００３１】
次いで、連鎖するＬ個の遺伝子座について定量的ＤＮＡタイピングを行い、各ＤＮＡプールに含まれるサンプルについて遺伝子型を決定する。したがって、ＤＮＡプール毎に決定された遺伝子型が遺伝子型プール情報となる。ここで、Ｌ個の遺伝子座は、二対立遺伝子であってもよいし、多対立遺伝子（ｍｕｌｔｉａｌｌｅｌｉｃ）であってもよい。各ＤＮＡプールについて定量的ＤＮＡタイピングを行うことによって、各ＤＮＡプールにおける各遺伝子座についての対立遺伝子コピー数が正確に決定される。
【００３２】
ここで、「対立遺伝子コピー」又は「ハプロタイプコピー」とは、所定の個体又は所定のプールにおける所定の遺伝子座の対立遺伝子又はハプロタイプを意味する。個体が所定の遺伝子座においてホモ接合体である場合、その個体は当該遺伝子座において１個の対立遺伝子（ただし２個の対立遺伝子コピー）を有すると判断される。すなわち、単一のプールにおいては、所定の遺伝子座で２Ｍ個の対立遺伝子コピー（２Ｍ個の対立遺伝子ではない）が存在することとなる。また、Ｍ＝１（単一の個体プール）である場合、単一のプールにおける定量的ＤＮＡタイピングは、各個体に対する一般的なＤＮＡタイピングと等しくなる。
【００３３】
以上のように、Ｎ個のＤＮＡプールについて遺伝子型プール情報を決定した後、ＥＭアルゴリズムを使用してＮ個の遺伝子型プール情報からハプロタイプ推定を行う。図５は、本発明において用いるＥＭアルゴリズムの各ステップを概略的に示すフローチャートである。以下、図５を参照しながら、ＥＭアルゴリズムによってＮ個の遺伝子型プール情報からハプロタイプ推定を行う方法を詳細に説明する。
【００３４】
第 １ ステップ：ハプロタイプ頻度への実数値の割り当て
Ａｉを第ｉの遺伝子座における対立遺伝子数とする。Ｌ個の遺伝子座に対して可能性があるハプロタイプ数は、
【００３５】
【数１】
【数２】
【００３６】
第 ２ ステップ：ハプロタイプの組合せ
ＤＮＡプールは、Ｍ個の個体からのサンプルを含む。従って、２Ｍ個のハプロタイプコピーがプール中に存在するはずである。２Ｍ個のハプロタイプコピーがＵ個のハプロタイプ全体から選択される（反復的サンプリングを許す）場合、少なくとも１つのハプロタイプコピーの組合せが、プールにおいてＬ個の遺伝子座全てで認められたプールされた遺伝子型データに一致するはずである。Ｃｊｍをハプロタイプコピーの第ｍの組合せとする。これは、Ｌ個の遺伝子座について第ｊのプールに対して認められた遺伝子型データと一致する（ここでは、ｍ ＝ １，２，．．，Ｑｊ）。Ｑｊは、第ｊのプールに対して認められた遺伝子型データに一致するハプロタイプコピーの組合せ数を表す（ステップＳ１０１）。
【００３７】
ステップ ３ ：尤度計算
Ｈａｒｄｙ−Ｗｅｉｎｂｅｒｇ平衡であると仮定すると、Ｃｊｍの事前確率（ｐｒｉｏｒ ｐｒｏｂａｂｉｌｉｔｙ）は、
【数３】
【数４】
【００３８】
【数５】
【００３９】
Ｎ個のプール全てに対する全体の尤度は、
【数６】
【００４０】
ステップ ４ ：期待値
第ｊのプールに対するＣｊｍの事後確率（ｐｏｓｔｅｒｉｏｒ ｐｒｏｂａｂｉｌｉｔｙ）は、以下のようにベイズ定理（Ｂａｙｅｓ’ ｔｈｅｏｒｅｍ）によって計算される。
【数７】
【数８】
【００４１】
ステップ ５ ：最大化
最大化は、全てのｉについてｐｉをＥｉ／（２ＭＮ）で置き換えることによって行われた（ステップＳ１０３）。
【００４２】
ステップ ６ ：反復
ステップ２〜５をＬａｌｌが収束するまで繰り返す（ステップＳ１０４）。Ｌｍａｘは、Ｌａｌｌが収束したときの値を表す。
【００４３】
反復の最後のステップ後のｐｉを、ｐｉの最尤推定値
【数９】
【００４４】
以上のように、ＥＭアルゴリズムを使用することによって遺伝子型プール情報からハプロタイプ頻度を推定することができる。すなわち、本アルゴリズムによれば、Ｍ個のサンプルを含むＮ個のＤＮＡプールから、Ｌ個の遺伝子座に関するハプロタイプ頻度を推定することができる。
【００４５】
また、以下に説明するように、いわゆるブートストラップ法を実行し、推定されたハプロタイプ頻度に関する標準偏差を求めることができる。
【００４６】
最尤推定値が与えられた Ｃ ｊｍ の事後確率の計算
第ｊのプールに対するＣｊｍ （ハプロタイプの集団頻度がｉ ＝ １，２，．．，Ｕに対して、
【数１０】
【００４７】
【数１１】
【００４８】
連鎖不平衡でないという仮定の下での尤度
連鎖不平衡でないという仮定の下でのデータの尤度計算は、以下のように行われた。
ｑｉｋを集団における第ｉの遺伝子座での第ｋの対立遺伝子の頻度とする。Ｖｉｊｋを第ｊのプールにおける第ｉの遺伝子座での第ｋの対立遺伝子のコピー数とする。またＷｉｊを第ｊのプールにおける第ｉの遺伝子座での異なる対立遺伝子数とする。任意のｉおよびｊに対して、
【数１２】
【００４９】
【数１３】
【数１４】
【数１５】
【００５０】
ｌｏｄ スコア
ｌｏｄスコアは、以下のようにして計算した。
【数１６】
連鎖不平衡でないという帰無仮説を排除するためのＰ値は、尤度比を等式（８）に組み込む、−ｌｎ（尤度比）とすることで計算された。なおこの統計値はχ２分布に漸近的に従うと仮定する。自由度は、
【数１７】
【００５２】
プール法によるハプロタイプ頻度の変動および連鎖不平衡の尺度
プールされた遺伝子型データから推定されたハプロタイプ頻度は、サンプルの様々な組合せによって変動する。そのような変動を調べるために、本発明者らは、様々な個体からのＤＮＡサンプルの様々な組合せを作製し、ハプロタイプ頻度を推定した。従って、合計ＭＮ個の個体が存在する場合、Ｍ個の異なる個体からのサンプルが各プール中に存在し、Ｎ個のプールが作製されるはずである。サンプルの異なる組合せは、
【数１８】
【００５３】
標準誤差を推定するための非パラメトリックブートストラップ法
非パラメトリックブートストラップ法を用いて、第ｉのハプロタイプ頻度の標準誤差、
【数１９】
【００５４】
最初のＤＮＡプールは、Ｎ個のプールから構成された。各プールはＭ個の個体からのＤＮＡを含んでいた。二重サンプリングを許すことで、最初のＮ個のプールから新しいセットのＮ個のプールを選び出すことによって、ブートストラップサンプルを構築した。次いで、新しいセットのプールにおけるデータを、ハプロタイプ頻度推定のアルゴリズム（すなわち、ｉ ＝ １，２，．．，Ｕに対するｐｉ）に当てはめた。
【００５５】
【数２０】
【００５６】
【数２１】
【００５７】
【数２２】
【数２３】
【００５８】
Ｄ 、 Ｄ’ およびρ ２ の推定
２つの二対立遺伝子座に対する連鎖不平衡の尺度、Ｄ、Ｄ’およびρ２を、以下のように２つの異なる方法で推定した。ｆｉｊを、第１遺伝子座および第２遺伝座のそれぞれにおける第ｉの対立遺伝子および第ｊの対立遺伝子を含むハプロタイプ頻度とする。Ｄ、Ｄ’およびρ２は、以下のようにｉ ＝ １，２およびｊ ＝ １，２に対する推定値、
【数２４】
【数２５】
第１の方法では、利用可能な全ての遺伝子座についてのデータを用いたハプロタイプ頻度の最尤推定値を用いて、
【数２６】
【数２７】
【００６０】
上述したアルゴリズムは、仮想的なデータではなく、検体から得た真のデータを用いて集団におけるハプロタイプ頻度を推定することができる。ここで、上述したアルゴリズムにおいては、上述する実施例において実証したように、ハプロタイプの相対頻度が０．１より大であり、かつＤＮＡプールに含まれる検体数が４以上である場合、特に、ハプロタイプの頻度の推定はかなり正確であることが示唆された。
【００６１】
また、本アルゴリズムは、コンピュータープログラムにおいて実行することができる。本アルゴリズムを実行できるコンピュータープログラムは、汎用コンピュータが読取り可能な記録媒体に記録することができる。記録媒体としては、可搬型或いは固定型のいずれであっても良く、例えば、ＭＯ、ＣＤ−ＲＯＭ、ＣＤ−ＲＡＭ、ＤＶＤ−ＲＯＭ及びＤＶＤ−ＲＡＭ等の光学式記録媒体、フレキシブルディスク及びハードディスク等の磁気記録媒体及び半導体メモリ等を挙げることができる。
【００６２】
また、本アルゴリズムを実行できるコンピュータープログラムは、ネットワークに接続されたコンピュータの記録装置に格納しておき、ネットワークを介して他のコンピュータに転送することもできる。本アルゴリズムを実行するコンピュータープログラムを提供する提供媒体としては、様々な形式のコンピュータに読み出し可能な媒体として頒布可能であって、特定のタイプの媒体に限定されるものではない。
【００６３】
一方、本アルゴリズムは、コンピュータとして実現することができる。この場合コンピュータは、本アルゴリズムを実行できるコンピュータープログラムを格納したハードディスク装置等の補助記憶装置と、サンプルの遺伝子型に関する情報を入力できる入力装置と、ハプロタイプ頻度を算出する演算装置と、算出したハプロタイプ頻度を表示する表示装置とを備える。
【００６４】
【実施例】
以下、実施例を用いて本発明を更に詳細に説明するが、本発明の技術的範囲はこれらの実施例に限定されるものではない。
【００６５】
〔実施例１〕
ＳＡＡ 遺伝子におけるハプロタイプ頻度の推定
ＳＡＡ遺伝子については、６つのＳＮＰ座位についての１５６人の被験者から得たハプロタイプデータが開示されている（ＭｏｒｉｇｕｃｈｉらＡｒｔｈｒｉｔｉｓ Ｒｈｅｕｍ ４４：１２６６−１２７２ （２００１））。当該論文において、全被験者についてディプロタイプ構成が決定されており、これらのデータが実際のデータを反映しているものと判断された。
【００６６】
そこで当該論文に開示されているデータを適宜用いて、上述したアルゴリズムを用いたハプロタイプ頻度の推定を行った。ハプロタイプ頻度の推定を行うに際して、当該論文に開示されているハプロタイプデータを混合するとともに、相データを除いた。
【００６７】
多重連鎖座位における相未知の遺伝子型データを用いて、上述したアルゴリズムを実装した「ｌｄｐｏｏｌｅｄ」プログラムによりハプロタイプの頻度推定を行った。
【００６８】
表１には、必要なＣＰＵ時間、収束までの反復試行数、ロッドスコア、χ２値、及びＰ値を示す。それぞれのＰ値は、６つの座位の全てが互いに独立であるという帰無仮説を棄却することの危険率（ｒｉｓｋ）を表す。
【００６９】
【表１】
なお、表１において、Ｌｏｄ値及びχ２値の計算方法は、上述した方法に準じた。自由度は、上記式（９）により算出した場合、５７となるはずである。各プールに含まれる被験者の数が増加するにつれて、ロッドスコア及びχ２値は低下したが、Ｐ値は増大した（表１）。これらの結果は、おそらくは、プーリングのために情報内容が減少することによって生じるものである。
【００７１】
表２には、１人、２人及び４人の被験者から得た遺伝子型データを含む各プールについて、異なる頻度推定プロトコールを用いて推定したハプロタイプ頻度の結果を示す。各プールが２人以上の被験者から得たＤＮＡを含む場合には、各被験者についてのオリジナルデータからサンプリングされる被験者の様々な組み合わせを利用して推定した相対頻度の、平均及び標準偏差も示した。
【００７２】
【表２】
なお、表２における「ハプロタイプ」には、少なくとも１種の推定方法において相対頻度が＞０．００１となるハプロタイプが含まれる。表２における「相対頻度」は１回の推定を行って得られた推定頻度である。表２における「平均」及び「標準偏差」は種々の無作為抽出法を用いてプールを作製して得られた推定頻度の平均及び標準誤差である。
【００７４】
被験者を一人しか含まれないプール（１個体／プール）についての頻度推定の結果は、ＤＮＡがプールされなかった既報論文中で示された結果（ＫｉｔａｍｕｒａらＡｎｎ． Ｈｕｍ． Ｇｅｎｅｔ． ６６：１８３−１９３（２００２））と同じであった。表２から判るように、プール中の被験者数の変化とともに推定したハプロタイプ頻度（相対頻度）が変化したが、プールに含まれる被験者の数が４人までの範囲では、主たるハプロタイプであるＡＣＴＧＣＣ、ＡＣＣＧＴＣ、ＡＧＣＧＣＴの相対頻度は依然として適正であった。
【００７５】
また、表２から判るように、上記の主たる３つのハプロタイプにおいては、プールを作成するために様々なランダムサンプリングを行うことによって推定した相対頻度に対する標準偏差が平均の１０％未満であった。しかしながら、マイナーなハプロタイプ（相対頻度が０．１未満）については、頻度推定は正確ではなかった。例えば、ハプロタイプＡＣＴＧＴＣの頻度は、４人の被験者を含むプールについての頻度推定によっては０．０と推定されたが、被験者を１人含むプールについての頻度推定によっては０．０１３であった。後者の頻度推定が正確であれば、このハプロタイプは、１５６人（すなわちハプロタイプコピーが３１２）の間で４回しか出現しないこととなる。
【００７６】
また、上述した方法に準じて、表２に示した相対頻度から算出したＤ値及びＤ’値を表３に示す。
【００７７】
【表３】
なお、表３においてＤ値及びＤ’値は１回の推定を行って得られた値であり、括弧内の値は無作為抽出法によって算出される「平均±標準偏差」を示す。この場合、ｆｉｊ値、すなわち２つの遺伝子座についてのハプロタイプの推定頻度は、上述した方法に記載したように、６つの遺伝子座についてのハプロタイプの推定頻度から計算した。ｆｉｊを２つの遺伝子座についての遺伝子型データから算出した場合、Ｄ値及びＤ’値が非常に類似している場合もあったが、その２つの方法が全く異なる値をもたらした場合もあった（データは示していない）。また、表３において、各プールが２人以上の被験者を含む場合は、異なる組み合わせの被験者を用いて算出した値の平均及び標準偏差を示した。この値は、異なる頻度推定プロトコールの間で異なっていたが、プールに含まれる被験者の数が最大４人までであって、かつ｜Ｄ｜値が０．１を上回る限りは、依然としてかなり一致した値を示した。しかしながら、いくつかの場合においては、プールを作成するために様々なランダムサンプリングを行うことによって算出したＤ及びＤ’の標準偏差が、平均の約５０％であり、異なる頻度推定プロトコールにより算出される値は互いに大きく異なっていた。｜Ｄ｜値は０．１を上回るが｜Ｄ’｜値は０．１を上回らない場合は、Ｄ’値がそれほど異なっていなかった点は興味深い。
【００７９】
一方、各ＤＮＡプールに含まれるハプロタイプコピーの組み合わせについて頻度推定を行った結果を表４に示す。各プロトコールについて、最初の１２人の被験者に対応するデータ部分だけを示した。
【００８０】
【表４】
なお、表４において「確率」は、プールに対するハプロタイプの各組合せの事後確率を意味する。この事後確率（表４における「確率」）によって推定ハプロタイプ頻度が示される。
【００８２】
表４から判るように、大部分のプールにおいて、最も確率の高いハプロタイプコピーの組み合わせについての事後確率が１又はほぼ１であったことを示している。ハプロタイプコピーの頻度推定を行った組み合わせについてその内容を比較すると、異なる頻度推定プロトコールの間で、多くの場合、その内容は互いに一致していた。例えば、表４から判るように、２人の被験者を含むプールを用いるプロトコールの場合のプール番号１の内容は、被験者を１人含むプールを用いた場合のプール番号１及びプール番号２の内容を組み合わせたものと同となっていた。なお、あるプロトコールによって頻度推定を行ったプールの内容は、別のプロトコールによって頻度推定を行ったプールの内容と一致しない場合もあった。
【００８３】
一方、上述したように推定したハプロタイプ頻度には誤差があるため、上述した方法に準じてブートストラップ法による標準誤差を算出した。上述したハプロタイプ頻度の推定の際と同様に、１５６人の被験者から得たＳＡＡ遺伝子についての遺伝子型データを用いて、本発明者らは上述の通りに１人、２人、又は４人の被験者を含むプールを作成した。その後、上述した方法に準じてブートストラップ法に適用した。結果を図１に示す。
【００８４】
なお、図１には、１人、２人又は４人の被験者を含むプールを用いて推定したハプロタイプ頻度の平均及び標準誤差を示す。ブートストラップサンプリングは、各ＤＮＡプールを用いたハプロタイプ頻度の推定について、それぞれ１０，０００回繰り返した。図１から判るように、各ＤＮＡプールに含まれる被験者の数に拘わらず、ハプロタイプＡＣＴＧＣＣ、ＡＣＣＧＴＣ及びＡＧＣＧＣＴについての推定頻度は、かなり安定している。また、図１に示すように、これらのハプロタイプに関する標準誤差を示すバーの長さは、平均を示すバーと比較してかなり短い。さらに、図１から判るように、各ＤＮＡプールについて推定したハプロタイプ頻度の平均は、同一のハプロタイプについてはほぼ同じ値を示した。しかしながら、マイナーなハプロタイプ（ｐｉ＜ ０．１）については、各ＤＮＡプールについて推定したハプロタイプ頻度の平均が有意に異なっており、前記誤差を示すバーは平均と比較するとかなり長かった。
【００８５】
なお、表１に示したように、各計算を行うのに必要な時間及びメモリーを記録した。１ＧＨｚのＰｅｎｔｉｕｍ ＩＩＩ（ＣＰＵ）及び１．５ＧＢメモリーを搭載したコンピュータを用いる場合、プールに含める被験者の数は、遺伝子座の数が６であれば最大で６人であった。遺伝子座の数が１３であれば、プールに含む被験者は最大でたった２人となる可能性がある。
【００８６】
一方で、複数のＤＮＡサンプルをプールしない場合（換言すれば１個体の被験者を含むＤＮＡプールの場合）には、２５個の遺伝子座について計算することが可能であった。これは、上述したアルゴリズムを実装したコンピュータープログラムが、各ＤＮＡプールに含まれるハプロタイプコピーの可能な組み合わせを用いるものであって、このステップが大量のメモリーを消費するためである。そのような組み合わせの数は、所定の遺伝子座に存在する対立遺伝子の数の検定力関数に従って増加し、かつ、ＤＮＡプールに含まれる被験者数の階乗として増加する。従って、コンピュータを用いた計算に要する時間及びメモリーがこれらの要因に依存するのは、上記のような組み合わせのための余裕が必要となるためである。
【００８７】
〔実施例２〕
実施例２では、メチレンテトラヒドロ葉酸レダクターゼ（ＭＴＨＦＲ）遺伝子における連鎖する２つの遺伝子座についてハプロタイプ頻度を推定した。ＭＴＨＦＲ遺伝子は、葉酸代謝に関係するＭＴＨＦＲ酵素をコードしている。ＭＴＨＦＲ遺伝子におけるこれら遺伝子座に関する情報は、計８０人の被験者に関する発表済データ（ＵｒａｎｏらＰｈａｒｍａｃｏｇｅｎｅｔｉｃｓ １２：１８３−１９０ （２００２））を用いた。これらの発表済データを用いてＤＮＡプールを作成し、プールされたデータからパラメーターを算出した。結果を表５に示す。
【００８８】
【表５】
なお、表５における「相対頻度」は１回の推定を行って得られた推定頻度である。表５における「平均」及び「標準偏差」は種々の無作為抽出法を用いてプールを作製して得られた推定頻度の平均及び標準誤差である。
【００９０】
表５に示したように、ＤＮＡプールに含まれる被験者の人数に拘わらず、推定されたハプロタイプ頻度が安定していた。言い換えると、推定されたハプロタイプ頻度は、４人の被験者を含むＤＮＡプールを使用した場合でさえかなり正確だったことを示唆している。これは、これらの遺伝子座に関する連鎖不平衡が非常に強いためであると考えられる。また、表５からは、相対頻度が０．１を上回るハプロタイプ同士については、頻度推定がかなり正確であったことも示している。
【００９１】
図２は、ＭＴＨＦＲ遺伝子のハプロタイプデータに対するブートストラップ法の結果を示す。これらのデータは、本アルゴリズムによれば、ハプロタイプの頻度がかなり高い場合にはＭＴＨＦＲ遺伝子のハプロタイプ頻度をかなり正確に推定できることを示唆している。
【００９２】
ＭＴＨＦＲ遺伝子について推定したハプロタイプ頻度から、実施例１と同様に計算したＤ値及びＤ’値を表６に示す。また、２人及び４人の被験者を含むプールを用いるプロトコールによって得たデータに関して、様々なランダムサンプリングを行ったものに対する平均及び標準偏差も示す。
【００９３】
【表６】
【００９４】
〔実施例３〕
実施例３では、Ｎ−アセチルトランスフェラーゼ２（ＮＡＴ２）遺伝子における連鎖する７つの遺伝子座についてハプロタイプ頻度を推定した。ＮＡＴ２遺伝子は、Ｎ−アセチル残基の転移に関与するＮＡＴ２酵素をコードしている。ＮＡＴ２遺伝子におけるこれら遺伝子座に関する情報は、計１１６人の被験者に関する発表済データ（ＴａｎａｋａらＪ． Ｒｈｅｕｍａｔｏｌ （Ｉｎ ｐｒｅｓｓ））を用いた。これらの発表済データを用いてＤＮＡプールを作成し、プールされたデータからパラメーターを算出した。結果を表７に示す。
【００９５】
【表７】
なお、表７における「ハプロタイプ」には、少なくとも１種の推定方法において相対頻度が＞０．００１となるハプロタイプが含まれる。表７における「相対頻度」は１回の推定を行って得られた推定頻度である。表７における「平均」及び「標準偏差」は種々の無作為抽出法を用いてプールを作製して得られた推定頻度の平均及び標準誤差である。
【００９７】
表７に示したように、ＤＮＡプールに含まれる被験者の人数に拘わらず、推定されたハプロタイプ頻度が安定していた。言い換えると、推定されたハプロタイプ頻度は、４人の被験者を含むＤＮＡプールを使用した場合でさえかなり正確だったことを示唆している。
【００９８】
図３は、ＮＡＴ２遺伝子のハプロタイプデータに対するブートストラップ法の結果を示す。これらのデータは、本アルゴリズムによれば、ハプロタイプの頻度がかなり高い場合にはＭＴＨＦＲ遺伝子のハプロタイプ頻度をかなり正確に推定できることを示唆している。
【００９９】
ＮＡＴ２遺伝子について推定したハプロタイプ頻度から、実施例１と同様に計算したＤ値及びＤ’値を表８に示す。
【０１００】
【表８】
【０１０１】
〔実施例４〕
実施例４では、スムーセリン（ｓｍｏｏｔｈｅｌｉｎ）遺伝子における連鎖する３６個の遺伝子座についてハプロタイプ頻度を推定した。ｓｍｏｏｔｈｅｌｉｎ遺伝子は、第２２染色体に存在し、平滑筋細胞骨格を構成するｓｍｏｏｔｈｅｌｉｎをコードする遺伝子である。ｓｍｏｏｔｈｅｌｉｎ遺伝子におけるこれら遺伝子座に関する情報は、計１０２人の日本人の被験者に関するデータを用いた。
【０１０２】
本例では、遺伝子座の数が３６箇所と多かったため、全遺伝子座を３つの領域に任意に分けてハプロタイプ頻度の推定を行った。すなわち、ｓｍｏｏｔｈｅｌｉｎ遺伝子における５’末端から３’末端に向かって、１３個の遺伝子座を含む第１の領域、１３個の遺伝子座を含む第２の領域、及び１０個の遺伝子座を含む第３の領域となるように全遺伝子座を分けた。
【０１０３】
これら各領域について、前記データを用いてＤＮＡプールを作成し、プールされたデータからハプロタイプの相対頻度を算出した。結果を表９に示す。
【０１０４】
【表９】
なお、表９において「ハプロタイプ」は、１個体を含むＤＮＡプールを用いた相対頻度が≧０．００１となるハプロタイプが含まれる。また、表９の「ハプロタイプ」欄における「Ｄ」は欠失を意味し、「Ｉ」は挿入を意味する。
【０１０６】
表９に示すように、１個体を含むＤＮＡプール及び２個体を含むＤＮＡプールのいずれを用いた場合でも、推定したハプロタイプの頻度（表９における相対頻度）は、そのハプロタイプ頻度が０．１を上回る限りにおいては、非常に近い値を示した。しかしながらマイナーなハプロタイプ（相対頻度が０．１未満）については、１個体を含むＤＮＡプールを用いた場合と２個体を含むＤＮＡプールを用いた場合とで類似しているとは言えなかった。
【０１０７】
次に、本例では、連鎖不平衡の強さを評価するため、上記方法に準じてペアワイズ方式で連鎖不平衡測定値ρ２を算出した。３６個の遺伝子座に関し（３６×３５）／２＝６３０ペアについて、算出したρ２を図４に示す。なお、本例では、ｆｉｊ値は２つの遺伝子座のみの遺伝子型データを用いて算出した。図４から判るように、多少の差異は観察されたものの、１個体を含むＤＮＡプールを用いて算出されたρ２値と２個体を含むＤＮＡプールを用いて算出されたρ２値とは、互いに非常に類似していた。様々なペアの遺伝子座について算出したρ２値の平均及び標準偏差は、１個体を含むＤＮＡプールを用いた場合はそれぞれ０．１１４及び０．２２４であり、２個体を含むＤＮＡプール用いた場合はそれぞれ０．１１８及び０．２２７であった。１個体を含むＤＮＡプールを用いて算出されたρ２値と２個体を含むＤＮＡプールを用いて算出されたρ２値との差の絶対値の平均及び標準偏差は、それぞれ０．０１１及び０．０２２であった。また、１個体を含むＤＮＡプールを用いて算出されたρ２値と２個体を含むＤＮＡプールを用いて算出されたρ２値との差の絶対値は、平均の約９．５％であった。
【０１０８】
【発明の効果】
以上、詳細に説明したように、本発明に係るハプロタイプ頻度推定方法によれば、複数のプールした塩基配列データを入力値として集団のハプロタイプ頻度を推定することができる。換言すると、本発明に係るハプロタイプ頻度推定方法は、プールすることなしに複数の塩基配列データをそれぞれ入力値として用いた場合と同様に、ハプロタイプ頻度を推定することができる。
【図面の簡単な説明】
【図１】ＳＡＡ遺伝子遺伝子についての遺伝子型データに関して１人、２人又は４人の被験者を含むプールを用いて、ブートストラップ法を適用して算出した、推定ハプロタイプ頻度の平均及び標準誤差を示す特性図である。
【図２】ＭＴＨＦＲ遺伝子についての遺伝子型データに関して１人、２人又は４人の被験者を含むプールを用いて、ブートストラップ法を適用して算出した、推定ハプロタイプ頻度の平均及び標準誤差を示す特性図である。
【図３】ＮＡＴ２遺伝子遺伝子についての遺伝子型データに関して１人、２人又は４人の被験者を含むプールを用いて、ブートストラップ法を適用して算出した、推定ハプロタイプ頻度の平均及び標準誤差を示す特性図である。
【図４】ｓｍｏｏｔｈｅｌｉｎ遺伝子についての遺伝子型データに関して、ペアワイズ方式で連鎖不平衡測定値ρ２を算出した結果を示す特性図である。
【図５】ＥＭアルゴリズムの各ステップを概略的に示すフローチャートである。
Claims (12)
- 集団に含まれる複数の検体に関する遺伝子型情報を集積してある遺伝子型プール情報を入力値として、期待値最大化アルゴリズムにより上記集団におけるハプロタイプ頻度を推定する方法。
- 上記複数の検体に由来するＤＮＡを含む複数のＤＮＡプールを調製し、当該複数のＤＮＡプールに含まれるＤＮＡにおける所定の遺伝子座に関するタイピングを行う工程を含むことを特徴とする請求項１記載のハプロタイプ頻度を推定する方法。
- ２以上６以下の検体に関する遺伝子型情報を集積して遺伝子型プール情報を作成することを特徴とする請求項１記載のハプロタイプ頻度を推定する方法。
- 前記期待値最大化アルゴリズムはＥＭアルゴリズムであり、前記遺伝子型プール情報に最尤推定を行うことによりハプロタイプ頻度を推定することを特徴とする請求項１記載のハプロタイプ頻度を推定する方法。
- 集団に含まれる複数の検体に関する遺伝子型情報を集積してある遺伝子型プール情報を入力値として、期待値最大化アルゴリズムにより上記集団におけるハプロタイプ頻度を推定するアルゴリズム。
- 上記複数の検体に由来するＤＮＡを含む複数のＤＮＡプールを調製し、当該複数のＤＮＡプールに含まれるＤＮＡにおける所定の遺伝子座に関するタイピングを行う工程を含むことを特徴とする請求項５記載のハプロタイプ頻度を推定するアルゴリズム。
- ２以上６以下の検体に関する遺伝子型情報を集積して遺伝子型プール情報を作成することを特徴とする請求項５記載のハプロタイプ頻度を推定するアルゴリズム。
- 前記期待値最大化アルゴリズムはＥＭアルゴリズムであり、前記遺伝子型プール情報に最尤推定を行うことによりハプロタイプ頻度を推定することを特徴とする請求項５記載のハプロタイプ頻度を推定するアルゴリズム。
- 集団に含まれる複数の検体に関する遺伝子型情報を集積してある遺伝子型プール情報を入力値として、期待値最大化アルゴリズムにより上記集団におけるハプロタイプ頻度を推定するプログラム。
- 上記複数の検体に由来するＤＮＡを含む複数のＤＮＡプールを調製し、当該複数のＤＮＡプールに含まれるＤＮＡにおける所定の遺伝子座に関するタイピングを行う工程を含むことを特徴とする請求項９記載のハプロタイプ頻度を推定するプログラム。
- ２以上６以下の検体に関する遺伝子型情報を集積して遺伝子型プール情報を作成することを特徴とする請求項９記載のハプロタイプ頻度を推定するプログラム。
- 前記期待値最大化アルゴリズムはＥＭアルゴリズムであり、前記遺伝子型プール情報に最尤推定を行うことによりハプロタイプ頻度を推定することを特徴とする請求項９記載のハプロタイプ頻度を推定するプログラム。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2002303395A JP2004192018A (ja) | 2002-10-16 | 2002-10-17 | Ｄｎａプールによるハプロタイプ頻度推定方法 |
Applications Claiming Priority (2)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2002302264 | 2002-10-16 | ||
| JP2002303395A JP2004192018A (ja) | 2002-10-16 | 2002-10-17 | Ｄｎａプールによるハプロタイプ頻度推定方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| JP2004192018A true JP2004192018A (ja) | 2004-07-08 |
Family
ID=32774396
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP2002303395A Pending JP2004192018A (ja) | 2002-10-16 | 2002-10-17 | Ｄｎａプールによるハプロタイプ頻度推定方法 |
Country Status (1)
| Country | Link |
|---|---|
| JP (1) | JP2004192018A (ja) |
Cited By (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2006059770A1 (ja) * | 2004-12-03 | 2006-06-08 | Riken | 確率算出方法 |
| JP2007233485A (ja) * | 2006-02-27 | 2007-09-13 | Fujitsu Ltd | 遺伝子多型解析支援プログラム、該プログラムを記録した記録媒体、遺伝子多型解析支援装置、および遺伝子多型解析支援方法 |
| WO2008136210A1 (ja) * | 2007-04-26 | 2008-11-13 | Riken | ハプロタイプ推定装置、および、プログラム |
| JP2009205551A (ja) * | 2008-02-28 | 2009-09-10 | Institute Of Physical & Chemical Research | ハプロタイプ推定装置、および、プログラム |
| CN108241792A (zh) * | 2016-12-23 | 2018-07-03 | 深圳华大基因科技服务有限公司 | 一种整合多平台基因分型结果的方法和装置 |
-
2002
- 2002-10-17 JP JP2002303395A patent/JP2004192018A/ja active Pending
Cited By (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2006059770A1 (ja) * | 2004-12-03 | 2006-06-08 | Riken | 確率算出方法 |
| JP4755601B2 (ja) * | 2004-12-03 | 2011-08-24 | 独立行政法人理化学研究所 | 確率算出方法 |
| JP2007233485A (ja) * | 2006-02-27 | 2007-09-13 | Fujitsu Ltd | 遺伝子多型解析支援プログラム、該プログラムを記録した記録媒体、遺伝子多型解析支援装置、および遺伝子多型解析支援方法 |
| WO2008136210A1 (ja) * | 2007-04-26 | 2008-11-13 | Riken | ハプロタイプ推定装置、および、プログラム |
| JP5129809B2 (ja) * | 2007-04-26 | 2013-01-30 | 独立行政法人理化学研究所 | ハプロタイプ推定装置、および、プログラム |
| US8401802B2 (en) | 2007-04-26 | 2013-03-19 | Riken | Haplotype estimating apparatus and program |
| JP2009205551A (ja) * | 2008-02-28 | 2009-09-10 | Institute Of Physical & Chemical Research | ハプロタイプ推定装置、および、プログラム |
| CN108241792A (zh) * | 2016-12-23 | 2018-07-03 | 深圳华大基因科技服务有限公司 | 一种整合多平台基因分型结果的方法和装置 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Choin et al. | Genomic insights into population history and biological adaptation in Oceania | |
| Weissensteiner et al. | Contamination detection in sequencing studies using the mitochondrial phylogeny | |
| US20250022535A1 (en) | Identifying Genetic Variants by Imputation | |
| Speidel et al. | A method for genome-wide genealogy estimation for thousands of samples | |
| Aulchenko et al. | Linkage disequilibrium in young genetically isolated Dutch population | |
| Rasmussen et al. | Genome-wide inference of ancestral recombination graphs | |
| Uitterlinden | An introduction to genome-wide association studies: GWAS for dummies | |
| Patin et al. | The impact of agricultural emergence on the genetic history of African rainforest hunter-gatherers and agriculturalists | |
| Tsai et al. | Comparison of three methods to estimate genetic ancestry and control for stratification in genetic association studies among admixed populations | |
| Zaitlen et al. | Leveraging population admixture to characterize the heritability of complex traits | |
| US20140067355A1 (en) | Using Haplotypes to Infer Ancestral Origins for Recently Admixed Individuals | |
| US20190338349A1 (en) | Methods and systems for high fidelity sequencing | |
| BR112016007401B1 (pt) | Método para determinar a presença ou ausência de uma aneuploidia cromossômica em uma amostra | |
| JP2003534560A (ja) | ハプロタイプ頻度の推定値を用いるｄｎａマーカーに基づく遺伝子解析の方法およびその使用 | |
| Li et al. | Single nucleotide polymorphism (SNP) detection and genotype calling from massively parallel sequencing (MPS) data | |
| Liu et al. | Systematic assessment of imputation performance using the 1000 Genomes reference panels | |
| Zhang et al. | Rapid and robust resampling-based multiple-testing correction with application in a genome-wide expression quantitative trait loci study | |
| JP2004192018A (ja) | Ｄｎａプールによるハプロタイプ頻度推定方法 | |
| Yu et al. | Detecting natural selection by empirical comparison to random regions of the genome | |
| CN120322562A (zh) | 基因组结构变异的检测方法、系统、设备及介质 | |
| Lin et al. | Smarter clustering methods for SNP genotype calling | |
| McCallum et al. | Empirical Bayes scan statistics for detecting clusters of disease risk variants in genetic studies | |
| Wang et al. | Estimating genome-wide copy number using allele-specific mixture models | |
| Zhu et al. | The analysis of ethnic mixtures | |
| Tanck et al. | Simultaneous estimation of gene‐gene and gene‐environment interactions for numerous loci using double penalized log–likelihood |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A621 | Written request for application examination |
Free format text: JAPANESE INTERMEDIATE CODE: A621 Effective date: 20051014 |
|
| A711 | Notification of change in applicant |
Free format text: JAPANESE INTERMEDIATE CODE: A711 Effective date: 20060104 |
|
| A521 | Written amendment |
Free format text: JAPANESE INTERMEDIATE CODE: A821 Effective date: 20060104 |
|
| A977 | Report on retrieval |
Free format text: JAPANESE INTERMEDIATE CODE: A971007 Effective date: 20070521 |
|
| A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20070605 |
|
| A521 | Written amendment |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20070806 |
|
| A02 | Decision of refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A02 Effective date: 20071120 |