JP5134397B2 - ハプロタイプ推定装置、および、プログラム - Google Patents

Description

本発明は、ハプロタイプ推定装置、および、プログラムに関する。

生活習慣病のような複雑な病気の原因となる遺伝子の解明や、さらには個別化医療の実現の為には、遺伝子型データなどの実験データから、ヒト等の個体のハプロタイプを推定する必要がある。

複数座位に渡るハプロタイプを推定する方法には、主に、（１）１座位当たりの遺伝子型データから複数座位に渡るハプロタイプとその頻度を推定する方法と、（２）塩基多型を考慮したコピー多型（Ｃｏｐｙ Ｎｕｍｂｅｒ Ｐｏｌｙｍｏｒｐｈｉｓｍ，あるいはＣｏｐｙ Ｎｕｍｂｅｒ Ｖａｒｉａｔｉｏｎ、本明細書中で「ＣＮＰ」と略す場合がある。）に関する実験データから一塩基多型のハプロタイプとその頻度を推定する方法と、が存在する。

まず、上述の（１）に示す、１座位当たりの遺伝子型データから複数座位に渡るハプロタイプとその頻度を推定する方法において、「１座位当たりの遺伝子型データ」とは、複数座位に渡る関連が分からない、座位毎に独立した（相（ｐｈａｓｅ）が未知の）遺伝子型データのことである。ここで、図１は、１座位あたりの遺伝子型データの一例を示す図である。図１において、Ｌは座位（Ｌｏｃｕｓ）、Ａはアレル（Ａｌｌｅｌｅ）を表している。

図１に示すように、１座位当たりの遺伝子型データは、各個体、各座位における、各アレルのカウント数データを有している。ここで、「カウント数データ」とは、個体において、各座位におけるアレルをカウントすることにより得られたカウント数のデータである。例えば、図１においては、個体１における座位（Ｌ_１）におけるアレル（Ａ_１）のカウント数は「１」であり、座位（Ｌ_３）におけるアレル（Ａ_１）のカウント数は「２」である。

ここで、上述のように、遺伝子型データは、直接的に相（ｐｈａｓｅ）を特定するものではなく、座位間のアレルの関連は未知であり、図１の例（例えば、個体１）でいえば、カウント数データから、座位Ｌ_１と座位Ｌ_２の間での相を特定することはできず、座位Ｌ_１におけるアレル（Ａ_１／Ａ_２）と座位Ｌ_２におけるアレル（Ａ_１／Ａ_３）の関連は未知である。そのため、ハプロタイプを推定する（相を特定する）ための手法が必要となる。

ここで、非特許文献１〜４に記載のハプロタイプ推定方法は、１座位あたりの遺伝子型データから、複数座位に渡るハプロタイプを推定する。ここで、「複数座位に渡るハプロタイプ」とは、複数座位に渡るアレルの組合せ（相を特定する組合せ）のことである。ここで、図２は、複数座位に渡るハプロタイプの組合せの一例を示す図である。図２において、Ａ（Ｌ）は、座位Ｌに対応するアレルＡを表している。

図２に示すように、例えば、ハプロタイプ１は、座位Ｌ_１においてＡ_１のアレルをもち、座位Ｌ_２においてＡ_１のアレルをもち、座位Ｌ_３においてＡ_１のアレルを持つことが特定されている。このように、従来のハプロタイプの推定方法においては、一般に２種類のアレルを想定し、一塩基多型（Ｓｉｎｇｌｅ Ｎｕｃｌｅｏｔｉｄｅ Ｐｏｌｙｍｏｒｐｈｉｓｍ、「ＳＮＰ」と略す。）などの遺伝子型データから、複数座位に渡るハプロタイプを推定している。また、この方法では、１座位当たりの遺伝子型データさえ分かれば、コピー数多型（Ｃｏｐｙ Ｎｕｍｂｅｒ Ｐｏｌｙｍｏｒｐｈｉｓｍ，またはＣｏｐｙ Ｎｕｍｂｅｒ Ｖａｒｉａｔｉｏｎ、以下「ＣＮＰ」と略す。）のデータからでもハプロタイプ推定を行うことができる。例えば、１座位あたりの遺伝子型データとして、コピーの数が多いアレル、少ないアレルの２種類を定義し（その数は曖昧ではあるが）、一塩基多型に対してと同じ方法で複数座位にわたるハプロタイプが推定されている（非特許文献５参照）。

次に、上述の（２）に示す、塩基多型を考慮したコピー多型に関する実験データから一塩基多型のハプロタイプとその頻度を推定する方法（例えば、本願出願人および発明者による特願２００７−２３７１３９（未公開であるため先行技術ではない。））に関して、コピー数多型は、ある区間の配列（「コピー単位」と呼ぶ。）が繰り返し現れることがあり、そのコピー数に個体差がある。ここで、「塩基多型を考慮したコピー数多型に関する実験データ」とは、蛍光色素プローブ等によって実験的に識別される異なる塩基の数のデータのことであり、そのデータは生物学的にはコピー単位上にある異なる塩基の数に相当し、例えば図３のように示される。ここで、図３は、塩基多型を考慮した一塩基多型の実験データの一例を示す図である。

図３において、Ｍは、蛍光色素プローブ等の標識で識別される個体間で違いのない配列部位（マーカー部位）、Ｆは、（異なる蛍光色素等で区別される）個体間で違いうる塩基（多型塩基）に、それぞれ対応する。換言すれば、図３は、各個体、各マーカー部位における、各多型塩基をカウントしたカウント数の表であり、図１の場合とは異なり、各マーカー部位におけるカウント総数が２以外の数値（０，１，２，３，４・・・）であることが想定されている。すなわち、多型塩基のカウント数は、個体における２つのハプロタイプ（すなわちディプロタイプ）のコピー数に依存する。ここで「コピー数多型のハプロタイプ」とは、相同染色体の１本上の区画における異なる塩基の組合せで表されたコピー単位の組合せのことであり、例えば、図４のように表現される。ここで、図４は、図３のデータから推定されるハプロタイプの表現例を示す図である。図４において、Ｆ（Ｍ）は、マーカー部位Ｍに対応する異なる塩基Ｆを表し、各コピー単位はＭに渡るＦ（Ｍ）の組合せで表され、異なるコピー単位は“，”で区切られている。

図４に示すように、ハプロタイプ表現形式は、塩基多型とコピー数多型とを同時に表現され、例えばハプロタイプ１ではコピー数が２であり、それぞれのコピー単位上ではＦ_１（Ｍ_１）Ｆ_１（Ｍ_２）Ｆ_１（Ｍ_３）とＦ_１（Ｍ_１）Ｆ_２（Ｍ_２）Ｆ_３（Ｍ_３）の多型塩基の組合せを有している。

チアンフア ニウ（Ｔｉａｎｈｕａ Ｎｉｕ）著 「アルゴリズムズ フォー インファリング ハプロタイプス（Ａｌｇｏｒｉｔｈｍｓ ｆｏｒ ｉｎｆｅｒｒｉｎｇ ｈａｐｌｏｔｙｐｅｓ）」 Ｇｅｎｅｔ Ｅｐｉｄｅｍｉｏｌ．、２００４年１２月、２７巻（４号）３３４−３４７頁 ジャオフイＳキン、チアンフア ニウ、ジュンＳリウ（Ｚｈａｏｈｕｉ Ｓ． Ｑｉｎ， Ｔｉａｎｈｕａ Ｎｉｕ， Ｊｕｎ Ｓ． Ｌｉｕ）著 「パーティション−ライゲーション−エクスペクテーション−マキシマイゼーション アルゴリズム フォー ハプロタイプ インフェアレンス ウィズ シングル−ヌクレオタイド ポリモーフィズムズ（Ｐａｒｔｉｔｉｏｎ−ｌｉｇａｔｉｏｎ−ｅｘｐｅｃｔａｔｉｏｎ−ｍａｘｉｍｉｚａｔｉｏｎ ａｌｇｏｒｉｔｈｍ ｆｏｒ ｈａｐｌｏｔｙｐｅ ｉｎｆｅｒｅｎｃｅ ｗｉｔｈ ｓｉｎｇｌｅ−ｎｕｃｌｅｏｔｉｄｅ ｐｏｌｙｍｏｒｐｈｉｓｍｓ）」Ａｍ Ｊ Ｈｕｍ Ｇｅｎｅｔ．、２００２年１１月、７１巻（５号）１２４２−１２４７頁 ローレント エクスコファー、モンゴメリー スラキン（Ｌａｕｒｅｎｔ Ｅｘｃｏｆｆｉｅｒ， Ｍｏｎｔｇｏｍｅｒｙ Ｓｌａｔｋｉｎ）著 「マキシマム−ライクリフッド エスティメーション オブ モレキュラー ハプロタイプ フロークェンシーズ イン ア ディプロイド ポピュレーション（Ｍａｘｉｍｕｍ−ｌｉｋｅｌｉｈｏｏｄ ｅｓｔｉｍａｔｉｏｎ ｏｆ ｍｏｌｅｃｕｌａｒ ｈａｐｌｏｔｙｐｅ ｆｒｅｑｕｅｎｃｉｅｓ ｉｎ ａ ｄｉｐｌｏｉｄ ｐｏｐｕｌａｔｉｏｎ）」 Ｍｏｌ Ｂｉｏｌ Ｅｖｏｌ．、１９９５年９月、１２巻（５号）、９２１−９２７頁 Ｍ．Ｅ．ホーリー、Ｋ．Ｋ．キッド（Ｍ． Ｅ． Ｈａｗｌｅｙ， Ｋ． Ｋ． Ｋｉｄｄ）著 「ハプロ： ア プログラム ユージング ザ ＥＭ アルゴリズム トゥー エスティメイト ザ フリークェンシーズ オブ マルチ−サイト ハプロタイプズ（ＨＡＰＬＯ： ａ ｐｒｏｇｒａｍ ｕｓｉｎｇ ｔｈｅ ＥＭ ａｌｇｏｒｉｔｈｍ ｔｏ ｅｓｔｉｍａｔｅ ｔｈｅ ｆｒｅｑｕｅｎｃｉｅｓ ｏｆ ｍｕｌｔｉ−ｓｉｔｅ ｈａｐｌｏｔｙｐｅｓ）」 Ｊ Ｈｅｒｅｄ．、１９９５年９−１０月、８６巻（５号）、４０９−４１１頁 リチャード レドン、シュンペイ イシカワ、カレンＲフィッチ他、全４３名（Ｒｉｃｈａｒｄ Ｒｅｄｏｎ， Ｓｈｕｍｐｅｉ Ｉｓｈｉｋａｗａ， Ｋａｒｅｎ Ｒ． Ｆｉｔｃｈ， ｅｔ ａｌ．）著 「グローバル バリエーション イン コピー ナンバー イン ザ ヒューマン ゲノム（Ｇｌｏｂａｌ ｖａｒｉａｔｉｏｎ ｉｎ ｃｏｐｙ ｎｕｍｂｅｒ ｉｎ ｔｈｅ ｈｕｍａｎ ｇｅｎｏｍｅ）」 ネーチャー（Ｎａｔｕｒｅ）出版 ２００６年１１月２３日、４４４巻（７１１８号）４４４−４５４頁

しかしながら、上述した方法によっては、塩基多型が考慮されていないコピー数多型に関するデータから、コピー数多型のハプロタイプとその頻度を推定することができないという問題があった。

例えば、コピー数多型の分析において、効率性を優先する実験プラットフォームによっては、背景技術で示したようなコピー単位上の多型塩基のデータが得られず、一座位当たりの二つの相同染色体に渡る合計のコピー数や曖昧なコピー数しか実験データとして得られない場合があり、上述した方法によってはハプロタイプとその頻度を推定することができないという問題があった。また、コピー数多型に関するデータに１座位当たりの遺伝子型データが加わったデータから、ハプロタイプとその頻度を推定する方法は開発されていないという問題があった。

