JPWO2019022018A1 - 多型検出法 - Google Patents

Abstract

本発明において、２つ以上の配列の間における多型を検出する方法が提供される。本発明の方法は、配列データにおける個々の配列の全長配列における位置を考慮することなく、複数の配列データ間での多型の検出を可能にする。本発明の方法は、配列データ中の個々の配列（例えば、次世代シーケンサーからのショートリード）を連結してより長い配列とすること（例えば、アセンブリ）を必要とせずに、多型を検出することができることを１つの特徴とする。

Description

本発明は、配列情報、とりわけ、ゲノム等の生体分子の配列情報の情報処理の分野に関する。

次世代シーケンサーの出現により、生物の全ゲノム配列情報が得られるようになった。次世代シーケンサーの配列情報から多型情報を得て、表現型との関連を調べることにより、その表現型の原因となる遺伝子の特定につながる。正確な多型情報の取得は、作物育種のみならず、ヒトの遺伝病の診断、生物種・品種等の特定等、幅広い分野で必要とされる基盤技術であり、これまでにない精度で多型情報が得られれば、そのインパクトは大きい。

次世代シーケンサーからの塩基配列データを用いた多型の検出は、まず最初に配列データをｂｗａ、またはｂｏｗｔｉｅのようなマッピングプログラムを用いてリファレンス配列上の位置情報とミスマッチの情報を得て、次に、ＳａｍｔｏｏｌｓやＧＡＴＫ等の多型抽出プログラムでＳＮＰやｉｎｄｅｌ等の多型情報を抽出するのが一般的である。

これらの方法では、多型の可能性のある部分は可能な限り出力するため、多くのノイズを含みこれらの技術のみでは、正確な多型解析が困難である。マイクロアレイやＤＮＡチップ等の別の技術を併用して用いられているというのが現状である。

本発明において、２つ以上の配列の間における多型を検出する方法が提供される。本発明の方法は、配列データにおける個々の配列の全長配列における位置を考慮することなく、複数の配列データ間での多型の検出を可能にする。本発明の方法は、配列データ中の個々の配列（例えば、次世代シーケンサーからのショートリード）を連結してより長い配列とすること（例えば、アセンブリ）を必要とせずに、多型を検出することができることを１つの特徴とする。

例えば、本発明は以下の項目を提供する。
（項目１） 対象配列データにおいてコントロール配列データに対する多型を検出する方法であって、
ａ）該対象配列データの長さｋの部分配列のサブセットの各々の部分配列の出現頻度を提供する工程と、
ｂ）該コントロール配列データの長さｋの部分配列のサブセットの各々の部分配列の出現頻度を提供する工程と、
ｃ）対象配列とコントロール配列とを比較し、該出現頻度の分布の比較に基づいて、多型を検出する工程と
を包含し、ここで、ｋは該対象配列および該コントロール配列のいずれか短いほうの全長以下の整数である、方法。
（項目２） 前記部分配列中の長さｋ−ｘの配列部分が共通する配列ごとに、長さｘの部分について出現頻度の分布を算出する工程をさらに含み、ここで、ｘはｋ未満の正の整数である、前記項目に記載の方法。
（項目３） 前記比較が、前記部分配列中の長さｋ−ｘの配列部分が共通する配列における、長さｘの部分の出現頻度の分布の差異の比較を含む、前記項目のいずれかに記載の方法。
（項目４） 前記部分配列中の長さｋ−ｘの配列部分を、ユニークな配列ごとにグルーピングする工程を含み、ここで、ｘはｋ未満の正の整数である、前記項目のいずれかに記載の方法。
（項目５） 前記長さｋ−ｘの配列部分をソートする工程を含む、前記項目のいずれかに記載の方法。
（項目６） 前記長さｋ−ｘの配列部分を文字列としてソートする工程を含む、前記項目のいずれかに記載の方法。
（項目７） 前記ｋが、前記対象配列における偶然同一を排除する長さである、前記項目のいずれかに記載の方法。
（項目８） 前記対象配列データおよび前記コントロール配列データが、生物のゲノムに由来する塩基配列データであり、前記ｋが、前記生物のゲノムにおいて、異なる箇所での偶然同一を排除する長さである、前記項目のいずれかに記載の方法。
（項目９） 長さｘが１〜２である、前記項目のいずれかに記載の方法。
（項目１０） 長さｘが１である、前記項目のいずれかに記載の方法。
（項目１１） 前記長さｘの部分が、前記部分配列の末端に存在する、前記項目のいずれかに記載の方法。
（項目１２） 前記対象配列データおよび前記コントロール配列データが塩基配列データであり、前記長さｘの部分が、前記部分配列の３’末端である、前記項目のいずれかに記載の方法。
（項目１３） 前記コントロール配列データのサブセットと前記対象配列データのサブセットとの間で、前記長さｘの部分の配列の出現頻度が異なる場合、該長さｘの部分の配列を、対象配列データにおけるコントロール配列データに対する多型として検出する、前記項目のいずれかに記載の方法。
（項目１４） 前記コントロール配列データのサブセットと前記対象配列データのサブセットとの間で、前記長さｘの部分の配列で最も高頻度のものが異なっている長さｋ−ｘの配列部分が存在する場合、該長さｘの部分の配列を、対象配列データにおける多型として検出する、前記項目のいずれかに記載の方法。
（項目１５） 前記対象配列データおよび前記コントロール配列データが塩基配列データであり、前記対象配列データのサブセットにおける前記長さｘの部分の配列で、前記コントロール配列データのサブセットにおける最も高頻度のものと同一の長さｘの部分の配列がノイズ以下のカウントしか存在しない長さｋ−ｘの配列部分が存在する場合、該長さｘの部分の配列を、対象配列データにおけるホモ多型として検出する、前記項目のいずれかに記載の方法。
（項目１６） 前記対象配列データおよび前記コントロール配列データが塩基配列データであり、対象サブセットにおける前記長さｘの部分の配列で、コントロール配列データのサブセットにおける最も高頻度のものと同一の長さｘの部分の配列が存在し、かつ、コントロール配列データのサブセットにおける最も高頻度のものと異なる長さｘの部分の配列が存在する長さｋ−ｘの配列部分が存在する場合、該長さｘの部分の配列を、対象配列データにおけるヘテロ多型として検出する、前記項目のいずれかに記載の方法。
（項目１７） 対象配列データ量から予測される出現頻度と比較して、前記出現頻度が少ない部分配列をノイズとする、前記項目のいずれかに記載の方法。
（項目１８） 前記対象配列データおよび前記コントロール配列データが塩基配列データであり、［（対象配列データ量）×（１−精度）］／（対象ゲノムサイズ）＋１未満の出現頻度の部分配列をノイズとする、前記項目のいずれかに記載の方法。
（項目１９） 前記対象配列データが、次世代シーケンシングによって得られた塩基配列データである、前記項目のいずれかに記載の方法。
（項目２０） 前記対象配列データが、個体から得られた配列データであり、前記コントロール配列データが、該個体と同種の別の個体、またはデータベースから得られた配列データである、前記項目のいずれかに記載の方法。
（項目２１） 前記対象配列データが、個体の組織試料から得られた配列データであり、前記コントロール配列データが、該個体の別の組織、またはデータベースから得られた配列データである、前記項目のいずれかに記載の方法。
（項目２２） 前記対象配列データが、細胞試料から得られた配列データであり、前記コントロール配列データが、別の細胞、またはデータベースから得られた配列データである、前記項目のいずれかに記載の方法。
（項目２３） 前記多型が、置換、挿入、欠失、コピー数多型（Ｃｏｐｙ Ｎｕｍｂｅｒ Ｖａｒｉａｔｉｏｎ， ＣＮＶ）、ＳＴＲＰ（ｓｈｏｒｔ ｔａｎｄｅｍ ｒｅｐｅａｔ ｐｏｌｙｍｏｒｐｈｉｓｍ）、逆位または転座である、前記項目のいずれかに記載の方法。
（項目２４） 前記多型が、置換である、前記項目のいずれかに記載の方法。
（項目２５） 前記対象配列に対するリファレンス配列における前記多型の位置を特定する工程をさらに含む、前記項目のいずれかに記載の方法。
（項目２６） 前記対象配列データおよび前記コントロール配列データが、生物のゲノムに由来する塩基配列データであり、前記多型のゲノム上の位置を特定する工程をさらに含む、前記項目のいずれかに記載の方法。
（項目２７） 検出された多型の部位について、リファレンス配列またはコントロール配列から作成したクエリ配列セットを用いて、対象配列データおよび／またはコントロール配列データとの比較を行い確認する工程をさらに含む、前記項目のいずれかに記載の方法。
（項目２８） 前記クエリ配列セットが、リファレンス配列またはコントロール配列において前記多型に該当する部位の文字を異なる文字に置換した変異型クエリ配列セットを含む、前記項目のいずれかに記載の方法。
（項目２９） 前記対象配列データおよび前記コントロール配列データが塩基配列データであり、検出された多型の部位について、対象配列データおよび／またはコントロール配列データの相補鎖の配列データに対して、リファレンス配列またはコントロール配列から作成したクエリ配列セットとの比較を行い確認する工程をさらに含む、前記項目のいずれかに記載の方法。
（項目３０） 前記対象配列データおよび前記コントロール配列データが塩基配列データであり、検出された多型の部位について、対象配列データおよび／またはコントロール配列データの変異型の塩基を有する配列データに対して、リファレンス配列またはコントロール配列から作成したクエリ配列セットとの比較を行い確認する工程をさらに含む、前記項目のいずれかに記載の方法。
（項目３１） 前記対象配列データおよび前記コントロール配列データが、生物のゲノムに由来する塩基配列データであり、前記ゲノムの配列が不明である、前記項目のいずれかに記載の方法。
（項目３２） 実験結果またはデータベースから対象配列データまたはコントロール配列データを取得する工程をさらに含む、前記項目のいずれかに記載の方法。
（項目Ｘ１） 対象配列データにおけるコントロール配列データに対する多型を含む部分配列中の多型ではない部分の少なくとも一部を含む配列を、該多型の識別子として割り当てることをさらに含む、前記項目のいずれか１項に記載の方法。
（項目Ｘ２） 前記多型の識別子をリファレンス配列にマッピングし、リファレンス上の該多型の位置を特定することを含む、前記項目のいずれかに記載の方法。
（項目３３） 対象配列データにおいてコントロール配列データに対する多型を検出する方法をコンピュータに実行させるためのプログラムであって、該方法は、
ａ）該対象配列データの長さｋの部分配列のサブセットをコンピュータに保存する工程であって、ここで、ｋは対象配列およびコントロール配列の全長以下の長さである、工程と、
ｂ）該対象配列データの長さｋのサブセットの各々の部分配列の出現頻度を算出する工程と、
ｃ）該コントロール配列データの長さｋの部分配列のサブセットにおける各々の部分配列の出現頻度をコンピュータに保存する工程と、
ｄ）対象配列とコントロール配列とを比較し、該出現頻度の分布の差異に基づいて、多型を検出する工程と
を包含する、プログラム。
（項目３３Ａ） 前記項目のいずれか１つまたは複数に記載される特徴を有する、前記項目に記載のプログラム。
（項目３４） 前記方法が、前記部分配列中の多型ではない部分の少なくとも一部を含む配列（前記部分配列全体であり得る。）を、検出された前記多型の名称として表示する工程をさらに含む、前記項目のいずれかに記載のプログラム。
（項目３５） 対象配列データにおいてコントロール配列データに対する多型を検出する方法をコンピュータに実行させるためのプログラムを格納する記録媒体であって、該方法は、
ａ）該対象配列データの長さｋの部分配列のサブセットをコンピュータに保存する工程であって、ここで、ｋは対象配列およびコントロール配列の全長以下の長さである、工程と、
ｂ）該対象配列データの長さｋのサブセットの各々の部分配列の出現頻度を算出する工程と、
ｃ）該コントロール配列データの長さｋの部分配列のサブセットにおける各々の部分配列の出現頻度をコンピュータに保存する工程と、
ｄ）対象配列とコントロール配列とを比較し、該出現頻度の分布の差異に基づいて、多型を検出する工程と
を包含する、記録媒体。
（項目３５Ａ） 前記項目のいずれか１つまたは複数に記載される特徴を有する、前記項目に記載の記録媒体。
（項目３６） 前記方法が、前記部分配列中の多型ではない部分の少なくとも一部を含む配列（前記部分配列全体であり得る。）を、検出された前記多型の名称として表示する工程をさらに含む、前記項目のいずれかに記載の記録媒体。
（項目３７） 対象配列データにおいてコントロール配列データに対する多型を検出するためのシステムであって、該システムは、
該対象配列データおよび該コントロール配列データの長さｋの部分配列のサブセットの各々の部分配列の出現頻度を提供するように構成された配列データ処理部であって、ここで、ｋは対象配列およびコントロール配列の全長以下の長さである、配列データ処理部と、
対象配列とコントロール配列とを比較し、該出現頻度の分布の差異に基づいて、多型を検出する工程とを行うように構成された、配列データ計算部と
を備える、システム。
（項目３７Ａ） 前記項目のいずれか１つまたは複数に記載される特徴を有する、前記項目に記載のシステム。
（項目３８） 前記システムが、前記部分配列中の多型ではない部分の少なくとも一部を含む配列（前記部分配列全体であり得る。）を、検出された前記多型の名称として表示する表示手段をさらに含む、前記項目のいずれかに記載のシステム。
（項目３９） 対象配列データにおいてコントロール配列データに対する多型を検出する方法であって、
（１）ａ）該対象配列データの長さｋの部分配列のサブセットの各々の部分配列の出現頻度を提供する工程と、
ｂ）該コントロール配列データの長さｋの部分配列のサブセットの各々の部分配列の出現頻度を提供する工程と、
ｃ）対象配列とコントロール配列とを比較し、該出現頻度の分布の比較に基づいて、多型を検出する工程と
によって、対象配列データにおける置換、コピー数多型、ＳＴＲＰ、挿入、欠失、逆位または転座を検出するプロセスと、
（２）ａ）該対象配列データの配列中の少なくとも２ヶ所の部分配列の、該コントロール配列上の位置を特定する工程と、
ｂ）対象配列データにおける該部分配列間の位置関係と、コントロール配列上の該部分配列間の位置関係とを比較する工程と、
ｃ）対象配列データにおける該部分配列間の位置関係と、コントロール配列上の該部分配列間の位置関係とが異なっている場合、目的とする多型があると判定し、該対象配列データにおける該部分配列部位間の文字を、対応するコントロール配列上の文字と、該部分配列部位を始点として順次比較して不一致となる部位を検出する工程と
によって、対象配列データにおける挿入、欠失、逆位、転座または置換を検出するプロセスと
を包含する、方法。
（項目３９Ａ） 前記項目のいずれか１つまたは複数に記載される特徴を有する、前記項目に記載の方法。
（項目４０） 対象配列データにおいてリファレンス配列データに対する多型を検出する方法であって、リファレンス配列データから、各々の位置情報と関連付けられたリファレンス配列のｋ長の部分配列セットを作成する工程を含み、さらに、
（Ａ１）該対象配列データの長さｋの部分配列のサブセットを生成し、ユニークな長さｋの部分配列の出現頻度を提供する工程と、
（Ａ２）該リファレンス配列のｋ長の部分配列セットの、ユニークな長さｋの部分配列の出現頻度を提供する工程と、
（Ａ３）該対象配列と該リファレンス配列とを比較し、該出現頻度の分布の比較に基づいて、挿入、欠失、置換、コピー数多型、ＳＴＲＰ、逆位または転座を検出する工程とを包含するプロセスと
（Ｂ１）該対象配列データの配列中の少なくとも２ヶ所のｋ長の部分配列をクエリとして、該リファレンス配列のｋ長の部分配列セットに対して検索を行い、該少なくとも２ヶ所の部分配列の、リファレンス配列上の位置を特定する工程と、
（Ｂ２）該対象配列データにおける該部分配列間の位置関係と、該リファレンス配列上の該部分配列間の位置関係とを比較する工程と、
（Ｂ３）該対象配列データにおける該部分配列間の位置関係と、該リファレンス配列上の該部分配列間の位置関係とが異なっている場合、挿入、欠失、逆位または転座があると判定し、該対象配列データにおける該部分配列部位間の文字を、対応するコントロール配列上の文字と、該部分配列部位を始点として順次比較して不一致となる部位を検出する工程を包含し、必要に応じて、
（Ｂ４）該位置関係が異ならない場合に、該対象配列データにおける該部分配列部位間の文字を、対応する前記コントロール配列上の文字と比較して不一致となる部位を検出する工程をさらに含み、不一致となる部位が存在する場合、置換が存在すると判定する工程をさらに含む、プロセスと、
を、同時に、並行して、または逐次的に行うことを特徴とする、方法。
（項目４０Ａ） 前記項目のいずれか１つまたは複数に記載される特徴を有する、前記項目に記載の方法。
（項目Ａ１） 対象配列データとコントロール配列データとの比較方法であって、
対象配列データにおけるコントロール配列データに対する多型を含む部分配列中の多型ではない部分の少なくとも一部を含む配列を、該多型の識別子として割り当てることを含む、方法。
（項目Ａ１Ａ） 前記項目のいずれか１つまたは複数に記載される特徴を有する、前記項目に記載の方法。
（項目Ａ２） 前記多型の識別子をリファレンス配列にマッピングし、リファレンス上の該多型の位置を特定することを含む、前記項目のいずれかに記載の方法。

本発明において、上記１または複数の特徴は、明示された組み合わせに加え、さらに組み合わせて提供され得ることが意図される。本発明のなおさらなる実施形態および利点は、必要に応じて以下の詳細な説明を読んで理解すれば、当業者に認識される。

本発明により、２つ以上の配列の間で、全長配列における位置を考慮する必要なく、正確に多型、特に置換を検出することができる。ｋ長配列を用いた置換変異の検出に関しては、ゲノムマッピングを行う前に多型検出ができることが１つの大きな特徴となる。そして、リファレンス配列が存在しない生物でも多型検出が可能であり、ｋ−ｍｅｒ自体を多型の名称として用いることが可能であるため、それゆえ連鎖解析等の遺伝解析で大きな変革をもたらす可能性がある。

