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

Description

本発明は、ハプロタイプ推定装置、および、プログラムに関し、特に、塩基多型を考慮したコピー数多型に関する実験データから、コンピュータによって効率的にコピー数多型のハプロタイプとその頻度を推定するハプロタイプ推定装置、および、プログラムに関する。

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

従来、１座位当たりの遺伝子型データから、複数座位に渡るハプロタイプを推定していた。ここで、１座位当たりの遺伝子型データとは、複数座位に渡る関連が分からない、座位毎に独立した（相（phase）が未知の）遺伝子型データのことである。ここで、図１は、１座位あたりの遺伝子型データの一例を示す図である。図１において、Ｌは座位（Locus）、Ａはアレル（Allele）を表している。

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

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

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

図２に示すように、例えば、ハプロタイプ１は、座位Ｌ₁においてＡ₁のアレルをもち、座位Ｌ₂においてＡ₁のアレルをもち、座位Ｌ₃においてＡ₁のアレルを持つことが特定されている。このように、従来のハプロタイプの推定方法においては、一般に２種類のアレルを想定し、一塩基多型（Single Nucleotide Polymorphism、「ＳＮＰ」と略す。）などの遺伝子型データから、複数座位に渡るハプロタイプを推定していた。

ここで、ＳＮＰなどの塩基多型の他に、コピー数多型（Copy Number Polymorphism, あるいはCopy Number Variation、本明細書中で「ＣＮＰ」と略す場合がある。）と呼ばれる多型が存在する。ここで、図３は、コピー数多型および塩基多型の一例を示す図である。図３において、Ｍは、蛍光色素プローブ等の標識で識別される個体間で違いのない配列部位（マーカー部位）、Ｆは、（異なる蛍光色素等で区別される）個体間で違いうる塩基（多型塩基）に、それぞれ対応する。

図３に示すように、コピー数多型においては、ある区間の配列（「コピー単位」と呼ぶ。）が繰り返し現れることがあり、そのコピー数に個体差がある。例えば、図３に示すように、互いに相同な染色体（染色体１〜４）において、染色体１ではコピー数が１であり、染色体２ではコピー数が０であり、染色体３ではコピー数が２であり、染色体４ではコピー数が３、とそれぞれ異なっている。

ここで、コピー単位上に多型塩基が存在する場合、コピー数多型がないゲノム領域における塩基多型とは異なり、その多型塩基のカウント数は、個体における２つのハプロタイプ（すなわちディプロタイプ）のコピー数に依存することになる。すなわち、コピー数多型がないゲノム領域における塩基多型では、カウント数は、基本的に０、１または２であるが（有性生殖をする個体は相同染色体として２本１組の染色体を持つため）、コピー数単位上に塩基多型がある場合、その多型塩基のカウント数は、０、１、２、３、４、５、・・・と、コピー数に依存して各個体で変化する。すなわち、コピー数単位上に塩基多型がある場合のカウント数は、従来技術におけるように１座位当たりの遺伝子型に直結するものではない。なお、コピー数多型におけるカウント数とは、コピー単位上の、標識によって特定されるマーカー部位に対応付けられた多型塩基をカウントすることにより得られたカウント数を意味する。

ここで、非特許文献５は、コピー数多型に対するハプロタイプ推定方法であり、アレルをコピー数が多いアレル、少ないアレルの２種類に分け、複数座位に渡るハプロタイプを推定する方法である。

また、コピー数多型に対するより一般的なハプロタイプ推定方法としては、コピー数をアレルの種類に対応させ、複数種類のアレルを想定し、複数座位に渡るハプロタイプを推定する方法がある。

また、特許文献１は、ハプロタイプの頻度を算出する方法として、ＥＭ（Expectation-Maximization）法を用いるものである。

特開２００４−１９２０１８号公報 チアンフア ニウ（Tianhua Niu）著 「アルゴリズムズ フォー インファリング ハプロタイプス（Algorithms for inferring haplotypes）」 Genet Epidemiol.、２００４年１２月、２７巻（４号）３３４−３４７頁 ジャオフイＳキン、チアンフア ニウ、ジュンＳリウ（Zhaohui S. Qin, Tianhua Niu, Jun S. Liu）著 「パーティション−ライゲーション−エクスペクテーション−マキシマイゼーション アルゴリズム フォー ハプロタイプ インフェアレンス ウィズ シングル−ヌクレオタイド ポリモーフィズムズ（Partition−ligation−expectation-maximization algorithm for haplotype inference with single-nucleotide polymorphisms）」Am J Hum Genet.、２００２年１１月、７１巻（５号）１２４２−１２４７頁 ローレント エクスコファー、モンゴメリー スラキン（Laurent Excoffier, Montgomery Slatkin）著 「マキシマム−ライクリフッド エスティメーション オブ モレキュラー ハプロタイプ フロークェンシーズ イン ア ディプロイド ポピュレーション（Maximum-likelihood estimation of molecular haplotype frequencies in a diploid population）」 Mol Biol Evol.、１９９５年９月、１２巻（５号）、９２１−９２７頁 Ｍ．Ｅ．ホーリー、Ｋ．Ｋ．キッド（M. E. Hawley, K. K. Kidd）著 「ハプロ： ア プログラム ユージング ザ ＥＭ アルゴリズム トゥー エスティメイト ザ フリークェンシーズ オブ マルチ−サイト ハプロタイプズ（HAPLO: a program using the EM algorithm to estimate the frequencies of multi-site haplotypes）」 J Hered.、１９９５年９−１０月、８６巻（５号）、４０９−４１１頁 リチャード レンドン、シュンペイ イシカワ、カレンＲフィッチ他、全４３名（Richard Redon, Shumpei Ishikawa, Karen R. Fitch, et al.）著 「グローバル バリエーション イン コピー ナンバー イン ザ ヒューマン ゲノム（Global variation in copy number in the human genome）」 ネーチャー（Nature）出版 ２００６年１１月２３日、４４４巻（７１１８号）４４４−４５４頁

しかしながら、従来、塩基多型を考慮したコピー数多型に関する実験データから、コンピュータによってコピー数多型のハプロタイプとその頻度を推定する方法は開発されていない、という問題があった。

すなわち、従来技術においては、コピー数をアレルの種類と仮定することにより、個体はコピー数１の染色体を２本持つものとしてハプロタイプの推定を行うので、コピー数の概念は捨象されてしまうという問題がある。

また、従来のコピー数多型に対するハプロタイプ推定方法においては、コピー数をアレルの種類に対応付けるため、原理的に、塩基多型を考慮したコピー数多型に関する実験データから、ハプロタイプを推定することができない、という問題がある。

また、特許文献１に記載のＥＭ法によるハプロタイプ頻度推定方法は、コピー数多型を考慮したハプロタイプの推定方法として応用することができない、という問題がある。

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

このような目的を達成するため、請求項１に記載のハプロタイプ推定装置は、集団における各個体のコピー数多型と塩基多型を含む遺伝子型データからハプロタイプを推定する、制御部と記憶部を少なくとも備えたハプロタイプ推定装置において、上記記憶部は、上記個体毎に、上記遺伝子型データを用いて、上記コピー数多型のコピー単位上の、標識によって特定されるマーカー部位に対応付けられた多型塩基をカウントすることにより得られたカウント数を、当該多型塩基の種類毎に記憶するカウント数テーブル、を備え、上記制御部は、上記カウント数テーブルに格納された上記個体の上記カウント数に基づいて、上記マーカー部位毎に上記カウント数の和を集計し、上記カウント数の和のうちの最大値を求める最大値算出手段と、上記多型塩基の種類に対応付けられた多型識別文字を、当該多型塩基の上記カウント数列挙する多型識別文字列挙手段と、上記コピー単位が上記最大値の回数繰り返される任意の文字列を作成するよう、列挙された上記多型識別文字を並び替える文字列作成手段と、上記文字列作成手段によって作成された上記文字列を、任意の上記コピー単位で２分し、ハプロタイプ文字列の組合せとして格納するハプロタイプ文字列格納手段と、上記集団において、同一である上記ハプロタイプ文字列の数を集計し、当該ハプロタイプ文字列の上記集団における頻度を求め、当該頻度が所定の条件を満たす上記各個体の上記ハプロタイプ文字列の上記組合せを、上記ハプロタイプの組合せとして推定するハプロタイプ推定手段と、を備えたことを特徴とする。

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

また、請求項３のハプロタイプ推定装置は、請求項１または２に記載のハプロタイプ推定装置において、上記多型識別文字列挙手段は、上記多型塩基の種類として塩基欠失を表す上記多型識別文字を、上記最大値から上記和を減じた数だけ列挙すること、を特徴とする。

また、請求項４のハプロタイプ推定装置は、請求項１乃至３のいずれか一つに記載のハプロタイプ推定装置において、上記文字列作成手段は、作成した上記文字列における、上記多型識別文字の文字数と、上記カウント数テーブルにおける、当該多型塩基の種類毎の上記カウント数と、が一致することを確認し、一致しない場合に当該文字列を除外すること、を特徴とする。

また、請求項５に記載のハプロタイプ推定装置は、請求項１乃至４のいずれか一つに記載のハプロタイプ推定装置において、上記制御部は、上記コピー単位を上記マーカー部位の単位で複数の区画に分け、上記区画に含まれる上記マーカー部位の上記塩基多型について、上記最大値算出手段、上記多型識別文字列挙手段、上記文字列作成手段、上記ハプロタイプ文字列格納手段、および、上記ハプロタイプ推定手段に、当該区画毎に処理を実行させるよう制御する区画化手段と、上記区画毎に上記ハプロタイプ推定手段により推定された上記ハプロタイプ文字列について、上記コピー単位の繰り返しパターンが同じ上記ハプロタイプ文字列同士で、上記コピー単位における上記各区画を互いに連結することにより、上記複数の上記区画から構成される上記ハプロタイプ文字列を再現する連結手段と、を備えたことを特徴とする。

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

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

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

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

また、請求項８に記載のハプロタイプ推定装置は、請求項１乃至７のいずれか一つに記載のハプロタイプ推定装置において、上記記憶部は、上記各個体における上記コピー数多型のコピー数の総和を記憶するコピー数総和記憶手段、を更に備え、上記制御部は、上記コピー数総和記憶手段に記憶された上記各個体における上記コピー数の総和を、任意の一つの上記多型塩基の上記カウント数とした上記カウント数テーブルを作成するカウント数テーブル作成手段と、上記ハプロタイプ推定手段により推定された上記ハプロタイプの組み合わせにおいて、上記ハプロタイプ文字列中の上記一つの多型塩基に対応付けられた上記多型識別文字の数を算出し、算出した数を上記ハプロタイプにおける上記コピー数として設定するコピー数設定手段と、を更に備えたことを特徴とする。

また、本発明のハプロタイプ推定装置は、請求項６または７に記載のハプロタイプ推定装置において、上記ハプロタイプ頻度算出手段は、上記Ｍステップにおいて算出された上記ハプロタイプ文字列の頻度と、前回の上記Ｍステップにおいて算出された上記ハプロタイプ文字列の頻度と、の対数尤度差を求め、上記対数尤度差が所定の閾値以下となった場合に、上記頻度の値が収束したと判定すること、を特徴とする。