本発明は、上記に鑑みてなされたもので、塩基多型が考慮されていないコピー数多型に関する実験データから、曖昧なコピー数の総和しか得られない場合であっても、ハプロタイプとその頻度を高精度で推定することができる、ハプロタイプ推定装置、および、プログラムを提供することを目的とする。また、コピー数多型に関するデータに１座位当たりの遺伝子型データが加わったデータであっても、ハプロタイプとその頻度を推定することができる、ハプロタイプ推定装置、および、プログラムを提供することを目的とする。

このような目的を達成するため、請求項１に記載のハプロタイプ推定装置は、集団における各個体の少なくともコピー数多型を含む実験データからハプロタイプを推定する、制御部と記憶部を少なくとも備えたハプロタイプ推定装置において、上記記憶部は、上記個体毎に、上記実験データから得られた、標識によって特定されるマーカー部位に対応付けられた上記コピー数多型のコピー数の総和Ｎに対する一または複数の条件式を、上記マーカー部位の種類毎に記憶する多型テーブル、を備え、上記制御部は、上記コピー数の総和Ｎについて、足し合わせた場合に、上記多型テーブルに記憶された上記コピー数の総和Ｎの上記条件式を満たす任意の２つの整数に分割するコピー数総和分割手段と、上記コピー数総和分割手段により分割された上記２つの整数をそれぞれ多型識別文字として表し、上記個体において上記多型識別文字を列挙したハプロタイプ文字列の組合せとして格納するハプロタイプ文字列格納手段と、上記集団において、同一である上記ハプロタイプ文字列の数を集計し、当該ハプロタイプ文字列の上記集団における頻度を求め、当該頻度が所定の条件を満たす上記各個体の上記ハプロタイプ文字列の上記組合せを、上記ハプロタイプの組合せとして推定するハプロタイプ推定手段と、を備えたことを特徴とする。

また、請求項２に記載のハプロタイプ推定装置は、請求項１に記載のハプロタイプ推定装置において、上記条件式は、上記コピー数の総和Ｎの条件を規定する一または複数の等式または不等式を含むこと、を特徴とする。

また、請求項３のハプロタイプ推定装置は、請求項１に記載のハプロタイプ推定装置において、上記多型テーブルは、上記個体毎に、塩基多型を更に含む上記実験データから得られた、上記標識によって特定される上記マーカー部位に対応付けられた多型塩基をカウントしたカウント数を、当該多型塩基の種類毎に更に記憶し、上記ハプロタイプ文字列格納手段は、上記ハプロタイプ文字列の組合せにおいて、上記多型塩基の種類に対応付けた上記多型識別文字を更に列挙すること、を特徴とする。

また、請求項４のハプロタイプ推定装置は、請求項１に記載のハプロタイプ推定装置において、上記ハプロタイプ推定手段は、上記ハプロタイプ文字列の頻度を、ハーディ・ワインバーグの法則に基づいて算出し、上記所定の条件を、上記集団におけるハーディ・ワインバーグ平衡とすること、を特徴とする。

また、請求項５に記載のハプロタイプ推定装置は、請求項１乃至４のいずれか一つに記載のハプロタイプ推定装置において、上記ハプロタイプ文字列格納手段は、作成した上記ハプロタイプ文字列の組合せにおける上記多型識別文字が表す上記整数の和が、上記多型テーブルにおける上記コピー数の総和Ｎの上記条件式を満たすか否か確認し、上記条件式を満たさない場合に当該ハプロタイプ文字列を除外し、および／または、作成した上記ハプロタイプ文字列の組合せが他の上記ハプロタイプ文字列の組合せと重複するか否か判断し、重複する場合に当該ハプロタイプ文字列の組合せを除外すること、を特徴とする。

また、請求項６に記載のハプロタイプ推定装置は、請求項１乃至５のいずれか一つに記載のハプロタイプ推定装置において、上記ハプロタイプ推定手段は、ＥＭ（Ｅｘｐｅｃｔａｔｉｏｎ−Ｍａｘｉｍｉｚａｔｉｏｎ）法を用いて、上記集団における上記ハプロタイプ文字列の頻度を、当該ハプロタイプ文字列を少なくとも一方に有する上記組合せの頻度により重み付けして算出するＭステップと、上記組合せの頻度を、当該組合せを構成する上記ハプロタイプ文字列の頻度の積により求め、当該組合せの頻度に基づいて上記重みを算出するＥステップと、を上記頻度の値が収束するまで交互に繰り返すハプロタイプ頻度算出手段、を更に備えたことを特徴とする。

また、請求項７に記載のハプロタイプ推定装置は、請求項６に記載のハプロタイプ推定装置において、上記Ｍステップは、下記の数式１に基づいて、上記ハプロタイプ文字列の頻度を算出し、上記Ｅステップは、下記の数式２に基づいて上記ハプロタイプ文字列の上記組合せの頻度を求め、上記重みとして当該ハプロタイプ文字列の上記組合せの頻度を上記集団における上記組合せの頻度の総和で除して算出すること、を特徴とする。

（ここで、Ｐ（ｈ_ｉ）は上記ハプロタイプ文字列の頻度を表し、ｈは上記ハプロタイプ文字列を表し、ｉは上記ハプロタイプ文字列のインデックスを表す。また、ｎは上記集団を構成する上記個体の数、ｊは上記多型テーブルにおける上記コピー数の総和Ｎの上記条件式および／または上記カウント数に関するパターンのインデックス、ｋは上記ハプロタイプ文字の組合せのインデックス、Ｎ（ｃ_ｊ）は上記パターンｊを持つ上記個体の数を表す。また、δ（ｈ_ｉ，ｄ_ｊｋ）は、上記組合せｄ_ｊｋが一方に当該ハプロタイプ文字列ｈ_ｉを有する場合に１を返し、両方に当該ハプロタイプ文字列ｈ_ｉを有する場合に２を返し、当該ハプロタイプ文字列ｈ_ｉを持たない場合に０を返す関数であり、ｄは上記ハプロタイプ文字列の上記組合せを表す。また、ｗ_ｊｋは上記ハプロタイプ文字列の上記組合せの頻度による上記重みである。）

（ここで、Ｐ（ｄ_ｊｋ）は、上記ハプロタイプ文字列の上記組合せの頻度を表す。また、ｈ_ｌおよびｈ_ｍは当該組合せを構成する２つの上記ハプロタイプ文字列を表し、Ｐ（ｈ_ｌ）およびＰ（ｈ_ｍ）は、当該２つの上記ハプロタイプ文字列の頻度をそれぞれ表す。）

また、請求項８に記載のプログラムは、集団における各個体の少なくともコピー数多型を含む実験データからハプロタイプを推定する、制御部と記憶部を少なくとも備えたハプロタイプ推定装置に実行させるためのプログラムであって、上記記憶部は、上記個体毎に、上記実験データから得られた、標識によって特定されるマーカー部位に対応付けられた上記コピー数多型のコピー数の総和Ｎに対する一または複数の条件式を、上記マーカー部位の種類毎に記憶する多型テーブル、を備えており、上記制御部において、上記コピー数の総和Ｎについて、足し合わせた場合に、上記多型テーブルに記憶された上記コピー数の総和Ｎの上記条件式を満たす任意の２つの整数に分割するコピー数総和分割ステップと、上記コピー数総和分割ステップにて分割された上記２つの整数をそれぞれ多型識別文字として表し、上記個体において上記多型識別文字を列挙したハプロタイプ文字列の組合せとして格納するハプロタイプ文字列格納ステップと、上記集団において、同一である上記ハプロタイプ文字列の数を集計し、当該ハプロタイプ文字列の上記集団における頻度を求め、当該頻度が所定の条件を満たす上記各個体の上記ハプロタイプ文字列の上記組合せを、上記ハプロタイプの組合せとして推定するハプロタイプ推定ステップと、を実行させることを特徴とする。

本発明によれば、塩基多型が考慮されていないコピー数多型に関する実験データから、曖昧なコピー数の総和しか得られない場合であっても、ハプロタイプとその頻度を高精度で推定することができる。

また、本発明によれば、コピー数多型に関するデータに１座位当たりの遺伝子型データが加わったデータであっても、ハプロタイプとその頻度を推定することができる。

また、本発明によれば、文字列による取り得るハプロタイプ文字列の組合せの推定結果が、実験データと一致するか否かを検証し、冗長な文字列を排除するので、精度よく取り得るハプロタイプの組合せを算出することができる。

また、本発明によれば、ＥＭ法を用いて、高精度にハプロタイプの尤度を検定することができる。

また、本発明によれば、個体におけるコピー数多型のコピー数の総和データを用いて、ハプロタイプにおけるコピー数とその頻度を高精度で推定することができる。

以下に、本発明にかかるハプロタイプ推定装置、ハプロタイプ推定方法、および、プログラム、並びに、記録媒体の実施の形態を図面に基づいて詳細に説明する。なお、この実施の形態によりこの発明が限定されるものではない。

［本発明の原理］

以下、本発明の原理および概要について説明し、その後、本発明の構成および処理等について詳細に説明する。ここで、図５は、１座位当たりの二つの相同染色体に渡る総コピー数（コピー数の総和）の多型テーブルの一例を表す図である。図５において、Ｌは座位を表す。

図５に示すように、実験結果によっては、各個体（個体１、個体２、・・・）、各座位（Ｌ_１、Ｌ_２、Ｌ_３、・・・）における、二つの相同染色体に渡る総コピー数の表（テーブル）が得られる場合（コピー数多型がある場合）があり、総コピー数は２に限られない。

また、実験によっては、総コピー数がはっきりと確定されず、例えば“０コピーあるいは１コピー”、“２よりは多いコピー数”など、コピー数の状態が唯一つに決まっていない総コピー数（すなわち曖昧な総コピー数）を含むデータが得られる場合がある。ここで、図６は、曖昧なコピー数を含む、１座位当たりの二つの相同染色体に渡る総コピー数のデータの一例を示す図である。

図６において、“０ｏｒ１”は、０コピーあるいは１コピーで表し、“＞２”は、２よりは多いコピー数を表している。実際に、実験で得られるデータは、このような曖昧なデータしか得られない場合があるので、このようなデータから、相同染色体１本上の一座位におけるコピー数と、そのコピー数をアレルと見た時の頻度、あるいはそのようなアレルの複数座位に渡る組合せとしてのハプロタイプとその頻度を推定する必要があった。

また、そのような総コピー数のデータに、それとは異なる座位における１座位当たりの遺伝子型データが組み合わさったデータから、ハプロタイプを推定する必要があった。ここで、図７は、１座位当たりの二つの相同染色体に渡る総コピー数のデータに、１座位当たりの遺伝子型データが組み合わさったデータの一例を示す図である。図７に示すように、座位Ｌ_１およびＬ_４では、１座位当たりの二つの相同染色体に渡る総コピー数を表しており、座位Ｌ_２、Ｌ_３、Ｌ_５では、１座位当たりの遺伝子型データが表現されている。

ここで、「ハプロタイプ」とは、二つの相同染色体に渡る総コピー数（曖昧な総コピー数も含む。）が得られた座位に対しては、相同染色体１本上におけるコピー数をアレルとし、遺伝子型が得られた座位に対してはそこで既に定義されたアレルを使った、複数座位に渡るアレルの組合せのことである。ここで、図８は、本発明におけるハプロタイプの表現例を示す図である。図８において、Ｎ（Ｌ）は、相同染色体１本上の座位Ｌにおけるコピー数Ｎ（整数値）を表しており、Ａ（Ｌ）は座位Ｌに対応するアレルＡを表している。

図７において示されるような、二つの相同染色体に渡る総コピー数のデータは、遺伝子型データではないため、従来技術（１）を適用することができないという問題があった。また、その総コピー数データは塩基多型を考慮したコピー数多型に関する実験データとも異なり、さらに今回取り扱うハプロタイプはコピー単位の組合せのことではない為、従来技術（２）も適用することができないという問題があった。