図１は、本発明の方法の１つの実施形態を例示するフロー図である。図１においては、本発明の方法において行われ得る主な工程の概要が示される。 図２は、本発明の方法の１つの実施形態を例示するフロー図である。図２においては、リファレンス配列の端からｋ−ｍｅｒ（ｋ＝２０）を順に得て、対象配列および参照配列のｋ−ｍｅｒ（ｋ＝２０）の頻度を位置情報と一緒に出力する工程の例が示される。このような出力を用いることにより、ＣＮＶや挿入・欠失、置換を検出することができる。 図３は、本発明の方法の１つの実施形態を例示するフロー図である。図３においては、リファレンス配列、コントロール配列（参照配列）および対象配列からのｋ長部分配列のサブセットおよび各々の部分配列の出現頻度を提供する工程の例が示される。 図４は、本発明の方法の１つの実施形態を例示するフロー図である。図４においては、対象配列とコントロール配列とを比較し、該出現頻度の分布の比較に基づいて、多型を検出する工程の例が示される。 図５は、本発明の方法の１つの実施形態を例示するフロー図である。図５においては、リファレンス配列における多型の位置を特定する工程の例が示される。 図６は、本発明の方法の１つの実施形態を例示するフロー図である。図６においては、検出された多型を確認する工程の例が示される。 図７は、コントロール配列（N1）と対象配列（N1S7）の間でのｋ−ｍｅｒ配列の比較結果の一部を示す。コントロール配列および対象配列のｋ長部分配列サブセットのそれぞれの部分配列と、当該部分配列の各々がコントロール配列および対象配列のｋ長部分配列サブセットに出現する頻度が示される。当該比較によって、コントロール配列由来の配列は対象配列では検出されず、逆に対象配列由来の配列は参照配列では検出されていないことから、多型を検出することができることが示される。下線の塩基がコントロールと対象とで異なっており、多型を示している。図７は、変異がホモ型の場合に観察される結果の一例である。図７における参照配列のｋ−ｍｅｒ部分配列は上からそれぞれ配列番号１〜６０に対応し、対象配列のｋ−ｍｅｒ部分は、上から配列番号１〜２０、６１〜８０、４０〜６０に対応する。 図８は、コントロール配列（N1）と対象配列（N1S7）の間でのｋ−ｍｅｒ配列の比較結果の一部を示す。コントロール配列および対象配列のｋ長部分配列サブセットのそれぞれの部分配列と、当該部分配列の各々がコントロール配列および対象配列のｋ長部分配列サブセットに出現する頻度が示される。当該比較によって、コントロール配列由来の配列は対象配列では検出されず、逆に対象配列由来の配列は参照配列では検出されていないことから、多型を検出することができることが示される。下線の塩基がコントロールと対象とで異なっており、多型を示している。図８は、変異がヘテロ型の場合に観察される結果の一例である。図８における参照配列のｋ−ｍｅｒ部分配列は上からそれぞれ配列番号８１〜１４０に対応し、対象配列のｋ−ｍｅｒ部分は、上から配列番号８１〜１００、１４１〜１６０、１２０〜１４０に対応する。 図９は、コントロール配列（N1）と対象配列（N1S5、N1S6、N1S7、N1S10）のｋ−ｍｅｒ配列の部分配列サブセットを整列させることによって、イネリファレンス配列の対応する位置から始まる配列と対応する配列の出現頻度を求めた結果を示す図である。ｋ−ｍｅｒ配列の出現頻度の変化から、N1S7においてヘテロ変異が生じ、N1S10においてホモ変異となったことを検出することができる。 図１０は、コントロール配列（N1）と対象配列（N1S5、N1S6、N1S7、N1S10）との間で、Ｐｏｌｙｍｏｒｐｈｉｃ Ｅｄｇｅ Ｄｅｔｅｃｔｉｏｎによって多型を検出・確認した結果の一部を示す。Ｃｈｒは染色体番号、Ｐｏｓは染色体上の位置を示し、Ｒｅｆはリファレンス配列での当該位置の塩基、Ａｌｔは変異が存在する場合の当該位置の塩基を示す。Ｐｏｌｙｍｏｒｐｈｉｃ Ｅｄｇｅ Ｄｅｔｅｃｔｉｏｎの欄では、Ｃｏｎｔｒｏｌ配列（N1）の２０ｍｅｒ部分配列における最終塩基の出現頻度と、各対象（Ｔａｒｇｅｔ）配列の２０ｍｅｒ部分配列における最終塩基の出現頻度が示される。最も高頻度の最終塩基が異なっている部分を、多型として検出している。Ｖｅｒｉｆｙの欄では、リファレンス配列から作成したクエリ配列セット（Ｒｅｆ：野生型クエリセット、Ａｌｔ：変異型クエリセット）に対するＣｏｎｔｒｏｌ配列データおよびＴａｒｇｅｔ配列データの出現頻度が示される。最終的に判定した各部位における各試料の遺伝子型を、Ｇｅｎｏｔｙｐｅの欄に示す（Ｍ：ホモ変異、Ｈ：ヘテロ変異、Ｗ：野生型）。 図１１は、コントロール配列（N1）と対象配列（N1S5、N1S6、N1S7、N1S10）との間で、Ｐｏｌｙｍｏｒｐｈｉｃ Ｅｄｇｅ Ｄｅｔｅｃｔｉｏｎによって多型を検出・確認した結果の一部を示す。Ｃｈｒは染色体番号、Ｐｏｓは染色体上の位置を示し、Ｒｅｆはリファレンス配列での当該位置の塩基、Ａｌｔは変異が存在する場合の当該位置の塩基を示す。Ｐｏｌｙｍｏｒｐｈｉｃ Ｅｄｇｅ Ｄｅｔｅｃｔｉｏｎの欄では、Ｃｏｎｔｒｏｌ配列（N1）の２０ｍｅｒ部分配列における最終塩基の出現頻度と、各対象（Ｔａｒｇｅｔ）配列の２０ｍｅｒ部分配列における最終塩基の出現頻度が示される。最も高頻度の最終塩基が異なっている部分を、多型として検出している。Ｖｅｒｉｆｙの欄では、リファレンス配列から作成したクエリ配列セットに対するＣｏｎｔｒｏｌ配列データおよびＴａｒｇｅｔ配列データの出現頻度が示される。最終的に判定した各部位における各試料の遺伝子型を、Ｇｅｎｏｔｙｐｅの欄に示す（Ｍ：ホモ変異、Ｈ：ヘテロ変異、Ｗ：野生型）。 図１２は、各世代のイネサンプルについてＰＣＲ法で対象となる多型部分を増幅して、サンガー法で塩基配列を確認した結果を示す図である。各試料の遺伝子型は、Ｍ：ホモ変異、Ｈ：ヘテロ変異、Ｗ：野生型として示される。Ｃｈｒは染色体番号、Ｐｏｓは染色体上の位置を示し、Ｒｅｆはリファレンス配列での当該位置の塩基、Ａｌｔは変異が存在する場合の当該位置の塩基を示す。N1、N1S1、N1S2、N1S3、N1S4、N1S5、N1S6、N1S7、N1S8、N1S9、N1S10と世代を重ねながら、ヘテロ変異が生じ、その後ホモ変異として定着する様子が観察されることが分かる。また、本発明の方法によって多型を検出した結果とよく一致しており、本発明の方法による多型検出が高い精度を有していることが理解される。 図１３は、ナイジェリアのヨルバ族男性（NA18507）の配列データを用いて、リファレンスゲノム配列データから作成された参照（コントロール）配列に対する多型を検出した解析結果の一部を示す。対象配列データは、すでにＩｌｌｕｍｉｎａ社の次世代シーケンサーで解析されてＮＣＢＩに登録・公開されたデータをダウンロードして用いた。当該塩基配列セットの実験ＩＤのＵＲＬは、https://www.ncbi.nlm.nih.gov/sra/SRX016231であり、配列のアクセッション番号は、SRR034939〜SRR034975の範囲であった。ｋ−１長の配列は、それぞれ上から配列番号１６１〜１９０に対応する。 図１４は、本発明において行われ得る確認工程における、クエリ配列セットの作成の模式図である。上部の配列がリファレンス配列であり、下線太字で示されるＴが検出された多型部位を示す。当該多型部位を含む部分配列のセットを生成し、クエリ配列セットとすることができる。各配列は、それぞれ上から配列番号２６７〜２７５に対応する。 図１５Ａは、本発明のシステムの実施形態を模式的に示した図である。 図１５Ｂは、本発明のシステムのさらなる実施形態を模式的に示した図である。 図１６は、本発明の方法の実施形態を模式的に示した図である。 図１７は、本発明の方法によるコピー数多型（ＣＮＶ）の検出の結果の一部を示す図である。イネ第７染色体の２６６９４７９５位置（図中、矢印で示す）からコピー数多型部位が開始されている。検出されたコピー数多型部位は、レトロトランスポゾンTos17に対応する。このトランスポゾンは４．１ｋｂあるため、図１７には最初のジャンクションの部分のみ示している。培養時間に応じた転移によるコピー数の増加が本発明の方法によって検出されていることが理解される。 図１８は、ｋ−ｍｅｒ配列の頻度を用いる多型検出フローと、部分配列の位置関係を用いる多型検出フローとを組み合わせて行う場合の一実施形態を示すフロー図である。

以下、本発明を最良の形態を示しながら説明する。本明細書の全体にわたり、単数形の表現は、特に言及しない限り、その複数形の概念をも含むことが理解されるべきである。従って、単数形の冠詞（例えば、英語の場合は「ａ」、「ａｎ」、「ｔｈｅ」など）は、特に言及しない限り、その複数形の概念をも含むことが理解されるべきである。また、本明細書において使用される用語は、特に言及しない限り、当該分野で通常用いられる意味で用いられることが理解されるべきである。したがって、他に定義されない限り、本明細書中で使用される全ての専門用語および科学技術用語は、本発明の属する分野の当業者によって一般的に理解されるのと同じ意味を有する。矛盾する場合、本明細書（定義を含めて）が優先する。

（定義）
以下に本明細書において特に使用される用語の定義および／または基本的技術内容を適宜説明する。

本明細書において、「配列」とは、各々が何らかの値を取る複数の変数であって、それら複数の変数の順序の情報をさらに含むものをいう。代表的には文字列で表示される。

本明細書において、「対象配列」とは、多型を検出しようとする任意の配列をいい、本明細書においては、「ターゲット」、「ターゲット配列」、「ｔａｒｇｅｔ」とも表記する場合がある。

本明細書において、「コントロール配列」とは、その配列との差異を多型として検出するための基準として用いられる任意の配列をいい、本明細書においては、「コントロール」、「参照配列」、「比較配列」、「ｃｏｎｔｒｏｌ」とも表記する場合がある。

本明細書において、「多型（ｐｏｌｙｍｏｒｐｈｉｓｍ）」とは、対象配列中においてコントロール配列と異なっている任意の部分を指す。本明細書において、「変異」も同様の意味で使用することができる。

本明細書において、「リファレンス（ｒｅｆｅｒｅｎｃｅ）配列」とは、対象配列および／またはコントロール配列の全長の配列として扱うことができる配列を指す。いかなる配列を全長配列とするかは、対象配列および／またはコントロール配列として用いる配列に応じて適宜決定されるものであり、例示されるものに限定されないが、例えば、ウェブ上のデータベース等に存在する、全ゲノム配列、染色体全長配列、遺伝子全長配列、プラスミド全長配列、エクソン全長配列、タンパク質全長配列などをリファレンス配列として用いることができる。

本明細書において、「配列データ」とは、ある配列についての情報を与えるデータをいう。代表的には、配列そのものも配列データということができ、また、配列の一部について情報を与えるデータ（例えば、ゲノム配列に対するシーケンシングによる解析データ）も配列データとして包含される。

本明細書において、ある配列の「部分配列」とは、その配列に含まれる任意の配列をいう。

本明細書において、「サブセット」とは、配列の集合と、それらの配列の部分配列の集合とを合わせた集合の任意の部分集合をいう。

本明細書において、「次世代シーケンシング」とは、配列決定プロセスを並列化し、一度のランで数千万から数億の配列データを生成するシーケンシング技法である。「次世代シーケンサー」とは、次世代シーケンシングを行うための機器を指す。

「偶然同一を排除する」とは、ある配列と、偶然に同一の配列が出現する期待値を１未満にすることをいう。

本明細書において、「カバレッジ」とは、配列データの量が、配列全長の何倍に相当しているかを指す。「カバー率」、「〜倍の読み」などと称される場合もある。

本明細書において、「配列構造体」とは、配列中における、物理的に分離された一連の配列をいう。例えば、ゲノム配列の文脈では、染色体のそれぞれは配列構造体ということができる。

本明細書において、「転座」とは、複数の配列構造体を有する配列中で、ある配列構造体上の部分配列が、他の配列構造体上に移動している多型をいう。

本明細書において、「ジャンクション」とは、一部が同一である２つの配列について、同一である部分と同一でない部分の境界を指す。

本明細書において、「識別子」とは、ある多型を他の多型と区別するために付される名称を指す。一般的には、多型の開始位置と型で記載されることが多いが、本明細書において記載される識別子を用いることができる。

本明細書において、「エッジ」とは、配列において多型を含む部分の末端をさす。

（好ましい実施形態）
以下に本発明の好ましい実施形態を説明する。以下に提供される実施形態は、本発明のよりよい理解のために提供されるものであり、本発明の範囲は以下の記載に限定されるべきでないことが理解される。従って、当業者は、本明細書中の記載を参酌して、本発明の範囲内で適宜改変を行うことができることは明らかである。また、本発明の以下の実施形態は単独でも使用されあるいはそれらを組み合わせて使用することができることが理解される。

なお、以下で説明する実施の形態は、いずれも包括的または具体的な例を示すものである。以下の実施の形態で示される数値、形状、材料、構成要素、構成要素の配置位置及び接続形態、ステップ、ステップの順序などは、一例であり、請求の範囲を限定する主旨ではない。また、以下の実施の形態における構成要素のうち、最上位概念を示す独立請求項に記載されていない構成要素については、任意の構成要素として説明される。

（本発明の多型検出の概要）
本発明は、対象配列データにおいてコントロール配列データに対する多型を検出する方法を提供する。この方法は、ａ）該対象配列データの長さｋの部分配列のサブセットの各々の部分配列の出現頻度を提供する工程と、ｂ）該コントロール配列データの長さｋの部分配列のサブセットの各々の部分配列の出現頻度を提供する工程と、ｃ）対象配列とコントロール配列とを比較し、該出現頻度の分布の比較に基づいて、多型を検出する工程とを包含し、ここで、ｋは該対象配列および該コントロール配列のいずれか短いほうの全長以下の整数である、方法を提供する。本発明の例示的なフローは図１６に説明されている。

本発明の方法は、対象配列データと、コントロール配列データとの２つの配列データ（例えば、次世代シーケンサー解析結果）からの直接比較により多型を検出することが可能である点で、従来法と根本的に異なる。とりわけ、塩基配列における多型検出において、ゲノム上の位置の考慮なしに直接比較する方法は、新規なものであると考えられる。

１つの実施形態において、本発明の方法は、配列データから、一定長（ｋ長）の部分配列のセットを得ることを１つの特徴とする。１つの実施形態において、本発明の方法は、配列データから、ｋ長の部分配列のセットにおける各部分配列の頻度分布を得ることを１つの特徴とする。一部の実施形態では、配列データから、配列を１つずつずらしながらｋ長の部分配列のセットを作成する。

１つの実施形態では、ｋ長の配列のうち、ｋ−ｘ（ｘは１など）の配列が同じであるデータをソートし、異なる部分（ｘ長の部分に該当）の頻度を検出する。Ｘ長の部分は、部分配列中での位置は限定されず、配列中の中央部でもよい。しかしながら、Ｘ長の部分を、部分配列の末端（例えば、塩基配列においては３’末端または５’末端）にすることは、ソートなどの処理を顕著に簡便化・高速化するため、好ましい。このように、ｋ長の配列のうち末端部で異なる部分（ｘ長の部分に該当）の頻度を検出する場合、本明細書において、「Ｐｏｌｙｍｏｒｐｈｉｃ Ｅｄｇｅ Ｄｅｔｅｃｔｉｏｎ（ＰＥＤ）」または「ｅｄｇｅ ｄｅｔｅｃｔｉｏｎ」と呼称される場合がある。

ここで、ｋの値としては、配列データの各配列（例えば、次世代シーケンサーの各々のショートリード）の長さを上限とした任意の値を挙げることができ、例えば、約５００、約４００、約３００、約２００、約１００、約５０、約４０、約３０、約２５、約２０、約１５等を挙げることができる。ｋが増加することにより、ｋ−ｍｅｒ配列のデータは指数的に増加する（例えば、塩基配列の場合、ｋが１塩基増えるごとに塩基の組み合わせは４倍になる）ため、例えば、塩基配列の場合、ｋ＝２０〜２５程度が好ましいが、理論上は、例えば、ｋ＝５００等でも用いることが可能であり、制限されるものではない。ヒトの場合であると、ｋ＝１７以下だと偶然一致が生じる確率が高くなるが、ゲノムサイズが小さな生物であれば、例えば、ｋ＝１５などより小さなｋ値を用いることも可能である。１つの実施形態ではｋ＝２０を用いる。

ｋ−ｘ長の配列が同じであるデータのｘの部分の文字が比較対象間で異なる場合、その文字に多型（変異）が含まれると考えられる。挿入・欠失変異も変異の末端文字を検出できる。例えば、ｋ−ｘ長の塩基配列が同じであるデータのｘの部分の塩基が比較対象間で異なる場合、その塩基に多型（変異）が含まれると考えられる。

一部の実施形態において、得られた配列セットのうち、同じ配列に関して出現回数で整理したデータを算出する。この工程は、計算機を使用して簡便に行うことができ、例えば、Ｕｎｉｘでは以下：

のように実装することができ、（文字列で）ソートされた配列と頻度を示す数値のデータを生成することができる。さらに、対象とコントロールの頻度データを同一ｋ−ｍｅｒでまとめる際に、例えば、Ｕｎｉｘではjoinコマンド等を使用して行うことが可能である。

配列出現回数を、配列データのカバレッジ（何倍読みか）と比較することによって、配列データにおける差異を評価することも可能である。例えば、ゲノム配列に対して４０倍のデータ量の配列解析からの配列データに対して、出現頻度が１となるようなものはノイズと考えることができる。

