(57)【要約】【課題】  観察された遺伝子型データと連続値を取る表現型データとに基づいて、ハプロタイプ頻度とディプロタイプ形に加え、ディプロタイプに基づいた量的表現型の分布を定める平均値及び標準偏差も同時に推定する。【解決手段】  遺伝子型データと連続値を取る表現型データとに基づいて、ハプロタイプ頻度と量的表現型の分布を定める平均値及び標準偏差とを母数とし、所定のハプロタイプを含むディプロタイプ型と所定の表現型とに関連性がないという仮説の下での尤度を最大化して得られる最大尤度（L0max）と、所定のハプロタイプを含むディプロタイプ形と所定の表現型とに関連性があるという仮説の下での尤度を最大化して得られる、ハプロタイプ頻度並びに浸透率の最尤推定値と最大尤度（Lmax）とを求めるステップaと、ステップaで求めた最尤推定値から量的表現型の分布を定める平均値及び標準偏差を求めるステップｂを含む。【選択図】  なし【特許請求の範囲】
【請求項１】
遺伝子型データと、連続値を取る表現型データとに基づいて、ハプロタイプ頻度と量的表現型の分布を定める平均値及び標準偏差とを母数とし、所定のハプロタイプを含むディプロタイプ形と連続値を取る表現型データ分布とに関連性がないという仮説の下での尤度を最大化して得られる最大尤度（L0max）と、所定のハプロタイプを含むディプロタイプ形と連続値を取る表現型データ分布とに関連性があるという仮説の下での尤度を最大化して得られる、ハプロタイプ頻度並びに量的表現型の分布を定める平均値及び標準偏差の最尤推定値と最大尤度（Lmax）とを求めるステップaと、
ステップaで求めた最尤推定値から量的表現型の分布を定める平均値及び標準偏差を求めるステップｂと
を含む量的表現型の分布を定める平均値及び標準偏差推定方法。
【請求項２】
上記ステップaでは、下記式（I）をΘ、
【数１】
2005129024.tif 000095
（なお
【数２】
2005129024.tif 000096
であり、可能な全てのディプロタイプ形に対する分布の平均値である）及びσ（標準偏差）の上で最大化させることで上記最大尤度（Lmax）を求め、
【数３】
2005129024.tif 000097
下記式（II）をΘ、
【数４】
2005129024.tif 000098
（なお
【数５】
2005129024.tif 000099
であり、可能な全てのディプロタイプ形に一定な平均値である）及びσ（標準偏差）の上で最大化させることで上記最大尤度（L0max）を求める
【数６】
2005129024.tif 000100
（上記式（I）及び（II）において、Θはハプロタイプ頻度のベクトルを示す。上記式（I）及び（II）において、
【数７】
2005129024.tif 000101
はi番目の個体がディプロタイプ形の実現値akを持つ確率である。diはi番目の個体のディプロタイプ形を表すランダム変数である。
上記式(I)において、
【数８】
2005129024.tif 000102
は、所定の表現型に関するハプロタイプの集合の要素を含むディプロタイプ形の集合をD+とし、N個体中i番目の個体のディプロタイプ形をdiとしたときに
【数９】
2005129024.tif 000103
の下で量的表現型xをきたす確率密度関数を意味する。
上記式(II)において、
【数１０】
2005129024.tif 000104
は
【数１１】
2005129024.tif 000105
の下で量的表現型xをきたす確率を意味する。また、Ψiはi番目の個体の確率変数としての表現型を示す。wiはi番目の個体の実測値としての表現型を示す。akは可能なk番目の個体のディプロタイプ形を示し、Aiはi番目の個体の遺伝子型データgiと合致するディプロタイプ形の集合を示す。）
ことを特徴とする請求項１記載の量的表現型の分布を定める平均値及び標準偏差推定方法。
【請求項３】
コホート研究又は臨床治験試験の結果から得られた、所定の集団において観察された遺伝子型データ及び表現型データを用いることを特徴とする請求項１記載の量的表現型の分布を定める平均値及び標準偏差推定方法。
【請求項４】
母数とした上記ハプロタイプ頻度に基づいて、個別のハプロタイプを区別する情報を与える座位及びこれらの座位から組み合わせによって重複する情報をもつ座位を特定し、特定した座位をマスクすることによって検定対象のハプロタイプ集合を限定して定義することを特徴とする請求項１記載の量的表現型の分布を定める平均値及び標準偏差推定方法。
【請求項５】
遺伝子型データと、連続値を取る表現型データとに基づいて、ハプロタイプ頻度と量的表現型の分布を定める平均値及び標準偏差とを母数とし、所定のハプロタイプを含むディプロタイプ形と連続値を取る表現型データ分布とに関連性がないという仮説の下での尤度を最大化して得られる最大尤度（L0max）と、所定のハプロタイプを含むディプロタイプ形と連続値を取る表現型データ分布とに関連性があるという仮説の下での尤度を最大化して得られる、ハプロタイプ頻度並びに量的表現型の分布を定める平均値及び標準偏差の最尤推定値と最大尤度（Lmax）とを求める工程aと、
工程aで求めた最尤推定値から量的表現型の分布を定める平均値及び標準偏差を求める工程ｂと
をコンピュータに実行させる量的表現型の分布を定める平均値及び標準偏差推定プログラム。
【請求項６】
上記工程aでは、下記式（I）をΘ、
【数１２】
2005129024.tif 000106
（なお
【数１３】
2005129024.tif 000107
であり、可能な全てのディプロタイプ形に対する分布の平均値である）及びσ（標準偏差）の上で最大化させることで上記最大尤度（Lmax）を求め、
【数１４】
2005129024.tif 000108
下記式（II）をΘ、
【数１５】
2005129024.tif 000109
（なお
【数１６】
2005129024.tif 000110
であり、可能な全てのディプロタイプ形に一定な平均値である）及びσ（標準偏差）の上で最大化させることで上記最大尤度（L0max）を求める
【数１７】
2005129024.tif 000111
（上記式（I）及び（II）において、Θはハプロタイプ頻度のベクトルを示す。上記式（I）及び（II）において、
【数１８】
2005129024.tif 000112
はi番目の個体がディプロタイプ形の実現値akを持つ確率である。diはi番目の個体のディプロタイプ形を表すランダム変数である。
上記式(I)において、
【数１９】
2005129024.tif 000113
は、所定の表現型に関するハプロタイプの集合の要素を含むディプロタイプ形の集合をD+とし、N個体中i番目の個体のディプロタイプ形をdiとしたときに
【数２０】
2005129024.tif 000114
の下で量的表現型xをきたす確率密度関数を意味する。
上記式(II)において、
【数２１】
2005129024.tif 000115
は
【数２２】
2005129024.tif 000116
の下で量的表現型xをきたす確率を意味する。また、Ψiはi番目の個体の確率変数としての表現型を示す。wiはi番目の個体の実測値としての表現型を示す。akは可能なk番目の個体のディプロタイプ形を示し、Aiはi番目の個体の遺伝子型データgiと合致するディプロタイプ形の集合を示す。）
ことを特徴とする請求項５記載の量的表現型の分布を定める平均値及び標準偏差推定プログラム。
【請求項７】
コホート研究又は臨床治験試験の結果から得られた、所定の集団において観察された遺伝子型データ及び表現型データを用いることを特徴とする請求項５記載の量的表現型の分布を定める平均値及び標準偏差推定プログラム。
【請求項８】
母数とした上記ハプロタイプ頻度に基づいて、個別のハプロタイプを区別する情報を与える座位及びこれらの座位から組み合わせによって重複する情報をもつ座位を特定し、特定した座位をマスクすることによって検定対象のハプロタイプ集合を限定して定義することを特徴とする請求項５記載の量的表現型の分布を定める平均値及び標準偏差推定プログラム。
【請求項９】
遺伝子型データと、連続値を取る表現型データとに基づいて、ハプロタイプ頻度と量的表現型の分布を定める平均値及び標準偏差とを母数とし、所定のハプロタイプを含むディプロタイプ形と連続値を取る表現型データ分布とに関連性がないという仮説の下での尤度を最大化して得られる最大尤度（L0max）と、所定のハプロタイプを含むディプロタイプ形と連続値を取る表現型データ分布とに関連性があるという仮説の下での尤度を最大化して得られる、ハプロタイプ頻度並びに量的表現型の分布を定める平均値及び標準偏差の最尤推定値と最大尤度（Lmax）とを求めるステップaと、
ステップaで求めた最大尤度（L0max）及び最大尤度（Lmax）から尤度比を求め、所定のハプロタイプを含むディプロタイプ形と所定の量的表現型とに関連性があるという仮説をχ2分布により検定するステップbと
を含むハプロタイプと量的表現型との関連性検定方法。
【請求項１０】
上記ステップaでは、下記式（I）をΘ、
【数２３】
2005129024.tif 000117
（なお
【数２４】
2005129024.tif 000118
であり、可能な全てのディプロタイプ形に対する分布の平均値である）
及びσ（標準偏差）の上で最大化させることで上記最大尤度（Lmax）を求め、
【数２５】
2005129024.tif 000119
下記式（II）をΘ、
【数２６】
2005129024.tif 000120
（なお
【数２７】
2005129024.tif 000121
であり、可能な全てのディプロタイプ形に一定な平均値である）及びσ（標準偏差）の上で最大化させることで上記最大尤度（L0max）を求める
【数２８】
2005129024.tif 000122
（上記式（I）及び（II）において、Θはハプロタイプ頻度のベクトルを示す。上記式（I）及び（II）において、
【数２９】
2005129024.tif 000123
はi番目の個体がディプロタイプ形の実現値akを持つ確率である。diはi番目の個体のディプロタイプ形を表すランダム変数である。
上記式(I)において、
【数３０】
2005129024.tif 000124
は、所定の表現型に関するハプロタイプの集合の要素を含むディプロタイプ形の集合をD+とし、N個体中i番目の個体のディプロタイプ形をdiとしたときに
【数３１】
2005129024.tif 000125
の下で量的表現型xをきたす確率密度関数を意味する。
上記式(II)において、
【数３２】
2005129024.tif 000126
は
【数３３】
2005129024.tif 000127
の下で量的表現型xをきたす確率を意味する。また、Ψiはi番目の個体の確率変数としての表現型を示す。wiはi番目の個体の実測値としての表現型を示す。akは可能なk番目の個体のディプロタイプ形を示し、Aiはi番目の個体の遺伝子型データgiと合致するディプロタイプ形の集合を示す。）
ことを特徴とする請求項９記載の関連性検定方法。
【請求項１１】
上記ステップｂでは、統計量として-2log（Lmax/L0max）を求め（ただしlogは自然対数を示す。）、所定のハプロタイプを含むディプロタイプ形と量的表現型とが無関係の場合には当該統計量が自由度１のχ2分布に漸近的に従うので、当該統計量が限界値χ2（ここで、限界値χ2は自由度１のχ2分布において、累積分布関数が1-αとなる値である。ここで、αは検定の危険率である。）を超えない場合は所定のハプロタイプを含むディプロタイプ形と所定の量的表現型とに関連性があるとは言えないと判断し、当該統計量が限界値χ2を超える場合は所定のハプロタイプと所定の量的表現型とに関連性があると判断することを特徴とする請求項９記載の関連性検定方法。
【請求項１２】
コホート研究又は臨床治験試験の結果から得られた、所定の集団において観察された遺伝子型データ及び連続値を取る表現型データを用いることを特徴とする請求項９記載の関連性検定方法。
【請求項１３】
母数とした上記ハプロタイプ頻度に基づいて、個別のハプロタイプを区別する情報を与える座位及びこれらの座位から組み合わせによって重複する情報をもつ座位を特定し、特定した座位をマスクすることによって検定対象のハプロタイプ集合を限定して定義することを特徴とする請求項９記載の関連性検定方法。
【請求項１４】
遺伝子型データと、連続値を取る表現型データとに基づいて、ハプロタイプ頻度と量的表現型の分布を定める平均値及び標準偏差とを母数とし、所定のハプロタイプを含むディプロタイプ形と連続値を取る表現型データ分布とに関連性がないという仮説の下での尤度を最大化して得られる最大尤度（L0max）と、所定のハプロタイプを含むディプロタイプ形と連続値を取る表現型データ分布とに関連性があるという仮説の下での尤度を最大化して得られる、ハプロタイプ頻度並びに量的表現型の分布を定める平均値及び標準偏差の最尤推定値と最大尤度（Lmax）とを求める工程aと、
ステップaで求めた最大尤度（L0max）及び最大尤度（Lmax）から尤度比を求め、所定のハプロタイプと所定の量的表現型とに関連性があるという仮説をχ2分布により検定する工程bと
をコンピュータに実行させるハプロタイプと量的表現型との関連性検定プログラム。
【請求項１５】
上記工程aでは、下記式（I）をΘ、
【数３４】
2005129024.tif 000128
（なお
【数３５】
2005129024.tif 000129
であり、可能な全てのディプロタイプ形に対する分布の平均値である）及びσ（標準偏差）の上で最大化させることで上記最大尤度（Lmax）を求め、
【数３６】
2005129024.tif 000130
下記式（II）をΘ、
【数３７】
2005129024.tif 000131
（なお
【数３８】
2005129024.tif 000132
であり、可能な全てのディプロタイプ形に一定な平均値である）及びσ（標準偏差）の上で最大化させることで上記最大尤度（L0max）を求める
【数３９】
2005129024.tif 000133
（上記式（I）及び（II）において、Θはハプロタイプ頻度のベクトルを示す。上記式（I）及び（II）において、
【数４０】
2005129024.tif 000134
はi番目の個体がディプロタイプ形の実現値akを持つ確率である。diはi番目の個体のディプロタイプ形を表すランダム変数である。
上記式(I)において、
【数４１】
2005129024.tif 000135
は、所定の表現型に関するハプロタイプの集合の要素を含むディプロタイプ形の集合をD+とし、N個体中i番目の個体のディプロタイプ形をdiとしたときに
【数４２】
2005129024.tif 000136
の下で量的表現型xをきたす確率密度関数を意味する。
上記式(II)において、
【数４３】
2005129024.tif 000137
は
【数４４】
2005129024.tif 000138
の下で量的表現型xをきたす確率を意味する。また、Ψiはi番目の個体の確率変数としての表現型を示す。wiはi番目の個体の実測値としての表現型を示す。akは可能なk番目の個体のディプロタイプ形を示し、Aiはi番目の個体の遺伝子型データgiと合致するディプロタイプ形の集合を示す。）
ことを特徴とする請求項１４記載の関連性検定プログラム。
【請求項１６】
上記工程ｂでは、統計量として-2log（Lmax/L0max）を求め（ただしlogは自然対数を示す。）、所定のハプロタイプを含むディプロタイプ形と量的表現型とが無関係の場合には当該統計量が自由度１のχ2分布に漸近的に従うので、当該統計量が限界値χ2（ここで、限界値χ2は自由度１のχ2分布において、累積分布関数が1-αとなる値である。ここで、αは検定の危険率である。）を超えない場合は所定のハプロタイプを含むディプロタイプ形と所定の量的表現型とに関連性があるとは言えないと判断し、当該統計量が限界値χ2を超える場合は所定のハプロタイプと所定の量的表現型とに関連性があると判断することを特徴とする請求項１４記載の関連性検定プログラム。
【請求項１７】
コホート研究又は臨床治験試験の結果から得られた、所定の集団において観察された遺伝子型データ及び連続値を取る表現型データを用いることを特徴とする請求項１４記載の関連性検定プログラム。
【請求項１８】
母数とした上記ハプロタイプ頻度に基づいて、個別のハプロタイプを区別する情報を与える座位及びこれらの座位から組み合わせによって重複する情報をもつ座位を特定し、特定した座位をマスクすることによって検定対象のハプロタイプ集合を限定して定義することを特徴とする請求項１４記載の関連性検定プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、コホート研究、臨床治験試験の結果から得られた、所定の集団において観察された遺伝子型データ及び連続値を取る表現型データを用いた量的表現型の分布を定める平均値及び標準偏差推定方法に関し、また、量的表現型の分布を定める平均値及び標準偏差推定方法を用いることによって得られた推定値を用いたハプロタイプと量的表現型との関連性検定方法に関する。
【背景技術】
【０００２】
連鎖不平衡（LD）とハプロタイプ構造に基づいた多型の解析は益々重要と成りつつある。ここで、一つの配偶子に存在する、連鎖する複数の対立遺伝子の組合せをハプロタイプと定義する。また、異なる対立遺伝子の非独立を連鎖不平衡と定義する。なお、対立遺伝子は、分子生物学的には片方の染色体における多型を意味する。
【０００３】
多数の個体における多数の連鎖する座位（対立遺伝子）に関する、最近のデータ解析によりヒトゲノムの集団遺伝学的視点からの構造が明らかになってきた。即ち、ヒトゲノムは、ハプロタイプブロック構造（またはLDブロック）を持っていることが明らかになった。一つのブロックの中ではLDは強く、主たるハプロタイプの数は限られている。このブロックの中のSNP (single nucleotide polymorphism, 一塩基置換)からタグSNPを抽出し、関連解析に用いることができる。即ち、ハプロタイプ構造は形質の詳細なマッピングに非常に有益である。
【０００４】
また、ハプロタイプは、ヒトの遺伝的情報をさまざまな表現型と関連付けるためにも有効である。一つの表現型はそれぞれのSNPと関連しているだけではなく、ハプロタイプと関連していることも多い。我々はハプロタイプを完全データ、遺伝子型（例えばSNP遺伝子型）を不完全データと解釈することができる。その理由はハプロタイプデータから遺伝子型データは完全に復活できるが逆は不可能であるからである。
【０００５】
実際、我々は一つの座位の一つのアレルを、その座位にそのアレルを持った複数のハプロタイプの集合と再定義することができる。従って、ハプロタイプやディプロタイプ形（ハプロタイプの組合せ）に基づいて、多型と表現型との関連を考察するほうが、アレルや遺伝子型に基づいて考察するよりもより一般的ものといえる。なお、最近の多くの研究によれば、例えば糖尿病や薬剤反応性などの表現型は、SNP等の遺伝子型ではなくハプロタイプやディプロタイプ形に関連していたと報告されている。
【０００６】
しかし、実験的に直接ハプロタイプを決めるという報告がなされているが、個人のハプロタイプやディプロタイプ形は容易に観察することは不可能である。ハプロタイプを直接観察する代わりに、多座位遺伝子型データから色々なアルゴリズムにより推定される。即ち、Clarkのアルゴリズム[非特許文献１参照]、EM(expectation-maximization)アルゴリズム[非特許文献２〜６参照]、PHASEアルゴリズム、PLアルゴリズム、およびPL-EMアルゴリズムなどが提案されている。
【０００７】
ハプロタイプに基づいたゲノムレベルの関連解析では、上記アルゴリズムのいずれかによりハプロタイプ頻度が推定される。そして推定された値に基づいて、ケースとコントロールのグループ間でハプロタイプ頻度を比較する。Zaykinらは、離散、あるいは連続の形質について一つに集中していないディプロタイプ形の存在下でも可能な、回帰に基づくハプロタイプ頻度と形質の関連解析のアルゴリズムを発表した[非特許文献７参照]。また、Fallinらも、ケース/コントロールのサンプルデザインに基づく、推定されたハプロタイプ頻度と形質の関連を検定する手法を提案した[非特許文献８参照]。
【０００８】
しかし、個体レベルの表現型に焦点を置く場合、それらはハプロタイプよりディプロタイプ形に基礎を置いているはずである。この関係は、一つの座位についての遺伝子型とアレルの関係と同様である。従って、表現型と遺伝情報の関連を調べるためには、場合によっては、ディプロタイプ形の異なる個体の間での罹患者の割合を比較するほうが、異なった表現型の間でのハプロタイプの頻度を比較するよりも有効である。
【０００９】
ディプロタイプ形の異なった個体間で罹患者の割合を比較するためには、それぞれの個体のディプロタイプ形が上記のアルゴリズムの一つを用いて推定され、すべての個体が一つのハプロタイプか一つのディプロタイプ形のあるなしに基づいて分類される。このように分類されたそれぞれの集団を罹患者、非罹患者でさらに分類した後、独立性の検定が行われる。
【００１０】
しかし問題は、少なくとも幾人かの個体についてはディプロタイプ形が完全に決まらないことである（あいまいなディプロタイプ形）。真のハプロタイプ情報ではなく、推定されたハプロタイプ情報により個人を分類することにより起きるタイプI過誤の程度は明らかではない。
【００１１】
【非特許文献１】Clark AG, Weiss KM, Nickerson DA et al (1998) Haplotype structure and population genetic inferences from nucleotide-sequence variation in human lipoprotein lipase、 Am. J. Hum. Genet. 63: 595-612
【非特許文献２】Excoffier L, Slatkin M (1995) Maximum-likelihood estimation of molecular haplotype frequencies in a diploid population、Mol. Biol. Evol. 12: 921-927
【非特許文献３】Hawley ME, Kidd KK (1995) HAPLO: a program using the EM algorithm to estimate the frequencies of multi-site haplotypes、J. Hered. 86: 409-411
【非特許文献４】Schneider S, Roessli D, Excoffier L (2000) Arlequin: a software for population genetics data analysis、Ver 2.000、Genetics and Biometry Laboratory、Department of Anthropology、University of Geneva、Geneva
【非特許文献５】Long JC, Williams RC, Urbanek M (1995) An E-M algorithm and testing strategy for multiple locus haplotypes、Am. J. Hum. Genet. 56: 799-810
【非特許文献６】Kitamura Y, Moriguchi M, Kaneko H, Morisaki H, Morisaki T, Toyama K, Kamatani N (2002)、Determination of probability distribution of diplotype configuration (diplotype distribution) for each subject from genotypic data using the EM algorithm、Ann. Hum. Genet. 66: 183-193
【非特許文献７】Zaykin DV, Westfall PH, Young SS, Karnoub MA, Wagner MJ, Ehm MG (2002)、Testing association of statistically inferred haplotypes with discrete and continuous traits in samples of unrelated individuals、Hum. Hered. 53: 79-91
【非特許文献８】Fallin D, Schork NJ (2000)、Accuracy of haplotype frequency estimation for biallelic loci, via the expectation-maximization algorithm for unphased diploid genotype data、Am. J. Hum. Genet. 67: 947-959
【発明の開示】
【発明が解決しようとする課題】
【００１２】
そこで、本発明は、上述した問題に鑑み、観察された遺伝子型データと、連続値を取る表現型データとに基づいて、ハプロタイプ頻度とディプロタイプ形に加え、ディプロタイプに基づいた量的表現型の分布を定める平均値及び標準偏差推定方法に関し、また、量的表現型の分布を定める平均値及び標準偏差推定方法を用いることによって得られた推定値を用いたハプロタイプと量的表現型との関連性検定方法を提供することを目的とする。
【課題を解決するための手段】
【００１３】
上述した目的を達成した本発明に係るアルゴリズムは、それぞれの個体のディプロタイプ形を疑問の余地無く決定する必要が無く、集団のハプロタイプ頻度、それぞれの個体のディプロタイプ分布（ディプロタイプ形の事後確率分布）、量的表現型の分布を定める平均値及び標準偏差を遺伝子型データと連続値を取る表現型データが与えられた下で最尤推定することができる。本発明に係るアルゴリズムを用いると、例えば、コホート研究や臨床治療試験から得られた遺伝子型データと表現型のデータとを用いてハプロタイプの存在と一つの表現型の関連を検定することができる。本発明者らは、このアルゴリズムの有効性をシミュレーションと現実データの両方を用いて検討し、コホート研究や臨床治療試験の遺伝子型データと連続値を取る表現型データの解析に非常に有効なことを見出し、本発明を完成するに至った。
すなわち、本発明は以下を包含する。
【００１４】
（１）  遺伝子型データと、連続値を取る表現型データとに基づいて、ハプロタイプ頻度と量的表現型の分布を定める平均値及び標準偏差とを母数とし、所定のハプロタイプを含むディプロタイプ形と連続値を取る表現型データ分布とに関連性がないという仮説の下での尤度を最大化して得られる最大尤度（L0max）と、所定のハプロタイプを含むディプロタイプ形と連続値を取る表現型データ分布とに関連性があるという仮説の下での尤度を最大化して得られる、ハプロタイプ頻度並びに量的表現型の分布を定める平均値及び標準偏差の最尤推定値と最大尤度（Lmax）とを求めるステップaと、
ステップaで求めた最尤推定値から量的表現型の分布を定める平均値及び標準偏差を求めるステップｂと
を含む量的表現型の分布を定める平均値及び標準偏差推定方法。
【００１５】
（２）  上記ステップaでは、下記式（I）をΘ、
【数１】
2005129024.tif 000002
（なお
【００１６】
【数２】
2005129024.tif 000003
であり、可能な全てのディプロタイプ形に対する分布の平均値である）及びσ（標準偏差）の上で最大化させることで上記最大尤度（Lmax）を求め、
【００１７】
【数３】
2005129024.tif 000004
下記式（II）をΘ、
【００１８】
【数４】
2005129024.tif 000005
（なお
【００１９】
【数５】
2005129024.tif 000006
であり、可能な全てのディプロタイプ形に一定な平均値である）及びσ（標準偏差）の上で最大化させることで上記最大尤度（L0max）を求める
【００２０】
【数６】
2005129024.tif 000007
（上記式（I）及び（II）において、Θはハプロタイプ頻度のベクトルを示す。上記式（I）及び（II）において、
【００２１】
【数７】
2005129024.tif 000008
はi番目の個体がディプロタイプ形の実現値akを持つ確率である。diはi番目の個体のディプロタイプ形を表すランダム変数である。
【００２２】
上記式(I)において、
【数８】
2005129024.tif 000009
は、所定の表現型に関するハプロタイプの集合の要素を含むディプロタイプ形の集合をD+とし、N個体中i番目の個体のディプロタイプ形をdiとしたときに
【００２３】
【数９】
2005129024.tif 000010
の下で量的表現型xをきたす確率密度関数を意味する。
【００２４】
上記式(II)において、
【数１０】
2005129024.tif 000011
は
【００２５】
【数１１】
2005129024.tif 000012
の下で量的表現型xをきたす確率を意味する。また、Ψiはi番目の個体の確率変数としての表現型を示す。wiはi番目の個体の実測値としての表現型を示す。akは可能なk番目の個体のディプロタイプ形を示し、Aiはi番目の個体の遺伝子型データgiと合致するディプロタイプ形の集合を示す。）
ことを特徴とする（１）記載の量的表現型の分布を定める平均値及び標準偏差推定方法。
【００２６】
（３）  コホート研究又は臨床治験試験の結果から得られた、所定の集団において観察された遺伝子型データ及び表現型データを用いることを特徴とする（１）記載の量的表現型の分布を定める平均値及び標準偏差推定方法。
【００２７】
（４）  母数とした上記ハプロタイプ頻度に基づいて、個別のハプロタイプを区別する情報を与える座位及びこれらの座位から組み合わせによって重複する情報をもつ座位を特定し、特定した座位をマスクすることによって検定対象のハプロタイプ集合を限定して定義することを特徴とする（１）記載の量的表現型の分布を定める平均値及び標準偏差推定方法。
【００２８】
（５）  遺伝子型データと、連続値を取る表現型データとに基づいて、ハプロタイプ頻度と量的表現型の分布を定める平均値及び標準偏差とを母数とし、所定のハプロタイプを含むディプロタイプ形と連続値を取る表現型データ分布とに関連性がないという仮説の下での尤度を最大化して得られる最大尤度（L0max）と、所定のハプロタイプを含むディプロタイプ形と連続値を取る表現型データ分布とに関連性があるという仮説の下での尤度を最大化して得られる、ハプロタイプ頻度並びに量的表現型の分布を定める平均値及び標準偏差の最尤推定値と最大尤度（Lmax）とを求めるステップaと、
ステップaで求めた最大尤度（L0max）及び最大尤度（Lmax）から尤度比を求め、所定のハプロタイプを含むディプロタイプ形と所定の量的表現型とに関連性があるという仮説をχ2分布により検定するステップbと
を含むハプロタイプと量的表現型との関連性検定方法。
【００２９】
（６）  上記ステップaでは、下記式（I）をΘ、
【数１２】
2005129024.tif 000013
（なお
【００３０】
【数１３】
2005129024.tif 000014
であり、可能な全てのディプロタイプ形に対する分布の平均値である）及びσ（標準偏差）の上で最大化させることで上記最大尤度（Lmax）を求め、
【００３１】
【数１４】
2005129024.tif 000015
下記式（II）をΘ、
【００３２】
【数１５】
2005129024.tif 000016
（なお
【００３３】
【数１６】
2005129024.tif 000017
であり、可能な全てのディプロタイプ形に一定な平均値である）及びσ（標準偏差）の上で最大化させることで上記最大尤度（L0max）を求める
【００３４】
【数１７】
2005129024.tif 000018
（上記式（I）及び（II）において、Θはハプロタイプ頻度のベクトルを示す。上記式（I）及び（II）において、
【００３５】
【数１８】
2005129024.tif 000019
はi番目の個体がディプロタイプ形の実現値akを持つ確率である。diはi番目の個体のディプロタイプ形を表すランダム変数である。
【００３６】
上記式(I)において、
【数１９】
2005129024.tif 000020
は、所定の表現型に関するハプロタイプの集合の要素を含むディプロタイプ形の集合をD+とし、N個体中i番目の個体のディプロタイプ形をdiとしたときに
【００３７】
【数２０】
2005129024.tif 000021
の下で量的表現型xをきたす確率密度関数を意味する。
【００３８】
上記式(II)において、
【数２１】
2005129024.tif 000022
は
【００３９】
【数２２】
2005129024.tif 000023
の下で量的表現型xをきたす確率を意味する。また、Ψiはi番目の個体の確率変数としての表現型を示す。wiはi番目の個体の実測値としての表現型を示す。akは可能なk番目の個体のディプロタイプ形を示し、Aiはi番目の個体の遺伝子型データgiと合致するディプロタイプ形の集合を示す。）
ことを特徴とする（５）記載の関連性検定方法。
【００４０】
（６）  上記ステップｂでは、統計量として-2log（Lmax/L0max）を求め（ただしlogは自然対数を示す。）、所定のハプロタイプを含むディプロタイプ形と量的表現型とが無関係の場合には当該統計量が自由度１のχ2分布に漸近的に従うので、当該統計量が限界値χ2（ここで、限界値χ2は自由度１のχ2分布において、累積分布関数が1-αとなる値である。ここで、αは検定の危険率である。）を超えない場合は所定のハプロタイプを含むディプロタイプ形と所定の量的表現型とに関連性があるとは言えないと判断し、当該統計量が限界値χ2を超える場合は所定のハプロタイプと所定の量的表現型とに関連性があると判断することを特徴とする（５）記載の関連性検定方法。
【００４１】
（８）  コホート研究又は臨床治験試験の結果から得られた、所定の集団において観察された遺伝子型データ及び連続値を取る表現型データを用いることを特徴とする（５）記載の関連性検定方法。
【００４２】
（９）  母数とした上記ハプロタイプ頻度に基づいて、個別のハプロタイプを区別する情報を与える座位及びこれらの座位から組み合わせによって重複する情報をもつ座位を特定し、特定した座位をマスクすることによって検定対象のハプロタイプ集合を限定して定義することを特徴とする（５）記載の関連性検定方法。
【００４３】
また、本発明は、上記（１）〜（９）いずれかの方法における各ステップを、コンピュータに実行させるためのプログラムである。
【発明の効果】
【００４４】
本発明によれば、遺伝子型データ及び連続値を取る表現型データを用いて、個体のディプロタイプ形を疑問の余地無く決定する必要が無く、集団のハプロタイプ頻度、それぞれの個体のディプロタイプ形の事後確率分布、量的表現型の分布を定める平均値及び標準偏差を最尤推定することができる。本発明に係るアルゴリズムを用いると、例えば、コホート研究や臨床治療試験から得られた遺伝子型データと連続値を取る表現型のデータとを用いてハプロタイプの存在と量的表現型が示す値の範囲との関連を検定することが可能となる。
【発明を実施するための最良の形態】
【００４５】
以下、本発明を詳細に説明する。
本発明に係る量的表現型の分布を定める平均値及び標準偏差推定方法及びディプロタイプ形と量的表現型との関連性検定方法は、以下に説明するアルゴリズムによって実現される。本発明で用いられるアルゴリズム（以下、「本アルゴリズム」と称する）においては、コホート研究又は臨床治療試験から得られた個人の遺伝子型データと連続値を取る表現型データとを用いて、量的表現型と所定のハプロタイプの存在の関連を検定し、ハプロタイプを基礎とした量的表現型の分布を定める平均値及び標準偏差を推定する。本アルゴリズムは、EM(expectation-maximization)アルゴリズムに基づいて作成されたものである。
【００４６】
本アルゴリズムは、集団のハプロタイプ頻度に加え、所定のハプロタイプを保有している個体と保有していない個体の間で異なる量的表現型の分布を定める平均値及び標準偏差を推定することができる。従って、本アルゴリズムによれば、相対危険も最尤推定することができる。
【００４７】
本発明において、「連続値を取る表現型データ」とは、例えば、血圧値等の臨床検査データ、薬剤の血中濃度、投与量、DNAマイクロアレイの遺伝子発現量、タンパク質発現量といった連続値を取る表現型を示すデータを意味する。当該表現型は、いわゆる量的形質座位（QTL; quantitative trait locus）における表現型を意味する。以下、連続値を取る表現型データを「QTL表現型」と称する。「量的表現型の分布を定める平均値及び標準偏差」とは、連続値を取る表現型データの分布に関するパラメータであって、「QTL表現型を取る値の範囲」と言い換えることができる。
【００４８】
具体的に、本アルゴリズムにおいては、先ず、連続値を取る表現型と関連するハプロタイプは無いと言う仮説の下での最大尤度(L0max)(帰無仮説)および、連続値を取る表現型と関連するハプロタイプはあると言う仮説の下での最大尤度(Lmax)（対立仮説）を計算する。次に、本アルゴリズムでは、統計量、例えば-2log(L0max/Lmax)（以下単に統計量と称する）を算出し、この統計量に基づいて、連続値を取る表現型とハプロタイプとの関連の検定を行う。
【００４９】
本アルゴリズムは、検体の遺伝子型情報から、量的表現型の分布を定める平均値及び標準偏差を推定する方法に適用することができる。すなわち、本アルゴリズムによれば、検体の遺伝子型情報に基づいて、当該検体において所定の量的表現型の分布を定める平均値及び標準偏差を推定することができる。これにより、本アルゴリズムは特に、遺伝的要因と、血圧値等の臨床検査データ、薬剤の血中濃度、投与量、DNAマイクロアレイの遺伝子発現量、タンパク質発現量といった連続値を取る個人におけるQTL表現型との関連の解析に有用である。
【００５０】
本アルゴリズムは、コンピュータソフトウェアQTLhaploに搭載することで、コンピュータ上で実行することができる。ここでコンピュータとは、動作を全て制御するCPUと、プログラムの実行指示等を入力できるキーボード及びマウス等の入力装置と、ディスプレイ装置等の表示装置と、一時的な情報及びプログラム等が記録されるメモリーと、各種データやプログラム等が格納されたハードディスク等の記憶装置とを備えるものである。また、コンピュータは、インターネット等の通信回線網を介して外部のデータベースや他のコンピュータ等と接続されたものであっても良い。
【００５１】
本アルゴリズムをコンピュータ上で実行する場合、コンピュータソフトウェアQTLhaploをインストールする。これにより、コンピュータは、コンピュータソフトウェアQTLhaploに従って、CPUの制御のもと本アルゴリズムを実行することができる。なお、遺伝子型データ及び連続値を取る表現型データは、例えばインターネット等の通信回線網を介して取得することもできる。
【００５２】
すなわち、本アルゴリズムは、コンピュータを、遺伝子型データ及び連続値を取る表現型データを取得する入力手段と、入力手段で取得した遺伝子型データ及び連続値を取る表現型データを用いてハプロタイプ頻度並びに並びに量的表現型の分布を定める平均値及び標準偏差の最尤推定値と最大尤度（Lmax）とを求める制御手段（演算手段）として機能させるものである。
【００５３】
また、本アルゴリズムは、当該制御手段によって、上記最大尤度（L0max）及び最大尤度（Lmax）から尤度比を求め、所定のディプロタイプ形と所定の量的表現型とに関連性があるという仮説をχ2分布により検定するようにコンピュータを機能させることができる。
【００５４】
ここで、遺伝子型データとは、ある個体に関して、いわゆるSNPタイピング等を実施した結果として得られる多型の位置を意味する情報と、多型の種類を意味する情報とを含むデータである。遺伝子型データは、個体を特定する情報を除くことにより匿名化されていてもよい。
【００５５】
連続値を取る表現型データとは、ある個体に関して、量的表現型の連続値における所定の数値を意味するデータ、或いは、量的表現型の連続値における所定の範囲に入るか否かを意味するデータである。
【００５６】
本アルゴリズム
以下、本アルゴリズムを詳細に説明する。
＜標本空間＞
l個の連鎖するSNP座位があるとする。すべての可能なハプロタイプの数はL=2lである。さらに、無限のハプロタイプコピーの集合を定義する。ここで、ハプロタイプの頻度はΘ＝(θ1,..,θJ,..,θL)であり、ここでθJはjのハプロタイプ頻度である。ただし、
【００５７】
【数２３】
2005129024.tif 000024
である。N人の個体のそれぞれに、ハプロタイプコピーの集合よりランダムに抽出して、二つのハプロタイプコピーを順番に与える。ディプロタイプ形は、順序付けられた二つのハプロタイプコピーと定義する。a1、a2、…、aＬ２を可能なディプロタイプ形とする。i番目の個体のディプロタイプ形がaｋである確率はＰ（di＝ak｜Θ）＝θlθmである。ここにdiはi番目の個体のディプロタイプ形を表すランダム変数であり、l及びmはakを構成するハプロタイプの順番である。これはハプロタイプレベルでのハーディーワインバーグ平衡が仮定されていることを意味している。
【００５８】
i番目の個体は、確率密度関数に従ってQTL表現型ψｉを発現する。QTL表現型は、diに依存した平均値と、diに依存しない固定の標準偏差に従う正規分布に従うものと仮定する。実験の結果の集合を（Θ,Ｄ,Ψ）と表し、Ｄ＝（d1,...,dN）は個体iのディプロタイプ形のベクトル、Ψ＝（ψ1,...,ψN）は個体iのQTL表現型のベクトルを表す。平均値μは、他と異なる影響を持つ表現型に関係したハプロタイプhbをdiが含むか否かに依存するものとする。Ｄ＋はハプロタイプhbを含むディプロタイプ形の集合を表すものとする。これにより、ディプロタイプ形に依存したQTL表現型に関して、ただ2つの正規分布を定義することとなる。一つは
【００５９】
【数２４】
2005129024.tif 000025
の場合の分布を与えるＮ（μ1,σ）であり、もう一つは、
【００６０】
【数２５】
2005129024.tif 000026
の場合の分布を与えるＮ（μ2,σ）である。
【００６１】
ここで、
【数２６】
2005129024.tif 000027
の場合に個体iがQTL表現型xを発現する確率密度関数をfμ1(x)とし、
【００６２】
【数２７】
2005129024.tif 000028
の場合に個体iがQTL表現型xを発現する確率密度関数をfμ2(x)とする。
【００６３】
ここで、ψiを個体iのQTL表現型とすると、
【数２８】
2005129024.tif 000029
及び
【００６４】
【数２９】
2005129024.tif 000030
である。ここで、Θと、fμ1(x)及びfμ2(x)は独立であり、ψiはdiの条件の下でΘと独立であることに注意する。
【００６５】
なお、本アルゴリズムにおいて、理論的にはすべてのディプロタイプ形に対して平均値及び標準偏差によって規定される量的表現型の分布関数を仮定することが可能である。しかし、すべてのディプロタイプ形に当該分布関数を対応付けることは現実的ではない。そこで、本アルゴリズムにおいては、上述したように、ただ二つの分布関数fμ1(x)及びfμ2(x)を仮定する。ここで、検定対象である表現型に関連したハプロタイプhbは、唯一つのハプロタイプとするばかりではなく、ハプロタイプの部分集合H+として定義することができる。即ちHallをすべてのハプロタイプの集合とし、H+をHallの部分集合で、H+の要素がディプロタイプ形に含まれていると他と異なった表現型をきたすハプロタイプの集合とする。典型的な例ではH+はただ一つのハプロタイプを含むが、複数のハプロタイプを要素として含むこともできる。H+は、本アルゴリズムにおける検定対象ハプロタイプの集合となる。もし、H+が特定の座位で特定のアレルを含むすべてのハプロタイプの集合と定義すれば、（ハプロタイプではなく）アレルと量的表現型との関連を検定することと同じになる。
【００６６】
また、Hallの部分集合としてはH+に限定されず、以下に説明するHlを定義してもよい。Hlは、EMアルゴリズムによって推定されたハプロタイプ分布に基づいて、個別のハプロタイプを区別する情報を与える座位及びこれらの座位から組み合わせによって重複する情報をもつ座位を特定し、特定した座位をマスクすることによってHallの部分集合として定義される。
【００６７】
ここでマスクとは、ハプロタイプを構成する複数の座位のうち、１以上の特定の座位については全ての多型が当てはまるものとして情報を隠蔽することを意味する。したがって、マスクによって一部が隠蔽された不完全ハプロタイプは、複数のハプロタイプを要素として含む集合Hlとして表現される概念である。その定義から明らかに
【００６８】
【数３０】
2005129024.tif 000031
となる。さらに、唯一つの座位の情報を用いて他の座位をマスクした不完全ハプロタイプを構築すると、使用した座位のSNP情報のみが用いられることからSNPと同義となる。これを集合HSNPとして定義すると、HSNPはHlの特別な場合であり、
【００６９】
【数３１】
2005129024.tif 000032
となる。
【００７０】
検定対象としてH+の代わりに不完全ハプロタイプHlを用いる場合の合理性は以下のとおりである。
【００７１】
１）ハプロタイプで遺伝子多型を表現する場合、多型の原因は塩基置換と組み換えによるものである。ある領域のハプロタイプがある量的表現型と関連する原因座位と関連するときに、複数のハプロタイプが対応付けられるならば、それは原因座位が発生した後に生じた突然変異、あるいは組み換えによるものである。不完全ハプロタイプを構築することにより、突然変異は特定の座位のマスク、組み換えは連続する座位のマスクとして表現することが可能である。
【００７２】
２）L座位のSNP情報を用いて不完全ハプロタイプを構築する際に、マスクする座位を0からL-1まで変化させることにより、全ての情報を用いたハプロタイプからSNPまでを本アルゴリズムの検定対象に含めることが可能となる。
【００７３】
ここで、L座位のSNP情報を用いて全ての不完全ハプロタイプHlを構築することを考えると、座位ごとに2つのアレルとマスク操作からなる3通りの情報を適用することとなるため、単純には3L-1とおりの組み合わせを考慮する必要が生じる。ハプロタイプ推定における組み合わせ数が2Lとおりであることと比較しても、膨大な組み合わせ数となる。しかしながら、連鎖不平衡が強い領域においては、現実にはハプロタイプは10にも満たない場合がほとんどであり、単純に全ての組み合わせを構築する必要はない。したがって、不完全ハプロタイプHlは以下の１〜３のアルゴリズムに従って構築することができる。
1.EMアルゴリズムによるハプロタイプ推定を行う。
2.推定されたハプロタイプ分布より、個別のハプロタイプを区別する情報を与える座位を抽出し、さらにそれらの座位から組み合わせによって重複する情報を持つものを削除することによって、ハプロタイプタグSNP（以下、htSNPとする）を定める。
3.htSNPの座位に対して、2つのアレルとマスク操作からなる3通りの情報を順次適用して、不完全ハプロタイプHlを構築する。さらに、異なるマスク方法でも同じHlが構築されるケースがあるため、これを削除する。
【００７４】
以下においては、H+を検定対象とする場合を説明するが、この場合と同様にして検定対象として不完全ハプロタイプHlとすることができる。
＜尤度関数＞
本アルゴリズムに用いる観察データは、個体の遺伝子型データ及びQTL表現型データである。ここで、遺伝子型データのベクトルをGobs=(g1,g2,..,gN)とし、表現型データのベクトルをΨobs=(w1,w2,..,wN)とする。ここでgiとwiは、それぞれi番目の個体の観察される遺伝子型、QTL表現型である。可能なすべてのディプロタイプ形に対する分布の平均を以下のように定義する。
【００７５】
【数３２】
2005129024.tif 000033
そうすると、尤度関数は次のようになる。
【００７６】
【数３３】
2005129024.tif 000034
ここでAiはi番目の個体についてgiに合致するakの集合である。diは
【００７７】
【数３４】
2005129024.tif 000035
及びσと独立であり、ψiはdiの条件の下でΘと独立であることから、
【００７８】
【数３５】
2005129024.tif 000036
が得られる。ここで、再び、Aiはgiに合致する個体iのディプロタイプ形の集合である。個体iは次式で与えられる確率密度関数によってQTL表現型xを発現する。
【００７９】
【数３６】
2005129024.tif 000037
帰無仮説の場合、すなわちQTL表現型の分布がディプロタイプ形と無関係であるときの尤度関数は、
【００８０】
【数３７】
2005129024.tif 000038
である。ここで、帰無仮説の場合はディプロタイプ形によらずQTL表現型の分布の平均値は一定であり、これを
【００８１】
【数３８】
2005129024.tif 000039
とする。また、Aiはgiに合致する個体iのディプロタイプ形の集合である。
＜EMアルゴリズム＞
本アルゴリズムにおいては、等式(I)をΘ、
【００８２】
【数３９】
2005129024.tif 000040
及びσの上で最大化し、得られた最大尤度をLmaxとして算出する。また、本アルゴリズムにおいては、等式(II)をΘ、
【００８３】
【数４０】
2005129024.tif 000041
及びσの上で最大化し、得られた最大尤度をL0maxとして算出する。次に、本アルゴリズムでは、尤度比L0max/Lmaxをハプロタイプの存在とQTL表現型との関連の検定に用いる。
【００８４】
Lmaxの最大化については、推定すべきパラメータはΘ=(θ1,θ2,...,θL)、
【数４１】
2005129024.tif 000042
及びσであるが、L0maxの最大化については推定すべきパラメータはΘ=(θ1,θ2,...,θL)、
【００８５】
【数４２】
2005129024.tif 000043
及びσである。帰無仮説の下では-2log(L0max/Lmax)は自由度1のχ2分布に従う。
【００８６】
もしd1, d2,...,dNとψ1,ψ2,...,ψNに関して完全データが得られるならば、θ1,θ2,...,θL、
【数４３】
2005129024.tif 000044
及びσ最尤推定量は、それぞれ
【００８７】
【数４４】
2005129024.tif 000045
のように得ることができる。
【００８８】
ただし、j＝1,2,...,Lはハプロタイプの番号、njは個体中のj番目のハプロタイプコピーの数、N＋はハプロタイプhbを保有する個体の数、N−はハプロタイプhbを保有しない個体の数を表す。しかしながら、完全データは得られず、単に個体の遺伝子型とQTL表現型を観察できるのみである。したがって、真の値の代りに
【００８９】
【数４５】
2005129024.tif 000046
の期待値を代入する以下のEMアルゴリズム（ステップ１〜９）を構築する。
【００９０】
ステップ１
n＝0について、初期値として
【数４６】
2005129024.tif 000047
を与える。但し、
【００９１】
【数４７】
2005129024.tif 000048
であり、
【００９２】
【数４８】
2005129024.tif 000049
である。
【００９３】
ステップ２
n＝0について、初期値として
【数４９】
2005129024.tif 000050
を与える。
【００９４】
ステップ３
n＝0について、標準偏差の初期値にσ(n)を与える。σはハプロタイプhbを保有するか否かに拘わらず、同じ値を取ると仮定する。
【００９５】
ステップ４
全ての個体iについて、giと合致する順列ディプロタイプ形akに対して、
【数５０】
2005129024.tif 000051
を計算する。ここでAiは、giに合致する個体iのディプロタイプ形amの集合である。giに合致するakのみについて計算することに注意する。さらに、diは
【００９６】
【数５１】
2005129024.tif 000052
及びσ(n)と独立であり、ψiはdiの条件の下でΘ(n)と独立なので、上記式(III)は、
【００９７】
【数５２】
2005129024.tif 000053
となる。
【００９８】
ステップ５
Ｎ個体中のj番目のハプロタイプコピーの数njはランダム変数であるので、j番目のハプロタイプコピーの数の期待値が次式のように定義できる。
【００９９】
【数５３】
2005129024.tif 000054
【０１００】
ここで、gj(ak)はakの中に含まれているj番目のハプロタイプコピーの数を表し、Aiは再びgjに合致する個体iのディプロタイプ形の集合を表す。gj(ak)は、0、1、2のいずれかであることを指摘しておく。期待値は全てのjについて計算する。
【０１０１】
ステップ６
ここで、
【数５４】
2005129024.tif 000055
及び
【０１０２】
【数５５】
2005129024.tif 000056
は、ランダム変数であるので、それぞれ期待値が定義できる。
【０１０３】
まず、
【数５６】
2005129024.tif 000057
については次式で定義される。
【０１０４】
【数５７】
2005129024.tif 000058
ここで分子と分母は、giに合致する個体の順列ディプロタイプ形の集合に対する表現型に関係したハプロタイプを含む集合の割合、すなわち、ub/uoによってそれぞれ重み付けられている。
【０１０５】
次に、
【数５８】
2005129024.tif 000059
については次式で定義される。
【０１０６】
【数５９】
2005129024.tif 000060
ここで分子と分母は、giに合致する個体iの順列ディプロタイプ形の集合に対する表現型に関係したハプロタイプを含まない集合の割合、すなわち、vb/voによってそれぞれ重み付けられている。
【０１０７】
さらに、
【数６０】
2005129024.tif 000061
については次式で定義される。
【０１０８】
【数６１】
2005129024.tif 000062
ここでσは、giに合致する個体iの順列ディプロタイプ形の集合に対する表現型に関係したハプロタイプを含む集合の割合（すなわちub/uo）と、giに合致する個体iの順列ディプロタイプ形の集合に対する表現型に関係したハプロタイプを含まない集合の割合（すなわちvb/vo）とによってそれぞれ重み付けられている。また、nは
【０１０９】
【数６２】
2005129024.tif 000063
で与えられる。
【０１１０】
ステップ７
ステップ５の計算結果を用いて、次ステップの計算を行うためΘを以下のように更新する。
【数６３】
2005129024.tif 000064
ステップ6の計算結果を用いて、次ステップの計算を行うため
【０１１１】
【数６４】
2005129024.tif 000065
及びσを以下のように更新する。
【０１１２】
【数６５】
2005129024.tif 000066
【０１１３】
ステップ８
ステップ４からステップ７までを値が収束するまで繰り返す。収束した時の値は最尤推定量
【数６６】
2005129024.tif 000067
及び
【０１１４】
【数６７】
2005129024.tif 000068
とみなされる。
【０１１５】
ステップ９
局所的な極大値への収束を避けるため、異なる初期値のセット
【数６８】
2005129024.tif 000069
及びσ(0)について繰返し収束計算を行ない、それらのうちで尤度関数が最大となるものを求める。ここで、
【０１１６】
【数６９】
2005129024.tif 000070
が対立仮説における最大尤度Lmaxである。もしも
【０１１７】
【数７０】
2005129024.tif 000071
を条件として、ステップ4からステップ7までを繰り返せば、帰無仮説に対する最大尤度L0maxが得られる。帰無仮説の下では、統計量-2log(L0max/Lmax)は自由度1のχ2分布に従うことが期待される。
【０１１８】
上記のEMアルゴリズム及び量的表現型の分布を定める平均値及び標準偏差推定アルゴリズム（本アルゴリズム）は、例えばコンピュータソフトウェアに搭載することができる。本アルゴリズムをコンピュータソフトウェアに搭載することによって、全ての計算をコンピュータにおいて行うことができる。
【０１１９】
また、本アルゴリズムを搭載したソフトウェアをインストールしたコンピュータは、通信回線網を介して外部ネットワークと接続されていてもよく、例えば、コホート研究又は臨床治療試験から得られた個人の遺伝子型データ及び連続値を取る表現型データを当該通信回線網を介して取得することができる。また、本アルゴリズムによって推定された量的表現型の分布を定める平均値及び標準偏差推定を、通信回線網を介して外部に出力することもできる。
【０１２０】
＜シミュレーション＞
分布パラメータ推定精度
本アルゴリズムの推定精度をシミュレーションによって検討した。採用したシミュレーションは、ハプロタイプ頻度を仮定した集団よりＮ個体の遺伝子型及び表現型を決定して標本を作成し、得られた標本に対して推定及び検定を行う方法である。集団のハプロタイプ頻度Θは、集団遺伝学的なモデルを仮定せず、過去に我々の研究によって得られたSAA(血清アミロイド A)遺伝子に関するハプロタイプ頻度を用いた[Moriguchi et al. (2001) A novel single-nucleotide polymorphism at the 5'-flanking region of SAA1 associated with risk of type AA amyloidosis secondary to rheumatoid arthritis. Arthritis Rheum 44:1266-1272]。SAA遺伝子に関するハプロタイプデータは6座位のSNPのデータを含んでいるが、6座位のうち、座位1と座位4は連鎖不平衡があまり強くない。このため、これらの座位を含めた6座位のハプロタイプ分布と、これらの座位を除いた4座位のハプロタイプ分布を集団のハプロタイプ分布と仮定した。使用したハプロタイプ分布を表1に示す。
【０１２１】
【表１】
2005129024.tif 000072
【０１２２】
仮定したハプロタイプ頻度に基づき、順番を付けた2つのハプロタイプコピーをＮ個体のそれぞれに対して抽出することによって遺伝子型を決めた。さらに、表1に示すハプロタイプの一つをN(μ1,σ)の分布に従うQTL表現型に関係したハプロタイプと決めた。表現型に関係したハプロタイプを保有している場合にランダムにN(μ1,σ)の分布に従うQTL表現型を与え、表現型に関係したハプロタイプを保有していない場合にランダムにN(μ2,σ)の分布に従うQTL表現型を与えた。その後、相を除去して上記のアルゴリズムを適用し、分布パラメータの推定を行った。
【０１２３】
4座位のハプロタイプ分布を用いたシミュレーションデータによる検証計算の結果を表2の推定・検定結果の欄に示す。なお、推定・検定結果の欄は、異なる標本の大きさＮに対して、帰無仮説に対応するμ1＝μ2＝160、σ＝5、対立仮説に対応するμ1≠μ2、σ＝5で行った。
【０１２４】
【表２】
2005129024.tif 000073
【０１２５】
本アルゴリズムの推定精度は、シミュレーションデータにおける真の解と比較をすることによって評価することができる。シミュレーションデータを生成した時の順列ディプロタイプ形によって、ハプロタイプCCTCを保有している個体と保有していない個体に分類し、それぞれの平均値μ1及びμ2を求め、μ1及びμ2に基づいた全体の標準偏差を求めると表2の標本の統計量の欄となる。標本から得られた統計量と本アルゴリズムによる推定結果を比較すると、4座位のシミュレーションデータの条件設定においては誤差は大きくても0.3%程度であり、本アルゴリズムは非常に良い推定精度を持っていると結論される。また、特定のハプロタイプと関連する表現型の平均値の推定精度が良いことから、少なくとも1000人規模のデータでは表現型と関連するハプロタイプの検出が十分に可能であることが分かる。
【０１２６】
次に、連鎖不平衡の弱い座位を含めた場合の推定精度を検討するため、6座位のハプロタイプ分布を用いたシミュレーションデータによる検証計算の結果を表3の推定・検定結果の欄に示す。なお、推定・検定結果の欄は、異なる標本の大きさＮに対して、帰無仮説に対応するμ1＝μ2＝160、σ＝5、対立仮説に対応するμ1≠μ2、σ＝5で行った。
【０１２７】
【表３】
2005129024.tif 000074
【０１２８】
4座位の場合と同様に、ハプロタイプACCGTCを持っている個体と持っていない個体に分類し、それぞれの平均値μ1及びμ2を求め、μ1及びμ2に基づいた全体の標準偏差を求めると表3の標本の統計量の欄となる。標本から得られた統計量と本アルゴリズムによる推定結果を比較すると、6座位のシミュレーションデータの条件設定においても、特定のハプロタイプと関連する表現型は分離可能であることは明らかである。しかしながら、
【０１２９】
【数７１】
2005129024.tif 000075
が大きい場合、ここではμ1＝180やμ1＝170の場合には、標準偏差の最尤推定値
【０１３０】
【数７２】
2005129024.tif 000076
が大きくなる傾向にある。これは、本アルゴリズムで採用している標準偏差の期待値を算出する式の性質によるものである。
【０１３１】
【数７３】
2005129024.tif 000077
が大きい時には、
【０１３２】
【数７４】
2005129024.tif 000078
と
【０１３３】
【数７５】
2005129024.tif 000079
の値には大きな差が生じる。例えばμ1が正解である場合、
【０１３４】
【数７６】
2005129024.tif 000080
は
【０１３５】
【数７７】
2005129024.tif 000081
よりも大きくなる確率が大きい。ここで、ディプロタイプ形が一意に決まらないとすれば、それぞれに分配される偏差の平方の和は正解と比較して大きくなる。表1の6座位の分布では連鎖不平衡が弱い座位も含まれており、比較的小さな頻度のハプロタイプも出現していることから、ディプロタイプ形が一意に決まらない個体もある。このため、推定精度が悪くなっている。
【０１３６】
本アルゴリズムにおける標準偏差の期待値の計算方法から、平均値の差が小さい場合とディプロタイプ形が一意に決まる個体が多い場合に推定精度が高くなる。したがって、ハプロタイプブロックを同定した上で本アルゴリズムを適用する方が良いと結論される。また、計算結果のディプロタイプ形の分布を見れば、分布に広がりがある場合には推定精度が悪くなっていると判断することが可能である。 
【０１３７】
ディプロタイプ形推定精度
表4は4座位の検証計算結果より、μ1＝165、μ2＝160、σ＝5.0、Ｎ＝1000のケースにおいて、最初10個体のハプロタイプ頻度分布の事後確率であるディプロタイプ分布を示している。
【０１３８】
【表４】
2005129024.tif 000082
【０１３９】
4座位のデータではディプロタイプ形が集中しているため、両者の差は小さい。ただし、詳細を見れば、EMアルゴリズムにより推定されたハプロタイプ頻度分布を用いた事後確率と、QTL表現型を考慮した事後確率を比較すると、個体2はハプロタイプCCTCを持つ個体2では後者が大きく計算されている。
【０１４０】
表5は6座位の検証計算結果より、上記の4座位と同じ条件であるμ1＝165、μ2＝160、σ＝5.0、Ｎ＝1000のケースにおいて、51番目から60番目の10個体のハプロタイプ頻度分布の事後確率であるディプロタイプ分布を示している。
【０１４１】
【表５】
2005129024.tif 000083
【０１４２】
EMアルゴリズムにより推定されたハプロタイプ頻度分布を用いた事後確率と、QTL表現型を考慮した事後確率を比較すると、ハプロタイプACCGTCを持つ個体55では後者が大きく計算されていることが判る。
【０１４３】
表４及び５に示したディプロタイプ形の分布推定の結果は、表現型と関連があるハプロタイプを持つ場合にディプロタイプ形の推定精度は高くなることを示している。これは、表現型とハプロタイプの関連という遺伝的な効果をモデル化した本アルゴリズムが、個体の遺伝情報の推定精度を高めるということを示すものである。
【０１４４】
帰無仮説の下での統計量-2log(L0max/Lmax)の経験的分布:
先ず、統計量-2log(L0max/Lmax)の経験的分布を帰無仮説の下でのシミュレーションにより検討した。帰無仮説はμ1＝μ2とすることに対応する。統計量-2log(L0max/Lmax)は、標本一つに対して一つの値が定まる。シミュレーションにより多数の標本を生成して、この統計量の分布を調べることにより、経験的な分布が得られる。
【０１４５】
4座位のハプロタイプ分布を用いたミュレーションによって得られた検定統計量-2log(L0max/Lmax)の分布に関するヒストグラムを図1に示す。図1は検定統計量がほぼ自由度1のχ２分布に従っていることを示している。なお、図１における左のヒストグラムはＮ＝100、標本数＝10000であり、右のヒストグラムはＮ＝1000、標本数＝10000である。図１中、ヒストグラムは棒グラフで示され、自由度1のχ２分布の確率密度関数は曲線で示されている。
【０１４６】
また、推定された分布パラメータとタイプIエラーを表6に示す。なお、ここでは帰無仮説に対応するμ1＝μ2＝160、σ＝5の条件下で、異なる標本の大きさＮ＝100又はＮ＝1000の標本を繰り返し抽出するシミュレーションを行った。それぞれのパラメータの組に対して、10000回の標本抽出を行った。
【０１４７】
【表６】
2005129024.tif 000084
【０１４８】
タイプIエラーは5%の有意水準に対し、Ｎ＝100、標本数10000のケースでは4.96%、Ｎ＝1000、標本数10000のケースでは5.14%といずれも5%に近い値となっている。また、6座位のハプロタイプ分布を用いたミュレーションによって得られた推定された分布パラメータとタイプIエラーも表6に示す。タイプIエラーは5%の有意水準に対し、Ｎ＝100、標本数10000のケースでは6.06%、Ｎ＝1000、標本数10000のケースでは5.41%といずれも5%に近いが、大きめの値となっている。ディプロタイプ形が集中しない場合には、検定に関してもタイプIエラーが僅かではあるが増加する傾向であることを示すものである。
【０１４９】
したがって、本アルゴリズムにおいては、検定を行う場合にもディプロタイプ形が一意に決まる個体が多い場合に精度が高くなると考えられる。このため、ハプロタイプブロックを同定した上で本アルゴリズムを適用する方が良いと結論される。
【０１５０】
検出力の評価
次に、対立仮説の下でのシミュレーションを行い、検出力の評価を行った。μ1＝μ2＋γσとしてγを変化させることにより、γに対するタイプIIエラー、したがって検出力が評価される。シミュレーションにより、それぞれのγに対して標本を生成し、検出力の評価を行った。
【０１５１】
すなわち、μ1＝μ2＋γσとしてγを変化させて、検定統計量-2log(L0max/Lmax)の分布を評価することにより、検出力を評価することができる。4座位のハプロタイプ分布を用いたミュレーションにおいて、γ＝（μ1−μ２）/σを変化させて検出力を評価した結果を表7に示す。
【０１５２】
【表７】
2005129024.tif 000085
【０１５３】
また、γ＝（μ1−μ２）/σに対する検出力を図2に示す。Ｎ＝100の場合（図２中実線で示す。）、μ1−μ２＝５、すなわちγ＝（μ1−μ２）/σ＝1.0であれば検出力はほぼ1.0である。Ｎ＝1000の場合（図２中破線で示す。）、μ1−μ２＝1.3、すなわちγ＝（μ1−μ２）/σ＝0.3であれば検出力はほぼ1.0である。したがって、本アルゴリズムの検出力は十分に高いと考えられる。
【０１５４】
＜現実データの分析＞
CAPN10遺伝子と糖尿病との関係についての研究(Horikawa Y, Oda N, Cox NJ, Li X, Orho-Melander M, Hara M, Hinokio Y et al (2000) Genetic variation in the gene encoding calpain-10 is associated with type 2 diabetes mellitus. Nat Genet 26:163-175. )によって得られている3SNPsの遺伝子型データと検査データを用いて、本アルゴリズムによる解析を実施した。CAPN10遺伝子は2型糖尿病に関連する遺伝子であり、BMI、血糖値など関連する検査データが取得されている。図3は、当該研究によって得られている検査データの分布を示している。それぞれ、正規分布で近似可能な分布となっている。
【０１５５】
解析の第一段階として、本アルゴリズムのハプロタイプ推定機能を用いてハプロタイプ頻度を求めた。表8はSNPs頻度とHWE検定のP値を示しており、表9はハプロタイプ頻度を示している。また、表10は各座位間の連鎖不平衡尺度を示している。
【０１５６】
【表８】
2005129024.tif 000086
【０１５７】
【表９】
2005129024.tif 000087
【０１５８】
【表１０】
2005129024.tif 000088
【０１５９】
次に、ハプロタイプ推定の結果より、表現型に関係したハプロタイプとして、121、112、111、221、122、122及び212を仮定し、本アルゴリズムを用いた解析を実施した。表11に示す解析結果より、BS 30'およびBS 60'でハプロタイプ112が有意であることが示されており、BS 0'でハプロタイプ122が有意であることが示されている。また、表12に示す推定結果は、BS 30'およびBS 60'とハプロタイプ112との組合せで
【０１６０】
【数７８】
2005129024.tif 000089
であることが示されており、ハプロタイプ112がBS 30'およびBS 60'を有意に上昇されるという結果となった。一方、BS 0'とハプロタイプ122の推定結果は、
【０１６１】
【数７９】
2005129024.tif 000090
となっている。
【０１６２】
【表１１】
2005129024.tif 000091
【０１６３】
【表１２】
2005129024.tif 000092
【０１６４】
もともとハプロタイプ122の頻度が小さいため、該当する標本の数が少なかったことが検定結果を有意とした原因と考えられる。
【０１６５】
以上、現実データを用いた分析の結果より、本アルゴリズムは、遺伝子型データと連続値を取る表現形データとを用いて、ハプロタイプ頻度に加えて、それぞれ個体のディプロタイプ形の事後確率分布並びにQTL表現型パラメータを最尤推定することができる。これにより、本アルゴリズムによれは、遺伝子型データと連続値を取る表現型データより、個体のレベルでQTL表現型とハプロタイプの関連を検定することができる。また、本アルゴリズムによれば、特定のハプロタイプのあるなしにより異なった範囲で量的表現型を取るという仮定の下で、当該範囲を最尤推定することもでき、更に、それらの範囲から当然相対危険の最尤推定値も得ることができる。
【０１６６】
一方、上述したように、不完全ハプロタイプHlを検定対象として、CAPN10遺伝子と糖尿病との関係についての研究によって得られている3SNPsの遺伝子型データと検査データを用いた解析を実施した。このデータは、検査データであるBS 30’及びBS 60’において、ハプロタイプ112が有意であることが既に示されている（表８乃至１２参照）。
【０１６７】
不完全ハプロタイプHlを用いたBS 30’に対する解析結果を表１３に示す。
【０１６８】
【表１３】
2005129024.tif 000093
【０１６９】
表１３に示すように、解析結果において、最も低いP値を示しているのはハプロタイプ「122」であり、このハプロタイプが有意であることが示されている。ただし、不完全ハプロタイプを構築するhtSNPを定める際に、頻度の小さい「122」ハプロタイプは除外しており、「122」ハプロタイプは検定対象から除外されている。このように、不完全ハプロタイプを用いることによって、「表現型に関連したハプロタイプ」を探索できることが可能となった。
【０１７０】
さらに、不完全ハプロタイプHlを用いた他の適用例として、肝機能の指標であるGPTが45を超えた症例を示した98人の関節リウマチ患者に関して、CHST3遺伝子に関する4つのSNPからなる遺伝子型データと、GPTが45を超えた時点までにメソトレキサートの累積投与量のデータに対して、検定対象に不完全ハプロタイプを用いた解析を実施した。その結果を表１４に示す。
【０１７１】
【表１４】
2005129024.tif 000094
【０１７２】
解析結果において、最も低いP値を示しているのはハプロタイプ「*AT*」であり、全ての座位の情報を用いたハプロタイプ「AATC」あるいは「TATT」よりもより有意であることが示されている。これは、第二座位と第三座位による不完全ハプロタイプ「AT」が、表現型をとる原因とより強い関連があることを示している。すなわち、原因座位は「AATC」及び「TATT」の両方のハプロタイプと関連しており、第一座位及び第四座位をマスクすることにより、表現型との関連が強く表われたものと理解することができる。この結果は、不完全ハプロタイプを用いることにより、複数のハプロタイプと関連する表現型に対しても検出することが可能となることを示すものである。
【図面の簡単な説明】
【０１７３】
【図１】4座位のハプロタイプ分布を用いたミュレーションによって得られた検定統計量-2log(L0max/Lmax)の分布に関するヒストグラムである。
【図２】γ＝（μ1−μ２）/σに対する検出力を示す特性図である。
【図３】CAPN10遺伝子と糖尿病との関係についての研究によって得られている検査データの分布を示す特性図である。
【図１】
2005129024.tif 000139
【図２】
2005129024.tif 000140
【図３】
2005129024.tif 000141