また、請求項９に記載のプログラムは、集団における各個体のコピー数多型と塩基多型を含む遺伝子型データからハプロタイプを推定する、制御部と記憶部を少なくとも備えたハプロタイプ推定装置にハプロタイプ推定方法を実行させるためのプログラムであって、上記記憶部は、上記個体毎に、上記遺伝子型データを用いて、上記コピー数多型のコピー単位上の、標識によって特定されるマーカー部位に対応付けられた多型塩基をカウントすることにより得られたカウント数を、当該多型塩基の種類毎に記憶するカウント数テーブル、を備えており、上記制御部において実行される、上記カウント数テーブルに格納された上記個体の上記カウント数に基づいて、上記マーカー部位毎に上記カウント数の和を集計し、上記カウント数の和のうちの最大値を求める最大値算出ステップと、上記多型塩基の種類に対応付けられた多型識別文字を、当該多型塩基の上記カウント数列挙する多型識別文字列挙ステップと、上記コピー単位が上記最大値の回数繰り返される任意の文字列を作成するよう、列挙された上記多型識別文字を並び替える文字列作成ステップと、上記文字列作成ステップにおいて作成された上記文字列を、任意の上記コピー単位で２分し、ハプロタイプ文字列の組合せとして格納するハプロタイプ文字列格納ステップと、上記集団において、同一である上記ハプロタイプ文字列の数を集計し、当該ハプロタイプ文字列の上記集団における頻度を求め、当該頻度が所定の条件を満たす上記各個体の上記ハプロタイプ文字列の上記組合せを、上記ハプロタイプの組合せとして推定するハプロタイプ推定ステップと、を含むハプロタイプ推定方法を実行させることを特徴とする。

また、本発明は、ハプロタイプ推定方法に関するものであり、集団における各個体のコピー数多型と塩基多型を含む遺伝子型データからハプロタイプを推定する、制御部と記憶部を少なくとも備えたハプロタイプ推定装置において実行されるハプロタイプ推定方法であって、上記記憶部は、上記個体毎に、上記遺伝子型データを用いて、上記コピー数多型のコピー単位上の、標識によって特定されるマーカー部位に対応付けられた多型塩基をカウントすることにより得られたカウント数を、当該多型塩基の種類毎に記憶するカウント数テーブル、を備えており、上記制御部において実行される、上記カウント数テーブルに格納された上記個体の上記カウント数に基づいて、上記マーカー部位毎に上記カウント数の和を集計し、上記カウント数の和のうちの最大値を求める最大値算出ステップと、上記多型塩基の種類に対応付けられた多型識別文字を、当該多型塩基の上記カウント数列挙する多型識別文字列挙ステップと、上記コピー単位が上記最大値の回数繰り返される任意の文字列を作成するよう、列挙された上記多型識別文字を並び替える文字列作成ステップと、上記文字列作成ステップにおいて作成された上記文字列を、任意の上記コピー単位で２分し、ハプロタイプ文字列の組合せとして格納するハプロタイプ文字列格納ステップと、上記集団において、同一である上記ハプロタイプ文字列の数を集計し、当該ハプロタイプ文字列の上記集団における頻度を求め、当該頻度が所定の条件を満たす上記各個体の上記ハプロタイプ文字列の上記組合せを、上記ハプロタイプの組合せとして推定するハプロタイプ推定ステップと、を含むことを特徴とする。

また、本発明は、上記ハプロタイプ推定方法において、上記ハプロタイプ推定ステップは、上記ハプロタイプ文字列の頻度を、ハーディ・ワインバーグの法則に基づいて算出し、上記所定の条件を、上記集団におけるハーディ・ワインバーグ平衡とすること、を特徴とする。

また、本発明は、上記ハプロタイプ推定方法において、上記多型識別文字列挙ステップは、上記多型塩基の種類として塩基欠失を表す上記多型識別文字を、上記最大値から上記和を減じた数だけ列挙すること、を特徴とする。

また、本発明は、上記ハプロタイプ推定方法において、上記文字列作成ステップは、作成した上記文字列における、上記多型識別文字の文字数と、上記カウント数テーブルにおける、当該多型塩基の種類毎の上記カウント数と、が一致することを確認し、一致しない場合に当該文字列を除外すること、を特徴とする。

また、本発明は、上記ハプロタイプ推定方法において、上記制御部において実行される、上記コピー単位を上記マーカー部位の単位で複数の区画に分け、上記区画に含まれる上記マーカー部位の上記塩基多型について、上記最大値算出ステップ、上記多型識別文字列挙ステップ、上記文字列作成ステップ、上記ハプロタイプ文字列格納ステップ、および、上記ハプロタイプ推定ステップにおいて、当該区画毎に処理が実行されるよう制御する区画化ステップと、上記区画毎に上記ハプロタイプ推定ステップにより推定された上記ハプロタイプ文字列について、上記コピー単位の繰り返しパターンが同じ上記ハプロタイプ文字列同士で、上記コピー単位における上記各区画を互いに連結することにより、上記複数の上記区画から構成される上記ハプロタイプ文字列を再現する連結ステップと、を含むことを特徴とする。

また、本発明は、上記ハプロタイプ推定方法において、上記ハプロタイプ推定ステップは、ＥＭ（Expectation-Maximization）法を用いて、上記集団における上記ハプロタイプ文字列の頻度を、当該ハプロタイプ文字列を少なくとも一方に有する上記組合せの頻度により重み付けして算出するＭステップと、上記組合せの頻度を、当該組合せを構成する上記ハプロタイプ文字列の頻度の積により求め、当該組合せの頻度に基づいて上記重みを算出するＥステップと、を上記頻度の値が収束するまで交互に繰り返すハプロタイプ頻度算出ステップ、を更に含むことを特徴とする。

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

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

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

また、本発明は、上記ハプロタイプ推定方法において、上記記憶部は、上記各個体における上記コピー数多型のコピー数の総和を記憶するコピー数総和記憶手段、を更に備えており、上記制御部において実行される、上記コピー数総和記憶手段に記憶された上記各個体における上記コピー数の総和を、任意の一つの上記多型塩基の上記カウント数とした上記カウント数テーブルを作成するカウント数テーブル作成ステップと、上記ハプロタイプ推定ステップにおいて推定された上記ハプロタイプの組み合わせにおいて、上記ハプロタイプ文字列中の上記一つの多型塩基に対応付けられた上記多型識別文字の数を算出し、算出した数を上記ハプロタイプにおける上記コピー数として設定するコピー数設定ステップと、を更に含むことを特徴とする。

また、本発明は、上記ハプロタイプ推定方法において、上記ハプロタイプ頻度算出ステップは、上記Ｍステップにおいて算出された上記ハプロタイプ文字列の頻度と、前回の上記Ｍステップにおいて算出された上記ハプロタイプ文字列の頻度と、の対数尤度差を求め、上記対数尤度差が所定の閾値以下となった場合に、上記頻度の値が収束したと判定すること、を特徴とする。

本発明によれば、１座位当たりの遺伝子型データに直結しない、塩基多型を考慮したコピー数多型に関する実験データから、ハプロタイプとその頻度を高精度で推定することができる。

また、本発明によれば、カウント数が各マーカー部位間で一致しない場合であっても、不足するカウント数を塩基欠失に起因する結果であるとして処理を行うので、実験により検出することができない塩基欠失の多型も考慮して、精度の高いハプロタイプ推定装置を提供できる。

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

また、本発明によれば、パーティション−ライゲーション法（ＰＬ法）を用いて、コピー単位を複数の区画に分けて処理した結果を統合化するので、マーカー部位の数が増えた場合に取り得るハプロタイプの数が飛躍的に増大する場合でも、計算時間を節約することができる。

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

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

図１は、１座位あたりの遺伝子型データの一例を示す図である。 図２は、複数座位に渡るハプロタイプの組合せの一例を示す図である。 図３は、コピー数多型および塩基多型の一例を示す図である。 図４は、塩基多型を考慮したコピー数多型の実験データの一例を示す図である。 図５は、本発明の概要を模式的に示したフロー図である。 図６は、本発明のハプロタイプ表現形式の一例を示す図である。 図７は、本発明が適用される本ハプロタイプ推定装置の構成の一例を示すブロック図である。 図８は、本ハプロタイプ推定装置１００のハプロタイプ推定処理の一例を示す図である。 図９は、カウント数データの全パターンにおいて、そのカウント数パターンと矛盾しないディプロタイプを作成する処理の一例を示すフローチャートである。 図１０は、コピー数多型を考慮しないハプロタイプ推定方法のＰＬ法における連結（ligation）処理の一例を示す図である。 図１１は、本実施の形態の、コピー数多型を考慮に入れたハプロタイプ推定方法のＰＬ法における連結（ligation）処理の一例を示す図である。 図１２は、連結ステップにおける合理的な連結の方法を示す図である。 図１３は、本実施の形態２における処理の一例を示すフローチャートである。 図１４は、各個体におけるコピー数の総和データと、作成されたカウント数テーブルの対応関係の一例を示す図である。 図１５は、ＥＭ法によるハプロタイプ推定処理の一例を示すフローチャートである。 図１６は、シミュレーションの枠組みを示すフローチャートである。 図１７は、読込むハプロタイプ推定の結果の一例を表すデータ図である。 図１８は、ハプロタイプとその確率と同様のデータを、さらに個体の数５００を読み込み、出力されたカウント数データを示す図である。 図１９は、ハプロタイプ推定装置１００による、ハプロタイプとその頻度に関し出力された結果を示す図である。 図２０は、本実装の枠組みを示すフローチャートである。 図２１は、ＳＦ−５の、より詳しい処理の一例を示すフローチャートである。 図２２は、シミュレーション実験に用いたハプロタイプとその頻度のデータの一例を示す図である。 図２３は、読込まれたハプロタイプとその確率についての、作成された個体の数５００のカウント数データを示す図である。 図２４は、本実施例を適用した結果を示す図である。 図２５は、本実施例において設定した、ハプロタイプに存在するコピー数（コピー単位の数）とそのハプロタイプの確率を示す図である。 図２６は、図２５に示すコピー数を仮の塩基Ｆ_１のカウント数として表現した図である。 図２７は、抽出された各個体におけるコピー数の総和を示す図である。 図２８は、図２７に示すコピー数の総和を仮の塩基Ｆ_１のカウント数として表現した図である。 図２９は、ハプロタイプ推定装置１００により推定されたハプロタイプ文字列とその頻度を示す図である。 図３０は、推定されたハプロタイプ文字列から算出された各コピー数のハプロタイプとその頻度を示す図である。

符号の説明

１００ ハプロタイプ推定装置
１０２ 制御部
１０２ａ 最大値算出部
１０２ｂ 多型識別文字列挙部
１０２ｃ 文字列作成部
１０２ｄ ハプロタイプ文字列格納部
１０２ｅ ハプロタイプ推定部
１０２ｆ ハプロタイプ頻度算出部
１０２ｇ 区画化部
１０２ｈ 連結部
１０２ｊ カウント数テーブル作成部
１０２ｋ コピー数設定部
１０４ 通信制御インターフェース部
１０６ 記憶部
１０６ａ カウント数テーブル
１０６ｂ ハプロタイプ文字列ファイル
１０６ｃ 遺伝子型データファイル
１０６ｄ コピー数総和ファイル
１０８ 入出力制御インターフェース部
１１２ 入力部
１１４ 出力部
２００ 外部システム
３００ ネットワーク

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