本発明は、特に、「置換」多型（長さが変わらない＝欠失挿入ではない）の検出には、極めて高い効果を発揮する。多型部位が多コピーの場合は位置の特定ができない場合があり得るが、それでもなお、多型自体の検出は可能であり、多型に名称を付して特定できる。そのため、例えば、検出された各多型を、形質との関係を調べる多型マーカーとして使用することが可能であり、診断・育種・鑑定・品質管理（例えば、ｉＰＳ細胞の品質管理）・分類・検査にも応用可能である。

本発明を、次世代シーケンサーから得られた塩基配列データから直接多型を検出する方法として用いることで、２種類のサンプル間、およびリファレンス配列とサンプルとの間の多型の検出が可能となる。また、長さｋの部分配列、ｋ−ｘの部分配列は重複のないユニークな配列なので、配列自体を多型の識別子（名称）として利用可能である。このため、リファレンスゲノム配列が決定されていないため多型のゲノム上の位置関係が判別できない場合においても、世界共通の一意の多型の識別子として利用できる。本発明の１つの実施形態において、対象配列データとコントロール配列データとの比較方法であって、対象配列データにおけるコントロール配列データに対する多型を含む部分配列中の多型ではない部分の少なくとも一部を含む配列を、多型の識別子として割り当てることを含む、方法が提供される。また、多型の識別子をリファレンス配列にマッピングし、リファレンス上の多型の位置を特定することが可能である。

識別子は、多型ではない部分の少なくとも一部に加えて、多型自体も含み得る。多型塩基を含めた識別子は、リファレンス配列へのマッピングは難しいものの、連鎖解析に用いることが可能である。
例えば、
ＡＡＡＣＣＡＣＴＴＣＡＣＧＴＴＴＣＣＡ Ａ
ＡＡＡＣＣＡＣＴＴＣＡＣＧＴＴＴＣＣＡ Ｇ
という多型の例では、
ＡＡＡＣＣＡＣＴＴＣＡＣＧＴＴＴＣＣＡのＡ型
ＡＡＡＣＣＡＣＴＴＣＡＣＧＴＴＴＣＣＡのＧ型
ＡＡＡＣＣＡＣＴＴＣＡＣＧＴＴＴＣＣＡのＡ／Ｇのヘテロ型
という表現が記載の一例である。
多型を含めた表記の仕方の例としては、
ＡＡＡＣＣＡＣＴＴＣＡＣＧＴＴＴＣＣＡＡ型、
ＡＡＡＣＣＡＣＴＴＣＡＣＧＴＴＴＣＣＡＧ型、
そして、ヘテロ型は、
ＡＡＡＣＣＡＣＴＴＣＡＣＧＴＴＴＣＣＡＡ／ＡＡＡＣＣＡＣＴＴＣＡＣＧＴＴＴＣＣＡＧ
のように、２つの型を併記することが可能である。

本発明の１つの実施形態は、対象配列データにおいてコントロール配列データに対する多型を検出する方法である。１つの実施形態において、当該方法は、該対象配列データの長さｋの部分配列のサブセットの各々の部分配列の出現頻度を提供する工程を含む。ｋは該対象配列および該コントロール配列のいずれか短いほうの全長以下の整数である。１つの実施形態において、当該方法は、該コントロール配列データの長さｋの部分配列のサブセットの各々の部分配列の出現頻度を提供する工程を含む。１つの実施形態において、当該方法は、対象配列とコントロール配列とを比較し、該出現頻度の分布の比較に基づいて、多型を検出する工程とを含む。このような工程によって、全長配列における位置を考慮せず、また、配列を連結することなく配列データを比較し、多型を検出することができる。

さらなる実施形態において、本発明の方法は、部分配列中の長さｋ−ｘの配列部分が共通する配列ごとに、長さｘの部分について出現頻度の分布を算出する工程をさらに含む。ｘはｋ未満の正の整数である。この場合において、出現頻度の分布の比較として、前記部分配列中の長さｋ−ｘの配列部分が共通する配列における、長さｘの部分の出現頻度の分布の差異の比較が含まれ得る。これにより、多型検出の処理を高速化することが可能である。

一部の実施形態において、本発明の方法は、前記部分配列中の長さｋ−ｘの配列部分を、ユニークな配列ごとにグルーピングする工程を含む。これには、例えば、前記長さｋ−ｘの配列部分をソートする工程（例えば、前記長さｋ−ｘの配列部分を文字列としてソートする工程）が含まれ得る。

一部の実施形態において、ｋの値は、前記対象配列データ等における偶然同一を排除する長さである。例えば、前記対象配列データおよび前記コントロール配列データが、生物のゲノムに由来する塩基配列データである場合、前記ｋは、前記生物のゲノムにおいて、異なる箇所での偶然同一を排除する長さであり得る。これにより、多型の検出をより正確なものとすることが可能である。

長さｘは、限定されるものではないが、好ましくは１〜３であり、さらに好ましくは１〜２であり、より好ましくは１である。

１つの実施形態では、前記長さｘの部分が、前記部分配列の末端に存在する。例えば、前記対象配列データおよび前記コントロール配列データが塩基配列データである場合、前記長さｘの部分は、前記部分配列の３’末端または５’末端であり得る。長さｘの部分を部分配列の末端にとることは、比較処理の高速化・簡便化にとって望ましい。

出現頻度の分布の差異の比較により、例えば、以下のような多型の検出が可能である。１つの実施形態では、前記コントロール配列データのサブセットと前記対象配列データのサブセットとの間で、前記長さｘの部分の配列の出現頻度が異なる場合、該長さｘの部分の配列を、対象配列データにおけるコントロール配列データに対する多型として検出する。１つの実施形態では、前記コントロール配列データのサブセットと前記対象配列データのサブセットとの間で、前記長さｘの部分の配列で最も高頻度のものが異なっている長さｋ−ｘの配列部分が存在する場合、該長さｘの部分の配列を、対象配列データにおける多型として検出する。

１つの実施形態では、前記対象配列データおよび前記コントロール配列データが塩基配列データであり、前記対象配列データのサブセットにおける前記長さｘの部分の配列で、前記コントロール配列データのサブセットにおける最も高頻度のものと同一の長さｘの部分の配列がノイズ以下のカウントしか存在しない長さｋ−ｘの配列部分が存在する場合、該長さｘの部分の配列を、対象配列データにおけるホモ多型として検出する。ノイズは、本明細書において後述されるような基準を用いて判定することができる。１つの実施形態では、前記対象配列データおよび前記コントロール配列データが塩基配列データであり、対象サブセットにおける前記長さｘの部分の配列で、コントロール配列データのサブセットにおける最も高頻度のものと同一の長さｘの部分の配列が存在し、かつ、コントロール配列データのサブセットにおける最も高頻度のものと異なる長さｘの部分の配列が存在する長さｋ−ｘの配列部分が存在する場合、該長さｘの部分の配列を、対象配列データにおけるヘテロ多型として検出する。

一部の実施形態において、本発明の方法において、予測される出現頻度と比較して、出現頻度が少ない部分配列をノイズとすることが可能である。これにより、機械的に生じた差異と、実際に生じている多型とを識別して検出することが可能である。例えば、前記対象配列データおよび前記コントロール配列データが塩基配列データである場合には、対象配列データ量、配列データにおける予想されるエラー率（例えば、マニュアル・機器の公称値）、対象配列全長サイズ（例えば、ゲノムサイズ）等を考慮して、予想されるノイズのカウント程度、または予想されるノイズのカウント以下の出現頻度の部分配列をノイズとすることができる。１つの予測としては、生じるノイズの頻度の期待値は、ｉｎｔ（ゲノムを何倍読んだか×（１−精度）＋１）」となる。エラーが存在しない場合、精度は１になる。１つの実施形態では、対象配列データおよび前記コントロール配列データが塩基配列データである場合に、［（対象配列データ量）×（１−精度）］／（対象ゲノムサイズ）＋１未満の出現頻度の部分配列をノイズとすることが可能である。

例えば、エラー率が０．００１の場合は着目した塩基に関して１カウントでも出現する確率は０．００１になるためほぼ０と考えられるが、１０００塩基の範囲で見れば、どれか１塩基にエラーが含まれる計算となる。そのため、繰り上げた整数値をノイズの予測値とすることができると考えられ、ｉｎｔ関数中で１を加えることによって繰り上げた整数値とすることができる。

あるいは、本発明の方法において、予測される出現頻度と比較して、出現頻度が多い部分配列を、リピート配列部位として除外することができる。例えば、対象配列データの対象配列全長のカバレッジ（カバー率）の２倍を超えるような部分配列を解析から除外することができる。

１つの実施形態では、前記対象配列データおよび／またはコントロール配列データが、次世代シーケンシングによって得られた塩基配列データである。次世代シーケンシングにおける多型の検出においては、従来、リファレンスへのマッピングや、配列のアセンブリが必要とされており、このような工程で生じる不確実性によって、多型の検出が大きく阻害されていたため、このような工程を必要としない本発明の方法を次世代シーケンシングから得られた配列データに対して用いることは特に有利なものであるということができる。

対象配列データおよびコントロール配列データは、限定されるものではないが、多型を検出する上では、一定の共通性を持つ配列についての配列データであることが望ましい。しかしながら、配列の取得方法については各々同一でも異なっていてもよく、シーケンシングによって得られたデータ間での比較を行うことも、データベース等から得られたデータ間での比較を行うことも、シーケンシングによって得られたデータとデータベース等から得られたデータとの間での比較を行うことも可能である。

１つの実施形態では、対象配列データが、個体から得られた配列データであり、コントロール配列データが、該個体と同種の別の個体、またはデータベースから得られた配列データである。１つの実施形態では、対象配列データが、個体の組織試料から得られた配列データであり、コントロール配列データが、該個体の別の組織、またはデータベースから得られた配列データである。１つの実施形態では、対象配列データが、細胞試料から得られた配列データであり、コントロール配列データが、別の細胞、またはデータベースから得られた配列データである。

本発明の方法は、全長配列の情報を必要としないため、データベース等において例えば全長配列が公知でない場合にも用いることができ、例えば、対象配列データおよびコントロール配列データは、生物のゲノムに由来する塩基配列データである場合、前記ゲノムの配列が不明であってもよい。

本発明の方法によって検出できる多型としては、置換、挿入、欠失、コピー数変異（Ｃｏｐｙ Ｎｕｍｂｅｒ Ｖａｒｉａｔｉｏｎ， ＣＮＶ）、ＳＴＲＰ（ｓｈｏｒｔ ｔａｎｄｅｍ ｒｅｐｅａｔ ｐｏｌｙｍｏｒｐｈｉｓｍ）、逆位および転座が挙げられる。１つの実施形態において、本発明の方法は、上記複数の多型の任意の組み合わせを同時に検出することも可能である。さらなる実施形態において、本発明の方法は、上記複数の多型の全てを同時に検出することも可能である。特に、多型が置換である場合には、本発明の方法は、非常に高い検出力を発揮することが可能である。

対象配列に対するリファレンス配列が存在する場合、本発明の方法は、前記対象配列に対するリファレンス配列における前記多型の位置を特定する工程をさらに含むことができる。例えば、対象配列データおよびコントロール配列データが、生物のゲノムに由来する塩基配列データである場合、多型のゲノム上の位置を特定する工程をさらに含むことができる。この位置の特定は、本発明の方法が、多型を周囲の配列と関連づけて検出する（例えば、ｘ長部分の多型がｋ−ｘ長の配列と関連付けられる）ことを可能にしているため、リファレンス配列に対して検索を行うことにより、簡便に行うことが可能である。

本発明の方法は、検出した多型について確認する工程をさらに含むことができる。確認は、例えば、検出された多型の部位について、リファレンス配列またはコントロール配列から作成したクエリ配列セットを用いて、対象配列データおよび／またはコントロール配列データとの比較を行うことによって行うことができる。クエリ配列セットは、リファレンス配列またはコントロール配列において前記多型に該当する部位の文字を異なる文字に置換した変異型クエリ配列セット、および／またはリファレンス配列またはコントロール配列において前記多型に該当する部位の文字を置換していない野生型クエリ配列セットを含み得る。

本発明の方法は、対象配列データおよびコントロール配列データが塩基配列データである場合、検出された多型の部位について、対象配列データおよび／またはコントロール配列データの相補鎖の配列データに対して、リファレンス配列またはコントロール配列から作成したクエリ配列セットとの比較を行い確認する工程をさらに含むことができる。本発明の方法は、対象配列データおよびコントロール配列データが塩基配列データである場合、検出された多型の部位について、対象配列データおよび／またはコントロール配列データの対立遺伝子の配列データに対して、リファレンス配列またはコントロール配列から作成したクエリ配列セットとの比較を行い確認する工程をさらに含むことができる。ここで、対立遺伝子の配列データとして、実際の遺伝子の存在の有無とは関係なく、野生型に対する変異型の塩基を有する配列データを用いることができる。

本発明の方法は、実験結果またはデータベースから対象配列データまたはコントロール配列データを取得する工程を含んでもよい。また、本発明の方法においては、必ずしも配列データそのものを取得する必要はなく、配列データのサブセット、および／または配列データもしくは配列データのサブセットにおける頻度分布のデータを取得して実行することも可能である。

１つの局面において、本発明は、本発明の多型を検出する方法をコンピュータに実施させるための方法を実装するプログラム、該プログラムを記録した記録媒体、およびこれを実現するためのシステムを提供する。ここで採用され得る任意の特徴は本明細書の多型を検出する方法の説明に記載される任意の特徴またはその組み合わせを採用することができる。

したがって、１つの実施形態において、対象配列データにおいてコントロール配列データに対する多型を検出する方法をコンピュータに実行させるためのプログラムであって、該方法は、
ａ）該対象配列データの長さｋの部分配列のサブセットをコンピュータに保存する工程であって、ここで、ｋは対象配列およびコントロール配列の全長以下の長さである、工程と、
ｂ）該対象配列データの長さｋのサブセットの各々の部分配列の出現頻度を算出する工程と、
ｃ）該コントロール配列データの長さｋの部分配列のサブセットにおける各々の部分配列の出現頻度をコンピュータに保存する工程と、
ｄ）対象配列とコントロール配列とを比較し、該出現頻度の分布の差異に基づいて、多型を検出する工程と
を包含する、プログラムが提供される。さらなる実施形態において、上記方法が、前記部分配列中の多型ではない部分の少なくとも一部を含む配列（前記部分配列全体であり得る）を、検出された前記多型の名称として表示する工程をさらに含む、プログラムが提供される。

別の実施形態において、対象配列データにおいてコントロール配列データに対する多型を検出する方法をコンピュータに実行させるためのプログラムを格納する記録媒体であって、該方法は、
ａ）該対象配列データの長さｋの部分配列のサブセットをコンピュータに保存する工程であって、ここで、ｋは対象配列およびコントロール配列の全長以下の長さである、工程と、
ｂ）該対象配列データの長さｋのサブセットの各々の部分配列の出現頻度を算出する工程と、
ｃ）該コントロール配列データの長さｋの部分配列のサブセットにおける各々の部分配列の出現頻度をコンピュータに保存する工程と、
ｄ）対象配列とコントロール配列とを比較し、該出現頻度の分布の差異に基づいて、多型を検出する工程と
を包含する、記録媒体が提供される。さらなる実施形態において、方法が、前記部分配列中の多型ではない部分の少なくとも一部を含む配列（前記部分配列全体であり得る）を、検出された前記多型の名称として表示する工程をさらに含む、記録媒体が提供される。

別の実施形態において、対象配列データにおいてコントロール配列データに対する多型を検出するためのシステムであって、該システムは、
該対象配列データおよび該コントロール配列データの長さｋの部分配列のサブセットの各々の部分配列の出現頻度を提供するように構成された配列データ処理部であって、ここで、ｋは対象配列およびコントロール配列の全長以下の長さである、配列データ処理部と、
対象配列とコントロール配列とを比較し、該出現頻度の分布の差異に基づいて、多型を検出する工程とを行うように構成された、配列データ計算部と
を備える、システムが提供される。さらなる実施形態において、前記部分配列中の多型ではない部分の少なくとも一部を含む配列前記部分配列全体であり得る）を、検出された前記多型の名称として表示する表示手段をさらに含む、システムが提供される。

（具体的な例）
本発明の例示的な実施形態は、以下のような工程による方法である。

１．配列データについて、長さｋの部分配列のサブセットを得る。例えば、塩基配列データを端から１塩基ずつずらしながら、ｋ長の配列のセットを得る。以下の例では、次世代シーケンサーで得られたイネ個体の配列データ（データ量はイネゲノムの４０倍、リード長は１００塩基）を用いて、ｋ長を２０塩基として得られた配列セットを用いて説明する。得られた配列セットに関して、比較対象のサンプル間で５’末端からｋ−１長の配列が同じであるデータの３’末端の塩基が比較対象の間で異なる場合、その塩基は多型である。

（上記例において、１行目から各々配列番号１９１〜２０１である。）

最初の行が次世代シーケンサーから得られた塩基配列であり、以下の行がｋ長の部分配列を示す。この場合ｋ長を２０塩基で配列セットを得ている。

２．得られた部分配列セットの出現頻度データを得る。すなわち、部分配列セットのうちの同じ配列に関して、出現回数で整理したデータを得る。

（上記例において、１行目から各々配列番号２０２〜２１１である。）

配列セットを降順に整列して出現回数を配列の右に表示している。この例では、ゲノムの４０倍の解析なので、出現回数が１回程度の配列はノイズであると判断することができる。４０〜５０回程度の配列は、ゲノム上のユニークな配列由来であり、８９回出現した配列は、ゲノム上２ヶ所存在していると考えられる。

３．部分配列中の長さｋ−ｘの配列部分が共通する配列ごとに、長さｘの部分について出現頻度の分布を算出する。この例では、得られたｋ塩基の部分配列（ｋ＝２０）の頻度データをもとに、最初の１９塩基（ｋ−ｘ、ｘ＝１）に対する最後の１塩基のＡＣＧＴそれぞれの塩基の出現頻度データに変換する。

（上記例において、１行目から各々配列番号２１２〜２２０である。）
（２０塩基の頻度データを最初の１９塩基と最後のＡＣＧＴの塩基の頻度一覧に変換する。）