本発明は、上述した種々の問題に鑑みて本発明者により鋭意検討された結果得られたものであり、一座位当たりの二つの相同染色体に渡る、コピー数の総和（以下、「総コピー数」と呼ぶ。）のデータを処理し、相同染色体１本上の一座位におけるコピー数と、そのコピー数をアレルと見たときの頻度、あるいは、そのようなアレルの複数座位に渡る組合せとしてのハプロタイプとその頻度を推定する。さらには、一座位当たりの二つの相同染色体に渡る総コピー数のデータに、それとは異なる座位における一座位当たりの遺伝子型データが組み合わさったデータから、上述した背景技術のようなハプロタイプとその頻度を推定する。

ここでは、本発明の原理を説明するために一例として、総コピー数と遺伝子型が組み合わさったデータからの推定について説明する。なぜなら、総コピー数と遺伝子型が組み合わさったデータからハプロタイプを推定することは、総コピー数データのみからハプロタイプを推定することを含むより一般的な推定であり、かつ、総コピー数データから複数座位に渡るハプロタイプとその頻度を推定することは、一座位におけるコピー数のアレルとその頻度を推定することを含むより一般的な推定であるからである。

与えられたデータからハプロタイプとその頻度を推定する原理は、各個体に対しデータと矛盾しないディプロタイプ（２つの相同染色体上それぞれにある２つのハプロタイプの組合せ）を見出し、データを利用してハプロタイプの頻度を計算することである。ここで、二つの相同染色体に渡る（確定された）総コピー数がデータとして得られた場合、「矛盾しない」とは、各座位に対しては、ディプロタイプから数えられるその座位のコピー数アレルのコピー数の和が、そのデータにおける総コピー数に一致するということである。例えば、上述の図７のデータの個体１、座位Ｌ_１に対しては、図８におけるハプロタイプ１，ハプロタイプ２からなるディプロタイプを考えたとき、Ｎ_１１（Ｌ_１）＋Ｎ_１２（Ｌ_１）がデータにおける総コピー数２に一致するということである。

また、総コピー数が曖昧な総コピー数である場合、「矛盾しない」とは、それが複数の総コピー数が提示されたもの（例えば、不等号を含まず“ｏｒ”で複数の総コピー数が連結されたもの）であるときは、ディプロタイプから数えられるその座位のコピー数アレルのコピー数の和が、その複数の総コピー数のいずれかに一致するということである。例えば図７のデータの個体１、座位Ｌ_４に対しては、図８におけるハプロタイプ１，ハプロタイプ２から成るディプロタイプを考えたとき、Ｎ_４１（Ｌ_４）＋Ｎ_４２（Ｌ_４）が、データ「０ｏｒ１」が提示する複数の総コピー数Ｎ＝０，１のいずれかに一致するということである。また、曖昧な総コピー数が“ある数より大きい”など不等号で表現される総コピー数であるとき、「矛盾しない」とは、ディプロタイプから数えられるその座位のコピー数アレルのコピー数の和が、不等号で表現される条件を満たすということである。例えば、図７のデータの個体２、座位Ｌ_１に対しては、図８におけるハプロタイプ１，ハプロタイプ２からなるディプロタイプを考えたとき、Ｎ_１１（Ｌ_１）＋Ｎ_１２（Ｌ_１）がデータ「＞２」によって表現された“２より大きい”を満たす、すなわちＮ_１１（Ｌ_１）＋Ｎ_１２（Ｌ_１）＞２を満たす、ということである。

また、遺伝子型がデータとして得られた各座位に対しては、上記「矛盾しない」とは、ディプロタイプから数えられるその座位の各アレルの数が、その遺伝子型データにおける各アレルの数（カウント数）と一致するということである。例えば、図７のデータの個体１、座位Ｌ_２に対しては、図８におけるハプロタイプ１，ハプロタイプ２からなるディプロタイプを考えたとき、ディプロタイプのＬ_２におけるＡ_１，Ａ_２，Ａ_３，…の数それぞれ１，１，０，…（ディプロタイプ表現のＡ_１（Ｌ_２），Ａ_２（Ｌ_２）から得られるカウント数）が、データにおけるＡ_１，Ａ_２，Ａ_３，…の数にそれぞれ一致するということである。

以上、まとめると、「矛盾しない」とは、コピー数多型におけるコピー数の総和（２本の染色体に渡るコピー単位の総数）の条件や、塩基多型におけるカウント数（２本の染色体に渡る多型塩基に特異的なマーカー部位の総数）の条件に反しないことである。

ここでは、本発明の原理を１個体に対して説明したが、データの各座位における総コピー数及び各アレルの数が同じ個体（すなわち、コピー数データやカウント数データのパターンが同じ個体）は同じように扱えるので、それらの数のパターンでデータをまとめてパターン毎にインデックス（ｊ）を付け、集団における特定のパターンを有する個体数Ｎ（ｃ_ｊ）を保持しておいてもよい。以降、この数のパターンを「カウントパターン」と呼ぶ。（各カウントパターンに対し）与えられたデータと矛盾しないディプロタイプが見出せれば、Ｅｘｐｅｃｔａｔｉｏｎ−Ｍａｘｉｍｉｚａｔｉｏｎ（ＥＭ）法やＧｉｂｂｓサンプリング法などによって、データからハプロタイプの頻度が計算できる。

各カウントパターンに対し、データと矛盾しないディプロタイプを見出す方法は種々考えられる。ここではデータと矛盾しない、あらゆる可能なディプロタイプを見出すことを想定して説明する。方法としては、例えば各カウントパターンに対し、まず一座位毎にデータと矛盾しないあらゆる可能な遺伝子型を見出し、それから全座位に渡ってあらゆる可能なディプロタイプを構成する方法が考えられる。ここで、図１０は、本発明の概要を模式的に示したフロー図である。

その第一段階では、二本の相同染色体に渡る（確定された）総コピー数がデータとして得られた座位に対しては、その総コピー数と和が等しい（０を含む正の）整数値２つ一組をあらゆる組合せで作り、２つの整数値をそれぞれコピー数を表すアレルで、２つ一組を遺伝子型とすればよい。例えば、総コピー数がＮのとき、［Ｎ／０］，［Ｎ−１／１］，［Ｎ−２／２］，…，［０／Ｎ］（ここで、［ ］は遺伝子型、あるいはディプロタイプを表し、“／”でアレルまたはハプロタイプを分ける位置を表す。）に対し、重複する冗長な遺伝子型を除いたものが、求める遺伝子型である。総コピー数が曖昧な、“ｏｒ”で連結された総コピー数である場合は、それぞれの総コピー数に対して、これと同じ方法で遺伝子型を列挙すれば求められる。例えば、０ｏｒ１ｏｒ３の場合、総コピー数０に対し上と同じ方法で遺伝子型を列挙し（例えば、［０／０］）、さらに総コピー数１に対し同じ方法で列挙し（例えば、［１／０］）、さらに総コピー数３に対し同じ方法で列挙し（例えば、［３／０］、［２／１］）、あらゆる組み合わせの遺伝子型（例えば、［０／０］、［１／０］、［３／０］、［２／１］）を求める。

総コピー数が曖昧な、不等号で表現された総コピー数であり、しかも“Ｎより多いコピー数”を表す「＞Ｎ」である場合は、アレルのコピー数にも“Ｎより多いコピー数”を表す「＞Ｎ」を用意して｛０，１，…，Ｎ，＞Ｎ｝の各要素をコピー数アレルとし、それら２つのコピー数の和が総コピー数の不等式条件（＞Ｎ、即ちＮより大きい）を満たすようなあらゆる２つ一組（同じコピー数アレルからなる組を含む。）を遺伝子型とすればよい（ここで、２つのコピー数の和を出す際、アレルコピー数＞Ｎには、どんな数あるいは＞Ｎが足されてもその総コピー数は＞Ｎとなり、不等式条件は満たされる）。ここで、図９は、２つのコピー数アレルからなる遺伝子型とそのコピー数の和を一例として示す図である。

例えば、図９のように、｛０，１，…，Ｎ，＞Ｎ｝から２つのコピー数アレルをあらゆる組合せで作って各組を遺伝子型とし、その２つのアレルのコピー数の和が不等式条件＞Ｎを満たす組だけ選択すればそれが求める遺伝子型である。図９において、（ ）内の数字は２つのアレルのコピー数の和を表す。なお、総コピー数が“Ｎより少ないコピー数”を表す＜Ｎである場合は、０ｏｒ１ｏｒ２…ｏｒＮ−１である場合と同じである。

遺伝子型がデータとして得られた座位に対しては、データにおける数の分だけアレルを表す文字（多型識別文字）を取り、これを２つに分けて組を作り、遺伝子型とすればよい。遺伝子型データの座位に関しては、一座位当たりのアレルの総数は必ず２であるので簡単にアレルを２つに分けられ、遺伝子型は唯一つ得られる。

第二段階においては、各座位の遺伝子型から、全座位に渡るあらゆる可能なディプロタイプを作る。これには、全座位に渡って一つずつ遺伝子型をあらゆる組合せで取り、その各組合せに対し、各座位の各遺伝子型が持つ２つのアレルから、全座位に渡るアレルの組合せとしてのハプロタイプ２つの組合せから構成されるディプロタイプを、２^Ｍ−１個（ここでＭは全座位数）のあらゆる組合せで作り、このディプロタイプの作成を遺伝子型の組合せ全部に対し行った後、最後に重複する冗長なディプロタイプを除けばよい。

例えば、座位Ｌ_１の遺伝子型が［Ｂ_１１／Ｂ_１２］，［Ｂ_１３／Ｂ_１４］、座位Ｌ_２，Ｌ_３の遺伝子型がそれぞれ［Ｂ_２１／Ｂ_２２］，［Ｂ_３１／Ｂ_３２］であるとき（ここでＢは、コピー数アレルあるいは遺伝子型データで定義されているアレルを表す。）、まず全座位に渡る遺伝子型の組合せ、［Ｂ_１１／Ｂ_１２］，［Ｂ_２１／Ｂ_２２］，［Ｂ_３１／Ｂ_３２］と［Ｂ_１３／Ｂ_１４］，［Ｂ_２１／Ｂ_２２］，［Ｂ_３１／Ｂ_３２］を取る。次に、最初の組合せ［Ｂ_１１／Ｂ_１２］，［Ｂ_２１／Ｂ_２２］，［Ｂ_３１／Ｂ_３２］に対しては、［Ｂ_１１Ｂ_２１Ｂ_３１／Ｂ_１２Ｂ_２２Ｂ_３２］，［Ｂ_１１Ｂ_２１Ｂ_３２／Ｂ_１２Ｂ_２２Ｂ_３１］，［Ｂ_１１Ｂ_２２Ｂ_３１／Ｂ_１２Ｂ_２１Ｂ_３２］，［Ｂ_１１Ｂ_２２Ｂ_３２／Ｂ_１２Ｂ_２１Ｂ_３１］の２^３−１個のあらゆる組合せを列挙する。さらに、［Ｂ_１３／Ｂ_１４］，［Ｂ_２１／Ｂ_２２］，［Ｂ_３１／Ｂ_３２］に対しても、同様に２^３−１個のあらゆる組合せを列挙する。

最後に重複する冗長なディプロタイプを除く。図１０に、あるカウントパターンｃ_１に対し、データと矛盾しないディプロタイプを見出す方法の例を示した。その図において、Ｌ_１，Ｌ_３は二つの相同染色体に渡る総コピー数がデータとして得られた座位であり、Ｌ_２は遺伝子型がデータとして得られた座位である。図１０＜ＳＡ−２＞の記号（○の中に×）は、全座位に渡って一つずつ遺伝子型をあらゆる組合せで取ることを示す。

