JPWO2016104688A1 - 特定遺伝子座群又は個別の遺伝子座の遺伝型の判定方法、判定用コンピュータシステム及び判定用プログラム - Google Patents

Abstract

ＭＨＣ等の特定の遺伝子座群に対してマッピングを行ったリード情報を、最適化する手段を確率統計処理的な観点から提供することを目的としてなされた発明である。本発明では、特定の遺伝子座群のアリルに対する各リードの対応が特定されたリード情報に対して、全てのリードにおける個々の特定遺伝子座のアリルに対する期待マッピング数の数値化が行われるステップ、合計期待マッピング数の数値化が行われるステップ、各アリルに割り当てられたリードの割合が算出されるステップ、を繰り返し行って、上記リード情報の最適化をコンピュータにおいて行い、当該最適化情報を基に、特定の遺伝子座群の遺伝型を簡便かつ正確に推定することが可能な方法、コンピュータシステム、及び、コンピュータプログラムが提供される。

Description

本発明は、核酸に基づく遺伝的解析の分野に関する発明であり、より詳しくはヒト遺伝子の高性能シークエンサから導出されるリードの情報に基づき、選択された遺伝子座群又は個別の遺伝子座の遺伝型遺伝型の判定を、簡便かつ高感度に行うための手段を提供する発明である。当該遺伝子座群又は個別の遺伝子座としては、主要組織適合遺伝子抗原（major histocompatibility complex：ＭＨＣ）の遺伝子座群又は個別の遺伝子座が好適例として上げられる。ヒトのＭＨＣは、ヒト白血球型抗原（Human Leucocyte Antigen：ＨＬＡ）である。

特定の遺伝子座群の遺伝型（アリル）は、生体における様々な疾病や器質に関連している。

例えば、ヒトのＭＨＣであるＨＬＡのアリルは、臓器提供の際の適合性判定に重要であるのみならず、薬の感受性、糖尿病、自己免疫疾患、ナルコレプシー等の数多くの疾患との関連が、近年明らかになって来ている。よって、ＨＬＡ遺伝子座群の各ＨＬＡ遺伝座の遺伝型の決定を適切に行うことは、医療現場において極めて重要な事項である。

現在主流のＨＬＡ遺伝型の決定方法は、ＨＬＡアリルの前後にプライマーを設計することでターゲット領域をＰＣＲ増幅し、シークエンスする手法（例えば、特許文献１）である。しかしながら、この手法では、増幅用プライマーの使用の際にオフターゲット増幅による影響等により、偽遺伝子領域が増幅される等、正確なＨＬＡアリルの推定が困難になるという問題点が存在する。特許文献１は、特定のプライマーセットを用いることでこの問題を解決することを目指しているが、対象となるＨＬＡ遺伝子座が実質的には限定されており、ＨＬＡ遺伝子の個人間の相違の多様性を鑑みると、ＰＣＲ増幅を用いている以上、相応の不正確性が生ずると思われる。

さらに、既知のＨＬＡアリル配列を基にしてｃＤＮＡプローブを設計することでＨＬＡ遺伝子座におけるＳＮＰタイピングによりＨＬＡ型を決定する手法が存在する（Itoh Y et al., Immunogenetics 2005）。しかしながら、この手法では新規のアリルを決定することが困難であり、また、アリルにおける個人のゲノムバリエーションが、ＳＮＰタイピング精度に影響を与えるという問題点が存在する。

ＷＯ ２０１４／０６５４１０ Ａ１パンフレット

上記した通りに、高性能シークエンサを何らかの形で用いるＭＨＣ等の遺伝子座群における遺伝型の決定をする際に、ＰＣＲ法等の遺伝子増幅法を行う段階で生ずるオフターゲット領域の増幅等、不正確性の問題は、どれほど適切な増幅用プライマーを用いても常につきまとう問題である。これはどの遺伝子座群についても、ゲノム中に類似の塩基配列を持つ座位が複数存在する、あるいは遺伝的多型が存在する場合、当該遺伝子座群の各遺伝子座の遺伝型の決定の過程で同様の問題が伴うことになる。

そこで、ＭＨＣ遺伝子座群等の、特に複雑な構造の遺伝子座群に対してマッピングを行った高性能シークエンサで得られたリード情報を、各遺伝子座について遺伝型の決定をするために最適化する手段が必要となる。

本発明は、この最適化手段を確率統計処理的な観点から提供することを課題とする。

本発明者らは、この課題の解決に向けて検討を行う過程で、遺伝子座群として、ヒトのＭＨＣ遺伝子座群であるＨＬＡ遺伝子座群を選択して、被験者の検体ＤＮＡから高性能シークエンサにより全ゲノムシークエンスリードデータを得、当該遺伝子座群のＨＬＡアリル参照配列に対するマッピングを行いリードデータとＨＬＡアリル参照配列のマッピング対応情報を得た。当該リードデータに対して各リードの各ＨＬＡアリル毎の期待マッピング数と、各ＨＬＡアリル毎のアリル頻度を求めるために最尤推定処理を行い、さらに好ましくはベイズ推定処理を行うことによって、極めて高精度のＨＬＡ型決定に用いることが可能なＨＬＡアリル毎の期待マッピング数の最適化がなされることを見出し、ゲノム中に類似の塩基配列を持つ座位が複数存在する、あるいは遺伝的多型が多数知られている遺伝子座群の各遺伝子座のアリルを正確に決定できる本発明を完成した。なお、本発明において「高性能シークエンサ」は、本発明の実施に用いることができる、大量のリード情報を比較的短期間で提供することができるシークエンサであり、いわゆる「次世代シークエンサ」を含むものである。本発明を行った時点における次世代シークエンサとしては、例えば、Genome Sequencer FLX（Roche(454)社）、Genome Analyzer IIx、HiSeq2000、HiSeq2500、MiSeq (共にIllumina社)、SOLiD (Applied Biosystem社)、PacBio RS II (Pacific Biosciences 社)等が挙げられるが、これらに限定されるものではなく、現在、将来において提供される高性能シークエンサの全てを含むものである。

本発明完成時点において、リード情報には、概ね、シングルエンドのリード情報と、ペアエンドのリード情報の２種類の形式が認められる。シングルエンドのリード情報とは、リードに対応するＤＮＡ断片の塩基配列の片端の一定長又は可変長（概ね５０〜３００bp程度）についてのリード情報であり、ペアエンドのリード情報とは、当該ＤＮＡ断片の両端の一定長又は可変長（概ね５０〜３００bp程度）についてのリード情報である。技術の進歩に応じてリード情報の内容も日進月歩であるが、本発明においては現在又は将来提供されるリード情報を適用させることが可能である。

選択される遺伝子座群は、上記したように類似の塩基配列を有する遺伝子や擬遺伝子が複数存在し、当該遺伝子座群の各遺伝子座における遺伝的多型（アリル）が多数知られている遺伝子座が、好適な遺伝子座である。当業者の技術常識として、このような遺伝子座群として、ＨＬＡの遺伝子座群等のＭＨＣ遺伝子座群の他に、例えば、シトクロムＰ４５０（Ｃｙｔｏｃｈｒｏｍｅ Ｐ４５０：ＣＹＰ）遺伝子座群、免疫グロブリンをコードする遺伝子座群、T細胞受容体をコードする遺伝子座群、嗅覚受容体をコードする遺伝子群等が知られている。シトクロムＰ４５０は、薬物代謝、解毒に関与する酸化還元酵素ファミリーに属する酵素の総称であり、ヒトでは５７個の機能遺伝子及び５８個の擬遺伝子が知られている。さらに、当該遺伝子座群の遺伝的多型（アリル）はこれまでに２０００種以上が知られており、これらに応じて各種薬物の代謝速度に個人差が現れることが知られている。本明細書に開示した実施例のＭＨＣ遺伝子座群の一つであるＨＬＡ遺伝子座群を、他の複雑な構造の遺伝子座、例えば、シトクロムＰ４５０遺伝子座群に置き換えても同様の良好な効果が得られることは、本発明完成時において明らかであった。

また、個別の遺伝子座とは、遺伝子座群を構成する個別の遺伝子座であり、例えば、ＨＬＡの遺伝子座群であれば、ＨＬＡ-Ａ、ＨＬＡ-Ｂ、ＨＬＡ-Ｃ等が挙げられる。

［Ａ］ 本発明の最適化方法
本発明は、選択された遺伝子座群又は個別の遺伝子座（以下、特定遺伝子座群又は個別の遺伝子座ともいう）のアリル由来のＤＮＡのリード情報が混在したデータのリードの塩基配列に対してマッピングを行うことにより得られる、当該遺伝子座群又は個別の遺伝子座のアリルに対する各リードのマッピング対応が特定されたリード情報（以下、特定遺伝子座群又は個別遺伝子座の対応リード情報ともいう）に対して、下記のステップ（１）〜（６）の全部又は一部が実行されることを特徴とする、遺伝子のリード情報の最適化方法（以下、本発明の最適化方法ともいう）を提供する。本発明の最適化方法は、コンピュータにおいて実行される方法である。

（１） 個々のリードにおける個々の特定遺伝子座群又は個別の遺伝子座のアリルに対する期待マッピング数の数値化が行われるステップ。
「期待マッピング数」は、リード毎に各アリルに対して定義され、後述する「合計期待マッピング数」は、アリル毎に定義され、さらに、「合計期待マッピング数の和」は、当該遺伝子座群又は個別の遺伝子座について定義されるマッピング数である。本発明においては、「マッピング」と「アライメント」は同意義である。

（２） ステップ（１）において数値化された期待マッピング数が当該遺伝子座群又は個別の遺伝子座のアリル毎に合算されて合計期待マッピング数が算出されるステップ。

当該ステップ（２）を数式で一例を示せば、例えば下記式（Ｉ）は、特定遺伝子座群又は個別の遺伝子座のアリルｔから生ずる合計期待マッピング数ｒ_tの、存在度パラメータの現在の推定値に基づいた算出式である。

［式中、Ｚ_ｎｔは、もしリードｎがアリルｔから生じるならば１をとり、それ以外の場合は０である指標変数であり、Ｅ_z［Ｚ_ｎｔ］は、Ｚ_ｎｔの期待値である。ここで、Ｚ_ｎｔの期待値は、Ｚ_nt＝１の事後確率と同値である。 ］

（３） ステップ（２）において算出された合計期待マッピング数が、それぞれ全ての当該遺伝子座群又は個別の遺伝子座のアリルにおける合計期待マッピング数の和で除されて、当該遺伝子座のアリルにマッピングされているリード総量に対して当該遺伝子座群又は個別の遺伝子座の各アリルに割り当てられたリードの割合が算出されるステップ。

当該ステップ（３）を数式で一例を示せば、例えば下記式（II）は、リード総量に対して各特定遺伝子座群又は個別の遺伝子座のアリルに割り当てられたリードの割合Ｅ_ｑ［θ_ｔ］の算出式である。

［式中、ｒ_ｔは特定遺伝子座のアリルt から生ずる合計期待マッピング数、θは特定遺伝子座群又は個別の遺伝子座のアリル上のリードの存在量の頻度を示すパラメータである。］

（４） ステップ（３）において得られたリードの割合が、頻度として個々の当該遺伝子座群又は個別の遺伝子座のアリルに対して割り当てられ、当該割り当て頻度を前提にして、再びステップ（１）により改めて得られた個々のリードにおける個々の当該遺伝子座群又は個別の遺伝子座のアリル毎の期待マッピング数が算出されるステップ。

（５） ステップ（４）において得られた新たな期待マッピング数に対して、再びステップ（２）又は（３）が実行され、当該遺伝子座群又は個別の遺伝子座のアリルにマッピングされているリード総量に対して当該遺伝子座群又は個別の遺伝子座の各アリルに割り当てられたリードの割合が新たに算出されるステップ。

（６） ステップ（４）と（５）が、ステップ（４）において算出されるリード毎の個々の当該遺伝子座群又は個別の遺伝子座のアリルに対する期待マッピング数と、前回のステップ（４）で算出される当該期待マッピング数との間における差が全てのリードについて認められなくなるか、又は、ステップ（５）において算出されるリードの割合の値と、前回のステップ（５）で算出される当該割合の値との間における差が当該遺伝子座群又は個別の遺伝子座の全てのアリルについて認められなくなるまで、繰り返し実行され、収束したリード毎の個々の当該遺伝子座群又は個別の遺伝子座のアリルに対する期待マッピング数、又は、収束した当該遺伝子座群又は個別の遺伝子座のアリル毎のリードの割合の値が、最適化されたデータとして認定されるステップ。

上記の本発明の最適化方法においては、ステップ（１）と（４）において推定される個々のリード毎の特定遺伝子座群又は個別の遺伝子座のアリル毎に対する期待マッピング数と、ステップ（３）と（５）において合算される特定遺伝子座群又は個別の遺伝子座のアリル毎のアリル頻度、のそれぞれが算出される際に、最尤推定を行うためにExpectation Maximizationアルゴリズム（ＥＭアルゴリズム）が使用され、更に好ましくは、ベイズ推定を行うために変分ベイズ法を使用される。

ここで、ＥＭアルゴリズムと変分ベイズ法等の最適化手段を行うための前提を説明する。これは本願の優先権の基礎となる先の出願（特願２０１４−２６５７０４号：以下、先の出願ともいう）における開示について、さらに詳細に解説するものである。本発明においては、ＥＭアルゴリズムと変分ベイズ法の他に、スパースベイズ法、Gibbs sampling法、ＭＣＭＣ法、ＥＰ法、ＰｏｗｅｒＥＰ法等が挙げられる。

本明細書、請求の範囲、及び、図面に記載されたパラメータや数等を示す記号は、発明の開示の便宜ための記号であり、本発明はそれらの記号の種類に全く限定されない。

本発明の推定方法における推定要素として、観測データ、目的パラメータ、及び、観測データと目的パラメータを結びつける潜在変数が挙げられ、例えば、以下のように設定できる。

［観測データ］
観測データ（以下、Ｒ_ｎとも表現される）は、上記の通りに「特定遺伝子座群又は個別の遺伝子座のアリル由来のＤＮＡのリード情報が混在したデータにおける、ＤＮＡリード（以下、リードｎとも表現される）の塩基配列」である。ＤＮＡ混在データは、検体のＤＮＡシークエンスによるＤＮＡのリードを、これとは別に、ヒトＭＨＣであるＨＬＡアリル等の特定遺伝子座群又は個別の遺伝子座のアリルの参照配列をリファレンス配列としてマッピングを行って得られたリード個別の情報の総和として提供されるデータである。当該リファレンス配列は、例えば特定遺伝子座群がＨＬＡ遺伝子座群である場合は、ＩＭＧＴ／ＨＬＡデータベース等から得られるが、当該遺伝子座群について、過去に別の検体でシークエンシング等により決定されたゲノム配列を使用することも可能である。なお、新規の特定遺伝子座の遺伝型、例えば、新規のＨＬＡ型が明らかになった場合には、当該新規遺伝型がデータベース等に逐次繰り入れられていることが好適である。本発明によって明らかになった新規の遺伝型も同様の繰り入れを行うことが好適である。本発明による新規の遺伝型の決定については後述する。

観測データＲ_ｎは、上記の通りにＮ個（Ｎは自然数：リードの本数換算）のＤＮＡ混在データのうちのｎ番目のリードデータにおける塩基配列である。これは、Ｎ個の独立した一様に分布したリードデータとして観測されると仮定される。シングルエンドリードの場合は一本のリードに対して一つの観測データＲ_ｎが当て嵌められるが、ペアエンドリードの場合には、一本の断片に対して両端の塩基配列に対応した２つの観測データが組として当て嵌められる。すなわち、ペアエンドリードの場合は、例えばＲ_ｎａとＲ_ｎｂの組を構成するが、これらは同じ断片から由来する塩基情報であるため、このリードの組を一つの単位として扱うことによって、シングルエンドリードの場合と同様の統計的なモデルのもとで扱うことが容易に可能である（Nariai et al., Bioinformatics:15;29(18) 2013）。具体的には、当該先行文献に基づいた後述の［ペアエンドモデルの完全尤度］の記載の通りに計算対応が可能である。

［目的パラメータ］
目的パラメータは、上記の観測データＲ_ｎを基に推定がなされるパラメータである。本発明は、1つの目的パラメータ（以下、θとも表現される）を伴っている。

目的パラメータθは、ヒトＭＨＣであるＨＬＡアリル等の特定遺伝子座群又は個別の遺伝子座のアリルの頻度を表すベクトルである。例えば、パラメータベクトルθ＝（θ_１，．．．，θ_Ｔ）’（以下、特に断らない限り、ベクトル又は行列における「’」は転置を示す。）として、各特定遺伝子座群又は個別の遺伝子座のアリルについての存在度の分数を

の制約の下で示すことができる。この場合、特定遺伝子座群又は個別の遺伝子座のアリルはＴ個（Ｔは自然数）存在すると仮定され、個々の特定遺伝子座群又は個別の遺伝子座のアリルはｔ（ｔは１以上の整数）としてカウントされる。

目的パラメータθを推定することにより、本発明の目的である特定遺伝子座群の各遺伝子座又は個別の遺伝子座のアリル頻度、あるいは遺伝型の推定を行うことができる。

［潜在変数］
潜在変数は、上記観測データＲ_ｎが、どの特定遺伝子座群又は個別の遺伝子座のアリルから生成されたか、特定遺伝子座群又は個別の遺伝子座のアリルのどの場所から生成されたかを記述するため繰り入れられる非観測変数である。本発明においては、上記の２種の潜在変数（Ｔ_ｎ，Ｓ_ｎ）の当該２種、又は、Ｔ_ｎを単独で繰り入れてパラメータθを算出推定することで、的確にこれらの目的変数の推定を行い、さらにヒトＭＨＣであるＨＬＡ等の特定遺伝子座群の各遺伝子座又は個別の遺伝子座の遺伝型の推定を行うことができる。これらの潜在変数を、上記観測データＲ_ｎからの目的パラメータθの推測工程に、観測データＲ_ｎが依存するように繰り入れて、パラメータθを算出推定することで、各遺伝子座の遺伝型の推定を的確に行うことができる。

上記の潜在変数Ｔ_ｎは、リードｎのヒトＭＨＣであるＨＬＡアリル等の特定遺伝子座群又は個別の遺伝子座のアリル選択に関する、上記θに依存する変数である。Ｔ_ｎ＝ｔは、リードｎが特定遺伝子座群又は個別の遺伝子座のアリルｔから発生することを意味している。

上記の潜在変数Ｓ_ｎは、リードｎの開始位置に関する、上記Ｔ_ｎに依存する変数である。Ｓ_ｎ＝ｓは、リードｎが、位置ｓ（１≦ｓ≦ｌ_ｔ−Ｌ＋１）（ｌ_ｔは、特定遺伝子座群又は個別の遺伝子座のアリルｔの長さであり、Ｌはリード長である）から発生していることを意味している。ここで、一般的にヒトＭＨＣであるＨＬＡ等の本発明の推定方法の対象となる特定遺伝子座群又は個別の遺伝子座のアリルの長さｌ_ｔは、数百塩基長から数万塩基長であり、リードの塩基長よりも長いことが一般的である。また、開始位置ｓが１とは、特定遺伝子座群又は個別の遺伝子座のアリルの最初の塩基からリードが読まれたことを意味する。言い換えればＳ_ｎは、リードをヒトＭＨＣであるＨＬＡ等の特定遺伝子座群又は個別の遺伝子座のアリルの参照配列にマッピングした際の、参照配列における開始位置のことを意味している。

後述するように、特にペアエンドモデルの場合は、例えば、リード間の塩基断片の長さを反映する潜在変数Ｆ_ｎ等を、上記Ｔ_ｎやＳ_ｎと共に遺伝型の推定計算に繰り入れることができる。

［本発明の推定方法の表現］
上記指標を用いた本発明の推定方法は、例えば、「選択された遺伝子座群又は個別の遺伝子座のアリル由来のＤＮＡのリード情報が混在したデータにおけるリード全体の塩基配列を観測データＲとして、個々のリードにおける個々の当該遺伝子座群又は個別の遺伝子座のアリルに対する期待マッピング数を求めるステップ、並びに、目的パラメータである当該遺伝子座群又は個別の遺伝子座のアリル頻度θ（θはＴ次元ベクトル、Ｔは当該遺伝子座群又は個別の遺伝子座のアリルの種類数）の推定値を求めるステップ、を含む被験者由来のＤＮＡのリード情報の当該遺伝子座群又は個別の遺伝子座のアリルへのマッピングを、コンピュータにより最適化する最適化方法において、
上記目的パラメータθ、及び、観測データＲを媒介する潜在変数である、（ａ）リードｎの当該遺伝子座群又は個別の遺伝子座のアリル選択に関するθに依存する変数Ｔ_ｎ、及び、（ｂ）リードｎの開始位置に関するＴ_ｎに依存するＳ_ｎについて、
リードｎの塩基配列を観測データＲ_ｎとして、観測データＲ_ｎからの目的パラメータθの推測工程において観測データＲ_ｎが依存するように、少なくとも（ｉ）変数Ｔ_ｎ及びＳ_ｎ、あるいは、(ii) 変数Ｔ_ｎ、を繰り入れて当該推定値を算出することを特徴とする、最適化方法。 」
として表現され得る。

［シングルエンドモデルの完全尤度］
上記の指標を用いた本発明の最適化方法のパラメータと変数同士の依存関係を反映させた完全尤度（事後同時分布）は、条件付き確率の積として分解される。具体的には、下記式（１）により表される。各記号は、特に断らない限り、前記した通りである。