４．部分配列中の長さｋ−ｘの配列部分が共通する配列における、長さｘの部分の出現頻度の分布の差異を比較する。例えば、コントロールと調べたい対象の配列データ（ターゲット）からの頻度データを１９ｍｅｒの配列でまとめた一覧を作成する。多型がない場合、最後の塩基の頻度は同じ塩基で最多となる。最後の塩基の頻度でコントロールと対象で塩基が異なる場合は、その塩基が多型である。

（上記例において、１行目から各々配列番号２１３、２１５、２１７、および２１８である。）

コントロールと調べたい対象の最後の塩基の出現頻度の一覧。この場合、ＡＡＡＡＧＡＴＣＴＡＴＧＡＧＣＡＣＴＣ（配列番号２１８）の次にはコントロールではＡのみであるが、対象ではＡとＧのヘテロザイガスであることがわかる。また、ホモザイガスの多型が生じる場合には、以下のように、最後の塩基として出現するものが異なり、検出することができる。

このように、ゲノム上の位置が不明でも、最初の１９塩基で表記される配列位置に続く塩基の多型を網羅することができる。ゲノムリファレンス配列が既知の場合、最初の１９塩基に対応するゲノム位置から多型の位置を決めることができる。順鎖、相補鎖ともに同じ位置に検出された場合は一塩基多型である。多型となる配列の末端部分を検出するため、挿入・欠失多型の検出も可能である。最初の１９塩基を多型の名称、最後の塩基を遺伝子型として表記するため、その多型を示す唯一の表記法として利用できる。あるいは、多型も含めてｋ−ｍｅｒ（例えば、２０塩基）全体で一つの名称とすることも可能である。例えばｋ＝５の場合には、ＡＣＧＴＡ型とＡＣＧＴＴ型といった表記が可能である。本発明の方法は、図１に示されるようなフローにしたがって、図１に示される工程を適宜採用することによって実行することが可能である。

（ｋ−ｍｅｒ）
１つの実施形態において、本発明の１つの特徴は、対象配列データにおけるコントロール配列データに対する多型の検出において、該対象配列データの長さｋの部分配列のサブセット、または該部分配列のサブセットの各々の部分配列の出現頻度を用いることである。ここで、ｋは対象配列およびコントロール配列のいずれか短いほうの全長以下の整数である。

長さｋの部分配列は、対象配列データ、コントロール配列データ、リファレンス配列データ等から切り出すことによって生成することができる。例えば、一定間隔でｋ長の配列を切り出すことによって生成することができ、１文字ごと、２文字ごと、３文字ごと、またはそれ以上の間隔で切り出して部分配列セットを生成することができる。あるいは、対象配列データ、コントロール配列データ、リファレンス配列データ等から、ランダムに位置を選択して切り出すことも可能である。全てのｋ長部分配列を生成することが望ましい場合には、１文字ごとに切り出し位置をずらしながらｋ長部分配列のサブセットを生成することができる。

長さｋは、対象配列、コントロール配列および／またはリファレンス配列における偶然同一を排除する長さであることが望ましい。偶然同一を排除することによって、異なる配列が対象配列の別の箇所に偶然含まれていたものをコントロール配列との差異として検出する可能性を低減し、より正確に解析することができる。ｋ長のある配列と、偶然に同一の配列が、対象配列、コントロール配列および／またはリファレンス配列に出現する期待値を１未満にすることが望ましい。

一般的に、ｖ：配列において各変数が取り得る値の種類、Ｌ：解析対象とする配列（対象配列、コントロール配列および／またはリファレンス配列）の全長（含まれる変数の数）として、ｖ＾ｋ＞Ｌとなる場合に、ｋが偶然同一を排除する長さであると考えられる。あるいは、配列全長が１０＾Ｌのようなオーダーで表される場合、両辺の対数をとり、ｋ＞Ｌ／ｌｏｇ（ｖ）を満たす場合に、ｋが偶然同一を排除する長さであると考えられる。

例えば、配列データが塩基配列データである場合、４つの文字が変数の値として考えられるため、ｖ＝４である。例えば、ヒトゲノムは３．１×１０の９乗のサイズであり、仮に１０の９乗長のアトランダムな塩基配列があった場合、９／ｌｏｇ（４）≒１５塩基が偶然一致を排除できる最小ｋ長と考えられる。ゲノムサイズが異なる生物に対しても、例えば、１０／ｌｏｇ（４）の場合はｋ＝１７で偶然一致を排除できると考えられる。長いほど偶然一致の配列が生じる可能性を減少させることができるが、データサイズがその分大きくなる。

（頻度）
１つの実施形態において、本発明の方法は、対象配列データおよび／またはコントロール配列データの長さｋの部分配列のサブセットの各々の部分配列の出現頻度を提供する工程を含む。

長さｋの部分配列のサブセットの各々の部分配列の出現頻度は、ユニークなｋ長の部分配列データに対する頻度がペアになったデータとして（例えば、１列目に部分配列、２列目にその頻度が提供される行列データとして）提供されてもよい。部分配列の出現頻度は、配列とその頻度とが関連付けられているものであれば、その形式は限定されるものではなく、対応する識別子を有する配列データと頻度データとを別個に出力することも可能である。また、部分配列の出現頻度は、既に存在するデータを取得することによって提供することも可能である。

長さｋの部分配列のサブセットの各々の部分配列の出現頻度を得る工程は、計算機を使用して簡便に行うことができ、例えば、Ｕｎｉｘでは以下：

のように実装することができる。

多型が存在している部位がゲノム上でユニークな配列部位であり、例えば１塩基置換が存在している場合、その置換を含むｋ−ｍｅｒ部分配列は対象配列データのサブセットには存在するが、コントロール配列データのサブセットには存在しないと考えられる。当該置換変異を含むｋ−ｍｅｒではすべてこの結果が得られると考えられ、結果的に当該置換変異部位を含む２×ｋ−１個のｋ−ｍｅｒ部分配列で出現頻度の差異が観察され、当該多型を検出することができる。

一部の実施形態において、本発明の方法において、予測される出現頻度と比較して、出現頻度が少ない部分配列をノイズとすることが可能である。これにより、機械的に生じた差異と、実際に生じている多型とを識別して検出することが可能である。例えば、前記対象配列データおよび前記コントロール配列データが塩基配列データである場合には、対象配列データ量、配列データにおける予想されるエラー率（例えば、マニュアル・機器の公称値）、対象配列全長サイズ（例えば、ゲノムサイズ）等を考慮して、予想されるノイズのカウント程度、または予想されるノイズのカウント以下の出現頻度の部分配列をノイズとすることができる。１つの予測としては、生じるノイズの頻度の期待値は、ｉｎｔ（ゲノムを何倍読んだか×（１−精度）＋１）」となる。

理論的には、シーケンサーの精度が９９％でゲノムの１００倍読みした場合は、１塩基のノイズが入ることになる。通常の解析はゲノムの４０倍読み程度であるため、ノイズは１塩基以下と考えられるが、カウントは整数値となるため、リードエラーは１とカウントされる。すなわち、２塩基以上のカウントはノイズ以外の何らかの要因がある可能性が高いと考えられる。ただし、低い確率だが、同一塩基に２回以上ノイズが入る可能性はある。予想されるノイズのカウントは、例えば、ｉｎｔ（ゲノムを何倍読んだか×（１−精度）＋１）である。この場合の精度は１００％の場合１、９９％の場合０．９９として計算できる。なお、ｉｎｔ（Ｘ）は、Ｘ未満の最大の整数を返す関数である。

エラー率（精度）が完全に０でない場合には、例えば０．００００１でも１０万か所のうち１ヶ所は１のカウントが生じると考えられる。そのため、この場合、結局１はノイズである可能性が高いと考えられるが、２はノイズにしては高すぎると言える。そのため、ｉｎｔ（ゲノムを何倍読んだか×（１−精度）＋１）で繰り上げて整数にすることにより、予測値を算出することが可能である。１つの実施形態では、カウントが１となる部分配列をノイズと判定する。例えば、対象配列データおよび前記コントロール配列データが塩基配列データである場合に、［（対象配列データ量）×（１−精度）］／（対象ゲノムサイズ）＋１未満の出現頻度の部分配列をノイズとすることにより、予想されるノイズのカウント以下の出現頻度の部分配列をノイズとして排除することができる。

頻度の算出から、コピー数変異（ＣＮＶ）の検出も可能である。例えば、コントロール配列と比較して約２倍以上の頻度が連続して検出された場合にＣＮＶと判定することができる。頻度の基準は、約２．５倍以上などと厳しくすることによって、誤検出を減少させることもできる。

あるいは、本発明の方法において、予測される出現頻度と比較して、出現頻度が多い部分配列を、リピート配列部位として除外することができる。例えば、対象配列データの対象配列全長のカバレッジ（カバー率）の２倍を超えるような部分配列を解析から除外することができる。

例えば、配列データ量の配列全長に対するカバレッジと比較し、同程度の出現頻度の部分配列は、対象配列の全長におけるユニークな配列に由来している配列と考えられる。さらに、配列データ量の配列全長に対するカバレッジと比較し、２倍程度の出現頻度の部分配列は、対象配列の全長において２箇所存在する配列に由来している配列と考えられる。

ｋが偶然同一を排除している場合、それを上回る場合には、部分配列が、リピート配列部位に由来しているものと判断することができる。例えば、ゲノムの４０倍のカバレッジのシーケンシングデータにおいて、出現頻度が４０〜５０回程度の配列は、ゲノム上のユニークな配列由来であり、８９回出現した配列は、ゲノム上２ヶ所存在していると考えられる。ｋ長配列のサブセットおよび／または各部分配列の出現頻度データの作成は、一例としては、図３に例示されるような工程を採用して行うことができる。

（ｋ−ｘ）
１つの実施形態では、本発明の方法は、ｋ長の部分配列中の長さｋ−ｘの配列部分が共通する配列ごとに、長さｘの部分について出現頻度の分布を算出することを特徴とする。ｘはｋ未満の正の整数である。長さｋの全長ではなく、ｋ−ｍｅｒの一部（ｘ長の部分）の文字の差異を見ることによって、計算量を顕著に減少させることが可能である。長さｘは、限定されるものではないが、好ましくは１〜２であり、より好ましくは１である。

一部の実施形態において、本発明の方法は、前記部分配列中の長さｋ−ｘの配列部分を、ユニークな配列ごとにグルーピングする工程を含む。これには、例えば、前記長さｋ−ｘの配列部分をソートする工程（例えば、前記長さｋ−ｘの配列部分を文字列としてソートする工程）が含まれ得る。

部分配列中の長さｋ−ｘの配列部分が共通する配列ごとの、長さｘの部分についての出現頻度の分布は、長さｋの部分配列の出現頻度から算出できる。配列の文字の種数をｖ（配列において各変数が取り得る値の種類）とした場合、長さｋの部分配列において、ｋ−ｘの配列部分が共通する配列が、ｋ−ｘの配列部分が共通する配列ごとにｖ＾ｘ種生じる。例えば、配列が塩基配列であり、ｘ＝１とした場合、長さｋの配列のセットの中には、ｋ−１の配列部分が共通する配列ごとに、ｘに対応する部分がＡ、Ｃ、Ｇ、Ｔである４種の配列が存在している。長さｋ−ｘの配列部分が共通する配列ごとの、長さｘの部分についての出現頻度は、それぞれに対応する長さｋの部分配列の頻度データに対応する。

１つの実施形態では、前記長さｘの部分が、前記部分配列の末端に存在する。例えば、前記対象配列データおよび前記コントロール配列データが塩基配列データである場合、前記長さｘの部分は、前記部分配列の３’末端または５’末端であり得る。長さｘの部分を部分配列の末端にとることは、比較処理の高速化・簡便化にとって望ましい。このように、ｋ長の配列のうち末端部で異なる部分（ｘ長の部分に該当）の頻度を検出する場合、対象となる配列の多型部位の「エッジ」（置換の場合はその位置そのものに該当するが、挿入・欠失変異の場合はその縁（エッジ）に該当する）を検出しているというように理解することができるため、本明細書において、「Ｐｏｌｙｍｏｒｐｈｉｃ Ｅｄｇｅ Ｄｅｔｅｃｔｉｏｎ（ＰＥＤ）」または「ｅｄｇｅ ｄｅｔｅｃｔｉｏｎ」と呼称される場合がある。

（比較・多型の検出）
出現頻度の分布の差異の比較により、例えば、以下のような多型の検出が可能である。

多型が存在している部位がゲノム上でユニークな配列部位であり、例えば置換が存在している場合、その置換を含むｋ−ｍｅｒ部分配列は対象配列データのサブセットには存在するが、コントロール配列データのサブセットには存在しないと考えられる。当該置換変異を含むｋ−ｍｅｒではすべてこの結果が得られると考えられ、結果的に当該置換変異部位を含む２×ｋ−１個のｋ−ｍｅｒ部分配列で出現頻度の差異が観察され、当該多型を検出することができる。そのような解析は、一例としては、図６に示されるような工程を採用して行うことが可能である。

１つの実施形態では、前記コントロール配列データのサブセットと前記対象配列データのサブセットとの間で、前記長さｘの部分の配列の出現頻度が異なる場合、該長さｘの部分の配列を、対象配列データにおけるコントロール配列データに対する多型として検出する。１つの実施形態では、前記コントロール配列データのサブセットと前記対象配列データのサブセットとの間で、前記長さｘの部分の配列で最も高頻度のものが異なっている長さｋ−ｘの配列部分が存在する場合、該長さｘの部分の配列を、対象配列データにおける多型として検出する。

１つの実施形態では、前記対象配列データおよび前記コントロール配列データが塩基配列データであり、前記対象配列データのサブセットにおける前記長さｘの部分の配列で、前記コントロール配列データのサブセットにおける最も高頻度のものと同一の長さｘの部分の配列がノイズ以下のカウントしか存在しない長さｋ−ｘの配列部分が存在する場合、該長さｘの部分の配列を、対象配列データにおけるホモ多型として検出する。１つの実施形態では、前記対象配列データおよび前記コントロール配列データが塩基配列データであり、対象サブセットにおける前記長さｘの部分の配列で、コントロール配列データのサブセットにおける最も高頻度のものと同一の長さｘの部分の配列が存在し、かつ、コントロール配列データのサブセットにおける最も高頻度のものと異なる長さｘの部分の配列が存在する長さｋ−ｘの配列部分が存在する場合、該長さｘの部分の配列を、対象配列データにおけるヘテロ多型として検出する。ｋ長配列のセットの頻度データの比較は、一例としては、図４に示されるような工程によって行うことが可能である。

このような比較は、１つの例としては、ｋ−１長配列と最後の塩基の頻度のファイル、controlとtargetを処理する場合、

のコマンドでｋ−１配列とcontrolとtargetを１行にまとめた頻度を出力することによって行うことができる。この出力結果の各行を調べる条件としては、限定されるものではないが、control、targetの両方でカウントが１塩基以下の塩基が２個以上存在し、controlあるいはtargetで１０以上のカウントを示した塩基に対応するtargetあるいはcontrolの塩基のカウントが１以下である事例が１ないし２回ある場合、多型の境界を検出したとすることができる。

すなわち、コントロール配列および／または対象配列の部分配列のサブセットにおける第１のカウントが第１の閾値を超えており、かつ、コントロール配列および／または対象配列の部分配列のサブセットにおける対応する第２のカウント（すなわち、第１のカウントがコントロール配列のものである場合、第２のカウントは対象配列のものであり、第１のカウントが対象配列のものである場合、第２のカウントはコントロール配列のものである）が第２の閾値を下回る場合、多型の境界を検出したとすることができる。

第１の閾値は、配列データのカバレッジによって変動するが、例えば、１０〜５０の範囲で設定することができる。第１の閾値は、例えば、１０〜４０、１０〜３０、１０〜２０、または１０〜１５の範囲で設定することができる。配列データのカバレッジが大きくなる場合には、第１の閾値も大きく設定することができ、例えば、ヒトゲノム解析では２００倍読み程度のデータも存在するが、この場合は、第１の閾値として２００を使用することができる。カバレッジを考慮して実際にその配列が存在しているといえるカウントを第１の閾値としてよく、例えば、カバレッジの約１００％、約９０％、約８０％、約７０％、約６０％、約４０％、約３０％、または約２０％等の値を用いることができる。

第２の閾値も同様に配列データのカバレッジによって変動するが、１〜７の範囲で設定することができる。第２の閾値は、例えば、１、２、３、４、５、６または７であり得る。

（配列）
本発明の対象配列、コントロール配列および／またはリファレンス配列としては、多型が生じ得る任意の配列を用いることができる。なお、コントロール配列として、リファレンス配列を用いることが可能である。代表的な実施形態では、対象配列、コントロール配列および／またはリファレンス配列は、生物学的配列であり、例えば、塩基配列（ＤＮＡ、ＲＮＡ、およびそれらのアナログ等の配列が包含される）、アミノ酸配列、または糖鎖配列等である。生物学的配列の例としては、例えば、ゲノム配列、染色体配列、遺伝子配列、プラスミド配列、エクソン配列、タンパク質配列等が挙げられる。

対象配列データおよびコントロール配列データは、限定されるものではないが、多型を検出する上では、一定の共通性を持つ配列についての配列データであることが望ましい。しかしながら、配列の取得方法については各々同一でも異なっていてもよく、シーケンシングによって得られたデータ間での比較を行うことも、データベース等から得られたデータ間での比較を行うことも、シーケンシングによって得られたデータとデータベース等から得られたデータとの間での比較を行うことも可能である。

１つの実施形態では、対象配列データが、個体から得られた配列データであり、コントロール配列データが、該個体と同種の別の個体、またはデータベースから得られた配列データである。１つの実施形態では、対象配列データが、個体の組織試料から得られた配列データであり、コントロール配列データが、該個体の別の組織、またはデータベースから得られた配列データである。１つの実施形態では、対象配列データが、細胞試料から得られた配列データであり、コントロール配列データが、別の細胞、またはデータベースから得られた配列データである。

本発明の方法は、全長配列の情報を必要としないため、データベース等において例えば全長配列が公知でない場合にも用いることができ、例えば、対象配列データおよびコントロール配列データは、生物のゲノムに由来する塩基配列データである場合、前記ゲノムの配列が不明であってもよい。