こうして得られたディプロタイプを用い、さらに与えられたデータを使って、ハプロタイプの頻度を求める。この方法も幾つか考えられるが、一例として、Ｅｘｐｅｃｔａｔｉｏｎ−Ｍａｘｉｍｉｚａｔｉｏｎ（ＥＭ）法を用いてもよい。これは、得られたディプロタイプに対し、その存在の重みを割り付け、そのディプロタイプが含むハプロタイプの個数を、その重み分を考慮して数え、ハプロタイプの頻度を計算し（Ｍステップ）、次にそのハプロタイプ頻度からハーディ・ワインバーグの法則を使って、ディプロタイプの存在の重みを更新し（Ｅステップ）、さらにその更新された重みから、Ｍステップ、次にＥステップ、さらにＭステップ、…と手続きを繰り返して頻度を更新していく方法である。例えば、以下の数式１のようなＭステップと、数式２−１の様なＥステップを交互に行って、ハプロタイプの頻度の更新していく。

上記数式２−１において、Ｐは頻度、ｎは個体の総数、ｉ，ｊ，ｋはハプロタイプ、カウントパターン、カウントパターン内でのディプロタイプのインデックス、ｈ，ｄはハプロタイプ、ディプロタイプ，Ｎ（ｃ_ｊ）はカウントパターンｃ_ｊを持つ個体の数、ｗはＥＭ法における、カウントパターン内でのディプロタイプの重みを表す。数２−１におけるディプロタイプの頻度は、ハーディ・ワインバーグの法則から計算する。ハーディ・ワインバーグの法則とは遺伝学における自然法則であり、この法則によって、ディプロタイプを構成する２つのハプロタイプとその確率（あるいは頻度）が分かったとき、そのディプロタイプの確率（あるいは頻度）が計算できる。この法則は、例えば、以下の数２−２の様に表現される。以上のような方法で、ハプロタイプとその頻度を推定する。ここで、頻度の非常に低いハプロタイプは存在しないと解釈してもよい。

［本発明の概要］
まず、本発明は、集団における各個体の少なくともコピー数多型を含む実験データからハプロタイプを推定する、制御部と記憶部を少なくとも備えたハプロタイプ推定装置において実行される。

そして、図１０に示すように、本ハプロタイプ推定装置は、個体毎に、実験データから得られた、標識によって特定されるマーカー部位に対応付けられたコピー数多型の総コピー数Ｎに対する一または複数の条件式を、マーカー部位の種類毎に記憶する多型テーブルを備える（ＳＡ−１）。ここで、コピー数Ｎに対する条件式は、コピー数の総和Ｎの条件を規定する一または複数の等式または不等式を含んでもよく、例えば、「Ｎ＝２」や「Ｎ＞０」や、「Ｎ＝１ｏｒ２」や「Ｎ≧２」等であってもよい。また、当該多型テーブルは、個体毎に、塩基多型を更に含む実験データから得られた、標識によって特定されるマーカー部位に対応付けられた多型塩基をカウントしたカウント数を、当該多型塩基の種類毎に更に格納されてもよい。

そして、本ハプロタイプ推定装置は、コピー数の総和Ｎについて、足し合わせた場合に、多型テーブルに記憶された総コピー数Ｎの条件式を満たす任意の２つの整数に分割する（ＳＡ−２）。例えば、本ハプロタイプ推定装置は、図５＜ＳＡ−２＞に示すように、総コピー数Ｎの条件式Ｎ＝２の場合に、［２／０］または［１／１］に分割する。また、例えば、総コピー数Ｎの条件式Ｎ≦２の場合は、［２／０］、［１／１］、［１／０］または［０／０］に分割する。

そして、本ハプロタイプ推定装置は、分割された２つの整数をそれぞれ多型識別文字（例えば、「２」や「３」など）として表し、個体において（複数の座位に渡って）多型識別文字を列挙したハプロタイプ文字列の組合せ（ディプロタイプを表す。）として格納する（ＳＡ−３〜４）。ここで、図１０＜ＳＡ−３＞に示すように、コピー数多型のコピー数を表す多型識別文字は、一例として、アラビア数字により表してもよい。ここで、図１０＜ＳＡ−４＞に示すように、多型テーブルに個体毎の標識によって特定されるマーカー部位に対応付けられた多型塩基の種類毎の多型塩基のデータがある場合は、ハプロタイプ文字列の組合せにおいて、多型塩基の種類に対応付けた多型識別文字（例えば、「Ａ_１」や「Ａ_２」など）を更に列挙してもよい（例えば、「２Ａ_１３／０Ａ_２０」）。ここで、本ハプロタイプ推定装置は、まず遺伝子型の組合せを求めてから、次に全ての可能なハプロタイプの組合せ（ディプロタイプ）を求めてもよい。

すなわち、本ハプロタイプ推定装置は、一例として、図１０に示すように、個体毎および座位（マーカー部位）毎の実験データ＜ＳＡ−１＞から、各座位における全ての可能な組合せを求め（ＳＡ−２）、次に、複数の座位に渡る全ての可能な遺伝子型の組合せを求め（ＳＡ−３）、最後に、遺伝子型の組合せから複数の座位に渡る全ての可能なハプロタイプの組合せ（ディプロタイプ）を表す文字列を求めてもよい（ＳＡ−４）。また、ここで、本ハプロタイプ推定装置は、作成したハプロタイプ文字列の組合せにおける多型識別文字が表す整数の和が、多型テーブルにおける総コピー数Ｎの条件式を満たすか否か（または多型塩基のカウント数と一致するか否か）を確認し、条件を満たさない場合に当該ハプロタイプ文字列を除外してもよい。また、ハプロタイプ推定装置は、作成したハプロタイプ文字列の組合せが他のハプロタイプ文字列の組合せと重複するか否か判断し、重複する場合に当該ハプロタイプ文字列の組合せを除外してもよい（ＳＡ−５）。

そして、本ハプロタイプ推定装置は、集団において、同一であるハプロタイプ文字列の数を集計し、ハプロタイプ文字列の集団における頻度を求め、当該頻度が所定の条件を満たす各個体のハプロタイプ文字列の組合せを、ハプロタイプの組合せとして推定する。ここで、本ハプロタイプ推定装置は、頻度をハーディ・ワインバーグ（Ｈａｒｄｙ−Ｗｅｉｎｂｅｒｇ）の法則に基づいて算出し、集団における頻度がハーディ・ワインバーグ平衡となる場合に所定の条件を満たすと判定してもよい。

また、本ハプロタイプ推定装置は、ハプロタイプの推定処理において、ＥＭ（Ｅｘｐｅｃｔａｔｉｏｎ−Ｍａｘｉｍｉｚａｔｉｏｎ）法を用いて、集団におけるハプロタイプ文字列の頻度を、当該ハプロタイプ文字列を少なくとも一方に有する組合せの頻度により重み付けして算出するＭステップと、組合せの頻度を、当該組合せを構成するハプロタイプ文字列の頻度の積により求め、当該組合せの頻度に基づいて重みを算出するＥステップと、を頻度の値が収束するまで交互に繰り返すことにより所定の条件を満たす頻度を算出してもよい。ＥＭ法の詳細な処理については、後述する。以上で、本発明の概要の説明を終える。

［ハプロタイプ推定装置の構成］
まず、本ハプロタイプ推定装置の構成について説明する。図１１は、本発明が適用される本ハプロタイプ推定装置の構成の一例を示すブロック図であり、該構成のうち本発明に関係する部分のみを概念的に示している。

図１１において、ハプロタイプ推定装置１００は、概略的に、ハプロタイプ推定装置１００の全体を統括的に制御するＣＰＵ等の制御部１０２、通信回線等に接続されるルータ等の通信装置（図示せず）に接続される通信制御インターフェース部１０４、入力部１１２や出力部１１４に接続される入出力制御インターフェース部１０８、および、各種のデータベースやテーブルなどを格納する記憶部１０６を備えて構成されており、これら各部は任意の通信路を介して通信可能に接続されている。

記憶部１０６に格納される各種のデータベースやテーブル（多型テーブル１０６ａ〜実験データファイル１０６ｃ）は、固定ディスク装置等のストレージ手段であり、各種処理に用いる各種のプログラムやテーブルやファイルやデータベース等を格納する。

これら記憶部１０６の各構成要素のうち、多型テーブル１０６ａは、個体毎に、実験データから得られた、標識によって特定されるマーカー部位に対応付けられたコピー数多型の総コピー数Ｎ（二本の相同染色体上のコピー単位の数の総和）に対する一または複数の条件式をマーカー部位の種類毎（すなわちコピー多型の種類毎）にコピー数データとして記憶する多型テーブルである。ここで、多型テーブル１０６ａは、個体毎に、塩基多型を更に含む実験データから得られた、標識によって特定されるマーカー部位に対応付けられた多型塩基をカウントしたカウント数を当該多型塩基の種類毎にカウント数データとして更に記憶してもよい。上述したように、図５や図６や図７や図１０＜ＳＡ−１＞は、多型テーブル１０６ａに格納されるコピー数データおよび／またはカウント数データの一例を示す。この多型テーブル１０６ａに格納される情報は、一例として図７等に示すように、個体毎、座位毎（コピー数多型の種類に対応したマーカー部位毎、および／または、各座位における多型塩基の種類毎）のコピー数（条件式を含んでもよい。）／カウント数を定義している。ここで、コピー数Ｎに対する条件式は、コピー数の総和Ｎの条件を規定する一または複数の等式または不等式を含んでもよく、例えば、「Ｎ＝３」や「Ｎ＞０」、「Ｎ＝１ｏｒ２」、「Ｎ≠２」、「０＜Ｎ≦２」等であってもよい。

また、ハプロタイプ文字列ファイル１０６ｂは、多型テーブル１０６ａに記憶された個体のコピー数データおよび／またはカウント数データに基づいて算出された、取り得るハプロタイプ文字列の組合せを記憶するハプロタイプ文字列記憶手段である。ハプロタイプ文字列ファイル１０６ｂは、一例として、ハプロタイプ文字列を、図８に示したようなハプロタイプ表現形式で記憶してもよい。

また、実験データファイル１０６ｃは、集団における各個体の少なくともコピー数多型を含む実験データを記憶する実験データ記憶手段である。ここで、実験データは、一塩基多型等の塩基多型のデータ（遺伝子型データ等）を含んでもよい。一例として、コピー数多型のコピー単位上の、標識によって特定されるマーカー部位に対応付けられた多型塩基を示す実験データ（ＤＮＡチップやＰＣＲ等による実験データなど）を記憶してもよい。また、標識としては、蛍光色素プローブの他、蛍光特性を持たない色素や、放射性同位体、ＧＦＰ・ＧＲＰなどのタンパク質、Ｈｉｓタグ、ビオチン化などによって識別可能なプローブ等を用いてもよい。

また、図１１において、通信制御インターフェース部１０４は、ハプロタイプ推定装置１００とネットワーク３００（またはルータ等の通信装置）との間における通信制御を行う。すなわち、通信制御インターフェース部１０４は、他の端末と通信回線を介してデータを通信する機能を有する。

また、図１１において、入出力制御インターフェース部１０８は、入力部１１２や出力部１１４の制御を行う。ここで、出力部１１４としては、モニタ（家庭用テレビを含む。）の他、スピーカ等を用いることができる。また、入力部１１２としては、キーボード、マウス、およびマイク等を用いることができる。

また、図１１において、制御部１０２は、ＯＳ（Ｏｐｅｒａｔｉｎｇ Ｓｙｓｔｅｍ）等の制御プログラム、各種の処理手順等を規定したプログラム、および所要データを格納するための内部メモリを有し、これらのプログラム等により、種々の処理を実行するための情報処理を行う。制御部１０２は、機能概念的に、コピー数総和分割部１０２ａ、ハプロタイプ文字列格納部１０２ｂ、ハプロタイプ推定部１０２ｃ、多型テーブル作成部１０２ｅを備えて構成されている。なお、理解の容易のために上述した記号を用いて説明することがある。

このうち、コピー数総和分割部１０２ａは、コピー数の総和Ｎについて、足し合わせた場合に、多型テーブル１０６ａに記憶された総コピー数Ｎの条件式を満たす任意の２つの整数に分割するコピー数総和分割手段である。ここで、コピー数総和分割部１０２ａは、総コピー数Ｎの条件式が曖昧な数を規定している場合（例えば、「Ｎ＝２」ではなく「Ｎ＞２」など）には、２つの整数の少なくとも一方を曖昧な整数（例えば、［＞２／０］）に分割してもよい。