［本発明の概要］
以下、本発明の概要について説明し、その後、本発明の構成および処理等について詳細に説明する。まず、コピー数多型がある場合とない場合で、ハプロタイプの推定の上でどのような違いがあるかについて説明する。

コピー数多型がないゲノム領域上の塩基多型では、カウント数は、基本的に０、１または２である。すなわち、上記図３の例でいえば、全ての個体のハプロタイプが、染色体１のようにコピー数を１しか持たない場合である。すなわち、コピー数多型がないゲノム領域上の塩基多型においては、各マーカー部位において、異なる各多型塩基の種類どうしのカウント数を総計すると、総数は必ず２となる。

例えば、図３における染色体１を２本持つ場合の例では、ハプロタイプは、[F₁(M₁)F₃(M₂)F₁(M₃)/F₁(M₁)F₃(M₂)F₁(M₃)]（"/"で２本の相同染色体を区別する。以下、(M)によるマーカー位置表示を省略する。）で表記され、M₁, M₂, M₃ に、それぞれ[F₁/F₁], [F₃/F₃], [F₁/F₁]の多型塩基が対応している（１座位当たりの遺伝子型に直結している）。すなわち、実験によって各マーカー部位、各塩基のカウント数が、M₁においてF₁, F₂, F₃ はそれぞれ2, 0, 0、M₂ においてそれぞれ0, 0, 2、M₃においてそれぞれ2, 0, 0 と識別され、総数はいずれのマーカー部位でも２(M₁ で2+0+0, M₂ で0+0+2, M₃で2+0+0)となる。

一方、コピー単位上に多型塩基が存在する場合では、１本の染色体に乗っているコピー単位の数は１以外の場合があるため、上記と異なり各マーカー部位において、各塩基のカウント数の総数が２以外となる場合がある。すなわち、カウント数の和は、個体における２つのハプロタイプのコピー数に依存し、０、１、２、３、４、５、・・・と、コピー数に依存的に各個体で変化する。このため、コピー数単位上に塩基多型がある場合のカウント数は、従来技術におけるように１座位当たりの遺伝子型に直結しないので、ハプロタイプの推定において工夫が必要となる。

例えば、図３における染色体１と染色体３とを相同染色体として持っている個体において、[F₁F₃F₁/F₁F₃F₁, F₂F₃F₁]（"," でコピー単位を区別する）が相同染色体に乗っており、M₁, M₂, M₃ に、それぞれ[F₁/F₁, F₂], [F₃/F₃, F₃], [F₁/F₁, F₁] が対応することになる。実験によって各マーカー部位、各塩基のカウント数が、M₁ においてF₁, F₂, F₃ はそれぞれ2, 1, 0、M₂ においてそれぞれ0, 0, 3、M₃ においてそれぞれ3, 0, 0 と識別され、総数は各マーカー部位で3 (M₁ で2+1+0, M₂ で0+0+3, M₃で3+0+0)となり、ハプロタイプのコピー数の和と一致する。

ここで、もし図３における染色体１と染色体２を持ったときは、[F₁F₃F₁/F₀F₀F₀] （F₀ は欠失を表し、実験的には塩基を識別できない。）となり、M₁, M₂, M₃に、それぞれ[F₁/F₀], [F₃/F₀], [F₁/F₀] が対応することになる。M₁においてF₁, F₂, F₃ はそれぞれ1, 0, 0、M₂ においてそれぞれ0, 0, 1、M₃ においてそれぞれ1, 0, 0 と識別され、総数は各マーカー部位で１となり、ハプロタイプのコピー数と一致する。

ここで、図４は、塩基多型を考慮したコピー数多型の実験データから得られたコピー数データの一例を示す図である。Ｍは、蛍光色素プローブ等の標識で識別される個体間で違いのない配列部位（マーカー部位）、Ｆは、（異なる蛍光色素等で区別される）個体間で違いうる塩基（多型塩基）に、それぞれ対応する。

図４において、個体１においては、M₁ で2+1+0, M₂ で1+1+1, M₃ で2+0+1で総数は３であり、上述の性質からコピー数の和は３と予測できる。一方、個体２においては、M₁で0+0+0, M₂ で0+3+0, M₃ で2+2+0であり、一致していないが、塩基欠失があることを鑑みるとコピー数の和は４（M₁では塩基欠失が４箇所、M₂では塩基欠失が１箇所と予想）であると考えられる。

以上をまとめると、塩基多型を考慮したコピー数多型に関する実験データは、蛍光色素プローブ等の標識によって実験的に識別される多型塩基のカウント数を記したカウント数データを含む。そして、当該カウント数データは、コピー数多型がないゲノム領域上の塩基多型とは異なり、１座位当たりの遺伝子型データに直結するものではないが、カウント数の総数は２つのハプロタイプのコピー数の和と一致する（塩基欠失の場合を除く）。本発明は、この性質を利用して、塩基多型を考慮したコピー数多型に関する実験データから、ハプロタイプを推定するものである。すなわち、本発明は、多型塩基をカウントし得られたカウント数データからハプロタイプとその頻度を推定するために、カウント数データと矛盾しないディプロタイプ（ハプロタイプの組合せ）を見出し、カウント数データを利用してハプロタイプの頻度を計算する。ここで、「矛盾しない」とは、ディプロタイプから数えられる各マーカー部位、各多型塩基のカウント数が、カウント数データにおける各マーカー部位、各塩基のカウント数と一致する、ということである。カウント数データと矛盾しないディプロタイプが見出せれば、ＥＭ（Expectation-Maximization ）法やGibbsサンプリング法などによって、カウント数データからハプロタイプの頻度が計算できる。以下に、本発明の概要を上記で各図の説明で用いた記号で補足しながら説明する。ここで、図５は、本発明の概要を模式的に示したフロー図である。

まず、本発明は、集団における各個体のコピー数多型と塩基多型を含む遺伝子型データからハプロタイプを推定する、制御部と記憶部を少なくとも備えたハプロタイプ推定装置または当該装置において実行される方法、そのプログラムである。

そして、本ハプロタイプ推定装置は、コピー数多型のコピー単位（ｕ）上の、標識によって特定されるマーカー部位（Ｍ）に対応付けられた多型塩基（Ｆ）をカウントすることにより得られたカウント数（ｃ₁）を、当該多型塩基の種類（Ｆ₁、Ｆ₂）毎、および、個体（個体１、個体２、・・・）毎に記憶するカウント数テーブルを備える。ここで、本ハプロタイプ推定装置は、記憶部に記憶された各個体におけるコピー数多型のコピー数の総和を、任意の一つの多型塩基のカウント数としたカウント数テーブルを作成してもよい。図５＜ＳＡ−１＞は、ある個体についてのカウント数について表されたカウント数テーブルの例である。

そして、本ハプロタイプ推定装置は、カウント数テーブルに格納された個体（例えば、個体１）のカウント数（c₁）に基づいて、マーカー部位（Ｍ₁、Ｍ₂、Ｍ₃）毎にカウント数の和（2+1、1+1、0+1）を集計し、カウント数の和のうちの最大値（この場合、３）を求める。

そして、本ハプロタイプ推定装置は、多型塩基の種類に対応付けられた多型識別文字（例えば、Ｆ₁、Ｆ₂）を、多型塩基のカウント数（例えば、Ｍ₁についてＦ₁は２個）ずつ列挙する（図５＜ＳＡ−２＞参照）。ここで、図５＜ＳＡ−２＞に示すように、本ハプロタイプ推定装置は、多型塩基の種類として塩基欠失を表す多型識別文字（Ｆ₀）を、最大値からカウント数の和を減じた数（例えば、Ｍ₂については3-2=１個）だけ列挙してもよい。

そして、本ハプロタイプ推定装置は、図５＜ＳＡ−３〜４＞に示すように、コピー単位が最大値の回数（この場合３回）繰り返される任意の文字列（＜ＳＡ−４＞参照）を作成するよう列挙された多型識別文字を並び替える。

すなわち、本ハプロタイプ推定装置は、一例として、図５に示すように、マーカー部位（Ｍ）毎に列挙された多型識別文字群を、まずコピー単位（ｕ）分、抽出し総当りの組合せを作ることにより、コピー単位（ｕ）の文字列を作成する（ＳＡ−３）。つぎに、コピー単位が最大値の回数繰り返される任意の文字列を作成するために、作成されたコピー単位（ｕ）の文字列を任意に組合せて、コピー単位（ｕ）の文字列の集合を作成する（ＳＡ−４の","で区切られた個々の文字列がコピー単位の文字列であり、｛ ｝で囲まれた文字列全体がコピー単位の文字列の集合である）。そして、本ハプロタイプ推定装置は、作成した文字列（コピー単位の文字列の集合）における、各多型塩基の種類を表す文字の文字数と、カウント数テーブルにおける、当該多型塩基の種類毎のカウント数と、が各マーカー部位において一致することを確認し、一致しない場合に当該文字列を除外する（ＳＡ−５）。

そして、本ハプロタイプ推定装置は、作成された文字列を、コピー単位で、任意の２つの文字列（「ハプロタイプ文字列」と呼ぶ。）に分離し、ハプロタイプ文字列の組合せとして格納する。例えば、図５＜ＳＡ−４＞最上段の文字列[F₁F₁F₂, F₁F₁F₀, F₁F₁F₂]を分離する場合、[/F₁F₁F₂, F₁F₁F₀, F₁F₁F₂][F₁F₁F₂/ F₁F₁F₀, F₁F₁F₂][F₁F₁F₂, F₁F₁F₀/ F₁F₁F₂][F₁F₁F₂, F₁F₁F₀, F₁F₁F₂/]の４通りのハプロタイプ文字列の組合せを作ることができる。ここで、例えば、[F₁F₃F₁/F₁F₃F₁, F₂F₃F₁] の"/"（スラッシュ）の両側の文字列は、[F₁F₃F₁]、[F₁F₃F₁, F₂F₃F₁]のように各個体における取り得る２つのハプロタイプを表現している。すなわち、２群１組で表された文字列（例えば[F₁F₃F₁/F₁F₃F₁, F₂F₃F₁]）は、１組のそれぞれのハプロタイプを表現すると同時にディプロタイプ（ハプロタイプ／ハプロタイプ）をも表現している。ここで、図６は、本発明により出力されるハプロタイプ表現形式の一例を示す図である。

図６に示すように、ハプロタイプ表現形式は、塩基多型とコピー数多型とを同時に表現する。すなわち、コピー数多型を考慮しない従来法でのハプロタイプ表現形式（図２参照）とは異なり、コピー数多型と塩基多型の両者を扱う本発明のハプロタイプ表現形式は、図６のようにコピー単位で表現される。

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

また、ここで、本ハプロタイプ推定装置は、コピー単位（ｕ）をマーカー部位（Ｍ）単位で複数の区画（パーティション）に分け、区画に含まれるマーカー部位の塩基多型について、上記の最大値の算出や、多型識別文字の列挙や、文字列の作成や、ハプロタイプ文字列の格納や、ハプロタイプの推定処理を、区画毎に実行させるよう制御してもよい。