例えば、次世代シーケンサーのリードデータで対象（ターゲット）とコントロール間の配列の直接比較による多型検出は既存の技術では不可能であった。リファレンス配列がある場合は、多型をゲノム上にマップすることができるが、リファレンスが存在しない場合（例えば、リファレンスゲノム配列が作成されていない生物）で、リファレンス情報を用いることなく対象とコントロール間での多型を検出できるというのは画期的である。リファレンスのない生物でのＦ２分離集団で表現型の分離と連鎖する多型を検出すれば、ゲノム上の位置が不明であったとしても、表現型に対応するＤＮＡマーカーを得ることができ、応用範囲は非常に広いと考えられる。実際、育種にはゲノム上の位置情報も重要ではあるが、位置情報が不明でも優良形質にリンクしているＤＮＡマーカーでの選抜ができれば、選抜育種への利用が十分に可能である。

同一個体（例えば、ヒト）の正常組織と変異を含み得る組織（例えば、癌組織）の配列を直接比較して多型を検出することも可能で、いったんリファレンスゲノムにマップしてから差を調べる方法に比べて、多型の捕捉率、捕捉精度ともに非常に高くなる。

１つの実施形態では、本発明の方法で用いる対象配列データおよび／またはコントロール配列データは、シーケンシングによって得られた塩基配列データである。シーケンシングの手法としては、サンガー法、マクサム・ギルバード法、単一分子リアルタイムシーケンシング（例えば、Ｐａｃｉｆｉｃ Ｂｉｏｓｃｉｅｎｃｅｓ、Ｍｅｎｌｏ Ｐａｒｋ、Ｃａｌｉｆｏｒｎｉａ）、イオン半導体シーケンシング（例えば、Ｉｏｎ Ｔｏｒｒｅｎｔ、Ｓｏｕｔｈ Ｓａｎ Ｆｒａｎｃｉｓｃｏ、Ｃａｌｉｆｏｒｎｉａ）、パイロシーケンシング（例えば、４５４、Ｂｒａｎｆｏｒｄ、Ｃｏｎｎｅｃｔｉｃｕｔ）、ライゲーションによるシーケンシング（例えば、Ｌｉｆｅ Ｔｅｃｈｎｏｌｏｇｉｅｓ、Ｃａｒｌｓｂａｄ、ＣａｌｉｆｏｒｎｉａのＳＯＬｉＤシーケンシング）、合成および可逆性ターミネーターによるシーケンシング（例えば、Ｉｌｌｕｍｉｎａ、Ｓａｎ Ｄｉｅｇｏ、Ｃａｌｉｆｏｒｎｉａ）、透過型電子顕微鏡法などの核酸イメージング技術、ナノポアシーケンシングなどがある。

１つの実施形態では、本発明の方法で用いる対象配列データおよび／またはコントロール配列データは、次世代シーケンシングによって得られた配列データであり得る。次世代シーケンシングとしては、シーケンシングバイシンセシス、パイロシーケンシング、ライゲーションによるシーケンシング、イオン半導体シーケンシング、ナノポアシーケンシング等が挙げられる。次世代シーケンシングデータからの多型の検出においては、リファレンスへのマッピングやアセンブリによって精度が制限されていたため、本発明の方法を用いた場合に大きな利益が得られると考えられる。

１つの実施形態では、本発明の方法で用いる対象配列データおよび／またはコントロール配列データは、ジニトロフェニル化法、ヒドラジン分解法、カルボキシペプチダーゼ法、エドマン法もしくはそれらの方法を自動化する装置（ペプチドシーケンサーあるいはプロテインシーケンサー）を用いる方法、質量分析（例えば、タンデム質量分析計（ＭＳ／ＭＳ））を用いた方法（例えば、シーケンスタグ法）等から得られたアミノ酸配列データである。

本発明の対象配列データおよび／またはコントロール配列データの由来となる生物種としては、生物学的配列を有するものである以上は何ら制限されない。一部を例示すると、動物としては、ヒトもしくは非ヒト哺乳動物（例えば、マウス、ラット、ウサギ、ヒツジ、ブタ、ウシ、ウマ、ネコ、イヌ、サル、チンパンジー）、鳥類、爬虫類、両生類、魚類等の脊椎動物、無脊椎動物、例えば、昆虫、線形動物などを挙げることができる。植物としては、イネ、コムギ、トウモロコシ、ジャガイモ、オオムギ、サツマイモ、ソバ、シロイヌナズナ、ミヤコグサ、トマト、キュウリ、キャベツ、白菜、ナス、サトウキビ、ソルガム、リンゴ、ミカン、バナナ、桃、ポプラ、松、杉、被子植物、裸子植物、シダ、コケ、藻類などを挙げることができる。その他、真菌、細菌、ウイルス等でもよい。

さらに、これらの生物の一部分、例えば、組織、細胞等に由来する対象配列データおよび／またはコントロール配列データを解析し、多型を検出することも可能である。

（変異）
本発明の方法は、例えば、置換、挿入、欠失、コピー数変異、ＳＴＲＰ（ｓｈｏｒｔ ｔａｎｄｅｍ ｒｅｐｅａｔ ｐｏｌｙｍｏｒｐｈｉｓｍ）、逆位または転座等の多型の検出に用いることができる。変異のエッジの部分が検出されるため、挿入・欠失の結果、長さｘの配列に違いがあればそのエッジ部分を検出することができる。ｋ−ｍｅｒ内に収まりきる場合であれば、ＳＴＲＰ（ｓｈｏｒｔ ｔａｎｄｅｍ ｒｅｐｅａｔ ｐｏｌｙｍｏｒｐｈｉｓｍ）を検出することも可能である。ＳＴＲ（ｓｈｏｒｔ ｔａｎｄｅｍ ｒｅｐｅａｔ）は、マイクロサテライトとも称され、２〜７塩基からなる配列が２〜数十回反復するもので、この回数に多型が見られる。部分配列の出現頻度によって、コピー数多型（ＣＮＶ）を検出することもできる。エッジ検出という観点からは、逆位、転座のエッジも検出することが可能である。

特に、多型が置換である場合には、本発明の方法は、非常に高い検出力を発揮することが可能である。

（位置の特定）
対象配列に対するリファレンス配列が存在する場合、本発明の方法は、対象配列に対するリファレンス配列における前記多型の位置を特定する工程をさらに含むことができる。例えば、対象配列データおよびコントロール配列データが、生物のゲノムに由来する塩基配列データである場合、多型のゲノム上の位置を特定する工程をさらに含むことができる。この位置の特定は、本発明の方法が、多型を周囲の配列と関連づけて検出する（例えば、ｘ長部分の多型がｋ−ｘ長の配列と関連付けられる）ことを可能にしているため、リファレンス配列に対して検索を行うことにより、簡便に行うことが可能である。

リファレンス配列に対する検索は、一例としては、図５に示されるような工程によって二分検索用リファレンスゲノム配列データを作成し、その後二分検索によって多型境界塩基のマッピングによって検索を行うことができる。

他の方法として、ｕｎｉｘのｊｏｉｎコマンドによって対象配列データ中の部分配列の、リファレンス配列における位置および向きを出力することによってマッピングを行う方法を用いることができる。より詳細には、対象配列のコントロール配列上の位置を決定する方法であって、ａ）コントロール配列中の複数のｋ長の部分配列について、配列およびコントロール配列中の位置と向きを出力する工程と、ｂ）対象配列中の複数のｋ長の部分配列について、配列および対象配列中の位置を出力する工程と、ｃ）ａ）およびｂ）で得られた配列を比較し、同一の部分配列に対応するコントロール配列中の位置と対象配列中の位置とを対応付ける工程とを含み、ここで、ｋは、対象配列の長さを超えない長さである、方法を使用することができる。当該方法については、本出願人により本出願と同日に出願された「挿入・欠失・逆位・転座・置換検出法」との名称の出願（整理番号ＮＧ０１２ＰＣＴ／Ｆ５−１８ＰＣＴ０７５）を参照することができる。

（確認（ｖｅｒｉｆｙ））
本発明の方法は、検出した多型について確認する工程をさらに含むことができる。確認は、例えば、検出された多型の部位について、リファレンス配列またはコントロール配列から作成したクエリ配列セットを用いて、対象配列データおよび／またはコントロール配列データとの比較を行うことによって行うことができる。クエリ配列セットは、リファレンス配列またはコントロール配列において前記多型に該当する部位の文字を異なる文字に置換した変異型クエリ配列セット、および／またはリファレンス配列またはコントロール配列において前記多型に該当する部位の文字を置換していない野生型クエリ配列セットを含み得る。

本発明の方法は、対象配列データおよびコントロール配列データが塩基配列データである場合、検出された多型の部位について、対象配列データおよび／またはコントロール配列データの相補鎖の配列データに対して、リファレンス配列またはコントロール配列から作成したクエリ配列セットとの比較を行い確認する工程をさらに含むことができる。本発明の方法は、対象配列データおよびコントロール配列データが塩基配列データである場合、検出された多型の部位について、対象配列データおよび／またはコントロール配列データの対立遺伝子の配列データに対して、リファレンス配列またはコントロール配列から作成したクエリ配列セットとの比較を行い確認する工程をさらに含むことができる。確認する工程は、一例としては、図６に示されるフローに従って、図６に示される工程を適宜採用して行うことができる。ここで、対立遺伝子の配列データとして、実際の遺伝子の存在の有無とは関係なく、野生型に対する変異型の塩基を有する配列データを用いることができる。

以下、コントロール配列が野生型（つまりリファレンスゲノム配列とほぼ同一である）から得られた場合、または、コントロール配列が、リファレンスゲノム配列から対象配列と同一長で作成された配列である場合について例示する。

次世代シーケンサーで読まれたショートリードの塩基配列長がＬの場合、リファレンスゲノム配列の対象となる多型塩基位置を起点として、Ｌ−１塩基前の位置からＬ−１塩基後の位置までの２Ｌ−１塩基長の配列を得て、多型塩基位置を推定された多型塩基に置換した置換配列とコントロールの非置換配列を作成する。置換配列、および、非置換配列それぞれを１塩基ずつずらしながらＬ長のクエリ配列セットを作成する（例えば、図１４に例示される）。個々のクエリ配列を用いて対象配列とコントロール配列に対して完全一致するカウントを取得する。ホモ型変異の場合は変異を導入したクエリ配列では対象配列のカウントが大部分になり、ヘテロ型変異の場合は、対象配列とコントロール配列への一致カウントが概ね半分ずつとなると考えられる。非置換配列をクエリにした場合は、基本的には大部分コントロール配列にヒットする。非置換配列が対象配列にヒットする場合は、多型ではないと判断し、除外することが可能である。

コントロール配列、対象配列の検索のため、コントロール配列、対象配列、それぞれ相補鎖の配列も合わせたのち辞書順にソートして同一配列は一つにまとめたデータセットを用いて、二分探索法によりクエリ配列を検索する。実施形態において、Ｆａｓｔｑファイルからの塩基配列データをソートしてユニークなデータセットにする時点で、つまり最初の段階で、各リードの相補鎖も一緒にソートしてユニーク処理することも可能である。

ｋ−ｍｅｒを切り出す前にsortとuniqの処理をすることは、シーケンス反応時にＰＣＲ増幅のステップが入る場合があり（入らないキットも存在する）、同一の配列が複数回リードデータの中に出てくる場合があることに対処する上で好ましい場合がある。そのまま処理するとｋ−ｍｅｒの分布が歪む可能性がある。同一配列ではあるがＮが含まれるようなリードも別物と認識されて歪む原因になり得るため、Ｎを含まない配列で、その配列とその相補鎖配列をsortしてuniq処理した配列からｋ−ｍｅｒを得ることが可能である。

歪むことを許容する場合、もとの配列の長さは揃っている必要はなく、サンガー法で得られたような長さが一定でない配列データでもｋ−ｍｅｒにして多型検出、マッピングまで行うことが可能である。

確認工程における、sort_uniq配列を二分検索して、変異型と野生型のリード数を調べる工程では、配列データの長さが揃っていることが好ましい。本発明者らの知見によれば、確認工程の前にｋ−１配列で多型をマップした段階で、順鎖と相補鎖の両方で多型が検出された場合は、ほぼ間違いなくＳＮＰであることがわかっている。順鎖、相補鎖のどちらかがリピート領域にあって片側しか検出できないものでも実際にＳＮＰである場合があるが、このような場合、sort_uniq配列を二分検索して確認することで、実際のＳＮＰかどうかの判断をすることができる場合があり得る。二分検索による確認を行って捕捉率を上げるためには、スタートのショートリード（次世代シーケンサーで得られた塩基配列データ）の長さが揃っていることが好ましい場合がある。対象とコントロールの配列長が同一である必要はなく、対象配列データとコントロール配列データそれぞれの中で長さが一定であれば、好適に二分検索によって確認を行うことができる。

あるいは、以下のように確認を行うことができる。リファレンス配列から対象配列のＬ長で変異部位を含むようにして部分配列を切り出し、変異に置換したセットと、置換しないセットを作り、位置関係、変異の有無等の記載と一緒にソートして出力する。このデータとソートした対象配列を、ｕｎｉｘコマンドのｊｏｉｎで処理（または適切な等価な処理）して、対象配列中に含まれる、野生型と変異型の配列を選び、配列数を変異部位ごとに調べる。選んだ配列を、ｓｏｒｔした後、ｕｎｉｑ−ｃのコマンドで配列数を数えることができる。同様の操作を、コントロール配列（Ｌ’長）に対しても行う。対象個体に対するコントロール個体がある場合は、この個体から得られたリード配列をコントロール配列として用いることができる。コントロール個体がない場合は、リファレンス配列から、Ｌ長で切り出して作った配列をコントロール配列として使うことが可能である。対象配列およびコントロール配列で長さが違う場合は、それぞれの長さに対応する変異型と野生型のデータセットを作り、対応する個数を調べることができる。当該方法については、本出願人により本出願と同日に出願された「挿入・欠失・逆位・転座・置換検出法」との名称の出願（整理番号ＮＧ０１２ＰＣＴ／Ｆ５−１８ＰＣＴ０７５）を参照することができる。

例えば、イルミナ社の次世代シーケンサー（例えば、ＨｉＳｅｑ）は同じ長さの配列を出力するため、対象配列データにおいて、特に長さをそろえる処理をする必要はない。この場合、対象配列データと同じ長さを有するクエリ配列のセットを作成することが可能であるため、直接二分検索を行うことができ、確実な結果を得る上で有利であり得る。

長さにばらつきがあるショートリードによる配列データあるいはショートリードの集合ではない配列データ（例えば、サンガー法で得られた配列）でも、本発明において適用可能である。確認工程においては、長さを揃える処理（例えば、リード中のクオリティスコアが最大となるＬ長配列を選択して切り出す、一端からＬ長配列を切り出す等）の処理をしたデータを用いて二分検索を行うことが可能である。あるいは、配列データ（リード）をＢＬＡＳＴのターゲット配列（データベース）にして、配列データとは長さの異なるクエリ配列を検索して数を数えることも可能である。

コンピュータで計算させる場合、sort_uniq配列や二分探索用のリファレンス配列に対して二分探索で完全一致の配列を探索する場合に、ファイルサイズが大きくなる場合が多いため、高速のシステムを利用するか、適切に高速化を図ることが好ましい。高速化の手法としては、すべてオンメモリで計算を行う、ＳＳＤ等高速ディスクにファイルを置く等、ハード的な手法が存在する。ソフト的には、単純にソートされたファイルに対して二分検索を行うのではなく、検索対象のファイルをＢｕｒｒｏｗｓ−Ｗｈｅｅｌｅｒ変換してより高速化を行うこともできる。

（プログラム、記録媒体およびシステム）
１つの局面において、本発明は、本発明の多型を検出する方法をコンピュータに実施させるための方法を実装するプログラム、該プログラムを記録した記録媒体、およびこれを実現するためのシステムを提供する。ここで採用され得る任意の特徴は本明細書の多型を検出する方法の説明に記載される任意の特徴またはその組み合わせを採用することができる。

したがって、１つの実施形態において、対象配列データにおいてコントロール配列データに対する多型を検出する方法をコンピュータに実行させるためのプログラムであって、該方法は、
ａ）該対象配列データの長さｋの部分配列のサブセットをコンピュータに保存する工程であって、ここで、ｋは対象配列およびコントロール配列の全長以下の長さである、工程と、
ｂ）該対象配列データの長さｋのサブセットの各々の部分配列の出現頻度を算出する工程と、
ｃ）該コントロール配列データの長さｋの部分配列のサブセットにおける各々の部分配列の出現頻度をコンピュータに保存する工程と、
ｄ）対象配列とコントロール配列とを比較し、該出現頻度の分布の差異に基づいて、多型を検出する工程と
を包含する、プログラムが提供される。さらなる実施形態において、方法が、前記部分配列中の多型ではない部分の少なくとも一部を含む配列（前記部分配列全体であり得る）を、検出された前記多型の名称として表示する工程をさらに含む、プログラムが提供される。プログラムはどのような言語で記述されてもよい。

別の実施形態において、対象配列データにおいてコントロール配列データに対する多型を検出する方法をコンピュータに実行させるためのプログラムを格納する記録媒体であって、該方法は、
ａ）該対象配列データの長さｋの部分配列のサブセットをコンピュータに保存する工程であって、ここで、ｋは対象配列およびコントロール配列の全長以下の長さである、工程と、
ｂ）該対象配列データの長さｋのサブセットの各々の部分配列の出現頻度を算出する工程と、
ｃ）該コントロール配列データの長さｋの部分配列のサブセットにおける各々の部分配列の出現頻度をコンピュータに保存する工程と、
ｄ）対象配列とコントロール配列とを比較し、該出現頻度の分布の差異に基づいて、多型を検出する工程と
を包含する、記録媒体が提供される。さらなる実施形態において、方法が、前記部分配列中の多型ではない部分の少なくとも一部を含む配列（前記部分配列全体であり得る）を、検出された前記多型の名称として表示する工程をさらに含む、記録媒体が提供される。プログラムはどのような言語で記述されてもよい。１つの実施形態では、記録媒体は、内部に格納され得るＲＯＭやＨＤＤ、磁気ディスク、ＵＳＢメモリ等のフラッシュメモリなどの外部記憶装置でありうる。

