遺伝子関連解析装置及び遺伝子関連解析プログラム
【課題】稀なマーカーの頻度差についても検討可能な小標本理論に基づいた検出方法を提供する。
【解決手段】２以上の座位に関する遺伝子型データを含む個体に関する遺伝子型データセット、並びに、当該個体に関する表現型データを入力する入力手段と、
サンプリング法により上記遺伝子型データセットから複数のサンプリングデータセットを作成し、作成されたサンプリングデータセットに含まれる遺伝子型データセットを用いて当該サンプリングデータセットにおけるハプロタイプ頻度を最尤推定し、得られたハプロタイプ頻度から各個体のディプロタイプ形事後確率を計算し、計算された事後確率からディプロタイプ形をサンプリングし、サンプリングされたディプロタイプ形からハプロタイプ及び/又はディプロタイプ形の分割表を作成し独立性検定を行うことでp値を求め、算出されたp値からp値の分布を算出する処理を行う演算手段と、
上記演算手段で算出されたp値の分布に基づいて当該分布の統計量を出力する出力手段とを備える。
【解決手段】２以上の座位に関する遺伝子型データを含む個体に関する遺伝子型データセット、並びに、当該個体に関する表現型データを入力する入力手段と、
サンプリング法により上記遺伝子型データセットから複数のサンプリングデータセットを作成し、作成されたサンプリングデータセットに含まれる遺伝子型データセットを用いて当該サンプリングデータセットにおけるハプロタイプ頻度を最尤推定し、得られたハプロタイプ頻度から各個体のディプロタイプ形事後確率を計算し、計算された事後確率からディプロタイプ形をサンプリングし、サンプリングされたディプロタイプ形からハプロタイプ及び/又はディプロタイプ形の分割表を作成し独立性検定を行うことでp値を求め、算出されたp値からp値の分布を算出する処理を行う演算手段と、
上記演算手段で算出されたp値の分布に基づいて当該分布の統計量を出力する出力手段とを備える。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、集団から取得した遺伝子型データ及び表現型データを用いて関連解析を行う遺伝子関連解析装置及び遺伝子関連解析プログラムに関する。
【背景技術】
【０００２】
表現型を記述するゲノム領域を特定するためにSNPsを用いた関連解析が数多く行われている。現在は個人当り55万SNPsを一度に測定する技術なども開発されており、大量の多型情報を用いた高密度SNPs解析が主流となっている。高密度SNPs解析ではSNPs間の遺伝的距離が短く、連鎖不平衡があるためにマーカーであるSNPsの情報が独立ではないという特徴がある。そのため実際の関連解析においてはこのような連鎖不平衡を考慮した解析が求められる。
【０００３】
一方でSNPは通常2つの対立因子(allele)しか持たない。そのため、ある座位に関するSNPを用いた2つの表現型間の独立性検定では、alleleやdominant model、recessive modelの元での2×2の分割表や、遺伝子型を使った2×3の分割表などが通常用いられる。表現型への寄与がある程度あり、比較的頻度の高いSNPsの検出には、このような解析は非常に有効であるが、
1.様々な因子の効果により1つの座位の影響が小さい場合
2.頻度の小さな原因多型
の場合にはこれらを検出することは困難であると考えられる。
【０００４】
1つの座位の影響が小さい場合には、環境要因や複数の遺伝子間の相互作用等が考えられる。これらの解析には環境要因を記述する変数の導入や、 複数のマーカーとの関連などを考慮した解析が必要になり、検出は非常に困難になることが予想される。
【０００５】
マーカーに用いるSNPsセットに常に目的の表現型の原因となるSNPsが含まれていると仮定することは出来ない。元々マーカーとしてのSNPsは、マーカーとしての有効性のために頻度の小さいSNPsは選ばれていない。そのため表現型の寄与度が大きくとも、頻度が低ければこれがSNPマーカーとなることはない。一方で頻度の小さいSNPは、たとえ周辺領域との連鎖不平衡が高くとも、頻度の高い周辺の1つのSNPで記述することは出来ない。このような状況の例としては、common disease/multiple rare variants(CD/MRVs)や、非常に稀な副作用を起こす薬剤に対してこの原因を記述するゲノム上の変異を検出する場合等が考えられる。このような非常に稀な多型(例えば頻度0.1%)を効率的に検出するためには、 頻度の低いマーカーが必要になる。頻度の低いマーカーとしては、稀なSNPsを 網羅的にタイピングすることや、マイクロサテライトマーカーを用いることが考えられるが、タイピングコストを考えると決して効率的な方法であるとは言えない。
【０００６】
そこで頻度の高いSNPsを複数個使用することで対立因子数を増やすことが考えられる。目的となる感受性変異がマーカーの変異よりも後で起った場合を想定すると、周辺のSNPや複数のSNPsを組み合わせて検出するよりもこの原因多型を挟んだ複数のマーカーによるハプロタイプやディプロタイプを用いた解析の方がより有意差を検出しやすく、有効な解析であると考えられる。
【０００７】
これまで大規模な家系を用いた遺伝統計学的手法では、家系情報からハプロタイプを特定し、これをマーカーとして独立性検定を行ってきた。しかしながら一般集団を用いた関連解析では通常各個人のハプロタイプを同定することは出来ない。そのため集団の情報を使って推定する方法論が盛んに研究されてきた。このような方法論の延長としてハプロタイプを用いた関連解析も研究されており、各表現型毎に最尤法によって求められたハプロタイプ頻度を使った尤度比検定などの検出方法論が研究されている。
【発明の開示】
【発明が解決しようとする課題】
【０００８】
しかしながら、尤度比検定は大標本理論に基づいており、稀な頻度差を検討するには不向きな方法論である。そこで、本発明は、稀なマーカーの頻度差についても検討可能な小標本理論に基づいた検出方法を提供することを目的とする。
【課題を解決するための手段】
【０００９】
上述した目的を達成するため、本発明者が鋭意検討した結果、小標本に適用可能な検定方法としては、分割表を用いたFisher's exact testを挙げられるが、ハプロタイプやディプロタイプは通常の遺伝子型データから直接求めることが出来ないため、直接Fisher's exact test を行うことは出来ないといった問題を見いだした。そして、この問題を解決するため、ハプロタイプやディプロタイプを推定するとともに、サンプリングによる揺らぎ及び相の推定による揺らぎを考慮した遺伝子関連解析方法を開発することに成功した。
【００１０】
本発明は以下を包含する。
（１）２以上の座位に関する遺伝子型データを含む個体に関する遺伝子型データセット、並びに、当該個体に関する表現型データを入力する入力手段と、サンプリング法により上記遺伝子型データセットから複数のサンプリングデータセットを作成し、作成されたサンプリングデータセットに含まれる遺伝子型データセットを用いて当該サンプリングデータセットにおけるハプロタイプ頻度を最尤推定し、得られたハプロタイプ頻度から各個体のディプロタイプ形事後確率を計算し、計算された事後確率からディプロタイプ形をサンプリングし、サンプリングされたディプロタイプ形からハプロタイプ及び/又はディプロタイプ形の分割表を作成し独立性検定を行うことでp値を求め、算出されたp値からp値の分布を算出する処理を行う演算手段と、上記演算手段で算出されたp値の分布に基づいて当該分布の統計量を出力する出力手段とを備える遺伝子関連解析装置。
【００１１】
（２）上記出力手段は、上記分布の統計量として、少なくともp値の信頼区間の上限値及び平均値を出力することを特徴とする（１）記載の遺伝子関連解析装置。
【００１２】
（３）上記演算手段は、上記サンプリング法としてブートストラップ法を適用してサンプリングデータセットを作成することを特徴とする（１）記載の遺伝子関連解析装置。
【００１３】
（４）上記演算手段は、上記独立性検定として少なくともフィッシャー正確性検定を適用してｐ値を算出することを特徴とする（１）記載の遺伝子関連解析装置。
【００１４】
（５）入力手段、演算手段及び出力手段を備えるコンピュータに、上記入力手段が、２以上の座位に関する遺伝子型データを含む個体に関する遺伝子型データセット、並びに、当該個体に関する表現型データを入力するステップと、上記演算手段が、サンプリング法により上記遺伝子型データセットから複数のサンプリングデータセットを作成し、作成されたサンプリングデータセットに含まれる遺伝子型データセットを用いて当該サンプリングデータセットにおけるハプロタイプ頻度を最尤推定し、得られたハプロタイプ頻度から各個体のディプロタイプ形事後確率を計算し、計算された事後確率からディプロタイプ形をサンプリングし、サンプリングされたディプロタイプ形からハプロタイプ及び/又はディプロタイプ形の分割表を作成し独立性検定を行うことでp値を求め、算出されたp値からp値の分布を算出する処理を行うステップと、上記出力手段が上記演算手段で算出されたp値の分布に基づいて当該分布の統計量を出力するステップとを実行させる遺伝子関連解析プログラム。
【００１５】
（６）上記出力手段は、上記分布の統計量として、少なくともp値の信頼区間の上限値及び平均値を出力することを特徴とする特徴とする（５）記載の遺伝子関連解析プログラム。
【００１６】
（７）上記演算手段は、上記サンプリング法としてブートストラップ法を適用してサンプリングデータセットを作成することを特徴とする（５）記載の遺伝子関連解析プログラム。
【００１７】
（８）上記演算手段は、上記独立性検定として少なくともフィッシャー正確性検定を適用してｐ値を算出することを特徴とする（５）記載の遺伝子関連解析プログラム。
【発明の効果】
【００１８】
本発明に係る遺伝子関連解析装置及び遺伝子関連解析プログラムによれば、原因となる多型の頻度が非常に低い場合であっても、表現型との関連を検出することができる。
【発明を実施するための最良の形態】
【００１９】
以下、本発明を詳細に説明する。
１. ハプロタイプやディプロタイプを用いた関連解析の有効性
１−１. CD/CV vs CD/MRVs
晩成型のcommon diseaseの原因においてSNPsの様なgerm-lineからの寄与が小さくないと考えられる場合には、大きくcommon disease/Common variant(CD/CV)仮説と common disease/multiple rare variants(CD/MRVs)仮説が考えられる。一般集団を用いたゲノムワイドなSNPs解析は、主にCD/CV仮説に従う変異を抽出することを目的としている。しかしながら、突然変異率や集団での遺伝的浮動を考えるとCD/CV仮説に従う変異はそれ程多く無く、むしろCD/MRVs仮説に従う多数の稀な変異が多い可能性を示唆している。このような稀な変異は、これまで主に伝統的な連鎖解析によって抽出されてきた。
疾患に寄与する変異を考える場合、表現型への寄与度と集団中での頻度とをパラメータとして考えることができる。各場合の検出方法論について表1にまとめた。
【００２０】
【表１】
【００２１】
頻度が小さく、寄与度も小さい原因変異を特定するには非常に多数のサンプルが必要になると考えられ、これを検出することは困難であると考えられる。一方でこれまで連鎖解析の対象であった頻度が低く、寄与度が大きい原因変異をcommon SNPsを用いたcase-control studyで検出することができる可能性がある。一般に表現型原因変異とマーカーとのr2(後述する参考例A参照)が小さくなると検出力は劇的に下がる。そこで頻度の小さい変異を検出するためには、
【００２２】
【数１】
となるような頻度の小さいマーカーが必要になる。Common SNPsから効率的に頻度の小さいマーカーを作成する方法としてハプロタイプを利用することが考えられる。
【００２３】
１−２.関連解析の検出力
本節ではマーカー密度が非常に高い場合における検出力について説明する。まずは原因変異とマーカーが一致していないために起る見掛けのオッズ比(マーカーオッズ比)と検出力について説明し、続いて大標本理論に基づく検出力の近似式(後述する参考例B参照)を用いて様々なパラメータに対する検出力の傾向を検討する。最後にsimulation dataを使った検出力について説明し、ハプロタイプを使った関連解析がどのような場合に有効であるのかを説明する。
【００２４】
１−２−１. 見掛けのオッズ比
表現型の原因多型が稀な場合、例えオッズ比が高くとも所謂common SNPsでは、これを検出することは非常に困難になる。例えばオッズ比が5となるような変異であっても集団中の頻度が1%であればcaseでの頻度は4.8%程度にしかならない。この3.8%の頻度差を頻度10%のSNPsで検出する場合、見掛けのオッズ比は最大(D'=1)でも
【００２５】
【数２】
となってしまう。よって、稀な変異を単座位のcommon SNPで解析する場合、影響力の大きい稀なSNPが原因であったとしても見掛けのオッズ比が小さくなり、これを検出するにはより多くのサンプルが必要となる。頻度の小さいマーカーで記述したときの分割表を表2に示し、common SNPで記述したときの分割表を表3に示す。上記の状況で表現型の原因となる稀なSNPをタイピングしていた場合には、表2のような分割表が期待されP=0.019となる。しかしながら、このSNPとD'=1となる様な頻度10%のSNPでも期待値としては表3の様な分割表となり、P=0.218となってしまう。
【００２６】
【表２】
【００２７】
【表３】
【００２８】
表３に示した分割表のように、稀な変異を単座位のcommon SNPで解析する場合にはp値が大きくなり帰無仮説を棄却できない。すなわち、単座位のcommon SNPと表現型との関連性を検出することができない。
【００２９】
１−２−２. D'=1モデル
表現型の原因となる変異の割合がcaseとcontrolの間で異なるモデルを考える。 ただしマーカーの密度は十分に高く、原因変異とD'=1となる場合のみを考える。この場合、原因変異とマーカーとの関係は図1に示す2通りが考えられる。ここで表現型の原因変異をD、対立因子をdと表し、 この変異とD'=1の関係にあるマーカーのalleleをM及びmで表す。また集団中でのこれらのalleleの頻度をそれぞれfD、fd、fM及びfmと表現することにし、ハプロタイプ頻度をそれぞれhDM、hDm、hdM及びhdmとする。以下では図１に示したタイプ１とタイプ２におけるマーカーの頻度を導き、このマーカーの検出力を説明する。
【００３０】
タイプ１
一般集団中では、マーカーの存在比はallele頻度に従いfM：fmであるが、caseにおいては原因変異Dの頻度の増加に伴いマーカーの頻度にも歪みが生じる。原因変異の割合の増分をδと書くと、原因変異の比率はfD＋δ：fdとなり、マーカーの存在比はfM＋δ：fmとなる。この時、caseにおけるマーカーMの頻度f'Mは規格化を考慮すると
【００３１】
【数３】
となり、原因変異によりもたらされるマーカー頻度の増分は
【００３２】
【数４】
と表されることになる。式(3)から、規格化を考慮するとマーカーの頻度差の増分は集団中(control)のマーカー頻度に依存することがわかる。式(3)では集団中のマーカーの頻度が増大するに従って頻度差が線形に減少している。
【００３３】
また式(3)は原因変異の頻度fDに依存していない。これは様々な効果をδに押し込めているためで、実際にはδはfdとオッズ比ORから決まるパラメータである。簡単な計算から
【００３４】
【数５】
となり、式(3)は、
【００３５】
【数６】
と表現される。
【００３６】
タイプ２
原因変異がcommonの場合には、マーカーよりも原因変異が古い場合も想定される。このとき、集団中では、fD：fdの比であった原因変異の割合が、case中ではfD＋δ：fdとなったとする。このδの増分がハプロタイプ頻度に従ってマーカーの頻度に分配されたとすれば、hDMの頻度は
【００３７】
【数７】
となる。よって、よってマーカーMのcaseにおける頻度は、hDM＝fMであるから
【００３８】
【数８】
となり、頻度差は
【００３９】
【数９】
となる。式(4)を用いれば式(8)は
【００４０】
【数１０】
となる。式(5)と式(9)を比較すると、タイプ１ではfMが増加すると頻度差が減少するが、タイプ２ではfDが増加すると頻度差が減少する。ただしtype 1では
【００４１】
【数１１】
であるが、タイプ２では
【００４２】
【数１２】
である。即ち両ケースともfM=fDが頻度差の上限となる。
【００４３】
１−２−３.D'=1モデルの検出力
検出力は原因変異の頻度とオッズ比、及びマーカー頻度に依存する。オッズ比をOR=1.5及びOR=3.0に固定した場合のタイプ1、タイプ2の場合における検出力を図2及び図3に示した。ただしこの検出力は２×２の分割表におけるχ2検定の近似式を用いた場合である(後述する参考例B参照)。fD=fMの条件の場合(即ち原因変異とマーカーとの間でr2=1)に検出力が最大になり、タイプ1の場合には
【００４４】
【数１３】
タイプ2の場合には
【００４５】
【数１４】
であるためこの曲線の内側に各検出力曲線が描かれる。なお、図2及び図3において、(A)は3次元プロットを示し、(B)は2次元プロットを示している。図2及び図3においては、原因頻度を0.01から0.05刻みで表示している(caseを100人、controlを100人とした)。
【００４６】
また、原因変異のオッズ比を3.0とした場合のマーカーMの見掛けのオッズ比とマーカーの頻度差f'M-fMをそれぞれ図4及び図5に示した。これらの曲線から原因変異とマーカーとの間の連鎖不平衡がD'=1であっても、頻度が異なれば急激に検出力が落ちることがわかる。これまでの知見通り、原因変異の頻度が0.5付近でfD=fMの場合に検出力は最大となるが、
【００４７】
【数１５】
の場合には必ずしもマーカーの頻度が0.5付近で良いと言う訳ではなく、r2=1の様なマーカーが最も検出力が高い。そのため様々な頻度のマーカーを用意することが望ましいが、現実的には難しい。このためcommon variantの検出においてもハプロタイプの方が有利である状況は十分に考えられ、ハプロタイプを使った関連解析はcommon variant及びrare variantsの検出を試みる上で検討する価値はあると考えられる。なお、図4及び図5において、(A)は3次元プロットを示し、(B)は2次元プロットを示している。図4及び図5においては、原因頻度を0.01から0.05刻みで表示している(caseを100人、controlを100人とした)。
【００４８】
１−２−４. 注意点
以上の説明からD'=1のマーカーであったとしても原因変異とマーカーの頻度が異なれば検出力は急激に落ちるため、ハプロタイプを使った解析は非常に有効である。しかしながら、ここでハプロタイプを使った場合に検出力を下げる要因を挙げておく。
1.前節の検出力の説明は、2×2の分割表における大標本理論を用いた場合の結果である。ハプロタイプを用いた場合には2×nの分割表になり、このため自由度が大きくなるために検出力が落ちる。
2.遺伝子型データでは相が分かっていないために生じる推定誤差により検出力が下がる。
3.分割表のカラムを幾つか統合検定やdominantやrecessiveのモデルを 導入した場合、多数の組み合わせが発生し新たな多重性が生じる。
4.必ずしも稀な原因変異が稀なハプロタイプに乗っているとは限らない。頻度1%のハプロタイプに原因変異がのっている確率はやはり1%であろう。この場合にはハプロタイプを構成する座位数を増やすことが考えられるが、alleleの数が多くなりすぎると自由度が大きくなり検出力は落ちる。また各個人のディプロタイプ形の推定も困難になる。よってどこまでハプロタイプを伸ばせば良いのかの議論は別に必要になる。この場合、連鎖不平衡マップなどが役に立つと考えられる。
【００４９】
これまでの解析は大標本理論に基づく２×２の分割表の結果である。そこでより現実的な２×ｎの分割表となるシミュレーションデータを使って検出力を調べた結果を次節に示した。
【００５０】
１−３.シミュレーションデータを使った検出力
１−３−１. Case及びcontrolにおける母数の作成
実際のSNPsデータから9座位分のハプロタイプ頻度を推定し、これをcontrolのハプロタイプ頻度母数とした(表4参照)。
【００５１】
【表４】
【００５２】
またこの時の連鎖不平衡指標r2及びD'を表5及び図6に示した。なお、図６において上三角部分がr2であり、下三角部分がD'である。
【００５３】
【表５】
【００５４】
図６から、3座位目から6座位目までの4座位間でD'の値が非常に高くブロックを形成していると考えられる。この4座位分のハプロタイプ頻度を表6に示した。
【００５５】
【表６】
【００５６】
この4座位のハプロタイプにおいて、2番目、3番目に頻度の高いハプロタイプにそれぞれδ=0.15を与えてcase-1、case-2のハプロタイプ頻度母数を作成した。また希なハプロタイプの影響を調べるために、表4における11番目のハプロタイプにやはりδ=0.15を与えてcase-3のハプロタイプ頻度母数を作成した。この結果得られたハプロタイプ頻度及びアレル頻度も表4に示し、これを図7に図示した。なお、図７(A)はハプロタイプ頻度母数を示し、図７(B)はアレル頻度母数を示している。
【００５７】
１−３−２.シミュレーションデータの作成
表4で示した母数を使って、それぞれ100人づつ100ケースランダムにサンプリングを行った。ただしcontrolに対してはcontrol対controlの比較を行うためにcontrol-1及びcontrol-2としてそれぞれ100ケース、計200ケースのサンプリングを行った。このサンプリング方法では、相乗モデル（multiplicative model）となり、caseでもハーディ・ワインベルグ平衡（Hardy-Weinberg equilibrium）が成立していることになる。尚、各ハプロタイプ頻度、アレル頻度の平均、分散は多項分布に従うため、分散はnpi(1-pi)となる(ここでnはサンプリング数であり、piは頻度を意味する)。
【００５８】
１−３−３.シミュレーションデータの独立性検定の結果
各ケースにおける独立性検定の結果を、ボックス図と有意となった数のヒストグラムとで示した。座位毎の特異性を見るためにアレルについては各座位毎の結果を示し、またハプロタイプは3〜9座位の長さで構成し、これを、1座位ずつwindowをずらして計算した結果を示した。
【００５９】
Control-1 vs control-2 (100人対100人)
Control-1対control-2で100ケースずつサンプリングした時の単座位のアレル、及び3〜9座位のハプロタイプによる分割表の独立性検定(Fisher's exact test)のp値の分布を図8-1及び図8-2にボックス図として示した。また、有意水準を1%及び5%とした時に有意となった数のヒストグラムを図9-1及び図9-2に示した。Control対controlの比較ではp値の分布は一様分布となるが、図8-1及び図8-2は大体対称に分布しており、また平均も0.5付近にある。偽陽性の数も多少の増減はあるがほぼ期待通りの結果となっている。図8-1及び図8-2において、(A)はアレルについてのp値を示し、（B）〜(H)はそれぞれ3〜9座位のハプロタイプについてのp値を示している。図9-1及び図9-2において、(A)はアレルについて有意となった数を示し、（B）〜(H)はそれぞれ3〜9座位のハプロタイプについて有意となった数を示している。
【００６０】
Case-1 vs control-1 (15% ハプロタイプ, 100人対100人)
上記と同様に、Case-1対control-1の場合におけるp値の分布を図10-1及び図10-2にボックス図として示した。また、有意水準を1%及び5%とした時に有意となった数のヒストグラムを図11-1及び図11-2に示した。図10-1及び図10-2において、(A)はアレルについてのp値を示し、（B）〜(H)はそれぞれ3〜9座位のハプロタイプについてのp値を示している。図11-1及び図11-2において、(A)はアレルについて有意となった数を示し、（B）〜(H)はそれぞれ3〜9座位のハプロタイプについて有意となった数を示している。
【００６１】
Case-2 vs control-1 (10% ハプロタイプ, 100人対100人)
上記と同様に、Case-2対control-1の場合におけるp値の分布を図12-1及び図12-2にボックス図として示した。また、有意水準を1%及び5%とした時に有意となった数のヒストグラムを図13-1及び図13-2に示した。図12-1及び図12-2において、(A)はアレルについてのp値を示し、（B）〜(H)はそれぞれ3〜9座位のハプロタイプについてのp値を示している。図13-1及び図13-2において、(A)はアレルについて有意となった数を示し、（B）〜(H)はそれぞれ3〜9座位のハプロタイプについて有意となった数を示している。
【００６２】
Case-3 vs control-1 (1.5% ハプロタイプ, 10人対100人)
上記と同様に、Case-3対control-1の場合におけるp値の分布を図14-1及び図14-2にボックス図として示した。また、有意水準を1%及び5%とした時に有意となった数のヒストグラムを図15-1及び図15-2に示した。図14-1及び図14-2において、(A)はアレルについてのp値を示し、（B）〜(H)はそれぞれ3〜9座位のハプロタイプについてのp値を示している。図15-1及び図15-2において、(A)はアレルについて有意となった数を示し、（B）〜(H)はそれぞれ3〜9座位のハプロタイプについて有意となった数を示している。
【００６３】
１−３−４.シミュレーションデータのまとめ
シミュレーションデータを使って単座位のみを用いた独立性検定とハプロタイプによる分割表の独立性検定のp値の分布について調べた（図8、10、12及び14）。case対controlの比較において有意となった数は、連鎖不平衡の構造に強く依存しており、頻度15%のハプロタイプに変異を与えた場合はハプロタイプの方が有意になった割合が多く、頻度10%では単座位の検定において突出して有意となった数が多い座位があった。これはcase-2の場合に、5番目の座位のアレル頻度差が4座位分のハプロタイプ頻度差と同じ頻度差(約12%)となっている特殊事情のためである(即ち、この1つのアレルがこのハプロタイプのtag SNPとなっている)。
一方でcase-3では単座位の検定で有意となった割合は非常に小さく、偽陽性の量と余り変わらない。しかしながらハプロタイプを用いると有意となる数が非常に多くなることがわかる。この結果は、希なハプロタイプ変異の特徴である。
【００６４】
検出力
次に、実際の検出過程を想定して、各データ毎に有意と判定された数の内訳を調べた（図9、11、13及び15）。ここでは、有意水準に対してBonferroni補正を行い、アレル及びハプロタイプを用いた検定で共に有意差が検出されなかった場合の数、アレルの検定のみで検出された数、ハプロタイプを用いた検定のみで検出された数、アレル及びハプロタイプを用いた検定の両方で有意となった数をそれぞれ表7に示した。
【００６５】
【表７】
【００６６】
表７に示すように、Case-1の場合において、有意水準5%でハプロタイプでのみ有意となった数がアレルのみで有意となった数を若干上回っているが、有意水準1%では単座位の検定の方がやや有利となっている。Case-2の場合において、上述の特殊事情により、単座位での検定の方が検出力が高い。前節の説明では頻度の小さいハプロタイプに揺らぎを与えた場合、common SNPsに対してハプロタイプを用いた検定が有利になるはずであったがこのケースでは逆の結果となった。しかしながらハプロタイプを用いて検出される数(表中のハプロタイプとbothを足した数)はcase-1に比べて増えており、頻度の小さなハプロタイプに頻度差がある場合の方が検出力が高くなる傾向は再現している。この様に、単座位、ハプロタイプを用いた検定の有利/不利は個別の事情を反映しており単純ではない。ただし、傾向としては頻度の小さなハプロタイプに揺らぎがある場合に、common SNPsに対してハプロタイプがより有利になることが多い。実際希なハプロタイプの場合であるcase-3では、圧倒的にハプロタイプが有利となっている。
【００６７】
ハプロタイプの座位数
Case-1及びcase-2では4座位分(3〜6座位)のハプロタイプにδ=0.15の揺らぎを与えてcaseのハプロタイプ頻度母数を作成したが、検出力は必ずしも4座位のハプロタイプが高い訳ではなかった。検出力はハプロタイプの長さに対して比較的ロバストであると考えられる。一方で、Case-3では9座位全体のハプロタイプに揺らぎを与えているが、図14及び図15に示すように、前の方(1〜6番目の座位程度まで)が小さなp値が多く、有意となる数も多い。Bonferroni補正後の検出された割合でみると、有意水準5%では9座位全体を使った方が検出力は一応高くなっている。ここでは、相が分かっているデータで調べているが、実際には各個人のディプロタイプ形を推定すること考えると、D'を使って求めたLD blockに合わせた長さのハプロタイプを用いることも有効であると考えられる。
【００６８】
１−４. まとめ
本章ではハプロタイプやディプロタイプを用いた関連解析の有効性について説明した。ハプロタイプやディプロタイプを用いると連鎖不平衡の情報を自然に取り入れることができ、単座位の解析よりも有利であることが期待されるが、検出力の説明はそれ程単純ではなく単座位の方が有利である状況も多い。ハプロタイプやディプロタイプを用いると原因変異をより良く記述できるマーカーが出現する可能性が高まるが、逆に独立性検定においては自由度が大きくなるために検出力が落ちることがある。検出力の大きさはこのトレードオフによって支配されている。ここで単座位よりハプロタイプやディプロタイプを用いた解析が有利である状況をまとめると、
1.稀なハプロタイプに表現型の原因変異がのっている。
2.原因変異は頻度の高いSNPsで記述されるハプロタイプ上にのっており、ハプロタイプを用いることによってこの頻度が著しく原因変異に近くなる。
3.一つの遺伝子上に同じ表現型に影響する複数の変異がある場合。
の様な状況が考えられる。
【００６９】
２．ディプロタイプ形推定による分割表のp値の区間推定
各個人のディプロタイプ形を推定し、これを用いて分割表による関連解析を行う場合、個人のサンプリングによる揺らぎの他に、推定による揺らぎが存在する。本章では、この推定による揺らぎの検討を行う。
【００７０】
２−１.サンプリング揺らぎと推定揺らぎ
統計学的にp値とは、帰無仮説の元でそのデータ(outcome)が出現する確率である。このp値がある有意水準αを下回った場合、「サンプリングによる揺らぎを考慮しても帰無仮説が成立しているとは考え難い」と解釈し帰無仮説を棄却する。ブートストラップ法は、得られているデータを母集団として復元抽出を行うノンパラメトリックブートストラップ法と、あるモデルのパラメータ推定を得られているデータから行いこれを母数としてランダムサンプリングを行うパラメトリックブートストラップ法に大別される。得られているデータを母集団とするサンプリング揺らぎから、母集団からのサンプリングの揺らぎが評価出来るという仮定の元で、母集団のサンプリングの揺らぎを評価する方法論である。
【００７１】
疾患感受性ゲノム領域を特定する方法論として、一般集団を用いたSNPsによる関連解析がサンプルを集めやすく且つ検出力の高い方法論として注目されている。通常、この解析方法では疾患群と健常群の遺伝子型データが取得される。このとき、SNPsの相に関する情報、即ちハプロタイプの情報は失われている。よってハプロタイプやディプロタイプ形を使った関連解析を、単点のアレルや遺伝子型解析と同様に「分割表の独立性検定」を使って行うにはまず各個人のディプロタイプ形を推定することが必要になる。しかしながら、ディプロタイプ形もあるモデルの元での推定値であり、推定の揺らぎが生じる。そこでハプロタイプやディプロタイプ形を使った「分割表の独立性検定」では、母集団からのサンプリングによる揺らぎとディプロタイプ形推定の揺らぎの2つを考慮することが必要となる。
【００７２】
２−２.p値の区間推定
サンプリング揺らぎと推定揺らぎを評価するために、ブートストラップ法によるp値の分布で真のデータ（outcome）のp値を評価することを考える。ディプロタイプ形推定において、ある与えられたハプロタイプ頻度の集合
【００７３】
【数１６】
の元でベイズの定理を用いて
【００７４】
【数１７】
と事後確率によって評価することを考える。ここで
【００７５】
【数１８】
は2つのハプロタイプの組み合わせで表現したディプロタイプ形であり、ハーディ・ワインベルグ平衡の元では
【００７６】
【数１９】
ここで、δqrはクロネッカーδであり、q=rの時のみ1となり他は0である。Viは個人iの遺伝子型データの集合を表し、このデータにおいて可能な全てのハプロタイプの組み合わせvpについて和
【００７７】
【数２０】
がこのデータが出現する確率(尤度)を表すことになる。式(10)の事後確率に従って各個人のディプロタイプ形をサンプリングすれば、与えられたハプロタイプ頻度の元で推定の不確実性を考慮することができる。しかしながら、通常はハプロタイプ頻度も推定量である。そこでこの推定量の揺らぎをブートストラップ法で評価する。ブートストラップ法を用いることにより、得られたデータを母集団としてハプロタイプ頻度の区間推定を行うことができる。この信頼区間はサンプルデータを母集団としているため、推定によるハプロタイプ頻度の揺らぎを表現していると考えられる。そこで、ブートストラップ法に従い、得られているデータからブートストラップサンプルを作成するステップ、それぞれのブートストラップサンプルを用いてハプロタイプ頻度を最尤推定するステップ、得られたハプロタイプ頻度から各個人のディプロタイプ形事後確率を求めるステップ、この事後確率に従ってディプロタイプ形をサンプリングするステップ、ディプロタイプ形からハプロタイプ及び/又はディプロタイプ形による分割表を作成してp値を計算するステップによって、推定揺らぎを考慮したp値の分布を得ることができる。
【００７８】
本来、ブートストラップ法は、サンプリングによる揺らぎを評価する方法論として発達してきたが、本発明においてはブートストラップサンプルの揺らぎから母集団のサンプリングによる揺らぎを推定するのではなく、直接得られているデータによるハプロタイプ頻度推定揺らぎを評価している。このようにしてデータ（outcome）のp値の信頼区間を構成する。
【００７９】
２−３.注意点
実際の計算上における注意点を下記にまとめた。
先ず、p値の信頼区間の構成においては、ブートストラップサンプルを作成し、ハプロタイプ頻度を求め、これを使ってディプロタイプ形の事後確率を求めている。この時、稀なハプロタイプを持っている個人の事後確率が0になることがある。ブートストラップ法によるサンプリングでこの稀なハプロタイプが含まれない場合には、頻度を0と推定するためこの個人出現確率が0になる。よってハプロタイプやディプロタイプ形の分割表を構成したときに合計が元データと同じ数になるとは限らなくなる。
【００８０】
次に、ブートストラップ法によるサンプリングとディプロタイプ形事後確率に従った2つのサンプリングを行っている。このとき、ディプロタイプ形事後確率に従ったサンプリング数は連鎖不平衡の強さに従って変更することが考えられる。連鎖不平衡が強い領域では殆んどのサンプルのディプロタイプ形事後確率が1つのハプロタイプペアのみに集中し、他は殆んど0になる。この場合にはディプロタイプ形事後確率に従ったサンプリングは1又は2とすることが好ましい。しかしながら、連鎖平衡が弱い領域ではディプロタイプ形事後確率が複数のハプロタイプペアへ割れることが予想され、この場合にはより多い数のサンプリング数とすることが考えられる。
【００８１】
なお、ハプロタイプ頻度の推定揺らぎの評価方法としてブートストラップ法を用いたが、本発明はこれに限定されない。すなわち、ブートストラップ法以外のサンプリング法によってサンプリングを行っても良い。例えば、尤度分布からハプロタイプ頻度をサンプリングする方法が挙げられる。この場合、Markov chain Monte Carlo(MCMC)法などを使ってハプロタイプ頻度に依存する尤度分布を作成することができる。
【００８２】
２−４.発明の構成
上述したp値の区間推定は、図１６に示すようなコンピュータに上述した演算処理を実行させることによって実行することができる。コンピュータは、CPU101（制御手段）、ROM102、RAM103、入力手段104、送信/受信手段105、出力手段106、ハードディスクドライブ(HDD)107及びCD-ROMドライブ108を備える。また、コンピュータは、ROM102、RAM103、HDD107及びパブリックデータベース110等に記録されたデータを検索する検索手段111、検索手段111で検索したデータや入力手段104で入力したデータを加算減算処理する演算手段112を備える。
【００８３】
ここで、本発明に係る遺伝子関連解析プログラムは、例えば、ROM102、RAM103又はHDD107に記憶される。そして、本プログラムに従ってCPU101がコンピュータの上記ハードウェアの駆動制御を行うことで、上記「２−２.p値の区間推定」に説明した情報処理を行いp値の区間推定を実行する。
【００８４】
CPU101は、コンピュータ全体を制御し、遺伝子関連解析処理を実行する。RAM103は、遺伝子関連解析処理を実行する上で必要なデータを一時的に格納する。入力手段104は、キーボードやマウス等であり、遺伝子関連解析処理を実行する上で必要な条件を入力するとき等に操作される。送信／受信手段105は、CPU101の命令に基づいて、通信回線を介してパブリックデータベース110等との間でデータの送受信処理を実行する。出力手段106は、表現型と関連付けられた遺伝子形データ、入力手段104から入力された各種条件、遺伝子関連解析結果等を、CPU101からの命令に基づいて表示処理を実行する。なお、出力手段106としては、コンピュータのディスプレイ、又はプリンターなどが例示される。HDD107は、遺伝子関連解析プログラム、表現型データと関連付けられた遺伝子型データ等を格納し、CPU101の命令に基づいて、格納しているプログラム又はデータ等を読み出し、例えばRAM103に格納する。CD-ROMドライブ108は、CPU101の指示に基づいて、CD-ROM109に格納されている遺伝子関連解析プログラムや、表現型データと関連付けられた遺伝子型データ等から、プログラム又はデータ等を読み出し、例えばRAM103に格納する。なお、表現型データと関連付けられた遺伝子型データとは、個体から採取された２以上の座位に関する遺伝子型データであって、当該個体の表現型データ（例えば、疾患の有無、薬剤抵抗性の有無、副作用の有無等の情報）と関連付けられたものである。表現型データと関連付けられた遺伝子型データは、コンピュータ内部の記憶装置（HDD107、RAM103及びROM102）に格納されていてもよいが、パブリックデータベース110に格納されていてもよい。
【００８５】
本発明に係る遺伝子関連解析プログラムをコンピュータにインストールすることによって、当該コンピュータは遺伝子関連解析を実行する遺伝子関連解析装置として機能することとなる。すなわち、本発明に係る遺伝子関連解析装置は、入力手段104により２以上の座位に関する遺伝子型データを含む個体に関する遺伝子型データセット、並びに、当該個体に関する表現型データを入力する。また、遺伝子関連解析装置において演算手段112は、サンプリング法（例えばブートストラップ法）により上記遺伝子型データセットから複数のサンプリングデータセットを作成し、作成されたサンプリングデータセットに含まれる遺伝子型データセットを用いて当該サンプリングデータセットにおけるハプロタイプ頻度を最尤推定し、得られたハプロタイプ頻度から各個体のディプロタイプ形事後確率を計算し、計算された事後確率からディプロタイプ形をサンプリングし、サンプリングされたディプロタイプ形からハプロタイプ及び/又はディプロタイプ形の分割表を作成し独立性検定を行うことでp値を求め、算出されたp値からp値の分布を算出する処理を行う。そして、遺伝子関連解析装置は、演算手段112で算出されたp値の分布に基づいて、p値の信頼区間及びp値の平均値等の統計量を出力手段106により出力する。
【００８６】
本発明に係る遺伝子関連解析装置は、p値の分布に関する統計量として、例えばp値の信頼区間を出力することが好ましい。また、p値の分布に関する統計量としては、p値の信頼区間の上限値、下限値、p値の平均値、p値の中央値、幾何平均値（logPの平均値のexp）、最頻値、トリム平均を挙げることができる。また、本発明に係る遺伝子関連解析装置において演算手段112は、演算手段112で算出されたp値の信頼区間及び/又はp値の平均値を有意水準αと比較し、当該有意水準αを上回った数をカウントする処理を行い、有意水準を上回った数を出力手段106で出力しても良い。
【００８７】
本発明に係る遺伝子関連解析装置によれば、遺伝子関連解析の結果として、表現型データ（caseとcontrol）に対応するようにp値の信頼区間及びp値の平均値といったp値分布の統計量が出力される。出力されたp値の信頼区間及びp値の平均値等の統計量に基づいて、表現型データを記述するハプロタイプを特定することができる。本発明に係る遺伝子関連解析装置によれば、サンプリングの揺らぎとハプロタイプの推定の揺らぎを考慮した出力結果を得られるため、従来の大標本理論に基づく解析結果よりも信頼性の高い解析を行うことができる。
【００８８】
なお、本発明に係る遺伝子関連解析装置は、上述したp値分布に関する統計量の他に、従前の解析方法に準じて演算した解析結果をあわせて出力するものであっても良い。すなわち、演算手段112は、入力手段104によって入力された上記遺伝子型データセットから、多座位を含むハプロタイプ頻度を最尤推定し、推定されたハプロタイプ頻度から上記遺伝子型データセットに含まれる各個体のディプロタイプ形事後確率を計算し、計算された事後確率からディプロタイプ形をサンプリングし、サンプリングされたディプロタイプ形からハプロタイプ及び/又はディプロタイプ形の分割表を作成し独立性検定を行うことでp値を求める処理を行う。そして、演算手段112によって算出されたp値を出力手段106で出力する。得られたp値は、サンプリングの揺らぎとハプロタイプの推定の揺らぎを考慮した出力結果ではなく、上述したように算出されたp値分布の統計量と比較することができる。
【００８９】
３.シミュレーションデータの解析
本節では、2節において相の情報がある時にハプロタイプやディプロタイプ形の検出力を評価したデータについて、2章の方法論を使ってp値の区間推定を行い、推定の揺らぎを考慮した場合のp値の推定精度、関連解析への影響、及び検出力について調べた結果をまとめた。
【００９０】
３−１. 計算の手続き
今回用いたブートストラップによるp値の区間推定方法は、以下の計算手順に従って実行した。
1.6座位分のハプロタイプ頻度を全サンプルを用いて最尤推定
2.推定されたハプロタイプ頻度を用いて各個人のディプロタイプ形事後確率を計算
3.事後確率から個人のディプロタイプ形をサンプリング
4.サンプリングされたディプロタイプ形からハプロタイプ及びディプロタイプ形の分割表を作成して独立性検定(p値を求める)
5.ノンパラメトリックブートストラップ法に従いサンプル集団から再サンプリングを行い、各再サンプリングデータのハプロタイプ頻度を最尤推定
6.上記5.で得られた各ハプロタイプ頻度に対して、上記2.から4.を繰り返しp値の分布を得る。
7.推定するハプロタイプ頻度を1座位分ずらし、上記を繰り返す。
となる。次節以降に示す計算では、5.のブートストラップ法におけるサンプリングを1000回とし、3.の事後確率からのディプロタイプ形サンプリングはそれぞれ2回ずつ行い、1つのcase-controlデータに対して2000個のハプロタイプ及びディプロタイプ形の分割表を作成してp値を求めた。
【００９１】
３−２. ｐ値の推定精度
表8及び表9に遺伝子型データからディプロタイプ形を推定して得られたハプロタイプ及びディプロタイプ形の分割表のp値とデータ（outcome）のp値との差(絶対値)の平均と分散を、それぞれ100人対100人データ（表8）、10人対100人データ（表9）について示した。ここでMLEとは最尤推定法による単点推定（従前の解析方法に準じた方法）であり、BS meanとは本発明に係る遺伝子関連解析により得られたp値分布における統計量（平均値）であり、intervalとは本発明に係る遺伝子関連解析により得られたp値分布における信頼区間の中に真のデータ（outcome）のp値が含まれていた数である。またCaseとはcase対control-1、Controlはcontrol-1対control-2の比較結果である。
【００９２】
【表８】
【００９３】
【表９】
【００９４】
表8及び表9から判るように、どのケースでも欠測率の増加に伴い平均誤差(mean)及び推定の揺らぎ(se)が増大している。またBS meanの方がMLEよりも推定の平均誤差及び推定の揺らぎとも小さくなっている。これによりハプロタイプ頻度推定の揺らぎを考慮したBS meanの方がMLEのみによる推定よりも精度が高く安定した推定ができていると考えることができる。
【００９５】
信頼区間の構成については、表8を見る限り非常に成功していると言えるだろう。10人対100人の比較である表9では若干カバー率が悪くなっているが、サンプル数が少ない場合には推定の揺らぎが大きくなるため1000回程度のブートストラップ法では高い精度で信頼区間を構成することは難しかったためと考えられる。よって、サンプル数が少ない場合には、ブートストラップ法によるサンプルの回数を適宜増加させることが好ましい。
【００９６】
一方、推定精度のp値依存性を見るために、横軸に真のハプロタイプの分割表のp値をとり、縦軸に推定されたp値(MLE又はBS mean)から真のp値を引いた値をとってプロットした結果を図17-1及び図17-2に示した。なお、図17-1には100人対100人の場合の結果を示し、図17-2には10人対100人の場合の結果を示した。図17-1において(A)及び（B）は、欠測率を0％とした時のcase及びcontrolをそれぞれ示している。図17-1において(C)及び（D）は、欠測率を5％とした時のcase及びcontrolをそれぞれ示している。図17-1において(E)及び（F）は、欠測率を10％とした時のcase及びcontrolをそれぞれ示している。図17-2において(G)及び（H）は、欠測率を0％とした時のcase及びcontrolをそれぞれ示している。図17-2において(I)及び（J）は、欠測率を5％とした時のcase及びcontrolをそれぞれ示している。図17-2において(K)及び（L）は、欠測率を10％とした時のcase及びcontrolをそれぞれ示している。
【００９７】
図17-1及び図17-2から、欠測率が高くなると推定精度が悪くなっていることが分かる。またp値は、
【００９８】
【数２１】
であるためp値が0に近いところでは上方に、p値が1.0に近いところでは下方に分布が広がっている。このため、ある有意水準を定めた場合の検定では、やや保守的になることが予想される。この問題はブートストラップ法で得られたp値の中央値や幾何平均を使うことによって補正される可能性がある。
【００９９】
３−３. 関連解析への影響
３−３−１. 有意となった数
表10及び表11にディプロタイプ形推定によるハプロタイプの分割表においてFisher's exact testのp値が有意となった数をそれぞれ100人対100人データ（表10）及び10人対100人データ（表11）に対して示した。前節と同様にcaseとはcase対control-1、controlはcontrol-1対control-2の比較であり、knownは相が分かっているときのp値を用いた場合、MLEは最尤推定されたハプロタイプ頻度のみを用いた推定、BS meanは本発明に係る遺伝子関連解析により得られたp値分布における統計量（平均値）であり、BS upperは同じく本発明に係る遺伝子関連解析により得られたp値分布における統計量（信頼区間の上限値）であり、missingは遺伝子型データの欠測率である。
【０１００】
【表１０】
【０１０１】
【表１１】
【０１０２】
表10及び表11から、全体的にknownに比べて有意と判定されている数は減っており、また欠測率の増加に伴い有意となる数が減っていることが分かる。これは前節の図17-1及び図17-2で示した様にp値が小さいところでは大きめの予測を行うバイアスのためと考えられる。また、BS upperは欠測率の影響を強く受けており、これはp値の区間推定における信頼区間の幅が広くなったためであろう。
【０１０３】
３−３−２. 判定の精度
有意水準毎にハプロタイプの相が既知の場合のp値(正解データ)と推定値された場合のp値に対して、正解及び不正解の比較を行った。表12及び表13にハプロタイプを使って分割表を作成した場合の結果をそれぞれ100人対100人データ（表12）、10人対100人データ（表13）についてまとめた。表12及び表13においてMLE、BS mean、BS upper（90％信頼区間)については前節と同様である。また、表12及び表13において、TPは真陽性、TNは真陰性、FPは偽陽性及びFNは偽陰性の数を示している。さらに、表12及び表13において、Accとして(TP+TN)/totalの式に従って算出された値も示した。この値は判定の精度を意味している。
【０１０４】
【表１２】
【０１０５】
【表１３】
【０１０６】
表12及び表13から、単点推定であるMLEではFPの数が本発明に係る方法と比べて多くなっており、control対controlの比較ではFPとFNは大体同数程度出現している。BS upperは保守的な指標であるためFNの数が非常に多く、欠測の影響も他の指標に比べて強く出ていることも確認できる。今回の結果からはBS meanが最も精度の高い指標となっており、10人対100人のcase対controlの結果を除いて最も精度が高い。またcontrol対controlの比較においてFPの数が少なく、やや保守的な結果をもたらしているが偽陽性をコントロールしやすい指標であると言える。これらの結果は、本発明に係る遺伝子関連解析装置によってハプロタイプ頻度の推定の揺らぎを考慮してp値分布の統計量を指標にすることが、遺伝子関連解析においてより有効であることを示していると考えられる。
【０１０７】
３−４. 検出力
単座位のアリル及び遺伝子型の分割表で有意となった割合、相が分かっている状態で構築したハプロタイプの分割表で有意となった割合、及び遺伝子型データからハプロタイプの分割表のp値を区間推定した場合に有意となった割合を、100人対100人、10人対100人のデータ毎に表14及び表15にまとめた。なお、有意水準は5%とした。
【０１０８】
【表１４】
【０１０９】
【表１５】
【０１１０】
表14及び表15から、単純に有意となった割合は、MLE、BS mean、BS upperの順に多く、MLEが最も検出力が高くなっているように見える。しかしながら、MLEでは偽陽性もある程度の割合で観察されておりデータの欠測率や連鎖不平衡の大きさによっては正確な検定を行えているとは言えない これに対して、表12におけるcase対controlの比較では、最も保守的なBS upperを検定統計量として採用すると欠測率が高いと単座位よりも検出力が低くなっており、推定の揺らぎによって検出力が下がっていると考えられる。しかしながらBS meanでは相の推定における揺らぎを考慮しても、単座位よりも検出力が高くなっていることが分かる。
【０１１１】
一方で、表12におけるcontrol対controlの比較では、相の推定における揺らぎを考慮することによるバイアスが存在し、やや保守的な結果となっている。逆にサンプル数の少ない表13のケースでは相の情報があるときよりも有意と判定されている割合が多くなっている。しかしながら、表13におけるcontrol対controlの相の推定による偽陽性の増大は僅かであり、今回の方法論ではBS meanを使うとやや保守的になるが、偽陽性の量をある程度コントロールすることが可能であると考えられる。また、単純な相加平均よりも頑健な中央値や幾何平均(logpの平均)を検定統計量として使用すれば、このバイアスが小さいなることが期待できる。
【０１１２】
３−５. 大標本理論に基づく検定統計量の影響
ハプロタイプやディプロタイプを用いた分割表の独立性検定では、単座位のSNPsによる分割表とは異なりカウント数の少ない升目が出現する。そのため大標本理論に基づくχ2検定は適さず、Fisher's exact testの様な小標本検定が必要になる。この様子を表16及び表17に、同じハプロタイプ頻度母数からランダムにそれぞれ100人対100人（表16）、10人対100人ずつ（表17）100ケースサンプリングして得られた相が既知の場合におけるcontrol対controlの比較の結果有意となった割合を示した。
【０１１３】
【表１６】
【０１１４】
【表１７】
【０１１５】
表16ではFisherに比べてχ2検定は保守的になっており、逆に表17では与えられた有意水準以上に有意となる割合が多くなっている。Fisher's exact testの結果は両方ともほぼ有意水準に近い値を示しており、帰無仮説の下で偽陽性をコントロールできていることが分かる。
【０１１６】
５．まとめ
ハプロタイプやディプロタイプを使った関連解析の有効性について主に検出力について説明した。大標本理論の元では近似的な検出力曲線を描くことが可能であり、この結果、原因変異とr2=1となるようなマーカーを用いることが非常に有効であることがわかった。対立因子が2であるSNPsに比べて、ハプロタイプでは複数の対立因子が出現するため、原因変異とr2が高いマーカーが出現する可能性が高まる。特にcommon SNPsマーカーは通常頻度5%以上であるため、原因変異が稀な場合にはハプロタイプを用いた解析が非常に有効であった。このような状況は、表現型データとして薬剤に対する副作用の有無などが考えられ、ハプロタイプやディプロタイプ形を使った独立性検定は、pharmacogenomics(PGx)の解析において特に有効であると考えられる。
【０１１７】
しかしながら、通常のSNPsデータには相の情報がなく直接ハプロタイプを使って関連解析を行うことはできない。そこで遺伝子型データからハプロタイプ頻度やディプロタイプ形を推定することになるが、推定量であるために推定の揺らぎが生じる。よって、実際のハプロタイプやディプロタイプ形を使った関連解析では、この推定の揺らぎを考慮しなければならない。推定の揺らぎを評価するために、一例としてブートストラップ法によるp値の分布を算出してp値の区間推定を行った。シミュレーションデータデータを使った評価では、p値の推定精度は欠損率に依存するが、統計的な有意差を検討することは十分可能であり、また指標によって偽陽性の量をコントロールすることも可能である。よって実際の遺伝子型データから単座位解析では検出することが難しい領域が、本発明に係る遺伝子関連解析装置によって検出することが可能となる。この場合、帰無仮説を棄却する根拠として、「サンプリングの揺らぎと推定の揺らぎを同時に考慮しても、帰無仮説が成立しているとは考えにくい。」とすることができる。
【０１１８】
〔参考例〕
Ａ． 連鎖不平衡指標（Linkage disequilibrium indices）
座位1における所定のアリルをa1、座位2における所定のアリルをa2とし、それぞれのアリル頻度をpa1及びpa2とする。また座位1におけるa1以外のアリルを
【０１１９】
【数２２】
座位2におけるa2以外のアリルを
【０１２０】
【数２３】
と表現する。更に、この2つのアリルのハプロタイプa1a2の頻度をha1a2とすると連鎖不平衡指標D'及びr2は、
【０１２１】
【数２４】
と定義される。
【０１２２】
Ｂ. 大標本理論に基づく検出力の近似式
簡単のため2×2の分割表におけるχ2検定について考える。今、2つの二項分布Bi（n1, p1）とBi（n2, p2）に従う確率変数をそれぞれX1及びX2とすると、
【０１２３】
【数２５】
の時X1及びX2はそれぞれN(n1p1, n1p1(1-p1))及びN(n2p2, n2p2(1-p2))に従う。従って
【０１２４】
【数２６】
は、N(p1, p1(1-p1)/n1)及びN(p2, p2(1-p2)/n2)に従い、
【０１２５】
【数２７】
となる。この時、母比率の検定を帰無仮説H0及び対立仮説H1をそれぞれ
【０１２６】
【数２８】
として検定を行うことを考える。帰無仮説H0の元では、検定統計量
【０１２７】
【数２９】
はN(0,1)に従う。ここで
【０１２８】
【数３０】
である。従ってZ2は自由度1のχ2分布に従い、実際Z2は2×2の分割表のχ2統計量と一致する。
【０１２９】
有意水準αにおける両側検定は式(18)を使うと
【数３１】
と表されるので、
【０１３０】
【数３２】
となる。ここで
【０１３１】
【数３３】
は 式(16)からN(0,1)に従うので、
【０１３２】
【数３４】
と近似する(即ち
【０１３３】
【数３５】
の揺らぎを無視)と、検出力β(n1,p1,n2,p2)(図18参照)は
【０１３４】
【数３６】
と近似される 。ただし、
【０１３５】
【数３７】
であり
【０１３６】
【数３８】
である。
【図面の簡単な説明】
【０１３７】
【図１】D'=1の場合における原因変異とマーカーとの関係を示す模式図である。
【図２】オッズ比をOR=1.5に固定した場合のタイプ1、タイプ2の場合における検出力を示す特性図である。
【図３】オッズ比をOR=3.0に固定した場合のタイプ1、タイプ2の場合における検出力を示す特性図である。
【図４】原因変異のオッズ比を3.0とした場合のマーカーMの見掛けのオッズ比を示す特性図である。
【図５】原因変異のオッズ比を3.0とした場合のマーカーMのマーカーの頻度差f'M-fを示す特性図である。
【図６】シミュレーションデータとして使用したSNPsデータにおける1〜9座位のLDマップ（上三角部分がr2、下三角部分がD'）を示す特性図である。
【図７】原因変異の割合の増分を与えたときのハプロタイプ頻度及びアレル頻度を示す特性図である。
【図８−１】Control-1対control-2で100ケースずつサンプリングした時の単座位のアレル、及び3〜9座位のハプロタイプによる分割表の独立性検定(Fisher's exact test)のp値の分布を示すボックス図である。
【図８−２】Control-1対control-2で100ケースずつサンプリングした時の単座位のアレル、及び3〜9座位のハプロタイプによる分割表の独立性検定(Fisher's exact test)のp値の分布を示すボックス図である。
【図９−１】Control-1対control-2で100ケースずつサンプリングした時の単座位のアレル、及び3〜9座位のハプロタイプによる分割表の独立性検定(Fisher's exact test)で得られたp値に対して有意水準を1%及び5%とした時に有意となった数を示すヒストグラムである。
【図９−２】Control-1対control-2で100ケースずつサンプリングした時の単座位のアレル、及び3〜9座位のハプロタイプによる分割表の独立性検定(Fisher's exact test)で得られたp値に対して有意水準を1%及び5%とした時に有意となった数を示すヒストグラムである。
【図１０−１】Case-1対control-1で100ケースずつサンプリングした時の単座位のアレル、及び3〜9座位のハプロタイプによる分割表の独立性検定(Fisher's exact test)のp値の分布を示すボックス図である。
【図１０−２】Case-1対control-1で100ケースずつサンプリングした時の単座位のアレル、及び3〜9座位のハプロタイプによる分割表の独立性検定(Fisher's exact test)のp値の分布を示すボックス図である。
【図１１−１】Case-1対control-1で100ケースずつサンプリングした時の単座位のアレル、及び3〜9座位のハプロタイプによる分割表の独立性検定(Fisher's exact test)で得られたp値に対して有意水準を1%及び5%とした時に有意となった数を示すヒストグラムである。
【図１１−２】Case-1対control-1で100ケースずつサンプリングした時の単座位のアレル、及び3〜9座位のハプロタイプによる分割表の独立性検定(Fisher's exact test)で得られたp値に対して有意水準を1%及び5%とした時に有意となった数を示すヒストグラムである。
【図１２−１】Case-2対control-1で100ケースずつサンプリングした時の単座位のアレル、及び3〜9座位のハプロタイプによる分割表の独立性検定(Fisher's exact test)のp値の分布を示すボックス図である。
【図１２−２】Case-2対control-1で100ケースずつサンプリングした時の単座位のアレル、及び3〜9座位のハプロタイプによる分割表の独立性検定(Fisher's exact test)のp値の分布を示すボックス図である。
【図１３−１】Case-2対control-1で100ケースずつサンプリングした時の単座位のアレル、及び3〜9座位のハプロタイプによる分割表の独立性検定(Fisher's exact test)で得られたp値に対して有意水準を1%及び5%とした時に有意となった数を示すヒストグラムである。
【図１３−２】Case-2対control-1で100ケースずつサンプリングした時の単座位のアレル、及び3〜9座位のハプロタイプによる分割表の独立性検定(Fisher's exact test)で得られたp値に対して有意水準を1%及び5%とした時に有意となった数を示すヒストグラムである。
【図１４−１】Case-3対control-1で10人対100人サンプリングした時の単座位のアレル、及び3〜9座位のハプロタイプによる分割表の独立性検定(Fisher's exact test)のp値の分布を示すボックス図である。
【図１４−２】Case-3対control-1で10人対100人サンプリングした時の単座位のアレル、及び3〜9座位のハプロタイプによる分割表の独立性検定(Fisher's exact test)のp値の分布を示すボックス図である。
【図１５−１】Case-3対control-1で10人対100人サンプリングした時の単座位のアレル、及び3〜9座位のハプロタイプによる分割表の独立性検定(Fisher's exact test)で得られたp値に対して有意水準を1%及び5%とした時に有意となった数を示すヒストグラムである。
【図１５−２】Case-3対control-1で10人対100人サンプリングした時の単座位のアレル、及び3〜9座位のハプロタイプによる分割表の独立性検定(Fisher's exact test)で得られたp値に対して有意水準を1%及び5%とした時に有意となった数を示すヒストグラムである。
【図１６】本発明を適用した遺伝子関連解析装置のハードウェア構成を示す模式図である。
【図１７−１】横軸に真のハプロタイプの分割表のp値をとり、縦軸に推定されたp値(MLE又はBS mean)から真のp値を引いた値をとってプロットした結果を示す特性図である。
【図１７−２】横軸に真のハプロタイプの分割表のp値をとり、縦軸に推定されたp値(MLE又はBS mean)から真のp値を引いた値をとってプロットした結果を示す特性図である。
【図１８】帰無仮説H0及び対立仮説H1における統計的検定の模式図である。
【符号の説明】
【０１３８】
101…CPU（制御手段）、102…ROM、103…RAM、104…入力手段、105…送信/受信手段、106…出力手段、107…ハードディスクドライブ(HDD)、108…CD-ROMドライブ、109…CD-ROM、110…パブリックデータベース、111…検索手段、112…演算手段
【技術分野】
【０００１】
本発明は、集団から取得した遺伝子型データ及び表現型データを用いて関連解析を行う遺伝子関連解析装置及び遺伝子関連解析プログラムに関する。
【背景技術】
【０００２】
表現型を記述するゲノム領域を特定するためにSNPsを用いた関連解析が数多く行われている。現在は個人当り55万SNPsを一度に測定する技術なども開発されており、大量の多型情報を用いた高密度SNPs解析が主流となっている。高密度SNPs解析ではSNPs間の遺伝的距離が短く、連鎖不平衡があるためにマーカーであるSNPsの情報が独立ではないという特徴がある。そのため実際の関連解析においてはこのような連鎖不平衡を考慮した解析が求められる。
【０００３】
一方でSNPは通常2つの対立因子(allele)しか持たない。そのため、ある座位に関するSNPを用いた2つの表現型間の独立性検定では、alleleやdominant model、recessive modelの元での2×2の分割表や、遺伝子型を使った2×3の分割表などが通常用いられる。表現型への寄与がある程度あり、比較的頻度の高いSNPsの検出には、このような解析は非常に有効であるが、
1.様々な因子の効果により1つの座位の影響が小さい場合
2.頻度の小さな原因多型
の場合にはこれらを検出することは困難であると考えられる。
【０００４】
1つの座位の影響が小さい場合には、環境要因や複数の遺伝子間の相互作用等が考えられる。これらの解析には環境要因を記述する変数の導入や、 複数のマーカーとの関連などを考慮した解析が必要になり、検出は非常に困難になることが予想される。
【０００５】
マーカーに用いるSNPsセットに常に目的の表現型の原因となるSNPsが含まれていると仮定することは出来ない。元々マーカーとしてのSNPsは、マーカーとしての有効性のために頻度の小さいSNPsは選ばれていない。そのため表現型の寄与度が大きくとも、頻度が低ければこれがSNPマーカーとなることはない。一方で頻度の小さいSNPは、たとえ周辺領域との連鎖不平衡が高くとも、頻度の高い周辺の1つのSNPで記述することは出来ない。このような状況の例としては、common disease/multiple rare variants(CD/MRVs)や、非常に稀な副作用を起こす薬剤に対してこの原因を記述するゲノム上の変異を検出する場合等が考えられる。このような非常に稀な多型(例えば頻度0.1%)を効率的に検出するためには、 頻度の低いマーカーが必要になる。頻度の低いマーカーとしては、稀なSNPsを 網羅的にタイピングすることや、マイクロサテライトマーカーを用いることが考えられるが、タイピングコストを考えると決して効率的な方法であるとは言えない。
【０００６】
そこで頻度の高いSNPsを複数個使用することで対立因子数を増やすことが考えられる。目的となる感受性変異がマーカーの変異よりも後で起った場合を想定すると、周辺のSNPや複数のSNPsを組み合わせて検出するよりもこの原因多型を挟んだ複数のマーカーによるハプロタイプやディプロタイプを用いた解析の方がより有意差を検出しやすく、有効な解析であると考えられる。
【０００７】
これまで大規模な家系を用いた遺伝統計学的手法では、家系情報からハプロタイプを特定し、これをマーカーとして独立性検定を行ってきた。しかしながら一般集団を用いた関連解析では通常各個人のハプロタイプを同定することは出来ない。そのため集団の情報を使って推定する方法論が盛んに研究されてきた。このような方法論の延長としてハプロタイプを用いた関連解析も研究されており、各表現型毎に最尤法によって求められたハプロタイプ頻度を使った尤度比検定などの検出方法論が研究されている。
【発明の開示】
【発明が解決しようとする課題】
【０００８】
しかしながら、尤度比検定は大標本理論に基づいており、稀な頻度差を検討するには不向きな方法論である。そこで、本発明は、稀なマーカーの頻度差についても検討可能な小標本理論に基づいた検出方法を提供することを目的とする。
【課題を解決するための手段】
【０００９】
上述した目的を達成するため、本発明者が鋭意検討した結果、小標本に適用可能な検定方法としては、分割表を用いたFisher's exact testを挙げられるが、ハプロタイプやディプロタイプは通常の遺伝子型データから直接求めることが出来ないため、直接Fisher's exact test を行うことは出来ないといった問題を見いだした。そして、この問題を解決するため、ハプロタイプやディプロタイプを推定するとともに、サンプリングによる揺らぎ及び相の推定による揺らぎを考慮した遺伝子関連解析方法を開発することに成功した。
【００１０】
本発明は以下を包含する。
（１）２以上の座位に関する遺伝子型データを含む個体に関する遺伝子型データセット、並びに、当該個体に関する表現型データを入力する入力手段と、サンプリング法により上記遺伝子型データセットから複数のサンプリングデータセットを作成し、作成されたサンプリングデータセットに含まれる遺伝子型データセットを用いて当該サンプリングデータセットにおけるハプロタイプ頻度を最尤推定し、得られたハプロタイプ頻度から各個体のディプロタイプ形事後確率を計算し、計算された事後確率からディプロタイプ形をサンプリングし、サンプリングされたディプロタイプ形からハプロタイプ及び/又はディプロタイプ形の分割表を作成し独立性検定を行うことでp値を求め、算出されたp値からp値の分布を算出する処理を行う演算手段と、上記演算手段で算出されたp値の分布に基づいて当該分布の統計量を出力する出力手段とを備える遺伝子関連解析装置。
【００１１】
（２）上記出力手段は、上記分布の統計量として、少なくともp値の信頼区間の上限値及び平均値を出力することを特徴とする（１）記載の遺伝子関連解析装置。
【００１２】
（３）上記演算手段は、上記サンプリング法としてブートストラップ法を適用してサンプリングデータセットを作成することを特徴とする（１）記載の遺伝子関連解析装置。
【００１３】
（４）上記演算手段は、上記独立性検定として少なくともフィッシャー正確性検定を適用してｐ値を算出することを特徴とする（１）記載の遺伝子関連解析装置。
【００１４】
（５）入力手段、演算手段及び出力手段を備えるコンピュータに、上記入力手段が、２以上の座位に関する遺伝子型データを含む個体に関する遺伝子型データセット、並びに、当該個体に関する表現型データを入力するステップと、上記演算手段が、サンプリング法により上記遺伝子型データセットから複数のサンプリングデータセットを作成し、作成されたサンプリングデータセットに含まれる遺伝子型データセットを用いて当該サンプリングデータセットにおけるハプロタイプ頻度を最尤推定し、得られたハプロタイプ頻度から各個体のディプロタイプ形事後確率を計算し、計算された事後確率からディプロタイプ形をサンプリングし、サンプリングされたディプロタイプ形からハプロタイプ及び/又はディプロタイプ形の分割表を作成し独立性検定を行うことでp値を求め、算出されたp値からp値の分布を算出する処理を行うステップと、上記出力手段が上記演算手段で算出されたp値の分布に基づいて当該分布の統計量を出力するステップとを実行させる遺伝子関連解析プログラム。
【００１５】
（６）上記出力手段は、上記分布の統計量として、少なくともp値の信頼区間の上限値及び平均値を出力することを特徴とする特徴とする（５）記載の遺伝子関連解析プログラム。
【００１６】
（７）上記演算手段は、上記サンプリング法としてブートストラップ法を適用してサンプリングデータセットを作成することを特徴とする（５）記載の遺伝子関連解析プログラム。
【００１７】
（８）上記演算手段は、上記独立性検定として少なくともフィッシャー正確性検定を適用してｐ値を算出することを特徴とする（５）記載の遺伝子関連解析プログラム。
【発明の効果】
【００１８】
本発明に係る遺伝子関連解析装置及び遺伝子関連解析プログラムによれば、原因となる多型の頻度が非常に低い場合であっても、表現型との関連を検出することができる。
【発明を実施するための最良の形態】
【００１９】
以下、本発明を詳細に説明する。
１. ハプロタイプやディプロタイプを用いた関連解析の有効性
１−１. CD/CV vs CD/MRVs
晩成型のcommon diseaseの原因においてSNPsの様なgerm-lineからの寄与が小さくないと考えられる場合には、大きくcommon disease/Common variant(CD/CV)仮説と common disease/multiple rare variants(CD/MRVs)仮説が考えられる。一般集団を用いたゲノムワイドなSNPs解析は、主にCD/CV仮説に従う変異を抽出することを目的としている。しかしながら、突然変異率や集団での遺伝的浮動を考えるとCD/CV仮説に従う変異はそれ程多く無く、むしろCD/MRVs仮説に従う多数の稀な変異が多い可能性を示唆している。このような稀な変異は、これまで主に伝統的な連鎖解析によって抽出されてきた。
疾患に寄与する変異を考える場合、表現型への寄与度と集団中での頻度とをパラメータとして考えることができる。各場合の検出方法論について表1にまとめた。
【００２０】
【表１】
【００２１】
頻度が小さく、寄与度も小さい原因変異を特定するには非常に多数のサンプルが必要になると考えられ、これを検出することは困難であると考えられる。一方でこれまで連鎖解析の対象であった頻度が低く、寄与度が大きい原因変異をcommon SNPsを用いたcase-control studyで検出することができる可能性がある。一般に表現型原因変異とマーカーとのr2(後述する参考例A参照)が小さくなると検出力は劇的に下がる。そこで頻度の小さい変異を検出するためには、
【００２２】
【数１】
となるような頻度の小さいマーカーが必要になる。Common SNPsから効率的に頻度の小さいマーカーを作成する方法としてハプロタイプを利用することが考えられる。
【００２３】
１−２.関連解析の検出力
本節ではマーカー密度が非常に高い場合における検出力について説明する。まずは原因変異とマーカーが一致していないために起る見掛けのオッズ比(マーカーオッズ比)と検出力について説明し、続いて大標本理論に基づく検出力の近似式(後述する参考例B参照)を用いて様々なパラメータに対する検出力の傾向を検討する。最後にsimulation dataを使った検出力について説明し、ハプロタイプを使った関連解析がどのような場合に有効であるのかを説明する。
【００２４】
１−２−１. 見掛けのオッズ比
表現型の原因多型が稀な場合、例えオッズ比が高くとも所謂common SNPsでは、これを検出することは非常に困難になる。例えばオッズ比が5となるような変異であっても集団中の頻度が1%であればcaseでの頻度は4.8%程度にしかならない。この3.8%の頻度差を頻度10%のSNPsで検出する場合、見掛けのオッズ比は最大(D'=1)でも
【００２５】
【数２】
となってしまう。よって、稀な変異を単座位のcommon SNPで解析する場合、影響力の大きい稀なSNPが原因であったとしても見掛けのオッズ比が小さくなり、これを検出するにはより多くのサンプルが必要となる。頻度の小さいマーカーで記述したときの分割表を表2に示し、common SNPで記述したときの分割表を表3に示す。上記の状況で表現型の原因となる稀なSNPをタイピングしていた場合には、表2のような分割表が期待されP=0.019となる。しかしながら、このSNPとD'=1となる様な頻度10%のSNPでも期待値としては表3の様な分割表となり、P=0.218となってしまう。
【００２６】
【表２】
【００２７】
【表３】
【００２８】
表３に示した分割表のように、稀な変異を単座位のcommon SNPで解析する場合にはp値が大きくなり帰無仮説を棄却できない。すなわち、単座位のcommon SNPと表現型との関連性を検出することができない。
【００２９】
１−２−２. D'=1モデル
表現型の原因となる変異の割合がcaseとcontrolの間で異なるモデルを考える。 ただしマーカーの密度は十分に高く、原因変異とD'=1となる場合のみを考える。この場合、原因変異とマーカーとの関係は図1に示す2通りが考えられる。ここで表現型の原因変異をD、対立因子をdと表し、 この変異とD'=1の関係にあるマーカーのalleleをM及びmで表す。また集団中でのこれらのalleleの頻度をそれぞれfD、fd、fM及びfmと表現することにし、ハプロタイプ頻度をそれぞれhDM、hDm、hdM及びhdmとする。以下では図１に示したタイプ１とタイプ２におけるマーカーの頻度を導き、このマーカーの検出力を説明する。
【００３０】
タイプ１
一般集団中では、マーカーの存在比はallele頻度に従いfM：fmであるが、caseにおいては原因変異Dの頻度の増加に伴いマーカーの頻度にも歪みが生じる。原因変異の割合の増分をδと書くと、原因変異の比率はfD＋δ：fdとなり、マーカーの存在比はfM＋δ：fmとなる。この時、caseにおけるマーカーMの頻度f'Mは規格化を考慮すると
【００３１】
【数３】
となり、原因変異によりもたらされるマーカー頻度の増分は
【００３２】
【数４】
と表されることになる。式(3)から、規格化を考慮するとマーカーの頻度差の増分は集団中(control)のマーカー頻度に依存することがわかる。式(3)では集団中のマーカーの頻度が増大するに従って頻度差が線形に減少している。
【００３３】
また式(3)は原因変異の頻度fDに依存していない。これは様々な効果をδに押し込めているためで、実際にはδはfdとオッズ比ORから決まるパラメータである。簡単な計算から
【００３４】
【数５】
となり、式(3)は、
【００３５】
【数６】
と表現される。
【００３６】
タイプ２
原因変異がcommonの場合には、マーカーよりも原因変異が古い場合も想定される。このとき、集団中では、fD：fdの比であった原因変異の割合が、case中ではfD＋δ：fdとなったとする。このδの増分がハプロタイプ頻度に従ってマーカーの頻度に分配されたとすれば、hDMの頻度は
【００３７】
【数７】
となる。よって、よってマーカーMのcaseにおける頻度は、hDM＝fMであるから
【００３８】
【数８】
となり、頻度差は
【００３９】
【数９】
となる。式(4)を用いれば式(8)は
【００４０】
【数１０】
となる。式(5)と式(9)を比較すると、タイプ１ではfMが増加すると頻度差が減少するが、タイプ２ではfDが増加すると頻度差が減少する。ただしtype 1では
【００４１】
【数１１】
であるが、タイプ２では
【００４２】
【数１２】
である。即ち両ケースともfM=fDが頻度差の上限となる。
【００４３】
１−２−３.D'=1モデルの検出力
検出力は原因変異の頻度とオッズ比、及びマーカー頻度に依存する。オッズ比をOR=1.5及びOR=3.0に固定した場合のタイプ1、タイプ2の場合における検出力を図2及び図3に示した。ただしこの検出力は２×２の分割表におけるχ2検定の近似式を用いた場合である(後述する参考例B参照)。fD=fMの条件の場合(即ち原因変異とマーカーとの間でr2=1)に検出力が最大になり、タイプ1の場合には
【００４４】
【数１３】
タイプ2の場合には
【００４５】
【数１４】
であるためこの曲線の内側に各検出力曲線が描かれる。なお、図2及び図3において、(A)は3次元プロットを示し、(B)は2次元プロットを示している。図2及び図3においては、原因頻度を0.01から0.05刻みで表示している(caseを100人、controlを100人とした)。
【００４６】
また、原因変異のオッズ比を3.0とした場合のマーカーMの見掛けのオッズ比とマーカーの頻度差f'M-fMをそれぞれ図4及び図5に示した。これらの曲線から原因変異とマーカーとの間の連鎖不平衡がD'=1であっても、頻度が異なれば急激に検出力が落ちることがわかる。これまでの知見通り、原因変異の頻度が0.5付近でfD=fMの場合に検出力は最大となるが、
【００４７】
【数１５】
の場合には必ずしもマーカーの頻度が0.5付近で良いと言う訳ではなく、r2=1の様なマーカーが最も検出力が高い。そのため様々な頻度のマーカーを用意することが望ましいが、現実的には難しい。このためcommon variantの検出においてもハプロタイプの方が有利である状況は十分に考えられ、ハプロタイプを使った関連解析はcommon variant及びrare variantsの検出を試みる上で検討する価値はあると考えられる。なお、図4及び図5において、(A)は3次元プロットを示し、(B)は2次元プロットを示している。図4及び図5においては、原因頻度を0.01から0.05刻みで表示している(caseを100人、controlを100人とした)。
【００４８】
１−２−４. 注意点
以上の説明からD'=1のマーカーであったとしても原因変異とマーカーの頻度が異なれば検出力は急激に落ちるため、ハプロタイプを使った解析は非常に有効である。しかしながら、ここでハプロタイプを使った場合に検出力を下げる要因を挙げておく。
1.前節の検出力の説明は、2×2の分割表における大標本理論を用いた場合の結果である。ハプロタイプを用いた場合には2×nの分割表になり、このため自由度が大きくなるために検出力が落ちる。
2.遺伝子型データでは相が分かっていないために生じる推定誤差により検出力が下がる。
3.分割表のカラムを幾つか統合検定やdominantやrecessiveのモデルを 導入した場合、多数の組み合わせが発生し新たな多重性が生じる。
4.必ずしも稀な原因変異が稀なハプロタイプに乗っているとは限らない。頻度1%のハプロタイプに原因変異がのっている確率はやはり1%であろう。この場合にはハプロタイプを構成する座位数を増やすことが考えられるが、alleleの数が多くなりすぎると自由度が大きくなり検出力は落ちる。また各個人のディプロタイプ形の推定も困難になる。よってどこまでハプロタイプを伸ばせば良いのかの議論は別に必要になる。この場合、連鎖不平衡マップなどが役に立つと考えられる。
【００４９】
これまでの解析は大標本理論に基づく２×２の分割表の結果である。そこでより現実的な２×ｎの分割表となるシミュレーションデータを使って検出力を調べた結果を次節に示した。
【００５０】
１−３.シミュレーションデータを使った検出力
１−３−１. Case及びcontrolにおける母数の作成
実際のSNPsデータから9座位分のハプロタイプ頻度を推定し、これをcontrolのハプロタイプ頻度母数とした(表4参照)。
【００５１】
【表４】
【００５２】
またこの時の連鎖不平衡指標r2及びD'を表5及び図6に示した。なお、図６において上三角部分がr2であり、下三角部分がD'である。
【００５３】
【表５】
【００５４】
図６から、3座位目から6座位目までの4座位間でD'の値が非常に高くブロックを形成していると考えられる。この4座位分のハプロタイプ頻度を表6に示した。
【００５５】
【表６】
【００５６】
この4座位のハプロタイプにおいて、2番目、3番目に頻度の高いハプロタイプにそれぞれδ=0.15を与えてcase-1、case-2のハプロタイプ頻度母数を作成した。また希なハプロタイプの影響を調べるために、表4における11番目のハプロタイプにやはりδ=0.15を与えてcase-3のハプロタイプ頻度母数を作成した。この結果得られたハプロタイプ頻度及びアレル頻度も表4に示し、これを図7に図示した。なお、図７(A)はハプロタイプ頻度母数を示し、図７(B)はアレル頻度母数を示している。
【００５７】
１−３−２.シミュレーションデータの作成
表4で示した母数を使って、それぞれ100人づつ100ケースランダムにサンプリングを行った。ただしcontrolに対してはcontrol対controlの比較を行うためにcontrol-1及びcontrol-2としてそれぞれ100ケース、計200ケースのサンプリングを行った。このサンプリング方法では、相乗モデル（multiplicative model）となり、caseでもハーディ・ワインベルグ平衡（Hardy-Weinberg equilibrium）が成立していることになる。尚、各ハプロタイプ頻度、アレル頻度の平均、分散は多項分布に従うため、分散はnpi(1-pi)となる(ここでnはサンプリング数であり、piは頻度を意味する)。
【００５８】
１−３−３.シミュレーションデータの独立性検定の結果
各ケースにおける独立性検定の結果を、ボックス図と有意となった数のヒストグラムとで示した。座位毎の特異性を見るためにアレルについては各座位毎の結果を示し、またハプロタイプは3〜9座位の長さで構成し、これを、1座位ずつwindowをずらして計算した結果を示した。
【００５９】
Control-1 vs control-2 (100人対100人)
Control-1対control-2で100ケースずつサンプリングした時の単座位のアレル、及び3〜9座位のハプロタイプによる分割表の独立性検定(Fisher's exact test)のp値の分布を図8-1及び図8-2にボックス図として示した。また、有意水準を1%及び5%とした時に有意となった数のヒストグラムを図9-1及び図9-2に示した。Control対controlの比較ではp値の分布は一様分布となるが、図8-1及び図8-2は大体対称に分布しており、また平均も0.5付近にある。偽陽性の数も多少の増減はあるがほぼ期待通りの結果となっている。図8-1及び図8-2において、(A)はアレルについてのp値を示し、（B）〜(H)はそれぞれ3〜9座位のハプロタイプについてのp値を示している。図9-1及び図9-2において、(A)はアレルについて有意となった数を示し、（B）〜(H)はそれぞれ3〜9座位のハプロタイプについて有意となった数を示している。
【００６０】
Case-1 vs control-1 (15% ハプロタイプ, 100人対100人)
上記と同様に、Case-1対control-1の場合におけるp値の分布を図10-1及び図10-2にボックス図として示した。また、有意水準を1%及び5%とした時に有意となった数のヒストグラムを図11-1及び図11-2に示した。図10-1及び図10-2において、(A)はアレルについてのp値を示し、（B）〜(H)はそれぞれ3〜9座位のハプロタイプについてのp値を示している。図11-1及び図11-2において、(A)はアレルについて有意となった数を示し、（B）〜(H)はそれぞれ3〜9座位のハプロタイプについて有意となった数を示している。
【００６１】
Case-2 vs control-1 (10% ハプロタイプ, 100人対100人)
上記と同様に、Case-2対control-1の場合におけるp値の分布を図12-1及び図12-2にボックス図として示した。また、有意水準を1%及び5%とした時に有意となった数のヒストグラムを図13-1及び図13-2に示した。図12-1及び図12-2において、(A)はアレルについてのp値を示し、（B）〜(H)はそれぞれ3〜9座位のハプロタイプについてのp値を示している。図13-1及び図13-2において、(A)はアレルについて有意となった数を示し、（B）〜(H)はそれぞれ3〜9座位のハプロタイプについて有意となった数を示している。
【００６２】
Case-3 vs control-1 (1.5% ハプロタイプ, 10人対100人)
上記と同様に、Case-3対control-1の場合におけるp値の分布を図14-1及び図14-2にボックス図として示した。また、有意水準を1%及び5%とした時に有意となった数のヒストグラムを図15-1及び図15-2に示した。図14-1及び図14-2において、(A)はアレルについてのp値を示し、（B）〜(H)はそれぞれ3〜9座位のハプロタイプについてのp値を示している。図15-1及び図15-2において、(A)はアレルについて有意となった数を示し、（B）〜(H)はそれぞれ3〜9座位のハプロタイプについて有意となった数を示している。
【００６３】
１−３−４.シミュレーションデータのまとめ
シミュレーションデータを使って単座位のみを用いた独立性検定とハプロタイプによる分割表の独立性検定のp値の分布について調べた（図8、10、12及び14）。case対controlの比較において有意となった数は、連鎖不平衡の構造に強く依存しており、頻度15%のハプロタイプに変異を与えた場合はハプロタイプの方が有意になった割合が多く、頻度10%では単座位の検定において突出して有意となった数が多い座位があった。これはcase-2の場合に、5番目の座位のアレル頻度差が4座位分のハプロタイプ頻度差と同じ頻度差(約12%)となっている特殊事情のためである(即ち、この1つのアレルがこのハプロタイプのtag SNPとなっている)。
一方でcase-3では単座位の検定で有意となった割合は非常に小さく、偽陽性の量と余り変わらない。しかしながらハプロタイプを用いると有意となる数が非常に多くなることがわかる。この結果は、希なハプロタイプ変異の特徴である。
【００６４】
検出力
次に、実際の検出過程を想定して、各データ毎に有意と判定された数の内訳を調べた（図9、11、13及び15）。ここでは、有意水準に対してBonferroni補正を行い、アレル及びハプロタイプを用いた検定で共に有意差が検出されなかった場合の数、アレルの検定のみで検出された数、ハプロタイプを用いた検定のみで検出された数、アレル及びハプロタイプを用いた検定の両方で有意となった数をそれぞれ表7に示した。
【００６５】
【表７】
【００６６】
表７に示すように、Case-1の場合において、有意水準5%でハプロタイプでのみ有意となった数がアレルのみで有意となった数を若干上回っているが、有意水準1%では単座位の検定の方がやや有利となっている。Case-2の場合において、上述の特殊事情により、単座位での検定の方が検出力が高い。前節の説明では頻度の小さいハプロタイプに揺らぎを与えた場合、common SNPsに対してハプロタイプを用いた検定が有利になるはずであったがこのケースでは逆の結果となった。しかしながらハプロタイプを用いて検出される数(表中のハプロタイプとbothを足した数)はcase-1に比べて増えており、頻度の小さなハプロタイプに頻度差がある場合の方が検出力が高くなる傾向は再現している。この様に、単座位、ハプロタイプを用いた検定の有利/不利は個別の事情を反映しており単純ではない。ただし、傾向としては頻度の小さなハプロタイプに揺らぎがある場合に、common SNPsに対してハプロタイプがより有利になることが多い。実際希なハプロタイプの場合であるcase-3では、圧倒的にハプロタイプが有利となっている。
【００６７】
ハプロタイプの座位数
Case-1及びcase-2では4座位分(3〜6座位)のハプロタイプにδ=0.15の揺らぎを与えてcaseのハプロタイプ頻度母数を作成したが、検出力は必ずしも4座位のハプロタイプが高い訳ではなかった。検出力はハプロタイプの長さに対して比較的ロバストであると考えられる。一方で、Case-3では9座位全体のハプロタイプに揺らぎを与えているが、図14及び図15に示すように、前の方(1〜6番目の座位程度まで)が小さなp値が多く、有意となる数も多い。Bonferroni補正後の検出された割合でみると、有意水準5%では9座位全体を使った方が検出力は一応高くなっている。ここでは、相が分かっているデータで調べているが、実際には各個人のディプロタイプ形を推定すること考えると、D'を使って求めたLD blockに合わせた長さのハプロタイプを用いることも有効であると考えられる。
【００６８】
１−４. まとめ
本章ではハプロタイプやディプロタイプを用いた関連解析の有効性について説明した。ハプロタイプやディプロタイプを用いると連鎖不平衡の情報を自然に取り入れることができ、単座位の解析よりも有利であることが期待されるが、検出力の説明はそれ程単純ではなく単座位の方が有利である状況も多い。ハプロタイプやディプロタイプを用いると原因変異をより良く記述できるマーカーが出現する可能性が高まるが、逆に独立性検定においては自由度が大きくなるために検出力が落ちることがある。検出力の大きさはこのトレードオフによって支配されている。ここで単座位よりハプロタイプやディプロタイプを用いた解析が有利である状況をまとめると、
1.稀なハプロタイプに表現型の原因変異がのっている。
2.原因変異は頻度の高いSNPsで記述されるハプロタイプ上にのっており、ハプロタイプを用いることによってこの頻度が著しく原因変異に近くなる。
3.一つの遺伝子上に同じ表現型に影響する複数の変異がある場合。
の様な状況が考えられる。
【００６９】
２．ディプロタイプ形推定による分割表のp値の区間推定
各個人のディプロタイプ形を推定し、これを用いて分割表による関連解析を行う場合、個人のサンプリングによる揺らぎの他に、推定による揺らぎが存在する。本章では、この推定による揺らぎの検討を行う。
【００７０】
２−１.サンプリング揺らぎと推定揺らぎ
統計学的にp値とは、帰無仮説の元でそのデータ(outcome)が出現する確率である。このp値がある有意水準αを下回った場合、「サンプリングによる揺らぎを考慮しても帰無仮説が成立しているとは考え難い」と解釈し帰無仮説を棄却する。ブートストラップ法は、得られているデータを母集団として復元抽出を行うノンパラメトリックブートストラップ法と、あるモデルのパラメータ推定を得られているデータから行いこれを母数としてランダムサンプリングを行うパラメトリックブートストラップ法に大別される。得られているデータを母集団とするサンプリング揺らぎから、母集団からのサンプリングの揺らぎが評価出来るという仮定の元で、母集団のサンプリングの揺らぎを評価する方法論である。
【００７１】
疾患感受性ゲノム領域を特定する方法論として、一般集団を用いたSNPsによる関連解析がサンプルを集めやすく且つ検出力の高い方法論として注目されている。通常、この解析方法では疾患群と健常群の遺伝子型データが取得される。このとき、SNPsの相に関する情報、即ちハプロタイプの情報は失われている。よってハプロタイプやディプロタイプ形を使った関連解析を、単点のアレルや遺伝子型解析と同様に「分割表の独立性検定」を使って行うにはまず各個人のディプロタイプ形を推定することが必要になる。しかしながら、ディプロタイプ形もあるモデルの元での推定値であり、推定の揺らぎが生じる。そこでハプロタイプやディプロタイプ形を使った「分割表の独立性検定」では、母集団からのサンプリングによる揺らぎとディプロタイプ形推定の揺らぎの2つを考慮することが必要となる。
【００７２】
２−２.p値の区間推定
サンプリング揺らぎと推定揺らぎを評価するために、ブートストラップ法によるp値の分布で真のデータ（outcome）のp値を評価することを考える。ディプロタイプ形推定において、ある与えられたハプロタイプ頻度の集合
【００７３】
【数１６】
の元でベイズの定理を用いて
【００７４】
【数１７】
と事後確率によって評価することを考える。ここで
【００７５】
【数１８】
は2つのハプロタイプの組み合わせで表現したディプロタイプ形であり、ハーディ・ワインベルグ平衡の元では
【００７６】
【数１９】
ここで、δqrはクロネッカーδであり、q=rの時のみ1となり他は0である。Viは個人iの遺伝子型データの集合を表し、このデータにおいて可能な全てのハプロタイプの組み合わせvpについて和
【００７７】
【数２０】
がこのデータが出現する確率(尤度)を表すことになる。式(10)の事後確率に従って各個人のディプロタイプ形をサンプリングすれば、与えられたハプロタイプ頻度の元で推定の不確実性を考慮することができる。しかしながら、通常はハプロタイプ頻度も推定量である。そこでこの推定量の揺らぎをブートストラップ法で評価する。ブートストラップ法を用いることにより、得られたデータを母集団としてハプロタイプ頻度の区間推定を行うことができる。この信頼区間はサンプルデータを母集団としているため、推定によるハプロタイプ頻度の揺らぎを表現していると考えられる。そこで、ブートストラップ法に従い、得られているデータからブートストラップサンプルを作成するステップ、それぞれのブートストラップサンプルを用いてハプロタイプ頻度を最尤推定するステップ、得られたハプロタイプ頻度から各個人のディプロタイプ形事後確率を求めるステップ、この事後確率に従ってディプロタイプ形をサンプリングするステップ、ディプロタイプ形からハプロタイプ及び/又はディプロタイプ形による分割表を作成してp値を計算するステップによって、推定揺らぎを考慮したp値の分布を得ることができる。
【００７８】
本来、ブートストラップ法は、サンプリングによる揺らぎを評価する方法論として発達してきたが、本発明においてはブートストラップサンプルの揺らぎから母集団のサンプリングによる揺らぎを推定するのではなく、直接得られているデータによるハプロタイプ頻度推定揺らぎを評価している。このようにしてデータ（outcome）のp値の信頼区間を構成する。
【００７９】
２−３.注意点
実際の計算上における注意点を下記にまとめた。
先ず、p値の信頼区間の構成においては、ブートストラップサンプルを作成し、ハプロタイプ頻度を求め、これを使ってディプロタイプ形の事後確率を求めている。この時、稀なハプロタイプを持っている個人の事後確率が0になることがある。ブートストラップ法によるサンプリングでこの稀なハプロタイプが含まれない場合には、頻度を0と推定するためこの個人出現確率が0になる。よってハプロタイプやディプロタイプ形の分割表を構成したときに合計が元データと同じ数になるとは限らなくなる。
【００８０】
次に、ブートストラップ法によるサンプリングとディプロタイプ形事後確率に従った2つのサンプリングを行っている。このとき、ディプロタイプ形事後確率に従ったサンプリング数は連鎖不平衡の強さに従って変更することが考えられる。連鎖不平衡が強い領域では殆んどのサンプルのディプロタイプ形事後確率が1つのハプロタイプペアのみに集中し、他は殆んど0になる。この場合にはディプロタイプ形事後確率に従ったサンプリングは1又は2とすることが好ましい。しかしながら、連鎖平衡が弱い領域ではディプロタイプ形事後確率が複数のハプロタイプペアへ割れることが予想され、この場合にはより多い数のサンプリング数とすることが考えられる。
【００８１】
なお、ハプロタイプ頻度の推定揺らぎの評価方法としてブートストラップ法を用いたが、本発明はこれに限定されない。すなわち、ブートストラップ法以外のサンプリング法によってサンプリングを行っても良い。例えば、尤度分布からハプロタイプ頻度をサンプリングする方法が挙げられる。この場合、Markov chain Monte Carlo(MCMC)法などを使ってハプロタイプ頻度に依存する尤度分布を作成することができる。
【００８２】
２−４.発明の構成
上述したp値の区間推定は、図１６に示すようなコンピュータに上述した演算処理を実行させることによって実行することができる。コンピュータは、CPU101（制御手段）、ROM102、RAM103、入力手段104、送信/受信手段105、出力手段106、ハードディスクドライブ(HDD)107及びCD-ROMドライブ108を備える。また、コンピュータは、ROM102、RAM103、HDD107及びパブリックデータベース110等に記録されたデータを検索する検索手段111、検索手段111で検索したデータや入力手段104で入力したデータを加算減算処理する演算手段112を備える。
【００８３】
ここで、本発明に係る遺伝子関連解析プログラムは、例えば、ROM102、RAM103又はHDD107に記憶される。そして、本プログラムに従ってCPU101がコンピュータの上記ハードウェアの駆動制御を行うことで、上記「２−２.p値の区間推定」に説明した情報処理を行いp値の区間推定を実行する。
【００８４】
CPU101は、コンピュータ全体を制御し、遺伝子関連解析処理を実行する。RAM103は、遺伝子関連解析処理を実行する上で必要なデータを一時的に格納する。入力手段104は、キーボードやマウス等であり、遺伝子関連解析処理を実行する上で必要な条件を入力するとき等に操作される。送信／受信手段105は、CPU101の命令に基づいて、通信回線を介してパブリックデータベース110等との間でデータの送受信処理を実行する。出力手段106は、表現型と関連付けられた遺伝子形データ、入力手段104から入力された各種条件、遺伝子関連解析結果等を、CPU101からの命令に基づいて表示処理を実行する。なお、出力手段106としては、コンピュータのディスプレイ、又はプリンターなどが例示される。HDD107は、遺伝子関連解析プログラム、表現型データと関連付けられた遺伝子型データ等を格納し、CPU101の命令に基づいて、格納しているプログラム又はデータ等を読み出し、例えばRAM103に格納する。CD-ROMドライブ108は、CPU101の指示に基づいて、CD-ROM109に格納されている遺伝子関連解析プログラムや、表現型データと関連付けられた遺伝子型データ等から、プログラム又はデータ等を読み出し、例えばRAM103に格納する。なお、表現型データと関連付けられた遺伝子型データとは、個体から採取された２以上の座位に関する遺伝子型データであって、当該個体の表現型データ（例えば、疾患の有無、薬剤抵抗性の有無、副作用の有無等の情報）と関連付けられたものである。表現型データと関連付けられた遺伝子型データは、コンピュータ内部の記憶装置（HDD107、RAM103及びROM102）に格納されていてもよいが、パブリックデータベース110に格納されていてもよい。
【００８５】
本発明に係る遺伝子関連解析プログラムをコンピュータにインストールすることによって、当該コンピュータは遺伝子関連解析を実行する遺伝子関連解析装置として機能することとなる。すなわち、本発明に係る遺伝子関連解析装置は、入力手段104により２以上の座位に関する遺伝子型データを含む個体に関する遺伝子型データセット、並びに、当該個体に関する表現型データを入力する。また、遺伝子関連解析装置において演算手段112は、サンプリング法（例えばブートストラップ法）により上記遺伝子型データセットから複数のサンプリングデータセットを作成し、作成されたサンプリングデータセットに含まれる遺伝子型データセットを用いて当該サンプリングデータセットにおけるハプロタイプ頻度を最尤推定し、得られたハプロタイプ頻度から各個体のディプロタイプ形事後確率を計算し、計算された事後確率からディプロタイプ形をサンプリングし、サンプリングされたディプロタイプ形からハプロタイプ及び/又はディプロタイプ形の分割表を作成し独立性検定を行うことでp値を求め、算出されたp値からp値の分布を算出する処理を行う。そして、遺伝子関連解析装置は、演算手段112で算出されたp値の分布に基づいて、p値の信頼区間及びp値の平均値等の統計量を出力手段106により出力する。
【００８６】
本発明に係る遺伝子関連解析装置は、p値の分布に関する統計量として、例えばp値の信頼区間を出力することが好ましい。また、p値の分布に関する統計量としては、p値の信頼区間の上限値、下限値、p値の平均値、p値の中央値、幾何平均値（logPの平均値のexp）、最頻値、トリム平均を挙げることができる。また、本発明に係る遺伝子関連解析装置において演算手段112は、演算手段112で算出されたp値の信頼区間及び/又はp値の平均値を有意水準αと比較し、当該有意水準αを上回った数をカウントする処理を行い、有意水準を上回った数を出力手段106で出力しても良い。
【００８７】
本発明に係る遺伝子関連解析装置によれば、遺伝子関連解析の結果として、表現型データ（caseとcontrol）に対応するようにp値の信頼区間及びp値の平均値といったp値分布の統計量が出力される。出力されたp値の信頼区間及びp値の平均値等の統計量に基づいて、表現型データを記述するハプロタイプを特定することができる。本発明に係る遺伝子関連解析装置によれば、サンプリングの揺らぎとハプロタイプの推定の揺らぎを考慮した出力結果を得られるため、従来の大標本理論に基づく解析結果よりも信頼性の高い解析を行うことができる。
【００８８】
なお、本発明に係る遺伝子関連解析装置は、上述したp値分布に関する統計量の他に、従前の解析方法に準じて演算した解析結果をあわせて出力するものであっても良い。すなわち、演算手段112は、入力手段104によって入力された上記遺伝子型データセットから、多座位を含むハプロタイプ頻度を最尤推定し、推定されたハプロタイプ頻度から上記遺伝子型データセットに含まれる各個体のディプロタイプ形事後確率を計算し、計算された事後確率からディプロタイプ形をサンプリングし、サンプリングされたディプロタイプ形からハプロタイプ及び/又はディプロタイプ形の分割表を作成し独立性検定を行うことでp値を求める処理を行う。そして、演算手段112によって算出されたp値を出力手段106で出力する。得られたp値は、サンプリングの揺らぎとハプロタイプの推定の揺らぎを考慮した出力結果ではなく、上述したように算出されたp値分布の統計量と比較することができる。
【００８９】
３.シミュレーションデータの解析
本節では、2節において相の情報がある時にハプロタイプやディプロタイプ形の検出力を評価したデータについて、2章の方法論を使ってp値の区間推定を行い、推定の揺らぎを考慮した場合のp値の推定精度、関連解析への影響、及び検出力について調べた結果をまとめた。
【００９０】
３−１. 計算の手続き
今回用いたブートストラップによるp値の区間推定方法は、以下の計算手順に従って実行した。
1.6座位分のハプロタイプ頻度を全サンプルを用いて最尤推定
2.推定されたハプロタイプ頻度を用いて各個人のディプロタイプ形事後確率を計算
3.事後確率から個人のディプロタイプ形をサンプリング
4.サンプリングされたディプロタイプ形からハプロタイプ及びディプロタイプ形の分割表を作成して独立性検定(p値を求める)
5.ノンパラメトリックブートストラップ法に従いサンプル集団から再サンプリングを行い、各再サンプリングデータのハプロタイプ頻度を最尤推定
6.上記5.で得られた各ハプロタイプ頻度に対して、上記2.から4.を繰り返しp値の分布を得る。
7.推定するハプロタイプ頻度を1座位分ずらし、上記を繰り返す。
となる。次節以降に示す計算では、5.のブートストラップ法におけるサンプリングを1000回とし、3.の事後確率からのディプロタイプ形サンプリングはそれぞれ2回ずつ行い、1つのcase-controlデータに対して2000個のハプロタイプ及びディプロタイプ形の分割表を作成してp値を求めた。
【００９１】
３−２. ｐ値の推定精度
表8及び表9に遺伝子型データからディプロタイプ形を推定して得られたハプロタイプ及びディプロタイプ形の分割表のp値とデータ（outcome）のp値との差(絶対値)の平均と分散を、それぞれ100人対100人データ（表8）、10人対100人データ（表9）について示した。ここでMLEとは最尤推定法による単点推定（従前の解析方法に準じた方法）であり、BS meanとは本発明に係る遺伝子関連解析により得られたp値分布における統計量（平均値）であり、intervalとは本発明に係る遺伝子関連解析により得られたp値分布における信頼区間の中に真のデータ（outcome）のp値が含まれていた数である。またCaseとはcase対control-1、Controlはcontrol-1対control-2の比較結果である。
【００９２】
【表８】
【００９３】
【表９】
【００９４】
表8及び表9から判るように、どのケースでも欠測率の増加に伴い平均誤差(mean)及び推定の揺らぎ(se)が増大している。またBS meanの方がMLEよりも推定の平均誤差及び推定の揺らぎとも小さくなっている。これによりハプロタイプ頻度推定の揺らぎを考慮したBS meanの方がMLEのみによる推定よりも精度が高く安定した推定ができていると考えることができる。
【００９５】
信頼区間の構成については、表8を見る限り非常に成功していると言えるだろう。10人対100人の比較である表9では若干カバー率が悪くなっているが、サンプル数が少ない場合には推定の揺らぎが大きくなるため1000回程度のブートストラップ法では高い精度で信頼区間を構成することは難しかったためと考えられる。よって、サンプル数が少ない場合には、ブートストラップ法によるサンプルの回数を適宜増加させることが好ましい。
【００９６】
一方、推定精度のp値依存性を見るために、横軸に真のハプロタイプの分割表のp値をとり、縦軸に推定されたp値(MLE又はBS mean)から真のp値を引いた値をとってプロットした結果を図17-1及び図17-2に示した。なお、図17-1には100人対100人の場合の結果を示し、図17-2には10人対100人の場合の結果を示した。図17-1において(A)及び（B）は、欠測率を0％とした時のcase及びcontrolをそれぞれ示している。図17-1において(C)及び（D）は、欠測率を5％とした時のcase及びcontrolをそれぞれ示している。図17-1において(E)及び（F）は、欠測率を10％とした時のcase及びcontrolをそれぞれ示している。図17-2において(G)及び（H）は、欠測率を0％とした時のcase及びcontrolをそれぞれ示している。図17-2において(I)及び（J）は、欠測率を5％とした時のcase及びcontrolをそれぞれ示している。図17-2において(K)及び（L）は、欠測率を10％とした時のcase及びcontrolをそれぞれ示している。
【００９７】
図17-1及び図17-2から、欠測率が高くなると推定精度が悪くなっていることが分かる。またp値は、
【００９８】
【数２１】
であるためp値が0に近いところでは上方に、p値が1.0に近いところでは下方に分布が広がっている。このため、ある有意水準を定めた場合の検定では、やや保守的になることが予想される。この問題はブートストラップ法で得られたp値の中央値や幾何平均を使うことによって補正される可能性がある。
【００９９】
３−３. 関連解析への影響
３−３−１. 有意となった数
表10及び表11にディプロタイプ形推定によるハプロタイプの分割表においてFisher's exact testのp値が有意となった数をそれぞれ100人対100人データ（表10）及び10人対100人データ（表11）に対して示した。前節と同様にcaseとはcase対control-1、controlはcontrol-1対control-2の比較であり、knownは相が分かっているときのp値を用いた場合、MLEは最尤推定されたハプロタイプ頻度のみを用いた推定、BS meanは本発明に係る遺伝子関連解析により得られたp値分布における統計量（平均値）であり、BS upperは同じく本発明に係る遺伝子関連解析により得られたp値分布における統計量（信頼区間の上限値）であり、missingは遺伝子型データの欠測率である。
【０１００】
【表１０】
【０１０１】
【表１１】
【０１０２】
表10及び表11から、全体的にknownに比べて有意と判定されている数は減っており、また欠測率の増加に伴い有意となる数が減っていることが分かる。これは前節の図17-1及び図17-2で示した様にp値が小さいところでは大きめの予測を行うバイアスのためと考えられる。また、BS upperは欠測率の影響を強く受けており、これはp値の区間推定における信頼区間の幅が広くなったためであろう。
【０１０３】
３−３−２. 判定の精度
有意水準毎にハプロタイプの相が既知の場合のp値(正解データ)と推定値された場合のp値に対して、正解及び不正解の比較を行った。表12及び表13にハプロタイプを使って分割表を作成した場合の結果をそれぞれ100人対100人データ（表12）、10人対100人データ（表13）についてまとめた。表12及び表13においてMLE、BS mean、BS upper（90％信頼区間)については前節と同様である。また、表12及び表13において、TPは真陽性、TNは真陰性、FPは偽陽性及びFNは偽陰性の数を示している。さらに、表12及び表13において、Accとして(TP+TN)/totalの式に従って算出された値も示した。この値は判定の精度を意味している。
【０１０４】
【表１２】
【０１０５】
【表１３】
【０１０６】
表12及び表13から、単点推定であるMLEではFPの数が本発明に係る方法と比べて多くなっており、control対controlの比較ではFPとFNは大体同数程度出現している。BS upperは保守的な指標であるためFNの数が非常に多く、欠測の影響も他の指標に比べて強く出ていることも確認できる。今回の結果からはBS meanが最も精度の高い指標となっており、10人対100人のcase対controlの結果を除いて最も精度が高い。またcontrol対controlの比較においてFPの数が少なく、やや保守的な結果をもたらしているが偽陽性をコントロールしやすい指標であると言える。これらの結果は、本発明に係る遺伝子関連解析装置によってハプロタイプ頻度の推定の揺らぎを考慮してp値分布の統計量を指標にすることが、遺伝子関連解析においてより有効であることを示していると考えられる。
【０１０７】
３−４. 検出力
単座位のアリル及び遺伝子型の分割表で有意となった割合、相が分かっている状態で構築したハプロタイプの分割表で有意となった割合、及び遺伝子型データからハプロタイプの分割表のp値を区間推定した場合に有意となった割合を、100人対100人、10人対100人のデータ毎に表14及び表15にまとめた。なお、有意水準は5%とした。
【０１０８】
【表１４】
【０１０９】
【表１５】
【０１１０】
表14及び表15から、単純に有意となった割合は、MLE、BS mean、BS upperの順に多く、MLEが最も検出力が高くなっているように見える。しかしながら、MLEでは偽陽性もある程度の割合で観察されておりデータの欠測率や連鎖不平衡の大きさによっては正確な検定を行えているとは言えない これに対して、表12におけるcase対controlの比較では、最も保守的なBS upperを検定統計量として採用すると欠測率が高いと単座位よりも検出力が低くなっており、推定の揺らぎによって検出力が下がっていると考えられる。しかしながらBS meanでは相の推定における揺らぎを考慮しても、単座位よりも検出力が高くなっていることが分かる。
【０１１１】
一方で、表12におけるcontrol対controlの比較では、相の推定における揺らぎを考慮することによるバイアスが存在し、やや保守的な結果となっている。逆にサンプル数の少ない表13のケースでは相の情報があるときよりも有意と判定されている割合が多くなっている。しかしながら、表13におけるcontrol対controlの相の推定による偽陽性の増大は僅かであり、今回の方法論ではBS meanを使うとやや保守的になるが、偽陽性の量をある程度コントロールすることが可能であると考えられる。また、単純な相加平均よりも頑健な中央値や幾何平均(logpの平均)を検定統計量として使用すれば、このバイアスが小さいなることが期待できる。
【０１１２】
３−５. 大標本理論に基づく検定統計量の影響
ハプロタイプやディプロタイプを用いた分割表の独立性検定では、単座位のSNPsによる分割表とは異なりカウント数の少ない升目が出現する。そのため大標本理論に基づくχ2検定は適さず、Fisher's exact testの様な小標本検定が必要になる。この様子を表16及び表17に、同じハプロタイプ頻度母数からランダムにそれぞれ100人対100人（表16）、10人対100人ずつ（表17）100ケースサンプリングして得られた相が既知の場合におけるcontrol対controlの比較の結果有意となった割合を示した。
【０１１３】
【表１６】
【０１１４】
【表１７】
【０１１５】
表16ではFisherに比べてχ2検定は保守的になっており、逆に表17では与えられた有意水準以上に有意となる割合が多くなっている。Fisher's exact testの結果は両方ともほぼ有意水準に近い値を示しており、帰無仮説の下で偽陽性をコントロールできていることが分かる。
【０１１６】
５．まとめ
ハプロタイプやディプロタイプを使った関連解析の有効性について主に検出力について説明した。大標本理論の元では近似的な検出力曲線を描くことが可能であり、この結果、原因変異とr2=1となるようなマーカーを用いることが非常に有効であることがわかった。対立因子が2であるSNPsに比べて、ハプロタイプでは複数の対立因子が出現するため、原因変異とr2が高いマーカーが出現する可能性が高まる。特にcommon SNPsマーカーは通常頻度5%以上であるため、原因変異が稀な場合にはハプロタイプを用いた解析が非常に有効であった。このような状況は、表現型データとして薬剤に対する副作用の有無などが考えられ、ハプロタイプやディプロタイプ形を使った独立性検定は、pharmacogenomics(PGx)の解析において特に有効であると考えられる。
【０１１７】
しかしながら、通常のSNPsデータには相の情報がなく直接ハプロタイプを使って関連解析を行うことはできない。そこで遺伝子型データからハプロタイプ頻度やディプロタイプ形を推定することになるが、推定量であるために推定の揺らぎが生じる。よって、実際のハプロタイプやディプロタイプ形を使った関連解析では、この推定の揺らぎを考慮しなければならない。推定の揺らぎを評価するために、一例としてブートストラップ法によるp値の分布を算出してp値の区間推定を行った。シミュレーションデータデータを使った評価では、p値の推定精度は欠損率に依存するが、統計的な有意差を検討することは十分可能であり、また指標によって偽陽性の量をコントロールすることも可能である。よって実際の遺伝子型データから単座位解析では検出することが難しい領域が、本発明に係る遺伝子関連解析装置によって検出することが可能となる。この場合、帰無仮説を棄却する根拠として、「サンプリングの揺らぎと推定の揺らぎを同時に考慮しても、帰無仮説が成立しているとは考えにくい。」とすることができる。
【０１１８】
〔参考例〕
Ａ． 連鎖不平衡指標（Linkage disequilibrium indices）
座位1における所定のアリルをa1、座位2における所定のアリルをa2とし、それぞれのアリル頻度をpa1及びpa2とする。また座位1におけるa1以外のアリルを
【０１１９】
【数２２】
座位2におけるa2以外のアリルを
【０１２０】
【数２３】
と表現する。更に、この2つのアリルのハプロタイプa1a2の頻度をha1a2とすると連鎖不平衡指標D'及びr2は、
【０１２１】
【数２４】
と定義される。
【０１２２】
Ｂ. 大標本理論に基づく検出力の近似式
簡単のため2×2の分割表におけるχ2検定について考える。今、2つの二項分布Bi（n1, p1）とBi（n2, p2）に従う確率変数をそれぞれX1及びX2とすると、
【０１２３】
【数２５】
の時X1及びX2はそれぞれN(n1p1, n1p1(1-p1))及びN(n2p2, n2p2(1-p2))に従う。従って
【０１２４】
【数２６】
は、N(p1, p1(1-p1)/n1)及びN(p2, p2(1-p2)/n2)に従い、
【０１２５】
【数２７】
となる。この時、母比率の検定を帰無仮説H0及び対立仮説H1をそれぞれ
【０１２６】
【数２８】
として検定を行うことを考える。帰無仮説H0の元では、検定統計量
【０１２７】
【数２９】
はN(0,1)に従う。ここで
【０１２８】
【数３０】
である。従ってZ2は自由度1のχ2分布に従い、実際Z2は2×2の分割表のχ2統計量と一致する。
【０１２９】
有意水準αにおける両側検定は式(18)を使うと
【数３１】
と表されるので、
【０１３０】
【数３２】
となる。ここで
【０１３１】
【数３３】
は 式(16)からN(0,1)に従うので、
【０１３２】
【数３４】
と近似する(即ち
【０１３３】
【数３５】
の揺らぎを無視)と、検出力β(n1,p1,n2,p2)(図18参照)は
【０１３４】
【数３６】
と近似される 。ただし、
【０１３５】
【数３７】
であり
【０１３６】
【数３８】
である。
【図面の簡単な説明】
【０１３７】
【図１】D'=1の場合における原因変異とマーカーとの関係を示す模式図である。
【図２】オッズ比をOR=1.5に固定した場合のタイプ1、タイプ2の場合における検出力を示す特性図である。
【図３】オッズ比をOR=3.0に固定した場合のタイプ1、タイプ2の場合における検出力を示す特性図である。
【図４】原因変異のオッズ比を3.0とした場合のマーカーMの見掛けのオッズ比を示す特性図である。
【図５】原因変異のオッズ比を3.0とした場合のマーカーMのマーカーの頻度差f'M-fを示す特性図である。
【図６】シミュレーションデータとして使用したSNPsデータにおける1〜9座位のLDマップ（上三角部分がr2、下三角部分がD'）を示す特性図である。
【図７】原因変異の割合の増分を与えたときのハプロタイプ頻度及びアレル頻度を示す特性図である。
【図８−１】Control-1対control-2で100ケースずつサンプリングした時の単座位のアレル、及び3〜9座位のハプロタイプによる分割表の独立性検定(Fisher's exact test)のp値の分布を示すボックス図である。
【図８−２】Control-1対control-2で100ケースずつサンプリングした時の単座位のアレル、及び3〜9座位のハプロタイプによる分割表の独立性検定(Fisher's exact test)のp値の分布を示すボックス図である。
【図９−１】Control-1対control-2で100ケースずつサンプリングした時の単座位のアレル、及び3〜9座位のハプロタイプによる分割表の独立性検定(Fisher's exact test)で得られたp値に対して有意水準を1%及び5%とした時に有意となった数を示すヒストグラムである。
【図９−２】Control-1対control-2で100ケースずつサンプリングした時の単座位のアレル、及び3〜9座位のハプロタイプによる分割表の独立性検定(Fisher's exact test)で得られたp値に対して有意水準を1%及び5%とした時に有意となった数を示すヒストグラムである。
【図１０−１】Case-1対control-1で100ケースずつサンプリングした時の単座位のアレル、及び3〜9座位のハプロタイプによる分割表の独立性検定(Fisher's exact test)のp値の分布を示すボックス図である。
【図１０−２】Case-1対control-1で100ケースずつサンプリングした時の単座位のアレル、及び3〜9座位のハプロタイプによる分割表の独立性検定(Fisher's exact test)のp値の分布を示すボックス図である。
【図１１−１】Case-1対control-1で100ケースずつサンプリングした時の単座位のアレル、及び3〜9座位のハプロタイプによる分割表の独立性検定(Fisher's exact test)で得られたp値に対して有意水準を1%及び5%とした時に有意となった数を示すヒストグラムである。
【図１１−２】Case-1対control-1で100ケースずつサンプリングした時の単座位のアレル、及び3〜9座位のハプロタイプによる分割表の独立性検定(Fisher's exact test)で得られたp値に対して有意水準を1%及び5%とした時に有意となった数を示すヒストグラムである。
【図１２−１】Case-2対control-1で100ケースずつサンプリングした時の単座位のアレル、及び3〜9座位のハプロタイプによる分割表の独立性検定(Fisher's exact test)のp値の分布を示すボックス図である。
【図１２−２】Case-2対control-1で100ケースずつサンプリングした時の単座位のアレル、及び3〜9座位のハプロタイプによる分割表の独立性検定(Fisher's exact test)のp値の分布を示すボックス図である。
【図１３−１】Case-2対control-1で100ケースずつサンプリングした時の単座位のアレル、及び3〜9座位のハプロタイプによる分割表の独立性検定(Fisher's exact test)で得られたp値に対して有意水準を1%及び5%とした時に有意となった数を示すヒストグラムである。
【図１３−２】Case-2対control-1で100ケースずつサンプリングした時の単座位のアレル、及び3〜9座位のハプロタイプによる分割表の独立性検定(Fisher's exact test)で得られたp値に対して有意水準を1%及び5%とした時に有意となった数を示すヒストグラムである。
【図１４−１】Case-3対control-1で10人対100人サンプリングした時の単座位のアレル、及び3〜9座位のハプロタイプによる分割表の独立性検定(Fisher's exact test)のp値の分布を示すボックス図である。
【図１４−２】Case-3対control-1で10人対100人サンプリングした時の単座位のアレル、及び3〜9座位のハプロタイプによる分割表の独立性検定(Fisher's exact test)のp値の分布を示すボックス図である。
【図１５−１】Case-3対control-1で10人対100人サンプリングした時の単座位のアレル、及び3〜9座位のハプロタイプによる分割表の独立性検定(Fisher's exact test)で得られたp値に対して有意水準を1%及び5%とした時に有意となった数を示すヒストグラムである。
【図１５−２】Case-3対control-1で10人対100人サンプリングした時の単座位のアレル、及び3〜9座位のハプロタイプによる分割表の独立性検定(Fisher's exact test)で得られたp値に対して有意水準を1%及び5%とした時に有意となった数を示すヒストグラムである。
【図１６】本発明を適用した遺伝子関連解析装置のハードウェア構成を示す模式図である。
【図１７−１】横軸に真のハプロタイプの分割表のp値をとり、縦軸に推定されたp値(MLE又はBS mean)から真のp値を引いた値をとってプロットした結果を示す特性図である。
【図１７−２】横軸に真のハプロタイプの分割表のp値をとり、縦軸に推定されたp値(MLE又はBS mean)から真のp値を引いた値をとってプロットした結果を示す特性図である。
【図１８】帰無仮説H0及び対立仮説H1における統計的検定の模式図である。
【符号の説明】
【０１３８】
101…CPU（制御手段）、102…ROM、103…RAM、104…入力手段、105…送信/受信手段、106…出力手段、107…ハードディスクドライブ(HDD)、108…CD-ROMドライブ、109…CD-ROM、110…パブリックデータベース、111…検索手段、112…演算手段
【特許請求の範囲】
【請求項１】
２以上の座位に関する遺伝子型データを含む個体に関する遺伝子型データセット、並びに、当該個体に関する表現型データを入力する入力手段と、
サンプリング法により上記遺伝子型データセットから複数のサンプリングデータセットを作成し、作成されたサンプリングデータセットに含まれる遺伝子型データセットを用いて当該サンプリングデータセットにおけるハプロタイプ頻度を最尤推定し、得られたハプロタイプ頻度から各個体のディプロタイプ形事後確率を計算し、計算された事後確率からディプロタイプ形をサンプリングし、サンプリングされたディプロタイプ形からハプロタイプ及び/又はディプロタイプ形の分割表を作成し独立性検定を行うことでp値を求め、算出されたp値からp値の分布を算出する処理を行う演算手段と、
上記演算手段で算出されたp値の分布に基づいて当該分布の統計量を出力する出力手段とを備える遺伝子関連解析装置。
【請求項２】
上記出力手段は、上記分布の統計量として、少なくともp値の信頼区間の上限値及び平均値を出力することを特徴とする請求項１記載の遺伝子関連解析装置。
【請求項３】
上記演算手段は、上記サンプリング法としてブートストラップ法を適用してサンプリングデータセットを作成することを特徴とする請求項１記載の遺伝子関連解析装置。
【請求項４】
上記演算手段は、上記独立性検定として少なくともフィッシャー正確性検定を適用してｐ値を算出することを特徴とする請求項１記載の遺伝子関連解析装置。
【請求項５】
入力手段、演算手段及び出力手段を備えるコンピュータに、
上記入力手段が、２以上の座位に関する遺伝子型データを含む個体に関する遺伝子型データセット、並びに、当該個体に関する表現型データを入力するステップと、
上記演算手段が、サンプリング法により上記遺伝子型データセットから複数のサンプリングデータセットを作成し、作成されたサンプリングデータセットに含まれる遺伝子型データセットを用いて当該サンプリングデータセットにおけるハプロタイプ頻度を最尤推定し、得られたハプロタイプ頻度から各個体のディプロタイプ形事後確率を計算し、計算された事後確率からディプロタイプ形をサンプリングし、サンプリングされたディプロタイプ形からハプロタイプ及び/又はディプロタイプ形の分割表を作成し独立性検定を行うことでp値を求め、算出されたp値からp値の分布を算出する処理を行うステップと、
上記出力手段が上記演算手段で算出されたp値の分布に基づいて当該分布の統計量を出力するステップと、
を実行させる遺伝子関連解析プログラム。
【請求項６】
上記出力手段は、上記分布の統計量として、少なくともp値の信頼区間の上限値及び平均値を出力することを特徴とする特徴とする請求項５記載の遺伝子関連解析プログラム。
【請求項７】
上記演算手段は、上記サンプリング法としてブートストラップ法を適用してサンプリングデータセットを作成することを特徴とする請求項５記載の遺伝子関連解析プログラム。
【請求項８】
上記演算手段は、上記独立性検定として少なくともフィッシャー正確性検定を適用してｐ値を算出することを特徴とする請求項５記載の遺伝子関連解析プログラム。
【請求項１】
２以上の座位に関する遺伝子型データを含む個体に関する遺伝子型データセット、並びに、当該個体に関する表現型データを入力する入力手段と、
サンプリング法により上記遺伝子型データセットから複数のサンプリングデータセットを作成し、作成されたサンプリングデータセットに含まれる遺伝子型データセットを用いて当該サンプリングデータセットにおけるハプロタイプ頻度を最尤推定し、得られたハプロタイプ頻度から各個体のディプロタイプ形事後確率を計算し、計算された事後確率からディプロタイプ形をサンプリングし、サンプリングされたディプロタイプ形からハプロタイプ及び/又はディプロタイプ形の分割表を作成し独立性検定を行うことでp値を求め、算出されたp値からp値の分布を算出する処理を行う演算手段と、
上記演算手段で算出されたp値の分布に基づいて当該分布の統計量を出力する出力手段とを備える遺伝子関連解析装置。
【請求項２】
上記出力手段は、上記分布の統計量として、少なくともp値の信頼区間の上限値及び平均値を出力することを特徴とする請求項１記載の遺伝子関連解析装置。
【請求項３】
上記演算手段は、上記サンプリング法としてブートストラップ法を適用してサンプリングデータセットを作成することを特徴とする請求項１記載の遺伝子関連解析装置。
【請求項４】
上記演算手段は、上記独立性検定として少なくともフィッシャー正確性検定を適用してｐ値を算出することを特徴とする請求項１記載の遺伝子関連解析装置。
【請求項５】
入力手段、演算手段及び出力手段を備えるコンピュータに、
上記入力手段が、２以上の座位に関する遺伝子型データを含む個体に関する遺伝子型データセット、並びに、当該個体に関する表現型データを入力するステップと、
上記演算手段が、サンプリング法により上記遺伝子型データセットから複数のサンプリングデータセットを作成し、作成されたサンプリングデータセットに含まれる遺伝子型データセットを用いて当該サンプリングデータセットにおけるハプロタイプ頻度を最尤推定し、得られたハプロタイプ頻度から各個体のディプロタイプ形事後確率を計算し、計算された事後確率からディプロタイプ形をサンプリングし、サンプリングされたディプロタイプ形からハプロタイプ及び/又はディプロタイプ形の分割表を作成し独立性検定を行うことでp値を求め、算出されたp値からp値の分布を算出する処理を行うステップと、
上記出力手段が上記演算手段で算出されたp値の分布に基づいて当該分布の統計量を出力するステップと、
を実行させる遺伝子関連解析プログラム。
【請求項６】
上記出力手段は、上記分布の統計量として、少なくともp値の信頼区間の上限値及び平均値を出力することを特徴とする特徴とする請求項５記載の遺伝子関連解析プログラム。
【請求項７】
上記演算手段は、上記サンプリング法としてブートストラップ法を適用してサンプリングデータセットを作成することを特徴とする請求項５記載の遺伝子関連解析プログラム。
【請求項８】
上記演算手段は、上記独立性検定として少なくともフィッシャー正確性検定を適用してｐ値を算出することを特徴とする請求項５記載の遺伝子関連解析プログラム。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８−１】
【図８−２】
【図９−１】
【図９−２】
【図１０−１】
【図１０−２】
【図１１−１】
【図１１−２】
【図１２−１】
【図１２−２】
【図１３−１】
【図１３−２】
【図１４−１】
【図１４−２】
【図１５−１】
【図１５−２】
【図１６】
【図１７−１】
【図１７−２】
【図１８】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８−１】
【図８−２】
【図９−１】
【図９−２】
【図１０−１】
【図１０−２】
【図１１−１】
【図１１−２】
【図１２−１】
【図１２−２】
【図１３−１】
【図１３−２】
【図１４−１】
【図１４−２】
【図１５−１】
【図１５−２】
【図１６】
【図１７−１】
【図１７−２】
【図１８】
【公開番号】特開２００９−６９９１１（Ｐ２００９−６９９１１Ａ）
【公開日】平成２１年４月２日（２００９．４．２）
【国際特許分類】
【出願番号】特願２００７−２３４６３３（Ｐ２００７−２３４６３３）
【出願日】平成１９年９月１０日（２００７．９．１０）
【出願人】（５９２１３１９０６）みずほ情報総研株式会社 (187)
【公開日】平成２１年４月２日（２００９．４．２）
【国際特許分類】
【出願日】平成１９年９月１０日（２００７．９．１０）
【出願人】（５９２１３１９０６）みずほ情報総研株式会社 (187)
[ Back to top ]