式（１）において、
ｐ（Ｔ_ｎ＝ｔ｜θ）は、θが所与のもと、リードｎが特定遺伝子座群又は個別の遺伝子座のアリルｔから発生する確率である。この確率は、ｐ（Ｔ_ｎ＝ｔ｜θ）＝θ_ｔとして計算され得る（（１）(ａ)）。

ｐ（Ｓ_ｎ＝ｓ｜Ｔ_ｎ＝ｔ）は、特定遺伝子座群又は個別の遺伝子座のアリルｔが所与のもと、リードｎが位置ｓから発生する確率である。この確率は、ｐ（Ｓ_ｎ＝ｓ｜Ｔ_ｎ＝ｔ）＝１／（ｌ_ｔ−Ｌ＋１）として計算され得る（（１）(ｂ)）。ｌ_tはアリルｔの参照配列の長さ、Ｌはリード長を表す。

ｐ（Ｒ_ｎ｜Ｔ_ｎ＝ｔ，Ｓ_ｎ＝ｓ）は、特定遺伝子座群又は個別の遺伝子座のアリル選択、及び、リードｎの開始位置が所与のもと、リードｎの塩基配列を観測する確率である。ここで、Ｔ_ｎ及びＳ_ｎを要約するための指標変数Ｚ_ｎｔｓ、又は、Ｔ_ｎを要約するための指標変数Ｚ_ｎｔ、を導入することが好適である（（１）(ｃ)）。

Ｚ_ｎｔｓは、（Ｔ_ｎ，Ｓ_ｎ）＝（ｔ，ｓ）の場合、１に等しく、さもなければゼロである。仮に、π_ｎを、リードｎの可能なマッピングについての全（ｔ，ｓ）組のセットとして、その時、各（ｔ，ｓ）∈π_ｎについて、下記式（２）：

（式中、subst ( , , )は、１からシークエンスの置換エラーを引いた数値を取るベースクオリティスコア依存置換確率関数であり、ｒ_ｎ［ｘ］は、リードｎの位置ｘの塩基文字であり、ｑ_ｎ［ｘ］は、リードｎの位置ｘのベースクオリティスコアであり、ｃ_ｔ［ｘ］は、特定遺伝子座群又は個別の遺伝子座のアリルｔの対応するＤＮＡ配列の位置ｘの塩基文字である）
によってリード配列の確率を計算することができる。ベースクオリティスコア置換確率関数,「subst ( , , )」 は、Phredベースクオリティスコアにしたがって決定することも可能であり、ＤＮＡ−Ｓｅｑデータからリードの参照ＤＮＡ配列に対する最も良いアラインメントから見積もることもできる。なお、Phredベースクオリティスコアは、高性能シークエンサからＦＡＳＴＱフォーマットとして出力される塩基配列情報と共に提供される塩基読み取り精度の目安となるスコア、すなわち、シークエンサが出力するエラー率を示すスコアである(Phred quality score)。具体的には、当該スコアＱは、
Ｑ＝−１０log_１０Ｙ（Ｙは、エラー率）、で表される。

一方、Ｚ_ｎｔは、Ｔ_ｎ＝ｔの場合、１に等しく、さもなければゼロである。Ｚ_ｎｔは、上記Ｚ_ｎｔｓを各ｔについて可能なｓについて全てを考慮したものであるため、（２）式から

によってリード配列の確率を計算することができる。

また、上記はシークエンスの置換エラーを考慮した計算式であるが、シークエンスの挿入・欠失エラーを考慮した計算式も、同様に容易に導出可能である（Nariai et al., Bioinformatics:15;29(18)）。

上記式（１）で示した、本発明の推定方法の完全尤度の数式は、以上のように解釈される。

この数式（１）は、シングルエンドモデルによる本発明の推定方法に係る潜在変数の事後確率や事後分布を求める基礎となるものである。

なお、潜在変数Ｚ_ｎｔについて、上記のようにＺ_ｎｔｓに基づく変数としてではなく、Ｚ_ｎｔｓとは独立した潜在変数として、すなわち上記のマッピングしたポジションを示す「ｓ」を全く考慮に入れない潜在変数として設定が可能である。例えば、リードｎがアリルｔのどこかにマッピングされていれば、ポジション「ｓ」は考慮に入れずに、マッピングツールが与えるマッピングスコアのみを利用して、期待マッピング数を算出する等も可能である。以下の開示における潜在変数Ｚ_ｎｔは、原則としてＺ_ｎｔｓに基づく変数として用いられているが、ここに示すＺ_ｎｔｓから独立した変数としてＺ_ｎｔを用いることも可能である。

［ペアエンドモデルの完全尤度］
ペアエンドデータの場合、上記の指標を用いた本発明の最適化方法のパラメータと変数同士の依存関係を反映させた完全尤度（事後同時分布）は、条件付き確率の積として分解される。具体的には、下記式（３）により表される。各記号は、特に断らない限り、前記した通りである。

ここで、Ｆ_ｎはペアエンドリードの組「Ｒ_ｎａとＲ_ｎｂ」のリファレンス配列へのマッピングから推測される塩基断片（フラグメント）の長さである。

式（３）右辺において、ｐ（Ｔ_ｎ＝ｔ｜θ）は、θが所与のもと、リードｎが特定遺伝子座群又は個別の遺伝子座のアリルｔから発生する確率である。この確率は、ｐ（Ｔ_ｎ＝ｔ｜θ）＝θ_ｔとして計算され得る（（３）（ａ））。

ｐ（Ｆ_ｎ＝ｆ｜Ｔ_ｎ＝ｔ）は、遺伝子座のアリルｔが所与のもと、塩基断片の長さｆが発生する確率である。ｄ_Ｆ（ｘ）を、事前に与えられている塩基断片の長さの分布とすると、この確率は、

として計算され得る。ここでｌ_ｔはアリルｔの参照配列の長さ、Ｌはリード長である。塩基断片の長さの分布ｄ_Ｆ（ｘ）は、例えば、平均μ_Ｆ、標準偏差σ_Ｆの正規分布として与える。平均μ_Ｆ、標準偏差σ_Ｆは、塩基断片を作成した際において塩基断片長の分布が実験的に分かっていればその値を指定しても良いが、事前に多数のペアエンドリードをアライメントした結果からこれらのパラメータを推定して指定しても良い。

ｐ（Ｓ_ｎ＝ｓ｜Ｔ_ｎ＝ｔ，Ｆ_ｎ＝ｆ）は、特定遺伝子座群又は個別の遺伝子座のアリルｔが所与のもと、断片長ｆのリードｎの組が位置ｓから発生する確率である。この確率は、ｐ（Ｓ_ｎ＝ｓ｜Ｔ_ｎ＝ｔ，Ｆ_ｎ＝ｆ）＝１／（ｌ_ｔ−ｆ＋１）として計算され得る（（３）（ｂ））。ｌ_ｔはアリルｔの参照配列の長さ、Ｌはリード長、ｆは塩基断片の長さを表す。

ｐ（Ｒ_ｎａ｜Ｔ_ｎ＝ｔ，Ｓ_ｎ＝ｓ）及びｐ（Ｒ_nb｜Ｔ_ｎ＝ｔ，Ｓ_ｎ＝ｓ，Ｆ_ｎ＝ｆ）は、特定遺伝子座群又は個別の遺伝子座のアリル選択、リードｎの組の開始位置、断片長が所与のもと、下記式（４−１，４−２）で計算される：

式中、subst ( , , )は、１からシークエンスの置換エラーを引いた数値を取るベースクオリティスコア依存置換確率関数であり、ｒ_ｎａ［ｘ］は、リードｎの組の一つ目の塩基配列Ｒ_ｎａの位置ｘの塩基文字であり、ｑ_ｎａ［ｘ］は、リードｎの組の一つ目の塩基配列位置ｘのベースクオリティスコアであり、ｃ_ｔａ［ｘ］は、特定遺伝子座群又は個別の遺伝子座のアリルｔのＤＮＡ配列について、リードｎの組の一つ目の塩基配列とマッピングされた位置ｘの塩基文字であり、ｒ_ｎｂ［ｘ］は、リードｎの組の二つ目の塩基配列Ｒ_ｎｂの位置ｘの塩基文字であり、ｑ_ｎｂ［ｘ］は、リードｎの組の一つ目の塩基配列位置ｘのベースクオリティスコアであり、ｃ_ｔｂ［ｘ］は、特定遺伝子座群又は個別の遺伝子座のアリルｔのＤＮＡ配列について、リードｎの組の二つ目の塩基配列とマッピングされた位置ｘの塩基文字である。

上記数式（３）は、ペアエンドモデルによる本発明の推定方法に係る潜在変数の事後確率や事後分布を求める基礎となるものである。

なお、上記のシングルエンドモデルにおける潜在変数Ｚ_ｎｔについての開示は、このペアエンドモデルにおいても適用することができる。

［ハイパーパラメータ］
特に、推定手段として変分ベイズ法等のベイズ推定法を行う際、ハイパーパラメータα_０（０＜α_０）が繰り入れられ計算されることが好適であり、特に０＜α_０≦０．１、もしくは対数尤度の下限を最大化する値であることが好適である。適切な値のハイパーパラメータα_０の繰り入れを行うことにより、はずれ値に強いロバスト性に優れたベイズ推定を行うことが可能となる。

これを、上記式（II）を基に数式にて表すと、例えば、下記式（II）’で表される。

［式中、α_ｔ＝α_０＋ｒ_ｔである。］

ハイパーパラメータα_０は、ベイズ推定における枠組みにおいて加味される定数である。すなわち、特定遺伝子座群又は個別の遺伝子座のアリル上のリードの存在量を示すパラメータθは、ベイズ推定における枠組みにおいては事後分布として推定することが可能であり、当該θの事前分布としてディリクレ分布（式（III））：

［式中、Ｃは定数であり、Π^Ｔ _ｔ＝１θ_ｔ＝１、Ｔは検討する特定遺伝子座群又は個別の遺伝子座のアリルの数であり、α_ｔはハイパーパラメータである。］
を仮定する。パラメータθの複雑さ（θ_ｔ＞０となる個数）をコントロールするハイパーパラメータα_０を、測定データの対数周辺尤度を最大化するように選択する。

そして、測定データを前提としてθの事後分布を予測することは、潜在変数に対する積分を必要とし、閉形式で計算し難い。そこで、潜在変数とパラメータθの因子分解を仮定することによって、事後確率分布の近似式を得て導出される式が、上記式（II）’である。

ハイパーパラメータについては、改めて記載を行う。

上記した本発明の最適化方法は、特定遺伝子座群又は個別遺伝子座の対応リード情報であれば特に限定されずに用いることができる。例えば、特定遺伝子座群のある遺伝子座又は個別の遺伝子座に対応したプライマーを用いて調製された遺伝子増幅産物に対して、高性能シークエンサによる処理を経て得られた当該遺伝子座対応リード情報であっても良いし、当該遺伝子座対応リード情報をさらに当該遺伝子座のアリルとのマッピングを行ったリード情報であっても良い。しかしながら、本発明の最適化方法は、このような事前の特定遺伝子座群のある遺伝子座又は個別の遺伝子座における被験者の遺伝子の増幅工程を行わずに、被験者の遺伝子検体を高性能シークエンサで処理を行って得られる全ゲノムリード情報に対して、当該遺伝子座群又は個別の遺伝子座のアリルに対するマッピングを行った特定遺伝子座群又は個別遺伝子座の対応リード情報に対しても適用することができる。

本発明においては、特定遺伝子座群又は個別遺伝子座の対応リード情報を取得するためのマッピングが、下記のステップ（ａ）及び（ｂ）により実行されることを特徴とする、本発明の最適化方法を提供する。以下、これらのステップを行うプロセスを、「本発明の特定遺伝子座群又は個別の遺伝子座のマッピングプロセス」ともいう。

（ａ） 遺伝子シークエンサにより得られた被験者のリードの塩基配列情報において、ヒト遺伝子の塩基配列に対するマッピングが行われ、特定遺伝子座群又は個別の遺伝子座のアリルにマッピングされたリードが抽出されるステップ、

この最初のマッピングの対象はヒト全ゲノム配列であることが好適な態様の一つであり、慣習的に特定の人物（解析対象として選択された人物、あるいは特定の人物の組み合わせ）のゲノム配列が対象になる。通常は、国際ゲノムコンソーシアム等の機関が決定したゲノム配列である。この１回目のマッピングによって、当該遺伝子座群又は個別の遺伝子座のアリルに関係の無いリードを除くことができる。なお、上記のヒト全ゲノム配列以外に、例えば、ターゲットシークエンス、Ｅｘоｍｅシークエンス、ＲＮＡシークエンス、ＰａｃＢｉｏ ＲＳ II、Ｏｘｆｏｒｄ Ｎａｎｏｐｏｒｅ等のロングリードシークエンスデータ等も上記マッピング対象配列として用いることができる。

（ｂ） ステップ（ａ）により抽出された特定遺伝子座群又は個別の遺伝子座のアリルにマッピングされたリードの配列情報を、データベースに登録されている当該遺伝子座群又は個別の遺伝子座のアリルの塩基配列に対してマッピングを行うことにより、マッピングされたリードが当該遺伝子座群又は個別の遺伝子座のアリル毎に抽出され、当該遺伝子座群又は個別の遺伝子座のアリルに対する各リードのマッピング対応が特定されたリード情報が得られるステップ、

この２回目のマッピングの対象は、データベースに登録されている特定遺伝子座群又は個別の遺伝子座の全てのアリルの遺伝子配列である。これにより、いわば仮の特定遺伝子座群又は個別の遺伝子座の対応リード情報を得ることができる。

本発明の特定遺伝子座群又は個別の遺伝子座のマッピングプロセスに対応する、上記のステップ（ｂ）におけるマッピングは、一つのリードが複数の特定遺伝子座群又は個別の遺伝子座のアリルに対してマッピングされることを許容するものであることが好適である。この時点で機械的にマッピング対象が絞り込まれると、導出される特定遺伝子座群又は個別遺伝子座の対応リード情報に関して不適切なバイアスを折り込んでしまう可能性が強くなる。

上述した通りに、上記ステップ（ａ）において用いる遺伝子シークエンサにより得られたリードの配列情報は、各ＤＮＡ断片の両端からの読み取り（それぞれ５０〜３００ｂｐ程度）を行うペアエンドの配列情報であってもよい。各ＤＮＡ断片の片側からの読み取り（５０〜３００ｂｐ程度）を行うシングルエンドの配列情報でも良いが、ペアエンドの方が、１本のＤＮＡフラグメント（通常、３００〜１０００ｂｐ程度）に対応した配列情報が当該リードの両端で結ばれる範囲で特定され、より精度の高いマッピングが期待され、結果としてより精度の高いアリルの推定が可能となる。

また、ステップ（ａ）の特定遺伝子座群又は個別の遺伝子座のアリルにマッピングされたリードに加えて、ヒト遺伝子に対してマッピングがなされなかったリードが併せて抽出され、これが（ｂ）ステップの再マッピングの対象とされることが好ましい。上記ステップ（ａ）のマッピングの対象ゲノムは、特定の人物もしくは複数の人物の組み合わせのゲノムであり、具体的検出対象とは合わない場合が想定されるからである。例えば、特定遺伝子座群がヒトＭＨＣであるＨＬＡ遺伝子座群である場合に、マッピングの対象ゲノムが西洋人のゲノムであり、被験者が日本人の場合には、ＨＬＡ遺伝子座群の配列がマッピング対象ゲノムとは大きく異なる可能性があり、当該被験者のＨＬＡ遺伝子座群由来のリードが当該対象ゲノムにマッピングされない可能性がある。これを担保するために上記の処理が行われる。

［Ｂ］ 本発明の判定方法
本発明の最適化方法で得られた特定遺伝子座群又は個別遺伝子座の対応リード情報を、そのまま特定遺伝子座群の各遺伝子座又は個別の遺伝子座の遺伝型の判定指標として用いることが可能である。特に、高機能シークエンサの使用の段階で、特定遺伝子座群のある遺伝子座又は個別の遺伝子座に対応した遺伝子増幅用プライマーを用いたＰＣＲ法等の遺伝子増幅操作により、当該遺伝子座のアリルに対応したリードに絞り込んでいる場合は、その傾向が認められる。その場合には、当該遺伝子座の対応リード情報におけるリードの割合から当該遺伝子座のアリル毎のリードの個別深度が算出され、当該個別深度の大きな当該遺伝子座のアリルから順に２個以内を、被験者の当該遺伝子座の遺伝型として決定を行うことができる（本発明の判定方法）。しかしながら、このような場合であっても、結果に対する再検討を行い、偽陽性の可能性を出来る限り除くことが好適である。

また、上記の本発明の特定遺伝子座群又は個別遺伝子座のマッピングプロセスを行った場合を含め、事前の当該遺伝子座群の遺伝子座に対する遺伝子増幅法を用いた絞り込みを行わなかった場合においては、結果に対する再検討を行う必要性はより高くなる。

事前の遺伝子座に対する遺伝子増幅法を用いた絞り込みを行わない手法の場合（例えば、全ゲノムシークエンスをリファレンス配列とする場合）は、下記の再検討プロセスを行うことがより望ましい。

この本発明の判定方法の好適な態様は、上述した本発明の最適化方法により得られた特定遺伝子座群又は個別遺伝子座の対応リード情報が、個別深度として算出されて再評価が行われることにより、極めて確度が高くなった特定遺伝子座群の各遺伝子座の又は個別の遺伝子座の遺伝型の判定方法である。

すなわち本発明により、本発明の最適化方法により得られた特定遺伝子座群又は個別の遺伝子座のアリルのリードの割合から当該遺伝子座群又は個別の遺伝子座のアリル毎のリードの個別深度が算出され、当該遺伝子座群の各遺伝子座又は個別の遺伝子座について当該個別深度の大きな当該遺伝子座のアリルから順に２個以内について選択され、当該遺伝子座の遺伝型の要素として決定がなされる特定遺伝子座群の各遺伝子座又は個別の遺伝子座の遺伝型の判定方法において、全リード深度の５〜５０％、好ましくは１０〜３０％のいずれかの頻度数が棄却閾値として設定され、当該閾値以下の個別深度の特定遺伝子座群又は個別の遺伝子座のアリルは当該遺伝子座群又は個別の遺伝子座の遺伝型決定の要素から除外されることを特徴とする特定遺伝子座群の各遺伝子座又は個別の遺伝子座の遺伝型の判定方法、が提供される。

「全リード深度」とは、高性能シークエンサによる遺伝子検体の処理に伴い算出される、塩基毎の対応するＤＮＡ断片の数のことを意味するもので、被験対象の全ゲノムに対してどの程度の高性能シークエンサにおける重複読み取りが行われたかの平均値を示す指標である。具体的には、シークエンサで読まれた全リードに含まれる総塩基数を、ゲノムの長さ（ヒトの全ゲノムは３０億塩基）で除した値である。例えば、１００ｂｐシングルエンドリードが９億リードあれば、全リード深度は「１００×９億／３０億 ＝３０」で、「３０×」となる。ここで「全ゲノムの長さ」は、全ゲノムシークエンスでは無い場合、例えば被験対象が「特定のＨＬＡ遺伝子座」である場合には、特定のＨＬＡ遺伝子座における全塩基数となる。

「個別深度」とは、本発明においては、特定遺伝子座群又は個別の遺伝子座のアリルの塩基に対応するリードの重なりの数を数えたときに、平均どれくらいのリードが重なっているかを示す指標であって、本発明の最適化方法により得られた特定遺伝子座群又は個別の遺伝子座のアリルのリードの割合から算出され得る。

具体的には、「特定遺伝子座群又は個別の遺伝子座のアリルに割り振られたリード数」は、「当該遺伝子座群又は個別の遺伝子座のアリルのリードの割合× 当該遺伝子座群又は個別の遺伝子座のアリルにマッピングされた総リード数 」で算出される。

上述したように、本発明の最適化方法に従って「特定遺伝子座群又は個別の遺伝子座のアリルのリードの割合」が算出される。また、「特定遺伝子座群又は個別の遺伝子座のアリルにマッピングされた総リード数」は、例えば、「本発明の特定遺伝子座群又は個別遺伝子座のマッピングプロセス」において導出される。よって、「各特定遺伝子座群又は個別の遺伝子座のアリルに割り振られたリード数」は、本発明の最適化方法が行われることによって算出される。

そして「個別深度」は、
「特定遺伝子座群又は個別の遺伝子座のアリルにマッピングされたリードの総塩基数 ／ 当該遺伝子座群又は個別の遺伝子座のアリルのリファレンス配列の塩基数 」で算出される。

「特定遺伝子座群又は個別の遺伝子座のアリルにマッピングされたリードの総塩基数」は、各リードの平均的な塩基数は既知数であるから、上述した「各特定遺伝子座群又は個別の遺伝子座のアリルに割り振られたリード数」に、当該平均的塩基数を乗ずることによって算出される。「特定遺伝子座群又は個別の遺伝子座のアリルのリファレンス配列の塩基数」は、特定遺伝子座群又は個別の遺伝子座のアリル毎に異なる既知数である。よって、上記特定遺伝子座群又は個別の遺伝子座のアリルの「個別深度」は、上記したように、本発明の最適化方法により得られた特定遺伝子座群又は個別の遺伝子座のアリルのリードの割合から算出される。

以上をまとめると、各特定遺伝子座群又は個別の遺伝子座のアリルにおける個別深度「ｄ_ｔ」は、次式（ＩＶ）によって算出される。