別の実施形態において、対象配列データにおいてコントロール配列データに対する多型を検出するためのシステムであって、該システムは、該対象配列データおよび該コントロール配列データの長さｋの部分配列のサブセットの各々の部分配列の出現頻度を提供するように構成された配列データ処理部であって、ここで、ｋは対象配列およびコントロール配列の全長以下の長さである、配列データ処理部と、対象配列とコントロール配列とを比較し、該出現頻度の分布の差異に基づいて、多型を検出する工程とを行うように構成された、配列データ計算部とを備える、システムが提供される。さらなる実施形態において、前記部分配列中の多型ではない部分の少なくとも一部を含む配列（前記部分配列全体であり得る。）を、検出された前記多型の名称として表示する表示手段をさらに含む、システムが提供される。

次に、図１５Ａの機能ブロック図を参照して、本発明のシステム１の構成を説明する。なお、本図においては、単一のシステムで実現した場合を示しているが、複数のシステムで実現される場合も本発明の範囲に包含されることが理解される。

本発明のシステム１０００は、コンピュータシステムに内蔵されたＣＰＵ１００１にシステムバス１０２０を介してＲＡＭ１００３、ＲＯＭやＨＤＤ、磁気ディスク、ＵＳＢメモリ等のフラッシュメモリなどの外部記憶装置１００５及び入出力インターフェース（Ｉ／Ｆ）１０２５が接続されて構成される。入出力Ｉ／Ｆ１０２５には、キーボードやマウスなどの入力装置１００９、ディスプレイなどの出力装置１００７、及びモデムなどの通信デバイス１０１１がそれぞれ接続されている。外部記憶装置１００５は、情報データベース格納部１０３０とプログラム格納部１０４０とを備えている。何れも、外部記憶装置１００５内に確保された一定の記憶領域である。

このようなハードウェア構成において、入力装置１００９を介して各種の指令（コマンド）が入力されることで、又は通信Ｉ／Ｆや通信デバイス１０１１等を介してコマンドを受信することで、この記憶装置１００５にインストールされたソフトウェアプログラムがＣＰＵ１００１によってＲＡＭ１００３上に呼び出されて展開され実行されることで、ＯＳ（オペレーションシステム）と協働して本発明の対象配列データにおいてコントロール配列データに対する多型を検出する方法の機能を奏するようになっている。もちろん、このような協働する場合以外の仕組みでも本発明を実装することは可能である。

本発明の実装において、対象配列データの長さｋの部分配列のサブセットをコンピュータに保存する工程であって、ここで、ｋは対象配列およびコントロール配列の全長以下の長さである、工程を行う際に、対象配列データおよび／または対象配列データの長さｋの部分配列のデータは、入力装置１００９を介して入力され、あるいは、通信Ｉ／Ｆや通信デバイス１０１１等を介して入力されるか、あるいは、データベース格納部１０３０に格納されたものであってもよい。次に、対象配列データの長さｋのサブセットの各々の部分配列の出現頻度を算出する工程は、プログラム格納部１０４０に格納されたプログラム、または、入力装置１００９を介して各種の指令（コマンド）が入力されることで、又は通信Ｉ／Ｆや通信デバイス１０１１等を介してコマンドを受信することで、この外部記憶装置１００５にインストールされたソフトウェアプログラムによって実行することができる。あるいは、あらかじめ算出されている出現頻度を、入力装置１００９を介して入力することができる。出現頻度データは、出力装置１００７を通じて出力されるかまたは情報データベース格納部１０３０等の外部記憶装置１００５に格納されてもよい。

次に、コントロール配列データの長さｋの部分配列のサブセットにおける各々の部分配列の出現頻度をコンピュータに保存する工程を行う差異、コントロール配列データ、コントロール配列データの長さｋの部分配列のサブセットのデータ、または部分配列の出現頻度のデータは、入力装置１００９を介して入力され、あるいは、通信Ｉ／Ｆや通信デバイス１０１１等を介して入力されるか、あるいは、データベース格納部１０３０に格納されたものであってもよく、プログラム格納部１０４０に格納されたプログラム、または、入力装置１００９を介して各種の指令（コマンド）が入力されることで、又は通信Ｉ／Ｆや通信デバイス１０１１等を介してコマンドを受信することで、この外部記憶装置１００５にインストールされたソフトウェアプログラムによってこれらのデータを処理してコントロール配列データの長さｋの部分配列のサブセットにおける各々の部分配列の出現頻度を提供してもよい。

対象配列とコントロール配列とを比較し、該出現頻度の分布の差異に基づいて、多型を検出する工程は、プログラム格納部１０４０に格納されたプログラム、または、入力装置１００９を介して各種の指令（コマンド）が入力されることで、又は通信Ｉ／Ｆや通信デバイス１０１１等を介してコマンドを受信することで、この外部記憶装置１００５にインストールされたソフトウェアプログラムによって実行することができる。

データベース格納部１０３０には、これらのデータや計算結果、もしくは通信デバイス１０１１等を介して取得した情報が随時書き込まれ、更新される。各入力配列セット中の各々の配列、参照データベースの各遺伝子情報ＩＤ等の情報を各マスタテーブルで管理することにより、蓄積対象となるサンプルに帰属する情報を、各マスタテーブルにおいて定義されたＩＤにより管理することが可能となる。

データベース格納部１０３０には、上記計算結果は、配列に関する情報、例えば、生物学的情報、生化学的情報、医学的情報、例えば疾患、障害、生体情報等の既知の情報と関連付けて格納されてもよい。このような関連付けは、ネットワーク（インターネット、イントラネット等）を通じて入手可能なデータをそのまままたはネットワークのリンクとしてなされてもよい。

また、プログラム格納部１０４０に格納されるコンピュータプログラムは、コンピュータを、上記した処理システム、例えば、配列データの提供、部分配列サブセットの提供、出現頻度データの算出、出現頻度データの比較、多型の検出、多型の確認などの処理を実施するシステムとして構成するものである。これらの各機能は、それぞれが独立したコンピュータプログラムやそのモジュール、ルーチンなどであり、上記ＣＰＵ１００１によって実行されることでコンピュータを各システムや装置として構成させるものである。なお、本発明の例示においては、それぞれのシステムにおける各機能が協働してそれぞれのシステムを構成しているものとするが、この処理のためのプログラムもまた、それぞれ外部記憶装置または通信デバイスまたは入力装置を介して提供されうる。

本発明がシステムとして構成される場合は、対象配列データおよび／またはコントロール配列データ、それらの長さｋのサブセットのデータ、ならびに／あるいはそれらの出現頻度データの提供は、まとめて配列データ処理部としてもよい。また、出現頻度の分布の比較および多型の検出は、配列データ計算部としてまとめてもよい。

また、図１５Ｂに示されるように、クラスター構造を有する計算システムによって本発明の方法を実装してもよい。１つの実施形態では、システムはクラスター構成であり、ヘッドとノードからなる。ノードは検索の高速化を図るため、主記憶装置にＳＳＤを用いることができる。１つの実施形態では、ヘッド１台に対して複数のノード（例えば１２台）で運用することができる。１つの実施形態では、計算システムはクラスター構造を持ち、主コンピュータ（クラスターヘッド）に大容量記憶装置（ＨＤＤ）を搭載して解析データおよび結果を保存する。クラスターヘッドより、分割したデータを各ノードに送り計算を実行し、結果をクラスターヘッドに集約する。クラスターヘッド、ノード共に、中央制御素子（ＣＰＵ）、メモリ（ＲＡＭ）を搭載し、通信インターフェース（ＮＩＣ）を介してデータの通信を行い得る。ノードには高速での検索処理をするため、ソリッドステートドライブ（ＳＳＤ）を主記憶装置とすることができる。各ノードに搭載されるＣＰＵ、ＲＡＭ、ＳＳＤ等は、他のノードと共有されてもよく、物理的に分離していてもよい。

（例示的計算フロー）
本発明は、例えば、以下のフローにおいて実施することができる。

１．次世代シーケンサーから得られたＦａｓｔｑ形式のファイルより、配列内にＮを含まない塩基配列を選んで、相補鎖の配列と合わせて、それぞれ出力する。

２．出力されたファイル（reads）を辞書順にソートし同一の配列は一つにまとめる。
ｕｎｉｘコマンドは以下の通り。

相同な配列がゲノム上に複数存在する部分に対しての解析を行う場合は、readsに対してsortのみを行い、uniq処理は行わない場合もあり得る。例えば、マルチコピー領域に対してのＣＮＶを検出する場合にはsortのみを行うことが有利であり得る。

３．ソートされた各塩基配列をそれぞれ５’末端から１塩基ずつずらしながら、ｋ−ｍｅｒ（ここでは２０塩基）の配列を対象の塩基配列の３’末端に到達するまで繰り返し出力（k-mer_file）する。

４．出力されたｋ−ｍｅｒの配列を辞書順にソートし、同一配列は一つにまとめて、出現回数を配列と共に表記したファイルを作成する。
ｕｎｉｘコマンドは以下の通り。

uniq -cコマンドなどのように、出力結果が頻度→配列の順になる場合、uniq-c等のコマンドの後に配列→頻度の順に出力されるフィルタープログラムを通してもよい。

５．k-mer_count_fileの各行のデータに対して、配列の５’末端よりｋ−１塩基の配列を得て、３’末端の塩基すなわちｋ番目の塩基をＡ、Ｃ、Ｇ、Ｔの出現回数として表記したデータに変換する。
ｋ−１ｍｅｒの配列 Ａの回数 Ｃの回数 Ｇの回数 Ｔの回数
という形式で出力される。

６．このようなデータを対象（target）と比較（control）のサンプルから得られた次世代シーケンサーの配列データからそれぞれ得る。

７．controlとtargetの５の方法で作成したデータをｋ−１ｍｅｒの配列でまとめたデータを作成する。
ｕｎｉｘコマンドは以下の通り

８．joinコマンドで得られたデータでcontrolとtargetで異なる塩基でカウントが示されるデータを選び出す。

この例では、controlでは、ＡＣＴＴＴＣＴＴＣＡＡＧＧＴＣＴＧＴＴ（配列番号２２５）に続く塩基はＧであるが、targetではＣである。つまり、ｋ−１ｍｅｒのユニークな識別子（名称）に続く塩基がＧ型、あるいはＣ型という表記で多型を表現することができる。それぞれの塩基に対応する数は、この多型が検出された独立のｒｅａｄの数である。

９．このｋ−１ｍｅｒの識別子をリファレンスゲノム上にマッピングすれば、容易に多型位置を決めることができる。

１０．本法は配列情報をリファレンスゲノムにマップする前にcontrolとtarget間の多型を検出するので、リファレンスゲノムが未知の生物種でも多型を検出できる。

１１．ゲノム上の場所が決められない場合は、多型の名称はｋ−１ｍｅｒの配列自身で表すことができ、遺伝子型はそれに続く多型塩基となる。このデータセットを用いて、表現型に対するアソシエーション解析が可能である。

１２．リファレンスゲノムへのマッピングには、リファレンスゲノムを３の方法と同様に各染色体を５’側から１塩基ずつずらしながらｋ−ｍｅｒで切り出し、ｋ−ｍｅｒ配列、染色体番号、位置、向きを同一行に表記したデータをｋ−ｍｅｒの配列で辞書順に整列したデータに対して、二分検索法によりｋ−１ｍｅｒの配列のゲノム上の位置を決めることができる。

１，２行目のように同一配列も複数行でそれぞれの位置がわかるので、対象配列がリピート領域に存在しても、検索では対応する候補領域が列挙されて出力される。

今回は、複数の位置情報が出力された場合は、位置不定として多型情報の出力から除外している。複数の位置のどれかという出力の仕方も可能である。

１３．リファレンスゲノムにマップされた多型の確認を行う。２で作成されたsort_uniq配列（ターゲット）と同様に、リファレンス配列も５’末端から２塩基ごとに２で作成された配列と同じ長さの配列を切り出し、相補鎖とともにsort,uniq処理をしたデータを作成する。このリファレンスとターゲットのsort_uniqデータが検索の対象となる。

次に、リファレンスゲノム配列から、多型位置を含む２で作成された配列と同じ長さの配列セット（検索を行うクエリセット）を切り出す。この配列セットは配列の５’末端から３’末端までのすべての位置で多型位置の塩基を含む配列セットである。配列セット数は配列長と同一となる。リファレンスの配列セット（リファレンスセット）に対して、多型位置の塩基を予想された多型塩基に置換した配列セット（ターゲットセット）を作成する。クエリ用配列セットの作成は、図１４に例示される。

１４．リファレンスセット、ターゲットセットをクエリにして、リファレンスゲノム、及び、２で作成されたsort_uniq配列をそれぞれ、検索（例えば、二分検索法またはjoinコマンドを用いた方法）し、それぞれのセットに対してsort_uniq配列が何個マッチするかを調べる。リファレンスセットではリファレンスのsort_uniq配列のみにヒットする。これに対して、ターゲットセットでは、ホモ型の変異の場合は、ターゲットのsort_uniqデータのみから検出され、ヘテロ型の変異の場合は、リファレンスと、ターゲットのsort_uniqデータの両方から検出されるはずである。このようにして、予想と一致する検索値を示した多型を抽出すると、以下のように結果を出力することが可能である。

（組み合わせ）
本明細書において、上述のとおり置換、コピー数多型、ＳＴＲＰ、挿入、欠失、逆位または転座を検出するのに有用な方法を記載しているが、かかるプロセスは、置換、挿入、欠失、逆位または転座を検出するのに有用な以下に記載するプロセスと組み合わせて行うことができる。例えば、組み合わせた方法は、図１８に示されるようなフローに従って実行することが可能である。

このようにプロセスを組み合わせることで、配列に存在し得る多くの種類の多型を網羅的に高い検出力で検出することが可能である。このようなプロセスの組み合わせは、例えば、複数のプロセスを同時に、並行して、または逐次的に行うことによって達成することが可能である。例えば、図１５Ｂに示されるようなクラスター構造を有する計算システムによって、異なるノードを用いてそれぞれのプロセスを行うことにより、組み合わせの方法を実装することができる。

置換、挿入、欠失、逆位または転座を検出するのに有用なプロセスとしては、対象配列データの配列中の少なくとも２ヶ所の部分配列の、コントロール配列上の位置を特定する工程を含むプロセスがある。ここで、部分配列は、ｋ長の部分配列を用いることができる。好ましくは、コントロール配列は、配列上の位置情報が特定できる配列であり、より好ましくは、コントロール配列はリファレンス配列である。

プロセスは、対象配列データにおける部分配列間の位置関係と、コントロール配列上の部分配列間の位置関係とを比較する工程を含み得る。ここで、対象配列データにおける部分配列間の位置関係と、コントロール配列上の部分配列間の位置関係とが異なっている場合、目的とする多型があると判定することができる。例えば、部分配列が、コントロール配列の異なる配列構造体上に存在する場合、転座が生じていると判定すること、部分配列が、コントロール配列の同一の配列構造体上に存在し、かつ、向きが対象配列データ上のものと異なっている場合、逆位が存在すると判定すること、部分配列が、コントロール配列の同一の配列構造体上に存在し、向きが対象配列データ上のものと同一であり、部分配列の距離が、コントロール配列上で対象配列データ上の距離より短い場合、欠失が存在すると判定すること、および／または部分配列が、コントロール配列の同一の配列構造体上に存在し、向きが対象配列データ上のものと同一であり、部分配列の距離が、コントロール配列上で対象配列データ上の距離より長い場合、挿入が存在すると判定することを含み得る。位置関係が異ならない場合、処理を終了してもよく、目的とする多型はないと判定してもよく、対象配列データにおける部分配列部位間の文字を、対応するコントロール配列上の文字と比較して不一致となる部位を検出する工程をさらに行い、不一致となる部位が存在する場合、置換が存在すると判定してもよい。

プロセスは、位置関係が異なっている場合、目的とする多型があると判定し、対象配列データにおける部分配列部位間の文字を、対応するコントロール配列上の文字と、部分配列部位を始点として順次比較して不一致となる部位を検出する工程を含み得る。かかる工程により、検出した多型の境界塩基を検出することができる。

例えば、このようなプロセスを組み合わせた場合、本発明の１つの実施形態では、
対象配列データにおいてコントロール配列データに対する多型を検出する方法であって、
（１）ａ）該対象配列データの長さｋの部分配列のサブセットの各々の部分配列の出現頻度を提供する工程と、
ｂ）該コントロール配列データの長さｋの部分配列のサブセットの各々の部分配列の出現頻度を提供する工程と、
ｃ）対象配列とコントロール配列とを比較し、該出現頻度の分布の比較に基づいて、多型を検出する工程と
によって、対象配列データにおける置換、コピー数多型、ＳＴＲＰ、挿入、欠失、逆位または転座を検出するプロセスと、
（２）ａ）該対象配列データの配列中の少なくとも２ヶ所の部分配列の、該コントロール配列上の位置を特定する工程と、
ｂ）対象配列データにおける該部分配列間の位置関係と、コントロール配列上の該部分配列間の位置関係とを比較する工程と、
ｃ）対象配列データにおける該部分配列間の位置関係と、コントロール配列上の該部分配列間の位置関係とが異なっている場合、目的とする多型があると判定し、該対象配列データにおける該部分配列部位間の文字を、対応するコントロール配列上の文字と、該部分配列部位を始点として順次比較して不一致となる部位を検出する工程と
によって、対象配列データにおける挿入、欠失、逆位、転座または置換を検出するプロセスと
を包含する、方法が提供される。

