(57)【要約】【課題】ユーザ指定の信頼度を反映した精度を実現しつつ処理を高速化させる。【解決手段】本方法は、ゲノム中の指定部分に存在する一塩基多型の集合についてジェノタイプとフェノタイプとの関係を表す統計量を観測データに基づいて算出する工程と、観測データにおいてフェノタイプの値を入れ替えた場合における統計量ｓを算出する処理を、算出した統計量が統計量ｓより高くなる回数の割合が第１の閾値を超えるか若しくは第２の閾値を下回るという事象を検出するまで、又は上記事象が検出されなければ予め指定された回数実行するまで、実行する処理実行工程と、上記割合を一塩基多型の集合に対するｐ値として出力する工程とを含む。ここで第１の閾値及び第２の閾値は、統計量ｓより高い統計量が算出されるケースの出現分布及び予め指定された信頼度に基づき、予め指定された有意水準を表す値から動的に決定されるものである。【選択図】図１
2019053340.tif 000002
【特許請求の範囲】
【請求項１】
ゲノム中の指定部分に存在する一塩基多型の集合についてジェノタイプとフェノタイプとの関係を表す統計量を観測データに基づいて算出するステップと、
前記観測データにおいて前記フェノタイプの値を入れ替えた場合における統計量を算出する処理を、算出した前記統計量が前記観測データに基づき算出した前記統計量より高くなる回数の割合が第１の閾値を超えるか若しくは第２の閾値を下回るという事象を検出するまで、又は前記事象が検出されなければ予め指定された回数実行するまで、実行する処理実行ステップと、
前記割合を前記一塩基多型の集合に対するｐ値として出力するステップと、
を、コンピュータに実行させ、
前記第１の閾値及び前記第２の閾値は、前記観測データに基づき算出した前記統計量より高い統計量が算出されるケースの出現分布及び予め指定された信頼度に基づき、予め指定された有意水準を表す値から動的に決定される
ことを特徴とするプログラム。
【請求項２】
前記第１の閾値が、前記予め指定された有意水準を表す値に対して、前記出現分布において前記予め指定された信頼度に対応する信頼区間の幅の１／２の値を加算した値であり、
前記第２の閾値が、前記予め指定された有意水準を表す値から、前記出現分布において前記予め指定された信頼度に対応する信頼区間の幅の１／２の値を引いた値である
請求項１記載のプログラム。
【請求項３】
前記処理実行ステップにおいて、
前記割合が０である場合には、（予め指定された第２の信頼度に対応する第２の有意水準）＞（算出した前記統計量が前記観測データに基づき算出した前記統計量より低くなる確率を（１−前記予め指定された有意水準）として、前記処理を現在の実行回数だけ実行すると、前記処理の現在の実行回数だけ、算出した前記統計量が前記観測データに基づき算出した前記統計量より低くなる確率）となる第２の事象、又は
前記割合が１である場合には、（前記第２の有意水準）＞（算出した前記統計量が前記観測データに基づき算出した前記統計量より高くなる確率を前記予め指定された有意水準を表す値に対して処理すべき前記一塩基多型の集合の数を乗じた値として、前記処理を現在の実行回数だけ実行すると、前記処理の現在の実行回数だけ、算出した前記統計量が前記観測データに基づき算出した前記統計量より高くなる確率）となる第３の事象
を検出するまで前記処理を実行して、
前記第２の事象を検出すると、１／前記処理の現在の実行回数をｐ値として出力し、
前記第３の事象を検出すると、１をｐ値として出力する
請求項１又は２記載のプログラム。
【請求項４】
前記ｐ値に基づき、前記一塩基多型の集合の採否を決定するステップ
をさらに前記コンピュータに実行させるための請求項１乃至３のいずれか１つ記載のプログラム。
【請求項５】
前記処理実行ステップが、
前記処理を所定回数実行する毎に、前記事象の発生の有無を判断するステップを
含む請求項１乃至４のいずれか１つ記載のプログラム。
【請求項６】
ゲノム中の指定部分に存在する一塩基多型の集合についてジェノタイプとフェノタイプとの関係を表す統計量を観測データに基づいて算出するステップと、
前記観測データにおいて前記フェノタイプの値を入れ替えた場合における統計量を算出する処理を、算出した前記統計量が前記観測データに基づき算出した前記統計量より高くなる回数の割合が第１の閾値を超えるか若しくは第２の閾値を下回るという事象を検出するまで、又は前記事象が検出されなければ予め指定された回数実行するまで、実行する処理実行ステップと、
前記割合を前記一塩基多型の集合に対するｐ値として出力するステップと、
を、コンピュータが実行し、
前記第１の閾値及び前記第２の閾値は、前記観測データに基づき算出した前記統計量より高い統計量が算出されるケースの出現分布及び予め指定された信頼度に基づき、予め指定された有意水準を表す値から動的に決定される
ことを特徴とする情報処理方法。
【請求項７】
ゲノム中の指定部分に存在する一塩基多型の集合についてジェノタイプとフェノタイプとの関係を表す統計量を観測データに基づいて算出する算出手段と、
前記観測データにおいて前記フェノタイプの値を入れ替えた場合における統計量を算出する処理を、算出した前記統計量が前記観測データに基づき算出した前記統計量より高くなる回数の割合が第１の閾値を超えるか若しくは第２の閾値を下回るという事象を検出するまで、又は前記事象が検出されなければ予め指定された回数実行するまで、実行する処理実行手段と、
前記割合を前記一塩基多型の集合に対するｐ値として出力する手段と、
を有し、
前記第１の閾値及び前記第２の閾値は、前記観測データに基づき算出した前記統計量より高い統計量が算出されるケースの出現分布及び予め指定された信頼度に基づき、予め指定された有意水準を表す値から動的に決定される
ことを特徴とする情報処理装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、Gene-based Testに関する。
【背景技術】
【０００２】
ＧＷＡＳ（Genome-Wide Association Studies）解析は、フェノタイプと変異の関係を調べる方法で、ＳＮＰ（一塩基多型。Single-Nucleotide Polymorphism）を中心に調べるSNP-based testと、特定遺伝子上に関連するＳＮＰが集合しているという仮定の下でフェノタイプ（形質）とＳＮＰ集合の関連性を検定するGene-based testとがある。
【０００３】
SNP-based testの場合、特定のフェノタイプとＳＮＰとの関連性を検出するためには、そのＳＮＰが集団中で共通であることが要件となる。しかしながら、集団中のレアなＳＮＰに関しては、実際に関連性があったとしても、統計的に有意な水準で関連性を検出することが出来ない。そのため、特定のフェノタイプとレアなＳＮＰとの関連性を明確に検出するためには、遺伝子上に存在するレアなＳＮＰを集合として観測することでフェノタイプとの関係を見るGene-based testで行うことがよいとされている。
【０００４】
しかしながら、既存のGene-based testに用いられる統計手法を用いて正確なｐ値（帰無仮説が正しいという条件の下で、実際のデータから計算された統計量よりも大きな統計量が観測される確率）を得るためには、permutation test（並べ替え検定とも呼ぶ）を行うことになる。このpermutation testには、多大な計算コストを要することが分かっている。
【０００５】
より具体的には、全ゲノム中の遺伝子部分（エクソン）に存在するレアバリアント（レアなＳＮＰ）とフェノタイプの関係を見ようとすると、通常の統計学的な立場からは、そのレアバリアントとフェノタイプとの関係を統計的に処理した場合に、所定の閾値（例えば０．０１や０．０５）よりも小さいｐ値を有するかどうかを判断することになる。
【０００６】
このｐ値を求める統計的処理を全ゲノムの中に存在する遺伝子部分毎に行うことになり、この遺伝子は全ゲノム中に２００００個ほど存在している。そのため、多重検定という概念を用いると、遺伝子毎の統計学的処理によりその関係性の確からしさを確認するためには、所定の閾値／遺伝子数を基準とすることになる。そうすると、Gene-based testを行う場合のようにpermutation testを実行する場合、有意であるかどうかを見極めるために要する試行回数は最低でも１／（所定の閾値／遺伝子数)となる。
【０００７】
例えば、所定の閾値＝０．０１、遺伝子数＝２００００であれば、permutation testをそれぞれの遺伝子に対して少なくとも２００万回行う必要があり、遺伝子全体ではその２００００倍である４００億回行う必要となる。さらに、複数のフェノタイプについて関連性を計算する場合には、その分permutation testの試行回数は増加することになる。
【０００８】
このため、permutation testの試行回数を削減するための技術が幾つか提案されている。
【０００９】
例えば、以下に示す非特許文献１では、adaptive permutation testとして、最終的に計算されるｐ値が有意ではないと分かったところでテストを中止することで試行回数を削減する技術が開示されている。
【００１０】
さらに、以下に示す非特許文献２では、Ｍ回実行するpermutation testの途中におけるｐ値が、閾値αを超える場合にはテストを途中で打ち切ること、また、ｐ値が、閾値α＋（ｐ値の分布の標準偏差σMの６倍）を超えるか又は閾値α−（ｐ値の分布の標準偏差σMの６倍）を下回る場合にテストを途中で打ち切ることが開示されている。
【００１１】
これらの技術でもある程度の高速化がなされるが、特に非特許文献２では、閾値αを中心として「標準偏差σMの６倍」の２倍という固定的な範囲をテスト継続の条件としているが、このような範囲が適切な範囲であるかどうかは不明である。よって、無駄にテストを継続してしまう若しくは中断してしまう場合もある。
【先行技術文献】
【非特許文献】
【００１２】
【非特許文献１】PLINK/SEQ, PSEQ association analysis [online] [2015年10月1日検索]  インターネット＜https://atgu.mgh.harvard.edu/plinkseq/assoc.shtml＞
【非特許文献２】SNP&Variation Suite v8.4.1. Manual Collapsing Methods for DNA-Sequence Analysis [Online] [2015年10月16日検索]  ＜http://doc.goldenhelix.com/SVS/latest/svsmanual/collapsing_methods.html＞
【発明の概要】
【発明が解決しようとする課題】
【００１３】
従って、本発明の目的は、一側面によれば、ユーザ指定の信頼度を反映した精度を実現しつつpermutation testを高速化させる技術を提供することである。
【課題を解決するための手段】
【００１４】
本発明に係る方法は、（Ａ）ゲノム中の指定部分に存在する一塩基多型の集合についてジェノタイプとフェノタイプとの関係を表す統計量を観測データに基づいて算出するステップと、（Ｂ）観測データにおいてフェノタイプの値を入れ替えた場合における統計量を算出する処理を、算出した統計量が観測データに基づき算出した統計量より高くなる回数の割合が第１の閾値を超えるか若しくは第２の閾値を下回るという事象を検出するまで、又は上記事象が検出されなければ予め指定された回数実行するまで、実行する処理実行ステップと、（Ｃ）上記割合を一塩基多型の集合に対するｐ値として出力するステップと、を含む。そして、第１の閾値及び第２の閾値は、観測データに基づき算出した統計量より高い統計量が算出されるケースの出現分布及び予め指定された信頼度に基づき、予め指定された有意水準を表す値から動的に決定されるものである。
【発明の効果】
【００１５】
一側面によれば、ユーザ指定の信頼度を反映した精度を実現しつつpermutation testを高速化させることができるようになる。
【図面の簡単な説明】
【００１６】
【図１】図１は、実施の形態に係る情報処理装置の構成例を示す図である。
【図２】図２は、本実施の形態に係る処理フローを示す図である。
【図３】図３は、観測データの一例を示す図である。
【図４】図４は、本実施の形態に係る処理フローを示す図である。
【図５】図５は、置換処理の処理フローを示す図である。
【図６】図６は、コンピュータの機能ブロック図である。
【発明を実施するための形態】
【００１７】
本発明の一実施の形態に係る情報処理装置の一例を図１に示す。
【００１８】
本実施の形態に係る情報処理装置１００は、入力部１１０と、観測データ格納部１２０と、第１データ格納部１３０と、ｐ値算出部１４０と、第２データ格納部１５０と、抽出部１６０と、出力部１７０とを有する。
【００１９】
入力部１１０は、ユーザから、各種パラメータ、処理すべき観測データ及びＳＮＰの集合（ＳＮＰｓｅｔと呼ぶ）の指定などを受け付け、第１データ格納部１３０などに格納する。入力部１１０は、情報処理装置１００がネットワークを介して接続されている他のコンピュータから、それらのデータを受信する場合もある。観測データ格納部１２０は、複数の被観測者の各々について、フェノタイプのデータと、ゲノムのデータとを対応付けて格納している。これらのデータについては予め格納されている場合もあれば、入力部１１０から入力される場合もある。
【００２０】
ｐ値算出部１４０は、各ＳＮＰｓｅｔについてｐ値を算出し、算出されたｐ値をＳＮＰｓｅｔに対応付けて第２データ格納部１５０に格納する。このため、ｐ値算出部１４０は、統計量算出部１４１と、置換処理部１４２とを有する。統計量算出部１４１は、予め定められた統計量を算出する。置換処理部１４２は、ブートストラップ法などの方法に従って観測データにおける各被観測者のフェノタイプの値を入れ替える処理等を実行する。
【００２１】
抽出部１６０は、第２データ格納部１５０から、所定の条件を満たすｐ値が算出されたＳＮＰｓｅｔを抽出し、出力部１７０に出力する。出力部１７０は、ユーザ又は他のコンピュータに対して、抽出部１６０の抽出結果を出力する。なお、出力部１７０は、例えば第１データ格納部１３０に格納されているデータに基づき、第２データ格納部１５０に格納されているデータ等の処理途中のデータを出力するようにしても良い。
【００２２】
次に、図２乃至図５を用いて、このような情報処理装置１００が行う処理について詳細に説明する。
【００２３】
入力部１１０は、例えばユーザから、各種パラメータの入力を受け付け、第１データ格納部１３０に格納する（図２：ステップＳ１）。本実施の形態では、１つのＳＮＰｓｅｔについて統計量の最大計算回数Ｂと、第１の信頼度α1と、第２の信頼度α2と、以下で説明する置換処理において算出する統計量の計算回数Ｔと、第１の信頼度α1に対応する有意水準ｐsigとを入力する。なお、ｐsigについては、多重検定を考慮して設定される。すなわち、ｐsig＝（１−α1）／ＳＮＰｓｅｔ数となる。また、Ｂは、少なくとも１／ｐsigとなる。さらに、第２の信頼度α2に対応する有意水準ｐa2を入力するようにしても良い。また、第１の信頼度α1と第２の信頼度α2については同一である場合もある。
【００２４】
また、入力部１１０は、例えばユーザから、処理すべき観測データ及びＳＮＰｓｅｔの指定を受け付け、第１データ格納部１３０又は観測データ格納部１２０に格納する（ステップＳ３）。例えば、観測データ格納部１２０において、どこからどこまでの観測データを用いるか、例えば遺伝子の位置及びその遺伝子におけるＳＮＰの位置等のデータを指定する。
【００２５】
観測データ格納部１２０に格納される観測データの一例を図３に示す。図３では、被観測者毎に、フェノタイプ（ここではｘ及びｙ。例えば身長体重など。）と、ゲノムのデータとが対応付けられて格納されている。例えば、ゲノムには、複数の遺伝子が含まれているが、ここでは遺伝子Ｈに着目しており、その中で変異がある複数の遺伝子座を処理すべきＳＮＰｓｅｔとして指定する。図３の例では変異がある２つのＳＮＰが指定される。なお、フェノタイプについても、多数の種類のフェノタイプが観測データに含まれる場合には、いずれのフェノタイプについて処理するのかについて指定する場合もある。
【００２６】
さらに、遺伝子毎にＳＮＰｓｅｔを指定するのではなく、ある程度の遺伝子のまとまりについてＳＮＰｓｅｔを指定するようにしても良い。また、ここでは、複数のＳＮＰｓｅｔが指定されるものとする。
【００２７】
そうすると、ｐ値算出部１４０は、処理すべきＳＮＰｓｅｔのうち、未処理のＳＮＰｓｅｔを１つ特定する（ステップＳ５）。例えば、処理すべきＳＮＰｓｅｔのリストを用意して上から順番に特定する。
【００２８】
また、ｐ値算出部１４０の統計量算出部１４１は、特定されたＳＮＰｓｅｔについて、観測データそれ自身に関する所定の統計量ｓを算出する（ステップＳ７）。本実施の形態に係る統計量ｓは、ゲノム中の指定遺伝子等に存在する一塩基多型の集合に関してジェノタイプとフェノタイプとの関係を表す統計量である。より具体的には、例えばＳＫＡＴ（Sequential Kernel Association Test）の統計量ｓを算出することとする。
【００２９】
ＳＫＡＴの統計量ｓは、以下のように表される。
ｓ＝（ｙ−μ）TＧＷＧT（ｙ−μ）    （１）
ここでｎを被観測者数、ｍをＳＮＰｓｅｔに含まれるＳＮＰの数として、ｙは、観測されたフェノタイプのｎ次元ベクトルを表し、μは、帰無仮説（すなわち着目しているフェノタイプはジェノタイプとの関係を有さないという仮説）の下におけるフェノタイプの予測平均のｎ次元ベクトルである。また、行列Ｇは、以下のように表される。
Ｇ＝（ｇ1，・・・，ｇi，・・・ｇm）T    （２）
ｇiは、ｉ番目のＳＮＰについてｎ人の被観測者のジェノタイプを含むｎ次元のベクトルである。なお、ｇiは、離散値（例えば０，１，２）でも連続値（例えば０以上２以下）でも良い。
行列Ｗは、以下のように表される。
Ｗ＝diag（ｗ1，・・・，ｗj，・・・ｗm）  （３）
Ｗは、ｍ×ｍの対角行列であり、ｗjは、ｊ番目のＳＮＰについての重み値を表している。
【００３０】
本実施の形態で採用可能な統計量は、ＳＫＡＴの統計量だけに限定されず、例えば、Burden Testの統計量やc-alpha Test統計量などを用いることができる。これらの統計量については、例えば、Seunggeng Lee et al, "Rare-Variant Association Analysis: Study Designs and Statistical Tests", The American journal of Human Genetics 95, 5-23, July 3, 2014、Benjamin et al. “Testing for an Unusual Distribution of Rare Variants” PLoS Genet 7(3), 2011などを参照のこと。
【００３１】
次に、ｐ値算出部１４０は、集合ＳBcにＮＵＬＬを代入し、統計量の計算回数を表す変数Ｂcに０を代入する（ステップＳ９）。
【００３２】
そして、ｐ値算出部１４０は、Ｂcがユーザ指定のＢより小さいか否かを判断する（ステップＳ１１）。この条件が満たされる場合には、処理は端子Ａを介して図４の処理に移行する。
【００３３】
一方、ステップＳ１１の条件が満たされない場合、すなわち計算回数Ｂcがユーザ指定の最大計算回数Ｂ以上となった場合には、ｐ値算出部１４０は、ステップＳ５で特定されたＳＮＰｓｅｔについての推定ｐ値として、図４の処理で計算されるｐBcを、当該ＳＮＰｓｅｔのデータ（例えばＳＮＰｓｅｔに含まれるＳＮＰの位置データ）に対応付けて第２データ格納部１５０に格納する（ステップＳ１３）。
【００３４】
その後、ｐ値算出部１４０は、未処理のＳＮＰｓｅｔが存在するか否かを判断する（ステップＳ１５）。存在する場合には、処理はステップＳ５に戻る。一方、未処理のＳＮＰｓｅｔが存在しない場合には、抽出部１６０は、第２データ格納部１５０において、推定ｐ値がユーザが指定したｐsig未満であるＳＮＰｓｅｔを有意水準を超えるＳＮＰｓｅｔとして抽出して、抽出したＳＮＰｓｅｔを出力部１７０に出力する。出力部１７０は、抽出部１６０の抽出結果を、出力装置（例えば表示装置、印刷装置、他のコンピュータなど）に出力する（ステップＳ１７）。
【００３５】
図２に示した処理は、従来とおおよそ変わらない処理である。処理の高速化については図４及び図５を用いて詳しく述べる。
【００３６】
端子Ａを介して図２のステップＳ１１から移行してくると、置換処理部１４２は、置換処理を実行する（図４：ステップＳ２１）。この置換処理については、図５を用いて詳しく述べる。
【００３７】
置換処理部１４２は、カウンタｔを０に初期化し、集合ＳTにＮＵＬＬを代入する（図５：ステップＳ４１）。そして、置換処理部１４２は、観測データ格納部１２０に格納されている観測データにおいてフェノタイプの入れ替えを実行する（ステップＳ４３）。各被観測者について、復元抽出又は非復元抽出によってフェノタイプの値を抽出して、抽出したフェノタイプの値に入れ替えて、フェノタイプ入れ替え後の観測データを生成する。すなわち、処理対象の被観測者のフェノタイプの値を母集団として、母集団から標本を選ぶ場合に１つの個体を選んだ後にその個体を母集団の中に戻す方法が復元抽出であり、一度選んだ個体を母集団の中に戻さずに標本を選ぶ方法が非復元抽出である。図３の例では、被観測者１００１のフェノタイプの値を、被観測者１００３のフェノタイプの値として設定するといった処理を行う。
【００３８】
そして、置換処理部１４２は、統計量算出部１４１に、フェノタイプ入れ替え後の観測データに関する所定の統計量ｓtを算出させる（ステップＳ４５）。統計量算出部１４１は、処理対象の観測データは異なるが、算出方法はステップＳ７で述べた方法と同じである。
【００３９】
その後、置換処理部１４２は、算出された統計量ｓtを、集合ＳTに入れる（ステップＳ４７）。また、置換処理部１４２は、カウンタｔを１インクリメントする（ステップＳ４９）。
【００４０】
そして、置換処理部１４２は、カウンタｔの値が、ユーザ指定の定数Ｔを超えたか否かを判断する（ステップＳ５１）。カウンタｔの値が、ユーザ指定の定数Ｔ以下であれば、処理はステップＳ４３に戻る。一方、カウンタｔの値が、ユーザ指定の定数Ｔを超えた場合には、処理は呼び出し元の処理に戻る。なお、Ｔは１であってもよいが、ある程度の数値である方が、処理の効率化が図られる。
【００４１】
このような処理を実行すれば、Ｔ個の統計量ｓtが集合ＳTに入れられることになる。
【００４２】
図４の処理の説明に戻って、ｐ値算出部１４０は、集合ＳBcに、集合ＳBc及び集合ＳTを代入する（ステップＳ２３）。また、ｐ値算出部１４０は、変数Ｂcに、集合ＳBcに含まれる統計量の数を代入する（ステップＳ２５）。さらに、ｐ値算出部１４０は、変数ｂに、集合ＳBc内において統計量ｓよりも高い統計量ｓtの数を計数して代入する（ステップＳ２７）。
【００４３】
そして、ｐ値算出部１４０は、ｂ／Ｂcを算出して、推定ｐ値を表す変数ｐBcに代入する（ステップＳ２９）。
【００４４】
さらに、ｐ値算出部１４０は、ガウス分布Ｎ（０，ｐBc（１−ｐBc）／Ｂc）から、信頼度α1の信頼区間の幅ｄaを算出する（ステップＳ３１）。
【００４５】
本実施の形態では、推定ｐ値がユーザが指定したｐsig未満（又は以下）であるＳＮＰｓｅｔを抽出することになるので、明らかにｐsigを下回るＳＮＰｓｅｔ及び明らかにｐsigを上回るようなＳＮＰｓｅｔについては、途中で処理を打ち切っても問題が生じない。しかしながら、どのような条件を満たした場合に、途中で処理を打ち切るかという問題が別に生ずる。
【００４６】
具体的には、ｐsigを中心として一定の幅に入るｐBcが計算された場合には処理を継続し、ｐsigを中心として一定の幅に入らないｐBcが計算された場合には処理を途中で打ち切ることとする。そこで、本実施の形態では、ユーザが指定した信頼度α1に応じた精度が出るように、当該信頼度α1に応じた幅を設定する。
【００４７】
より具体的には、統計量ｓ未満の統計量が算出されると１となり、それ以外であれば０となるバイナリのランダムな変数Ｘを考え、計算された統計量についてのＸの期待値Ｅｘｐ［Ｘ］及び分散Ｖａｒ［Ｘ］は、以下のように表される。
Ｅｘｐ［Ｘ］＝ｐBc＝ｂ／Ｂc
Ｖａｒ［Ｘ］＝Ｅｘｐ［Ｘ2］−Ｅｘｐ［Ｘ］2＝ｐBc（１−ｐBc）
【００４８】
大数の法則から期待値の分散は以下のようになる。
Ｖａｒ［ｐBc］＝ｐBc（１−ｐBc）／Ｂc    （４）
【００４９】
さらに中心極限定理から、ｐBcはガウス分布Ｎ（０，ｐBc（１−ｐBc）／Ｂc）（平均０で、分散は（４）式のガウス分布）に対応するものと考える。そうすると、このガウス分布における信頼度α1に対応する信頼区間の幅ｄaを、信頼度α1に応じた幅として用いることが適当であることが分かる。
【００５０】
ガウス分布Ｎ（０，ｐBc（１−ｐBc）／Ｂc）は、統計量ｓtの計算回数が増加すれば動的に変化する。そして、この分布Ｎ（０，ｐBc（１−ｐBc）／Ｂc）が分かれば、ユーザ指定の信頼度α1に対応する信頼区間を算出する関数は既知であるから、その幅ｄaの詳細な計算方法については説明を省略する。
【００５１】
一方、一時的にｂ＝０又はｂ＝Ｂcとなる状態が生じ得る。このような場合には、分散が０になるので、上で述べたような基準にて処理を途中で終了させるのは妥当ではない。そこで、本実施の形態では、負の二項分布を採用する。ここで、Ｙを試行回数を表す正の整数のランダム変数として、ｐsig_single＝ｐsig＊ＳＮＰｓｅｔ数とする。そうすると、真のｐ値は、最も高くてｐsigであり、最も低くてｐsig_singleとなる。
【００５２】
よって、ｂ＝０の場合には、Ｂcの確率は負の二項分布ＮＢ（Ｙ＝Ｂc；Ｂc，１−ｐsig）と表され、この確率が、最も高くてｐBc＝ｐsigで信頼度α2となる有意水準ｐa2（＝１−α2）より小さければ、処理を途中で停止して、ｐBc＝１／Ｂcと設定する。
【００５３】
すなわち、ｐBcが０である場合には、以下の事象を検出すると、処理を停止する。
予め指定された信頼度α2に対応する有意水準ｐa2＞（算出した統計量ｓtが観測データに基づき算出した統計量ｓより低くなる確率を（１−ｐsig）として、Ｂc回だけ実行すると、Ｂc回だけ、算出した統計量ｓtが観測データに基づき算出した統計量ｓより低くなる確率）
【００５４】
また、ｂ＝Ｂcの場合には、Ｂcの確率は負の二項分布ＮＢ（Ｙ＝Ｂc；Ｂc，ｐsig_single）と表され、この確率が、最も低くてｐBc＝ｐsig_singleで有意水準ｐa2より小さければ、処理を途中で停止して、ｐBc＝１と設定する。
【００５５】
すなわち、ｐBc＝１である場合は、以下の事象を検出すると、処理を停止する。
予め指定された信頼度α2に対応する有意水準ｐa2＞（算出した統計量ｓtが観測データに基づき算出した統計量ｓより高くなる確率をｐsig_singleとして、Ｂc回だけ実行すると、Ｂc回だけ、算出した統計量ｓtが観測データに基づき算出した統計量ｓより高くなる確率）
【００５６】
このような考え方から、ｐ値算出部１４０は、ｐBc＝０且つＮＢ（Ｙ＝Ｂc；Ｂc，１−ｐsig）＜ｐa2という条件を満たすか否かを判断する（ステップＳ３３）。このような条件を満たす場合には、ｐBc＝１／Ｂcとして処理を途中で停止しても良いので、ｐ値算出部１４０は、ｐBcに対して１／Ｂcと設定して（ステップＳ３４）、処理は端子Ｃを介して図２のステップＳ１３に戻る。
【００５７】
一方、ステップＳ３３の条件を満たさない場合には、ｐ値算出部１４０は、ｐBc＝１且つＮＢ（Ｙ＝Ｂc；Ｂc，ｐsig_single＝ｐsig×ＳＮＰｓｅｔ数）＜ｐa2という条件を満たすか否かを判断する（ステップＳ３５）。このような条件を満たす場合には、ｐBc＝１として処理を途中で停止しても良い。よって、処理は端子Ｃを介して図２のステップＳ１３に戻る。
【００５８】
一方、ステップＳ３５の条件を満たさない場合には、ｐ値算出部１４０は、ｐBc＞ｐsig＋ｄa／２又はｐBc＜ｐsig−ｄa／２という条件を満たしているか否かを判断する（ステップＳ３９）。すなわち、ｐBcが、ｐsigを中心に、信頼区間の幅ｄaに入るか否かを判断している。この条件を満たす場合には、処理を途中で停止して良いので、端子Ｃを介して図２のステップＳ１３に戻る。なお、この条件に含まれる不等式は、その関係を維持しつつ変形することが可能である。一方、この条件を満たさない場合には、処理は端子Ｂを介して図２のステップＳ１１に戻る。
【００５９】
このように、ｐBcが、ｐsigを中心に、ユーザ指定の信頼度が反映されて決定される信頼区間の幅ｄaに入る場合には、入念に置換処理を実行して正確なｐ値を算出するようにする。ｄaは、統計量ｓtの算出状況に応じて変化する分布Ｎ（０，ｐBc（１−ｐBc）／Ｂc）に基づき動的に変化する。
【００６０】
ステップＳ３３及びＳ３５の例外を除けば、ｐBcが、ｐsigを中心に、信頼区間の幅ｄaに入らなければ、その時点で次のＳＮＰｓｅｔについての処理に移行するので、それだけ置換処理の実行回数を削減できる。従って、全体の処理速度を高速化できるようになる。
【００６１】
以上述べた実施の形態は一例であって、本発明はこれに限定されるものではない。例えば、図１に示した機能ブロック図は一例であって、プログラムモジュール構成とは一致しない場合もある。処理フローについては、処理結果が変わらない限り、処理順番を入れ替えたり、処理を並列に実行するようにしても良い。本実施の形態の場合、特にＳＮＰｓｅｔ毎に、異なるコンピュータ又はプロセッサコアで並列化しても良い。
【００６２】
また、重要なＳＮＰｓｅｔについては、ｐBc＜ｐsigが成り立つ限り置換処理を繰り返し実行することで正確なｐ値を算出するようにしても良い。重要か否かについては、フェノタイプの入れ替え無しで算出された統計量ｓの値に基づき決定するようにしても良い。特に、フェノタイプの入れ替え無しで算出された統計量ｓの値が、所定の閾値以上となる場合には、本処理の対象外とするようにしても良い。
【００６３】
さらに、図１の情報処理装置１００のように１台のコンピュータで実装するのではなく、クライアント端末から１又は複数の情報処理装置１００にネットワークを介してアクセスして処理を行わせるようにしても良い。複数の情報処理装置１００を用いる場合には、上で述べた機能を分担すればよい。
【００６４】
なお、上で述べた情報処理装置１００は、コンピュータ装置であって、図６に示すように、メモリ２５０１とＣＰＵ（Central Processing Unit）群２５０３とハードディスク・ドライブ（ＨＤＤ：Hard Disk Drive）２５０５と表示装置２５０９に接続される表示制御部２５０７とリムーバブル・ディスク２５１１用のドライブ装置２５１３と入力装置２５１５とネットワークに接続するための通信制御部２５１７とがバス２５１９で接続されている。オペレーティング・システム（ＯＳ：Operating System）及び本実施例における処理を実施するためのアプリケーション・プログラムは、ＨＤＤ２５０５に格納されており、ＣＰＵ群２５０３により実行される際にはＨＤＤ２５０５からメモリ２５０１に読み出される。ＣＰＵ群２５０３は、アプリケーション・プログラムの処理内容に応じて表示制御部２５０７、通信制御部２５１７、ドライブ装置２５１３を制御して、所定の動作を行わせる。また、処理途中のデータについては、主としてメモリ２５０１に格納されるが、ＨＤＤ２５０５に格納されるようにしてもよい。本技術の実施例では、上で述べた処理を実施するためのアプリケーション・プログラムはコンピュータ読み取り可能なリムーバブル・ディスク２５１１に格納されて頒布され、ドライブ装置２５１３からＨＤＤ２５０５にインストールされる。インターネットなどのネットワーク及び通信制御部２５１７を経由して、ＨＤＤ２５０５にインストールされる場合もある。このようなコンピュータ装置は、上で述べたＣＰＵ群２５０３、メモリ２５０１などのハードウエアとＯＳ及びアプリケーション・プログラムなどのプログラムとが有機的に協働することにより、上で述べたような各種機能を実現する。
【００６５】
以上述べた本実施の形態をまとめると、以下のようになる。
【００６６】
本実施の形態に係る方法（検定方法、統計処理方法、テスト方法等とも呼ぶ）は、（Ａ）ゲノム中の指定部分に存在する一塩基多型の集合についてジェノタイプとフェノタイプとの関係を表す統計量を観測データに基づいて算出するステップと、（Ｂ）観測データにおいてフェノタイプの値を入れ替えた場合における統計量を算出する処理を、算出した統計量が観測データに基づき算出した統計量より高くなる回数の割合が第１の閾値を超えるか若しくは第２の閾値を下回るという事象を検出するまで、又は事象が検出されなければ予め指定された回数実行するまで、実行する処理実行ステップと、（Ｃ）上記割合を一塩基多型の集合に対するｐ値として出力するステップとを含む。そして、上で述べた第１の閾値及び第２の閾値は、観測データに基づき算出した統計量より高い統計量が算出されるケースの出現分布及び予め指定された信頼度に基づき、予め指定された有意水準を表す値から動的に決定されるものである。
【００６７】
このように、予め指定された信頼度を反映させた精度が実現されるように第１及び第２の閾値が設定され、第１の閾値を超えるか又は第２の閾値を下回る場合には処理が途中で打ち切られるので、permutation testの高速化も実現されている。さらに、第１及び第２の閾値は、観測データに基づき算出した統計量より高い統計量が算出されるケースの出現分布にも基づくものであるため、第１及び第２の閾値は統計量の算出状況に応じて動的に変化する。すなわち、本実施の形態は、統計量の算出方法に依存せずに適用可能である。
【００６８】
なお、上で述べた第１の閾値が、予め指定された有意水準を表す値に対して、上記出現分布において予め指定された信頼度に対応する信頼区間の幅の１／２の値を加算した値であり、上で述べた第２の閾値が、予め指定された有意水準を表す値から、上記出現分布において予め指定された信頼度に対応する信頼区間の幅の１／２の値を引いた値である場合もある。このように、予め指定された有意水準を表す値を中心に、「予め指定された信頼度に対応する信頼区間の幅」を設定することで、精度からして処理継続が適切な一塩基多型の集合を適切に特定できるようになる。
【００６９】
上で述べた処理実行ステップにおいて、（ｂ１）上記割合が０である場合には、（予め指定された第２の信頼度に対応する第２の有意水準）＞（算出した統計量が観測データに基づき算出した統計量より低くなる確率を（１−予め指定された有意水準）として、上記処理を現在の実行回数だけ実行すると、上記処理の現在の実行回数だけ、算出した統計量が観測データに基づき算出した統計量より低くなる確率）となる第２の事象、又は（ｂ２）上記割合が１である場合には、（上記第２の有意水準）＞（算出した統計量が観測データに基づき算出した統計量より高くなる確率を予め指定された有意水準を表す値に対して処理すべき一塩基多型の集合の数を乗じた値として、上記処理を現在の実行回数だけ実行すると、上記処理の現在の実行回数だけ、算出した統計量が観測データに基づき算出した統計量より高くなる確率）となる第３の事象を検出するまで上記処理を実行するようにしても良い。この場合、上記第２の事象を検出すると、１／処理の現在の実行回数をｐ値として出力し、上記第３の事象を検出すると、１をｐ値として出力するようにしても良い。上記割合が０である場合や１である場合等、特別なケースにも対処できる。
【００７０】
また、本実施の形態に係る方法は、（Ｄ）上記ｐ値に基づき、一塩基多型の集合の採否を決定するステップをさらに含むようにしても良い。これによって有意な一塩基多型の集合を得ることができるようになる。
【００７１】
また、上で述べた処理実行ステップが、上記処理を所定回数実行する毎に、上記事象の発生の有無を判断するステップを含むようにしても良い。このようにすれば、より処理の効率化が図られる。
【００７２】
なお、上で述べたような処理をプロセッサ又はコンピュータに実行させるためのプログラムを作成することができ、当該プログラムは、例えばフレキシブル・ディスク、ＣＤ−ＲＯＭなどの光ディスク、光磁気ディスク、半導体メモリ（例えばＲＯＭ）、ハードディスク等のコンピュータ読み取り可能な記憶媒体又は記憶装置に格納される。なお、処理途中のデータについては、ＲＡＭ等の記憶装置に一時保管される。
【符号の説明】
【００７３】
１００  情報処理装置
１１０  入力部
１２０  観測データ格納部
１３０  第１データ格納部
１４０  ｐ値算出部    １４１  統計量算出部    １４２  置換処理部
１５０  第２データ格納部
１６０  抽出部
１７０  出力部
【図１】
2019053340.tif 000003
【図２】
2019053340.tif 000004
【図３】
2019053340.tif 000005
【図４】
2019053340.tif 000006
【図５】
2019053340.tif 000007
【図６】
2019053340.tif 000008