［式中、Ｎは総リード数であり、ｃ_ｎはリードｎが含む塩基数であり、Ｅ［Ｚ_ｎｔ］はリードｎの遺伝子座のアリルｔへの期待マッピング数であり、ｌ_ｔは各遺伝子座のアリルの参照配列の長さ（塩基数）である。ｔは１からＴ（遺伝子座におけるアリルの総数）まで、ｎは１からＮまでを採ることができる。］

なお、ここに示したＥ［Ｚ_ｎｔ］は、上述したＥ［Ｚ_ｎｔｓ］を、各アリルｔについて可能なｓについて全て足し合わせたものとしても計算できる。

上記のように「棄却深度」が全リード深度を基に設定されることにより、個別深度が小さい特定遺伝子座群又は個別の遺伝子座のアリルが遺伝型の決定要素から除外され、いわば偽陽性にあたる特定遺伝子座群又は個別の遺伝子座のアリル候補の除外によって本発明の判定方法の確度を高めることができる。

棄却深度は、上述のように全リード深度の５〜５０％、好ましくは同１０〜３０％、のいずれかの頻度数として選択することが可能である。棄却深度が小さければ、特定遺伝子座群又は個別の遺伝子座のアリルが、本発明の判定方法における特定遺伝子座群又は個別の遺伝子座の遺伝型の決定要素の候補となる機会が多くなるが、偽陽性を拾ってしまう危険性も増すことになる。逆に棄却深度が大きければ、偽陽性を拾う可能性は少なくなるけれども、真に被験者の特定遺伝子座群又は個別の遺伝子座の遺伝型を示すアリルを棄却してしまう可能性が高くなる。

「個別深度の大きな特定遺伝子座群の各遺伝子座又は個別の遺伝子座のアリルから順に２個以内」、すなわち特定遺伝子座群の各遺伝子座又は個別の遺伝子座のアリルの個数の最大値を「２個」としたのは、異なる２個の遺伝子座のアリルがこれらのヘテロ接合と決定され、これを超える個数分はノイズとして排除されるという意味である。そして、当該個数が「１個」であれば当該遺伝子座は当該アリルのホモ接合、もしくはヘテロ接合かつもう一方のアリルは未知、と決定され、「０個」であれば該当する当該遺伝子座のアリルは未知と決定されるものである。

さらに具体的な本発明の判定方法の態様を挙げれば、上述した特定遺伝子座群又は個別の遺伝子座の遺伝型決定の要素からの除外が行われた後、下記（ｉ）又は（ii）の決定がなされることが好適である。
（ｉ） 特定遺伝子座群の各遺伝子座又は個別の遺伝子座について遺伝型決定の対象が１個のアリルについては、当該１個のアリルの個別深度が前記棄却閾値の２倍以上の場合には、当該１個のアリルはホモ接合と決定がなされ、若しくは、前記棄却閾値の２倍より小さい場合はヘテロ接合であると決定がなされ、
（ii） 特定遺伝子座群の各遺伝子座又は個別の遺伝子座について遺伝型決定の対象が２個のアリルについては、個別深度が大きな方が小さい方の２倍未満である場合には、両アリルはヘテロ接合であるとの決定がなされ、若しくは、個別深度が大きな方が小さい方の２倍以上である場合には、大きな方のアリルはホモ接合であるとの決定がなされる。

本発明の判定方法をこのような態様として行うことにより、一層的確な被験者の特定遺伝子座群の各遺伝子座又は個別の遺伝子座のアリルの決定が可能となる。

なお、決定されるべきアリルが新規である場合には、当該新規アリルと最も近い既知のアリルがまず決定され、当該既知アリルの塩基配列と決定されるべき新規アリルの塩基配列の置換、挿入、欠失等による差分を認識することで、当該新規アリルの塩基配列決定をすることができる。当該新規アリルの塩基配列は、新しい遺伝型として、対象データベース等に逐次登録を行うことが好適である。

以下、本発明の最適化方法と判定方法を、「本発明の方法」と総称することもある。

［Ｃ］ 本発明のコンピュータシステム
本発明のコンピュータシステムは、上述した本発明の方法を行う手段となるシステムであり、特に断らない限りは同一の用語は概念として重複する。「アルゴリズム」とは、コンピュータ分野の一般的な概念と同じく、問題を解くための手順を定式化した形で表現したものを意味する。

本発明のコンピュータシステムは、通常のコンピュータシステムに関わるハードウエアを備えることができる。すなわち、通常ハードディスクドライブに該当する「記録部」、ＣＰＵに相当する「演算処理部」の他、例えば、ＲＡＭに相当する「一時記憶部」、キーボード、マウス、タッチパネル等に相当する「操作部」、ディスプレイに相当する「表示部」、操作部に応じたシリアル又はパラレルインターフェース等に相当する「出入力インターフェース（ＩＦ）部」、ビデオメモリとＤ／Ａ変換部を備え、表示部のビデオ方式に応じたアナログ信号を出力する「通信インターフェース（ＩＦ）部」を備えている。当該通信ＩＦ部では、外部の情報、特に、ヒトゲノムデータベース等のヒトゲノム情報とデータ交換を行うことができる。

以下においては特に断らない限り、本発明のコンピュータシステムの「演算処理部」が行う処理として説明する。「演算処理部」は、「操作部」が操作されて「通信ＩＦ部」を介して、特にヒトゲノムデータベースのデータを取得して「記録部」に記録し、適宜当該「記録部」からデータを「一時記憶部」に読み出し、所定の処理を行った後、その結果を再度「記録部」に記録する。当該「演算処理部」は、「操作部」の操作を促す画面データや処理結果を表示する画面データを作成し、入力ＩＦ部のビデオＲＡＭを介して、これらの画像を「表示部」に表示する。本発明のプログラムは、用時又は予め「記録部」に記録、あるいは、外部のハードウエア資源に記録されており、必要に応じて「演算処理部」において、記載されたアルゴリズムに従った演算処理が行われる。

本発明のコンピュータシステムは、特定遺伝子座群又は個別遺伝子座の対応リード情報を最適化するコンピュータシステムであって、記録部と演算処理部を備え、下記の処理（Ａ）〜（Ｇ）の全て又は一部；
（Ａ） 当該記録部には、被験者由来のＤＮＡのリード情報が、リードの配列及びリードのマッピング先である当該遺伝子座群又は個別の遺伝子座のアリルのデータとして記録されており、
（Ｂ） 当該演算処理部では、前記記録部の情報に基づいて、個々のリードにおける個々の当該遺伝子座群又は個別の遺伝子座のアリルに対する期待マッピング数の数値化処理が実行され、
（Ｃ） 上記処理（Ｂ）において数値化された期待マッピング数が当該遺伝子座群又は個別の遺伝子座のアリル毎に合算されて合計期待マッピング数が算出され、
（Ｄ） 上記処理（Ｃ）において算出された合計期待マッピング数が、それぞれ全ての当該遺伝子座群又は個別の遺伝子座のアリルにおける合計期待マッピング数の和で除されて、当該遺伝子座群又は個別の遺伝子座のアリルにマッピングされているリード総量に対して当該遺伝子座群又は個別の遺伝子座の各アリルに割り当てられたリードの割合が算出される処理が実行され、
（Ｅ） 上記処理（Ｃ）において算出されたリードの割合が、頻度として個々の当該遺伝子座群又は個別の遺伝子座のアリルに対して割り当てられ、当該割り当て頻度を前提にして、再び上記処理（Ｂ）により改めて算出された個々のリードにおける個々の当該遺伝子座群又は個別の遺伝子座のアリル毎の期待マッピング数が算出される処理が実行され、
（Ｆ） 上記処理（Ｅ）により算出された新たな期待マッピング数に対して、再び上記処理（Ｃ）又は（Ｄ）が実行されて、当該遺伝子座群又は個別の遺伝子座のアリルにマッピングされているリード総量に対して当該遺伝子座群又は個別の遺伝子座の各アリルに割り当てられたリードの割合が新たに算出される処理が実行され、
（Ｇ） 上記処理（Ｅ）と（Ｆ）が、処理（Ｅ）において算出されるリード毎の個々の当該遺伝子座群又は個別の遺伝子座のアリルに対する期待マッピング数と、前回の処理（Ｅ）において算出される当該期待マッピング数との間における差が全てのリードについて認められなくなるか、又は、処理（Ｆ）において算出されるリードの割合の値と、前回の処理（Ｆ）で算出される当該割合の値との差が当該遺伝子座群又は個別の遺伝子座の全てのアリルについて認められなくなるまで、繰り返し実行され、収束したリード毎の個々の当該遺伝子座群又は個別の遺伝子座のアリルに対する期待マッピング数、又は、収束した当該遺伝子座群又は個別の遺伝子座のアリル毎のリードの割合の値（アリル頻度）が、最適化されたデータとして認定が行われる；
が実行されることを特徴とするコンピュータシステムである。

上記処理（Ｂ）〜（Ｇ）は、全ての特定遺伝子座群又は個別の遺伝子座のアリルに対して包括的に実行されることが好適である。この包括的な実行は、全ての当該遺伝子座群又は個別の遺伝子座のアリルに関して、一緒に全リードのマッピングの最適化アルゴリズムを実行することを意味するものである。

さらに本発明のコンピュータシステムにおいて、被験者のデータベースに登録されている遺伝子のリード情報の当該遺伝子座群又は個別の遺伝子座のアリルに対するマッピングは、下記の（ａ）及び（ｂ）の処理により実行され得る。下記の処理（ａ）及び（ｂ）で示される過程は、特定遺伝子座群を対象とする場合は、その全ての遺伝子座について同時に行われることが好適である。

（ａ） 遺伝子シークエンサにより得られた被験者のリードの配列情報に対して、ヒト全遺伝子の参照塩基配列に対するマッピングの後、当該遺伝子座群の各遺伝子座又は個別の遺伝子座にマッピングされたリードが抽出される処理。

（ｂ） 前記（ａ）の処理により抽出された当該遺伝子座群又は個別の遺伝子座にマッピングされたリードの配列情報に対して、データベースに登録されている当該遺伝子座群又は個別の遺伝子座のアリルの塩基配列とのマッピングの後、マッピングされたリードが当該遺伝子座群又は個別の遺伝子座の各アリルに対する各リードのマッピング対応及びマッピング状態、すなわちリード配列のリファレンス配列におけるマッピング位置、リード配列とリファレンス配列の差異、及びマッピングスコアが特定されたリード情報が得られる処理。

処理（ｂ）において実行されるマッピングは、一つのリードが複数の特定遺伝子座群又は個別の遺伝子座のアリルに対してマッピングされることを許容することが好適である。

処理（ａ）の特定遺伝子座群又は個別の遺伝子座のアリルにマッピングされたリードに加えて、ヒト遺伝子に対してマッピングがなされなかったリードが併せて抽出処理され、これが（ｂ）処理の再マッピングの対象とされることが好適である。

さらに本発明は、被験者の特定遺伝子座群の各遺伝子座の遺伝型の判定を行うコンピュータシステムであって、記録部と演算処理部を備え、下記（α）〜（δ）の処理の全部又は一部；
（α） 当該記録部には、本発明の最適化方法により得られた、被験者の当該遺伝子座群又は個別の遺伝子座のアリル頻度、及び、全リード深度、が少なくとも記録されており；
（β） 当該演算処理部では、前記記録部の当該遺伝子座群又は個別の遺伝子座のアリル頻度を基とする、当該遺伝子座群又は個別の遺伝子座のアリル毎の個別深度への算出処理、及び、個々の当該遺伝子座群又は個別の遺伝子座のアリルに対する算出された当該個別深度の割り振り処理が実行され、
（γ） 棄却閾値として設定されている、全リード深度の平均の５〜５０％、好ましくは１０〜３０％のいずれかの頻度数に対して、当該数値以下の個別深度の当該遺伝子座群又は個別の遺伝子座のアリルは当該遺伝子座群の各遺伝子座又は個別の遺伝子座の遺伝型決定の要素から除外される処理が実行され、
（δ）：
（δ）−１ （γ）の除外処理の実行の後、当該遺伝子座群の各遺伝子座又は個別の遺伝子座の遺伝型決定の対象が１個のアリルについては、当該１個のアリルの個別深度が前記棄却閾値の２倍以上である場合には、当該アリルはホモ接合と決定がなされる処理が実行され、又は、前記棄却閾値の２倍より小さい場合はヘテロ接合であると決定がなされる処理が実行され、
（δ）−２ （γ）の除外処理の実行の後、当該遺伝子座群の各遺伝子座又は個別の遺伝子座の遺伝型決定の対象が２個のアリルについては、個別深度が大きな方が小さい方の２倍未満である場合には、両アリルはヘテロ接合であるとの決定がなされる処理が実行され、又は、個別深度が大きな方が小さい方の２倍以上である場合には、大きな方のアリルはホモ接合であるとの決定がなされる処理が実行される、
ことを特徴とするコンピュータシステムを提供する。

上記に加えて、例えば（γ）の除外処理の実行の後、特定遺伝子座群の各遺伝子座又は個別の遺伝子座の遺伝型決定の対象が０個のアリルについては、当該アリルに基づく当該遺伝子座の遺伝型の決定はなされない処理が実行される設定を、本発明のコンピュータシステムにおいて行うことも可能である。

なお、これらのコンピュータシステムのカテゴリーは「物」であり、「装置」として置き換えることも可能である。

［Ｄ］ 本発明のプログラム
本発明のプログラムは、本発明のコンピュータシステムに本発明の方法を実行させるためのアルゴリズムを備えたコンピュータプログラムであり、特に断らない限りは同一の用語は概念として重複する。

本発明のプログラムは、特定遺伝子座群又は個別遺伝子座の対応リード情報を最適化するコンピュータプログラムであって、コンピュータに下記の第１の機能〜第７の機能の全て又は一部；
（Ａ） 被験者由来のＤＮＡのリード情報が、リードの配列及びリードのマッピング先である当該遺伝子座群又は個別の遺伝子座のアリルのデータとして記録されている記録部から、当該リード情報を読み出す、第１の機能、
（Ｂ） 上記第１の機能により読み出したリード情報に基づいて、個々のリードにおける個々の当該遺伝子座群又は個別の遺伝子座のアリルに対する期待マッピング数の数値化処理を実行する、第２の機能、
（Ｃ） 上記第２の機能により数値化した期待マッピング数が、当該遺伝子座群又は個別の遺伝子座のアリル毎に合算されて合計期待マッピング数を算出する、第３の機能、
（Ｄ） 上記第３の機能により算出した合計期待マッピング数を、それぞれ全ての当該遺伝子座群又は個別の遺伝子座のアリルにおける合計期待マッピング数の和で除して、当該遺伝子座のアリルにマッピングされているリード総量に対して当該遺伝子座群又は個別の遺伝子座の各アリルに割り当てられたリードの割合を算出する、第４の機能、
（Ｅ） 上記第４の機能により算出したリードの割合を、頻度として個々の当該遺伝子座群又は個別の遺伝子座のアリルに対して割り当て、当該割り当て頻度を前提にして、再び第２の機能で改めて算出した、個々のリードにおける個々の当該遺伝子座群又は個別の遺伝子座のアリル毎の期待マッピング数を算出する、第５の機能、
（Ｆ） 上記第５の機能により算出した新たな期待マッピング数に対して、再び上記第３の機能又は第４の機能を実行して、当該遺伝子座群又は個別の遺伝子座のアリルにマッピングされているリード総量に対して当該遺伝子座群又は個別の遺伝子座の各アリルに割り当てられたリードの割合を新たに算出する、第６の機能、
（Ｇ） 上記第５の機能と第６の機能を、第５の機能の実行により算出するリード毎の個々の当該遺伝子座群又は個別の遺伝子座のアリルに対する期待マッピング数と、前回の第５の機能の実行により算出する当該期待マッピング数との間における差が全てのリードについて認められなくなるか、又は、上記第６の機能の実行により算出するリードの割合の値と、前回の第６の機能の実行により算出する当該割合の値との差が全ての当該遺伝子座群又は個別の遺伝子座のアリルについて認められなくなるまで、繰り返し実行し、収束したリード毎の個々の当該遺伝子座群又は個別の遺伝子座のアリルに対する期待マッピング数、又は、収束した当該遺伝子座群又は個別の遺伝子座のアリル毎のリードの割合の値を、最適化されたデータとして認定する、第７の機能；
を実現させるアルゴリズムが含まれることを特徴とする、コンピュータプログラムである。

さらに本発明は、被験者のデータベースに登録されている遺伝子のリード情報の特定遺伝子座群又は個別の遺伝子座のアリルに対するマッピングを、下記（ａ）及び（ｂ）に従って行う機能をコンピュータにおいて実現するアルゴリズムが含まれることを特徴とする、本発明のプログラムを提供する。
（ａ） 遺伝子シークエンサにより得られた被験者のリードの配列情報に対して、ヒト遺伝子の塩基配列に対するマッピングの後、当該遺伝子座群又は個別の遺伝子座のアリルにマッピングされたリードを抽出する機能。
（ｂ） 機能（ａ）により抽出された当該遺伝子座群又は個別の遺伝子座のアリルにマッピングされたリードの配列情報に対して、データベースに登録されている当該遺伝子座群又は個別の遺伝子座のアリルの塩基配列とのマッピングの後、当該遺伝子座群又は個別の遺伝子座のアリルに対する各リードのマッピング対応及びマッピング状態、すなわちリード配列のリファレンス配列におけるマッピング位置、リード配列とリファレンス配列の差異、及びマッピングスコアが特定されたリード情報を得る機能。

さらに本発明は、被験者の特定遺伝子座群の各遺伝子座又は個別の遺伝子座の遺伝型の判定を行うコンピュータプログラムであって、下記（α）〜（δ）の機能をコンピュータに実現させるためのアルゴリズムが含まれることを特徴とする、コンピュータプログラムを提供する。
（α） 前記のコンピュータプログラムの実行により得られた、当該遺伝子座群又は個別の遺伝子座のアリルのリードの割合、及び、全リード深度、を少なくとも読み出す、機能α。
（β） 前記機能αの実行により読み出した当該遺伝子座群又は個別の遺伝子座のアリルのリードの割合から、当該遺伝子座群又は個別の遺伝子座のアリル毎の個別深度への算出処理を実行し、個々の当該遺伝子座群又は個別の遺伝子座のアリルに対して算出された当該個別深度を割り振る処理を実行する、機能β。
（γ） 棄却閾値として全リード深度の５〜５０％、好ましくは１０〜３０％のいずれかの頻度数を設定し、前記機能Ｂの実行により特定された当該数値以下の個別深度の当該遺伝子座群又は個別の遺伝子座のアリルを、当該遺伝子座群の各遺伝子座又は個別の遺伝子座の遺伝型決定の要素から除外する処理を実行する、機能γ。
（δ） 下記（δ）−１及び（δ）−２に示す機能δ。
（δ）−１ 前記機能γの除外処理の実行の後、当該遺伝子座群の各遺伝子座又は個別の遺伝子座の遺伝型決定の対象が１個のアリルについては、当該１個のアリルの個別深度が前記棄却閾値の２倍以上である場合には、このアリルをホモ接合と決定し、又は、前記棄却閾値の２倍より小さい場合はヘテロ接合であると決定する処理を実行し、
（δ）−２ 前記機能γの除外処理の実行の後、当該遺伝子座群の各遺伝子座又は個別の遺伝子座の遺伝型決定の対象が２個のアリルについては、個別深度が大きな方が小さい方の２倍未満である場合には、両アリルはヘテロ接合であると決定し、又は、個別深度が大きな方が小さい方の２倍以上である場合には、大きな方のアリルはホモ接合であると決定する処理を実行する。

上記の処理に、例えば前記機能γの除外処理の実行の後、特定遺伝子座群の各遺伝子座又は個別の遺伝子座の遺伝型決定の対象が０個のアリルについては、当該アリルに基づく遺伝型の決定を行わない処理の実行をする機能を積極的に加えることもできる。

本発明のコンピュータプログラムは、例えば、Ｃ言語、Java（登録商標）、Perl、Python等で記載することが可能である。

本発明はさらに、本発明のプログラムが記録されていることを特徴とする、コンピュータにおいて読み取り可能な記録媒体又はコンピュータに接続し得る記録媒体（以下、本発明の記録媒体ともいう）を提供する。これらの記録媒体としては、フレキシブルディスク、フラッシュメモリ、ハードディスク等の磁気的媒体、ＣＤ、ＤＶＤ、ＢＤ等の光学的媒体、ＭＯ、ＭＤ等の磁気光学的媒体等が挙げられ、特に限定されるものではない。本発明のコンピュータシステムの典型は、本発明のプログラムを実行することを特徴とするものである。

本発明により、高性能シークエンサのリード情報に由来する特定遺伝子座群又は個別の遺伝子座のマッピング情報を、当該遺伝子座群の各遺伝型又は個別の遺伝子座の遺伝型の判定に向けて最適化する方法、最適化コンピュータシステム、及び、最適化用コンピュータプログラム、並びに当該最適化手段を用いた当該遺伝子座群の各遺伝子座の遺伝型又は個別の遺伝子座の遺伝型の判定方法、判定用コンピュータシステム、及び、判定用コンピュータプログラムが提供される。これらの本発明により、極めて的確かつ簡便にゲノム中に類似の塩基配列を持つ座位が複数存在する、あるいは遺伝的多型が多数知られている遺伝子座群の各遺伝子座又は個別の遺伝子座の遺伝型の判定を行うことが可能である。当該遺伝子座群としては、例えば、ヒトＭＨＣであるＨＬＡの遺伝子座群、シトクロムＰ４５０遺伝子座群等、免疫グロブリンをコードする遺伝子座群、T細胞受容体をコードする遺伝子座群、嗅覚受容体をコードする遺伝子座群が挙げられる。