本発明のさらなる実施形態では、例えば、
対象配列データにおいてリファレンス配列データに対する多型を検出する方法であって、リファレンス配列データから、各々の位置情報と関連付けられたリファレンス配列のｋ長の部分配列セットを作成する工程を含み、さらに、
（Ａ１）該対象配列データの長さｋの部分配列のサブセットを生成し、ユニークな長さｋの部分配列の出現頻度を提供する工程と、
（Ａ２）該リファレンス配列のｋ長の部分配列セットの、ユニークな長さｋの部分配列の出現頻度を提供する工程と、
（Ａ３）該対象配列と該リファレンス配列とを比較し、該出現頻度の分布の比較に基づいて、挿入、欠失、置換、コピー数多型、ＳＴＲＰ、逆位または転座を検出する工程とを包含するプロセスと
（Ｂ１）該対象配列データの配列中の少なくとも２ヶ所のｋ長の部分配列をクエリとして、該リファレンス配列のｋ長の部分配列セットに対して二分検索を行い、該少なくとも２ヶ所の部分配列の、リファレンス配列上の位置を特定する工程と、
（Ｂ２）該対象配列データにおける該部分配列間の位置関係と、該リファレンス配列上の該部分配列間の位置関係とを比較する工程と、
（Ｂ３）該対象配列データにおける該部分配列間の位置関係と、該リファレンス配列上の該部分配列間の位置関係とが異なっている場合、挿入、欠失、逆位または転座があると判定し、該対象配列データにおける該部分配列部位間の文字を、対応するコントロール配列上の文字と、該部分配列部位を始点として順次比較して不一致となる部位を検出する工程を包含し、必要に応じて、
（Ｂ４）該位置関係が異ならない場合に、該対象配列データにおける該部分配列部位間の文字を、対応する前記コントロール配列上の文字と比較して不一致となる部位を検出する工程をさらに含み、不一致となる部位が存在する場合、置換が存在すると判定する工程をさらに含む、プロセスと、
を、同時に、並行して、または逐次的に行うことを特徴とする、方法が提供される。

（一般技術）
本明細書において用いられる分子生物学的手法、生化学的手法、微生物学的手法、バイオインフォマティクスは、当該分野において公知であり、周知でありまたは慣用される任意のものが使用され得る。

本明細書において「または」は、文章中に列挙されている事項の「少なくとも１つ以上」を採用できるときに使用される。「もしくは」も同様である。本明細書において「２つの値」の「範囲内」と明記した場合、その範囲には２つの値自体も含む。

本明細書において引用された、科学文献、特許、特許出願などの参考文献は、その全体が、各々具体的に記載されたのと同じ程度に本明細書において参考として援用される。

以上、本発明の理解を容易にするために好ましい実施形態を示して説明してきた。以下に、実施例に基づいて本発明を説明するが、上述の説明および以下の実施例は、例示の目的のみに提供され、本発明を限定する目的で提供したのではない。従って、本発明の範囲は、本明細書に具体的に記載された実施形態にも実施例にも限定されず、特許請求の範囲によってのみ限定される。

（実施例１：イネＳＮＰ検出および検証）
（概要）
本発明の方法を用いて、以下のイネ配列データを用いてイネのＳＮＰの検出および検証を行った。
対象：N1S5、N1S6、N1S7、N1S10
コントロール：N1
リファレンス：イネリファレンスゲノム（IRGSP1.0）
なお、本明細書において、サンプルの表記は、大文字で表記しても小文字で表記しても同じものを指すことに留意されたい。

（材料および方法）
（サンプル）
イネ品種日本晴の種子１粒をN1個体とし、発芽させて育てて葉をサンプリングした。N1個体に実った種子をN1S1とし、その種子１粒を発芽させて育てて葉をサンプリングした。N1S1個体に実った種子をN1S2とし、同様にして、N1S10世代まで、葉と種子をサンプリングした。

（次世代シーケンシング）
N1、N1S5、N1S6、N1S7、N1S10の葉のサンプルについて、次世代シーケンサーでの全ゲノム解析を行った。シーケンサーは、Ｉｌｌｕｍｉｎａ社のＨｉＳｅｑ２０００を用い、ペアエンドで解析した。リード長はN1、N1S5、N1S6、N1S7は１００塩基、N1S10のみ１０１塩基であった。

シーケンスライブラリーを用いて、シーケンスの鋳型となるクラスターを形成し、鋳型ＤＮＡの塩基配列を取得した。シーケンスデータの解析は付属のソフトウェアを使用しベースコールを行い、ｆａｓｔｑ形式ファイルとして出力した。

N1、N1S6の解析は以下のとおり製造業者のマニュアルに従った条件で行った。
表５：各作業に使用したマニュアル名、バージョン番号

表６：クラスター形成、シーケンス及びシーケンス解析に使用した装置、試薬、ソフトウェア

N1S5、N1S7、N1S10の解析は以下のとおり製造業者のマニュアルに従った条件で行った。
表７：各作業に使用したマニュアル名、バージョン番号

表８：クラスター形成、シーケンス及びシーケンス解析に使用した装置、試薬、ソフトウェア

（多型の検出）
各サンプルについて次世代シーケンサーから得られたＦａｓｔｑ形式のファイルより、配列内にＮを含まない塩基配列を選んで、相補鎖の配列と合わせて、それぞれ出力した。

出力されたファイル（reads）を辞書順にソートし同一の配列は一つにまとめた。ｕｎｉｘコマンドは以下の通りであった。

ソートされた各塩基配列をそれぞれ５’末端から１塩基ずつずらしながら、ｋ−ｍｅｒ（本実施例では２０塩基）の配列を対象の塩基配列の３’末端に到達するまで繰り返し出力（k-mer_file）した。

出力されたｋ−ｍｅｒの配列を辞書順にソートし、同一配列は一つにまとめて、出現回数を配列と共に表記したファイルを作成した。ｕｎｉｘコマンドは以下の通りであった。

ここで、コントロール配列：N1、対象配列：N1S7について、ｋ−ｍｅｒ配列を整列させ、多型の検出を行った。多型が存在している部位が、ゲノム上でユニークな配列部位であり、そこに例えば一塩基置換が存在している場合、その置換を含むｋ−ｍｅｒ配列は対象配列では存在するが、コントロール配列では存在しないと考えられる。したがって、当該置換変異を含むｋ−ｍｅｒ配列では、対象配列では存在し、コントロール配列で存在しないため、結果的に当該置換変異部位を含むｋ個のｋ−ｍｅｒ配列で出現頻度の差異が見られると考えられる。ｋ−ｍｅｒ配列の整列は、リファレンス配列に沿う形で整列させた。

k-mer_count_fileの各行のデータに対して、配列の５’末端よりｋ−１塩基の配列を得て、３’末端の塩基すなわちｋ番目の塩基をＡ、Ｃ、Ｇ、Ｔの出現回数として表記したデータに変換した。「ｋ−１ｍｅｒの配列 Ａの回数 Ｃの回数 Ｇの回数 Ｔの回数」という形式で出力した。

コントロール配列と対象配列について、上記工程で出力したデータをｋ−１ｍｅｒの配列でまとめたデータを作成した。ｕｎｉｘコマンドは以下の通りであった。

さらに、joinコマンドで得られたデータでcontrolとtargetで異なる塩基でカウントが示されるデータを検出した。本実施例においては、イネゲノムの４０〜５０倍読みの次世代シーケンシングデータを対象にしたため、塩基の頻度が１００を超える場合はリピート配列部位として除外した。次にcontrol、targetの両方でカウントが１塩基以下の塩基が２個以上存在し、controlあるいはtargetで１０以上のカウントを示した塩基に対応するtargetあるいはcontrolの塩基のカウントが１以下である事例が１ないし２回ある場合、多型の境界を検出したものとし、当該事例が生じているｋ−１ｍｅｒの配列を、多型部位を含むものとして検出した。

リファレンスゲノムへのマッピングのために、リファレンスゲノム（イネリファレンスゲノム（IRGSP1.0））を上記工程と同様に各染色体を５’側から１塩基ずつずらしながらｋ−ｍｅｒで切り出し、ｋ−ｍｅｒ配列、染色体番号、位置、向きを同一行に表記したデータを作成し、ｋ−ｍｅｒの配列で辞書順に整列させた。当該データに対して、二分検索法により多型ｋ−１ｍｅｒの配列のゲノム上の位置を決定した。

リファレンスゲノムにマップされた多型の確認を行った。上記工程で作成されたsort_uniq配列（ターゲット）と同様に、リファレンス配列も５’末端から２塩基ごとにリード長と同じ長さの配列を切り出し、相補鎖とともにsort、uniq処理をしたデータを作成した。このリファレンスとターゲットのsort_uniqデータを検索の対象とした。

次に、リファレンスゲノム配列から、多型位置を含むリード長と同じ長さの配列セット（検索を行うクエリセット）を切り出した。この配列セットは配列の５’末端から３’末端までのすべての位置で多型位置の塩基を含む配列セットである。配列セット数は配列長と同一となる。リファレンスの配列セット（リファレンスセット）に対して、多型位置の塩基を予想された多型塩基に置換した配列セット（ターゲットセット）を作成した。

リファレンスセット、ターゲットセットをクエリにして、リファレンスゲノム、及び、次世代シーケンサーから得られたＦａｓｔｑ形式のファイルより出力されたファイル（reads）を辞書順にソートし同一の配列を一つにまとめたsort_uniq配列をそれぞれ、二分検索法で検索し、それぞれのセットに対してsort_uniq配列が何個マッチするかを調べた。リファレンスセットではリファレンスのsort_uniq配列のみにヒットする。これに対して、ターゲットセットでは、ホモ型の変異の場合は、ターゲットのsort_uniqデータのみから検出され、ヘテロ型の変異の場合は、リファレンスと、ターゲットのsort_uniqデータの両方から検出されるはずである。このようにして、予想と一致する検索値を示した多型を抽出した。

（ジェノタイプの確認）
ジェノタイプの確認は、ＰＣＲで当該領域を増幅し、サンガー法で決定した。

各領域の増幅に用いたプライマー配列は以下のとおりであった。

ＰＣＲ反応の反応条件は以下のとおりであった。

反応サイクルは、９４℃ ０．５分、６０℃ ０．５分、７２℃ １分を３０サイクルとして行った。

増幅したＤＮＡ断片を１％アガロースゲル電気泳動で分離し、０．５μｇ／ｍｌの濃度のエチジウムブロマイドで染色して、長波長紫外線ランプ（３６５ｎｍ）で蛍光を発するバンドを切り出し、Ｐｒｏｍｅｇａ社のＷｉｚａｒｄ（登録商標） ＳＶ Ｇｅｌ ａｎｄ ＰＣＲ Ｃｌｅａｎ−Ｕｐ Ｓｙｓｔｅｍ（Ｃａｔ．＃ Ａ９２８２）で断片の精製を行った。

サンガー法でのＳＮＰを含む塩基配列の確認
精製された断片をＢｉｇＤｙｅ（登録商標） Ｔｅｒｍｉｎａｔｏｒ ｖ３．１ Ｃｙｃｌｅ Ｓｅｑｕｅｎｃｉｎｇ Ｋｉｔ（Ｔｈｅｒｍｏ Ｆｉｓｈｅｒ Ｓｃｉｅｎｔｉｆｉｃ Ｃａｔ．＃ ４３３７４５５）で反応を行い、ＤＮＡシーケンサー ＡＢＩ ＰＲＩＳＭ ３１３０ｘｌで塩基配列の確認を行った。

（結果）
（シーケンシング）
イネ個体（N1、N1S5、N1S6、N1S7、N1S10）を次世代シーケンサーで解析した結果の塩基配列データは、ＤＤＢＪに送信されており、以下のアクセッション番号で登録されている。

各サンプルについてのリード数（総データ数）は、以下：

のとおりであった。

出力されたＦａｓｔｑファイルを処理したsort_uniqのデータ数は以下：

のとおりであった。sort_uniqは、Ｎを含まないリードとその相補鎖のデータをsortしたのちuniqで同一配列を一つにまとめたものである。このデータは、ｋ−ｍｅｒのデータと異なり配列のみのデータであり、頻度の数値データは含まない。

（多型の検出）
コントロール配列：N1、対象配列：N1S7について、ｋ−ｍｅｒ配列をリファレンス配列に沿って整列させ、多型の検出を行った結果は、図７および８に示される。下線を付された塩基がコントロールと対象で異なっており、多型が検出されたことが示されている。図９においては、コントロール配列（N1）と対象配列（N1S5、N1S6、N1S7、N1S10）のｋ−ｍｅｒ配列の部分配列サブセットを整列させることによって、イネリファレンス配列の対応する位置から始まる配列と対応する配列の出現頻度を求めた結果が示される。染色体番号、染色体の位置に続いて、N1、N1S5、N1S6、N1S7、N1S10の２０−ｍｅｒの頻度が示される。N1S7でヘテロ、N1S10でミュータントホモになり、野生型の２０−ｍｅｒがゼロになっていることがわかる。すなわち、ｋ−ｍｅｒ配列の出現頻度の変化から、N1S7においてヘテロ変異が生じ、N1S10においてホモ変異となったことを検出することができた。

さらに、最終的に上記手順によって、コントロール配列（N1）と対象配列（N1S5、N1S6、N1S7、N1S10）との間で検出された多型の一部を、図１０および１１に示した。これらの結果は、リファレンス配列を用いた確認と一致していた。Ｗが野生型、Ｈがヘテロ型、Ｍがミュータントホモ型を示す。

サンガー法により確認した各サンプルにおける多型を、図１２に示す。N1、N1S1、N1S2、N1S3、N1S4、N1S5、N1S6、N1S7、N1S8、N1S9、N1S10と世代を重ねながら、ヘテロ変異が生じ、その後ホモ変異として定着する様子が観察でき、この結果は、本発明の方法によって検出された多型とよく一致していた（図１２）。

したがって、本発明の方法によって、世代間に生じた多型を詳細に検出できることが示された。また、Ｐｏｌｙｍｏｒｐｈｉｃ Ｅｄｇｅ Ｄｅｔｅｃｔｉｏｎによって検出された多型が、サンガー法によっても確認されていることから、リファレンス配列（ゲノムリファレンス配列）を必要とせずに、配列データ間での多型の検出を行うことができることが実証された。

（実施例２：ナイジェリアのヨルバ族男性（NA18507）の配列解析）
（材料および方法）
コントロール配列データとして、ヒトゲノムリファレンスhg38を用いた。配列は、ftp://hgdownload.cse.ucsc.edu/goldenPath/hg38/chromosomes/よりchr1〜chr22とchrX、chrY、chrMの染色体データをダウンロードして用いた。alt、v1等のファイル名にコメントが付いているデータは除外した。

対象配列データとしては、ヒトゲノムの次世代シーケンス配列データNA18507をダウンロードして用いた。この配列データは、Ｉｌｌｕｍｉｎａ社の次世代シーケンサーで解析が行われたものであり、ＮＣＢＩに登録・公開されており、当該データをダウンロードして用いた。塩基配列セットの実験ＩＤのＵＲＬは、https://www.ncbi.nlm.nih.gov/sra/SRX016231であり、配列のアクセッション番号は、SRR034939〜SRR034975の範囲であった。

情報処理は、実施例１の（多型の検出）と同様に行った。

（結果）
結果の一部を、図１３に示す。相補鎖（ｒ）で検出された野生型・変異型の塩基は順鎖に変換して表示している。ヘテロ型の場合はそれぞれの塩基を並べて示している。参照、対象の配列ｘのＡ、Ｃ、Ｇ、Ｔの数はｋ−１配列に続く配列ｘの各塩基の頻度を示している。Ｐ〜Ｑ列に関しては、野生型あるいは変異型の塩基をもつ参照あるいは対象配列の数を示している。

順鎖、逆鎖の双方で同じ変異が検出された場合は、ほぼ間違いなくＳＮＰである。片側の鎖のみで検出される原因は、逆の鎖の配列がゲノム上でユニークでなく、一意で検出できなかったため、あるいは、挿入、欠失、転座、逆位の境界塩基を検出したためであると考えられる。

本発明の方法は、ヒトゲノムについても、正確に多型を検出できることが実証された。また、コントロール配列として、データベースから取得したリファレンス配列を用いることができることも示される。

（実施例３：同一個体の組織間での多型の検出）
（概要）
本発明の方法により、同一個体の組織間での多型の検出が可能であることを実証する。

（材料および方法）
ＮＣＢＩのＳＲＡよりfastq-dumpを用いて配列データを取得した。本データは、Texas Cancer Research Biobank Open Access Data Sharing: Genome Projectが登録したデータであり、詳細データについて、以下のＵＲＬ：https://www.ncbi.nlm.nih.gov/Traces/study/?acc=SRP060654において提供されている（Becnel, L. et al. NCBI Sequence Read Archive PRJNA284598 (2015).）。本発明の方法により、前記配列データにおいて多型を検出し、同一個体の正常組織と腫瘍組織との間の多型を本発明の方法によって検出できるかを検証した。

配列データの起源のサンプルは２種類で、配列データ名とサンプルの内容は以下の通りであった。
SRR2096532 コントロール血液 （Ｎｏｒｍａｌ）
SRR2096535 濾胞性リンパ腫 （９６９０／３： Ｆｏｌｌｉｃｕｌａｒ ｌｙｍｐｈｏｍａ）
リード数 （配列長１０１塩基）
SRR2096532 1300353764
SRR2096535 1339310760
sort_uniqの配列数
SRR2096532 2056683322
SRR2096535 2181081390

SRR2096532（正常組織）をコントロール配列データとして用い、SRR2096535（腫瘍組織）を対象配列データとして用いて解析を行った。

情報処理は、実施例２の（多型の検出）と同様に行った。

なお、本実施例においては、確認（ｖｅｒｉｆｙ）工程で、検出されたそれぞれの変異に対して、ターゲットでは、変異型が５リード以上、野生型が１リード以下、コントロールでは変異型が１リード以下、野生型が５リード以上の場合、ホモ型変異（Ｍ）とマークし、ターゲットのリード総数に対する変異型の割合が０．３より大きいか、０．７より小さく、かつ、コントロールで変異型が１リード以下、かつ、ターゲットで野生型リードが５以上の場合、ヘテロ型変異（Ｈ）とした。

［従来法］
この分野で一般的に広く用いられている、Ｓａｍｔｏｏｌｓを用いて同じデータを処理した。
従来法による解析は、以下の工程によって行った。工程０は準備で１回のみ実施し、工程１から５はサンプル毎に実施した。
０．準備：リファレンス配列にインデックス付加
１．ショートリードデータのマッピング
２．ＳＡＭ形式をＢＡＭに変換（マッピング位置でソートも）
３．Ｓａｍｔｏｏｌｓで多型部位の検出

リファレンス配列データとして、ヒトゲノムリファレンスhg38を用いた。配列は、ftp://hgdownload.cse.ucsc.edu/goldenPath/hg38/chromosomes/よりchr1〜chr22とchrX、chrY、chrMの染色体データをダウンロードして用いた。alt、v1等のファイル名にコメントが付いているデータは除外した。

（結果）
本発明の方法（ＰＥＤ）および従来法（ｂｗａ＋Ｓａｍｔｏｏｌｓ）で検出された多型の数は以下のとおりであった。