また、ハプロタイプ文字列格納部１０２ｂは、コピー数総和分割部１０２ａにより分割された２つの整数をそれぞれ多型識別文字（例えば、アラビア数字の「０」や「１」）として表し、個体において多型識別文字を列挙したハプロタイプ文字列の組合せ（例えば、「３０２／０２１」）としてハプロタイプ文字列ファイル１０６ｂに格納するハプロタイプ文字列格納手段である。ここで、ハプロタイプ文字列格納部１０２ｂは、多型テーブル１０６ａに、個体毎に、塩基多型を更に含む実験データ（遺伝子型データ）から得られた、標識によって特定されるマーカー部位に対応付けられた多型塩基が、当該多型塩基の種類毎に格納されている場合（すなわち、カウント数データが格納されている場合）に、ハプロタイプ文字列の組合せにおいて、多型塩基の種類に対応付けた多型識別文字を更に列挙してもよい（例えば、「３０２Ａ_１／０２１Ａ_２」）。ここで、ハプロタイプ文字列格納部１０２ｂは、作成したハプロタイプ文字列の組合せにおける対応する多型識別文字が表す整数の和（二本の相同染色体上のコピー多型のコピー数の和を表している。）が、多型テーブル１０６ａにおける総コピー数Ｎの条件式を満たすか否か確認し、条件式を満たさない場合に当該ハプロタイプ文字列を除外してもよい。また、ハプロタイプ文字列格納部１０２ｂは、作成したハプロタイプ文字列の組合せにおける多型塩基を表す多型識別文字の数が、多型テーブル１０６ａに格納されたカウント数データと一致するか否か確認し、一致しない場合に当該ハプロタイプ文字列を除外してもよい。また、ハプロタイプ文字列格納部１０２ｂは、作成したハプロタイプ文字列の組合せが他のハプロタイプ文字列と重複するか否か判断し、重複する場合に当該ハプロタイプ文字列を除外してもよい。

また、ハプロタイプ推定部１０２ｃは、ハプロタイプ文字列ファイル１０６ｂを参照して、集団において、同一であるハプロタイプ文字列の数を集計し、当該ハプロタイプ文字列の集団における頻度を求め、当該頻度が所定の条件を満たす各個体のハプロタイプ文字列の組合せを、ハプロタイプの組合せとして推定するハプロタイプ推定手段である。ここで、ハプロタイプ推定部１０２ｃは、ハプロタイプ文字列の頻度を、ハーディ・ワインバーグの法則に基づいて算出し、所定の条件を、集団におけるハーディ・ワインバーグ平衡としてもよい。

ここで、ハプロタイプ推定部１０２ｃは、図１１に示すように、ハプロタイプ頻度算出部１０２ｄを備えて構成される。ハプロタイプ頻度算出部１０２ｄは、ＥＭ（Ｅｘｐｅｃｔａｔｉｏｎ−Ｍａｘｉｍｉｚａｔｉｏｎ）法を用いて、集団におけるハプロタイプ文字列の頻度を、当該ハプロタイプ文字列を少なくとも一方に有する組合せの頻度により重み付けして算出するＭステップと、ハプロタイプの組合せの頻度を、当該組合せを構成する２つのハプロタイプ文字列の頻度の積により求め、当該組合せの頻度に基づいて重みを算出するＥステップと、を頻度の値が収束するまで交互に繰り返すハプロタイプ頻度算出手段である。ここで、ハプロタイプ頻度算出部１０２ｄは、Ｍステップにおいて算出されたハプロタイプ文字列の頻度と、前回のＭステップにおいて算出されたハプロタイプ文字列の頻度と、の対数尤度差を求め、対数尤度差が所定の閾値以下となった場合に、頻度の値が収束したと判定してもよい。これにより、ＥＭ法の収束条件を適切に設定することができるので、精度を保証しながら計算時間を節約することができる。

また、ここで、ハプロタイプ頻度算出部１０２ｄは、Ｍステップにおいて、下記の数式１に基づいて、ハプロタイプ文字列の頻度Ｐ（ｈ_ｉ）を算出し、Ｅステップにおいて、下記の数式２に基づいてハプロタイプ文字列の組合せの頻度Ｐ（ｄ_ｊｋ）を求め、重みｗ_ｊｋとして当該ハプロタイプ文字列の組合せの頻度を集団における組合せの頻度の総和で除して算出してもよい。

ここで、Ｐ（ｈ_ｉ）はハプロタイプ文字列の頻度を表し、ｈはハプロタイプ文字列を表し、ｉはハプロタイプ文字列のインデックスを表す。また、ｎは集団を構成する個体の数、ｊは、多型テーブル１０６ａにおけるコピー数データおよび／またはカウント数データのパターン（以下「カウントパターン」と呼ぶ。）のインデックス、ｋはハプロタイプ文字列の組合せのインデックスである。また、Ｎ（ｃ_ｊ）は、カウントパターンｊ（「ｃ_ｊ」と表記する。）を持つ個体の数を表す。また、δ（ｈ_ｉ，ｄ_ｊｋ）は、ハプロタイプ文字列の組合せｄ_ｊｋが一方に当該ハプロタイプ文字列ｈ_ｉを有する場合に１を返し、両方に当該ハプロタイプ文字列ｈ_ｉを有する場合に２を返し、当該ハプロタイプ文字列ｈ_ｉを持たない場合に０を返す関数であり、ｄはハプロタイプ文字列の組合せを表す。また、ｗ_ｊｋはハプロタイプ文字列の組合せの頻度による重み（当該ハプロタイプ文字列の組合せｄ_ｊｋの頻度を、組合せｄ_ｊｋの頻度の総和で除した数）である。

ここで、Ｐ（ｄ_ｊｋ）は、ハプロタイプ文字列の組合せの頻度を表す。また、ｈ_ｌおよびｈ_ｍは当該組合せを構成する２つのハプロタイプ文字列を表し、Ｐ（ｈ_ｌ）およびＰ（ｈ_ｍ）は、当該２つのハプロタイプ文字列の頻度をそれぞれ表す。

また、多型テーブル作成部１０２ｅは、実験データファイル１０６ｃに記憶された個体毎の実験データを用いて、標識によって特定されるマーカー部位に対応付けられた総コピー数Ｎ（二本の染色体に渡るコピー単位の数の総和）をコピー数多型の種類毎にコピー数データとして、および／または、標識によって特定されるマーカー部位に対応付けられた多型塩基をカウントしたカウント数を多型塩基の種類毎にカウント数データとして、多型テーブル１０６ａに格納する多型テーブル作成手段である。ここで、格納されるカウント数およびコピー数は、一意に特定される数に限られず、曖昧な数（例えば、「＞５」、「１ｏｒ２」、「≠０」）であってもよい。

ここで、本ハプロタイプ推定装置１００は、ルータ等の通信装置および専用線等の有線または無線の通信回線を介して、ネットワーク３００に通信可能に接続されてもよい。この場合、本システムは、概略的にハプロタイプ推定装置１００と、コピー数データやカウント数データ等に関する外部データベースやハプロタイプ推定プログラム等の外部プログラム等を提供する外部システム２００とを、ネットワーク３００を介して通信可能に接続して構成される。ここで、図１１において、ネットワーク３００は、ハプロタイプ推定装置１００と外部システム２００とを相互に接続する機能を有し、例えば、インターネット等である。

ここで、外部システム２００は、ネットワーク３００を介して、ハプロタイプ推定装置１００と相互に接続され、利用者に対してコピー数データやカウント数データ等に関する外部データベースやハプロタイプ推定プログラム等の外部プログラム等を実行するウェブサイトを提供する機能を有する。ここで、外部システム２００は、ＷＥＢサーバやＡＳＰサーバ等として構成していてもよく、そのハードウェア構成は、一般に市販されるワークステーション、パーソナルコンピュータ等の情報処理装置およびその付属装置により構成していてもよい。また、外部システム２００の各機能は、外部システム２００のハードウェア構成中のＣＰＵ、ディスク装置、メモリ装置、入力装置、出力装置、通信制御装置等およびそれらを制御するプログラム等により実現される。以上で、本ハプロタイプ推定装置１００の構成の説明を終える。

［本ハプロタイプ推定装置１００の処理］
次に、このように構成された実施の形態における本ハプロタイプ推定装置１００の処理の一例について、以下に図１２〜図１４を参照して詳細に説明する。

［ハプロタイプ推定処理］
本実施の形態におけるハプロタイプ推定処理の一例について、以下に図１２を参照して説明する。ここで、図１２は、本ハプロタイプ推定装置１００のハプロタイプ推定処理の一例を示す図である。

図１２に示すように、コピー数総和分割部１０２ａは、コピー数の総和Ｎについて、足し合わせた場合に、多型テーブルに記憶された総コピー数Ｎの条件式を満たす任意の２つの整数に分割する（ＳＢ−１）。ここで、コピー数総和分割部１０２ａは、総コピー数Ｎの条件式が曖昧な数を規定している場合（例えば、「Ｎ＞２」）には、２つの整数を曖昧な整数（例えば、［＞２／０］）で分割してもよい。

そして、ハプロタイプ文字列格納部１０２ｂは、コピー数総和分割部１０２ａにより分割された２つの整数をそれぞれ多型識別文字として表し、個体における多型識別文字を複数の座位に渡って列挙したハプロタイプ文字列の組合せを、とりうる全ての組み合せで算出してハプロタイプ文字列ファイル１０６ｂに格納する（ＳＢ−２）。なお、多型テーブル１０６ａに、遺伝子型データから得られた、標識によって特定されるマーカー部位に対応付けられた多型塩基のカウント数が、多型塩基の種類毎に格納されている場合には、ハプロタイプ文字列格納部１０２ｂは、多型塩基の種類に対応付けた多型識別文字を加えて、ハプロタイプ文字列の組合せに列挙してもよい。ここで、ハプロタイプ文字列格納部１０２ｂは、作成したハプロタイプ文字列の組合せにおける対応する多型識別文字が表す整数の和（二本の相同染色体上のコピー多型のコピー数の和を表している。）が、多型テーブル１０６ａにおける総コピー数Ｎの条件式を満たすか否か確認し、条件式を満たさない場合に当該ハプロタイプ文字列を除外してもよい。また、ハプロタイプ文字列格納部１０２ｂは、作成したハプロタイプ文字列の組合せが他のハプロタイプ文字列と重複するか否か判断し、重複する場合に当該ハプロタイプ文字列を除外してもよい。

そして、ハプロタイプ推定部１０２ｃは、ハプロタイプ文字列格納部１０２ｂにより個体毎に格納されたハプロタイプ文字列の組合せを記憶するハプロタイプ文字列ファイルハプロタイプ文字列ファイル１０６ｂを参照して、集団において、同一であるハプロタイプ文字列の数を集計する（ＳＢ−３）。

そして、ハプロタイプ推定部１０２ｃは、ハプロタイプ文字列の集団における頻度を計算する（ＳＢ−４）。ここで、ハプロタイプ推定部１０２ｃは、ハプロタイプ文字列の頻度をハーディ・ワインバーグの法則により算出してもよい。

そして、ハプロタイプ推定部１０２ｃは、ハプロタイプ文字列の頻度が所定の条件を満たす各個体のハプロタイプ文字列の組合せを抽出し、ハプロタイプの組合せとして推定する（ＳＢ−５）。ここで、ハプロタイプ推定部１０２ｃは、所定の条件として、集団におけるハーディ・ワインバーグ平衡を設定してもよい。以上で、本実施の形態におけるハプロタイプ推定処理を終える。