本発明の処理の流れを示した図面である。 ＨＬＡ対応リード情報の最適化工程の一態様を、潜在変数Ｚ_ｎｔを用いて示すフローシートである。 ＨＬＡ対応リード情報の最適化工程の一態様を、潜在変数Ｚ_ｎｔｓを用いて示す、ＥＭアルゴリズムのフローシートである。 ＨＬＡ対応リード情報の最適化工程の一態様を、潜在変数Ｚ_ｎｔｓを用いて示す、変分ベイズ法のフローシートである。 上記の最適化工程を経て得られた「各ＨＬＡアリルに割り当てられたリードの割合」を基に行う、ＨＬＡタイピング工程の一態様を示したフローシートである。 シミュレーションデータの範囲の種々の深度で、４ケタ解像度でのＨＬＡ型の予測精度を示した図面である。 本発明のシステムと他の方法によってＨＬＡ-Ａのアリル頻度を検討した結果を示すグラフである。 本発明のシステムと他の方法によってＨＬＡ-Ｂのアリル頻度を検討した結果を示すグラフである。 本発明のシステムと他の方法によってＨＬＡ-Ｃのアリル頻度を検討した結果を示すグラフである。

［Ａ］ ＥＭアルゴリズムと変分ベイズ法を用いた本発明の推定方法
ここでは、先の出願において開示されたＥＭアルゴリズムと変分ベイズ法の内容についてさらに詳細に説明を行う。用いられる記号は、特に断らない限り前述した通りである。

（１）ＥＭ（expectation-maximization)アルゴリズム
ＥＭアルゴリズムは、潜在変数が存在する確率モデルにおいてパラメータの最尤推定値を求めるための手法である。本発明の推定方法においては、例えば、下記の内容で行うことができる。

すなわちＥＭ法による本発明の推定方法は、下記（１）〜（５）のステップを行うことを特徴とする本発明の推定方法として例示される（この方法を「本発明のＥＭ法」ともいう）。

（ａ）所与されたθ_ｔの初期値に基づいて、Ｚ_ｎｔｓ＝１又はＺ_ｎｔ＝１の事後確率の第１の更新値を算出し、さらに、当該Ｚ_ｎｔｓ＝１又はＺ_ｎｔ＝１の事後確率の第１の更新値に基づいてθ_ｔの最尤推定値の第１の更新値を算出するステップ、あるいは、（ｂ）所与されたＺ_ｎｔｓ＝１又はＺ_ｎｔ＝１の事後確率の初期値に基づいて、θ_ｔの最尤推定値の第１の更新値を算出するステップ、
（ｃ） 直前のステップ（ｄ）（ただし、初回はステップ（ａ）又はステップ（ｂ）である）により算出されたθ_ｔの最尤推定値の更新値に基づいて、新たにＺ_ｎｔｓ＝１又はＺ_ｎｔ＝１の事後確率の更新値を算出するステップ、
（ｄ） ステップ（ｃ）において算出されたＺ_ｎｔｓ＝１又はＺ_ｎｔ＝１の事後確率の更新値に基づいて、θ_ｔの最尤推定値の更新値を新たに算出するステップ、
（ｅ） （ｉ）ステップ（ｃ）において算出されたＺ_ｎｔｓ＝１又はＺ_ｎｔ＝１、及び、ステップ（ｄ）において算出されたθ_ｔに基づいて対数尤度を計算して、対数尤度の収束性を評価するステップ、
（ii）ステップ（ｃ）で算出されたＺ_ｎｔｓ＝１又はＺ_ｎｔ＝１の事後確率の更新値の収束性を評価するステップ、あるいは、
（iii）ステップ（ｄ）で算出されたθ_ｔの最尤推定値の更新値の収束性を評価するステップであって、
収束が認められれば、それぞれのステップにおけるθ_ｔを最終推定値として決定し、収束が認められなければ、ステップ（ｃ）、（ｄ）、及び、（ｅ）の繰り返しを決定する。

上記の（ａ）ステップは、繰り返し最適化を行う本法の初期値を定義するステップである。

Ｚ_ｎｔｓ又はＺ_ｎｔの初期値は、ｎについては１からＮの各々の数字であり、ｔについては１からＴについての各々の数字である。Ｚ_ｎｔｓにおいて用いられるｓについては、１≦ｓ≦ｌ_ｔ−Ｌ＋１、ｌ_ｔは、特定遺伝子座のアリルｔの長さ、Ｌはリード長である。そして、あるリードが複数箇所にマッピングしている場合、同様に確からしいと仮定して、Ｚ_ｎｔｓ＝１、Ｚ_ｎｔ＝１の事後確率の初期値は１／［リードｎがマッピングしている合計箇所］とすることが好ましい。θ_ｔの初期値は、各々の特定遺伝子座のアリル頻度は等しいと仮定して、１／Ｔとすることが好ましい。なお、被験者の母集団が特定の人種の集団である等、特定遺伝型のアリル頻度の事前情報が既知の場合には、これをθ_ｔの初期値として用いることも可能である。

（ｂ）ステップは、現在の推定値（最初は（ａ）ステップで得られた初期化の値）に基づいて、Ｚ_ｎｔｓ＝１の事後確率ｒ_ｎｔｓ又はＺ_ｎｔ＝１の事後確率ｒ_ｎｔを算出するステップである。

すなわち、Ｐ（Ｚ_ｎｔｓ＝１｜θ^＊ _ｔ）又はＰ（Ｚ_ｎｔ＝１｜θ^＊ _ｔ）を算出するステップであり（＊は更新されたことを表している。以下、同様。）、

として計算される。

ここで、シングルエンドリードの場合、

（ｌｏｇρ_ｎｔｓを算出するための右辺第１項は上記式（１）（ａ）に、第２項は上記式（１）（ｃ）に、及び、第３項は上記式（２）について、対数計算することにより算出することができる）

また、ペアエンドリードの場合、

（ｌｏｇρ_ｎｔｓを算出するための右辺中括弧内第１項は上記式（３）（ａ）に、第２項は上記式（３）（ａ’）に、第３項は上記式（３）（ｂ）に、第４項は上記式（４−１）に、及び、第５項は上記式（４−２）について、対数計算することにより算出することができる）

一方Ｚ_ｎｔは、上記Ｚ_ｎｔｓを各ｔについて可能なｓについて全てを考慮したものであるため、

として計算される。

ここまでが本発明のＥＭ法のＥステップに相当する。

（ｃ）ステップは、本発明のＥＭ法のＭステップに相当するステップであり、（ｂ）ステップで得られたＺ_ｎｔｓ＝１又はＺ_ｎｔ＝１の事後確率であるｒ_ｎｔｓ又はｒ_ｎｔに基づき、θ_ｔの最尤推定値を求めるステップである。以下の通り、対数尤度を最大にするθ_ｔを算出する。

又は、

（ｄ）ステップと（ｅ）ステップは、上記の通りである。収束基準の好適な一例として、存在量パラメータθ_ｔ＞１０^−７である特定遺伝子座群又は個別の遺伝子座のアリル頻度についての、相対的変化１０^−３が収束基準として使用されるが、異なる収束基準を用いることも可能である。

なお、通常は（ｅ）ステップの「ステップ（ｃ）、（ｄ）、及び、（ｅ）の繰り返し」が、１回以上は行われるが、この繰り返し工程を１回も行わない、又は十分に収束していない状態でステップを終了して、特定遺伝子座群又は個別の遺伝子座のアリルのタイピングを行う事も可能である。

（２）変分ベイズ法
変分ベイズ法は、ベイズ推定法においてパラメータの事後確率分布を推定することによって、よりノイズに強い安定した推定を行うための方法である。

式（１）で示した、本発明の推定方法の完全尤度（事後同時分布）の数式は、全ての可能な潜在変数Ｚにわたる積分を伴い、解析的に解けない。そのため、当該完全尤度（事後同時分布）において、潜在的変数及びモデルパラメータの因子分解を、

と仮定することによって、近似値を求めるものである。

上述したように、θの事前分布について、本発明では、ディリクレ分布、

［式中、Π^Ｔ _ｔ＝１θ_ｔ＝１、Ｔは検討する特定遺伝子座群のアリルの数であり、α_ｔはハイパーパラメータであり、

そしてг(・)は、ガンマ関数である。］
を用いる。本発明において、ヒトＭＨＣであるＨＬＡ等の特定遺伝子座群又は個別の遺伝子座のアリル頻度の相対的差異について予備知識はないという仮定に基づいて、全ての特定遺伝子座群又は個別の遺伝子座のアリルについて単一のハイパーパラメータα_０が用いられることが好適である。当該単一ハイパーパラメータα_０は、目的パラメータの複雑さ、すなわち θ_ｔ ＞ ０ となる個数を制御する。α_０ ≧１の時、α_０−１は、特定遺伝子座のアリルに割り当てられるリードの事前カウントとして解釈され得、α_０＜１の時、当該事前分布は、特定遺伝子座群又は個別の遺伝子座のアリルのいくつかがゼロの傾向を与える。具体的には、０＜α_０≦０．１、もしくは対数周辺尤度の下限を最大化するα_０を設定することが好適である。アリル頻度の事前情報が既知の場合は、α_ｔをそれぞれのｔについてアリル頻度の事前情報の高いものから順に大きい値が与えられるようにα_ｔ ＞ ０を重み付けして設定することも可能である。ハイパーパラメータα_０が所与のもと、対数周辺尤度の下限は、変分ベイズ推定アルゴリズムによって繰り返しにより最大化される。

変分ベイズ法による本発明の推定方法は、下記（ａ）〜（ｄ）のステップを行うことを特徴とする本発明の推定方法として例示される（この方法を「本発明の変分ベイズ法」ともいう）。

（ａ） 所与された特定遺伝子座群又は個別の遺伝子座のアリルｔのアリル頻度についての予備知識の分布を示すハイパーパラメータα_ｔの初期値に基づくθ_ｔの事後分布の更新値に基づいてＺ_ｎｔｓ又はＺ_ｎｔの事後分布を算出し、さらに、当該Ｚ_ｎｔｓ又はＺ_ｎｔの事後分布に基づいてθ_ｔの第１の更新事後分布の更新値を算出するステップ、あるいは、（ｂ）所与されたＺ_ｎｔｓ又はＺ_ｎｔの初期分布に基づいてθ_ｔの第１の事後分布の更新値を算出するステップ、
（ｃ） 直前のステップ（ｄ）（ただし、初回はステップ（ａ）又はステップ（ｂ）である）により算出されたθ_ｔに基づいて、新たにＺ_ｎｔｓ又はＺ_ｎｔの事後分布を算出するステップ、
（ｄ） ステップ（ｃ）において算出された Ｚ_ｎｔｓ又はＺ_ｎｔの事後分布を基にして、θ_ｔの事後分布を新たに算出して更新するステップ、
（ｅ） ステップ（ｄ）において得られたθ_ｔの事後分布の期待値の収束性を評価するステップであって、当該期待値における収束が認められれば、当該収束期待値をθ_ｔの推定値として決定し、収束が認められなければ、ステップ（ｃ）、（ｄ）、及び、（ｅ）の繰り返しを決定する。

上記の（ａ）ステップは、繰り返し最適化を行う本法の初期値を定義するステップである。具体的には、各特定遺伝子座群又は個別の遺伝子座のアリルｔについて、ｑ^＊（θ）のパラメータであるα^＊ _ｔについて、各々の特定遺伝子座群又は個別の遺伝子座のアリル頻度は等しいと仮定して、α^＊ _ｔ＝（Ｎ／Ｔ＋ α_０） ／ Σ_ｔ（Ｎ／Ｔ ＋ α_０）を設定することが好ましい。

（ｂ）ステップは、ｑ^＊（θ）の現在の推定値が所与のもと、Ｅ_ｚ［Ｚ_ｎｔｓ］又はＥ_ｚ［Ｚ_ｎｔ］を計算するステップである。Ｅ_ｚ［Ｚ_ｎｔｓ］の値は、ｎについては１からＮの各々の数字であり、ｔについては１からＴについての各々の数字であり、ｓについては１≦ｓ≦ｌ_ｔ−Ｌ＋１、ｌ_ｔは、特定遺伝子座群又は個別の遺伝子座のアリルｔの長さ、Ｌはリード長である。Ｅ_ｚ［Ｚ_ｎｔｓ］は、ｑ^＊（θ）の現在の推定値に基づいて、

ここで、シングルエンドリードの場合、

（ｌｏｇρ_ｎｔｓを算出するための右辺第１項の変数は上記式（１）（ａ）の対数に期待値を取ったもの、第２項は上記式（１）（ｃ）に、及び、第３項は上記式（２）について、対数計算することにより算出することができる）

また、ペアエンドリードの場合、

（ｌｏｇρ_ｎｔｓを算出するための右辺中括弧内第１項は上記式（３）（ａ）の対数に期待値を取ったもの、第２項は上記式（３）（ａ’）に、第３項は上記式（３）（ｂ）に、第４項は上記式（４−１）に、及び、第５項は上記式（４−２）について、対数計算することにより算出することができる）

一方Ｚ_ｎｔは、上記Ｚ_ｎｔｓを各ｔについて可能なｓについて全てを考慮したものであるため、

として計算される。

以下、Ｚ_ｎｔｓの代わりに、Ｚ_ｎｔを用いて一連の計算を行い、以降のステップを実行することも可能である。

ここで、

（式中、ψ(・)はディガンマ関数である）
である。

（ｃ）ステップは、ｑ^＊（Ｚ）の現在の推定値が所与のもと、Ｅ_θ［θ_ｔ］を計算するステップである。

Ｅ_θ［θ_ｔ］は、ｑ^＊（Ｚ）の現在の推定値に基づいて、

として計算される。

従って、ｑ^＊（θ）もまたディリクレ分布であり、事前分布ｐ（θ）は共役事前分布である。

（ｄ）ステップにおいて、収束基準の好適な一例として、存在量パラメータの期待値Ｅ_θ［θ_ｔ］＞１０^−７である特定遺伝子座のアリル頻度についての、相対的変化１０^−３が収束基準として使用されるが、異なる収束基準を用いることも可能である。

ここで変分の下限についての検討を行い、本発明の変分ベイズ法の更新により得られる収束値が、真の目的変数θを近似するものであることを示す。

本発明の変分ベイズ法における対数周辺尤度は、

（式中、

として分解することができる。

KL(q||p)は、q(θ, Z)とp(θ, Z |R)の間のカルバック・ライブラー距離であるので、対数周辺尤度は、L(q)によって下限を与えられる。すなわち、カルバック・ライブラー距離は常に０以上であるため、L(q)が当該周辺尤度の下限を構成する。上記の（ｂ）ステップと（ｃ）ステップを繰り返し更新する度に、L(q)（対数周辺尤度の下限）を増加させることが一般的に示される（Bishop CM. Pattern Recognition and Machine Learning. Springer Science:Business Media, LLC, New York, NY, USA; 2006）。

上記の因子分解仮定

を用いて、

（式中、

である。

なお、通常は（ｅ）ステップの「ステップ（ｃ）、（ｄ）、及び、（ｅ）の繰り返し」が、１回以上は行われるが、この繰り返し工程を１回も行わない、又は収束する以前でステップを終了する場合も想定される。

（３）ＥＭアルゴリズムを用いるコンピュータシステムとコンピュータプログラム

（ｉ）上述したＥＭアルゴリズムを用いるコンピュータシステムとして、例えば、下記のコンピュータシステムが挙げられる。この本発明のコンピュータシステムの態様は、上記のＥＭアルゴリズムを用いた推定方法をコンピュータにおいて実行することを特徴とするコンピュータシステムである。

当該コンピュータシステムは、特定遺伝子座群又は個別遺伝子座の対応リード情報が混在したデータについて、個々のリードにおける個々の当該遺伝子座群のアリルに対する期待マッピング数を最適化するコンピュータシステムであって、記録部と演算処理部を具え、下記の処理（Ａ）〜（Ｅ）の全部又は一部；
（Ａ） 当該記録部には、被験者由来のＤＮＡのリード情報が、リードの配列及びリードのマッピング先である当該遺伝子座群又は個別の遺伝子座のアリルの観測データとして記録されており、
（Ｂ） 当該演算処理部では読み出された前記観測データに基づき、下記の初期化処理（Ｂ）−１及び（Ｂ）−２のいずれかが実行され、
（Ｂ）−１：当該遺伝子座群又は個別の遺伝子座のアリル頻度に関する変数θの初期値の算出処理、
（Ｂ）−２：上記変数θ及び、観測データである被験者のＤＮＡのリード情報が混在したデータにおけるリードの塩基配列を媒介する２種の潜在変数としての下記（ａ）及び（ｂ）：
（ａ）リードｎの当該遺伝子座群又は個別の遺伝子座のアリル選択に関する、θに依存する変数Ｔ_ｎ、
（ｂ）リードｎの開始位置に関する、Ｔ_ｎに依存するＳ_ｎ、
が要約された、指標変数Ｚ_ｎｔｓ（Ｚ_ｎｔｓは、（Ｔ_ｎ，Ｓ_ｎ）＝（ｔ，ｓ）の場合１であり、それ以外は０である。）＝１、又は、潜在変数Ｔ_ｎが要約された、Ｚ_ｎｔ（Ｚ_ｎｔは、Ｔ_ｎ＝ｔの場合１であり、それ以外は０である。）＝１、の事後確率の初期値の算出処理、
（Ｃ） 当該演算処理部において、上記処理（Ｂ）−１で算出された変数θに基づき、当該指標変数Ｚ_ｎｔｓ又はＺ_ｎｔ＝１の事後確率の算出処理がなされ、
（Ｄ） 当該演算処理部において、上記処理（Ｂ）−２、又は、処理（Ｃ）で算出された当該指標変数Ｚ_ｎｔｓ又はＺ_ｎｔ＝１の事後確率に基づいて変数θの最尤推定値の第１の更新値が算出され、
（Ｅ） 当該演算処理部において、上記処理（Ｄ）で算出された変数θの最尤推定値の第１の更新値に基づいて上記処理（Ｃ）と処理（Ｄ）が再度実行され、さらに、変数θの第２の更新値が算出されるループ処理が、新たな更新値と前回の更新値との間における差異が実質的に認められなくなるまで繰り返し実行されて、収束した変数θが最適化されたθとして、上記記録部に記録がなされる；
処理が実行されることを特徴とするコンピュータシステムである。

（ii）上述したＥＭアルゴリズムを用いるコンピュータプログラムとして、例えば、下記のコンピュータプログラムが挙げられる。この本発明のコンピュータプログラムの態様は、上記のＥＭアルゴリズムを用いた推定方法をコンピュータにおいて実現することを特徴とするコンピュータプログラムである。

当該コンピュータプログラムは、特定遺伝子座群又は個別遺伝子座の対応リード情報が混在したデータについて、個々のリードにおける個々の当該遺伝子座群又は個別の遺伝子座のアリルに対する期待マッピング数を最適化するコンピュータプログラムであって、コンピュータに下記の第１の機能〜第５の機能の全て又は一部；
（Ａ） 被験者由来のＤＮＡのリード情報が、リードの配列及びリードのマッピング先である当該遺伝子座群又は個別の遺伝子座のアリルの観測データとして記録されている記録部から当該データを読み出す、第１の機能、
（Ｂ） 上記第１の機能によって読み出した前記観測データに基づき、下記の初期化処理（Ｂ）−１及び（Ｂ）−２のいずれかを実行する、第２の機能、
（Ｂ）−１：当該遺伝子座群又は個別の遺伝子座のアリル頻度に関する変数θの初期値の算出処理、
（Ｂ）−２：上記変数θ及び、観測データである被験者のＤＮＡのリード情報が混在したデータにおけるリードの塩基配列を媒介する２種の潜在変数としての下記（ａ）及び（ｂ）：
（ａ）リードｎの当該遺伝子座群又は個別の遺伝子座のアリル選択に関する、θに依存する変数Ｔ_ｎ、
（ｂ）リードｎの開始位置に関する、Ｔ_ｎに依存するＳ_ｎ、
が要約された、指標変数Ｚ_ｎｔｓ（Ｚ_ｎｔｓは、（Ｔ_ｎ，Ｓ_ｎ）＝（ｔ，ｓ）の場合１であり、それ以外は０である。）＝１、又は、潜在変数Ｔ_ｎが要約された、Ｚ_ｎｔ（Ｚ_ｎｔは、Ｔ_ｎ＝ｔの場合１であり、それ以外は０である。）＝１、の事後確率の初期値の算出処理、
（Ｃ） 上記第２の機能の（Ｂ）−１で算出した変数θに基づき、当該指標変数Ｚ_ｎｔｓ又はＺ_ｎｔの算出処理を行う、第３の機能、
（Ｄ） 上記第２の機能の（Ｂ）−２、又は、第３の機能により算出した当該指標変数Ｚ_ｎｔｓ＝１又はＺ_ｎｔ＝１基づいて変数θの最尤推定値の第１の更新値を算出する、第４の機能、
（Ｅ） 上記第４の機能で算出した変数θの最尤推定値の第１の更新値に基づいて上記第３の機能と第４の機能を再度実行し、さらに、変数θの第２の更新値を算出するループ処理を、新たな更新値と前回の更新値との間における差異が実質的に認められなくなるまで繰り返し実行し、収束した変数θが最適化されたθとして、上記記録部に記録を行う、第５の機能；
を実現させるアルゴリズムが含まれることを特徴とするコンピュータプログラムである。