この場合、本ハプロタイプ推定装置は、区画毎に算出されたハプロタイプ文字列について、コピー単位の繰り返しパターンが同じハプロタイプ文字列同士で、コピー単位における各区画を互いに連結（ライゲーション）することにより、複数の上記区画から構成されるハプロタイプ文字列を再現する。これにより、計算時間を短縮することができる。

また、本ハプロタイプ推定装置は、ハプロタイプの推定処理において、ＥＭ（Expectation-Maximization）法を用いて、集団におけるハプロタイプ文字列の頻度を、当該ハプロタイプ文字列を少なくとも一方に有する組合せの頻度により重み付けして算出するＭステップと、組合せの頻度を、当該組合せを構成するハプロタイプ文字列の頻度の積により求め、当該組合せの頻度に基づいて重みを算出するＥステップと、を頻度の値が収束するまで交互に繰り返すことにより所定の条件を満たす頻度を算出してもよい。ＥＭ法の詳細な処理については、後述する。また、本ハプロタイプ推定装置は、コピー数の総和が任意の一つの多型塩基のカウント数としたカウント数テーブルを用いてハプロタイプを推定する場合に、推定されたハプロタイプの組み合わせにおいて、ハプロタイプ文字列中の当該一つの多型塩基に対応付けられた多型識別文字の数を算出し、算出した数をハプロタイプにおけるコピー数として設定してもよい。以上で、本発明の概要の説明を終える。

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

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

記憶部１０６に格納される各種のデータベースやテーブル（カウント数テーブル１０６ａ〜コピー数総和ファイル１０６ｄ）は、固定ディスク装置等のストレージ手段であり、各種処理に用いる各種のプログラムやテーブルやファイルやデータベースやウェブページ等を格納する。

これら記憶部１０６の各構成要素のうち、カウント数テーブル１０６ａは、個体毎に、遺伝子型データを用いて、コピー数多型のコピー単位上の、標識によって特定されるマーカー部位に対応付けられた多型塩基をカウントすることにより得られたカウント数を、当該多型塩基の種類毎に記憶するカウント数テーブルである。上述したように、図４や図５＜ＳＡ−１＞は、カウント数テーブル１０６ａに格納されるカウント数データの一例を示す。このカウント数テーブル１０６ａに格納される情報は、図４に示すように、個体毎、マーカー部位毎、および、多型塩基の種類毎のカウント数を定義している。

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

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

また、コピー数総和ファイル１０６ｄは、各個体におけるコピー数多型のコピー数の総和（二本の相同染色体上のコピー単位の数の総和）を記憶するコピー数総和記憶手段である。

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

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

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

このうち、最大値算出部１０２ａは、カウント数テーブル１０６ａに記憶された個体のカウント数データに基づいて、マーカー部位（Ｍ）毎にカウント数（ｃ）の和を集計し、カウント数（ｃ）の和のうちの最大値を求める最大値算出手段である。

また、多型識別文字列挙部１０２ｂは、多型塩基の種類に対応付けられた多型識別文字（Ｆ₁、Ｆ₂、Ｆ₃など）を、当該多型塩基のカウント数（ｃ）列挙する多型識別文字列挙手段である。ここで、多型識別文字列挙部１０２ｂは、多型塩基の種類として塩基欠失を表す多型識別文字（一例として、Ｆ₀）を、最大値算出部１０２ａで算出された、最大値から和を減じた数だけ列挙してもよい。

また、文字列作成部１０２ｃは、コピー単位（ｕ）が最大値の回数繰り返される任意の文字列を作成するよう、多型識別文字列挙部１０２ｂによって列挙された多型識別文字を並び替える文字列作成手段である。ここで、文字列作成部１０２ｃは、作成した文字列における多型識別文字の文字数と、カウント数テーブル１０６ａにおける当該多型塩基の種類毎のカウント数（ｃ）と、が一致することを確認し、一致しない場合に当該文字列を除外してもよい。

また、ハプロタイプ文字列格納部１０２ｄは、文字列作成部１０２ｃによって作成された文字列を、任意のコピー単位で２分し、ハプロタイプ文字列の組合せとしてハプロタイプ文字列ファイル１０６ｂに格納するハプロタイプ文字列格納手段である。

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

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

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

ここで、Ｐ（ｈ_i）はハプロタイプ文字列の頻度を表し、ｈはハプロタイプ文字列を表し、ｉはハプロタイプ文字列のインデックスを表す。また、ｎは集団を構成する個体の数、ｊはカウント数のパターンのインデックス、ｋはハプロタイプ文字列の組合せのインデックス、Ｎ（ｃ_j）はカウント数のパターンｊを持つ個体の数を表す。また、δ(h_i,d_jk)は、組合せd_jkが一方に当該ハプロタイプ文字列h_iを有する場合に１を返し、両方に当該ハプロタイプ文字列h_iを有する場合に２を返し、当該ハプロタイプ文字列h_iを持たない場合に０を返す関数であり、ｄはハプロタイプ文字列の組合せを表す。また、ｗ_jkはハプロタイプ文字列の組合せの頻度による重みである。

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

再び図７に戻り、区画化部１０２ｇは、コピー単位（ｕ）をマーカー部位（Ｍ）の単位で複数の区画（パーティション）に分け、区画毎に含まれるマーカー部位の塩基多型について、最大値算出部１０２ａ、多型識別文字列挙部１０２ｂ、文字列作成部１０２ｃ、ハプロタイプ文字列格納部１０２ｄ、および、ハプロタイプ推定部１０２ｅに、処理を実行させるよう制御する区画化手段である。すなわち、区画化部１０２ｇは、ＰＬ（Partition-Ligation）法に基づくパーティション（Partition）ステップを実行する。

また、連結部１０２ｈは、区画毎にハプロタイプ推定部１０２ｅにより推定されたハプロタイプ文字列について、コピー単位（ｕ）の繰り返しパターンが同じハプロタイプ文字列同士で、コピー単位（ｕ）における各区画を互いに連結（ライゲーション）することにより、複数の区画から構成されるハプロタイプ文字列を再現する連結手段である。すなわち、連結部１０２ｈは、ＰＬ法に基づくライゲーション（Ligation）ステップを実行する。以上が本ハプロタイプ推定装置１００の内部構成である。

また、カウント数テーブル作成部１０２ｊは、遺伝子型データファイル１０６ｃに記憶された個体毎の遺伝子型データを用いて、マーカー部位に対応付けられた多型塩基をカウントし、カウント数を多型塩基の種類毎に、カウント数テーブル１０６ａに対応付けて格納するカウント数テーブル作成手段である。ここで、カウント数テーブル作成部１０２ｊは、コピー数総和ファイル１０６ｄに記憶された各個体におけるコピー数の総和を、任意の一つの多型塩基のカウント数としたカウント数テーブル１０６ａを作成してもよい。

また、コピー数設定部１０２ｋは、ハプロタイプ推定部１０２ｅにより推定されたハプロタイプの組み合わせにおいて、ハプロタイプ文字列中の任意の一つの多型塩基に対応付けられた多型識別文字の数を算出し、算出した数をハプロタイプにおけるコピー数として設定するコピー数設定手段である。

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

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

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

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

図８に示すように、最大値算出部１０２ａは、カウント数テーブル１０６ａに記憶された個体のカウント数（ｃ）データに基づいて、マーカー部位（Ｍ）毎にカウント数（ｃ）の和を集計し、カウント数（ｃ）の和のうちの最大値を求める（ＳＢ−１）。ここで、カウント数テーブル作成部１０２ｊは、予め、遺伝子型データファイル１０６ｃに記憶された各個体の遺伝子型データに基づいて、マーカー部位（Ｍ）に対応付けられた多型塩基をカウントし、カウント数（ｃ）を多型塩基の種類毎に、カウント数テーブル１０６ａに対応付けて格納してもよい。

そして、多型識別文字列挙部１０２ｂは、多型塩基の種類に対応付けられた多型識別文字（Ｆ₁、Ｆ₂、Ｆ₃など）を、当該多型塩基のカウント数（ｃ）ずつ列挙する（ＳＢ−２）。

そして、多型識別文字列挙部１０２ｂは、塩基欠失を表す多型識別文字（Ｆ₀）を、最大値算出部１０２ａで算出された最大値から和を減じた数だけ列挙する（ＳＢ−３）。このとき、列挙された多型識別文字は、コピー単位の数に最大値を乗じた数になる。

そして、文字列作成部１０２ｃは、コピー単位（ｕ）が最大値の回数繰り返される任意の文字列を作成するよう、多型識別文字列挙部１０２ｂによって列挙された多型識別文字を並び替える（ＳＢ−４）。

そして、文字列作成部１０２ｃは、作成した文字列における各多型識別文字の文字数と、カウント数テーブル１０６ａにおける当該多型塩基の種類毎のカウント数（ｃ）と、が一致することを確認し、一致しない場合に当該文字列を除外する（ＳＢ−５）。このとき、カウント数テーブル１０６ａには塩基欠失の多型の数のデータを記憶していないので、当然のことながら、塩基欠失を表す多型識別文字（Ｆ₀）の数との一致を確認しなくともよい。

そして、ハプロタイプ文字列格納部１０２ｄは、文字列作成部１０２ｃによって作成された文字列を、任意のコピー単位で２分し、ハプロタイプ文字列の組合せとしてハプロタイプ文字列ファイル１０６ｂに格納する（ＳＢ−６）。例えば、コピー単位を３つ持つ文字列では４つの組合せを算出することができ、コピー単位を４つもつ文字列では５つの組合せを算出することができ、コピー単位をＮ個もつ文字列では（Ｎ＋１）個の組合せを算出することができる。

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

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

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

［最大値算出処理〜ハプロタイプ文字列格納処理］
取り得るハプロタイプの組合せを求めるための、最大値算出処理〜ハプロタイプ文字列格納処理の一例について図９を参照しながら詳細に説明する。

ここで、取り得るハプロタイプの組合せを過不足なく求める方法、すなわち、カウント数データと矛盾しないディプロタイプ（ハプロタイプの組合せ）を見出す方法の一例としては、まず、カウント数データと矛盾しないコピー単位（ｕ）の集合を見出し（多型識別文字列挙部１０２ｂによる多型識別文字列挙手段に相当する）、つぎに、ディプロタイプ（ハプロタイプの組合せ）を構成する（ハプロタイプ文字列格納部１０２ｄによるハプロタイプ文字列格納手段に相当する）方法がある。ここで、カウント数データと矛盾しないコピー単位の集合とは、文字の集合（文字列）に属する全コピー単位（ｕ）に渡って、コピー単位（ｕ）の各マーカー部位（Ｍ）における各多型塩基のカウント数（ｃ）を総計したとき、その各マーカー部位（Ｍ）、各多型塩基の総カウント数（ｃ）が、カウント数テーブル１０６ａに記憶されたカウント数データにおける各マーカー部位、各塩基のカウント数と一致する、ということである。ここで、図９は、カウント数データの全パターンにおいて、そのカウント数パターンと矛盾しないディプロタイプを作成する処理の一例を示すフローチャートである。

図９に示すように、まず、文字列作成部１０２ｃは、カウント数テーブル１０６ａに記憶されたカウント数データから一意なカウント数パターン（ｃ_i）を算出する（ＳＣ−２１）。