［ＥＭ法による処理］
ハプロタイプ推定部１０２ｃの処理によるＥＭ法の詳細な処理の一例について説明する。すなわち、ハプロタイプ推定部１０２ｃは、ハプロタイプ頻度算出部１０２ｄの処理により、集団におけるハプロタイプ文字列の頻度を、ＥＭ（Ｅｘｐｅｃｔａｔｉｏｎ−Ｍａｘｉｍｉｚａｔｉｏｎ）法を用いて、効率よく計算する。

ここで、ＥＭ法とは、得られたディプロタイプ（本実施の形態においては、ハプロタイプ文字列の組合せとして表す。）に対し、その存在の重みを割り付け、そのディプロタイプが含むハプロタイプの個数を、重み分を考慮して数え、ハプロタイプの頻度を計算し（Ｍステップ）、次にそのハプロタイプ頻度からハーディ・ワインバーグの法則を使って、ディプロタイプの存在の重みを更新し（Ｅステップ）、さらにその更新された重みから、Ｍステップ、次にＥステップ、さらにＭステップ、・・・と処理を繰り返して、頻度を更新していく方法である。例えば、下記の数式１に基づくＭステップと、下記の数式２（数式２−１および数式２−２）に基づくＥステップを交互に行うことにより、ハプロタイプの頻度の更新していく。

また、Ｐ（ｈ_ｉ）はハプロタイプ文字列の頻度を表し、ｈはハプロタイプ文字列を表し、ｉはハプロタイプ文字列のインデックスを表す。また、ｎは集団を構成する個体の数、ｊはカウントパターン（多型テーブル１０６ａにおけるコピー数データおよび／またはカウント数データのパターン）のインデックス、ｋはハプロタイプ文字の組合せのインデックス、Ｎ（ｃ_ｊ）はカウントパターン（ｃ_ｊ）を持つ個体の数を表す。また、δ（ｈ_ｉ，ｄ_ｊｋ）は、ハプロタイプ文字列の組合せｄ_ｊｋが一方に当該ハプロタイプ文字列ｈ_ｉを有する場合に１を返し、両方に当該ハプロタイプ文字列ｈ_ｉを有する場合に２を返し、当該ハプロタイプ文字列ｈ_ｉを持たない場合に０を返す関数であり、ｄはハプロタイプ文字列の組合せを表す。また、ｗ_ｊｋは下記の数式２−１に基づくハプロタイプ文字列の組合せの頻度による重みである。

ここでＰ（ｄ_ｊｋ）は、ハーディ・ワインバーグの法則を表す下記の数式２−２に基づいて計算する。

ここで、上記の数式２−２は、ハーディ・ワインバーグの法則を示している。ハーディ・ワインバーグの法則とは遺伝学における自然法則であり、この法則によって、ディプロタイプを構成する２つのハプロタイプとその確率（あるいは頻度）が分かった時、そのディプロタイプの確率（あるいは頻度）が計算できる。なお、上記ＥＭ法においては、頻度の非常に低いハプロタイプは存在しないと解釈される。このように、ハプロタイプ推定部１０２ｃは、ハプロタイプ頻度算出部１０２ｄの処理により、一例として上記のＥＭ法を用いて、ハプロタイプの組合せとその頻度を推定する。ここで、ハプロタイプ頻度算出部１０２ｄは、下記の数式に基づいて、Ｍステップにおいて算出されたハプロタイプ文字列の頻度と、前回のＭステップにおいて算出されたハプロタイプ文字列の頻度と、の対数尤度差を求め、対数尤度差が所定の閾値以下となった場合に、頻度の値が収束したと判定してもよい。

以上で、ＥＭ法による処理の一例の説明を終える。

［実施例］
本実施の形態を、プログラミング言語Ｐｅｒｌで実装した実施例について、以下に図１３〜図１４を参照して説明する。図１３は、本実装の枠組みを示すフローチャートである。図１４は、コピー数総和分割処理およびハプロタイプ文字列格納処理の一例を示すフローチャートである。なお、以下の説明において、説明の簡単のために、本実施の形態におけるハプロタイプ文字列を単に「ハプロタイプ」と、ハプロタイプ文字列の組合せを「ディプロタイプ」と述べる場合がある。

図１３に示すように、まず、ハプロタイプ推定装置１００は、多型テーブル１０６ａを参照してデータ（例えば、図５、図６、図７等に示した多型テーブルのデータ）を読み込む（ＳＣ−１）。ここで、ハプロタイプ推定装置は、ＥＭ法において収束判定に使われる数値も読み込んでもよい。

次に、ハプロタイプ推定装置１００は、コピー数総和分割処理およびハプロタイプ文字列格納処理として、データと矛盾しないディプロタイプを全て算出する（ＳＣ−２）。具体的には、以下のＳＣ−２１〜ＳＣ−２５を行う。

すなわち、取り得るハプロタイプの組合せを過不足なく求めるため（多型テーブル１０６ａのコピー数データおよび／またはカウント数データと矛盾しないディプロタイプを見出すため）、ハプロタイプ推定装置は、ディプロタイプ（ハプロタイプの組合せ）を構成する。ここで、データと矛盾しないディプロタイプとは、作成したディプロタイプ（２本の相同染色体に渡る文字の集合）において、コピー数多型におけるコピー数の総和（２本の染色体に渡るコピー単位の総数）の条件や、塩基多型におけるカウント数（２本の染色体に渡る多型塩基に特異的なマーカー部位の総数）の条件に反しないことである。具体的には、多型識別文字が表す２本の染色体上のコピー数を足し合わせたとき、多型テーブル１０６ａに記憶されたコピー数データの条件に反しない（各コピー多型の総コピー数Ｎの条件を満たす）ことである。また、多型識別文字が表す２本の染色体上の多型塩基をカウントしたとき、多型テーブル１０６ａに記憶されたカウント数データのカウント数と一致するということである。ここで、図１４は、コピー数データおよび／またはカウント数データの全パターンにおいて、そのカウントパターンと矛盾しないディプロタイプを作成する処理の一例を示すフローチャートである。

図１４に示すように、まず、ハプロタイプ推定装置１００は、多型テーブル１０６ａに記憶されたコピー数データおよび／またはカウント数データから全カウントパターンを算出する（読み出す）（ＳＣ−２１）。

つぎに、ハプロタイプ推定装置１００は、カウントパターン（ｃ_ｉ）のイテレーション（ｉイテレーション）に入る。ここでｉイテレーションは、最初ｉを１に初期化し、１イテレーション毎にｉを１ずつ増加させ、“ｉ＜＝カウントパターンの個数”である限り繰り返すこととする。

ｉイテレーション内において、ハプロタイプ推定装置１００は、コピー数総和分割部１０２ａおよびハプロタイプ文字列格納部１０２ｂの処理により、各座位に渡って遺伝子型を算出する（ＳＣ−２２）。コピー数総和分割部１０２ａは、二つの相同染色体に渡る（確定された）総コピー数の座位に対しては、その総コピー数と和が等しい（０を含む正の）整数値２つ一組をあらゆる組合せで作る。ハプロタイプ文字列格納部１０２ｂは、その２つの整数値を、アレルを表す文字（多型識別文字）として扱う。多型テーブル１０６ａに格納された総コピー数Ｎが曖昧な、“ｏｒ”で連結された総コピー数である場合は、コピー数総和分割部１０２ａは、それぞれの総コピー数に対して、上記と同じ方法でそれぞれ２つ一組を作る。多型テーブル１０６ａに格納された総コピー数が曖昧な、整数Ｎより少ないコピー数である事を表す“＜Ｎ”である場合は、コピー数総和分割部１０２ａは、“＜Ｎ”を“０ｏｒ１ｏｒ２…ｏｒＮ−１”に変換して、それぞれについて上と同じ方法で２つ一組を作る。多型テーブル１０６ａに格納された総コピー数が曖昧な、Ｎより多いコピー数である事を表す“＞Ｎ”である場合は、｛０，１，…，Ｎ，＞Ｎ｝の各要素をアレルを表す文字（多型識別文字）として、図９に示すように、あらゆる２つ一組を列挙してその２つの和を出し、それが総コピー数の不等式条件（＞Ｎ、即ちＮより大きい）を満たす組だけ選択する。ただし和を出す際、アレルコピー数“＞Ｎ”には、どんな数あるいは“＞Ｎ”が足されてもその総コピー数は“＞Ｎ”となり、不等式条件は満たされる。遺伝子型の座位に対しては、データにおける数の分だけアレルを表す文字（多型識別文字）を取り、２つ一組を作る（その他、適宜、特願２００７−２３７１３９号を参照）。

そして、ハプロタイプ推定装置１００は、各座位に渡って１つずつ遺伝子型を総当たりで取って、全座位に渡る遺伝子型のあらゆる組合せを作る（ＳＣ−２３）。

そして、ハプロタイプ推定装置１００は、作成した遺伝子型の組合せに関するイテレーション（ｊイテレーション）に入る。ここで、ｊイテレーションは、最初ｊを１に初期化し、１イテレーション毎に１ずつ増加させ、“ｊ＜＝遺伝子の組合せの個数”である限り繰り返すこととする。

ｊイテレーション内において、ハプロタイプ文字列格納部１０２ｂは、遺伝子型の組合せｊに対し、全座位に渡って各遺伝子型からアレルを１つずつ取って、座位の順に多型識別文字を並べた文字列を作り、その文字列をハプロタイプ文字列とする。このとき、ハプロタイプ文字列格納部１０２ｂは、各遺伝子型からアレルを１つずつ取った際残った方のアレルからも座位の順にアレル文字を並べた文字列を作り、その文字列をペア（組合せ）となるもう一つのハプロタイプ文字列として、ディプロタイプ（ハプロタイプ文字列の組合せ）を作る（ＳＣ−２４）。この、全座位に渡って各遺伝子型からアレルを一つずつ取ることは総当たりで行って、あらゆる可能な２^Ｍ−１個（ここでＭは全座位数）のディプロタイプを得る。ハプロタイプ推定装置１００は、このディプロタイプの作成を遺伝子型の組合せ全部に対して行い、ｊイテレーションを終了する。

そして、ハプロタイプ文字列格納部１０２ｂは、得られたハプロタイプ文字列の組合せ（ディプロタイプ）から、重複する冗長なディプロタイプを除外する（ＳＣ−２５）。すなわち、ハプロタイプ文字列格納部１０２ｂは、作成したハプロタイプ文字列の組合せが他のハプロタイプ文字列の組合せと重複するか否か判断し、重複する場合に当該ハプロタイプ文字列の組合せを除外する。ここで、ハプロタイプ文字列格納部１０２ｂは、作成したハプロタイプ文字列の組合せにおいて、対応する多型識別文字が表す整数の和が、多型テーブル１０６ａにおけるコピー数の総和Ｎの条件式を満たすか否か確認し、条件式を満たさない場合に当該ハプロタイプ文字列を除外してもよい。

ハプロタイプ推定装置１００は、上記の処理を各カウントパターンに対して行い、ｉイテレーションを終え、最終的に、データと矛盾しない全ディプロタイプを算出し、ハプロタイプ文字列ファイル１０６ｂに格納する。以上が、取り得るハプロタイプの組合せを求めるための処理である。

再び、図１３に戻り、ハプロタイプ推定部１０２ｃは、ハプロタイプ頻度算出部１０２ｄの処理により、ＥＭ法を用いて、ハプロタイプ文字列ファイル１０６ｂに記憶されたハプロタイプ文字列の組合せに基づいて、ハプロタイプ文字列の頻度を計算する。

すなわち、まず、ハプロタイプ頻度算出部１０２ｄは、各カウントパターン（ｃ_ｉ）に対し、各ハプロタイプの組合せの存在の重みを初期化する（ＳＣ−３）。本実施例においては、初期値として重みを平等に割り付けた。すなわち、初期値の重みは、ｗ_ｊｋ＝１／ｎ_ｊである。ここでｗは重み、ｊはカウントパターンのインデックス、ｋはカウントパターン内でのディプロタイプのインデックスであり、ｎ_ｊはカウントパターン内でのディプロタイプの総数である。