ＰＥＤでＶｅｒｉｆｙ（確認工程）にかけた座位数は２２６０１で、そのうち順鎖、相補鎖共に置換変異が検出されたのは、５１４座位であった。Ｖｅｒｉｆｙの結果片側鎖のみで検出されたものの、ポジティブと判断された座位を合わせて１０４２座位に置換変異があると予想された。ＰＥＤで検出された１０４２座位はすべてヘテロ型であった。このことは、同一人物の血液と癌組織からの配列の比較であるため、原理的にホモ型変異が起こる可能性は非常に低いと考えられることと整合する。ＰＥＤではマッピングする前に直接対象とコントロールの比較を行っているため、対象・コントロールのＳＮＰ数は出力されない。

ＰＥＤで検出されたが、ｂｗａ＋Ｓａｍｔｏｏｌｓで検出できなかったＳＮＰは２０であった。ｂｗａ＋Ｓａｍｔｏｏｌｓではホモ型変異も多く検出されており、結果に非常に多くのノイズを含んでいると考えられる。

（考察）
同一人物からの組織間の比較なので、ミュータントホモはあり得ないと予想されていたところ、予想通りにヘテロ型のみが検出された。この精度で多型を検出できる系は、他にはなく、本発明の方法は従来技術に対して顕著に有利なものであると考えられる。

上記の結果から、同一個体の正常組織と癌組織の配列を直接比較して多型を検出することも可能であることが示された。本発明の方法は、いったんリファレンスゲノムにマップしてから差を調べる方法に比べて、多型の捕捉率、捕捉精度ともに非常に高くなる。

ＰＥＤではマッピングの前に、リード同士の比較からＳＮＰを検出するので確度の高いＳＮＰが得られる。本実施例で得られた結果を、挿入・欠失を検出するプログラムの結果と併せて考慮すると、癌細胞では、置換変異よりむしろ、二重鎖切断の後の除去修復に起因すると考えられる挿入・欠失変異の頻度が高いと考えられる。

（実施例４：コピー数多型の検出）
（材料および方法）
実施例１におけるN1S6と同じ世代の種子１粒から細胞培養を行い、１ヶ月、３ヶ月、５ヶ月後に再分化してイネの個体にした葉よりＤＮＡを抽出し、それぞれ1M1、3M1、5M1のサンプルとして用いた。実施例１におけるN1種子と同じ世代の種子を５ヶ月培養して再分化した個体を４世代自殖した個体からＤＮＡを抽出し、TTM2とTTM5のサンプルとして用いた。抽出したＤＮＡから次世代シーケンサーによって配列データを取得した。シーケンシングのプロトコルは実施例１におけるものと同様であった。これらの配列データおよび実施例１のN1S5、N1S6、N1S7、N1S10の配列データを対象配列データとして、N1を参照配列データとして用いた。これらの配列データについてのアクセッション番号、リード数、sort_uniqの配列数は以下の表に示される。なお、TTM5のデータは、SRR556174とSRR556175の２つのアクセッション番号に分割されている。TTM5のsort_uniqは２つのリードを合わせて一つのファイルとして作成した。

情報処理は、実施例１の（多型の検出）と同様に行い、参照配列データと、対象配列データの間で出現頻度が大きく異なったｋ−ｍｅｒを検出した。コントロールのN1より２倍程度以上の頻度が連続して検出された場合にＣＮＶと判定した。誤検出も出てくるので、２．５倍以上と厳しくすることもできる。

（結果）
結果の一部を図１７に示した。第７染色体の２６６９４７９５位置（図７中矢印で示される）に対応するｋ−ｍｅｒ配列から、それまでの位置に対応するｋ−ｍｅｒ配列の出現頻度の４倍程度の値が、N1〜N1S10で現れていることが分かる。

この位置からレトロトランスポゾンTos17の配列が始まる。このトランスポゾンはゲノムに２コピー存在しており、それぞれのトランスポゾンの末端部分は同じ配列（Ｌｏｎｇ Ｔｅｒｍｉｎａｌ Ｒｅｐｅａｔ、ＬＴＲ）を有する。そのため、第７染色体の２６６９４７９５より、それまでの４倍程度の値が、N1〜N1S10で現れたと考えられる。このトランスポゾンの全長は４．１ｋｂあるため、図１７には最初のジャンクションの部分のみ示している。

1M1、3M1、5M1はN1S6と同じ世代の種子１粒から細胞培養を行い、それぞれ１ヶ月、３ヶ月、５ヶ月後に再分化してイネの個体にした葉よりＤＮＡを抽出して解析したものであり、図１７に示される結果から、培養時間に応じてカウントが増えていることがわかる。これは、培養時間に比例してトランスポゾンが転移してコピー数が増加し、コピー数多型（ＣＮＶ）が生じたためと考えられる。TTM2とTTM5はN1種子と同じ世代の種子を５ヶ月培養して再分化した個体を４世代自殖した個体のＤＮＡであり、こちらも同様にコピー数の増加が認められる。

（考察）
Tos17は培養時のみ活性化されて転移するイネのトランスポゾンとして知られている。Tos17はレトロトランスポゾンなので、オリジナルは切り出されて転移することなく、Tos17のコピーの転移によってゲノム上のTos17のコピー数が増加する。そのため、Tos17は培養で転移してコピー数が増えることが以前から知られている。

本発明の方法によって、コピー数の変異を検出した結果、理論とよく一致して、培養時間が長くなるにつれて、Tos17のコピー数が増加していることが図１７に示されるカウントから理解される。

したがって、本発明の方法によって、コピー数多型を検出することが可能であることが実証される。

例えば、このようなコピー数多型の検出は、培養細胞（例えば、ｉＰＳ細胞等）において、品質の管理に用いることが可能であると考えられる。上記イネの培養細胞と同様にトランスポゾン等によるコピー数変異が、例えば、ヒトの培養細胞（ｉＰＳ等）で観察されている場合には、例えば、治療に用いるのは危険である可能性が高いという判断に用いることが可能である。

（関連出願）
本出願は、２０１７年７月２４日に出願された特願２０１７−１４２７８１号の優先権の利益を主張し、当該出願は、全ての目的において、その開示全体が本明細書において参考として援用される。さらに、本明細書において、本出願人により本出願と同日に出願された「挿入・欠失・逆位・転座・置換検出法」との名称の出願（整理番号ＮＧ０１２ＰＣＴ／Ｆ５−１８ＰＣＴ０７５）およびその基礎出願である２０１７年７月２４日に出願された特願２０１７−１４２７８２号（整理番号Ｊ１−１７３６９１６２）は、全ての目的において、その開示全体が本明細書において参考として援用される。

塩基配列解析で多型を検出するすべての分野で利用が可能で、ＤＮＡ育種利用の他、臨床検査、ｉＰＳ細胞の検査、メタゲノム解析、発現解析等、幅広い分野で利用することができる。

配列番号１〜６０：図７のｋ−ｍｅｒ参照配列
配列番号６１〜８０：図７のｋ−ｍｅｒ対象配列（変異が存在する部分）
配列番号８１〜１４０：図８のｋ−ｍｅｒ参照配列
配列番号１４１〜１６０：図８のｋ−ｍｅｒ対象配列（変異が存在する部分）
配列番号１６１〜１９０：図１３のｋ−１（ｋ＝２０）配列
配列番号１９１〜２２１：（具体的な例）で用いられた配列
配列番号２２２〜２３２：（例示的計算フロー）で用いられた配列
配列番号２３３〜２６６：実施例１で用いられたプライマーの配列
配列番号２６７〜２７５：図１１の配列

Claims (42)

1. 対象配列データにおいてコントロール配列データに対する多型を検出する方法であって、
   ａ）該対象配列データの長さｋの部分配列のサブセットの各々の部分配列の出現頻度を提供する工程と、
   ｂ）該コントロール配列データの長さｋの部分配列のサブセットの各々の部分配列の出現頻度を提供する工程と、
   ｃ）対象配列とコントロール配列とを比較し、該出現頻度の分布の比較に基づいて、多型を検出する工程と
   を包含し、ここで、ｋは該対象配列および該コントロール配列のいずれか短いほうの全長以下の整数である、方法。
2. 前記部分配列中の長さｋ−ｘの配列部分が共通する配列ごとに、長さｘの部分について出現頻度の分布を算出する工程をさらに含み、ここで、ｘはｋ未満の正の整数である、請求項１に記載の方法。
3. 前記比較が、前記部分配列中の長さｋ−ｘの配列部分が共通する配列における、長さｘの部分の出現頻度の分布の差異の比較を含む、請求項２に記載の方法。
4. 前記部分配列中の長さｋ−ｘの配列部分を、ユニークな配列ごとにグルーピングする工程を含み、ここで、ｘはｋ未満の正の整数である、請求項１〜３のいずれか１項に記載の方法。
5. 前記長さｋ−ｘの配列部分をソートする工程を含む、請求項４に記載の方法。
6. 前記長さｋ−ｘの配列部分を文字列としてソートする工程を含む、請求項５に記載の方法。
7. 前記ｋが、前記対象配列における偶然同一を排除する長さである、請求項１〜６のいずれか１項に記載の方法。
8. 前記対象配列データおよび前記コントロール配列データが、生物のゲノムに由来する塩基配列データであり、前記ｋが、前記生物のゲノムにおいて、異なる箇所での偶然同一を排除する長さである、請求項１〜７のいずれか１項に記載の方法。
9. 長さｘが１〜２である、請求項２〜８のいずれか１項に記載の方法。
10. 長さｘが１である、請求項９に記載の方法。
11. 前記長さｘの部分が、前記部分配列の末端に存在する、請求項２〜１０のいずれか１項に記載の方法。
12. 前記対象配列データおよび前記コントロール配列データが塩基配列データであり、前記長さｘの部分が、前記部分配列の３’末端である、請求項９に記載の方法。
13. 前記コントロール配列データのサブセットと前記対象配列データのサブセットとの間で、前記長さｘの部分の配列の出現頻度が異なる場合、該長さｘの部分の配列を、対象配列データにおけるコントロール配列データに対する多型として検出する、請求項２〜１２のいずれか１項に記載の方法。
14. 前記コントロール配列データのサブセットと前記対象配列データのサブセットとの間で、前記長さｘの部分の配列で最も高頻度のものが異なっている長さｋ−ｘの配列部分が存在する場合、該長さｘの部分の配列を、対象配列データにおける多型として検出する、請求項２〜１３のいずれか１項に記載の方法。
15. 前記対象配列データおよび前記コントロール配列データが塩基配列データであり、前記対象配列データのサブセットにおける前記長さｘの部分の配列で、前記コントロール配列データのサブセットにおける最も高頻度のものと同一の長さｘの部分の配列がノイズ以下のカウントしか存在しない長さｋ−ｘの配列部分が存在する場合、該長さｘの部分の配列を、対象配列データにおけるホモ多型として検出する、請求項２〜１４のいずれか１項に記載の方法。
16. 前記対象配列データおよび前記コントロール配列データが塩基配列データであり、対象サブセットにおける前記長さｘの部分の配列で、コントロール配列データのサブセットにおける最も高頻度のものと同一の長さｘの部分の配列が存在し、かつ、コントロール配列データのサブセットにおける最も高頻度のものと異なる長さｘの部分の配列が存在する長さｋ−ｘの配列部分が存在する場合、該長さｘの部分の配列を、対象配列データにおけるヘテロ多型として検出する、請求項２〜１５のいずれか１項に記載の方法。
17. 対象配列データ量から予測される出現頻度と比較して、前記出現頻度が少ない部分配列をノイズとする、請求項１〜１６のいずれか１項に記載の方法。
18. 前記対象配列データおよび前記コントロール配列データが塩基配列データであり、［（対象配列データ量）×（１−精度）］／（対象ゲノムサイズ）＋１未満の出現頻度の部分配列をノイズとする、請求項１７に記載の方法。
19. 前記対象配列データが、次世代シーケンシングによって得られた塩基配列データである、請求項１〜１８のいずれか１項に記載の方法。
20. 前記対象配列データが、個体から得られた配列データであり、前記コントロール配列データが、該個体と同種の別の個体、またはデータベースから得られた配列データである、請求項１〜１９のいずれか１項に記載の方法。
21. 前記対象配列データが、個体の組織試料から得られた配列データであり、前記コントロール配列データが、該個体の別の組織、またはデータベースから得られた配列データである、請求項１〜２０のいずれか１項に記載の方法。
22. 前記対象配列データが、細胞試料から得られた配列データであり、前記コントロール配列データが、別の細胞、またはデータベースから得られた配列データである、請求項１〜２１のいずれか１項に記載の方法。
23. 前記多型が、置換、挿入、欠失、コピー数多型（Ｃｏｐｙ Ｎｕｍｂｅｒ Ｖａｒｉａｔｉｏｎ， ＣＮＶ）、ＳＴＲＰ（ｓｈｏｒｔ ｔａｎｄｅｍ ｒｅｐｅａｔ ｐｏｌｙｍｏｒｐｈｉｓｍ）、逆位または転座である、請求項１〜２２のいずれか１項に記載の方法。
24. 前記多型が、置換である、請求項２３に記載の方法。
25. 前記対象配列に対するリファレンス配列における前記多型の位置を特定する工程をさらに含む、請求項１〜２４のいずれか１項に記載の方法。
26. 前記対象配列データおよび前記コントロール配列データが、生物のゲノムに由来する塩基配列データであり、前記多型のゲノム上の位置を特定する工程をさらに含む、請求項１〜２５のいずれか１項に記載の方法。
27. 検出された多型の部位について、リファレンス配列またはコントロール配列から作成したクエリ配列セットを用いて、対象配列データおよび／またはコントロール配列データとの比較を行い確認する工程をさらに含む、請求項２５または２６に記載の方法。
28. 前記クエリ配列セットが、リファレンス配列またはコントロール配列において前記多型に該当する部位の文字を異なる文字に置換した変異型クエリ配列セットを含む、請求項２７に記載の方法。
29. 前記対象配列データおよび前記コントロール配列データが塩基配列データであり、検出された多型の部位について、対象配列データおよび／またはコントロール配列データの相補鎖の配列データに対して、リファレンス配列またはコントロール配列から作成したクエリ配列セットとの比較を行い確認する工程をさらに含む、請求項２７または２８に記載の方法。
30. 前記対象配列データおよび前記コントロール配列データが塩基配列データであり、検出された多型の部位について、対象配列データおよび／またはコントロール配列データの変異型の塩基を有する配列データに対して、リファレンス配列またはコントロール配列から作成したクエリ配列セットとの比較を行い確認する工程をさらに含む、請求項２７〜２９のいずれか１項に記載の方法。
31. 前記対象配列データおよび前記コントロール配列データが、生物のゲノムに由来する塩基配列データであり、前記ゲノムの配列が不明である、請求項１〜３０のいずれか１項に記載の方法。
32. 実験結果またはデータベースから対象配列データまたはコントロール配列データを取得する工程をさらに含む、請求項１〜３１のいずれか１項に記載の方法。
33. 対象配列データにおけるコントロール配列データに対する多型を含む部分配列中の多型ではない部分の少なくとも一部を含む配列を、該多型の識別子として割り当てることをさらに含む、請求項１〜３２のいずれか１項に記載の方法。
34. 前記多型の識別子をリファレンス配列にマッピングし、リファレンス上の該多型の位置を特定することを含む、請求項３３に記載の方法。
35. 対象配列データにおいてコントロール配列データに対する多型を検出する方法をコンピュータに実行させるためのプログラムであって、該方法は、
    ａ）該対象配列データの長さｋの部分配列のサブセットをコンピュータに保存する工程であって、ここで、ｋは対象配列およびコントロール配列の全長以下の長さである、工程と、
    ｂ）該対象配列データの長さｋのサブセットの各々の部分配列の出現頻度を算出する工程と、
    ｃ）該コントロール配列データの長さｋの部分配列のサブセットにおける各々の部分配列の出現頻度をコンピュータに保存する工程と、
    ｄ）対象配列とコントロール配列とを比較し、該出現頻度の分布の差異に基づいて、多型を検出する工程と
    を包含する、プログラム。
36. 前記方法が、前記部分配列中の多型ではない部分の少なくとも一部を含む配列を、検出された前記多型の名称として表示する工程をさらに含む、請求項３５に記載のプログラム。
37. 対象配列データにおいてコントロール配列データに対する多型を検出する方法をコンピュータに実行させるためのプログラムを格納する記録媒体であって、該方法は、
    ａ）該対象配列データの長さｋの部分配列のサブセットをコンピュータに保存する工程であって、ここで、ｋは対象配列およびコントロール配列の全長以下の長さである、工程と、
    ｂ）該対象配列データの長さｋのサブセットの各々の部分配列の出現頻度を算出する工程と、
    ｃ）該コントロール配列データの長さｋの部分配列のサブセットにおける各々の部分配列の出現頻度をコンピュータに保存する工程と、
    ｄ）対象配列とコントロール配列とを比較し、該出現頻度の分布の差異に基づいて、多型を検出する工程と
    を包含する、記録媒体。
38. 前記方法が、前記部分配列中の多型ではない部分の少なくとも一部を含む配列を、検出された前記多型の名称として表示する工程をさらに含む、請求項３７に記載の記録媒体。
39. 対象配列データにおいてコントロール配列データに対する多型を検出するためのシステムであって、該システムは、
    該対象配列データおよび該コントロール配列データの長さｋの部分配列のサブセットの各々の部分配列の出現頻度を提供するように構成された配列データ処理部であって、ここで、ｋは対象配列およびコントロール配列の全長以下の長さである、配列データ処理部と、
    対象配列とコントロール配列とを比較し、該出現頻度の分布の差異に基づいて、多型を検出する工程とを行うように構成された、配列データ計算部と
    を備える、システム。
40. 前記システムが、前記部分配列中の多型ではない部分の少なくとも一部を含む配列を、検出された前記多型の名称として表示する表示手段をさらに含む、請求項３９に記載のシステム。
41. 対象配列データとコントロール配列データとの比較方法であって、
    対象配列データにおけるコントロール配列データに対する多型を含む部分配列中の多型ではない部分の少なくとも一部を含む配列を、該多型の識別子として割り当てることを含む、方法。
42. 前記多型の識別子をリファレンス配列にマッピングし、リファレンス上の該多型の位置を特定することを含む、請求項４１に記載の方法。

Images