（４）変分ベイズ法を用いるコンピュータシステムとコンピュータプログラム

（ｉ）上述した変分ベイズ法を用いるコンピュータシステムとして、例えば、下記のコンピュータシステムが挙げられる。この本発明のコンピュータシステムの態様は、上記の変分ベイズ法を用いた推定方法をコンピュータにおいて実行することを特徴とするコンピュータシステムである。

当該コンピュータシステムは、特定遺伝子座群又は個別遺伝子座の対応リード情報が混在したデータについて、個々のリードにおける個々のアリルに対する期待マッピング数を最適化するコンピュータシステムであって、記録部と演算処理部を具え、下記処理（Ａ）〜（Ｅ）の全部又は一部；
（Ａ） 当該記録部には、被験者由来のＤＮＡのリード情報が、リードの配列及びリードのマッピング先である当該遺伝子座群又は個別の遺伝子座のアリルの観測データとして記録されており、
（Ｂ） 当該演算処理部では読み出された前記観測データに基づき、下記の初期化処理（Ｂ）−１及び（Ｂ）−２のいずれかが実行され、
（Ｂ）−１：当該遺伝子座群又は個別の遺伝子座のアリル頻度についての予備知識の分布を示すハイパーパラメータα_ｔの初期値に基づくθ_ｔの事後分布の更新値の算出処理、
（Ｂ）−２：上記θの分布、及び、観測データである被験者由来のＤＮＡのリード情報が混在したデータにおけるリードの塩基配列を媒介する２種の潜在変数としての下記（ａ）及び（ｂ）：
（ａ）リードｎの当該遺伝子座群又は個別の遺伝子座のアリル選択に関する、θに依存する変数Ｔ_ｎ、
（ｂ）リードｎの開始位置に関する、Ｔ_ｎに依存するＳ_ｎ、
が要約された、指標変数Ｚ_ｎｔｓ（Ｚ_ｎｔｓは、（Ｔ_ｎ，Ｓ_ｎ）＝（ｔ，ｓ）の場合１であり、それ以外は０である。）、又は、潜在変数Ｔ_ｎが要約された、Ｚ_ｎｔ（Ｚ_ｎｔは、Ｔ_ｎ＝ｔの場合１であり、それ以外は０である。）の事後分布の初期分布の算出処理、
（Ｃ） 当該演算処理部において、上記処理（Ｂ）−１で算出された変数θの分布に基づき、当該指標変数Ｚ_ｎｔｓ又はＺ_ｎｔの事後分布の算出処理がなされ、
（Ｄ） 当該演算処理部において、上記処理（Ｂ）−２、又は、処理（Ｃ）で算出された当該指標変数Ｚ_ｎｔｓ又はＺ_ｎｔの事後分布の更新値に基づいて変数θの第１更新事後分布の更新値が算出され、
（Ｅ） 当該演算処理部において、上記処理（Ｄ）で算出された変数θの第１の更新事後分布に基づいて上記処理（Ｃ）と処理（Ｄ）が再度実行され、さらに、変数θの第２の更新事後分布が算出されるループ処理が、新たに更新された事後分布の期待値と前回に更新された事後分布の期待値との間における差異が実質的に認められなくなるまで繰り返し実行されて、収束したθの期待値が最適化されたθのデータとして、上記記録部に記録がなされる；
が実行されることを特徴とするコンピュータシステムである。

（ii）上述した変分ベイズ法を用いるコンピュータプログラムとして、例えば、下記のコンピュータプログラムが挙げられる。この本発明のコンピュータプログラムの態様は、上記の変分ベイズ法を用いた推定方法をコンピュータにおいて実現することを特徴とするコンピュータプログラムである。

当該コンピュータプログラムは、特定遺伝子座群又は個別遺伝子座の対応リード情報が混在したデータについて、個々のリードにおける個々の当該遺伝子座群又は個別の遺伝子座のアリルに対する期待マッピング数を最適化するコンピュータプログラムであって、コンピュータに下記の第１の機能〜第５の機能の全て又は一部；
（Ａ） 被験者由来のＤＮＡのリード情報が、リードの配列及びリードのマッピング先である当該遺伝子座群又は個別の遺伝子座のアリルの観測データとして記録されている記録部から当該データを読み出す、第１の機能、
（Ｂ） 上記第１の機能によって読み出した前記観測データに基づき、下記の初期化処理（Ｂ）−１及び（Ｂ）−２のいずれかを実行する、第２の機能、
（Ｂ）−１：当該遺伝子座群又は個別の遺伝子座のアリル頻度についての予備知識の分布を示すハイパーパラメータα_ｔの初期値に基づくθ_ｔの事後分布の更新値の算出処理、
（Ｂ）−２：上記θの分布、及び、観測データである被験者由来のＤＮＡのリード情報が混在したデータにおけるリードの塩基配列を媒介する２種の潜在変数としての下記（ａ）及び（ｂ）：
（ａ）リードｎの当該遺伝子座のアリル選択に関する、θに依存する変数Ｔ_ｎ、
（ｂ）リードｎの開始位置に関する、Ｔ_ｎに依存するＳ_ｎ、
が要約された、指標変数Ｚ_ｎｔｓ（Ｚ_ｎｔｓは、（Ｔ_ｎ，Ｓ_ｎ）＝（ｔ，ｓ）の場合１であり、それ以外は０である。）、又は、潜在変数Ｔ_ｎが要約された、Ｚ_ｎｔ（Ｚ_ｎｔは、Ｔ_ｎ＝ｔの場合１であり、それ以外は０である。）、の事後分布の初期分布の算出処理、
（Ｃ） 上記第２の機能の（Ｂ）−１で算出した変数θの分布に基づき、当該指標変数Ｚ_ｎｔｓ又はＺ_ｎｔの事後分布の算出処理を行う、第３の機能、
（Ｄ） 上記第２の機能の（Ｂ）−２、又は、第３の機能で算出した当該指標変数Ｚ_ｎｔｓ又はＺ_ｎｔの事後分布の更新値に基づいて変数θの第１更新事後分布の更新値を算出する、第４の機能、
（Ｅ） 上記第４の機能で算出した変数θの第１の更新事後分布に基づいて上記第３の機能と第４の機能を再度実行し、さらに、変数θの第２の更新事後分布を算出するループ処理を、新たに更新された事後分布の期待値と前回に更新された事後分布の期待値との間における差異が実質的に認められなくなるまで繰り返し実行して、収束したθの期待値を最適化されたθのデータとして、上記記録部に記録を行う第５の機能；
を実現させるアルゴリズムが含まれることを特徴とするコンピュータプログラムである。

［Ｂ］ 本発明の方法のさらに具体的な形態
ここでは、特定遺伝子座群をヒトＭＨＣであるＨＬＡ遺伝子座群として用いた例を開示するが、ここに開示したアルゴリズムは、他の遺伝子座、例えば、シトクロムＰ４５０遺伝子座群、免疫グロブリンをコードする遺伝子座群、T細胞受容体をコードする遺伝子座群、嗅覚受容体をコードする遺伝子座群を対象にしてもよい。上記したように、本発明の方法は最適化方法と判定方法の総称として用いる。

図１は、被験者のＨＬＡ遺伝子座の遺伝型を予測するための本発明の方法の処理の流れである。この処理の流れの各々の要素は、全てコンピュータにおいて、必要に応じてコンピュータネットワークやデータベースを介して、コンピュータソフトウエアのアルゴリズムを実現する命令に基づき電子的に行われるものである。

ボックスＢ１−１は、リードデータの存在を示しており、ボックスＢ１−２は、ヒトゲノム参照配列の存在を示している。上述した通り「リードデータ」（ボックスＢ１−１）は、被験者の遺伝子検体を高性能シークエンサで処理することにより得られる、個々のＤＮＡ断片（リード）の一部又は全部の塩基配列の電子情報であり、通常は読み取り精度を示す情報も付加されている（ＦＡＳＴＱ：ＤＮＡの塩基配列を表すＦＡＳＴＡフォーマットに由来する用語）。また、この段階で「全リード深度」が確定している。「ヒトゲノム参照配列」（ボックスＢ１−２）は、上述した通りに一個人もしくは複数のゲノム配列情報に帰着するものであり、その提供源は特に限定されない。被験者と当該ゲノム配列情報の人種が異なっており、ＨＬＡ遺伝子配列に関して異質性を含んでいることも想定される。参照配列のバージョン変更に本発明の実質的な効果が依存するものではなく、常にこれに対応させて被験者のＨＬＡ型を決定することが可能であり、現在又は将来提供される最新の参照配列を用いることができる。ヒトゲノム参照配列として、例えば、ＵＣＳＣ(University of California Santa Cruz)が管理しているヒトゲノム参照配列（ｈｇ１９等）、Genome Reference Consortiumが管理しているヒトゲノム参照配列（ＧＲＣｈ３７等）等が挙げられるが、これらに限定されるものではない。また、上述したように参照配列として、例えば、ターゲットシークエンス、Ｅｘоｍｅシークエンス、ＲＮＡシークエンス、ＰａｃＢｉｏ ＲＳ II、Ｏｘｆｏｒｄ Ｎａｎｏｐｏｒｅ等のロングリードシークエンスデータ等も用いることができる。