次にＳＣ−４に移り、ハプロタイプ推定部１０２ｃは、ＥＭ法のＭステップを数式１に基づいて行う。そして、対数尤度を下記の数式に従って計算し、記憶部１０６に保存する。下記の数式において、記号は上述と同様である。ここで対数尤度とは、計算されたディプロタイプの頻度がどれくらいデータを説明しているかの指標であり、本実施例においては、これをＥＭ法の収束の判定に用いる。ハプロタイプ推定部１０２ｃは、収束判定において、前回のＳＣ−４イテレーションで保存された対数尤度と今回のＳＣ−４イテレーションで計算された対数尤度との差を計算し、その差が一定値以内ならば、もはや対数尤度は改善されないと判定し、ＳＣ−５に処理を移す。そうでなければ、ＳＣ−４内にとどまり、ハプロタイプ推定部１０２ｃは、数式２に従って、ＥＭ法のＥステップを行う。そしてＳＣ−４のイテレーションを繰り返す。

そして、ハプロタイプ推定部１０２ｃは、上記条件によってＳＣ−５に処理を移した場合（対数尤度差一定値以内の場合）、ハプロタイプ及びその頻度を結果ファイルに（出力部１１４などに）出力する（ＳＣ−５）。以上で、本実施例の説明を終える。

［実証シミュレーション］
本実施の形態がコピー数多型及び一座位当たりの遺伝子型に関するデータから、ハプロタイプとその頻度を推定できるかどうかを確かめる為、シミュレーション実験を行った。実験は四つのタイプからハプロタイプとその頻度を推定できるかどうかシミュレーションを行った。すなわち、（１）一座位における二つの相同染色体に渡る総コピー数のデータから、その座位における相同染色体１本上のコピー数と、そのコピー数をアレルと見た時の頻度を推定できるかどうか、（２）（複数座位に渡る）一座位当たりの二つの相同染色体に渡る総コピー数のデータから、ハプロタイプとその頻度を推定できるかどうか、（３）一座位当たりの二つの相同染色体に渡る総コピー数のデータに、それとは異なる座位において、一座位当たりのＳＮＰの遺伝子型データが組み合わさったデータから、ハプロタイプとその頻度を推定できるかどうか、（４）一座位当たりの二つの相同染色体に渡る総コピー数が曖昧な総コピー数を含む場合のデータに、それとは異なる座位において、一座位当たりのＳＮＰ の遺伝子型データが組み合わさったデータから、ハプロタイプとその頻度を推定できるかどうか、についてシミュレーション実験を行った。図１５は、シミュレーションの枠組みを示すフローチャートである。

図１５に示すように、まず、ＳＤ−１において、ハプロタイプ推定装置は、ハプロタイプ（あるいはコピー数アレル）とその確率が書かれたファイルを読み込む。また、ハプロタイプ推定装置は、以下で用いる個体の数も読み込む。曖昧な総コピー数を含むデータを作成する場合は、曖昧でない総コピー数を曖昧な総コピー数に変換する表（例えば総コピー数０と１はいずれも０ｏｒ１に変換する、３以上の総コピー数はどれも＞２に変換する、等の対応づけが定義された表）が書かれたファイルも読み込む。

そして、ＳＤ−２において、ハプロタイプ推定装置は、読み込んだハプロタイプから構成しうるあらゆる２つ１組（ハプロタイプ２つから構成される１組）を作成し、１組を１つのディプロタイプとして、ハーディ・ワインバーグの法則（数式２−２参照）に基づいて、読み込んだハプロタイプの確率から全ディプロタイプの確率を計算する。

そして、ＳＤ−３において、ハプロタイプ推定装置は、ディプロタイプとその確率から、ディプロタイプの多項分布を構成し、与えられた個体の数だけディプロタイプをランダム抽出する。これは、例えば、Ｒ言語を使って、ｒｍｕｌｔｉｎｏｍ（１，ｓｉｚｅ＝個体の数，ｐｒｏｂ＝ｃ（ディプロタイプ１の確率、ディプロタイプ２の確率、ディプロタイプ３の確率、…））の様なコマンドで簡単に実施できる。抽出されたディプロタイプの１つ１つが、１つ１つの個体に相当する。

そして、ＳＤ−４において、ハプロタイプ推定装置は、各個体が持つディプロタイプを構成する２つのハプロタイプの各座位に対し、コピー数の座位の時は２つのハプロタイプに渡る総コピー数、ＳＮＰの座位の時は各塩基の数を数え、データを作成する。例えば、個体が持つディプロタイプが［１Ａ_２１／１Ａ_１２］であったら、Ｌ_１，Ｌ_３はコピー数の座位であるのでそれぞれ総コピー数は２（１＋１），３（１＋２）、Ｌ_２はＳＮＰの座位であるので、Ａ_１，Ａ_２の数を数えそれぞれ１，１となる。曖昧な総コピー数を含むデータを作成する場合は、入力で読み込んだ変換表に基づいて総コピー数を曖昧な総コピー数に変換する。

最後にＳＤ−５において、ハプロタイプ推定装置は、データをファイルに出力する。

以上で説明したシミュレーションの枠組みを使って、まず上述の実験（１）の場合について、図１６で示されたコピー数アレルとその確率、さらに個体の数５００を読み込み、データ図１７を作成した。次に、本実施の形態が図１７のデータだけから、図１６で示されたコピー数アレルとその確率（頻度）を推定（再現）できるか試験した。これは言い換えれば、不完全な観測データからの母集団比率の推定問題である。本実施の形態適用の際、ＥＭ法において収束判定に使われる対数尤度差は０．００１未満とした。適用の結果、コピー数アレルとその頻度に関し、図１８のような結果を得た。図１８（再現データ）にあって図１６（元データ）にないコピー数アレルの頻度は全て低い。上述のように、これらは存在しないと解釈される。図１６にあるコピー数アレルは全て図１８に現れており、かつ、それら推定頻度も正解頻度とほぼ等しい。よって、本実施の形態はコピー数アレルとその頻度を推定出来ることが実証された。

同様に、上述の実験（２）の場合についても、図１９で示されたハプロタイプとその確率、さらに個体の数５００を読み込み、データ図２０を作成した。次に本実施の形態が図２０のデータだけから、図１９で示されたハプロタイプとその確率を推定できるか試験した。本実施の形態適用の際、ＥＭ法において収束判定に使われる対数尤度差は上と同じとした。適用の結果、ハプロタイプとその頻度に関し、図２１のような結果を得た。図２１（再現データ）にあって図１９（元データ）にないハプロタイプの頻度は全て低い。上述のように、これらは存在しないと解釈される。図１９にあるハプロタイプは全て図２１に現れており、かつ、それら推定頻度も正解頻度とほぼ等しい。よって、本実施の形態はコピー数アレルの組合せとしてのハプロタイプとその頻度を推定出来ることが実証された。

同様に、上述の実験（３）の場合についても、図２２で示されたハプロタイプとその確率、さらに個体の数５００を読み込み、データ図２３を作成した。次に本実施の形態が図２３のデータだけから、図２２で示されたハプロタイプとその確率を推定できるか試験した。本実施の形態適用の際、ＥＭ法において収束判定に使われる対数尤度差は上と同じとした。適用の結果、ハプロタイプとその頻度に関し、図２４のような結果を得た。上述の実験結果と同様、この結果が示す通り、本実施の形態はコピー数アレルとＳＮＰアレルの組合せとしてのハプロタイプとその頻度を推定出来ることが実証された。

同様に、上述の実験（４）の場合についても、図２５で示されたハプロタイプとその確率、さらに個体の数５００を読み込み、データ図２６を作成した。この際、総コピー数０と１はいずれも曖昧な総コピー数０ｏｒ１に変換し、３以上の総コピー数はどれも曖昧な総コピー数＞２に変換した。次に、本実施の形態が図２６のデータだけから、図２５で示されたハプロタイプとその確率を推定できるか試験した。本実施の形態適用の際、ＥＭ法において収束判定に使われる対数尤度差は上と同じとした。適用の結果、ハプロタイプとその頻度に関し、図２７のような結果を得た。上述の実験結果と同様、この結果が示す通り、曖昧な総コピー数を含むデータの場合においても、本実施の形態はコピー数アレルとＳＮＰアレルの組合せとしてのハプロタイプとその頻度を推定出来ることが実証された。

［他の実施の形態］
さて、これまで本発明の実施の形態について説明したが、本発明は、上述した実施の形態以外にも、上記特許請求の範囲に記載した技術的思想の範囲内において種々の異なる実施の形態にて実施されてよいものである。

特に、上述の実施の形態において、上述した構成・処理以外にも、適宜、特願２００７−２３７１３９に記載の事項を補ってもよいものである。

また、ハプロタイプ推定装置１００がスタンドアローンの形態で処理を行う場合を一例に説明したが、ハプロタイプ推定装置１００とは別筐体で構成されるクライアント端末からの要求に応じて処理を行い、その処理結果を当該クライアント端末に返却するように構成してもよい。

また、実施の形態において説明した各処理のうち、自動的に行われるものとして説明した処理の全部または一部を手動的に行うこともでき、あるいは、手動的に行われるものとして説明した処理の全部または一部を公知の方法で自動的に行うこともできる。

このほか、上記文献中や図面中で示した処理手順、制御手順、具体的名称、各処理の登録データ等のパラメータを含む情報、データベース構成については、特記する場合を除いて任意に変更することができる。

また、ハプロタイプ推定装置１００に関して、図示の各構成要素は機能概略的なものであり、必ずしも物理的に図示の如く構成されていることを要しない。

例えば、ハプロタイプ推定装置１００の各装置が備える処理機能、特に制御部１０２にて行われる各処理機能については、その全部または任意の一部を、ＣＰＵ（Ｃｅｎｔｒａｌ Ｐｒｏｃｅｓｓｉｎｇ Ｕｎｉｔ）および当該ＣＰＵにて解釈実行されるプログラムにて実現することができ、あるいは、ワイヤードロジックによるハードウェアとして実現することも可能である。尚、プログラムは、後述する記録媒体に記録されており、必要に応じてハプロタイプ推定装置１００に機械的に読み取られる。すなわち、ＲＯＭまたはＨＤなどの記憶部１０６などは、ＯＳ（Ｏｐｅｒａｔｉｎｇ Ｓｙｓｔｅｍ）として協働してＣＰＵに命令を与え、各種処理を行うためのコンピュータプログラムが記録されている。このコンピュータプログラムは、ＲＡＭにロードされることによって実行され、ＣＰＵと協働して制御部を構成する。

また、このコンピュータプログラムは、ハプロタイプ推定装置１００に対して任意のネットワーク３００を介して接続されたアプリケーションプログラムサーバに記憶されていてもよく、必要に応じてその全部または一部をダウンロードすることも可能である。

また、本発明に係るプログラムを、コンピュータ読み取り可能な記録媒体に格納することもできる。ここで、この「記録媒体」とは、フレキシブルディスク、光磁気ディスク、ＲＯＭ、ＥＰＲＯＭ、ＥＥＰＲＯＭ、ＣＤ−ＲＯＭ、ＭＯ、ＤＶＤ等の任意の「可搬用の物理媒体」、あるいは、ＬＡＮ、ＷＡＮ、インターネットに代表されるネットワークを介してプログラムを送信する場合の通信回線や搬送波のように、短期にプログラムを保持する「通信媒体」を含むものとする。

また、「プログラム」とは、任意の言語や記述方法にて記述されたデータ処理方法であり、ソースコードやバイナリコード等の形式を問わない。なお、「プログラム」は必ずしも単一的に構成されるものに限られず、複数のモジュールやライブラリとして分散構成されるものや、ＯＳ（Ｏｐｅｒａｔｉｎｇ Ｓｙｓｔｅｍ）に代表される別個のプログラムと協働してその機能を達成するものをも含む。なお、実施の形態に示した各装置において記録媒体を読み取るための具体的な構成、読み取り手順、あるいは、読み取り後のインストール手順等については、周知の構成や手順を用いることができる。