つぎに、ハプロタイプ推定装置１００は、カウント数パターン（ｉ）のイテレーション（ｉイテレーション）に入る。ここでｉイテレーションは、最初ｉを１に初期化し、１イテレーション毎にｉを１ずつ増加させ、"i<=カウント数パターンの個数"である限り繰り返すこととする。ｉイテレーション内において、最大値算出部１０２ａは、カウント総数（和）s_j=Σ_kc_i(M_j, F_k) に対し、jに関する最大値を示すカウント総数（和）を見出し、これを最大コピー数ｚ_iとする（ＳＣ−２２）。

そして、ハプロタイプ推定装置１００は、マーカー部位（Ｍ）に関するイテレーション（ｊイテレーション）に入る（ＳＣ−２３）。ここで、ｊイテレーションは、最初ｊを１に初期化し、１イテレーション毎に１ずつ増加させ、"j<=マーカー部位の個数"である限り繰り返すこととする。ｊイテレーション内において、多型識別文字列挙部１０２ｂは、各マーカー部位Ｍ_jに対し、カウント数c_i(M_j, F_k)の個数分の多型識別文字Ｆ_kを全ｋに渡って列挙する。ここで、多型識別文字列挙部１０２ｂは、最大コピー数ｚ_iよりカウント総数s_jが少ない場合に、塩基欠損を表す多型識別文字Ｆ₀を不足分補う。そして、文字列作成部１０２ｃは、多型識別文字列挙部１０２ｂにより列挙された多型識別文字F_kの群から、あるマーカー部位（Ｍ）に含まれる多型識別文字の群を表すベクトルＡ_jを作成する（例えば、図６に示したＭ₁においては、Ａ₁=(F₁, F₁, F₂)である）。このように、ハプロタイプ推定装置１００は、ｊイテレーションを終了する。

次に、文字列作成部１０２ｃは、Ａ_jにある多型識別文字群を、各ｊに渡って１文字ずつそれぞれ総当たりで拾い、j=1, 2, 3, …の順番で文字列を、総当たり分作る（ＳＣ−２４）。例えば、ベクトルＡ₁=(a₁, a₂), Ａ₂=(b₁, b₂), Ａ₃=(c₁, c₂) である場合、a₁b₁c₁, a₁b₁c₂, a₁b₂c₁, a₁b₂c₂, a₂b₁c₁, a₂b₁c₂, a₂b₂c₁, a₂b₂c₂のように、総当たり分文字列を作る。ここで文字列作成部１０２ｃにより作成された文字列１つは、コピー単位１つに相当する。

そして、文字列作成部１０２ｃは、得られたコピー単位を一意化する。ここで、一意化とは、重複する余分なものを取り除くことである（ＳＣ−２５）。

そして、文字列作成部１０２ｃは、得られたコピー単位から、最大コピー数z_i個分コピー単位をあらゆる組合せで抽出して、コピー単位の集合をあらゆる組合せで作る（ＳＣ−２６）。これにより、コピー単位が最大値の回数（最大コピー数）繰り返された任意の文字列が得られる。

そして、文字列作成部１０２ｃは、得られた文字列（コピー単位の集合）を一意化する（ＳＣ−２７）。

そして、文字列作成部１０２ｃは、得られた文字列から、カウント数データと矛盾しない集合だけを、下記の数式を使って選択する（ＳＣ−２８）。すなわち、文字列作成部１０２ｃは、作成した文字列における多型識別文字の文字数と、カウント数テーブル１０６ａにおける当該多型塩基の種類毎のカウント数と、が一致することを下記の数式（あるいはそれと数学的に同値なもの）を用いて確認し、一致しない場合に当該文字列を除外する。

ここで、上記数式において、ｃ_iは、カウント数パターンを表す。ここでｉはパターンを指すインデックスであり、個体（個人）を指すインデックスではない（異なる個体（個人）が同じカウント数パターンを持っている場合は、同一のｃ_iで表される）。また、c_i(M_j, F_k) は、マーカー部位Ｍ_jにおける多型塩基Ｆ_kのカウント数を表す（例えば、上述した図４の個体１のパターンにおいては、c₁(M₁, F₁)=2, c₁(M₂, F₁)=1である）。ｕ_lはコピー単位を表し、ｕ_l(M_j)は、コピー単位のマーカー部位Ｍ_jにおける多型塩基（塩基欠失を意味するF₀を含む）を表すとする（例えば、上述した図３においては、u₁(M₁)=F₁, u₂(M₂)=F₃, u₃(M₃)=F₀, u₀(M₁)= F₀ である）。

そして、ハプロタイプ文字列格納部１０２ｄは、得られた各文字列を、任意のコピー単位で、２つに分けて、２つのハプロタイプ文字列を１組として全文字列について作成し、ハプロタイプ文字列ファイル１０６ｂに格納する（ＳＣ−２９）。例えば、"/"で２分するとして、[F₁F₁F₂, F₁F₁F₀, F₁F₂F₂] に関しては、[/F₁F₁F₂, F₁F₁F₀, F₁F₂F₂]，[F₁F₁F₂/ F₁F₁F₀, F₁F₂F₂], [F₁F₁F₂, F₁F₁F₀/ F₁F₂F₂], [F₁F₁F₂, F₁F₂F₂, F₁F₁F₀/] の４つの組を作る。ここで得られた２ハプロタイプ文字列１組が１ディプロタイプに相当する。ここで、ハプロタイプ文字列格納部１０２ｄは、特殊な組合せとして、全ての多型識別文字がF₀である文字列と全コピー単位との組を、追加してもよい（一例として、[F₀F₀F₀/ F₁F₁F₂, F₁F₁F₀, F₁F₂F₂]）。

ハプロタイプ推定装置１００は、上記の処理を各カウント数パターンに対して行い、ｉイテレーションを終え、最終的に、カウント数データと矛盾しない全ディプロタイプを算出する。以上が、取り得るハプロタイプの組合せを求めるための、最大値算出処理〜ハプロタイプ文字列格納処理の一例についての詳細な説明である。上記の最大値算出処理〜ハプロタイプ文字列格納処理により得られたハプロタイプ文字列の組合せを用いて、ハプロタイプ推定部１０２ｅは、集団におけるハプロタイプ文字列の頻度を求めることとなる。

［ＥＭ法による処理］
ハプロタイプ推定部１０２ｅの処理によるＥＭ法の詳細な処理の一例について説明する。すなわち、ハプロタイプ推定部１０２ｅは、ハプロタイプ頻度算出部１０２ｆの処理により、集団におけるハプロタイプ文字列の頻度を、ＥＭ（Expectation-Maximization）法を用いて、効率よく計算する。

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

また、Ｐ（ｈ_i）はハプロタイプ文字列の頻度を表し、ｈはハプロタイプ文字列を表し、ｉはハプロタイプ文字列のインデックスを表す。また、ｎは集団を構成する個体の数、ｊはカウント数のパターンのインデックス、ｋはハプロタイプ文字列の組合せのインデックス、Ｎ（ｃ_j）はカウント数のパターンｊを持つ個体の数を表す。また、δ(h_i,d_jk)は、組合せd_jkが一方に当該ハプロタイプ文字列h_iを有する場合に１を返し、両方に当該ハプロタイプ文字列h_iを有する場合に２を返し、当該ハプロタイプ文字列h_iを持たない場合に０を返す関数であり、ｄはハプロタイプ文字列の組合せを表す。また、ｗ_jkはハプロタイプ文字列の組合せの頻度による重みである。

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

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

［ＰＬ法による処理］
ＰＬ法による処理の一例について図１０〜図１２を参照しながら説明する。

最大値算出部１０２ａ〜ハプロタイプ推定部１０２ｅの処理により、カウント数テーブル１０６ａに記憶されたカウント数データから、ハプロタイプ文字列の組合せとその頻度を合理的理由の下で推定できるが、マーカー部位（Ｍ）の数が増えると、取り得るハプロタイプの数が飛躍的に増大するため、計算時間が多大になる。このような場合、計算時間を節約する方法として、ＰＬ（Partition-Ligation）法が用いられる。ここで、図１０は、コピー数多型を考慮しないハプロタイプ推定方法のＰＬ法におけるligation の一例を示す図である。

図１０に示すように、ＰＬ法は、複数の座位（L）を、より少数の区画（図１０において、Partitioned Loci）に分割し（パーティション（partition）ステップ）、その区画（パーティション）内で区画毎にハプロタイプの組合せとその頻度を推定し、頻度が低いハプロタイプを取り除く（頻度以外の別の基準でハプロタイプを復帰させることもある）。そして、残ったハプロタイプに関し、区画（一般に隣り合った区画）間でハプロタイプをつなぎ（ライゲーション（ligation）ステップ、図１０においてLigated loci）、つなぎ合わされたハプロタイプに対しつなぎ合わされた区画の遺伝子型データを使って、つなぎ合わされたハプロタイプの頻度の推定を行う、という方法である。ここで「ハプロタイプ同士をつなげる」とは、ハプロタイプを文字列（本実施の形態における「ハプロタイプ文字列」）と見たとき、文字列同士をつないでより長い文字列を作成することに相当する。なお、このＰＬ法を、直接、本実施の形態のコピー数多型と塩基多型を含む遺伝子型データに適用することはできない。すなわち、ligation ステップにおいて、座位を隣り合う座位でつなげることは、本実施の形態においては、各マーカー部位がコピー単位で区切られているので、つなぎ合わせることが出来ないからである。ここで、図１１は、本実施の形態の、コピー数多型を考慮に入れたハプロタイプ推定方法のＰＬ法における連結（ligation）処理の一例を示す図である。ここで、ｈ₀はコピー単位の数が０の欠失ハプロタイプを表す。

図１１に示すように、本実施の形態では、区画化部１０２ｇは、コピー単位（ｕ）をマーカー部位（Ｍ）の単位で複数の区画（パーティション）に分け（Partitioned markers）、区画に含まれるマーカー部位（Ｍ）の塩基多型について、最大値算出部１０２ａ〜ハプロタイプ頻度算出部１０２ｆに対し当該区画毎に処理を実行させるよう制御する（区画化ステップ）。そして、連結部１０２ｈは、区画毎に算出されたハプロタイプ文字列について、コピー単位の繰り返しパターンが同じハプロタイプ文字列同士で（例えば、図１１において、ｈ１とｈ４同士、ｈ２とｈ５同士）、コピー単位における各区画を互いに連結（ライゲーション）させることにより、複数の区画から構成されるハプロタイプ文字列を再現する（連結ステップ）。ここで、図１２は、連結ステップにおける合理的な連結の方法を示す図である。

すなわち、図１２に示すように、連結部１０２ｈは、異なる区画間でハプロタイプ文字列の連結ステップを行う際、同数のコピー単位（ｕ）を持つハプロタイプ文字列に対し、コピー単位同士を余すことなくつなげ、新しいハプロタイプを作成する。ここで「コピー単位同士をつなげる」とは、コピー単位を文字列と見たとき、文字列同士をつないでより長い文字列を作成することに相当する。連結部１０２ｈは、ハプロタイプが持つコピー単位同士を余すことなくつなげる処理を、コピー単位同士のあらゆる組合せで行う（例えば、図１２における染色体３を参照）。これには合理的理由があり、区画化ステップにおいて複数のマーカー部位（Ｍ）の単位で分割するとは、上述のコピー数多型を考慮しないＰＬ法におけるようなハプロタイプを分割するというのではなく、本実施の形態においてはコピー単位を分割することに相当するからである。