ステップＳ１は、１回目のマッピングを行うステップであり、上記ボックスＢ１−１の「リードデータ」を、ボックスＢ１−２の「ヒトゲノム参照配列」に対してマッピングを行うステップである。このマッピングに関しては、当業者であればコンピュータで当該マッピングを実現可能なアルゴリズムを含むコンピュータプログラムを作出して行うことも可能であるが、既に提供されているソフトウエアを用いることも可能である。既存のマッピング用のソフトウエアとして、例えば、ＢＷＡ−ＭＥＭ(http://bio-bwa.source.net/)、Bowtie2(http://bowtie-bio.sourceforge.net/bowtie2/index.shtml)、Novoalign (http://www.novocraft.com/products/novoalign/) 等が挙げられる。が挙げられる。また、上記リードデータがペアエンドシークエンスデータの場合は、リードの両端の塩基配列（ペア）の一方がＨＬＡ遺伝子座にマッピングされており、他方がマッピングされていない場合は、当該ペアの両方のリードを抽出して以下の工程に用いることが好適である。

なお、上記のヒトゲノム参照配列には「おとり配列」が含まれることが好適である。「おとり配列」とは、予め用意された参照配列には存在しないゲノム配列である。おとり配列を用いない場合、ヒトゲノム参照配列に登録されている配列由来ではないリードが、既存のヒトゲノム参照配列のどこかにマッピングされてしまう危険性が高くなり、マッピングの精度が低下する可能性があるからである。おとり配列としては、ｈｓ３７ｄ５（ftp: //ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/phase2_reference_assembly_sequence/hs37d5.fa.gz）等が例示される。

ボックスＢ２は、上記Ｓ１における「１回目のマッピングの結果」の存在を示すボックスであり、マッピングの結果は、ＳＡＭ形式又はＢＡＭ形式等のフォーマットで存在している。

ステップＳ２−１は、上記ボックスＢ２の「１回目のマッピングの結果」から、「ＨＬＡ遺伝子座群にマッピングされたリードの抽出」を行うステップであり、ステップＳ２−２は、ヒトゲノム参照配列のどの箇所にもマッピングがなされなかった「非マッピングリードの抽出」を行うステップである。これらの抽出の実行に関しては、当業者であればコンピュータで当該抽出を実現可能なアルゴリズムを含むコンピュータプログラムを作出して行うことも可能であるが、既に提供されているソフトウエアを用いることも可能である。既存の抽出用のソフトウエアとして、例えば、ＳＡＭtools(http://samtools.sourceforge.net/)が挙げられる。この抽出ステップにおいては必要な処理な処理が、少なくとも全てのＨＬＡ遺伝子座にわたるので数が多く、可能な限り効率的な演算を行うことが好適である。よって、全てのＨＬＡ遺伝子座（現状では、ＨＬＡ−Ａ、−Ｂ、−Ｃ、−ＤＭ、−ＤＯ、−ＤＰ、−ＤＱ、−ＤＲ、−Ｅ、−Ｆ、−Ｇ、−Ｈ、−Ｊ、−Ｋ、−Ｌ、−Ｐ、−Ｖ、−ＭＩＣ、及び−ＴＡＰがＩＭＧＴ／ＨＬＡデータベースに登録されている）に関する抽出が行われるように、上記コンピュータプログラムが設計されていることが好適である。上記ＳＡＭtoolsにおいては、容易に抽出処理が可能である。ステップＳ２−２の「非マッピングリードの抽出」は、上述した通りに参照ゲノム配列と被験者のゲノム配列の根源的なズレによってマッピングされなかったＨＬＡ遺伝子座群由来のリードの喪失を防ぐためにステップとして設定されていることが好ましい。

これらのステップＳ２により、実質的にＨＬＡ遺伝子座以外にマッピングされたリードの情報が除かれる。

ボックスＢ３は、上記ステップＳ２により抽出されたリード情報（抽出済みリード）の存在を示すボックスであり、ボックスＢ４は、「ＨＬＡアリルの参照配列」の情報の存在を示すボックスである。ＨＬＡアリルの参照配列情報は、ＨＬＡアリルのゲノム情報が格納されているデータベースに由来する電子情報である。当該データベースとしては、例えば、ＩＭＧＴ／ＨＬＡ(http://www.ebi.ac.uk/ipd/imgt/hla/)、等が挙げられる。これらのデータベース情報は、可能な限り最新の更新情報として取得することが好ましい。また、当該ＨＬＡアリルのゲノム情報には、実際にはタンパク質をコードしない遺伝子である「偽遺伝子」のＨＬＡアリルも含めることが好適である。

ステップＳ３は、前記ボックスＢ３の「抽出済みリード情報」の、ボックスＢ４の「ＨＬＡアリルの参照配列」に対する「２回目のマッピング」が実行されるステップである。この２回目のマッピングの実行形式は、一つのリードが複数のＨＬＡアリルに対してマッピングされることを許容するように行われるようにマッピング用のソフトウエアが設計されていることが好適であり、上記のＢＷＡ−ＭＥＭにおいては「−ａオプション」、Ｂｏｗｔｉｅ２においては「−ｋオプション」を指定して実行がなされることにより、この複数のＨＬＡアリルに対するマッピングをコンピュータにおいて実現させることが可能である。

ボックスＢ５は、上記ステップＳ３における「２回目のマッピングの結果」、すなわち、「ＨＬＡ対応リード情報」の存在を示すボックスであり、当該情報はＳＡＭ形式又はＢＡＭ形式等のフォーマットで存在している。

図１のうち、ここまでが上述した「本発明のマッピングプロセス」を示している。

以下、ステップＳ４は上記「ＨＬＡ対応リード情報」を最適化するステップであり、その詳細をフローシートとして図２（図２−１、図２−２、図２−３）に示す。図２−１は、潜在変数Ｚ_ｎｔを用いた本発明の最適化方法の一態様（変分ベイズ法とＥＭアルゴリズム）を示すものであり、図２−２は、潜在変数Ｚ_ｎｔｓを用いた本発明の最適化方法におけるＥＭアルゴリズムを用いた一態様を示すものであり、図２−３は、潜在変数Ｚ_ｎｔｓを用いた本発明の最適化方法における変分ベイズ法を用いた一態様を示すものである。これらの態様に伴う詳細な計算工程は、特定遺伝子座群のアリルに対して前述した通りである。さらに、ボックスＢ６は所望する「ＨＬＡ型の決定」がなされた電子情報の存在を示すボックスであるが、上記ステップＳ４からボックスＢ６に至る過程において、図３（後述）において示した棄却閾値を用いたＨＬＡタイピング工程が好適に行われる。

［潜在変数Ｚ_ｎｔを用いた最適化法の態様］
図２−１は、上記のように「ＨＬＡ対応リード情報」の潜在変数Ｚ_ｎｔを用いた最適化工程の一態様を示すフローシートである。このフローシートは、直接的にはＨＬＡアリル上の期待されるリードカウントについて、ベイズ推定を行うために、変分ベイズ法による予測を行う態様であるが、下記のハイパーパラメータα_０の繰り入れを行わない、最尤推定を行うためにＥＭアルゴリズムによる予測を行うことも可能である。

ステップＳ４−１１は、上記ボックスＢ４に示される「ＨＬＡアリルの参照配列」の読み出し、及び、ボックスＢ５に示される「ＨＬＡ対応リード情報」の読み出しを行うステップである。

ステップＳ４−１２は、本コンピュータプログラムの実行による最適化処理に際して用いられるリードｎの各ＨＬＡアリルへの期待マッピング数Ｅ［Ｚ_ｎｔ］と、リード総量に対して各ＨＬＡアリルに割り当てられた合計期待マッピング数の割合の期待値Ｅ［θ_ｔ］の初期化を行うためのステップである。

ステップＳ４−１３は、上記ステップＳ４−１２により読み出した初期化値に基づき、全てのリードにおける個々のＨＬＡアリルに対する合計期待マッピング数の数値化処理を行うステップである。ここに記載された数式の意義は、前述した通りである。ｒ_ｔが「ＨＬＡアリルｔに割り当てられた数値化された合計期待マッピング数」である。

ステップＳ４−１４は、上記ステップＳ４−１３によって数値化された「合計期待マッピング数」を、それぞれ全てのＨＬＡアリルにおける合計期待マッピング数の和で除して、リード総量に対して各ＨＬＡアリルに割り当てられた合計期待マッピング数の割合の期待値Ｅ［θ_ｔ］を算出する処理を行うステップである。本ステップＳ４−１４において、上述のようにハイパーパラメータα_０が繰り入れられており、ここでは変分ベイズ法の態様が記述されている。当該ハイパーパラメータα_０は、例えば、事前知識が何も無いと仮定して、一様分布を与えるα_０＝１等が用いられるが、対数尤度の下限を最大化するα_０、例えば、α_０＝０．０１を用いることも好ましい。さらに前述したように、このハイパーパラメータの繰り入れを行わずに、パラメータθ_ｔの最尤推定をおこなうＥＭアルゴリズム態様のフローシートとして用いることも可能である。

ステップＳ４−１５は、ループ処理に係わるステップであり、
（ｉ）上記ステップＳ４−１４により算出した合計期待マッピング数の割合の分布を、頻度分布として個々のＨＬＡアリルに対して割り当て、当該割り当て頻度分布を前提にして、再び上記Ｓ４−１３で改めて個々のリードにおける個々のＨＬＡアリル毎の期待マッピング数を算出する機能、
（ii）上記（ｉ）により算出した新たなリード毎の個々のＨＬＡアリル毎の期待マッピング数に対して、再び上記ステップＳ４−１３又はＳ４−１４を実行して、ＨＬＡアリルにマッピングされているリード総量に対して各ＨＬＡアリルに割り当てられたリードの割合の分布を新たに算出する機能、
（iii）上記（ｉ）と（ii）における算出処理を、今回の（ｉ）の実行により算出するリード毎の個々のＨＬＡアリルに対する期待マッピング数と、前回の（ｉ）の実行により算出する当該期待マッピング数との間における差が全てのリードについて認められなくなるか、又は、今回の（ii）の実行により算出する合計期待マッピング数の割合の分布の期待値と、前回の（ii）の実行により算出する当該割合の分布の期待値との差が全てのＨＬＡアリルについて認められなくなるまで、繰り返し実行し、収束したリード毎の個々のＨＬＡアリルに対する期待マッピング数、又は、収束したＨＬＡアリル毎のリードの割合の分布の期待値を、「収束値」として最適化されたデータとして認定する機能；
が、このステップＳ４−１５の記述に含まれる。

［潜在変数Ｚ _ｎｔｓ を用いた最適化法の態様（１）］
図２−２は、上記のように「ＨＬＡ対応リード情報」の潜在変数Ｚ_ｎｔｓを用いた最適化工程の一態様を示すフローシートである。このフローシートは、ＨＬＡアリル上の期待されるリードカウントについて、最尤推定を行うためにＥＭアルゴリズムによる予測を行うためのものである。

ステップＳ４−２１は上記ステップＳ４−１１と同様に、「ＨＬＡ対応リード情報」の読み出しを行うステップである。

ステップＳ４−２２は、本コンピュータプログラムの実行による最適化処理に際して用いられる潜在変数Ｚ_ｎｔｓ（リードｎがＨＬＡアリルｔの位置ｓから生成されている場合は１を取る潜在変数）と、θ_ｔ（リード総量に対して各ＨＬＡアリルに割り当てられたリードの割合）の初期化を行うためのステップである。潜在変数Ｚ_ｎｔｓの期待値は、Ｅ［Ｚ_ｎｔｓ］と記述されている。ここでＮは「リード総数」であり、Ｔは「ＨＬＡアリルの種類の総数」である。ｌ_ｔは、ＨＬＡアリルｔの長さであり、Ｌはリード長である。Ｍ_ｎは、観測されたリードｎのマッピング先の総数である。Ｅ［Ｚ_ｎｔｓ］は、事前の無情報・一様分布を前提として、１／Ｍ_ｎとして初期設定がなされている。θ_ｔは、事前の無情報・一様分布を前提として、１／Ｔが初期値として与えられている。

ステップＳ４−２３は、上記ステップＳ４−２２によって読み出した初期化情報に基づき、全てのリードにおける個々のＨＬＡアリルに対する期待マッピング数の数値化処理を行うステップである。すなわち、直前に更新された（初回は上記ステップＳ４−２２の初期化値）θ_tに基づいて、リードｎのＨＬＡアリルｔの位置ｓへの割り当て分を示す潜在変数Ｚ_ｎｔｓの期待値（Ｅ_ｚ［Ｚ_ｎｔｓ］）、すなわち、「期待マッピング数」を再計算するステップであり、ＥＭ法のＥステップに相当する。

ステップＳ４−２４は、上記ステップＳ４−２３によって数値化された「期待マッピング数」を、ＨＬＡアリル毎に合算して「合計期待マッピング数」（ｒ_ｔ＝Σ_{ｎ’，ｔ’＝ｔ，ｓ’} Ｅ_ｚ［Ｚ_{ｎ’ｔ’ｓ’}］）を算出し、さらに当該合計期待マッピング数を、それぞれ全てのＨＬＡアリルにおける合計期待マッピング数の和で除して、リード総量に対して各ＨＬＡアリルに割り当てられたリードの割合「θ_ｔ」の最尤推定値を算出する処理（θ_ｔ＝ｒ_ｔ／Σ_ｔ’ｒ_ｔ’）を行うステップであり、ＥＭ法のＭステップに相当するステップである。このステップは、対数尤度を局所最大化するθ_ｔを算出することに対応している。

ステップＳ４−２５は、ループ処理に係わるステップである。すなわち、上記ステップＳ４−２４で算出したθ_ｔの最尤推定値の第１の更新値に基づいて上記ステップＳ４−２３とステップＳ４−２４を再度実行し、さらに、変数θの第２の更新値を算出するループ処理を、新たな更新値と前回の更新値との間における差異が全てのＨＬＡアリルについて実質的に認められなくなるまで繰り返し実行し、収束した変数θが最適化されたθとして、上記記録部に記録を行う。

［潜在変数Ｚ _ｎｔｓ を用いた最適化法の態様（２）］
図２−３は、上記のように「ＨＬＡ対応リード情報」の潜在変数Ｚ_ｎｔｓを用いた最適化工程の一態様を示すフローシートである。このフローシートは、ＨＬＡアリル上の期待されるリードカウント、及びＨＬＡアリル頻度について、変分ベイズ法による予測を行うためのものである。

ステップＳ４−３１は上記ステップＳ４−１１と同様に、「ＨＬＡ対応リード情報」の読み出しを行うステップである。

ステップＳ４−３２は、本コンピュータプログラムの実行による最適化処理に際して用いられる潜在変数Ｚ_ｎｔｓ（リードｎがＨＬＡアリルｔの位置ｓから生成された場合は１、そうではない場合は０を取る指標変数）の事後分布と、θ_ｔ（リード総量に対してＨＬＡアリルに割り当てられたリードの割合）の事後分布の初期化を行うためのステップである。潜在変数Ｚ_ｎｔｓの事後分布の期待値は、Ｅ［Ｚ_ｎｔｓ］と記述されている。ここでＮは「リード総数」であり、Ｔは「ＨＬＡアリルの種類の総数」である。ｌ_ｔは、ＨＬＡアリルの長さであり、Ｌはリード長である。Ｍ_ｎは、観測されたリードｎのマッピング先の総数である。Ｅ［Ｚ_ｎｔｓ］は、事前の無情報・一様分布を前提として、１／Ｍ_ｎとして初期設定がなされている。θ_ｔの事後分布の期待値であるＥ_θ［θ_ｔ］は、ディリクレ分布を事前分布として、α^＊ｔ／Σ_ｔ’α^＊ｔ’として与えられる。ここでα^＊ _ｔ＝α_０＋Ｎ／Ｔである。前述したように、α_０はハイパーパラメータであり、例えば、事前知識が何も無いと仮定して、一様分布を与えるα_０＝１等が用いられるが、対数尤度の下限を最大化するα_０、例えばα_０＝０．０１を用いることも好ましい。

ステップＳ４−３３は、上記ステップＳ４−３２によって読み出した初期化情報に基づき、全てのリードにおける個々のＨＬＡアリルに対する期待マッピング数の数値化処理を行うステップである。すなわち、直前に更新された（初回は上記ステップＳ４−３２の初期化値）θ_tの事後分布に基づいて、リードｎのＨＬＡアリルｔの位置ｓへの割り当て分を示す潜在変数Ｚ_ｎｔｓの事後分布の期待値（Ｅ_ｚ［Ｚ_ｎｔｓ］）を再計算するステップであり、変分ベイズ法のＥステップに相当する。

ステップＳ４−３４は、上記ステップＳ４−３３によって数値化された「期待マッピング数」（Ｅ_ｚ［Ｚ_ｎｔｓ］）によりθ_ｔの事後分布を算出する、変分ベイズ法のＭステップに相当するステップであり、θ_ｔの事後分布の期待値Ｅ_θ［θ_ｔ］は、直前のＳ４−３３で算出されたＥ_ｚ［Ｚ_ｎｔｓ］によって、

として計算される。

ステップＳ４−３５は、上記ステップＳ４−３３とＳ４−３４を繰り返すループ処理を行うか、当該ループ処理を終了するかの選択を行うステップである。

すなわち、上記ステップＳ４−２４で算出した変数θ_ｔの第１の更新事後分布に基づいて上記ステップＳ４−３３とＳ４−３４を再度実行し、さらに、変数θの第２の更新事後分布を算出するループ処理を、新たに更新された事後分布の期待値と前回に更新された事後分布の期待値との間における差異が実質的に認められなくなるまで繰り返し実行して、収束したθの期待値を最適化されたθのデータとして、上記記録部に記録を行う。

このようにして得られた「ｒ_ｔ」又は「θ_ｔ」は、そのままＨＬＡタイピングの指標とすることも可能であるが、好適には後述のＨＬＡタイピング工程が施されることが好適である。

図３は、上記の最適化工程を経て得られた「各ＨＬＡアリルに割り当てられたリードの割合『θ_ｔ』、もしくはリードの割合『θ_ｔ』の期待値」を基に、全リード深度を基準とする個別深度に対する棄却閾値を用いた当該ＨＬＡタイピング工程の一例を示したフローシートである。

ステップＳ５−１は、上記の「各ＨＬＡアリルに割り当てられたリードの割合」から、各ＨＬＡアリルにおける「個別深度」を算出する処理を行うステップである。変数「ｔ」と「ｎ」の初期化を行い、Ｓ５−１ボックス内に示された数式を用いて当該個別深度「ｄ_ｔ」の算出処理が行われる。ここの算出過程の詳細については既に述べた。

ステップＳ５−２は、ボックス内に記されているように、「各ＨＬＡ遺伝子座について、ＨＬＡアリルのうち、個別深度「ｄ_ｔ」が大きなものから２つを選択し、最も大きい個別深度のアリルと、２番目に大きい個別深度のアリルを、それぞれ指定する処理を行うステップである。

ステップＳ５−３は、上記の最も大きいアリルの個別深度「ｄ_{ｆｉｒｓｔ}」が、棄却深度「Ｄ」よりも小さいか否かの選択処理を行うステップである。棄却深度「Ｄ」は、上述したように「全リード深度」の５〜５０％、好適には１０〜３０％の間で選択される数値である。もしも、個別深度「ｄ_{ｆｉｒｓｔ}」が棄却深度「Ｄ」よりも小さい場合には、「ＨＬＡ型は決定されない」（結論Ｄ５−１）との判断がなされ、小さくない場合には、次の選択ステップＳ５−４に移行する判断がなされる。

ステップＳ５−４は、上記の２番目に大きいアリルの個別深度「ｄ_{ｓｅｃｏｎｄ}」が、棄却深度「Ｄ」よりも小さいか否かの選択処理を行うステップである。もしも、個別深度「ｄ_{ｓｅｃｏｎｄ}」が棄却深度「Ｄ」よりも小さい場合には、選択ステップＳ５−５に移行する判断がなされ、小さくない場合には、選択ステップＳ５−６に移行する判断がなされる。

ステップＳ５−５は、個別深度「ｄ_{ｓｅｃｏｎｄ}」が棄却深度「Ｄ」よりも小さい場合に適用される選択ステップである。すなわち、上記「ｄ_{ｆｉｒｓｔ}」が個別深度「Ｄ」の２倍より大きい場合には、「ＨＬＡ型は、個別深度が最大のＨＬＡアリルのホモ接合である」（結論Ｄ５−２）と決定する判断処理がなされ、２倍よりも大きくない場合には、「ＨＬＡ型は、個別深度が最大のＨＬＡアリルのヘテロ接合と決定し、ヘテロのもう一方のアリルは決定されない」（結論Ｄ５−３）とする判断処理がなされる。

ステップＳ５−６は、ステップＳ５−４において「ｄ_{ｓｅｃｏｎｄ}＜Ｄ」ではない場合に適用される、さらなる選択ステップであって、上記「ｄ_{ｆｉｒｓｔ}」が個別深度「Ｄ」の２倍より大きい場合には、「ＨＬＡ型は、個別深度が最大のＨＬＡアリルのホモ接合である」（結論Ｄ５−４）と決定する判断処理がなされ、２倍よりも大きくない場合には、「ＨＬＡ型は、個別深度が最大と２番目のＨＬＡアリルのヘテロ接合である」（結論Ｄ５−５）とする判断処理がなされる。

以上のようにして本発明のシステムは、好適な形でのＨＬＡ型の決定（図１のボックスＢ６）を行うことができる。

以下に、対象をヒトＭＨＣであるＨＬＡとして行った本発明の実施例を開示する。

［Ａ］ 実施例において用いたソース
本実施例において用いたコンピュータシステムとコンピュータプログラム（以下、「本発明のシステム」と総称する）の実行に際して用いられたソースは、以下に示す通りである。

（１）本発明のＨＬＡマッピングプロセス
図１のボックスＢ１−１の「リード情報」を提供する次世代シークエンサは、HiSeq2000（Illumina社）を用いた。リード情報は、ＦＡＳＴＡＱフォーマットで提供され、引き続くリード情報も同様にＦＡＳＴＡＱフォーマットである。

ボックスＢ１−２の「ヒトゲノム参照配列」としては、ｈｇ１９（ＵＣＳＣ）又はＧＲＣｈ３７(Genome Reference Consortium)を、おとり配列（ｈｓ３７ｄ５）と併せて用いた。

ステップＳ１／Ｓ３のマッピングは、ＢＷＡ−ＭＥＭにより実行された。ステップS３については、オプション “-a” が指定され、各リードについて全ての可能なマッピング先が出力されるように実行された。

ボックスＢ２の「１回目のマッピング結果」、及び、ボックスＢ５の「２回目のマッピング結果」は、ＢＡＭ形式で用いられた。

ステップＳ２のマッピング結果の抽出用のソフトウエアとして、「ＳＡＭtools」を用いた。

ボックスＢ４の「ＨＬＡアリルの参照配列」は、ＩＭＧＴ／ＨＬＡデータベースからＦＡＳＴＡフォーマットで提供された。

（２）最適化プロセス
ステップＳ４の最適化プロセスは、図２−１の変分ベイズ法を用いたフローシートのアルゴリズムをペアエンドデータに適用することにより実行され、棄却深度は全リード深度の２０％と設定がなされて、図３のフローシートのアルゴリズムにより棄却プロセスが実行された。

［Ｂ］ＨＬＡタイピングの性能測定
（１）シミュレーション試験のあらまし
本発明のシステムの予測性能を、予測精度の観点で評価した。予測精度は、真のＨＬＡ型の中の真の陽性予測の分数として定義される。このシミュレーション実験において、６つのＨＬＡ遺伝子座（ＨＬＡ−Ａ、−Ｂ、−Ｃ、−ＤＱＡ１、−ＤＱＢ１、−ＤＲＢ１）に関する２つのＨＬＡアリル（ヘテロ接合又はホモ接合いずれか）が各個体で評価された。予測性能は、各方法について、２ケタ、４ケタ、６ケタ及び８ケタ解像度で別々に評価された。

（２）シミュレーションデータ分析
シミュレーションデータ分析を用いて、本発明のシステムがＨＬＡ型を予測する性能を、他のシステムと比べて評価した。比較システムとして（ａ）ＰＨＬＡＴ（Bai et al., BMC genomics, 15:325 (2014)）、及び、（ｂ）ＨＬＡminer（Warren et al., Genome medicine, 4(12):102 (2013)）を用いた。これらの比較システムは、ＨＬＡクラスＩ遺伝子座（ＨＬＡ−Ａ、−Ｂ、−Ｃ）、及び、クラスＩＩ遺伝子座（ＨＬＡ−ＤＱＡ１、−ＤＱＢ１及び−ＤＲＢ１）を、全ヒトゲノム配列データから４ケタ解像度で分類することができることが知られている。

まず、１０００個体分のヒト試料のシミュレーションデータを準備した。６種類のＨＬＡ遺伝子座のためのそのＨＬＡディプロタイプは、ＩＭＧＴ／ＨＬＡデータベースリリース３．１５．０に登録されたＨＬＡアリルからランダムに選択された。ＨＬＡ型が各個体について確定されると、１０００bpについて１つのＳＮＰがヒトゲノムにおける平均塩基多様性に基づいて各個体のＨＬＡアリルにおいて組み入れられた。次いで、１００bpのペアエンドのリードデータ（平均深度が５×、１０×、２０×及び３０×であって、その断片長分布の平均偏差及び標準偏差は、それぞれ３００bp及び４０bpとして設定された）が、０．１％の置換、削除及び挿入エラーを持って、ＨＬＡアリル配列の全域で均一に生成された。

表１は、本発明のシステム及び上記既存ツールの予測精度を、３０×シミュレーションデータ分析において示したものである。

表１において、本発明のシステムは既存のＰＨＬＡＴやＨＬＡminerでは実施出来ない８ケタ解像度でＨＬＡ型を９９．９４％という特に高い精度で推定することが出来た。既存のＰＨＬＡＴやＨＬＡminerと推定精度を比較した場合、４ケタ解像度及び６ケタ解像度おいてはこれらの既存のソフトウエアよりも良好な解像性能を示すことが明らかになった。図４は、シミュレーションデータの範囲の種々の深度で、４ケタ解像度でのＨＬＡ型の予測精度を示した図面である。図４において、本発明のシステムを用いた場合の予測精度は、検討を行った深度全体にわたって、ＰＨＬＡＴ及びＨＬＡminerでの精度よりも常に良好であった。

特に本発明のシステムでは、平均深度５×シミュレーションデータから、４ケタ解像度で９９．３６％の精度でＨＬＡ型の予測がなされた。ＰＨＬＡＴは、尤度を算するためにＳＮＰ部位だけの検討が行われるが、その限られた情報のみでは、他の多型部位（例えば欠失又は挿入）がＨＬＡ型を決定するために重要な場合には有効とはいえない。ＰＨＬＡＴのもう一つのあり得る欠点は、当該システムがＨＬＡアリル頻度について事前情報を必要とすることである。しかし、ＨＬＡアリル頻度はヒト集団の中でも特に多様であり、提供された遺伝子検体の人種的なルーツを推測することは、常に可能であるとは限らない。それに対して本発明のシステムではＨＬＡアリル頻度に関する事前情報は必要としない。

（３）現実データの分析（１）
本発明のシステムを、国際ＨａｐＭａｐプロジェクトで使用されたＣＥＵトリオ試料（ＮＡ１２８７８（子供）、ＮＡ１２８９１（父親）、及び、ＮＡ１２８９２（母親））の、ＰＣＲ法による増幅がなされていない全ゲノムシークエンスデータに適用した。１００bpのペアエンドデータを、HiSeq2000を用いて導出し、その平均挿入長は３００bpであり、当該範囲における深度は各試料について４５×であった（データ全てはIllumina社により提供された）。表２は、上記ＣＥＵトリオにおいて予測されたＨＬＡ型を示している。

表２より、本発明のシステムを用いることにより予測されたクラスＩ遺伝子座（ＨＬＡ−Ａ、−Ｂ、及び、−Ｃ）に関するＨＬＡ型は、４ケタ解像度で実験的に検証されたＨＬＡ型と一致した。これらを表２中の太文字テキストで示した。ＨＬＡ型の多くが、本発明のシステムにより、８ケタ解像度で予測された。本発明のシステムで予測されたＨＬＡ−Ａ、−Ｂ及び−Ｃ遺伝子座のＨＬＡ型は、「Ｂ＊０７：０２：０１」（ＮＡ１２８９１中の一つのアリル）を除いて、６ケタ解像度でＰＨＬＡＴにより予測されたものと一致していた。他の文献（Major et al., PloS one, 8(11):e78410 (2013)）もまた、「Ｂ＊０７：０２：０１」としてＮＡ１２８９１の当該ＨＬＡ−Ｂアリルの一つを報告している。これに対してＰＨＬＡＴは、当該ＨＬＡ型を、「Ｂ＊０７：０２：２９」として予測した。全体的に本発明のシステムで予測されたトリオ（子供、父親及び母親）のＨＬＡ型は、ＰＨＬＡＴで予測されたものよりも一致していた。

本発明のシステムで予測されたＨＬＡ−ＤＱＡ１、−ＤＱＢ１及び−ＤＲＢ１遺伝子座のＨＬＡ型は、ＤＱＡ１＊０１：０１：０２（ＮＡ１２８７８中の一つのアリル及びＮＡ１２８９２中の２つのアリル）を除いて、６ケタ解像度でＰＨＬＡＴにおいて予測されたものと一致していた。ＰＨＬＡＴは、本発明のシステムではＨＬＡ型「ＤＱＡ１＊０１：０１：０２」と予測したＨＬＡアリルを「ＤＱＡ１＊０１：０１：０１」として予測した。しかしながらそのゲノム配列自体が、ＩＭＧＴデータベースリリース３．１５．０において欠落しており、従って、両システム間の当否を判断することはできなかった。

（４）現実データの分析（２）
本発明のシステムを、１ＫＪＰＮ集団（東北メディカル・メガバンク計画のコホート調査に参加した健常な日本人1,070人）に適用し、ＨＬＡ−Ａ、ＨＬＡ−Ｂ、及び、ＨＬＡ−Ｃ遺伝子座のＨＬＡアリルを推定した。本例は、ＩＭＧＴ／ＨＬＡデータベースに登録されたゲノムＨＬＡアリル配列に対する配列リードのマッピングに基づいている。これにより、２１４０アリルの中から２０６３アリルについてのＨＬＡ−Ａアリルを、フル解像度（ＨＬＡ命名規則における８桁）で、分類が可能であることを確認した。

また、この際に予測された４桁解像度（アミノ酸配列を変更するヌクレオチドの相違）での、ＨＬＡ−Ａ、ＨＬＡ−Ｂ、及び、ＨＬＡ−Ｃのアリル頻度について、ＰＣＲ−ＳＳＯＰ(Itoh,Y. et al.,Immunogenetics 57, 717-729(2005))を用いて別の日本人１０１８人の集団において決定された４桁解像度の頻度に非常に似ていることを確認した。比較された二つの日本人集団はどちらも１０００人以上の十分な数のサンプルを含んでいるため、実際の日本人集団のＨＬＡアリル頻度に近い分布をしていると考えられる。本発明のシステムの推定結果とＰＣＲ−ＳＳＯＰの推定結果が非常に似ているという結果は、どちらの手法を使用しても、日本人集団のアリル頻度を４桁解像度で正しく推定できたということを示唆している。

図５−１はＨＬＡ−Ａについて、図５−２はＨＬＡ−Ｂについて、図５−３はＨＬＡ−Ｃについての両解析により算出されたＨＬＡアリル頻度を検討した結果を示すグラフである。それぞれ、縦軸はアリル頻度、横軸は４桁解像度のＨＬＡアリル型を示し、同じＨＬＡアリル型の左側のグラフバーは今回の１ＫＪＰＮ集団において本発明のシステムを用いて推定した解析結果であり、右側のグラフバーは別の日本人１０１８人の集団においてＰＣＲ−ＳＳＯＰが適用された解析結果を示している。

本発明のシステムは、ヒトＭＨＣであるＨＬＡ遺伝子座群やシトクロムＰ４５０遺伝子座群、免疫グロブリンをコードする遺伝子座群、T細胞受容体をコードする遺伝子座群、嗅覚受容体をコードする遺伝子座群等のゲノム中に類似の塩基配列を持つ座位が複数存在する、あるいは遺伝的多型が多数知られている遺伝子座群について、かならずしも遺伝子座特異的なプライマーデザイン又は当該遺伝子座群のアリル頻度の予備的知識の必要なく、全ゲノムシークエンスデータを用いる有効かつ正確な遺伝型のタイピングを実現する。個別、あるいは集団規模のシークエンスデータに関わらず、選択された任意の遺伝子座群での各遺伝子座の遺伝型のタイピングの際に、本発明のシステムを容易に適用することが可能であり、それは、遺伝型と表現型の関連を特定するための研究に、また、臓器移植の際のドナーとレシピエントのＨＬＡ型マッチングのような臨床業務に、有用である。また、ＨＬＡ以外のＭＨＣ、例えば、マウスＭＨＣであるＨ−２(histocompatibility-２)等の哺乳類のＭＨＣや、ニワトリのＢ遺伝子座等の鳥類のＭＨＣについて、本発明を適用して遺伝型を推定することにより、より的確かつ詳細な品種の鑑別、さらにペットや保護動物の疾病の治療指針を立てる基礎知見の提供手段として、本発明を用いることが可能である。

Ｂ１−１・・・リードデータの存在を示すボックス
Ｂ１−２・・・ヒトゲノム参照配列の存在を示すボックス
Ｓ１・・・１回目のマッピングを行うステップ
Ｂ２・・・１回目のマッピングの結果の存在を示すステップ
Ｓ２−１・・・１回目のマッピングの結果からリードの抽出を行うステップ
Ｓ２−２・・・非マッピングリードの抽出を行うステップ
Ｂ３・・・Ｓ２により抽出されたリード情報の存在を示すボックス
Ｂ４・・・ＨＬＡアリルの参照配列の存在を示すボックス
Ｓ３・・・２回目のマッピングを行うステップ
Ｂ５・・・ＨＬＡ対応リード情報の存在を示すボックス
Ｂ６・・・ＨＬＡ型の決定がなされた電子情報の存在を示すボックス
Ｓ４−１１・・・最適化を行うための読み出し機能を記述するステップ
Ｓ４−１２・・・最適化処理を行うためのパラメータの初期化を行うステップ
Ｓ４−１３・・・最適化を行うためのＥステップの機能を記述するステップ
Ｓ４−１４・・・最適化を行うためのＭステップの機能を記述するステップ
Ｓ４−１５・・・最適化を行うためのループ・収束機能を記述するステップ
Ｓ４−２１・・・最適化を行うための読み出し機能を記述するステップ
Ｓ４−２２・・・最適化処理を行うためのパラメータの初期化を行うステップ
Ｓ４−２３・・・最適化を行うためのＥステップの機能を記述するステップ
Ｓ４−２４・・・最適化を行うためのＭステップの機能を記述するステップ
Ｓ４−２５・・・最適化を行うためのループ・収束機能を記述するステップ
Ｓ４−３１・・・最適化を行うための読み出し機能を記述するステップ
Ｓ４−３２・・・最適化処理を行うためのパラメータの初期化を行うステップ
Ｓ４−３３・・・最適化を行うためのＥステップの機能を記述するステップ
Ｓ４−３４・・・最適化を行うためのＭステップの機能を記述するステップ
Ｓ４−３５・・・最適化を行うためのループ・収束機能を記述するステップ
Ｓ５−１・・・個別深度を算出する処理を行うステップ
Ｓ５−２・・・個別深度の大きな２つを選択する処理を行うステップ
Ｓ５−３・・・棄却深度と最大の個別深度の大小によって選択処理を行うステップ
Ｄ５−１・・・ＨＬＡ型が決定されない結論を示すボックス
Ｓ５−４・・・棄却深度と２番目の個別深度の大小によって選択処理を行うステップ
Ｓ５−５・・・２番目の個別深度が棄却深度よりも小さい場合に行われるステップ
Ｄ５−２・・・最大の個別深度のＨＬＡアリルのホモ接合との結論を示すボックス
Ｄ５−３・・・最大の個別深度のＨＬＡアリルのヘテロ接合であり、もう一方のアリルは決定しないとの結論を示すボックス
Ｓ５−６・・・２番目の個別深度が棄却深度よりも小さくない場合に行われるステップ
Ｄ５−４・・・最大の個別深度のＨＬＡアリルのホモ接合との結論を示すボックス
Ｄ５−５・・・最大と２番目の個別深度のＨＬＡアリルのヘテロ接合との結論を示すボックス

Claims (35)

1. ゲノム中に類似の塩基配列を持つ座位が複数存在する、若しくは遺伝的多型が多数知られている、選択された遺伝子座群又は個別の遺伝子座のアリル由来のＤＮＡのリード情報が混在したデータのリードの塩基配列に対してマッピングを行うことにより得られる、当該遺伝子座群又は個別の遺伝子座のアリルに対する各リードのマッピング対応が特定されたリード情報に対して、下記のステップ（１）〜（６）の全部又は一部が実行されることを特徴とする、遺伝子のリード情報の最適化方法。
   （１） 個々のリードにおける個々の当該遺伝子座群又は個別の遺伝子座のアリルに対する期待マッピング数の数値化が行われるステップ；
   （２） ステップ（１）において数値化された期待マッピング数が当該遺伝子座群又は個別の遺伝子座のアリル毎に合算されて合計期待マッピング数が算出されるステップ；
   （３） ステップ（２）において算出された合計期待マッピング数が、それぞれ個々の当該遺伝子座群又は個別の遺伝子座のアリルにおける合計期待マッピング数の和で除されて、当該遺伝子座群又は個別の遺伝子座のアリルにマッピングされているリード総量に対して当該遺伝子座群又は個別の遺伝子座の各アリルに割り当てられたリードの割合が算出されるステップ；
   （４） ステップ（３）において得られたリードの割合が、頻度として個々の当該遺伝子座群又は個別の遺伝子座のアリルに対して割り当てられ、当該割り当て頻度を前提にして、再びステップ（１）により改めて得られた個々のリードにおける個々の当該遺伝子座群又は個別の遺伝子座のアリル毎の期待マッピング数が算出されるステップ；
   （５） ステップ（４）において得られた新たな期待マッピング数に対して、再びステップ（２）又は（３）が実行され、当該遺伝子座群又は個別の遺伝子座のアリルにマッピングされているリード総量に対して当該遺伝子座群又は個別の遺伝子座の各アリルに割り当てられたリードの割合が新たに算出されるステップ；
   （６） ステップ（４）と（５）が、ステップ（４）において算出されるリード毎の個々の当該遺伝子座群又は個別の遺伝子座のアリルに対する期待マッピング数と、前回のステップ（４）で算出される当該期待マッピング数との間における差が全てのリードについて認められなくなるか、又は、ステップ（５）において算出されるリードの割合の値と、前回のステップ（５）で算出される当該割合の値との間における差が全ての当該遺伝子座群又は個別の遺伝子座のアリルについて認められなくなるまで、繰り返し実行され、収束したリード毎の個々の当該遺伝子座群又は個別の遺伝子座のアリルに対する期待マッピング数、あるいは、収束した当該遺伝子座群又は個別の遺伝子座のアリル毎のリードの割合の値が、最適化されたデータとして認定されるステップ。
2. 選択された遺伝子座群又は個別の遺伝子座のアリル由来のＤＮＡのリード情報が混在したデータにおけるリード全体の塩基配列を観測データＲとして、個々のリードにおける個々の当該遺伝子座群又は個別の遺伝子座のアリルに対する期待マッピング数を求めるステップ、並びに、目的パラメータである当該遺伝子座群又は個別の遺伝子座のアリル頻度θ（θはＴ次元ベクトル、Ｔは当該遺伝子座群又は個別の遺伝子座のアリルの種類数）の推定値を求めるステップ、を含む被験者由来のＤＮＡのリード情報の当該遺伝子座群又は個別の遺伝子座のアリルへのマッピングを、コンピュータにより最適化する最適化方法において、
   上記目的パラメータθ及び、観測データＲを媒介する潜在変数である、（ａ）リードｎの当該遺伝子座群又は個別の遺伝子座のアリル選択に関するθに依存する変数Ｔ_ｎ、及び、（ｂ）リードｎの開始位置に関するＴ_ｎに依存するＳ_ｎについて、
   リードｎの塩基配列を観測データＲ_ｎとして、観測データＲ_ｎからの目的パラメータθの推測工程において観測データＲ_ｎが依存するように、少なくとも（ｉ）変数Ｔ_ｎ及びＳ_ｎ、あるいは、(ii) 変数Ｔ_ｎ、を繰り入れて当該推定値を算出することを特徴とする、最適化方法。
3. 最尤推定法、又は、ベイズ推定法に基づくステップの実行により、個々のリードにおける個々の選択された遺伝子座群又は個別の遺伝子座のアリルに対する期待マッピング数、目的パラメータである当該遺伝子座群又は個別の遺伝子座のアリル頻度θの推定値を算出することを特徴とする、請求項１又は２に記載の最適化方法。
4. 潜在変数Ｔ_ｎ及びＳ_ｎが要約された、指標変数Ｚ_ｎｔｓ（Ｚ_ｎｔｓは、（Ｔ_ｎ，Ｓ_ｎ）＝（ｔ，ｓ）の場合１であり、それ以外は０である。）、又は、潜在変数Ｔ_ｎが要約された、Ｚ_ｎｔ（Ｚ_ｎｔは、Ｔ_ｎ＝ｔの場合１であり、それ以外は０である。）、を潜在変数として用いることを特徴とする、請求項２又は３に記載の最適化方法。
5. 下記（１）〜（５）のステップを行うことを特徴とする、請求項４に記載の最適化方法。
   （１）所与されたθ_ｔの初期値に基づいて、Ｚ_ｎｔｓ＝１又はＺ_ｎｔ＝１の事後確率の第１の更新値を算出し、さらに、当該Ｚ_ｎｔｓ＝１又はＺ_ｎｔ＝１の事後確率の第１の更新値に基づいてθ_ｔの最尤推定値の第１の更新値を算出するステップ、あるいは、（２）所与されたＺ_ｎｔｓ＝１又はＺ_ｎｔ＝１の事後確率の初期値に基づいて、θ_ｔの最尤推定値の第１の更新値を算出するステップ、
   （３） 直前のステップ（４）（ただし、初回はステップ（１）又はステップ（２）である）により算出されたθ_ｔの最尤推定値の更新値に基づいて、新たにＺ_ｎｔｓ＝１又はＺ_ｎｔ＝１の事後確率の更新値を算出するステップ、
   （４） ステップ（３）において算出されたＺ_ｎｔｓ＝１又はＺ_ｎｔ＝１の事後確率の更新値に基づいて、θ_ｔの最尤推定値の更新値を新たに算出するステップ、
   （５） （ｉ）ステップ（３）において算出されたＺ_ｎｔｓ＝１又はＺ_ｎｔ＝１、及び、ステップ（４）において算出されたθ_ｔに基づいて対数尤度を計算して、対数尤度の収束性を評価するステップ、
   （ii）ステップ（３）で算出されたＺ_ｎｔｓ＝１又はＺ_ｎｔ＝１の事後確率の更新値の収束性を評価するステップ、あるいは、
   （iii）ステップ（４）で算出されたθ_ｔの最尤推定値の更新値の収束性を評価するステップであって、
   収束が認められれば、それぞれのステップにおけるθ_ｔを最終推定値として決定し、収束が認められなければ、ステップ（３）、（４）、及び、（５）の繰り返しを決定する。
6. 下記（１）〜（５）のステップを行うことを特徴とする、請求項４に記載の最適化方法。
   （１） 所与された選択された遺伝子座群又は個別の遺伝子座のアリルｔのアリル頻度についての予備知識の分布を示すハイパーパラメータα_ｔの初期値に基づくθ_ｔの事後分布の更新値に基づいてＺ_ｎｔｓ又はＺ_ｎｔの事後分布を算出し、さらに、当該Ｚ_ｎｔｓ又はＺ_ｎｔの事後分布に基づいてθ_ｔの第１の更新事後分布の更新値を算出するステップ、あるいは、（２）所与されたＺ_ｎｔｓ又はＺ_ｎｔの初期分布に基づいてθ_ｔの第１の事後分布の更新値を算出するステップ、
   （３） 直前のステップ（４）（ただし、初回はステップ（１）又はステップ（２）である）により算出されたθ_ｔに基づいて、新たにＺ_ｎｔｓ又はＺ_ｎｔの事後分布を算出するステップ、
   （４） ステップ（３）において算出された Ｚ_ｎｔｓ又はＺ_ｎｔの事後分布を基にして、θ_ｔの事後分布を新たに算出して更新するステップ、
   （５） ステップ（４）において得られたθ_ｔの事後分布の期待値の収束性を評価するステップであって、当該期待値における収束が認められれば、当該収束期待値をθ_ｔの推定値として決定し、収束が認められなければ、ステップ（３）、（４）、及び、（５）の繰り返しを決定する。
7. 選択された遺伝子座群又は個別の遺伝子座のアリル由来のＤＮＡのリード情報が混在したデータは、被験者のリード情報をデータベースに登録されている当該遺伝子座群又は個別の遺伝子座のアリルの塩基配列に対してマッピングをすることにより得られる、当該遺伝子座群又は個別の遺伝子座のアリルに対する各リードのマッピング対応が特定されたリード情報であって、当該マッピングは、下記のステップ（ａ）及び（ｂ）により実行されることを特徴とする、請求項１〜６のいずれか１項に記載の最適化方法。
   （ａ） 被験者のリードの塩基配列情報において、ヒト遺伝子の塩基配列に対するマッピングが行われ、当該遺伝子座群又は個別の遺伝子座のアリルにマッピングされたリードが抽出されるステップ；
   （ｂ） ステップ（ａ）により抽出された当該遺伝子座群又は個別の遺伝子座のアリルにマッピングされたリードの配列情報に対して、データベースに登録されている当該遺伝子座群又は個別の遺伝子座のアリルの塩基配列とのマッピングが行われ、マッピングされたリードが当該遺伝子座群又は個別の遺伝子座のアリル毎に抽出され、当該遺伝子座群又は個別の遺伝子座のアリルに対する各リードのマッピング対応が特定されたリード情報が得られるステップ。
8. ステップ（ａ）及び（ｂ）において実行されるマッピングは、一つのリードが複数の選択された遺伝子座群又は個別の遺伝子座のアリルに対してマッピングされることを許容することを特徴とする、請求項７に記載の最適化方法。
9. ステップ（ａ）の選択された遺伝子座群又は個別の遺伝子座のアリルにマッピングされたリードに加えて、ヒト遺伝子に対してマッピングがなされなかったリードが併せて抽出され、これが（ｂ）ステップの再マッピングの対象とされることを特徴とする、請求項７又は８に記載の最適化方法。
10. 選択された遺伝子座群又は個別の遺伝子座は、ＭＨＣの遺伝子座群又は個別の遺伝子座であることを特徴とする、請求項１〜９のいずれか１項に記載の最適化方法。
11. ＭＨＣはＨＬＡであることを特徴とする、請求項１０に記載の最適化方法。
12. 請求項１〜１１のいずれか１項に記載の最適化方法により得られた選択された遺伝子座群又は個別の遺伝子座のアリル頻度から当該遺伝子座群又は個別の遺伝子座のアリル毎のリードの個別深度が算出され、 当該遺伝子座群の各遺伝子座又は個別の遺伝子座について当該個別深度の大きなアリルから順に２個以内について選択され、当該遺伝子座群の各遺伝子座又は個別の遺伝子座の遺伝型の要素として決定がなされることを特徴とする、選択された遺伝子座群又は個別の遺伝子座の遺伝型の決定方法。
13. 請求項１〜１２のいずれか１項に記載の最適化方法により得られた選択された遺伝子座群又は個別の遺伝子座のアリル頻度から当該遺伝子座群又は個別の遺伝子座のアリル毎のリードの個別深度が算出され、当該遺伝子座群の各遺伝子座又は個別の遺伝子座について当該個別深度の大きなアリルから順に２個以内について選択され、当該遺伝子座群の各遺伝子座又は個別の遺伝子座の遺伝型の要素として決定がなされる遺伝型の判定方法において、全リード深度の５〜５０％のいずれかの頻度数が棄却閾値として設定され、当該閾値以下の個別深度の当該遺伝子座群又は個別の遺伝子座のアリルは遺伝型決定の要素から除外されることを特徴とする、選択された遺伝子座群の各遺伝子座又は個別の遺伝子座の遺伝型の決定方法。
14. 選択された遺伝子座群の各遺伝子座又は個別の遺伝子座の遺伝型決定の要素からの除外が行われた後、下記（ｉ）又は（ii）の決定がなされることを特徴とする、請求項１３に記載の遺伝型の決定方法。
    （ｉ） 当該遺伝子座群の各遺伝子座又は個別の遺伝子座の遺伝型決定の対象が１個の当該遺伝子座のアリルについては、当該１個のアリルの個別深度が前記棄却閾値の２倍以上の場合には、当該アリルはホモ接合と決定がなされ、又は、前記棄却閾値の２倍より小さい場合はヘテロ接合であると決定がなされる。
    （ii） 当該遺伝子座群の各遺伝子座又は個別の遺伝子座の遺伝型決定の対象が２個の当該遺伝子座のアリルについては、個別深度が大きな方が小さい方の２倍未満である場合には、両アリルはヘテロ接合であるとの決定がなされ、又は、個別深度が大きな方が小さい方の２倍以上である場合には、大きな方のアリルはホモ接合であるとの決定がなされる。
15. 選択された遺伝子座群又は個別の遺伝子座は、ＭＨＣの遺伝子座群又は個別の遺伝子座であることを特徴とする、請求項１２〜１４のいずれか１項に記載の遺伝型の決定方法。
16. ＭＨＣはＨＬＡであることを特徴とする、請求項１５に記載の遺伝型の決定方法。
17. 最適化対象とされた遺伝子座群又は個別の遺伝子座のアリル由来のＤＮＡのリード情報が混在したデータのリードの塩基配列に対してマッピングを行うことにより得られる、当該遺伝子座群又は個別の遺伝子座のアリルに対する各リードのマッピング対応が特定されたリード情報、を最適化するコンピュータシステムであって、記録部と演算処理部を備え、下記の処理（Ａ）〜（Ｇ）の全て又は一部；
    （Ａ） 当該記録部には、被験者由来のＤＮＡのリード情報が、リードの配列及びリードのマッピング先である当該遺伝子座群又は個別の遺伝子座のアリルのデータとして記録されており、
    （Ｂ） 当該演算処理部では、前記記録部の情報に基づいて、個々のリードにおける個々の当該遺伝子座群又は個別の遺伝子座のアリルに対する期待マッピング数の数値化処理が実行され、
    （Ｃ） 上記処理（Ｂ）において数値化された期待マッピング数が当該遺伝子座群又は個別の遺伝子座のアリル毎に合算されて合計期待マッピング数が算出され、
    （Ｄ） 上記処理（Ｃ）において算出された合計期待マッピング数が、それぞれ個々の当該遺伝子座群又は個別の遺伝子座のアリルにおける合計期待マッピング数の和で除されて、当該遺伝子座群又は個別の遺伝子座のアリルにマッピングされているリード総量に対して当該遺伝子座群又は個別の遺伝子座の各アリルに割り当てられたリードの割合が算出される処理が実行され、
    （Ｅ） 上記処理（Ｃ）において算出されたリードの割合が、頻度として個々の当該遺伝子座群又は個別の遺伝子座のアリルに対して割り当てられ、当該割り当て頻度を前提にして、再び上記処理（Ｂ）により改めて算出された個々のリードにおける個々の当該遺伝子座群又は個別の遺伝子座のアリル毎の期待マッピング数が算出される処理が実行され、
    （Ｆ） 上記処理（Ｅ）により算出された新たな期待マッピング数に対して、再び上記処理（Ｃ）又は（Ｄ）が実行されて、当該遺伝子座群又は個別の遺伝子座のアリルにマッピングされているリード総量に対して当該遺伝子座群又は個別の遺伝子座の各アリルに割り当てられたリードの割合が新たに算出される処理が実行され、
    （Ｇ） 上記処理（Ｅ）と（Ｆ）の処理が、処理（Ｅ）において算出されるリード毎の個々の当該遺伝子座群又は個別の遺伝子座のアリルに対する期待マッピング数と、前回の処理（Ｅ）において算出される当該期待マッピング数との間における差が全てのリードについて認められなくなるか、又は、処理（Ｆ）において算出されるリードの割合の値と、前回の処理（Ｆ）で算出される当該割合の値との差が全ての当該遺伝子座群又は個別の遺伝子座のアリルについて認められなくなるまで、繰り返し実行され、収束したリード毎の個々の当該遺伝子座群又は個別の遺伝子座のアリルに対する期待マッピング数、又は、収束した当該遺伝子座群又は個別の遺伝子座のアリル毎のリードの割合の値が、最適化されたデータとして認定が行われる；
    が実行されることを特徴とするコンピュータシステム。
18. 選択された遺伝子座群又は個別の遺伝子座のアリル由来のＤＮＡのリード情報が混在したデータについて、個々のリードにおける個々のアリルに対する期待マッピング数を最適化するコンピュータシステムであって、記録部と演算処理部を具え、下記の処理（Ａ）〜（Ｅ）の全部又は一部；
    （Ａ） 当該記録部には、被験者由来のＤＮＡのリード情報が、リードの配列及びリードのマッピング先である当該遺伝子座群又は個別の遺伝子座のアリルの観測データとして記録されており、
    （Ｂ） 当該演算処理部では読み出された前記観測データに基づき、下記の初期化処理（Ｂ）−１及び（Ｂ）−２のいずれかが実行され、
    （Ｂ）−１：当該遺伝子座群又は個別の遺伝子座のアリル頻度に関する変数θの初期値の算出処理、
    （Ｂ）−２：上記変数θ及び、観測データである被験者のＤＮＡのリード情報が混在したデータにおけるリードの塩基配列を媒介する２種の潜在変数としての下記（ａ）及び（ｂ）：
    （ａ）リードｎの当該遺伝子座群又は個別の遺伝子座のアリル選択に関する、θに依存する変数Ｔ_ｎ、
    （ｂ）リードｎの開始位置に関する、Ｔ_ｎに依存するＳ_ｎ、
    が要約された、指標変数Ｚ_ｎｔｓ（Ｚ_ｎｔｓは、（Ｔ_ｎ，Ｓ_ｎ）＝（ｔ，ｓ）の場合１であり、それ以外は０である。）、又は、潜在変数Ｔ_ｎが要約された、Ｚ_ｎｔ（Ｚ_ｎｔは、Ｔ_ｎ＝ｔの場合１であり、それ以外は０である。）＝１の事後確率の初期値の算出処理、
    （Ｃ） 当該演算処理部において、上記処理（Ｂ）−１で算出された変数θに基づき、当該指標変数Ｚ_ｎｔｓ又はＺ_ｎｔ＝１の事後確率の算出処理がなされ、
    （Ｄ） 当該演算処理部において、上記処理（Ｂ）−２、又は、処理（Ｃ）で算出された当該指標変数Ｚ_ｎｔｓ又はＺ_ｎｔ＝１の事後確率に基づいて変数θの最尤推定値の第１の更新値が算出され、
    （Ｅ） 当該演算処理部において、上記処理（Ｄ）で算出された変数θの最尤推定値の第１の更新値に基づいて上記処理（Ｃ）と処理（Ｄ）が再度実行され、さらに、変数θの第２の更新値が算出されるループ処理が、新たな更新値と前回の更新値との間における差異が実質的に認められなくなるまで繰り返し実行されて、収束した変数θが最適化されたθとして、上記記録部に記録がなされる；
    処理が実行されることを特徴とするコンピュータシステム。
19. 選択された遺伝子座群又は個別の遺伝子座のアリル由来のＤＮＡのリード情報が混在したデータについて、個々のリードにおける個々の当該遺伝子座群又は個別の遺伝子座のアリルに対する期待マッピング数を最適化するコンピュータシステムであって、記録部と演算処理部を具え、下記処理（Ａ）〜（Ｅ）の全部又は一部；
    （Ａ） 当該記録部には、被験者由来のＤＮＡのリード情報が、リードの配列及びリードのマッピング先である当該遺伝子座群又は個別の遺伝子座のアリルの観測データとして記録されており、
    （Ｂ） 当該演算処理部では読み出された前記観測データに基づき、下記の初期化処理（Ｂ）−１及び（Ｂ）−２のいずれかが実行され、
    （Ｂ）−１：当該遺伝子座群又は個別の遺伝子座のアリル頻度についての予備知識の分布を示すハイパーパラメータα_ｔの初期値に基づくθ_ｔの事後分布の更新値の算出処理、
    （Ｂ）−２：上記θの分布、及び、観測データである被験者由来のＤＮＡのリード情報が混在したデータにおけるリードの塩基配列を媒介する２種の潜在変数としての下記（ａ）及び（ｂ）：
    （ａ）リードｎの当該遺伝子座群又は個別の遺伝子座のアリル選択に関する、θに依存する変数Ｔ_ｎ、
    （ｂ）リードｎの開始位置に関する、Ｔ_ｎに依存するＳ_ｎ、
    が要約された、指標変数Ｚ_ｎｔｓ（Ｚ_ｎｔｓは、（Ｔ_ｎ，Ｓ_ｎ）＝（ｔ，ｓ）の場合１であり、それ以外は０である。）、又は、潜在変数Ｔ_ｎが要約された、Ｚ_ｎｔ（Ｚ_ｎｔは、Ｔ_ｎ＝ｔの場合１であり、それ以外は０である。）の分布の初期分布の算出処理、
    （Ｃ） 当該演算処理部において、上記処理（Ｂ）−１で算出された変数θの分布に基づき、当該指標変数Ｚ_ｎｔｓ又はＺ_ｎｔの事後分布の算出処理がなされ、
    （Ｄ） 当該演算処理部において、上記処理（Ｂ）−２、又は、処理（Ｃ）で算出された当該指標変数Ｚ_ｎｔｓ又はＺ_ｎｔの事後分布の更新値に基づいて変数θの第１更新事後分布の更新値が算出され、
    （Ｅ） 当該演算処理部において、上記処理（Ｄ）で算出された変数θの第１の更新事後分布に基づいて上記処理（Ｃ）と処理（Ｄ）が再度実行され、さらに、変数θの第２の更新事後分布が算出されるループ処理が、新たに更新された事後分布の期待値と前回に更新された事後分布の期待値との間における差異が実質的に認められなくなるまで繰り返し実行されて、収束したθの期待値が最適化された当該遺伝子座群又は個別の遺伝子座のアリル頻度のデータとして、上記記録部に記録がなされる；
    が実行されることを特徴とするコンピュータシステム。
20. 選択された遺伝子座群又は個別の遺伝子座のアリル由来のＤＮＡのリード情報が混在したデータは、被験者のリード情報をデータベースに登録されている当該遺伝子座群又は個別の遺伝子座のアリルの塩基配列に対してマッピングをすることにより得られる、当該遺伝子座群又は個別の遺伝子座のアリルに対する各リードのマッピング対応が特定されたリード情報であって、当該マッピングは、下記の処理（ａ）及び（ｂ）により実行されることを特徴とする、請求項１７〜１９のいずれか１項に記載のコンピュータシステム。
    （ａ） 被験者のリードの塩基配列情報において、ヒト遺伝子の塩基配列に対するマッピングが行われ、当該遺伝子座群又は個別の遺伝子座のアリルにマッピングされたリードが抽出される処理；
    （ｂ） 処理（ａ）により抽出された当該遺伝子座群又は個別の遺伝子座のアリルにマッピングされたリードの配列情報に対して、データベースに登録されている当該遺伝子座群又は個別の遺伝子座のアリルの塩基配列とのマッピングが行われ、当該遺伝子座群又は個別の遺伝子座のアリルに対する各リードのマッピング対応及びマッピング状態が特定されたリード情報が得られる処理。
21. 処理（ａ）及び（ｂ）において実行されるマッピングは、一つのリードが複数の選択された遺伝子座群又は個別の遺伝子座のアリルに対してマッピングされることを許容することを特徴とする、請求項２０に記載のコンピュータシステム。
22. 処理（ａ）の選択された遺伝子座群又は個別の遺伝子座のアリルにマッピングされたリードに加えて、ヒト遺伝子に対してマッピングがなされなかったリードが併せて抽出処理され、これが（ｂ）処理の再マッピングの対象とされることを特徴とする、請求項２０又は２１に記載のコンピュータシステム。
23. 選択された被験者の遺伝子座群又は個別の遺伝子座の遺伝型遺伝型の判定を行うコンピュータシステムであって、記録部と演算処理部を備え、下記（α）〜（δ）の処理の全部又は一部；
    （α） 当該記録部には、請求項１〜１１のいずれかに記載の最適化方法により得られた、被験者の当該遺伝子座群又は個別の遺伝子座のアリル頻度、及び、全リード深度、が少なくとも記録されており；
    （β） 当該演算処理部では、前記記録部の当該遺伝子座群又は個別の遺伝子座のアリル頻度を基とする、当該遺伝子座群又は個別の遺伝子座のアリル毎の個別深度への算出処理、及び、個々の当該遺伝子座群又は個別の遺伝子座のアリルに対する算出された当該個別深度の割り振り処理が実行され；
    （γ） 棄却閾値として設定されている、全リード深度の平均の５〜５０％のいずれかの頻度数に対して、当該数値以下の個別深度の当該遺伝子座群又は個別の遺伝子座のアリルは当該遺伝子座群又は個別の遺伝子座の遺伝型決定の要素から除外される処理が実行され；
    （δ）：
    （δ）−１ （γ）の除外処理の実行の後、当該遺伝子座群の各遺伝子座又は個別の遺伝子座の遺伝型決定の対象が１個の当該遺伝子座のアリルについては、当該１個のアリルの個別深度が前記棄却閾値の２倍以上である場合には、当該アリルはホモ接合と決定がなされる処理が実行され、又は、前記棄却閾値の２倍より小さい場合はヘテロ接合であると決定がなされる処理が実行され；
    （δ）−２ （γ）の除外処理の実行の後、当該遺伝子座群の各遺伝子座又は個別の遺伝子座の遺伝型決定の対象が２個の当該遺伝子座のアリルについては、個別深度が大きな方が小さい方の２倍未満である場合には、両アリルはヘテロ接合であるとの決定がなされる処理が実行され、又は、個別深度が大きな方が小さい方の２倍以上である場合には、大きな方のアリルはホモ接合であるとの決定がなされる処理が実行される；
    が実行されることを特徴とするコンピュータシステム。
24. 選択された遺伝子座群又は個別の遺伝子座は、ＭＨＣの遺伝子座群又は個別の遺伝子座であることを特徴とする、請求項１７〜２３のいずれか１項に記載のコンピュータシステム。
25. ＭＨＣはＨＬＡであることを特徴とする、請求項２４に記載のコンピュータシステム。
26. 選択された遺伝子座群又は個別の遺伝子座のアリル由来のＤＮＡのリード情報が混在したデータのリードの塩基配列に対してマッピングを行うことにより得られる、群のアリルに対する各リードのマッピング対応が特定されたリード情報、を最適化するコンピュータプログラムであって、コンピュータに下記の第１の機能〜第７の機能の全て又は一部；
    （Ａ） 被験者由来のＤＮＡのリード情報が、リードの配列及びリードのマッピング先である当該遺伝子座群又は個別の遺伝子座のアリルのデータとして記録されている記録部から、当該リード情報を読み出す、第１の機能、
    （Ｂ） 上記第１の機能により読み出したリード情報に基づいて、個々のリードにおける個々の当該遺伝子座群又は個別の遺伝子座のアリルに対する期待マッピング数の数値化処理を実行する、第２の機能、
    （Ｃ） 上記第２の機能により数値化した期待マッピング数が、当該遺伝子座群又は個別の遺伝子座のアリル毎に合算されて合計期待マッピング数を算出する、第３の機能、
    （Ｄ） 上記第３の機能により算出した合計期待マッピング数を、それぞれ全ての当該遺伝子座群又は個別の遺伝子座のアリルにおける合計期待マッピング数の和で除して、当該遺伝子座群又は個別の遺伝子座のアリルにマッピングされているリード総量に対して当該遺伝子座群又は個別の遺伝子座の各アリルに割り当てられたリードの割合を算出する、第４の機能、
    （Ｅ） 上記第４の機能により算出したリードの割合を、頻度として個々の当該遺伝子座群又は個別の遺伝子座のアリルに対して割り当て、当該割り当て頻度を前提にして、再び第２の機能で改めて算出した、個々のリードにおける個々の当該遺伝子座群又は個別の遺伝子座のアリル毎の期待マッピング数を算出する、第５の機能、
    （Ｆ） 上記第５の機能により算出した新たな期待マッピング数に対して、再び上記第３の機能又は第４の機能を実行して、当該遺伝子座群又は個別の遺伝子座のアリルにマッピングされているリード総量に対して当該遺伝子座群又は個別の遺伝子座の各アリルに割り当てられたリードの割合を新たに算出する、第６の機能、
    （Ｇ） 上記第５の機能と第６の機能を、第５の機能の実行により算出するリード毎の個々の当該遺伝子座群又は個別の遺伝子座のアリルに対する期待マッピング数と、前回の第５の機能の実行により算出する当該期待マッピング数との間における差が全てのリードについて認められなくなるか、あるいは、上記第６の機能の実行により算出するリードの割合の値と、前回の第６の機能の実行により算出する当該割合の値との差が当該遺伝子座群又は個別の遺伝子座の全てのアリルについて認められなくなるまで、繰り返し実行し、収束したリード毎の個々の当該遺伝子座群又は個別の遺伝子座のアリルに対する期待マッピング数、又は、収束した当該遺伝子座群又は個別の遺伝子座のアリル毎のリードの割合の値を、最適化されたデータとして認定する、第７の機能；
    を実現させるアルゴリズムが含まれることを特徴とする、コンピュータプログラム。
27. 選択された遺伝子座群又は個別の遺伝子座のアリル由来のＤＮＡのリード情報が混在したデータについて、個々のリードにおける個々の当該遺伝子座群又は個別の遺伝子座のアリルに対する期待マッピング数を最適化するコンピュータプログラムであって、コンピュータに下記の第１の機能〜第５の機能の全て又は一部；
    （Ａ） 被験者由来のＤＮＡのリード情報が、リードの配列及びリードのマッピング先である当該遺伝子座群又は個別の遺伝子座のアリルの観測データとして記録されている記録部から当該データを読み出す、第１の機能、
    （Ｂ） 上記第１の機能によって読み出した前記観測データに基づき、下記の初期化処理（Ｂ）−１及び（Ｂ）−２のいずれかを実行する、第２の機能、
    （Ｂ）−１：当該遺伝子座群又は個別の遺伝子座のアリル頻度に関する変数θの初期値の算出処理、
    （Ｂ）−２：上記変数θ及び、観測データである被験者のＤＮＡのリード情報が混在したデータにおけるリードの塩基配列を媒介する２種の潜在変数としての下記（ａ）及び（ｂ）：
    （ａ）リードｎの当該遺伝子座群又は個別の遺伝子座のアリル選択に関する、θに依存する変数Ｔ_ｎ、
    （ｂ）リードｎの開始位置に関する、Ｔ_ｎに依存するＳ_ｎ、
    が要約された、指標変数Ｚ_ｎｔｓ（Ｚ_ｎｔｓは、（Ｔ_ｎ，Ｓ_ｎ）＝（ｔ，ｓ）の場合１であり、それ以外は０である。）、又は、潜在変数Ｔ_ｎが要約された、Ｚ_ｎｔ（Ｚ_ｎｔは、Ｔ_ｎ＝ｔの場合１であり、それ以外は０である。）＝１の事後確率の初期値の算出処理、
    （Ｃ） 上記第２の機能の（Ｂ）−１で算出した変数θに基づき、当該指標変数Ｚ_ｎｔｓ又はＺ_ｎｔ＝１の事後確率の算出処理を行う、第３の機能、
    （Ｄ） 上記第２の機能の（Ｂ）−２、又は、第３の機能により算出した当該指標変数Ｚ_ｎｔｓ又はＺ_ｎｔ＝１の事後確率に基づいて変数θの最尤推定値の第１の更新値を算出する、第４の機能、
    （Ｅ） 上記第４の機能で算出した変数θの最尤推定値の第１の更新値に基づいて上記第３の機能と第４の機能を再度実行し、さらに、変数θの第２の更新値を算出するループ処理を、新たな更新値と前回の更新値との間における差異が実質的に認められなくなるまで繰り返し実行し、収束した変数θが最適化されたθとして、上記記録部に記録を行う、第５の機能；
    を実現させるアルゴリズムが含まれることを特徴とするコンピュータプログラム。
28. 選択された遺伝子座群又は個別の遺伝子座のアリル由来のＤＮＡのリード情報が混在したデータについて、個々のリードにおける個々の当該遺伝子座群又は個別の遺伝子座のアリルに対する期待マッピング数を最適化するコンピュータプログラムであって、コンピュータに下記の第１の機能〜第５の機能の全て又は一部；
    （Ａ） 被験者由来のＤＮＡのリード情報が、リードの配列及びリードのマッピング先である当該遺伝子座群又は個別の遺伝子座のアリルの観測データとして記録されている記録部から当該データを読み出す、第１の機能、
    （Ｂ） 上記第１の機能によって読み出した前記観測データに基づき、下記の初期化処理（Ｂ）−１及び（Ｂ）−２のいずれかを実行する、第２の機能、
    （Ｂ）−１：当該遺伝子座群又は個別の遺伝子座のアリル頻度についての予備知識の分布を示すハイパーパラメータα_ｔの初期値に基づくθ_ｔの事後分布の更新値の算出処理、
    （Ｂ）−２：上記θの分布、及び、観測データである被験者由来のＤＮＡのリード情報が混在したデータにおけるリードの塩基配列を媒介する２種の潜在変数としての下記（ａ）及び（ｂ）：
    （ａ）リードｎの当該遺伝子座群又は個別の遺伝子座のアリル選択に関する、θに依存する変数Ｔ_ｎ、
    （ｂ）リードｎの開始位置に関する、Ｔ_ｎに依存するＳ_ｎ、
    が要約された、指標変数Ｚ_ｎｔｓ（Ｚ_ｎｔｓは、（Ｔ_ｎ，Ｓ_ｎ）＝（ｔ，ｓ）の場合１であり、それ以外は０である。）、又は、潜在変数Ｔ_ｎが要約された、Ｚ_ｎｔ（Ｚ_ｎｔは、Ｔ_ｎ＝ｔの場合１であり、それ以外は０である。）の事後分布の初期分布の算出処理、
    （Ｃ） 上記第２の機能の（Ｂ）−１で算出した変数θの分布に基づき、当該指標変数Ｚ_ｎｔｓ又はＺ_ｎｔの事後分布の算出処理を行う、第３の機能、
    （Ｄ） 上記第２の機能の（Ｂ）−２、又は、第３の機能で算出した当該指標変数Ｚ_ｎｔｓ又はＺ_ｎｔの事後分布の更新値に基づいて変数θの第１更新事後分布の更新値を算出する、第４の機能、
    （Ｅ） 上記第４の機能で算出した変数θの第１の更新事後分布に基づいて上記第３の機能と第４の機能を再度実行し、さらに、変数θの第２の更新事後分布を算出するループ処理を、新たに更新された事後分布の期待値と前回に更新された事後分布の期待値との間における差異が実質的に認められなくなるまで繰り返し実行して、収束したθの期待値を最適化された当該遺伝子座群又は個別の遺伝子座のアリル頻度のデータとして、上記記録部に記録を行う第５の機能；
    を実現させるアルゴリズムが含まれることを特徴とするコンピュータプログラム。
29. 上記コンピュータプログラムにおいて、選択された遺伝子座群又は個別の遺伝子座のアリル由来のＤＮＡのリード情報が混在したデータは、被験者のリード情報をデータベースに登録されている当該遺伝子座群又は個別の遺伝子座のアリルの塩基配列に対してマッピングをすることにより得られる、当該遺伝子座群又は個別の遺伝子座のアリルに対する各リードのマッピング対応が特定されたリード情報であって、当該マッピングは、下記（ａ）及び（ｂ）に従って行う機能をコンピュータにおいて実現するアルゴリズムが含まれることを特徴とする、請求項２６〜２８のいずれか１項に記載のコンピュータプログラム。
    （ａ） 被験者のリードの配列情報に対して、ヒト遺伝子の塩基配列に対するマッピングの後、当該遺伝子座群又は個別の遺伝子座のアリルにマッピングされたリードを抽出する機能；
    （ｂ） 機能（ａ）により抽出された当該遺伝子座群又は個別の遺伝子座のアリルにマッピングされたリードの配列情報に対して、データベースに登録されている当該遺伝子座群又は個別の遺伝子座のアリルの塩基配列とのマッピングの後、マッピングされたリードを当該遺伝子座群又は個別の遺伝子座のアリル毎に抽出を行い、当該遺伝子座群又は個別の遺伝子座のアリルに対する各リードのマッピング対応が特定されたリード情報を得る機能。
30. 上記機能（ａ）及び（ｂ）におけるマッピングは、一つのリードが複数の選択された遺伝子座群又は個別の遺伝子座のアリルに対してマッピングされることを許容するマッピングであることを特徴とする、請求項２９に記載のコンピュータプログラム。
31. 機能（ａ）における選択された遺伝子座群又は個別の遺伝子座のアリルにマッピングされたリードに加えて、ヒト遺伝子に対してマッピングがなされなかったリードを併せて抽出処理し、これを機能（ｂ）の再マッピングの対象に含めることを特徴とする、請求項２９又は３０に記載のコンピュータプログラム。
32. 被験者の選択された遺伝子座群又は個別の遺伝子座の遺伝型の判定を行うコンピュータプログラムであって、下記（α）〜（δ）の機能をコンピュータに実現させるためのアルゴリズムが含まれることを特徴とする、コンピュータプログラム。
    （α） 請求項２６〜３１のいずれか１項に記載のコンピュータプログラムの実行により得られた、当該遺伝子座群又は個別の遺伝子座のアリル頻度、及び、全リード深度、を少なくとも読み出す、機能α；
    （β） 前記機能αの実行により読み出した当該遺伝子座群又は個別の遺伝子座のアリル頻度から、当該遺伝子座群又は個別の遺伝子座のアリル毎の個別深度への算出処理を実行し、個々のアリルに対して算出された当該個別深度を割り振る処理を実行する、機能β；
    （γ） 棄却閾値として全リード深度の５〜５０％のいずれかの頻度数を設定し、前記機能βの実行により特定された当該数値以下の個別深度の当該遺伝子座群又は個別の遺伝子座のアリルを、当該遺伝子座群の各遺伝子座又は個別の遺伝子座の遺伝型決定の要素から除外する処理を実行する、機能γ；
    （δ） 下記（δ）−１及び（δ）−２に示す機能δ：
    （δ）−１ 前記機能γの除外処理の実行の後、当該遺伝子座群の各遺伝子座又は個別の遺伝子座の遺伝型決定の対象が１個のアリルについては、当該１個のアリルの個別深度が前記棄却閾値の２倍以上である場合には、このアリルをホモ接合と決定し、又は、前記棄却閾値の２倍より小さい場合はヘテロ接合であると決定する処理を実行し；
    （δ）−２ 前記機能γの除外処理の実行の後、当該遺伝子座群の各遺伝子座又は個別の遺伝子座の遺伝型決定の対象が２個の当該遺伝子座群又は個別の遺伝子座のアリルについては、個別深度が大きな方が小さい方の２倍未満である場合には、両アリルはヘテロ接合であると決定し、又は、個別深度が大きな方が小さい方の２倍以上である場合には、大きな方のアリルはホモ接合であると決定する処理を実行する。
33. 選択された遺伝子座群又は個別の遺伝子座は、ＭＨＣの遺伝子座群又は個別の遺伝子座であることを特徴とする、請求項２６〜３２のいずれか１項に記載のコンピュータプログラム。
34. ＭＨＣはＨＬＡであることを特徴とする、請求項３３に記載のコンピュータプログラム。
35. 請求項２６〜３４のいずれか１項に記載のコンピュータプログラムが記録されていることを特徴とする、コンピュータにおいて読み取り可能な記録媒体。

Images