記憶部１０６に格納される各種のデータベース等（多型テーブル１０６ａ〜実験データファイル１０６ｃ）は、ＲＡＭ、ＲＯＭ等のメモリ装置、ハードディスク等の固定ディスク装置、フレキシブルディスク、光ディスク等のストレージ手段であり、各種処理やウェブサイト提供に用いる各種のプログラムやテーブルやデータベースやウェブページ用ファイル等を格納する。

また、ハプロタイプ推定装置１００は、既知のパーソナルコンピュータ、ワークステーション等の情報処理装置を接続し、該情報処理装置に本発明の方法を実現させるソフトウェア（プログラム、データ等を含む）を実装することにより実現してもよい。

更に、装置の分散・統合の具体的形態は図示するものに限られず、その全部または一部を、各種の付加等に応じた任意の単位で、機能的または物理的に分散・統合して構成することができる。

以上詳述に説明したように、本発明によれば、塩基多型が考慮されていないコピー数多型に関する実験データから、曖昧なコピー数の総和しか得られない場合であっても、ハプロタイプとその頻度を高精度で推定することができる、ハプロタイプ推定装置、および、プログラムを提供することが可能である。

また、コピー数多型に関するデータに１座位当たりの遺伝子型データが加わったデータであっても、ハプロタイプとその頻度を推定することができる、ハプロタイプ推定装置、および、プログラムを提供することが可能である。

１座位あたりの遺伝子型データの一例を示す図である。 複数座位に渡るハプロタイプの組合せの一例を示す図である。 塩基多型を考慮した一塩基多型の実験データの一例を示す図である。 図３のデータから推定されるハプロタイプの表現例を示す図である。 １座位当たりの二つの相同染色体に渡る総コピー数（コピー数の総和）の多型テーブルの一例を表す図である。 曖昧なコピー数を含む、１座位当たりの二つの相同染色体に渡る総コピー数のデータの一例を示す図である。 １座位当たりの二つの相同染色体に渡る総コピー数のデータに、１座位当たりの遺伝子型データが組み合わさったデータの一例を示す図である。 本発明におけるハプロタイプの表現例を示す図である。 ２つのコピー数アレルからなる遺伝子型とそのコピー数の和を一例として示す図である。 本発明の概要を模式的に示したフロー図である。 本発明が適用される本ハプロタイプ推定装置の構成の一例を示すブロック図である。 本ハプロタイプ推定装置１００のハプロタイプ推定処理の一例を示す図である。 本実装の枠組みを示すフローチャートである。 コピー数総和分割処理およびハプロタイプ文字列格納処理の一例を示すフローチャートである。 シミュレーションの枠組みを示すフローチャートである。 実験（１）の場合について、シミュレーション実験に用いたハプロタイプとその頻度のデータ（元データ）の一例を示す図である。 実験（１）の場合について、読込まれたハプロタイプとその確率についての、作成された個体の数５００のコピー数データを示す図である。 実験（１）の場合について、本実施の形態を適用したシミュレーション結果（再現データ）を示す図である。 実験（２）の場合について、シミュレーション実験に用いたハプロタイプとその頻度のデータ（元データ）の一例を示す図である。 実験（２）の場合について、読込まれたハプロタイプとその確率についての、作成された個体の数５００のコピー数データを示す図である。 実験（２）の場合について、本実施の形態を適用したシミュレーション結果（再現データ）を示す図である。 実験（３）の場合について、シミュレーション実験に用いたハプロタイプとその頻度のデータ（元データ）の一例を示す図である。 実験（３）の場合について、読込まれたハプロタイプとその確率についての、作成された個体の数５００のコピー数データおよびカウント数データを示す図である。 実験（３）の場合について、本実施の形態を適用したシミュレーション結果（再現データ）を示す図である。 実験（４）の場合について、シミュレーション実験に用いたハプロタイプとその頻度のデータ（元データ）の一例を示す図である。 実験（４）の場合について、読込まれたハプロタイプとその確率についての、作成された個体の数５００のコピー数データおよびカウント数データを示す図である。 実験（４）の場合について、本実施の形態を適用したシミュレーション結果（再現データ）を示す図である。

符号の説明

１００ ハプロタイプ推定装置
１０２ 制御部
１０２ａ コピー数総和分割部
１０２ｂ ハプロタイプ文字列格納部
１０２ｃ ハプロタイプ推定部
１０２ｄ ハプロタイプ頻度算出部
１０２ｅ 多型テーブル作成部
１０４ 通信制御インターフェース部
１０６ 記憶部
１０６ａ 多型テーブル
１０６ｂ ハプロタイプ文字列ファイル
１０６ｃ 実験データファイル
１０８ 入出力制御インターフェース部
１１２ 入力部
１１４ 出力部
２００ 外部システム
３００ ネットワーク

Claims (8)

1. 集団における各個体の少なくともコピー数多型を含む実験データからハプロタイプを推定する、制御部と記憶部を少なくとも備えたハプロタイプ推定装置において、
   上記記憶部は、
   上記個体毎に、上記実験データから得られた、標識によって特定されるマーカー部位に対応付けられた上記コピー数多型のコピー数の総和Ｎに対する一または複数の条件式を、上記マーカー部位の種類毎に記憶する多型テーブル、
   を備え、
   上記制御部は、
   上記コピー数の総和Ｎについて、足し合わせた場合に、上記多型テーブルに記憶された上記コピー数の総和Ｎの上記条件式を満たす任意の２つの整数に分割するコピー数総和分割手段と、
   上記コピー数総和分割手段により分割された上記２つの整数をそれぞれ多型識別文字として表し、上記個体において上記多型識別文字を列挙したハプロタイプ文字列の組合せとして格納するハプロタイプ文字列格納手段と、
   上記集団において、同一である上記ハプロタイプ文字列の数を集計し、当該ハプロタイプ文字列の上記集団における頻度を求め、当該頻度が所定の条件を満たす上記各個体の上記ハプロタイプ文字列の上記組合せを、上記ハプロタイプの組合せとして推定するハプロタイプ推定手段と、
   を備えたことを特徴とするハプロタイプ推定装置。
2. 請求項１に記載のハプロタイプ推定装置において、
   上記条件式は、上記コピー数の総和Ｎの条件を規定する一または複数の等式または不等式を含むこと、
   を特徴とするハプロタイプ推定装置。
3. 請求項１に記載のハプロタイプ推定装置において、
   上記多型テーブルは、
   上記個体毎に、塩基多型を更に含む上記実験データから得られた、上記標識によって特定される上記マーカー部位に対応付けられた多型塩基をカウントしたカウント数を、当該多型塩基の種類毎に更に記憶し、
   上記ハプロタイプ文字列格納手段は、
   上記ハプロタイプ文字列の組合せにおいて、上記多型塩基の種類に対応付けた上記多型識別文字を更に列挙すること、
   を特徴とするハプロタイプ推定装置。
4. 請求項１に記載のハプロタイプ推定装置において、
   上記ハプロタイプ推定手段は、
   上記ハプロタイプ文字列の頻度を、ハーディ・ワインバーグの法則に基づいて算出し、
   上記所定の条件を、上記集団におけるハーディ・ワインバーグ平衡とすること、
   を特徴とするハプロタイプ推定装置。
5. 請求項１乃至４のいずれか一つに記載のハプロタイプ推定装置において、
   上記ハプロタイプ文字列格納手段は、
   作成した上記ハプロタイプ文字列の組合せにおける上記多型識別文字が表す上記整数の和が、上記多型テーブルにおける上記コピー数の総和Ｎの上記条件式を満たすか否か確認し、上記条件式を満たさない場合に当該ハプロタイプ文字列を除外し、および／または、作成した上記ハプロタイプ文字列の組合せが他の上記ハプロタイプ文字列の組合せと重複するか否か判断し、重複する場合に当該ハプロタイプ文字列の組合せを除外すること、
   を特徴とするハプロタイプ推定装置。
6. 請求項１乃至５のいずれか一つに記載のハプロタイプ推定装置において、
   上記ハプロタイプ推定手段は、
   ＥＭ（Ｅｘｐｅｃｔａｔｉｏｎ−Ｍａｘｉｍｉｚａｔｉｏｎ）法を用いて、
   上記集団における上記ハプロタイプ文字列の頻度を、当該ハプロタイプ文字列を少なくとも一方に有する上記組合せの頻度により重み付けして算出するＭステップと、
   上記組合せの頻度を、当該組合せを構成する上記ハプロタイプ文字列の頻度の積により求め、当該組合せの頻度に基づいて上記重みを算出するＥステップと、
   を上記頻度の値が収束するまで交互に繰り返すハプロタイプ頻度算出手段、
   を更に備えたことを特徴とするハプロタイプ推定装置。
7. 請求項６に記載のハプロタイプ推定装置において、
   上記Ｍステップは、
   下記の数式１に基づいて、上記ハプロタイプ文字列の頻度を算出し、
   

   

   （ここで、Ｐ（ｈ_ｉ）は上記ハプロタイプ文字列の頻度を表し、ｈは上記ハプロタイプ文字列を表し、ｉは上記ハプロタイプ文字列のインデックスを表す。また、ｎは上記集団を構成する上記個体の数、ｊは上記多型テーブルにおける上記コピー数の総和Ｎの上記条件式および／または上記カウント数に関するパターンのインデックス、ｋは上記ハプロタイプ文字の組合せのインデックス、Ｎ（ｃ_ｊ）は上記パターンｊを持つ上記個体の数を表す。また、δ（ｈ_ｉ，ｄ_ｊｋ）は、上記組合せｄ_ｊｋが一方に当該ハプロタイプ文字列ｈ_ｉを有する場合に１を返し、両方に当該ハプロタイプ文字列ｈ_ｉを有する場合に２を返し、当該ハプロタイプ文字列ｈ_ｉを持たない場合に０を返す関数であり、ｄは上記ハプロタイプ文字列の上記組合せを表す。また、ｗ_ｊｋは上記ハプロタイプ文字列の上記組合せの頻度による上記重みである。）
   上記Ｅステップは、
   下記の数式２に基づいて上記ハプロタイプ文字列の上記組合せの頻度を求め、上記重みとして当該ハプロタイプ文字列の上記組合せの頻度を上記集団における上記組合せの頻度の総和で除して算出すること、
   

   

   （ここで、Ｐ（ｄ_ｊｋ）は、上記ハプロタイプ文字列の上記組合せの頻度を表す。また、ｈ_ｌおよびｈ_ｍは当該組合せを構成する２つの上記ハプロタイプ文字列を表し、Ｐ（ｈ_ｌ）およびＰ（ｈ_ｍ）は、当該２つの上記ハプロタイプ文字列の頻度をそれぞれ表す。）
   を特徴とするハプロタイプ推定装置。
8. 集団における各個体の少なくともコピー数多型を含む実験データからハプロタイプを推定する、制御部と記憶部を少なくとも備えたハプロタイプ推定装置に実行させるためのプログラムであって、
   上記記憶部は、
   上記個体毎に、上記実験データから得られた、標識によって特定されるマーカー部位に対応付けられた上記コピー数多型のコピー数の総和Ｎに対する一または複数の条件式を、上記マーカー部位の種類毎に記憶する多型テーブル、
   を備えており、
   上記制御部において、
   上記コピー数の総和Ｎについて、足し合わせた場合に、上記多型テーブルに記憶された上記コピー数の総和Ｎの上記条件式を満たす任意の２つの整数に分割するコピー数総和分割ステップと、
   上記コピー数総和分割ステップにて分割された上記２つの整数をそれぞれ多型識別文字として表し、上記個体において上記多型識別文字を列挙したハプロタイプ文字列の組合せとして格納するハプロタイプ文字列格納ステップと、
   上記集団において、同一である上記ハプロタイプ文字列の数を集計し、当該ハプロタイプ文字列の上記集団における頻度を求め、当該頻度が所定の条件を満たす上記各個体の上記ハプロタイプ文字列の上記組合せを、上記ハプロタイプの組合せとして推定するハプロタイプ推定ステップと、
   を実行させるためのプログラム。

Images