［実施の形態２］
各個体におけるコピー数多型のコピー数の総和を用いた実施の形態２について、図１３および図１４を参照して説明する。

本実施の形態２においては、実施の形態１のように蛍光色素プローブ等の標識によって識別される多型塩基のカウント数データからハプロタイプを推定するだけではなく、何らかの実験手段によって得られるコピー数多型のコピー数の総和（二つの相同染色体上に渡る総和）データからもハプロタイプを推定する。すなわち、本ハプロタイプ推定装置１００は、コピー数総和ファイル１０６ｄに記憶された各個体におけるコピー数多型のコピー数の総和データを用いて、ハプロタイプに存在するコピー数とそのハプロタイプの頻度を推定する。なお、本実施の形態２においては、ハプロタイプ推定の対象であるコピー単位上には必ずしも多型塩基が存在しなくてもよい。ここで、図１３は、本実施の形態２における処理の一例を示すフローチャートである。

図１３に示すように、コピー数総和ファイル１０６ｄは、各個体におけるコピー数多型のコピー数の総和を記憶する（ＳＤ−１）。ここで、「コピー数の総和」とは、二本の相同染色体にわたるコピー単位の数の総和のことである。

つぎに、カウント数テーブル作成部１０２ｊは、コピー数総和ファイル１０６ｄに記憶された各個体におけるコピー数の総和を、任意の一つの多型塩基のカウント数としたカウント数テーブル１０６ａを作成する（ＳＤ−２）。ここで、図１４は、各個体におけるコピー数の総和データと、作成されたカウント数テーブルの対応関係の一例を示す図である。

図１４左図に示すように、カウント数テーブル作成部１０２ｊは、コピー数総和ファイル１０６ｄに記憶されたあるコピー数多型（ＣＮＰ）領域における各個体１，２，３・・・のコピー数の総和データに基づいて、図１４右図に示すように、コピー数の総和２，３，０・・・を、暫定的にマーカー部位（Ｍ_１）における１種類の多型塩基（Ｆ_１）のカウント数として扱い、カウント数テーブル１０６ａを作成する。すなわち、図１４の例では、カウント数テーブル作成部１０２ｊは、個体１のコピー数の総和が「２」、個体２のコピー数の総和が「３」であるので、個体１についてはマーカー部位Ｍ_１における仮の多型塩基Ｆ_１のカウント数「２」とし、個体２についてはマーカー部位Ｍ_１における仮の多型塩基Ｆ₁のカウント数「３」としたカウント数テーブル１０６ａが作成されている。

そして、ハプロタイプ推定装置１００は、作成されたカウント数テーブル１０６ａに基づいて、実施の形態１において図８を用いて説明したハプロタイプ推定処理（ＳＢ−１〜９）を実行する（ＳＤ−３）。例えば、ハプロタイプ推定装置１００は、図１４の例では、個体１については、ハプロタイプ文字列の組合せとして[F₀/F₁, F₁]と[F₁/F₁]を作成し、個体２については、ハプロタイプ文字列の組合せとして[F₀/F₁, F₁, F₁]と[F₁/F₁, F₁] を作成し、ハプロタイプ推定処理を行う。ここで、Ｆ_１は、上述の任意の一つの多型塩基に対応付けられた多型識別文字であり、Ｆ_０は、コピー数が０であることを表す多型識別文字である。

そして、コピー数設定部１０２ｋは、ハプロタイプ推定部１０２ｅにより推定されたハプロタイプの組み合わせにおいて、ハプロタイプ文字列中の上述の任意の一つの多型塩基に対応付けられた多型識別文字（Ｆ_１）の数を算出し、算出した数をハプロタイプにおけるコピー数として設定する（ＳＤ−４）。例えば、コピー数設定部１０２ｋは、ハプロタイプ推定部１０２ｅの処理により推定されたハプロタイプの組み合わせにおいて、ハプロタイプ文字列[F₀]、[F₁]、[F₁, F₁]、[F₁, F₁, F₁]における多型塩基Ｆ_１を数え、それぞれコピー数「0」,「1」,「2」,「3」として設定する。これにより、目的であるハプロタイプに存在するコピー数（コピー単位の数）とそのハプロタイプの頻度を取得することができる。すなわち、例えば、コピー数が0, 1, 2, 3 であるコピー数多型のハプロタイプとその頻度が得られる。以上で、本実施の形態２の説明を終える。

［実施例］
上述の実施の形態をプログラミング言語Perlにて実装した実施例について、以下に図１５〜図３０を参照しながら説明する。

［ＥＭ法の実施例］
本実施例におけるＥＭ法によるハプロタイプ推定処理の一例について説明する。ここで、図１５は、ＥＭ法によるハプロタイプ推定処理の一例を示すフローチャートである。なお、以下の説明において、本実施の形態におけるハプロタイプ文字列を単に「ハプロタイプ」と、ハプロタイプ文字列の組合せを「ディプロタイプ」と述べる場合がある。

図１５に示すように、まず、最大値算出部１０２ａは、カウント数テーブル１０６ａを参照してカウント数データを読み込む（ＳＣ−１）。このとき、ハプロタイプ頻度算出部１０２ｆの処理による収束判定の数値を読み込んでもよい。

次に、ハプロタイプ推定装置１００は、［最大値算出処理〜ハプロタイプ文字列格納処理］の項で図９を用いて説明した処理（ＳＣ−２１〜ＳＣ−２９）を行う（ＳＣ−２）。

そして、以降のステップにおいて、ハプロタイプ推定部１０２ｅは、ハプロタイプ頻度算出部１０２ｆの処理により、ＥＭ法を用いて、ハプロタイプ文字列ファイル１０６ｂに記憶されたハプロタイプ文字列の組合せとカウント数テーブル１０６ａに記憶されたカウント数データに基づいて、ハプロタイプ文字列の頻度を計算する。

すなわち、まず、ハプロタイプ頻度算出部１０２ｆは、各カウント数パターン（ｃ_i）に対し、各ハプロタイプの組合せの存在の重みを初期化する（ＳＣ−３）。本実施例においては、初期値として重みを平等に割り付けた。すなわち、初期値の重みは、w_jk=1/n_jkである。ここでw, j, k は（発明の構成における）上述と同じであり（wは重み、jはカウント数パターン、kはカウント数パターン内でのディプロタイプ（ハプロタイプ文字列の組合せ）のインデックス）、n_jkはカウント数パターン内でのディプロタイプの総数である。

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

そして、ハプロタイプ推定部１０２ｅは、上記条件によってＳＣ−５に処理を移した場合、ハプロタイプ及びその頻度を結果ファイルに出力する。

[ＥＭ法の実証シミュレーション]
本実施の形態により塩基多型を考慮したコピー数多型（ＣＮＰ）に関するカウント数データから、コピー数多型のハプロタイプとその頻度を推定できるかどうかを確かめる為、シミュレーション実験を行った。すなわち、以下のシミュレーション実験においては、ハプロタイプ推定装置１００は、ハプロタイプ推定処理とは逆の処理を行い、結果（ハプロタイプと確率（頻度））からカウント数データを導き出し、そのカウント数データを用いて、ハプロタイプ推定処理を行い元の結果を再現できるか否かを実証する。ここで、図１６は、シミュレーションの枠組みを示すフローチャートである。

図１６に示すように、ハプロタイプ推定装置１００は、ＣＮＰのハプロタイプとその確率（頻度）が書かれたファイルを読み込む（ＳＥ−１）。これは、例えば図１７に示すようなデータである（すなわち、本来、ハプロタイプ推定の結果を表すデータである）。また、ハプロタイプ推定装置１００は、以降の処理で用いる個体の数も読み込む。

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

そして、ハプロタイプ推定装置１００は、ディプロタイプとその確率から、ディプロタイプの多項分布を構成し、与えられた個体の数だけディプロタイプをランダムに抽出する（ＳＥ−３）。これは、例えば、Ｒ言語を使って、rmultinom(1, size=個体の数, prob=c(ディプロタイプ１の確率、ディプロタイプ２の確率、ディプロタイプ３の確率、…))のようなコマンドで簡単に実施できる。ここにおいて抽出されたディプロタイプの１つ１つが、１つ１つの個体（のディプロタイプ）に相当する。

そして、ハプロタイプ推定装置１００は、各個体が持つディプロタイプからそれを構成するコピー単位（ｕ）を全て取り出す（ＳＥ−４）。すなわち、ハプロタイプ推定装置１００は、それらコピー単位（ｕ）の各マーカー部位（Ｍ）における各Ｆ_k（k=1, 2, 3, …）の（全コピー単位に渡る）カウント総数を数え、カウント数データを作成する。例えば、ハプロタイプ推定装置１００は、個体の持つディプロタイプが[F₂F₁/F₂F₂, F₂F₂]である場合に、全コピー単位[F₂F₁, F₂F₂, F₂F₂]を取り出し、各M₁, M₂ において、各F₁, F₂ のカウント総数を数え、この場合、M₁ におけるF₁, F₂ の総数はそれぞれ0, 3となり、M₂ におけるF₁, F₂ の総数はそれぞれ1, 2となる。

最後に、ハプロタイプ推定装置１００は、カウント数データをファイル（カウント数テーブル１０６ａ）に出力する（ＳＥ−５）。

以上説明したシミュレーションの枠組みにおいて、ハプロタイプ推定装置１００は、図１７で示したハプロタイプとその確率と同様のデータを、さらに個体の数５００を読み込み、カウント数データ（図１８）を作成した。

つづいて、本実施の形態のハプロタイプ推定装置１００が、出力した図１８のカウント数データだけから、図１７で示された元のハプロタイプとその確率を推定（再現）できるか試験した。これは言い換えれば、不完全な観測データからの母集団比率の推定問題である。本実施の形態を適用する際、ＥＭ法において収束判定に使われる対数尤度差の閾値は0.001未満とした。

そして本実施の形態を適用した結果、ハプロタイプ推定装置１００は、ハプロタイプとその頻度に関し、図１９に示す結果を出力した。図１９（再現データ）にあって図１７（元のデータ）にないハプロタイプの頻度は全て低い。すなわち、これらは存在しないと解釈される。逆に、図１７にあるハプロタイプは全て図１９に現れており、かつ、それら推定頻度も正解頻度とほぼ等しい。よって、本実施の形態はＣＮＰのハプロタイプとその頻度を精度よく推定することが実証された。

[ＰＬ法の実施例]
本実施の形態のＰＬ法を、階層型タイプのハプロタイプバッファ付ＰＬ法として、プログラミング言語Perlで実装した実施例について、以下に図２０〜図２２を参照しながら説明する。図２０は、本実装の枠組みを示すフローチャートである。

図２０に示すように、まず、区画化部１０２ｇは、カウント数データを読み込む（ＳＦ−１）。ここで、区画化部１０２ｇは、ＰＬ法において使われる区画（パーティション）内のマーカー数、区画内において取り除くべきハプロタイプの頻度（の閾値）、ハプロタイプ・バッファサイズも読み込む。また、区画化部１０２ｇは、ＥＭ法の実行に必要な上述のパラメータも読み込む。

そして、区画化部１０２ｇは、読込んだカウント数データを、読み込んだ区画内のマーカー部位（Ｍ）単位の個数ごとに、左側からｘ個の区画に分割する（ＳＦ−２）。この際、割り切れなかった最右端の区画は、この個数以下のマーカー部位（Ｍ）で構成される。

そして、区画化部１０２ｇは、イテレーション（ｉを１に初期化、１イテレーション毎にｉを１ずつ増加させ、"i<=x"である限り繰り返す）に入る。イテレーション内で、区画化部１０２ｇは、各区画内でのカウント数データに対し、ハプロタイプ頻度算出部１０２ｆの処理による上述のＥＭ法によってＣＮＰのハプロタイプとその頻度を計算する（ＳＦ−３）。なお、以降のＰＬ法による処理のため、ハプロタイプ頻度算出部１０２ｆの処理実行時に、各区画における各ハプロタイプの累積平均頻度を記録しておく。ここで、累積平均頻度とは、ＥＭ法での各イテレーションにおけるハプロタイプ頻度の、イテレーション開始時から収束時までの（イテレーションに渡る）平均値をいう。

そして、ハプロタイプ推定装置１００は、連結部１０２ｈの処理による連結（ライゲーション）ステップに入る。これはＳＦ−４、５、６とそれを取り囲むイテレーションである。まず、イテレーションの設定について説明する。

最初のイテレーション（図２０におけるｉイテレーション）は、階層型タイプのＰＬ法における、層のイテレーションを指す。このｉイテレーションにおいては、各層（即ち任意のｉ）において、次に深いイテレーション（図２０におけるｊイテレーション）での操作を、最終的に区画の数が１となるまで繰り返す。連結部１０２ｈは、次に深いイテレーション（ｊイテレーション）において、層にある区画を順次、連結（ライゲーション）する。連結し終えると、新たな区画の数はroundup(k/2)となる（ここで、ｋは連結する前の区画数であり、roundup()は切り上げを意味する）。すなわち、区画数は、連結し終えると、前の区画数が偶数ならば、連結し終えた後の区画数は前の1/2個となり、奇数ならば前の1/2+1個となる。このように、連結し終えると、区画の数は減る。すなわち、このｉイテレーションでは、層において区画を順次、連結し、連結し終えると、区画の数は減り、次の層へ行く。次の層でも区画を順次、連結し、連結し終えるとまた区画の数は減り、次の層へ行く。連結部１０２ｈは、この処理を、最終的に区画の数が１になるまで繰り返す。

次に深いイテレーション（図２０のｊイテレーション）では、前述のように層にある全ての区画を連結する。連結は区画ｊと区画ｊ＋１の隣接する区画間で行う。ここで、連結部１０２ｈは、最初の層では、ｊを１から始め、イテレーション毎順次２ずつ増やし、j<iである限り繰り返す。一方、次の層では、連結部１０２ｈは、ｊを、層における全区画（全パーティション）の数−１（i-1）から始め、順次２ずつ減らし、j>=1である限り繰り返す。このように、最初の層では左から、次の層では右から、その次の層では左から、さらにその次の層では右から、...、と層が進む度（ｉイテレーションが繰り返される度）、連結の開始を左、右と交互に設定し、上記のようにｊイテレーションを繰り返す。なお、区画数が奇数個の場合は、最後の区画は連結をせず、そのまま、次の層に持ち越す。

つぎに、イテレーション内部の処理の詳細について以下に説明する。

連結部１０２ｈは、区画ｊと区画ｊ＋１の両方の区画において、頻度が低いハプロタイプの除去と別条件によるハプロタイプの復帰を行う（ＳＦ−４）。別条件によるハプロタイプの復帰とは、除去後残ったハプロタイプの数が、与えられたハプロタイプ・バッファサイズ未満であった場合、ＥＭ法におけるハプロタイプの累積平均頻度が上位のものから順に、除去したハプロタイプをそのバッファ数まで復帰させる処理である。

そして、連結部１０２ｈは、区画ｊと区画ｊ＋１を合わせたより大きな区画に渡ってハプロタイプの連結処理を実行し、そのハプロタイプから構成され、かつカウント数データと矛盾しないディプロタイプ（ハプロタイプの組合せ）を作成する（ＳＦ−５）。ここで、図２１は、ＳＦ−５の、より詳しい処理の一例を示すフローチャートである。

図２１に示すように、まず、連結部１０２ｈは、区画ｊと区画ｊ＋１における全ハプロタイプのコピー数を数え、最大のコピー数を出す（ＳＦ−５１）。そして、連結部１０２ｈは、コピー数ｋを１にセットし、この最大のコピー数に達するまで、ｋイテレーションを実行する。このイテレーションの中のＳＦ−５１において、連結部１０２ｈは、区画ｊと区画ｊ＋１にあるコピー数ｋのハプロタイプを選択する（両方の区画からハプロタイプが得られない場合は、すぐ次のｋに行く）。そして、連結部１０２ｈは、それぞれの区画においてハプロタイプの個数を数え、ｌイテレーションとｍイテレーションを、その個数分、実行する。ここで、ｌイテレーションは、最初ｌを１に初期化し、１イテレーション毎にｌを１ずつ増加させ、"l<=区画ｊにおけるコピー数ｋのハプロタイプの数"である限り繰り返すこととする。また、ｍイテレーションは、最初ｍを１に初期化し、１イテレーション毎にｌを１ずつ増加させ、"m<=区画ｊ＋１におけるコピー数ｋのハプロタイプの数"である限り繰り返すこととする。このイテレーションの中で、実施の形態において上述したように、同数のコピー単位を持つハプロタイプに対し、コピー単位同士を余すことなく連結し、新しいハプロタイプをあらゆる組合せで作成する。

すなわち、連結部１０２ｈは、ハプロタイプｍをコピー単位という要素を持つベクトルと見なし、順列ベクトルを全て出す（ＳＦ−５２）。ここで、順列ベクトルとは、要素を順列的に入れ替えたベクトルのことである。順列ベクトルを全て出すのは、あらゆる組合せで新しいハプロタイプを作り出すためである。そして、連結部１０２ｈは、ハプロタイプｌをコピー単位という要素を持つベクトルと見なし、ハプロタイプｌのコピー単位と、上で出した各順列ベクトルにおけるコピー単位との、ベクトル的な文字列結合を行って、新しい文字列ベクトルを出す（ＳＦ−５２）。ベクトル的な文字列結合とは、実施の形態において上述したように、コピー単位を文字列と見なし、各要素独立に文字列を結合することである。こうして得られた文字列ベクトルを、次に（順序関係のない）文字列集合と見なす（本実施の形態におけるハプロタイプにおいてはコピー単位内に順序関係はない）。こうして、得られた文字列集合が、連結されたハプロタイプになる（ＳＦ−５３）。そして、ｍ，ｌイテレーションを終える。このようにして、それぞれのコピー数ごとに、コピー単位同士を余すことなくライゲートしたハプロタイプを、あらゆる組合せで作成する。

ｋイテレーションを終了すると、連結部１０２ｈは、連結されたハプロタイプを一意化する（重複する余分なハプロタイプを取り除く）（ＳＦ−５４）。

そして、ハプロタイプ推定装置１００は、各カウント数パターンに対するイテレーションに入る。すなわち、連結部１０２ｈは、連結されたハプロタイプから、２つ１組のディプロタイプをあらゆる組合せで作成する（ＳＦ−５５）。ただし、２つのハプロタイプのコピー数を足すと、カウント数パターンｃ_kの総コピー数となるように作成する。総コピー数とは、カウント総数Σnc_k(M_p, F_n) (n=1, 2, 3, …; p=1, 2, 3, …)の（p に関する）最大のカウント総数のことであり、区画化前（即ちＳＦ−２の処理前）に最大値算出部１０２ａにより算出されているものである。例えば、総コピー数が３の場合、コピー数が０と３の２つのハプロタイプ、あるいはコピー数が１と２の２つのハプロタイプから、ディプロタイプを作成する。そして、連結部１０２ｈは、作成された全ディプロタイプから、区画ｊと区画ｊ＋１を合わせた区画内でのカウント数データと、矛盾しないディプロタイプのみを選択する（ＳＦ−５６）。

再び図２０の説明に戻り、連結部１０２ｈは、区画ｊと区画ｊ＋１を合わせた区画において、ハプロタイプ頻度算出部１０２ｆの処理によるＥＭ法によるハプロタイプ頻度を計算させる（ＳＦ−６）。すなわち、ハプロタイプ頻度算出部１０２ｆは、区画ｊと区画ｊ＋１を合わせた区画におけるカウント数データ、ＳＦ−５で得られたディプロタイプを用い、ＭステップとＥステップを実行してハプロタイプ頻度を計算する。このとき、ハプロタイプ頻度算出部１０２ｆは、各ハプロタイプの累積平均頻度を記録する。

そして、ハプロタイプ推定装置１００は、ｉイテレーションを終了した後は、ハプロタイプ及びその頻度を結果ファイルに出力する（ＳＦ−７）。

［ＰＬ法の実証シミュレーション]
以上の本実施例により、区画化部１０２ｇおよび連結部１０２ｈの処理によるＰＬ法により同一の結果（コピー数多型のハプロタイプとその頻度）を推定でき、且つ計算処理の短縮を図ることができるか否かを確かめる為、シミュレーション実験を行った。すなわち、以下のシミュレーション実験においては、ハプロタイプ推定装置１００は、区画化処理および連結処理を行う場合と行わない場合の２通りのハプロタイプ推定処理を行い、区画化・連結処理を行わない場合と同様のハプロタイプ推定結果を再現でき、計算速度を速めることができることを実証する。図２２は、シミュレーション実験に用いたハプロタイプとその頻度のデータの一例を示す図である。図２３は、読込まれたハプロタイプとその確率についての、作成された個体の数５００のカウント数データである。

本実施例により、図２３のカウント数データだけから、図２２で示されたハプロタイプとその確率を短時間で推定できるか、ＰＬ法によって試験した。比較のため、ＰＬ法を用いないブルートフォース法（上述のＥＭ法の実施例と同じ方法）でも行った。ブルートフォース法及びＰＬ法で使うＥＭ法の設定は、上述の実施例の設定と同じとした。ＰＬ法において、マーカー部位（Ｍ）は、隣り合う４マーカー部位毎に分割した。区画内において取り除くべきハプロタイプの頻度は、10^-5 以下とした。ハプロタイプ・バッファサイズは５０とした。

本実施例を適用した結果、図２４のような結果を得た。この結果において、ブルートフォース法とＰＬ法の、計算時間とハプロタイプの推定頻度を示す。ブルートフォース法とＰＬ法いずれも推定頻度はほぼ同じであって、それは図２２の正解に近い。即ち、図２４にあって図２２にないハプロタイプの頻度は全て低く、図２２にあるハプロタイプは全て図２４に現れており、それら推定頻度は正解頻度とほぼ等しい。そして、ＰＬ法における計算時間はブルートフォース法より大幅に短い。これより、本発明のＰＬ法によって、ＣＮＰのハプロタイプとその頻度をより短い時間で推定できることが実証された。

［コピー単位の数による実施形態２の実証シミュレーション］
各個体におけるコピー数多型のコピー数の総和を用いた実施の形態２が適用される実施例について、図２５〜図３０を参照して説明する。

上述の実施の形態２を適用した本実施例では、コピー数多型のコピー数の総和を用いたカウント数データから、コピー数多型のハプロタイプとその頻度を推定できるかどうかを確かめる為、シミュレーション実験を行った。すなわち、本実施例では、予め、設定した集団におけるコピー数多型のハプロタイプの確率（頻度）に基づいて、個体におけるコピー数の総和データを作成する。そして、このコピー数総和データに基づいて、実施の形態２の処理を実行し、元のハプロタイプの確率（頻度）を推定できるか試験を行った。ここで、図２５は、本実施例において設定した、ハプロタイプに存在するコピー数（コピー単位の数）とそのハプロタイプの確率を示す図である。また、図２６は、図２５に示すコピー数を仮の塩基Ｆ_１のカウント数として表現した図である。

すなわち、シミュレーション実験の設定として、図２５に示すように、ハプロタイプのコピー数が０，１，２，３である場合の確率をそれぞれ０．２５として設定した。なお、このコピー数は、図２６に示すように、実施の形態２における処理において、仮の塩基Ｆ_１のカウント数として扱われる。

つぎに、当該設定した確率に基づいて、ハプロタイプのコピー数の組合せを任意に５００個体の抽出し、コピー数の総和（二本の相同染色体にわたるコピー数の総和）をコピー数総和ファイル１０６ｄに格納した。ここで、図２７は、抽出された各個体におけるコピー数の総和を示す図である。また、図２８は、図２７に示すコピー数の総和を仮の塩基Ｆ_１のカウント数として表現した図である。

つづいて、図２７に示すコピー数の総和データを用いて、ハプロタイプ推定装置１００が図２５に示したようなハプロタイプの確率（頻度）を推定できるか試験を行った。すなわち、まず、カウント数テーブル作成部１０２ｊは、図２７に示すようにコピー数総和ファイル１０６ｄに格納されたコピー数の総和を、図２８に示すように仮の塩基Ｆ_１のカウント数としたカウント数テーブルを作成した。

そして、ハプロタイプ推定装置１００は、以降のハプロタイプ推定処理を実行し、ハプロタイプとその頻度を推定した。なお、ハプロタイプ推定処理実行の際、ＥＭ法の設定は、上述の設定と同じとした。ここで、図２９は、ハプロタイプ推定装置１００により推定されたハプロタイプ文字列とその頻度を示す図である。また、図３０は、推定されたハプロタイプ文字列から算出された各コピー数のハプロタイプとその頻度を示す図である。

そこで、コピー数設定部１０２ｋは、図２９のように推定されたハプロタイプにおいて、ハプロタイプ文字列中の仮の塩基に対応付けられた多型識別文字（Ｆ_１）の数を算出し、図３０に示すように算出した数をハプロタイプにおけるコピー数として設定した。

図２５で示したハプロタイプの設定確率と図３０に示す推定結果であるハプロタイプの推定頻度を比較すると、コピー数が０，１，２，３であるハプロタイプの設定確率０．２５とほぼ同等のハプロタイプ推定頻度（約０．２４５〜０．２６０）が得られたのに対し、設定確率０としたコピー数４、５、６のハプロタイプ推定頻度は、ほぼ０に近い値（１．０５Ｅ−０５以下）を示した。すなわち、推定頻度は、正解頻度（設定確率）とほぼ等しく、コピー数の総和データだけから、ハプロタイプに存在するコピー単位の数とそのハプロタイプの頻度を推定できることが実証された。以上で、コピー単位の数による実施形態２の実証シミュレーションの説明を終える。

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

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

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

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

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

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

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

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

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

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

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

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

以上詳述に説明したように、本発明によれば、塩基多型を考慮したコピー数多型に関する実験データから、ハプロタイプとその頻度を高精度で推定することができる、ハプロタイプ推定装置、ハプロタイプ推定方法、および、プログラムを提供することができ、バイオテクノロジー分野や医療分野において利用可能である。

Claims (9)

1. 集団における各個体のコピー数多型と塩基多型を含む遺伝子型データからハプロタイプを推定する、制御部と記憶部を少なくとも備えたハプロタイプ推定装置において、
   上記記憶部は、
   上記個体毎に、上記遺伝子型データを用いて、上記コピー数多型のコピー単位上の、標識によって特定されるマーカー部位に対応付けられた多型塩基をカウントすることにより得られたカウント数を、当該多型塩基の種類毎に記憶するカウント数テーブル、
   を備え、
   上記制御部は、
   上記カウント数テーブルに格納された上記個体の上記カウント数に基づいて、上記マーカー部位毎に上記カウント数の和を集計し、上記カウント数の和のうちの最大値を求める最大値算出手段と、
   上記多型塩基の種類に対応付けられた多型識別文字を、当該多型塩基の上記カウント数列挙する多型識別文字列挙手段と、
   上記コピー単位が上記最大値の回数繰り返される任意の文字列を作成するよう、列挙された上記多型識別文字を並び替える文字列作成手段と、
   上記文字列作成手段によって作成された上記文字列を、任意の上記コピー単位で２分し、ハプロタイプ文字列の組合せとして格納するハプロタイプ文字列格納手段と、
   上記集団において、同一である上記ハプロタイプ文字列の数を集計し、当該ハプロタイプ文字列の上記集団における頻度を求め、当該頻度が所定の条件を満たす上記各個体の上記ハプロタイプ文字列の上記組合せを、上記ハプロタイプの組合せとして推定するハプロタイプ推定手段と、
   を備えたことを特徴とするハプロタイプ推定装置。
2. 請求項１に記載のハプロタイプ推定装置において、
   上記ハプロタイプ推定手段は、
   上記ハプロタイプ文字列の頻度を、ハーディ・ワインバーグの法則に基づいて算出し、
   上記所定の条件を、上記集団におけるハーディ・ワインバーグ平衡とすること、
   を特徴とするハプロタイプ推定装置。
3. 請求項１または２に記載のハプロタイプ推定装置において、
   上記多型識別文字列挙手段は、
   上記多型塩基の種類として塩基欠失を表す上記多型識別文字を、上記最大値から上記和を減じた数だけ列挙すること、
   を特徴とするハプロタイプ推定装置。
4. 請求項１乃至３のいずれか一つに記載のハプロタイプ推定装置において、
   上記文字列作成手段は、
   作成した上記文字列における、上記多型識別文字の文字数と、上記カウント数テーブルにおける、当該多型塩基の種類毎の上記カウント数と、が一致することを確認し、一致しない場合に当該文字列を除外すること、
   を特徴とするハプロタイプ推定装置。
5. 請求項１乃至４のいずれか一つに記載のハプロタイプ推定装置において、
   上記制御部は、
   上記コピー単位を上記マーカー部位の単位で複数の区画に分け、上記区画に含まれる上記マーカー部位の上記塩基多型について、上記最大値算出手段、上記多型識別文字列挙手段、上記文字列作成手段、上記ハプロタイプ文字列格納手段、および、上記ハプロタイプ推定手段に、当該区画毎に処理を実行させるよう制御する区画化手段と、
   上記区画毎に上記ハプロタイプ推定手段により推定された上記ハプロタイプ文字列について、上記コピー単位の繰り返しパターンが同じ上記ハプロタイプ文字列同士で、上記コピー単位における上記各区画を互いに連結することにより、上記複数の上記区画から構成される上記ハプロタイプ文字列を再現する連結手段と、
   を備えたことを特徴とするハプロタイプ推定装置。
6. 請求項１乃至５のいずれか一つに記載のハプロタイプ推定装置において、
   上記ハプロタイプ推定手段は、
   ＥＭ（Expectation-Maximization）法を用いて、
   上記集団における上記ハプロタイプ文字列の頻度を、当該ハプロタイプ文字列を少なくとも一方に有する上記組合せの頻度により重み付けして算出するＭステップと、
   上記組合せの頻度を、当該組合せを構成する上記ハプロタイプ文字列の頻度の積により求め、当該組合せの頻度に基づいて上記重みを算出するＥステップと、
   を上記頻度の値が収束するまで交互に繰り返すハプロタイプ頻度算出手段、
   を更に備えたことを特徴とするハプロタイプ推定装置。
7. 請求項６に記載のハプロタイプ推定装置において、
   上記Ｍステップは、
   下記の数式１に基づいて、上記ハプロタイプ文字列の頻度を算出し、
   

   

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

   

   （ここで、Ｐ（ｄ_jk）は、上記ハプロタイプ文字列の上記組合せの頻度を表す。また、ｈ_lおよびｈ_mは当該組合せを構成する２つの上記ハプロタイプ文字列を表し、Ｐ（ｈ_l）およびＰ（ｈ_m）は、当該２つの上記ハプロタイプ文字列の頻度をそれぞれ表す。）
   を特徴とするハプロタイプ推定装置。
8. 請求項１乃至７のいずれか一つに記載のハプロタイプ推定装置において、
   上記記憶部は、
   上記各個体における上記コピー数多型のコピー数の総和を記憶するコピー数総和記憶手段、
   を更に備え、
   上記制御部は、
   上記コピー数総和記憶手段に記憶された上記各個体における上記コピー数の総和を、任意の一つの上記多型塩基の上記カウント数とした上記カウント数テーブルを作成するカウント数テーブル作成手段と、
   上記ハプロタイプ推定手段により推定された上記ハプロタイプの組み合わせにおいて、上記ハプロタイプ文字列中の上記一つの多型塩基に対応付けられた上記多型識別文字の数を算出し、算出した数を上記ハプロタイプにおける上記コピー数として設定するコピー数設定手段と、
   を更に備えたことを特徴とするハプロタイプ推定装置。
9. 集団における各個体のコピー数多型と塩基多型を含む遺伝子型データからハプロタイプを推定する、制御部と記憶部を少なくとも備えたハプロタイプ推定装置にハプロタイプ推定方法を実行させるためのプログラムであって、
   上記記憶部は、
   上記個体毎に、上記遺伝子型データを用いて、上記コピー数多型のコピー単位上の、標識によって特定されるマーカー部位に対応付けられた多型塩基をカウントすることにより得られたカウント数を、当該多型塩基の種類毎に記憶するカウント数テーブル、
   を備えており、
   上記制御部において実行される、
   上記カウント数テーブルに格納された上記個体の上記カウント数に基づいて、上記マーカー部位毎に上記カウント数の和を集計し、上記カウント数の和のうちの最大値を求める最大値算出ステップと、
   上記多型塩基の種類に対応付けられた多型識別文字を、当該多型塩基の上記カウント数列挙する多型識別文字列挙ステップと、
   上記コピー単位が上記最大値の回数繰り返される任意の文字列を作成するよう、列挙された上記多型識別文字を並び替える文字列作成ステップと、
   上記文字列作成ステップにおいて作成された上記文字列を、任意の上記コピー単位で２分し、ハプロタイプ文字列の組合せとして格納するハプロタイプ文字列格納ステップと、
   上記集団において、同一である上記ハプロタイプ文字列の数を集計し、当該ハプロタイプ文字列の上記集団における頻度を求め、当該頻度が所定の条件を満たす上記各個体の上記ハプロタイプ文字列の上記組合せを、上記ハプロタイプの組合せとして推定するハプロタイプ推定ステップと、
   を含むハプロタイプ推定方法を実行させることを特徴とするプログラム。

Images