JP7166638B2 - 多型検出法 - Google Patents

多型検出法 Download PDF

Info

Publication number
JP7166638B2
JP7166638B2 JP2019532603A JP2019532603A JP7166638B2 JP 7166638 B2 JP7166638 B2 JP 7166638B2 JP 2019532603 A JP2019532603 A JP 2019532603A JP 2019532603 A JP2019532603 A JP 2019532603A JP 7166638 B2 JP7166638 B2 JP 7166638B2
Authority
JP
Japan
Prior art keywords
sequence
sequence data
length
data
polymorphism
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
JP2019532603A
Other languages
English (en)
Other versions
JPWO2019022018A1 (ja
Inventor
安藝雄 宮尾
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
National Agriculture and Food Research Organization
Original Assignee
National Agriculture and Food Research Organization
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by National Agriculture and Food Research Organization filed Critical National Agriculture and Food Research Organization
Publication of JPWO2019022018A1 publication Critical patent/JPWO2019022018A1/ja
Application granted granted Critical
Publication of JP7166638B2 publication Critical patent/JP7166638B2/ja
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G16INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16BBIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B20/00ICT specially adapted for functional genomics or proteomics, e.g. genotype-phenotype associations
    • G16B20/20Allele or variant detection, e.g. single nucleotide polymorphism [SNP] detection
    • GPHYSICS
    • G16INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16BBIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B25/00ICT specially adapted for hybridisation; ICT specially adapted for gene or protein expression
    • GPHYSICS
    • G16INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16BBIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B30/00ICT specially adapted for sequence analysis involving nucleotides or amino acids
    • G16B30/10Sequence alignment; Homology search

Landscapes

  • Bioinformatics & Cheminformatics (AREA)
  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Proteomics, Peptides & Aminoacids (AREA)
  • Biotechnology (AREA)
  • Biophysics (AREA)
  • Genetics & Genomics (AREA)
  • Molecular Biology (AREA)
  • Chemical & Material Sciences (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Analytical Chemistry (AREA)
  • Evolutionary Biology (AREA)
  • General Health & Medical Sciences (AREA)
  • Medical Informatics (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)
  • Information Retrieval, Db Structures And Fs Structures Therefor (AREA)
  • Apparatus Associated With Microorganisms And Enzymes (AREA)

Description

本発明は、配列情報、とりわけ、ゲノム等の生体分子の配列情報の情報処理の分野に関する。
次世代シーケンサーの出現により、生物の全ゲノム配列情報が得られるようになった。次世代シーケンサーの配列情報から多型情報を得て、表現型との関連を調べることにより、その表現型の原因となる遺伝子の特定につながる。正確な多型情報の取得は、作物育種のみならず、ヒトの遺伝病の診断、生物種・品種等の特定等、幅広い分野で必要とされる基盤技術であり、これまでにない精度で多型情報が得られれば、そのインパクトは大きい。
次世代シーケンサーからの塩基配列データを用いた多型の検出は、まず最初に配列データをbwa、またはbowtieのようなマッピングプログラムを用いてリファレンス配列上の位置情報とミスマッチの情報を得て、次に、SamtoolsやGATK等の多型抽出プログラムでSNPやindel等の多型情報を抽出するのが一般的である。
これらの方法では、多型の可能性のある部分は可能な限り出力するため、多くのノイズを含みこれらの技術のみでは、正確な多型解析が困難である。マイクロアレイやDNAチップ等の別の技術を併用して用いられているというのが現状である。
本発明において、2つ以上の配列の間における多型を検出する方法が提供される。本発明の方法は、配列データにおける個々の配列の全長配列における位置を考慮することなく、複数の配列データ間での多型の検出を可能にする。本発明の方法は、配列データ中の個々の配列(例えば、次世代シーケンサーからのショートリード)を連結してより長い配列とすること(例えば、アセンブリ)を必要とせずに、多型を検出することができることを1つの特徴とする。
例えば、本発明は以下の項目を提供する。
(項目1) 対象配列データにおいてコントロール配列データに対する多型を検出する方法であって、
a)該対象配列データの長さkの部分配列のサブセットの各々の部分配列の出現頻度を提供する工程と、
b)該コントロール配列データの長さkの部分配列のサブセットの各々の部分配列の出現頻度を提供する工程と、
c)対象配列とコントロール配列とを比較し、該出現頻度の分布の比較に基づいて、多型を検出する工程と
を包含し、ここで、kは該対象配列および該コントロール配列のいずれか短いほうの全長以下の整数である、方法。
(項目2) 前記部分配列中の長さk-xの配列部分が共通する配列ごとに、長さxの部分について出現頻度の分布を算出する工程をさらに含み、ここで、xはk未満の正の整数である、前記項目に記載の方法。
(項目3) 前記比較が、前記部分配列中の長さk-xの配列部分が共通する配列における、長さxの部分の出現頻度の分布の差異の比較を含む、前記項目のいずれかに記載の方法。
(項目4) 前記部分配列中の長さk-xの配列部分を、ユニークな配列ごとにグルーピングする工程を含み、ここで、xはk未満の正の整数である、前記項目のいずれかに記載の方法。
(項目5) 前記長さk-xの配列部分をソートする工程を含む、前記項目のいずれかに記載の方法。
(項目6) 前記長さk-xの配列部分を文字列としてソートする工程を含む、前記項目のいずれかに記載の方法。
(項目7) 前記kが、前記対象配列における偶然同一を排除する長さである、前記項目のいずれかに記載の方法。
(項目8) 前記対象配列データおよび前記コントロール配列データが、生物のゲノムに由来する塩基配列データであり、前記kが、前記生物のゲノムにおいて、異なる箇所での偶然同一を排除する長さである、前記項目のいずれかに記載の方法。
(項目9) 長さxが1~2である、前記項目のいずれかに記載の方法。
(項目10) 長さxが1である、前記項目のいずれかに記載の方法。
(項目11) 前記長さxの部分が、前記部分配列の末端に存在する、前記項目のいずれかに記載の方法。
(項目12) 前記対象配列データおよび前記コントロール配列データが塩基配列データであり、前記長さxの部分が、前記部分配列の3’末端である、前記項目のいずれかに記載の方法。
(項目13) 前記コントロール配列データのサブセットと前記対象配列データのサブセットとの間で、前記長さxの部分の配列の出現頻度が異なる場合、該長さxの部分の配列を、対象配列データにおけるコントロール配列データに対する多型として検出する、前記項目のいずれかに記載の方法。
(項目14) 前記コントロール配列データのサブセットと前記対象配列データのサブセットとの間で、前記長さxの部分の配列で最も高頻度のものが異なっている長さk-xの配列部分が存在する場合、該長さxの部分の配列を、対象配列データにおける多型として検出する、前記項目のいずれかに記載の方法。
(項目15) 前記対象配列データおよび前記コントロール配列データが塩基配列データであり、前記対象配列データのサブセットにおける前記長さxの部分の配列で、前記コントロール配列データのサブセットにおける最も高頻度のものと同一の長さxの部分の配列がノイズ以下のカウントしか存在しない長さk-xの配列部分が存在する場合、該長さxの部分の配列を、対象配列データにおけるホモ多型として検出する、前記項目のいずれかに記載の方法。
(項目16) 前記対象配列データおよび前記コントロール配列データが塩基配列データであり、対象サブセットにおける前記長さxの部分の配列で、コントロール配列データのサブセットにおける最も高頻度のものと同一の長さxの部分の配列が存在し、かつ、コントロール配列データのサブセットにおける最も高頻度のものと異なる長さxの部分の配列が存在する長さk-xの配列部分が存在する場合、該長さxの部分の配列を、対象配列データにおけるヘテロ多型として検出する、前記項目のいずれかに記載の方法。
(項目17) 対象配列データ量から予測される出現頻度と比較して、前記出現頻度が少ない部分配列をノイズとする、前記項目のいずれかに記載の方法。
(項目18) 前記対象配列データおよび前記コントロール配列データが塩基配列データであり、[(対象配列データ量)×(1-精度)]/(対象ゲノムサイズ)+1未満の出現頻度の部分配列をノイズとする、前記項目のいずれかに記載の方法。
(項目19) 前記対象配列データが、次世代シーケンシングによって得られた塩基配列データである、前記項目のいずれかに記載の方法。
(項目20) 前記対象配列データが、個体から得られた配列データであり、前記コントロール配列データが、該個体と同種の別の個体、またはデータベースから得られた配列データである、前記項目のいずれかに記載の方法。
(項目21) 前記対象配列データが、個体の組織試料から得られた配列データであり、前記コントロール配列データが、該個体の別の組織、またはデータベースから得られた配列データである、前記項目のいずれかに記載の方法。
(項目22) 前記対象配列データが、細胞試料から得られた配列データであり、前記コントロール配列データが、別の細胞、またはデータベースから得られた配列データである、前記項目のいずれかに記載の方法。
(項目23) 前記多型が、置換、挿入、欠失、コピー数多型(Copy Number Variation, CNV)、STRP(short tandem repeat polymorphism)、逆位または転座である、前記項目のいずれかに記載の方法。
(項目24) 前記多型が、置換である、前記項目のいずれかに記載の方法。
(項目25) 前記対象配列に対するリファレンス配列における前記多型の位置を特定する工程をさらに含む、前記項目のいずれかに記載の方法。
(項目26) 前記対象配列データおよび前記コントロール配列データが、生物のゲノムに由来する塩基配列データであり、前記多型のゲノム上の位置を特定する工程をさらに含む、前記項目のいずれかに記載の方法。
(項目27) 検出された多型の部位について、リファレンス配列またはコントロール配列から作成したクエリ配列セットを用いて、対象配列データおよび/またはコントロール配列データとの比較を行い確認する工程をさらに含む、前記項目のいずれかに記載の方法。
(項目28) 前記クエリ配列セットが、リファレンス配列またはコントロール配列において前記多型に該当する部位の文字を異なる文字に置換した変異型クエリ配列セットを含む、前記項目のいずれかに記載の方法。
(項目29) 前記対象配列データおよび前記コントロール配列データが塩基配列データであり、検出された多型の部位について、対象配列データおよび/またはコントロール配列データの相補鎖の配列データに対して、リファレンス配列またはコントロール配列から作成したクエリ配列セットとの比較を行い確認する工程をさらに含む、前記項目のいずれかに記載の方法。
(項目30) 前記対象配列データおよび前記コントロール配列データが塩基配列データであり、検出された多型の部位について、対象配列データおよび/またはコントロール配列データの変異型の塩基を有する配列データに対して、リファレンス配列またはコントロール配列から作成したクエリ配列セットとの比較を行い確認する工程をさらに含む、前記項目のいずれかに記載の方法。
(項目31) 前記対象配列データおよび前記コントロール配列データが、生物のゲノムに由来する塩基配列データであり、前記ゲノムの配列が不明である、前記項目のいずれかに記載の方法。
(項目32) 実験結果またはデータベースから対象配列データまたはコントロール配列データを取得する工程をさらに含む、前記項目のいずれかに記載の方法。
(項目X1) 対象配列データにおけるコントロール配列データに対する多型を含む部分配列中の多型ではない部分の少なくとも一部を含む配列を、該多型の識別子として割り当てることをさらに含む、前記項目のいずれか1項に記載の方法。
(項目X2) 前記多型の識別子をリファレンス配列にマッピングし、リファレンス上の該多型の位置を特定することを含む、前記項目のいずれかに記載の方法。
(項目33) 対象配列データにおいてコントロール配列データに対する多型を検出する方法をコンピュータに実行させるためのプログラムであって、該方法は、
a)該対象配列データの長さkの部分配列のサブセットをコンピュータに保存する工程であって、ここで、kは対象配列およびコントロール配列の全長以下の長さである、工程と、
b)該対象配列データの長さkのサブセットの各々の部分配列の出現頻度を算出する工程と、
c)該コントロール配列データの長さkの部分配列のサブセットにおける各々の部分配列の出現頻度をコンピュータに保存する工程と、
d)対象配列とコントロール配列とを比較し、該出現頻度の分布の差異に基づいて、多型を検出する工程と
を包含する、プログラム。
(項目33A) 前記項目のいずれか1つまたは複数に記載される特徴を有する、前記項目に記載のプログラム。
(項目34) 前記方法が、前記部分配列中の多型ではない部分の少なくとも一部を含む配列(前記部分配列全体であり得る。)を、検出された前記多型の名称として表示する工程をさらに含む、前記項目のいずれかに記載のプログラム。
(項目35) 対象配列データにおいてコントロール配列データに対する多型を検出する方法をコンピュータに実行させるためのプログラムを格納する記録媒体であって、該方法は、
a)該対象配列データの長さkの部分配列のサブセットをコンピュータに保存する工程であって、ここで、kは対象配列およびコントロール配列の全長以下の長さである、工程と、
b)該対象配列データの長さkのサブセットの各々の部分配列の出現頻度を算出する工程と、
c)該コントロール配列データの長さkの部分配列のサブセットにおける各々の部分配列の出現頻度をコンピュータに保存する工程と、
d)対象配列とコントロール配列とを比較し、該出現頻度の分布の差異に基づいて、多型を検出する工程と
を包含する、記録媒体。
(項目35A) 前記項目のいずれか1つまたは複数に記載される特徴を有する、前記項目に記載の記録媒体。
(項目36) 前記方法が、前記部分配列中の多型ではない部分の少なくとも一部を含む配列(前記部分配列全体であり得る。)を、検出された前記多型の名称として表示する工程をさらに含む、前記項目のいずれかに記載の記録媒体。
(項目37) 対象配列データにおいてコントロール配列データに対する多型を検出するためのシステムであって、該システムは、
該対象配列データおよび該コントロール配列データの長さkの部分配列のサブセットの各々の部分配列の出現頻度を提供するように構成された配列データ処理部であって、ここで、kは対象配列およびコントロール配列の全長以下の長さである、配列データ処理部と、
対象配列とコントロール配列とを比較し、該出現頻度の分布の差異に基づいて、多型を検出する工程とを行うように構成された、配列データ計算部と
を備える、システム。
(項目37A) 前記項目のいずれか1つまたは複数に記載される特徴を有する、前記項目に記載のシステム。
(項目38) 前記システムが、前記部分配列中の多型ではない部分の少なくとも一部を含む配列(前記部分配列全体であり得る。)を、検出された前記多型の名称として表示する表示手段をさらに含む、前記項目のいずれかに記載のシステム。
(項目39) 対象配列データにおいてコントロール配列データに対する多型を検出する方法であって、
(1)a)該対象配列データの長さkの部分配列のサブセットの各々の部分配列の出現頻度を提供する工程と、
b)該コントロール配列データの長さkの部分配列のサブセットの各々の部分配列の出現頻度を提供する工程と、
c)対象配列とコントロール配列とを比較し、該出現頻度の分布の比較に基づいて、多型を検出する工程と
によって、対象配列データにおける置換、コピー数多型、STRP、挿入、欠失、逆位または転座を検出するプロセスと、
(2)a)該対象配列データの配列中の少なくとも2ヶ所の部分配列の、該コントロール配列上の位置を特定する工程と、
b)対象配列データにおける該部分配列間の位置関係と、コントロール配列上の該部分配列間の位置関係とを比較する工程と、
c)対象配列データにおける該部分配列間の位置関係と、コントロール配列上の該部分配列間の位置関係とが異なっている場合、目的とする多型があると判定し、該対象配列データにおける該部分配列部位間の文字を、対応するコントロール配列上の文字と、該部分配列部位を始点として順次比較して不一致となる部位を検出する工程と
によって、対象配列データにおける挿入、欠失、逆位、転座または置換を検出するプロセスと
を包含する、方法。
(項目39A) 前記項目のいずれか1つまたは複数に記載される特徴を有する、前記項目に記載の方法。
(項目40) 対象配列データにおいてリファレンス配列データに対する多型を検出する方法であって、リファレンス配列データから、各々の位置情報と関連付けられたリファレンス配列のk長の部分配列セットを作成する工程を含み、さらに、
(A1)該対象配列データの長さkの部分配列のサブセットを生成し、ユニークな長さkの部分配列の出現頻度を提供する工程と、
(A2)該リファレンス配列のk長の部分配列セットの、ユニークな長さkの部分配列の出現頻度を提供する工程と、
(A3)該対象配列と該リファレンス配列とを比較し、該出現頻度の分布の比較に基づいて、挿入、欠失、置換、コピー数多型、STRP、逆位または転座を検出する工程とを包含するプロセスと
(B1)該対象配列データの配列中の少なくとも2ヶ所のk長の部分配列をクエリとして、該リファレンス配列のk長の部分配列セットに対して検索を行い、該少なくとも2ヶ所の部分配列の、リファレンス配列上の位置を特定する工程と、
(B2)該対象配列データにおける該部分配列間の位置関係と、該リファレンス配列上の該部分配列間の位置関係とを比較する工程と、
(B3)該対象配列データにおける該部分配列間の位置関係と、該リファレンス配列上の該部分配列間の位置関係とが異なっている場合、挿入、欠失、逆位または転座があると判定し、該対象配列データにおける該部分配列部位間の文字を、対応するコントロール配列上の文字と、該部分配列部位を始点として順次比較して不一致となる部位を検出する工程を包含し、必要に応じて、
(B4)該位置関係が異ならない場合に、該対象配列データにおける該部分配列部位間の文字を、対応する前記コントロール配列上の文字と比較して不一致となる部位を検出する工程をさらに含み、不一致となる部位が存在する場合、置換が存在すると判定する工程をさらに含む、プロセスと、
を、同時に、並行して、または逐次的に行うことを特徴とする、方法。
(項目40A) 前記項目のいずれか1つまたは複数に記載される特徴を有する、前記項目に記載の方法。
(項目A1) 対象配列データとコントロール配列データとの比較方法であって、
対象配列データにおけるコントロール配列データに対する多型を含む部分配列中の多型ではない部分の少なくとも一部を含む配列を、該多型の識別子として割り当てることを含む、方法。
(項目A1A) 前記項目のいずれか1つまたは複数に記載される特徴を有する、前記項目に記載の方法。
(項目A2) 前記多型の識別子をリファレンス配列にマッピングし、リファレンス上の該多型の位置を特定することを含む、前記項目のいずれかに記載の方法。
本発明において、上記1または複数の特徴は、明示された組み合わせに加え、さらに組み合わせて提供され得ることが意図される。本発明のなおさらなる実施形態および利点は、必要に応じて以下の詳細な説明を読んで理解すれば、当業者に認識される。
本発明により、2つ以上の配列の間で、全長配列における位置を考慮する必要なく、正確に多型、特に置換を検出することができる。k長配列を用いた置換変異の検出に関しては、ゲノムマッピングを行う前に多型検出ができることが1つの大きな特徴となる。そして、リファレンス配列が存在しない生物でも多型検出が可能であり、k-mer自体を多型の名称として用いることが可能であるため、それゆえ連鎖解析等の遺伝解析で大きな変革をもたらす可能性がある。
図1は、本発明の方法の1つの実施形態を例示するフロー図である。図1においては、本発明の方法において行われ得る主な工程の概要が示される。 図2は、本発明の方法の1つの実施形態を例示するフロー図である。図2においては、リファレンス配列の端からk-mer(k=20)を順に得て、対象配列および参照配列のk-mer(k=20)の頻度を位置情報と一緒に出力する工程の例が示される。このような出力を用いることにより、CNVや挿入・欠失、置換を検出することができる。 図3は、本発明の方法の1つの実施形態を例示するフロー図である。図3においては、リファレンス配列、コントロール配列(参照配列)および対象配列からのk長部分配列のサブセットおよび各々の部分配列の出現頻度を提供する工程の例が示される。 図4は、本発明の方法の1つの実施形態を例示するフロー図である。図4においては、対象配列とコントロール配列とを比較し、該出現頻度の分布の比較に基づいて、多型を検出する工程の例が示される。 図5は、本発明の方法の1つの実施形態を例示するフロー図である。図5においては、リファレンス配列における多型の位置を特定する工程の例が示される。 図6は、本発明の方法の1つの実施形態を例示するフロー図である。図6においては、検出された多型を確認する工程の例が示される。 図7は、コントロール配列(N1)と対象配列(N1S7)の間でのk-mer配列の比較結果の一部を示す。コントロール配列および対象配列のk長部分配列サブセットのそれぞれの部分配列と、当該部分配列の各々がコントロール配列および対象配列のk長部分配列サブセットに出現する頻度が示される。当該比較によって、コントロール配列由来の配列は対象配列では検出されず、逆に対象配列由来の配列は参照配列では検出されていないことから、多型を検出することができることが示される。下線の塩基がコントロールと対象とで異なっており、多型を示している。図7は、変異がホモ型の場合に観察される結果の一例である。図7における参照配列のk-mer部分配列は上からそれぞれ配列番号1~60に対応し、対象配列のk-mer部分は、上から配列番号1~20、61~80、40~60に対応する。 図8は、コントロール配列(N1)と対象配列(N1S7)の間でのk-mer配列の比較結果の一部を示す。コントロール配列および対象配列のk長部分配列サブセットのそれぞれの部分配列と、当該部分配列の各々がコントロール配列および対象配列のk長部分配列サブセットに出現する頻度が示される。当該比較によって、コントロール配列由来の配列は対象配列では検出されず、逆に対象配列由来の配列は参照配列では検出されていないことから、多型を検出することができることが示される。下線の塩基がコントロールと対象とで異なっており、多型を示している。図8は、変異がヘテロ型の場合に観察される結果の一例である。図8における参照配列のk-mer部分配列は上からそれぞれ配列番号81~140に対応し、対象配列のk-mer部分は、上から配列番号81~100、141~160、120~140に対応する。 図9は、コントロール配列(N1)と対象配列(N1S5、N1S6、N1S7、N1S10)のk-mer配列の部分配列サブセットを整列させることによって、イネリファレンス配列の対応する位置から始まる配列と対応する配列の出現頻度を求めた結果を示す図である。k-mer配列の出現頻度の変化から、N1S7においてヘテロ変異が生じ、N1S10においてホモ変異となったことを検出することができる。 図10は、コントロール配列(N1)と対象配列(N1S5、N1S6、N1S7、N1S10)との間で、Polymorphic Edge Detectionによって多型を検出・確認した結果の一部を示す。Chrは染色体番号、Posは染色体上の位置を示し、Refはリファレンス配列での当該位置の塩基、Altは変異が存在する場合の当該位置の塩基を示す。Polymorphic Edge Detectionの欄では、Control配列(N1)の20mer部分配列における最終塩基の出現頻度と、各対象(Target)配列の20mer部分配列における最終塩基の出現頻度が示される。最も高頻度の最終塩基が異なっている部分を、多型として検出している。Verifyの欄では、リファレンス配列から作成したクエリ配列セット(Ref:野生型クエリセット、Alt:変異型クエリセット)に対するControl配列データおよびTarget配列データの出現頻度が示される。最終的に判定した各部位における各試料の遺伝子型を、Genotypeの欄に示す(M:ホモ変異、H:ヘテロ変異、W:野生型)。 図11は、コントロール配列(N1)と対象配列(N1S5、N1S6、N1S7、N1S10)との間で、Polymorphic Edge Detectionによって多型を検出・確認した結果の一部を示す。Chrは染色体番号、Posは染色体上の位置を示し、Refはリファレンス配列での当該位置の塩基、Altは変異が存在する場合の当該位置の塩基を示す。Polymorphic Edge Detectionの欄では、Control配列(N1)の20mer部分配列における最終塩基の出現頻度と、各対象(Target)配列の20mer部分配列における最終塩基の出現頻度が示される。最も高頻度の最終塩基が異なっている部分を、多型として検出している。Verifyの欄では、リファレンス配列から作成したクエリ配列セットに対するControl配列データおよびTarget配列データの出現頻度が示される。最終的に判定した各部位における各試料の遺伝子型を、Genotypeの欄に示す(M:ホモ変異、H:ヘテロ変異、W:野生型)。 図12は、各世代のイネサンプルについてPCR法で対象となる多型部分を増幅して、サンガー法で塩基配列を確認した結果を示す図である。各試料の遺伝子型は、M:ホモ変異、H:ヘテロ変異、W:野生型として示される。Chrは染色体番号、Posは染色体上の位置を示し、Refはリファレンス配列での当該位置の塩基、Altは変異が存在する場合の当該位置の塩基を示す。N1、N1S1、N1S2、N1S3、N1S4、N1S5、N1S6、N1S7、N1S8、N1S9、N1S10と世代を重ねながら、ヘテロ変異が生じ、その後ホモ変異として定着する様子が観察されることが分かる。また、本発明の方法によって多型を検出した結果とよく一致しており、本発明の方法による多型検出が高い精度を有していることが理解される。 図13は、ナイジェリアのヨルバ族男性(NA18507)の配列データを用いて、リファレンスゲノム配列データから作成された参照(コントロール)配列に対する多型を検出した解析結果の一部を示す。対象配列データは、すでにIllumina社の次世代シーケンサーで解析されてNCBIに登録・公開されたデータをダウンロードして用いた。当該塩基配列セットの実験IDのURLは、https://www.ncbi.nlm.nih.gov/sra/SRX016231であり、配列のアクセッション番号は、SRR034939~SRR034975の範囲であった。k-1長の配列は、それぞれ上から配列番号161~190に対応する。 図14は、本発明において行われ得る確認工程における、クエリ配列セットの作成の模式図である。上部の配列がリファレンス配列であり、下線太字で示されるTが検出された多型部位を示す。当該多型部位を含む部分配列のセットを生成し、クエリ配列セットとすることができる。各配列は、それぞれ上から配列番号267~275に対応する。 図15Aは、本発明のシステムの実施形態を模式的に示した図である。 図15Bは、本発明のシステムのさらなる実施形態を模式的に示した図である。 図16は、本発明の方法の実施形態を模式的に示した図である。 図17は、本発明の方法によるコピー数多型(CNV)の検出の結果の一部を示す図である。イネ第7染色体の26694795位置(図中、矢印で示す)からコピー数多型部位が開始されている。検出されたコピー数多型部位は、レトロトランスポゾンTos17に対応する。このトランスポゾンは4.1kbあるため、図17には最初のジャンクションの部分のみ示している。培養時間に応じた転移によるコピー数の増加が本発明の方法によって検出されていることが理解される。 図18は、k-mer配列の頻度を用いる多型検出フローと、部分配列の位置関係を用いる多型検出フローとを組み合わせて行う場合の一実施形態を示すフロー図である。
以下、本発明を最良の形態を示しながら説明する。本明細書の全体にわたり、単数形の表現は、特に言及しない限り、その複数形の概念をも含むことが理解されるべきである。従って、単数形の冠詞(例えば、英語の場合は「a」、「an」、「the」など)は、特に言及しない限り、その複数形の概念をも含むことが理解されるべきである。また、本明細書において使用される用語は、特に言及しない限り、当該分野で通常用いられる意味で用いられることが理解されるべきである。したがって、他に定義されない限り、本明細書中で使用される全ての専門用語および科学技術用語は、本発明の属する分野の当業者によって一般的に理解されるのと同じ意味を有する。矛盾する場合、本明細書(定義を含めて)が優先する。
(定義)
以下に本明細書において特に使用される用語の定義および/または基本的技術内容を適宜説明する。
本明細書において、「配列」とは、各々が何らかの値を取る複数の変数であって、それら複数の変数の順序の情報をさらに含むものをいう。代表的には文字列で表示される。
本明細書において、「対象配列」とは、多型を検出しようとする任意の配列をいい、本明細書においては、「ターゲット」、「ターゲット配列」、「target」とも表記する場合がある。
本明細書において、「コントロール配列」とは、その配列との差異を多型として検出するための基準として用いられる任意の配列をいい、本明細書においては、「コントロール」、「参照配列」、「比較配列」、「control」とも表記する場合がある。
本明細書において、「多型(polymorphism)」とは、対象配列中においてコントロール配列と異なっている任意の部分を指す。本明細書において、「変異」も同様の意味で使用することができる。
本明細書において、「リファレンス(reference)配列」とは、対象配列および/またはコントロール配列の全長の配列として扱うことができる配列を指す。いかなる配列を全長配列とするかは、対象配列および/またはコントロール配列として用いる配列に応じて適宜決定されるものであり、例示されるものに限定されないが、例えば、ウェブ上のデータベース等に存在する、全ゲノム配列、染色体全長配列、遺伝子全長配列、プラスミド全長配列、エクソン全長配列、タンパク質全長配列などをリファレンス配列として用いることができる。
本明細書において、「配列データ」とは、ある配列についての情報を与えるデータをいう。代表的には、配列そのものも配列データということができ、また、配列の一部について情報を与えるデータ(例えば、ゲノム配列に対するシーケンシングによる解析データ)も配列データとして包含される。
本明細書において、ある配列の「部分配列」とは、その配列に含まれる任意の配列をいう。
本明細書において、「サブセット」とは、配列の集合と、それらの配列の部分配列の集合とを合わせた集合の任意の部分集合をいう。
本明細書において、「次世代シーケンシング」とは、配列決定プロセスを並列化し、一度のランで数千万から数億の配列データを生成するシーケンシング技法である。「次世代シーケンサー」とは、次世代シーケンシングを行うための機器を指す。
「偶然同一を排除する」とは、ある配列と、偶然に同一の配列が出現する期待値を1未満にすることをいう。
本明細書において、「カバレッジ」とは、配列データの量が、配列全長の何倍に相当しているかを指す。「カバー率」、「~倍の読み」などと称される場合もある。
本明細書において、「配列構造体」とは、配列中における、物理的に分離された一連の配列をいう。例えば、ゲノム配列の文脈では、染色体のそれぞれは配列構造体ということができる。
本明細書において、「転座」とは、複数の配列構造体を有する配列中で、ある配列構造体上の部分配列が、他の配列構造体上に移動している多型をいう。
本明細書において、「ジャンクション」とは、一部が同一である2つの配列について、同一である部分と同一でない部分の境界を指す。
本明細書において、「識別子」とは、ある多型を他の多型と区別するために付される名称を指す。一般的には、多型の開始位置と型で記載されることが多いが、本明細書において記載される識別子を用いることができる。
本明細書において、「エッジ」とは、配列において多型を含む部分の末端をさす。
(好ましい実施形態)
以下に本発明の好ましい実施形態を説明する。以下に提供される実施形態は、本発明のよりよい理解のために提供されるものであり、本発明の範囲は以下の記載に限定されるべきでないことが理解される。従って、当業者は、本明細書中の記載を参酌して、本発明の範囲内で適宜改変を行うことができることは明らかである。また、本発明の以下の実施形態は単独でも使用されあるいはそれらを組み合わせて使用することができることが理解される。
なお、以下で説明する実施の形態は、いずれも包括的または具体的な例を示すものである。以下の実施の形態で示される数値、形状、材料、構成要素、構成要素の配置位置及び接続形態、ステップ、ステップの順序などは、一例であり、請求の範囲を限定する主旨ではない。また、以下の実施の形態における構成要素のうち、最上位概念を示す独立請求項に記載されていない構成要素については、任意の構成要素として説明される。
(本発明の多型検出の概要)
本発明は、対象配列データにおいてコントロール配列データに対する多型を検出する方法を提供する。この方法は、a)該対象配列データの長さkの部分配列のサブセットの各々の部分配列の出現頻度を提供する工程と、b)該コントロール配列データの長さkの部分配列のサブセットの各々の部分配列の出現頻度を提供する工程と、c)対象配列とコントロール配列とを比較し、該出現頻度の分布の比較に基づいて、多型を検出する工程とを包含し、ここで、kは該対象配列および該コントロール配列のいずれか短いほうの全長以下の整数である、方法を提供する。本発明の例示的なフローは図16に説明されている。
本発明の方法は、対象配列データと、コントロール配列データとの2つの配列データ(例えば、次世代シーケンサー解析結果)からの直接比較により多型を検出することが可能である点で、従来法と根本的に異なる。とりわけ、塩基配列における多型検出において、ゲノム上の位置の考慮なしに直接比較する方法は、新規なものであると考えられる。
1つの実施形態において、本発明の方法は、配列データから、一定長(k長)の部分配列のセットを得ることを1つの特徴とする。1つの実施形態において、本発明の方法は、配列データから、k長の部分配列のセットにおける各部分配列の頻度分布を得ることを1つの特徴とする。一部の実施形態では、配列データから、配列を1つずつずらしながらk長の部分配列のセットを作成する。
1つの実施形態では、k長の配列のうち、k-x(xは1など)の配列が同じであるデータをソートし、異なる部分(x長の部分に該当)の頻度を検出する。X長の部分は、部分配列中での位置は限定されず、配列中の中央部でもよい。しかしながら、X長の部分を、部分配列の末端(例えば、塩基配列においては3’末端または5’末端)にすることは、ソートなどの処理を顕著に簡便化・高速化するため、好ましい。このように、k長の配列のうち末端部で異なる部分(x長の部分に該当)の頻度を検出する場合、本明細書において、「Polymorphic Edge Detection(PED)」または「edge detection」と呼称される場合がある。
ここで、kの値としては、配列データの各配列(例えば、次世代シーケンサーの各々のショートリード)の長さを上限とした任意の値を挙げることができ、例えば、約500、約400、約300、約200、約100、約50、約40、約30、約25、約20、約15等を挙げることができる。kが増加することにより、k-mer配列のデータは指数的に増加する(例えば、塩基配列の場合、kが1塩基増えるごとに塩基の組み合わせは4倍になる)ため、例えば、塩基配列の場合、k=20~25程度が好ましいが、理論上は、例えば、k=500等でも用いることが可能であり、制限されるものではない。ヒトの場合であると、k=17以下だと偶然一致が生じる確率が高くなるが、ゲノムサイズが小さな生物であれば、例えば、k=15などより小さなk値を用いることも可能である。1つの実施形態ではk=20を用いる。
k-x長の配列が同じであるデータのxの部分の文字が比較対象間で異なる場合、その文字に多型(変異)が含まれると考えられる。挿入・欠失変異も変異の末端文字を検出できる。例えば、k-x長の塩基配列が同じであるデータのxの部分の塩基が比較対象間で異なる場合、その塩基に多型(変異)が含まれると考えられる。
一部の実施形態において、得られた配列セットのうち、同じ配列に関して出現回数で整理したデータを算出する。この工程は、計算機を使用して簡便に行うことができ、例えば、Unixでは以下:
Figure 0007166638000001

のように実装することができ、(文字列で)ソートされた配列と頻度を示す数値のデータを生成することができる。さらに、対象とコントロールの頻度データを同一k-merでまとめる際に、例えば、Unixではjoinコマンド等を使用して行うことが可能である。
配列出現回数を、配列データのカバレッジ(何倍読みか)と比較することによって、配列データにおける差異を評価することも可能である。例えば、ゲノム配列に対して40倍のデータ量の配列解析からの配列データに対して、出現頻度が1となるようなものはノイズと考えることができる。
本発明は、特に、「置換」多型(長さが変わらない=欠失挿入ではない)の検出には、極めて高い効果を発揮する。多型部位が多コピーの場合は位置の特定ができない場合があり得るが、それでもなお、多型自体の検出は可能であり、多型に名称を付して特定できる。そのため、例えば、検出された各多型を、形質との関係を調べる多型マーカーとして使用することが可能であり、診断・育種・鑑定・品質管理(例えば、iPS細胞の品質管理)・分類・検査にも応用可能である。
本発明を、次世代シーケンサーから得られた塩基配列データから直接多型を検出する方法として用いることで、2種類のサンプル間、およびリファレンス配列とサンプルとの間の多型の検出が可能となる。また、長さkの部分配列、k-xの部分配列は重複のないユニークな配列なので、配列自体を多型の識別子(名称)として利用可能である。このため、リファレンスゲノム配列が決定されていないため多型のゲノム上の位置関係が判別できない場合においても、世界共通の一意の多型の識別子として利用できる。本発明の1つの実施形態において、対象配列データとコントロール配列データとの比較方法であって、対象配列データにおけるコントロール配列データに対する多型を含む部分配列中の多型ではない部分の少なくとも一部を含む配列を、多型の識別子として割り当てることを含む、方法が提供される。また、多型の識別子をリファレンス配列にマッピングし、リファレンス上の多型の位置を特定することが可能である。
識別子は、多型ではない部分の少なくとも一部に加えて、多型自体も含み得る。多型塩基を含めた識別子は、リファレンス配列へのマッピングは難しいものの、連鎖解析に用いることが可能である。
例えば、
AAACCACTTCACGTTTCCA A
AAACCACTTCACGTTTCCA G
という多型の例では、
AAACCACTTCACGTTTCCAのA型
AAACCACTTCACGTTTCCAのG型
AAACCACTTCACGTTTCCAのA/Gのヘテロ型
という表現が記載の一例である。
多型を含めた表記の仕方の例としては、
AAACCACTTCACGTTTCCAA型、
AAACCACTTCACGTTTCCAG型、
そして、ヘテロ型は、
AAACCACTTCACGTTTCCAA/AAACCACTTCACGTTTCCAG
のように、2つの型を併記することが可能である。
本発明の1つの実施形態は、対象配列データにおいてコントロール配列データに対する多型を検出する方法である。1つの実施形態において、当該方法は、該対象配列データの長さkの部分配列のサブセットの各々の部分配列の出現頻度を提供する工程を含む。kは該対象配列および該コントロール配列のいずれか短いほうの全長以下の整数である。1つの実施形態において、当該方法は、該コントロール配列データの長さkの部分配列のサブセットの各々の部分配列の出現頻度を提供する工程を含む。1つの実施形態において、当該方法は、対象配列とコントロール配列とを比較し、該出現頻度の分布の比較に基づいて、多型を検出する工程とを含む。このような工程によって、全長配列における位置を考慮せず、また、配列を連結することなく配列データを比較し、多型を検出することができる。
さらなる実施形態において、本発明の方法は、部分配列中の長さk-xの配列部分が共通する配列ごとに、長さxの部分について出現頻度の分布を算出する工程をさらに含む。xはk未満の正の整数である。この場合において、出現頻度の分布の比較として、前記部分配列中の長さk-xの配列部分が共通する配列における、長さxの部分の出現頻度の分布の差異の比較が含まれ得る。これにより、多型検出の処理を高速化することが可能である。
一部の実施形態において、本発明の方法は、前記部分配列中の長さk-xの配列部分を、ユニークな配列ごとにグルーピングする工程を含む。これには、例えば、前記長さk-xの配列部分をソートする工程(例えば、前記長さk-xの配列部分を文字列としてソートする工程)が含まれ得る。
一部の実施形態において、kの値は、前記対象配列データ等における偶然同一を排除する長さである。例えば、前記対象配列データおよび前記コントロール配列データが、生物のゲノムに由来する塩基配列データである場合、前記kは、前記生物のゲノムにおいて、異なる箇所での偶然同一を排除する長さであり得る。これにより、多型の検出をより正確なものとすることが可能である。
長さxは、限定されるものではないが、好ましくは1~3であり、さらに好ましくは1~2であり、より好ましくは1である。
1つの実施形態では、前記長さxの部分が、前記部分配列の末端に存在する。例えば、前記対象配列データおよび前記コントロール配列データが塩基配列データである場合、前記長さxの部分は、前記部分配列の3’末端または5’末端であり得る。長さxの部分を部分配列の末端にとることは、比較処理の高速化・簡便化にとって望ましい。
出現頻度の分布の差異の比較により、例えば、以下のような多型の検出が可能である。1つの実施形態では、前記コントロール配列データのサブセットと前記対象配列データのサブセットとの間で、前記長さxの部分の配列の出現頻度が異なる場合、該長さxの部分の配列を、対象配列データにおけるコントロール配列データに対する多型として検出する。1つの実施形態では、前記コントロール配列データのサブセットと前記対象配列データのサブセットとの間で、前記長さxの部分の配列で最も高頻度のものが異なっている長さk-xの配列部分が存在する場合、該長さxの部分の配列を、対象配列データにおける多型として検出する。
1つの実施形態では、前記対象配列データおよび前記コントロール配列データが塩基配列データであり、前記対象配列データのサブセットにおける前記長さxの部分の配列で、前記コントロール配列データのサブセットにおける最も高頻度のものと同一の長さxの部分の配列がノイズ以下のカウントしか存在しない長さk-xの配列部分が存在する場合、該長さxの部分の配列を、対象配列データにおけるホモ多型として検出する。ノイズは、本明細書において後述されるような基準を用いて判定することができる。1つの実施形態では、前記対象配列データおよび前記コントロール配列データが塩基配列データであり、対象サブセットにおける前記長さxの部分の配列で、コントロール配列データのサブセットにおける最も高頻度のものと同一の長さxの部分の配列が存在し、かつ、コントロール配列データのサブセットにおける最も高頻度のものと異なる長さxの部分の配列が存在する長さk-xの配列部分が存在する場合、該長さxの部分の配列を、対象配列データにおけるヘテロ多型として検出する。
一部の実施形態において、本発明の方法において、予測される出現頻度と比較して、出現頻度が少ない部分配列をノイズとすることが可能である。これにより、機械的に生じた差異と、実際に生じている多型とを識別して検出することが可能である。例えば、前記対象配列データおよび前記コントロール配列データが塩基配列データである場合には、対象配列データ量、配列データにおける予想されるエラー率(例えば、マニュアル・機器の公称値)、対象配列全長サイズ(例えば、ゲノムサイズ)等を考慮して、予想されるノイズのカウント程度、または予想されるノイズのカウント以下の出現頻度の部分配列をノイズとすることができる。1つの予測としては、生じるノイズの頻度の期待値は、int(ゲノムを何倍読んだか×(1-精度)+1)」となる。エラーが存在しない場合、精度は1になる。1つの実施形態では、対象配列データおよび前記コントロール配列データが塩基配列データである場合に、[(対象配列データ量)×(1-精度)]/(対象ゲノムサイズ)+1未満の出現頻度の部分配列をノイズとすることが可能である。
例えば、エラー率が0.001の場合は着目した塩基に関して1カウントでも出現する確率は0.001になるためほぼ0と考えられるが、1000塩基の範囲で見れば、どれか1塩基にエラーが含まれる計算となる。そのため、繰り上げた整数値をノイズの予測値とすることができると考えられ、int関数中で1を加えることによって繰り上げた整数値とすることができる。
あるいは、本発明の方法において、予測される出現頻度と比較して、出現頻度が多い部分配列を、リピート配列部位として除外することができる。例えば、対象配列データの対象配列全長のカバレッジ(カバー率)の2倍を超えるような部分配列を解析から除外することができる。
1つの実施形態では、前記対象配列データおよび/またはコントロール配列データが、次世代シーケンシングによって得られた塩基配列データである。次世代シーケンシングにおける多型の検出においては、従来、リファレンスへのマッピングや、配列のアセンブリが必要とされており、このような工程で生じる不確実性によって、多型の検出が大きく阻害されていたため、このような工程を必要としない本発明の方法を次世代シーケンシングから得られた配列データに対して用いることは特に有利なものであるということができる。
対象配列データおよびコントロール配列データは、限定されるものではないが、多型を検出する上では、一定の共通性を持つ配列についての配列データであることが望ましい。しかしながら、配列の取得方法については各々同一でも異なっていてもよく、シーケンシングによって得られたデータ間での比較を行うことも、データベース等から得られたデータ間での比較を行うことも、シーケンシングによって得られたデータとデータベース等から得られたデータとの間での比較を行うことも可能である。
1つの実施形態では、対象配列データが、個体から得られた配列データであり、コントロール配列データが、該個体と同種の別の個体、またはデータベースから得られた配列データである。1つの実施形態では、対象配列データが、個体の組織試料から得られた配列データであり、コントロール配列データが、該個体の別の組織、またはデータベースから得られた配列データである。1つの実施形態では、対象配列データが、細胞試料から得られた配列データであり、コントロール配列データが、別の細胞、またはデータベースから得られた配列データである。
本発明の方法は、全長配列の情報を必要としないため、データベース等において例えば全長配列が公知でない場合にも用いることができ、例えば、対象配列データおよびコントロール配列データは、生物のゲノムに由来する塩基配列データである場合、前記ゲノムの配列が不明であってもよい。
本発明の方法によって検出できる多型としては、置換、挿入、欠失、コピー数変異(Copy Number Variation, CNV)、STRP(short tandem repeat polymorphism)、逆位および転座が挙げられる。1つの実施形態において、本発明の方法は、上記複数の多型の任意の組み合わせを同時に検出することも可能である。さらなる実施形態において、本発明の方法は、上記複数の多型の全てを同時に検出することも可能である。特に、多型が置換である場合には、本発明の方法は、非常に高い検出力を発揮することが可能である。
対象配列に対するリファレンス配列が存在する場合、本発明の方法は、前記対象配列に対するリファレンス配列における前記多型の位置を特定する工程をさらに含むことができる。例えば、対象配列データおよびコントロール配列データが、生物のゲノムに由来する塩基配列データである場合、多型のゲノム上の位置を特定する工程をさらに含むことができる。この位置の特定は、本発明の方法が、多型を周囲の配列と関連づけて検出する(例えば、x長部分の多型がk-x長の配列と関連付けられる)ことを可能にしているため、リファレンス配列に対して検索を行うことにより、簡便に行うことが可能である。
本発明の方法は、検出した多型について確認する工程をさらに含むことができる。確認は、例えば、検出された多型の部位について、リファレンス配列またはコントロール配列から作成したクエリ配列セットを用いて、対象配列データおよび/またはコントロール配列データとの比較を行うことによって行うことができる。クエリ配列セットは、リファレンス配列またはコントロール配列において前記多型に該当する部位の文字を異なる文字に置換した変異型クエリ配列セット、および/またはリファレンス配列またはコントロール配列において前記多型に該当する部位の文字を置換していない野生型クエリ配列セットを含み得る。
本発明の方法は、対象配列データおよびコントロール配列データが塩基配列データである場合、検出された多型の部位について、対象配列データおよび/またはコントロール配列データの相補鎖の配列データに対して、リファレンス配列またはコントロール配列から作成したクエリ配列セットとの比較を行い確認する工程をさらに含むことができる。本発明の方法は、対象配列データおよびコントロール配列データが塩基配列データである場合、検出された多型の部位について、対象配列データおよび/またはコントロール配列データの対立遺伝子の配列データに対して、リファレンス配列またはコントロール配列から作成したクエリ配列セットとの比較を行い確認する工程をさらに含むことができる。ここで、対立遺伝子の配列データとして、実際の遺伝子の存在の有無とは関係なく、野生型に対する変異型の塩基を有する配列データを用いることができる。
本発明の方法は、実験結果またはデータベースから対象配列データまたはコントロール配列データを取得する工程を含んでもよい。また、本発明の方法においては、必ずしも配列データそのものを取得する必要はなく、配列データのサブセット、および/または配列データもしくは配列データのサブセットにおける頻度分布のデータを取得して実行することも可能である。
1つの局面において、本発明は、本発明の多型を検出する方法をコンピュータに実施させるための方法を実装するプログラム、該プログラムを記録した記録媒体、およびこれを実現するためのシステムを提供する。ここで採用され得る任意の特徴は本明細書の多型を検出する方法の説明に記載される任意の特徴またはその組み合わせを採用することができる。
したがって、1つの実施形態において、対象配列データにおいてコントロール配列データに対する多型を検出する方法をコンピュータに実行させるためのプログラムであって、該方法は、
a)該対象配列データの長さkの部分配列のサブセットをコンピュータに保存する工程であって、ここで、kは対象配列およびコントロール配列の全長以下の長さである、工程と、
b)該対象配列データの長さkのサブセットの各々の部分配列の出現頻度を算出する工程と、
c)該コントロール配列データの長さkの部分配列のサブセットにおける各々の部分配列の出現頻度をコンピュータに保存する工程と、
d)対象配列とコントロール配列とを比較し、該出現頻度の分布の差異に基づいて、多型を検出する工程と
を包含する、プログラムが提供される。さらなる実施形態において、上記方法が、前記部分配列中の多型ではない部分の少なくとも一部を含む配列(前記部分配列全体であり得る)を、検出された前記多型の名称として表示する工程をさらに含む、プログラムが提供される。
別の実施形態において、対象配列データにおいてコントロール配列データに対する多型を検出する方法をコンピュータに実行させるためのプログラムを格納する記録媒体であって、該方法は、
a)該対象配列データの長さkの部分配列のサブセットをコンピュータに保存する工程であって、ここで、kは対象配列およびコントロール配列の全長以下の長さである、工程と、
b)該対象配列データの長さkのサブセットの各々の部分配列の出現頻度を算出する工程と、
c)該コントロール配列データの長さkの部分配列のサブセットにおける各々の部分配列の出現頻度をコンピュータに保存する工程と、
d)対象配列とコントロール配列とを比較し、該出現頻度の分布の差異に基づいて、多型を検出する工程と
を包含する、記録媒体が提供される。さらなる実施形態において、方法が、前記部分配列中の多型ではない部分の少なくとも一部を含む配列(前記部分配列全体であり得る)を、検出された前記多型の名称として表示する工程をさらに含む、記録媒体が提供される。
別の実施形態において、対象配列データにおいてコントロール配列データに対する多型を検出するためのシステムであって、該システムは、
該対象配列データおよび該コントロール配列データの長さkの部分配列のサブセットの各々の部分配列の出現頻度を提供するように構成された配列データ処理部であって、ここで、kは対象配列およびコントロール配列の全長以下の長さである、配列データ処理部と、
対象配列とコントロール配列とを比較し、該出現頻度の分布の差異に基づいて、多型を検出する工程とを行うように構成された、配列データ計算部と
を備える、システムが提供される。さらなる実施形態において、前記部分配列中の多型ではない部分の少なくとも一部を含む配列前記部分配列全体であり得る)を、検出された前記多型の名称として表示する表示手段をさらに含む、システムが提供される。
(具体的な例)
本発明の例示的な実施形態は、以下のような工程による方法である。
1.配列データについて、長さkの部分配列のサブセットを得る。例えば、塩基配列データを端から1塩基ずつずらしながら、k長の配列のセットを得る。以下の例では、次世代シーケンサーで得られたイネ個体の配列データ(データ量はイネゲノムの40倍、リード長は100塩基)を用いて、k長を20塩基として得られた配列セットを用いて説明する。得られた配列セットに関して、比較対象のサンプル間で5’末端からk-1長の配列が同じであるデータの3’末端の塩基が比較対象の間で異なる場合、その塩基は多型である。
Figure 0007166638000002

(上記例において、1行目から各々配列番号191~201である。)
最初の行が次世代シーケンサーから得られた塩基配列であり、以下の行がk長の部分配列を示す。この場合k長を20塩基で配列セットを得ている。
2.得られた部分配列セットの出現頻度データを得る。すなわち、部分配列セットのうちの同じ配列に関して、出現回数で整理したデータを得る。
Figure 0007166638000003

(上記例において、1行目から各々配列番号202~211である。)
配列セットを降順に整列して出現回数を配列の右に表示している。この例では、ゲノムの40倍の解析なので、出現回数が1回程度の配列はノイズであると判断することができる。40~50回程度の配列は、ゲノム上のユニークな配列由来であり、89回出現した配列は、ゲノム上2ヶ所存在していると考えられる。
3.部分配列中の長さk-xの配列部分が共通する配列ごとに、長さxの部分について出現頻度の分布を算出する。この例では、得られたk塩基の部分配列(k=20)の頻度データをもとに、最初の19塩基(k-x、x=1)に対する最後の1塩基のACGTそれぞれの塩基の出現頻度データに変換する。
Figure 0007166638000004

(上記例において、1行目から各々配列番号212~220である。)
(20塩基の頻度データを最初の19塩基と最後のACGTの塩基の頻度一覧に変換する。)
4.部分配列中の長さk-xの配列部分が共通する配列における、長さxの部分の出現頻度の分布の差異を比較する。例えば、コントロールと調べたい対象の配列データ(ターゲット)からの頻度データを19merの配列でまとめた一覧を作成する。多型がない場合、最後の塩基の頻度は同じ塩基で最多となる。最後の塩基の頻度でコントロールと対象で塩基が異なる場合は、その塩基が多型である。
Figure 0007166638000005

(上記例において、1行目から各々配列番号213、215、217、および218である。)
コントロールと調べたい対象の最後の塩基の出現頻度の一覧。この場合、AAAAGATCTATGAGCACTC(配列番号218)の次にはコントロールではAのみであるが、対象ではAとGのヘテロザイガスであることがわかる。また、ホモザイガスの多型が生じる場合には、以下のように、最後の塩基として出現するものが異なり、検出することができる。
Figure 0007166638000006
このように、ゲノム上の位置が不明でも、最初の19塩基で表記される配列位置に続く塩基の多型を網羅することができる。ゲノムリファレンス配列が既知の場合、最初の19塩基に対応するゲノム位置から多型の位置を決めることができる。順鎖、相補鎖ともに同じ位置に検出された場合は一塩基多型である。多型となる配列の末端部分を検出するため、挿入・欠失多型の検出も可能である。最初の19塩基を多型の名称、最後の塩基を遺伝子型として表記するため、その多型を示す唯一の表記法として利用できる。あるいは、多型も含めてk-mer(例えば、20塩基)全体で一つの名称とすることも可能である。例えばk=5の場合には、ACGTA型とACGTT型といった表記が可能である。本発明の方法は、図1に示されるようなフローにしたがって、図1に示される工程を適宜採用することによって実行することが可能である。
(k-mer)
1つの実施形態において、本発明の1つの特徴は、対象配列データにおけるコントロール配列データに対する多型の検出において、該対象配列データの長さkの部分配列のサブセット、または該部分配列のサブセットの各々の部分配列の出現頻度を用いることである。ここで、kは対象配列およびコントロール配列のいずれか短いほうの全長以下の整数である。
長さkの部分配列は、対象配列データ、コントロール配列データ、リファレンス配列データ等から切り出すことによって生成することができる。例えば、一定間隔でk長の配列を切り出すことによって生成することができ、1文字ごと、2文字ごと、3文字ごと、またはそれ以上の間隔で切り出して部分配列セットを生成することができる。あるいは、対象配列データ、コントロール配列データ、リファレンス配列データ等から、ランダムに位置を選択して切り出すことも可能である。全てのk長部分配列を生成することが望ましい場合には、1文字ごとに切り出し位置をずらしながらk長部分配列のサブセットを生成することができる。
長さkは、対象配列、コントロール配列および/またはリファレンス配列における偶然同一を排除する長さであることが望ましい。偶然同一を排除することによって、異なる配列が対象配列の別の箇所に偶然含まれていたものをコントロール配列との差異として検出する可能性を低減し、より正確に解析することができる。k長のある配列と、偶然に同一の配列が、対象配列、コントロール配列および/またはリファレンス配列に出現する期待値を1未満にすることが望ましい。
一般的に、v:配列において各変数が取り得る値の種類、L:解析対象とする配列(対象配列、コントロール配列および/またはリファレンス配列)の全長(含まれる変数の数)として、v^k>Lとなる場合に、kが偶然同一を排除する長さであると考えられる。あるいは、配列全長が10^Lのようなオーダーで表される場合、両辺の対数をとり、k>L/log(v)を満たす場合に、kが偶然同一を排除する長さであると考えられる。
例えば、配列データが塩基配列データである場合、4つの文字が変数の値として考えられるため、v=4である。例えば、ヒトゲノムは3.1×10の9乗のサイズであり、仮に10の9乗長のアトランダムな塩基配列があった場合、9/log(4)≒15塩基が偶然一致を排除できる最小k長と考えられる。ゲノムサイズが異なる生物に対しても、例えば、10/log(4)の場合はk=17で偶然一致を排除できると考えられる。長いほど偶然一致の配列が生じる可能性を減少させることができるが、データサイズがその分大きくなる。
(頻度)
1つの実施形態において、本発明の方法は、対象配列データおよび/またはコントロール配列データの長さkの部分配列のサブセットの各々の部分配列の出現頻度を提供する工程を含む。
長さkの部分配列のサブセットの各々の部分配列の出現頻度は、ユニークなk長の部分配列データに対する頻度がペアになったデータとして(例えば、1列目に部分配列、2列目にその頻度が提供される行列データとして)提供されてもよい。部分配列の出現頻度は、配列とその頻度とが関連付けられているものであれば、その形式は限定されるものではなく、対応する識別子を有する配列データと頻度データとを別個に出力することも可能である。また、部分配列の出現頻度は、既に存在するデータを取得することによって提供することも可能である。
長さkの部分配列のサブセットの各々の部分配列の出現頻度を得る工程は、計算機を使用して簡便に行うことができ、例えば、Unixでは以下:
Figure 0007166638000007

のように実装することができる。
多型が存在している部位がゲノム上でユニークな配列部位であり、例えば1塩基置換が存在している場合、その置換を含むk-mer部分配列は対象配列データのサブセットには存在するが、コントロール配列データのサブセットには存在しないと考えられる。当該置換変異を含むk-merではすべてこの結果が得られると考えられ、結果的に当該置換変異部位を含む2×k-1個のk-mer部分配列で出現頻度の差異が観察され、当該多型を検出することができる。
一部の実施形態において、本発明の方法において、予測される出現頻度と比較して、出現頻度が少ない部分配列をノイズとすることが可能である。これにより、機械的に生じた差異と、実際に生じている多型とを識別して検出することが可能である。例えば、前記対象配列データおよび前記コントロール配列データが塩基配列データである場合には、対象配列データ量、配列データにおける予想されるエラー率(例えば、マニュアル・機器の公称値)、対象配列全長サイズ(例えば、ゲノムサイズ)等を考慮して、予想されるノイズのカウント程度、または予想されるノイズのカウント以下の出現頻度の部分配列をノイズとすることができる。1つの予測としては、生じるノイズの頻度の期待値は、int(ゲノムを何倍読んだか×(1-精度)+1)」となる。
理論的には、シーケンサーの精度が99%でゲノムの100倍読みした場合は、1塩基のノイズが入ることになる。通常の解析はゲノムの40倍読み程度であるため、ノイズは1塩基以下と考えられるが、カウントは整数値となるため、リードエラーは1とカウントされる。すなわち、2塩基以上のカウントはノイズ以外の何らかの要因がある可能性が高いと考えられる。ただし、低い確率だが、同一塩基に2回以上ノイズが入る可能性はある。予想されるノイズのカウントは、例えば、int(ゲノムを何倍読んだか×(1-精度)+1)である。この場合の精度は100%の場合1、99%の場合0.99として計算できる。なお、int(X)は、X未満の最大の整数を返す関数である。
エラー率(精度)が完全に0でない場合には、例えば0.00001でも10万か所のうち1ヶ所は1のカウントが生じると考えられる。そのため、この場合、結局1はノイズである可能性が高いと考えられるが、2はノイズにしては高すぎると言える。そのため、int(ゲノムを何倍読んだか×(1-精度)+1)で繰り上げて整数にすることにより、予測値を算出することが可能である。1つの実施形態では、カウントが1となる部分配列をノイズと判定する。例えば、対象配列データおよび前記コントロール配列データが塩基配列データである場合に、[(対象配列データ量)×(1-精度)]/(対象ゲノムサイズ)+1未満の出現頻度の部分配列をノイズとすることにより、予想されるノイズのカウント以下の出現頻度の部分配列をノイズとして排除することができる。
頻度の算出から、コピー数変異(CNV)の検出も可能である。例えば、コントロール配列と比較して約2倍以上の頻度が連続して検出された場合にCNVと判定することができる。頻度の基準は、約2.5倍以上などと厳しくすることによって、誤検出を減少させることもできる。
あるいは、本発明の方法において、予測される出現頻度と比較して、出現頻度が多い部分配列を、リピート配列部位として除外することができる。例えば、対象配列データの対象配列全長のカバレッジ(カバー率)の2倍を超えるような部分配列を解析から除外することができる。
例えば、配列データ量の配列全長に対するカバレッジと比較し、同程度の出現頻度の部分配列は、対象配列の全長におけるユニークな配列に由来している配列と考えられる。さらに、配列データ量の配列全長に対するカバレッジと比較し、2倍程度の出現頻度の部分配列は、対象配列の全長において2箇所存在する配列に由来している配列と考えられる。
kが偶然同一を排除している場合、それを上回る場合には、部分配列が、リピート配列部位に由来しているものと判断することができる。例えば、ゲノムの40倍のカバレッジのシーケンシングデータにおいて、出現頻度が40~50回程度の配列は、ゲノム上のユニークな配列由来であり、89回出現した配列は、ゲノム上2ヶ所存在していると考えられる。k長配列のサブセットおよび/または各部分配列の出現頻度データの作成は、一例としては、図3に例示されるような工程を採用して行うことができる。
(k-x)
1つの実施形態では、本発明の方法は、k長の部分配列中の長さk-xの配列部分が共通する配列ごとに、長さxの部分について出現頻度の分布を算出することを特徴とする。xはk未満の正の整数である。長さkの全長ではなく、k-merの一部(x長の部分)の文字の差異を見ることによって、計算量を顕著に減少させることが可能である。長さxは、限定されるものではないが、好ましくは1~2であり、より好ましくは1である。
一部の実施形態において、本発明の方法は、前記部分配列中の長さk-xの配列部分を、ユニークな配列ごとにグルーピングする工程を含む。これには、例えば、前記長さk-xの配列部分をソートする工程(例えば、前記長さk-xの配列部分を文字列としてソートする工程)が含まれ得る。
部分配列中の長さk-xの配列部分が共通する配列ごとの、長さxの部分についての出現頻度の分布は、長さkの部分配列の出現頻度から算出できる。配列の文字の種数をv(配列において各変数が取り得る値の種類)とした場合、長さkの部分配列において、k-xの配列部分が共通する配列が、k-xの配列部分が共通する配列ごとにv^x種生じる。例えば、配列が塩基配列であり、x=1とした場合、長さkの配列のセットの中には、k-1の配列部分が共通する配列ごとに、xに対応する部分がA、C、G、Tである4種の配列が存在している。長さk-xの配列部分が共通する配列ごとの、長さxの部分についての出現頻度は、それぞれに対応する長さkの部分配列の頻度データに対応する。
1つの実施形態では、前記長さxの部分が、前記部分配列の末端に存在する。例えば、前記対象配列データおよび前記コントロール配列データが塩基配列データである場合、前記長さxの部分は、前記部分配列の3’末端または5’末端であり得る。長さxの部分を部分配列の末端にとることは、比較処理の高速化・簡便化にとって望ましい。このように、k長の配列のうち末端部で異なる部分(x長の部分に該当)の頻度を検出する場合、対象となる配列の多型部位の「エッジ」(置換の場合はその位置そのものに該当するが、挿入・欠失変異の場合はその縁(エッジ)に該当する)を検出しているというように理解することができるため、本明細書において、「Polymorphic Edge Detection(PED)」または「edge detection」と呼称される場合がある。
(比較・多型の検出)
出現頻度の分布の差異の比較により、例えば、以下のような多型の検出が可能である。
多型が存在している部位がゲノム上でユニークな配列部位であり、例えば置換が存在している場合、その置換を含むk-mer部分配列は対象配列データのサブセットには存在するが、コントロール配列データのサブセットには存在しないと考えられる。当該置換変異を含むk-merではすべてこの結果が得られると考えられ、結果的に当該置換変異部位を含む2×k-1個のk-mer部分配列で出現頻度の差異が観察され、当該多型を検出することができる。そのような解析は、一例としては、図6に示されるような工程を採用して行うことが可能である。
1つの実施形態では、前記コントロール配列データのサブセットと前記対象配列データのサブセットとの間で、前記長さxの部分の配列の出現頻度が異なる場合、該長さxの部分の配列を、対象配列データにおけるコントロール配列データに対する多型として検出する。1つの実施形態では、前記コントロール配列データのサブセットと前記対象配列データのサブセットとの間で、前記長さxの部分の配列で最も高頻度のものが異なっている長さk-xの配列部分が存在する場合、該長さxの部分の配列を、対象配列データにおける多型として検出する。
1つの実施形態では、前記対象配列データおよび前記コントロール配列データが塩基配列データであり、前記対象配列データのサブセットにおける前記長さxの部分の配列で、前記コントロール配列データのサブセットにおける最も高頻度のものと同一の長さxの部分の配列がノイズ以下のカウントしか存在しない長さk-xの配列部分が存在する場合、該長さxの部分の配列を、対象配列データにおけるホモ多型として検出する。1つの実施形態では、前記対象配列データおよび前記コントロール配列データが塩基配列データであり、対象サブセットにおける前記長さxの部分の配列で、コントロール配列データのサブセットにおける最も高頻度のものと同一の長さxの部分の配列が存在し、かつ、コントロール配列データのサブセットにおける最も高頻度のものと異なる長さxの部分の配列が存在する長さk-xの配列部分が存在する場合、該長さxの部分の配列を、対象配列データにおけるヘテロ多型として検出する。k長配列のセットの頻度データの比較は、一例としては、図4に示されるような工程によって行うことが可能である。
このような比較は、1つの例としては、k-1長配列と最後の塩基の頻度のファイル、controlとtargetを処理する場合、
Figure 0007166638000008

のコマンドでk-1配列とcontrolとtargetを1行にまとめた頻度を出力することによって行うことができる。この出力結果の各行を調べる条件としては、限定されるものではないが、control、targetの両方でカウントが1塩基以下の塩基が2個以上存在し、controlあるいはtargetで10以上のカウントを示した塩基に対応するtargetあるいはcontrolの塩基のカウントが1以下である事例が1ないし2回ある場合、多型の境界を検出したとすることができる。
すなわち、コントロール配列および/または対象配列の部分配列のサブセットにおける第1のカウントが第1の閾値を超えており、かつ、コントロール配列および/または対象配列の部分配列のサブセットにおける対応する第2のカウント(すなわち、第1のカウントがコントロール配列のものである場合、第2のカウントは対象配列のものであり、第1のカウントが対象配列のものである場合、第2のカウントはコントロール配列のものである)が第2の閾値を下回る場合、多型の境界を検出したとすることができる。
第1の閾値は、配列データのカバレッジによって変動するが、例えば、10~50の範囲で設定することができる。第1の閾値は、例えば、10~40、10~30、10~20、または10~15の範囲で設定することができる。配列データのカバレッジが大きくなる場合には、第1の閾値も大きく設定することができ、例えば、ヒトゲノム解析では200倍読み程度のデータも存在するが、この場合は、第1の閾値として200を使用することができる。カバレッジを考慮して実際にその配列が存在しているといえるカウントを第1の閾値としてよく、例えば、カバレッジの約100%、約90%、約80%、約70%、約60%、約40%、約30%、または約20%等の値を用いることができる。
第2の閾値も同様に配列データのカバレッジによって変動するが、1~7の範囲で設定することができる。第2の閾値は、例えば、1、2、3、4、5、6または7であり得る。
(配列)
本発明の対象配列、コントロール配列および/またはリファレンス配列としては、多型が生じ得る任意の配列を用いることができる。なお、コントロール配列として、リファレンス配列を用いることが可能である。代表的な実施形態では、対象配列、コントロール配列および/またはリファレンス配列は、生物学的配列であり、例えば、塩基配列(DNA、RNA、およびそれらのアナログ等の配列が包含される)、アミノ酸配列、または糖鎖配列等である。生物学的配列の例としては、例えば、ゲノム配列、染色体配列、遺伝子配列、プラスミド配列、エクソン配列、タンパク質配列等が挙げられる。
対象配列データおよびコントロール配列データは、限定されるものではないが、多型を検出する上では、一定の共通性を持つ配列についての配列データであることが望ましい。しかしながら、配列の取得方法については各々同一でも異なっていてもよく、シーケンシングによって得られたデータ間での比較を行うことも、データベース等から得られたデータ間での比較を行うことも、シーケンシングによって得られたデータとデータベース等から得られたデータとの間での比較を行うことも可能である。
1つの実施形態では、対象配列データが、個体から得られた配列データであり、コントロール配列データが、該個体と同種の別の個体、またはデータベースから得られた配列データである。1つの実施形態では、対象配列データが、個体の組織試料から得られた配列データであり、コントロール配列データが、該個体の別の組織、またはデータベースから得られた配列データである。1つの実施形態では、対象配列データが、細胞試料から得られた配列データであり、コントロール配列データが、別の細胞、またはデータベースから得られた配列データである。
本発明の方法は、全長配列の情報を必要としないため、データベース等において例えば全長配列が公知でない場合にも用いることができ、例えば、対象配列データおよびコントロール配列データは、生物のゲノムに由来する塩基配列データである場合、前記ゲノムの配列が不明であってもよい。
例えば、次世代シーケンサーのリードデータで対象(ターゲット)とコントロール間の配列の直接比較による多型検出は既存の技術では不可能であった。リファレンス配列がある場合は、多型をゲノム上にマップすることができるが、リファレンスが存在しない場合(例えば、リファレンスゲノム配列が作成されていない生物)で、リファレンス情報を用いることなく対象とコントロール間での多型を検出できるというのは画期的である。リファレンスのない生物でのF2分離集団で表現型の分離と連鎖する多型を検出すれば、ゲノム上の位置が不明であったとしても、表現型に対応するDNAマーカーを得ることができ、応用範囲は非常に広いと考えられる。実際、育種にはゲノム上の位置情報も重要ではあるが、位置情報が不明でも優良形質にリンクしているDNAマーカーでの選抜ができれば、選抜育種への利用が十分に可能である。
同一個体(例えば、ヒト)の正常組織と変異を含み得る組織(例えば、癌組織)の配列を直接比較して多型を検出することも可能で、いったんリファレンスゲノムにマップしてから差を調べる方法に比べて、多型の捕捉率、捕捉精度ともに非常に高くなる。
1つの実施形態では、本発明の方法で用いる対象配列データおよび/またはコントロール配列データは、シーケンシングによって得られた塩基配列データである。シーケンシングの手法としては、サンガー法、マクサム・ギルバード法、単一分子リアルタイムシーケンシング(例えば、Pacific Biosciences、Menlo Park、California)、イオン半導体シーケンシング(例えば、Ion Torrent、South San Francisco、California)、パイロシーケンシング(例えば、454、Branford、Connecticut)、ライゲーションによるシーケンシング(例えば、Life Technologies、Carlsbad、CaliforniaのSOLiDシーケンシング)、合成および可逆性ターミネーターによるシーケンシング(例えば、Illumina、San Diego、California)、透過型電子顕微鏡法などの核酸イメージング技術、ナノポアシーケンシングなどがある。
1つの実施形態では、本発明の方法で用いる対象配列データおよび/またはコントロール配列データは、次世代シーケンシングによって得られた配列データであり得る。次世代シーケンシングとしては、シーケンシングバイシンセシス、パイロシーケンシング、ライゲーションによるシーケンシング、イオン半導体シーケンシング、ナノポアシーケンシング等が挙げられる。次世代シーケンシングデータからの多型の検出においては、リファレンスへのマッピングやアセンブリによって精度が制限されていたため、本発明の方法を用いた場合に大きな利益が得られると考えられる。
1つの実施形態では、本発明の方法で用いる対象配列データおよび/またはコントロール配列データは、ジニトロフェニル化法、ヒドラジン分解法、カルボキシペプチダーゼ法、エドマン法もしくはそれらの方法を自動化する装置(ペプチドシーケンサーあるいはプロテインシーケンサー)を用いる方法、質量分析(例えば、タンデム質量分析計(MS/MS))を用いた方法(例えば、シーケンスタグ法)等から得られたアミノ酸配列データである。
本発明の対象配列データおよび/またはコントロール配列データの由来となる生物種としては、生物学的配列を有するものである以上は何ら制限されない。一部を例示すると、動物としては、ヒトもしくは非ヒト哺乳動物(例えば、マウス、ラット、ウサギ、ヒツジ、ブタ、ウシ、ウマ、ネコ、イヌ、サル、チンパンジー)、鳥類、爬虫類、両生類、魚類等の脊椎動物、無脊椎動物、例えば、昆虫、線形動物などを挙げることができる。植物としては、イネ、コムギ、トウモロコシ、ジャガイモ、オオムギ、サツマイモ、ソバ、シロイヌナズナ、ミヤコグサ、トマト、キュウリ、キャベツ、白菜、ナス、サトウキビ、ソルガム、リンゴ、ミカン、バナナ、桃、ポプラ、松、杉、被子植物、裸子植物、シダ、コケ、藻類などを挙げることができる。その他、真菌、細菌、ウイルス等でもよい。
さらに、これらの生物の一部分、例えば、組織、細胞等に由来する対象配列データおよび/またはコントロール配列データを解析し、多型を検出することも可能である。
(変異)
本発明の方法は、例えば、置換、挿入、欠失、コピー数変異、STRP(short tandem repeat polymorphism)、逆位または転座等の多型の検出に用いることができる。変異のエッジの部分が検出されるため、挿入・欠失の結果、長さxの配列に違いがあればそのエッジ部分を検出することができる。k-mer内に収まりきる場合であれば、STRP(short tandem repeat polymorphism)を検出することも可能である。STR(short tandem repeat)は、マイクロサテライトとも称され、2~7塩基からなる配列が2~数十回反復するもので、この回数に多型が見られる。部分配列の出現頻度によって、コピー数多型(CNV)を検出することもできる。エッジ検出という観点からは、逆位、転座のエッジも検出することが可能である。
特に、多型が置換である場合には、本発明の方法は、非常に高い検出力を発揮することが可能である。
(位置の特定)
対象配列に対するリファレンス配列が存在する場合、本発明の方法は、対象配列に対するリファレンス配列における前記多型の位置を特定する工程をさらに含むことができる。例えば、対象配列データおよびコントロール配列データが、生物のゲノムに由来する塩基配列データである場合、多型のゲノム上の位置を特定する工程をさらに含むことができる。この位置の特定は、本発明の方法が、多型を周囲の配列と関連づけて検出する(例えば、x長部分の多型がk-x長の配列と関連付けられる)ことを可能にしているため、リファレンス配列に対して検索を行うことにより、簡便に行うことが可能である。
リファレンス配列に対する検索は、一例としては、図5に示されるような工程によって二分検索用リファレンスゲノム配列データを作成し、その後二分検索によって多型境界塩基のマッピングによって検索を行うことができる。
他の方法として、unixのjoinコマンドによって対象配列データ中の部分配列の、リファレンス配列における位置および向きを出力することによってマッピングを行う方法を用いることができる。より詳細には、対象配列のコントロール配列上の位置を決定する方法であって、a)コントロール配列中の複数のk長の部分配列について、配列およびコントロール配列中の位置と向きを出力する工程と、b)対象配列中の複数のk長の部分配列について、配列および対象配列中の位置を出力する工程と、c)a)およびb)で得られた配列を比較し、同一の部分配列に対応するコントロール配列中の位置と対象配列中の位置とを対応付ける工程とを含み、ここで、kは、対象配列の長さを超えない長さである、方法を使用することができる。当該方法については、本出願人により本出願と同日に出願された「挿入・欠失・逆位・転座・置換検出法」との名称の出願(整理番号NG012PCT/F5-18PCT075)を参照することができる。
(確認(verify))
本発明の方法は、検出した多型について確認する工程をさらに含むことができる。確認は、例えば、検出された多型の部位について、リファレンス配列またはコントロール配列から作成したクエリ配列セットを用いて、対象配列データおよび/またはコントロール配列データとの比較を行うことによって行うことができる。クエリ配列セットは、リファレンス配列またはコントロール配列において前記多型に該当する部位の文字を異なる文字に置換した変異型クエリ配列セット、および/またはリファレンス配列またはコントロール配列において前記多型に該当する部位の文字を置換していない野生型クエリ配列セットを含み得る。
本発明の方法は、対象配列データおよびコントロール配列データが塩基配列データである場合、検出された多型の部位について、対象配列データおよび/またはコントロール配列データの相補鎖の配列データに対して、リファレンス配列またはコントロール配列から作成したクエリ配列セットとの比較を行い確認する工程をさらに含むことができる。本発明の方法は、対象配列データおよびコントロール配列データが塩基配列データである場合、検出された多型の部位について、対象配列データおよび/またはコントロール配列データの対立遺伝子の配列データに対して、リファレンス配列またはコントロール配列から作成したクエリ配列セットとの比較を行い確認する工程をさらに含むことができる。確認する工程は、一例としては、図6に示されるフローに従って、図6に示される工程を適宜採用して行うことができる。ここで、対立遺伝子の配列データとして、実際の遺伝子の存在の有無とは関係なく、野生型に対する変異型の塩基を有する配列データを用いることができる。
以下、コントロール配列が野生型(つまりリファレンスゲノム配列とほぼ同一である)から得られた場合、または、コントロール配列が、リファレンスゲノム配列から対象配列と同一長で作成された配列である場合について例示する。
次世代シーケンサーで読まれたショートリードの塩基配列長がLの場合、リファレンスゲノム配列の対象となる多型塩基位置を起点として、L-1塩基前の位置からL-1塩基後の位置までの2L-1塩基長の配列を得て、多型塩基位置を推定された多型塩基に置換した置換配列とコントロールの非置換配列を作成する。置換配列、および、非置換配列それぞれを1塩基ずつずらしながらL長のクエリ配列セットを作成する(例えば、図14に例示される)。個々のクエリ配列を用いて対象配列とコントロール配列に対して完全一致するカウントを取得する。ホモ型変異の場合は変異を導入したクエリ配列では対象配列のカウントが大部分になり、ヘテロ型変異の場合は、対象配列とコントロール配列への一致カウントが概ね半分ずつとなると考えられる。非置換配列をクエリにした場合は、基本的には大部分コントロール配列にヒットする。非置換配列が対象配列にヒットする場合は、多型ではないと判断し、除外することが可能である。
コントロール配列、対象配列の検索のため、コントロール配列、対象配列、それぞれ相補鎖の配列も合わせたのち辞書順にソートして同一配列は一つにまとめたデータセットを用いて、二分探索法によりクエリ配列を検索する。実施形態において、Fastqファイルからの塩基配列データをソートしてユニークなデータセットにする時点で、つまり最初の段階で、各リードの相補鎖も一緒にソートしてユニーク処理することも可能である。
k-merを切り出す前にsortとuniqの処理をすることは、シーケンス反応時にPCR増幅のステップが入る場合があり(入らないキットも存在する)、同一の配列が複数回リードデータの中に出てくる場合があることに対処する上で好ましい場合がある。そのまま処理するとk-merの分布が歪む可能性がある。同一配列ではあるがNが含まれるようなリードも別物と認識されて歪む原因になり得るため、Nを含まない配列で、その配列とその相補鎖配列をsortしてuniq処理した配列からk-merを得ることが可能である。
歪むことを許容する場合、もとの配列の長さは揃っている必要はなく、サンガー法で得られたような長さが一定でない配列データでもk-merにして多型検出、マッピングまで行うことが可能である。
確認工程における、sort_uniq配列を二分検索して、変異型と野生型のリード数を調べる工程では、配列データの長さが揃っていることが好ましい。本発明者らの知見によれば、確認工程の前にk-1配列で多型をマップした段階で、順鎖と相補鎖の両方で多型が検出された場合は、ほぼ間違いなくSNPであることがわかっている。順鎖、相補鎖のどちらかがリピート領域にあって片側しか検出できないものでも実際にSNPである場合があるが、このような場合、sort_uniq配列を二分検索して確認することで、実際のSNPかどうかの判断をすることができる場合があり得る。二分検索による確認を行って捕捉率を上げるためには、スタートのショートリード(次世代シーケンサーで得られた塩基配列データ)の長さが揃っていることが好ましい場合がある。対象とコントロールの配列長が同一である必要はなく、対象配列データとコントロール配列データそれぞれの中で長さが一定であれば、好適に二分検索によって確認を行うことができる。
あるいは、以下のように確認を行うことができる。リファレンス配列から対象配列のL長で変異部位を含むようにして部分配列を切り出し、変異に置換したセットと、置換しないセットを作り、位置関係、変異の有無等の記載と一緒にソートして出力する。このデータとソートした対象配列を、unixコマンドのjoinで処理(または適切な等価な処理)して、対象配列中に含まれる、野生型と変異型の配列を選び、配列数を変異部位ごとに調べる。選んだ配列を、sortした後、uniq-cのコマンドで配列数を数えることができる。同様の操作を、コントロール配列(L’長)に対しても行う。対象個体に対するコントロール個体がある場合は、この個体から得られたリード配列をコントロール配列として用いることができる。コントロール個体がない場合は、リファレンス配列から、L長で切り出して作った配列をコントロール配列として使うことが可能である。対象配列およびコントロール配列で長さが違う場合は、それぞれの長さに対応する変異型と野生型のデータセットを作り、対応する個数を調べることができる。当該方法については、本出願人により本出願と同日に出願された「挿入・欠失・逆位・転座・置換検出法」との名称の出願(整理番号NG012PCT/F5-18PCT075)を参照することができる。
例えば、イルミナ社の次世代シーケンサー(例えば、HiSeq)は同じ長さの配列を出力するため、対象配列データにおいて、特に長さをそろえる処理をする必要はない。この場合、対象配列データと同じ長さを有するクエリ配列のセットを作成することが可能であるため、直接二分検索を行うことができ、確実な結果を得る上で有利であり得る。
長さにばらつきがあるショートリードによる配列データあるいはショートリードの集合ではない配列データ(例えば、サンガー法で得られた配列)でも、本発明において適用可能である。確認工程においては、長さを揃える処理(例えば、リード中のクオリティスコアが最大となるL長配列を選択して切り出す、一端からL長配列を切り出す等)の処理をしたデータを用いて二分検索を行うことが可能である。あるいは、配列データ(リード)をBLASTのターゲット配列(データベース)にして、配列データとは長さの異なるクエリ配列を検索して数を数えることも可能である。
コンピュータで計算させる場合、sort_uniq配列や二分探索用のリファレンス配列に対して二分探索で完全一致の配列を探索する場合に、ファイルサイズが大きくなる場合が多いため、高速のシステムを利用するか、適切に高速化を図ることが好ましい。高速化の手法としては、すべてオンメモリで計算を行う、SSD等高速ディスクにファイルを置く等、ハード的な手法が存在する。ソフト的には、単純にソートされたファイルに対して二分検索を行うのではなく、検索対象のファイルをBurrows-Wheeler変換してより高速化を行うこともできる。
(プログラム、記録媒体およびシステム)
1つの局面において、本発明は、本発明の多型を検出する方法をコンピュータに実施させるための方法を実装するプログラム、該プログラムを記録した記録媒体、およびこれを実現するためのシステムを提供する。ここで採用され得る任意の特徴は本明細書の多型を検出する方法の説明に記載される任意の特徴またはその組み合わせを採用することができる。
したがって、1つの実施形態において、対象配列データにおいてコントロール配列データに対する多型を検出する方法をコンピュータに実行させるためのプログラムであって、該方法は、
a)該対象配列データの長さkの部分配列のサブセットをコンピュータに保存する工程であって、ここで、kは対象配列およびコントロール配列の全長以下の長さである、工程と、
b)該対象配列データの長さkのサブセットの各々の部分配列の出現頻度を算出する工程と、
c)該コントロール配列データの長さkの部分配列のサブセットにおける各々の部分配列の出現頻度をコンピュータに保存する工程と、
d)対象配列とコントロール配列とを比較し、該出現頻度の分布の差異に基づいて、多型を検出する工程と
を包含する、プログラムが提供される。さらなる実施形態において、方法が、前記部分配列中の多型ではない部分の少なくとも一部を含む配列(前記部分配列全体であり得る)を、検出された前記多型の名称として表示する工程をさらに含む、プログラムが提供される。プログラムはどのような言語で記述されてもよい。
別の実施形態において、対象配列データにおいてコントロール配列データに対する多型を検出する方法をコンピュータに実行させるためのプログラムを格納する記録媒体であって、該方法は、
a)該対象配列データの長さkの部分配列のサブセットをコンピュータに保存する工程であって、ここで、kは対象配列およびコントロール配列の全長以下の長さである、工程と、
b)該対象配列データの長さkのサブセットの各々の部分配列の出現頻度を算出する工程と、
c)該コントロール配列データの長さkの部分配列のサブセットにおける各々の部分配列の出現頻度をコンピュータに保存する工程と、
d)対象配列とコントロール配列とを比較し、該出現頻度の分布の差異に基づいて、多型を検出する工程と
を包含する、記録媒体が提供される。さらなる実施形態において、方法が、前記部分配列中の多型ではない部分の少なくとも一部を含む配列(前記部分配列全体であり得る)を、検出された前記多型の名称として表示する工程をさらに含む、記録媒体が提供される。プログラムはどのような言語で記述されてもよい。1つの実施形態では、記録媒体は、内部に格納され得るROMやHDD、磁気ディスク、USBメモリ等のフラッシュメモリなどの外部記憶装置でありうる。
別の実施形態において、対象配列データにおいてコントロール配列データに対する多型を検出するためのシステムであって、該システムは、該対象配列データおよび該コントロール配列データの長さkの部分配列のサブセットの各々の部分配列の出現頻度を提供するように構成された配列データ処理部であって、ここで、kは対象配列およびコントロール配列の全長以下の長さである、配列データ処理部と、対象配列とコントロール配列とを比較し、該出現頻度の分布の差異に基づいて、多型を検出する工程とを行うように構成された、配列データ計算部とを備える、システムが提供される。さらなる実施形態において、前記部分配列中の多型ではない部分の少なくとも一部を含む配列(前記部分配列全体であり得る。)を、検出された前記多型の名称として表示する表示手段をさらに含む、システムが提供される。
次に、図15Aの機能ブロック図を参照して、本発明のシステム1の構成を説明する。なお、本図においては、単一のシステムで実現した場合を示しているが、複数のシステムで実現される場合も本発明の範囲に包含されることが理解される。
本発明のシステム1000は、コンピュータシステムに内蔵されたCPU1001にシステムバス1020を介してRAM1003、ROMやHDD、磁気ディスク、USBメモリ等のフラッシュメモリなどの外部記憶装置1005及び入出力インターフェース(I/F)1025が接続されて構成される。入出力I/F1025には、キーボードやマウスなどの入力装置1009、ディスプレイなどの出力装置1007、及びモデムなどの通信デバイス1011がそれぞれ接続されている。外部記憶装置1005は、情報データベース格納部1030とプログラム格納部1040とを備えている。何れも、外部記憶装置1005内に確保された一定の記憶領域である。
このようなハードウェア構成において、入力装置1009を介して各種の指令(コマンド)が入力されることで、又は通信I/Fや通信デバイス1011等を介してコマンドを受信することで、この記憶装置1005にインストールされたソフトウェアプログラムがCPU1001によってRAM1003上に呼び出されて展開され実行されることで、OS(オペレーションシステム)と協働して本発明の対象配列データにおいてコントロール配列データに対する多型を検出する方法の機能を奏するようになっている。もちろん、このような協働する場合以外の仕組みでも本発明を実装することは可能である。
本発明の実装において、対象配列データの長さkの部分配列のサブセットをコンピュータに保存する工程であって、ここで、kは対象配列およびコントロール配列の全長以下の長さである、工程を行う際に、対象配列データおよび/または対象配列データの長さkの部分配列のデータは、入力装置1009を介して入力され、あるいは、通信I/Fや通信デバイス1011等を介して入力されるか、あるいは、データベース格納部1030に格納されたものであってもよい。次に、対象配列データの長さkのサブセットの各々の部分配列の出現頻度を算出する工程は、プログラム格納部1040に格納されたプログラム、または、入力装置1009を介して各種の指令(コマンド)が入力されることで、又は通信I/Fや通信デバイス1011等を介してコマンドを受信することで、この外部記憶装置1005にインストールされたソフトウェアプログラムによって実行することができる。あるいは、あらかじめ算出されている出現頻度を、入力装置1009を介して入力することができる。出現頻度データは、出力装置1007を通じて出力されるかまたは情報データベース格納部1030等の外部記憶装置1005に格納されてもよい。
次に、コントロール配列データの長さkの部分配列のサブセットにおける各々の部分配列の出現頻度をコンピュータに保存する工程を行う差異、コントロール配列データ、コントロール配列データの長さkの部分配列のサブセットのデータ、または部分配列の出現頻度のデータは、入力装置1009を介して入力され、あるいは、通信I/Fや通信デバイス1011等を介して入力されるか、あるいは、データベース格納部1030に格納されたものであってもよく、プログラム格納部1040に格納されたプログラム、または、入力装置1009を介して各種の指令(コマンド)が入力されることで、又は通信I/Fや通信デバイス1011等を介してコマンドを受信することで、この外部記憶装置1005にインストールされたソフトウェアプログラムによってこれらのデータを処理してコントロール配列データの長さkの部分配列のサブセットにおける各々の部分配列の出現頻度を提供してもよい。
対象配列とコントロール配列とを比較し、該出現頻度の分布の差異に基づいて、多型を検出する工程は、プログラム格納部1040に格納されたプログラム、または、入力装置1009を介して各種の指令(コマンド)が入力されることで、又は通信I/Fや通信デバイス1011等を介してコマンドを受信することで、この外部記憶装置1005にインストールされたソフトウェアプログラムによって実行することができる。
データベース格納部1030には、これらのデータや計算結果、もしくは通信デバイス1011等を介して取得した情報が随時書き込まれ、更新される。各入力配列セット中の各々の配列、参照データベースの各遺伝子情報ID等の情報を各マスタテーブルで管理することにより、蓄積対象となるサンプルに帰属する情報を、各マスタテーブルにおいて定義されたIDにより管理することが可能となる。
データベース格納部1030には、上記計算結果は、配列に関する情報、例えば、生物学的情報、生化学的情報、医学的情報、例えば疾患、障害、生体情報等の既知の情報と関連付けて格納されてもよい。このような関連付けは、ネットワーク(インターネット、イントラネット等)を通じて入手可能なデータをそのまままたはネットワークのリンクとしてなされてもよい。
また、プログラム格納部1040に格納されるコンピュータプログラムは、コンピュータを、上記した処理システム、例えば、配列データの提供、部分配列サブセットの提供、出現頻度データの算出、出現頻度データの比較、多型の検出、多型の確認などの処理を実施するシステムとして構成するものである。これらの各機能は、それぞれが独立したコンピュータプログラムやそのモジュール、ルーチンなどであり、上記CPU1001によって実行されることでコンピュータを各システムや装置として構成させるものである。なお、本発明の例示においては、それぞれのシステムにおける各機能が協働してそれぞれのシステムを構成しているものとするが、この処理のためのプログラムもまた、それぞれ外部記憶装置または通信デバイスまたは入力装置を介して提供されうる。
本発明がシステムとして構成される場合は、対象配列データおよび/またはコントロール配列データ、それらの長さkのサブセットのデータ、ならびに/あるいはそれらの出現頻度データの提供は、まとめて配列データ処理部としてもよい。また、出現頻度の分布の比較および多型の検出は、配列データ計算部としてまとめてもよい。
また、図15Bに示されるように、クラスター構造を有する計算システムによって本発明の方法を実装してもよい。1つの実施形態では、システムはクラスター構成であり、ヘッドとノードからなる。ノードは検索の高速化を図るため、主記憶装置にSSDを用いることができる。1つの実施形態では、ヘッド1台に対して複数のノード(例えば12台)で運用することができる。1つの実施形態では、計算システムはクラスター構造を持ち、主コンピュータ(クラスターヘッド)に大容量記憶装置(HDD)を搭載して解析データおよび結果を保存する。クラスターヘッドより、分割したデータを各ノードに送り計算を実行し、結果をクラスターヘッドに集約する。クラスターヘッド、ノード共に、中央制御素子(CPU)、メモリ(RAM)を搭載し、通信インターフェース(NIC)を介してデータの通信を行い得る。ノードには高速での検索処理をするため、ソリッドステートドライブ(SSD)を主記憶装置とすることができる。各ノードに搭載されるCPU、RAM、SSD等は、他のノードと共有されてもよく、物理的に分離していてもよい。
(例示的計算フロー)
本発明は、例えば、以下のフローにおいて実施することができる。
1.次世代シーケンサーから得られたFastq形式のファイルより、配列内にNを含まない塩基配列を選んで、相補鎖の配列と合わせて、それぞれ出力する。
2.出力されたファイル(reads)を辞書順にソートし同一の配列は一つにまとめる。
unixコマンドは以下の通り。
Figure 0007166638000009
相同な配列がゲノム上に複数存在する部分に対しての解析を行う場合は、readsに対してsortのみを行い、uniq処理は行わない場合もあり得る。例えば、マルチコピー領域に対してのCNVを検出する場合にはsortのみを行うことが有利であり得る。
3.ソートされた各塩基配列をそれぞれ5’末端から1塩基ずつずらしながら、k-mer(ここでは20塩基)の配列を対象の塩基配列の3’末端に到達するまで繰り返し出力(k-mer_file)する。
4.出力されたk-merの配列を辞書順にソートし、同一配列は一つにまとめて、出現回数を配列と共に表記したファイルを作成する。
unixコマンドは以下の通り。
Figure 0007166638000010
uniq -cコマンドなどのように、出力結果が頻度→配列の順になる場合、uniq-c等のコマンドの後に配列→頻度の順に出力されるフィルタープログラムを通してもよい。
5.k-mer_count_fileの各行のデータに対して、配列の5’末端よりk-1塩基の配列を得て、3’末端の塩基すなわちk番目の塩基をA、C、G、Tの出現回数として表記したデータに変換する。
k-1merの配列 Aの回数 Cの回数 Gの回数 Tの回数
という形式で出力される。
Figure 0007166638000011
6.このようなデータを対象(target)と比較(control)のサンプルから得られた次世代シーケンサーの配列データからそれぞれ得る。
7.controlとtargetの5の方法で作成したデータをk-1merの配列でまとめたデータを作成する。
unixコマンドは以下の通り
Figure 0007166638000012
8.joinコマンドで得られたデータでcontrolとtargetで異なる塩基でカウントが示されるデータを選び出す。
Figure 0007166638000013
この例では、controlでは、ACTTTCTTCAAGGTCTGTT(配列番号225)に続く塩基はGであるが、targetではCである。つまり、k-1merのユニークな識別子(名称)に続く塩基がG型、あるいはC型という表記で多型を表現することができる。それぞれの塩基に対応する数は、この多型が検出された独立のreadの数である。
9.このk-1merの識別子をリファレンスゲノム上にマッピングすれば、容易に多型位置を決めることができる。
10.本法は配列情報をリファレンスゲノムにマップする前にcontrolとtarget間の多型を検出するので、リファレンスゲノムが未知の生物種でも多型を検出できる。
11.ゲノム上の場所が決められない場合は、多型の名称はk-1merの配列自身で表すことができ、遺伝子型はそれに続く多型塩基となる。このデータセットを用いて、表現型に対するアソシエーション解析が可能である。
12.リファレンスゲノムへのマッピングには、リファレンスゲノムを3の方法と同様に各染色体を5’側から1塩基ずつずらしながらk-merで切り出し、k-mer配列、染色体番号、位置、向きを同一行に表記したデータをk-merの配列で辞書順に整列したデータに対して、二分検索法によりk-1merの配列のゲノム上の位置を決めることができる。
Figure 0007166638000014
1,2行目のように同一配列も複数行でそれぞれの位置がわかるので、対象配列がリピート領域に存在しても、検索では対応する候補領域が列挙されて出力される。
今回は、複数の位置情報が出力された場合は、位置不定として多型情報の出力から除外している。複数の位置のどれかという出力の仕方も可能である。
13.リファレンスゲノムにマップされた多型の確認を行う。2で作成されたsort_uniq配列(ターゲット)と同様に、リファレンス配列も5’末端から2塩基ごとに2で作成された配列と同じ長さの配列を切り出し、相補鎖とともにsort,uniq処理をしたデータを作成する。このリファレンスとターゲットのsort_uniqデータが検索の対象となる。
次に、リファレンスゲノム配列から、多型位置を含む2で作成された配列と同じ長さの配列セット(検索を行うクエリセット)を切り出す。この配列セットは配列の5’末端から3’末端までのすべての位置で多型位置の塩基を含む配列セットである。配列セット数は配列長と同一となる。リファレンスの配列セット(リファレンスセット)に対して、多型位置の塩基を予想された多型塩基に置換した配列セット(ターゲットセット)を作成する。クエリ用配列セットの作成は、図14に例示される。
14.リファレンスセット、ターゲットセットをクエリにして、リファレンスゲノム、及び、2で作成されたsort_uniq配列をそれぞれ、検索(例えば、二分検索法またはjoinコマンドを用いた方法)し、それぞれのセットに対してsort_uniq配列が何個マッチするかを調べる。リファレンスセットではリファレンスのsort_uniq配列のみにヒットする。これに対して、ターゲットセットでは、ホモ型の変異の場合は、ターゲットのsort_uniqデータのみから検出され、ヘテロ型の変異の場合は、リファレンスと、ターゲットのsort_uniqデータの両方から検出されるはずである。このようにして、予想と一致する検索値を示した多型を抽出すると、以下のように結果を出力することが可能である。
Figure 0007166638000015
(組み合わせ)
本明細書において、上述のとおり置換、コピー数多型、STRP、挿入、欠失、逆位または転座を検出するのに有用な方法を記載しているが、かかるプロセスは、置換、挿入、欠失、逆位または転座を検出するのに有用な以下に記載するプロセスと組み合わせて行うことができる。例えば、組み合わせた方法は、図18に示されるようなフローに従って実行することが可能である。
このようにプロセスを組み合わせることで、配列に存在し得る多くの種類の多型を網羅的に高い検出力で検出することが可能である。このようなプロセスの組み合わせは、例えば、複数のプロセスを同時に、並行して、または逐次的に行うことによって達成することが可能である。例えば、図15Bに示されるようなクラスター構造を有する計算システムによって、異なるノードを用いてそれぞれのプロセスを行うことにより、組み合わせの方法を実装することができる。
置換、挿入、欠失、逆位または転座を検出するのに有用なプロセスとしては、対象配列データの配列中の少なくとも2ヶ所の部分配列の、コントロール配列上の位置を特定する工程を含むプロセスがある。ここで、部分配列は、k長の部分配列を用いることができる。好ましくは、コントロール配列は、配列上の位置情報が特定できる配列であり、より好ましくは、コントロール配列はリファレンス配列である。
プロセスは、対象配列データにおける部分配列間の位置関係と、コントロール配列上の部分配列間の位置関係とを比較する工程を含み得る。ここで、対象配列データにおける部分配列間の位置関係と、コントロール配列上の部分配列間の位置関係とが異なっている場合、目的とする多型があると判定することができる。例えば、部分配列が、コントロール配列の異なる配列構造体上に存在する場合、転座が生じていると判定すること、部分配列が、コントロール配列の同一の配列構造体上に存在し、かつ、向きが対象配列データ上のものと異なっている場合、逆位が存在すると判定すること、部分配列が、コントロール配列の同一の配列構造体上に存在し、向きが対象配列データ上のものと同一であり、部分配列の距離が、コントロール配列上で対象配列データ上の距離より短い場合、欠失が存在すると判定すること、および/または部分配列が、コントロール配列の同一の配列構造体上に存在し、向きが対象配列データ上のものと同一であり、部分配列の距離が、コントロール配列上で対象配列データ上の距離より長い場合、挿入が存在すると判定することを含み得る。位置関係が異ならない場合、処理を終了してもよく、目的とする多型はないと判定してもよく、対象配列データにおける部分配列部位間の文字を、対応するコントロール配列上の文字と比較して不一致となる部位を検出する工程をさらに行い、不一致となる部位が存在する場合、置換が存在すると判定してもよい。
プロセスは、位置関係が異なっている場合、目的とする多型があると判定し、対象配列データにおける部分配列部位間の文字を、対応するコントロール配列上の文字と、部分配列部位を始点として順次比較して不一致となる部位を検出する工程を含み得る。かかる工程により、検出した多型の境界塩基を検出することができる。
例えば、このようなプロセスを組み合わせた場合、本発明の1つの実施形態では、
対象配列データにおいてコントロール配列データに対する多型を検出する方法であって、
(1)a)該対象配列データの長さkの部分配列のサブセットの各々の部分配列の出現頻度を提供する工程と、
b)該コントロール配列データの長さkの部分配列のサブセットの各々の部分配列の出現頻度を提供する工程と、
c)対象配列とコントロール配列とを比較し、該出現頻度の分布の比較に基づいて、多型を検出する工程と
によって、対象配列データにおける置換、コピー数多型、STRP、挿入、欠失、逆位または転座を検出するプロセスと、
(2)a)該対象配列データの配列中の少なくとも2ヶ所の部分配列の、該コントロール配列上の位置を特定する工程と、
b)対象配列データにおける該部分配列間の位置関係と、コントロール配列上の該部分配列間の位置関係とを比較する工程と、
c)対象配列データにおける該部分配列間の位置関係と、コントロール配列上の該部分配列間の位置関係とが異なっている場合、目的とする多型があると判定し、該対象配列データにおける該部分配列部位間の文字を、対応するコントロール配列上の文字と、該部分配列部位を始点として順次比較して不一致となる部位を検出する工程と
によって、対象配列データにおける挿入、欠失、逆位、転座または置換を検出するプロセスと
を包含する、方法が提供される。
本発明のさらなる実施形態では、例えば、
対象配列データにおいてリファレンス配列データに対する多型を検出する方法であって、リファレンス配列データから、各々の位置情報と関連付けられたリファレンス配列のk長の部分配列セットを作成する工程を含み、さらに、
(A1)該対象配列データの長さkの部分配列のサブセットを生成し、ユニークな長さkの部分配列の出現頻度を提供する工程と、
(A2)該リファレンス配列のk長の部分配列セットの、ユニークな長さkの部分配列の出現頻度を提供する工程と、
(A3)該対象配列と該リファレンス配列とを比較し、該出現頻度の分布の比較に基づいて、挿入、欠失、置換、コピー数多型、STRP、逆位または転座を検出する工程とを包含するプロセスと
(B1)該対象配列データの配列中の少なくとも2ヶ所のk長の部分配列をクエリとして、該リファレンス配列のk長の部分配列セットに対して二分検索を行い、該少なくとも2ヶ所の部分配列の、リファレンス配列上の位置を特定する工程と、
(B2)該対象配列データにおける該部分配列間の位置関係と、該リファレンス配列上の該部分配列間の位置関係とを比較する工程と、
(B3)該対象配列データにおける該部分配列間の位置関係と、該リファレンス配列上の該部分配列間の位置関係とが異なっている場合、挿入、欠失、逆位または転座があると判定し、該対象配列データにおける該部分配列部位間の文字を、対応するコントロール配列上の文字と、該部分配列部位を始点として順次比較して不一致となる部位を検出する工程を包含し、必要に応じて、
(B4)該位置関係が異ならない場合に、該対象配列データにおける該部分配列部位間の文字を、対応する前記コントロール配列上の文字と比較して不一致となる部位を検出する工程をさらに含み、不一致となる部位が存在する場合、置換が存在すると判定する工程をさらに含む、プロセスと、
を、同時に、並行して、または逐次的に行うことを特徴とする、方法が提供される。
(一般技術)
本明細書において用いられる分子生物学的手法、生化学的手法、微生物学的手法、バイオインフォマティクスは、当該分野において公知であり、周知でありまたは慣用される任意のものが使用され得る。
本明細書において「または」は、文章中に列挙されている事項の「少なくとも1つ以上」を採用できるときに使用される。「もしくは」も同様である。本明細書において「2つの値」の「範囲内」と明記した場合、その範囲には2つの値自体も含む。
本明細書において引用された、科学文献、特許、特許出願などの参考文献は、その全体が、各々具体的に記載されたのと同じ程度に本明細書において参考として援用される。
以上、本発明の理解を容易にするために好ましい実施形態を示して説明してきた。以下に、実施例に基づいて本発明を説明するが、上述の説明および以下の実施例は、例示の目的のみに提供され、本発明を限定する目的で提供したのではない。従って、本発明の範囲は、本明細書に具体的に記載された実施形態にも実施例にも限定されず、特許請求の範囲によってのみ限定される。
(実施例1:イネSNP検出および検証)
(概要)
本発明の方法を用いて、以下のイネ配列データを用いてイネのSNPの検出および検証を行った。
対象:N1S5、N1S6、N1S7、N1S10
コントロール:N1
リファレンス:イネリファレンスゲノム(IRGSP1.0)
なお、本明細書において、サンプルの表記は、大文字で表記しても小文字で表記しても同じものを指すことに留意されたい。
(材料および方法)
(サンプル)
イネ品種日本晴の種子1粒をN1個体とし、発芽させて育てて葉をサンプリングした。N1個体に実った種子をN1S1とし、その種子1粒を発芽させて育てて葉をサンプリングした。N1S1個体に実った種子をN1S2とし、同様にして、N1S10世代まで、葉と種子をサンプリングした。
(次世代シーケンシング)
N1、N1S5、N1S6、N1S7、N1S10の葉のサンプルについて、次世代シーケンサーでの全ゲノム解析を行った。シーケンサーは、Illumina社のHiSeq2000を用い、ペアエンドで解析した。リード長はN1、N1S5、N1S6、N1S7は100塩基、N1S10のみ101塩基であった。
シーケンスライブラリーを用いて、シーケンスの鋳型となるクラスターを形成し、鋳型DNAの塩基配列を取得した。シーケンスデータの解析は付属のソフトウェアを使用しベースコールを行い、fastq形式ファイルとして出力した。
N1、N1S6の解析は以下のとおり製造業者のマニュアルに従った条件で行った。
表5:各作業に使用したマニュアル名、バージョン番号
Figure 0007166638000016

表6:クラスター形成、シーケンス及びシーケンス解析に使用した装置、試薬、ソフトウェア
Figure 0007166638000017
N1S5、N1S7、N1S10の解析は以下のとおり製造業者のマニュアルに従った条件で行った。
表7:各作業に使用したマニュアル名、バージョン番号
Figure 0007166638000018

表8:クラスター形成、シーケンス及びシーケンス解析に使用した装置、試薬、ソフトウェア
Figure 0007166638000019
(多型の検出)
各サンプルについて次世代シーケンサーから得られたFastq形式のファイルより、配列内にNを含まない塩基配列を選んで、相補鎖の配列と合わせて、それぞれ出力した。
出力されたファイル(reads)を辞書順にソートし同一の配列は一つにまとめた。unixコマンドは以下の通りであった。
Figure 0007166638000020
ソートされた各塩基配列をそれぞれ5’末端から1塩基ずつずらしながら、k-mer(本実施例では20塩基)の配列を対象の塩基配列の3’末端に到達するまで繰り返し出力(k-mer_file)した。
出力されたk-merの配列を辞書順にソートし、同一配列は一つにまとめて、出現回数を配列と共に表記したファイルを作成した。unixコマンドは以下の通りであった。
Figure 0007166638000021
ここで、コントロール配列:N1、対象配列:N1S7について、k-mer配列を整列させ、多型の検出を行った。多型が存在している部位が、ゲノム上でユニークな配列部位であり、そこに例えば一塩基置換が存在している場合、その置換を含むk-mer配列は対象配列では存在するが、コントロール配列では存在しないと考えられる。したがって、当該置換変異を含むk-mer配列では、対象配列では存在し、コントロール配列で存在しないため、結果的に当該置換変異部位を含むk個のk-mer配列で出現頻度の差異が見られると考えられる。k-mer配列の整列は、リファレンス配列に沿う形で整列させた。
k-mer_count_fileの各行のデータに対して、配列の5’末端よりk-1塩基の配列を得て、3’末端の塩基すなわちk番目の塩基をA、C、G、Tの出現回数として表記したデータに変換した。「k-1merの配列 Aの回数 Cの回数 Gの回数 Tの回数」という形式で出力した。
コントロール配列と対象配列について、上記工程で出力したデータをk-1merの配列でまとめたデータを作成した。unixコマンドは以下の通りであった。
Figure 0007166638000022
さらに、joinコマンドで得られたデータでcontrolとtargetで異なる塩基でカウントが示されるデータを検出した。本実施例においては、イネゲノムの40~50倍読みの次世代シーケンシングデータを対象にしたため、塩基の頻度が100を超える場合はリピート配列部位として除外した。次にcontrol、targetの両方でカウントが1塩基以下の塩基が2個以上存在し、controlあるいはtargetで10以上のカウントを示した塩基に対応するtargetあるいはcontrolの塩基のカウントが1以下である事例が1ないし2回ある場合、多型の境界を検出したものとし、当該事例が生じているk-1merの配列を、多型部位を含むものとして検出した。
リファレンスゲノムへのマッピングのために、リファレンスゲノム(イネリファレンスゲノム(IRGSP1.0))を上記工程と同様に各染色体を5’側から1塩基ずつずらしながらk-merで切り出し、k-mer配列、染色体番号、位置、向きを同一行に表記したデータを作成し、k-merの配列で辞書順に整列させた。当該データに対して、二分検索法により多型k-1merの配列のゲノム上の位置を決定した。
リファレンスゲノムにマップされた多型の確認を行った。上記工程で作成されたsort_uniq配列(ターゲット)と同様に、リファレンス配列も5’末端から2塩基ごとにリード長と同じ長さの配列を切り出し、相補鎖とともにsort、uniq処理をしたデータを作成した。このリファレンスとターゲットのsort_uniqデータを検索の対象とした。
次に、リファレンスゲノム配列から、多型位置を含むリード長と同じ長さの配列セット(検索を行うクエリセット)を切り出した。この配列セットは配列の5’末端から3’末端までのすべての位置で多型位置の塩基を含む配列セットである。配列セット数は配列長と同一となる。リファレンスの配列セット(リファレンスセット)に対して、多型位置の塩基を予想された多型塩基に置換した配列セット(ターゲットセット)を作成した。
リファレンスセット、ターゲットセットをクエリにして、リファレンスゲノム、及び、次世代シーケンサーから得られたFastq形式のファイルより出力されたファイル(reads)を辞書順にソートし同一の配列を一つにまとめたsort_uniq配列をそれぞれ、二分検索法で検索し、それぞれのセットに対してsort_uniq配列が何個マッチするかを調べた。リファレンスセットではリファレンスのsort_uniq配列のみにヒットする。これに対して、ターゲットセットでは、ホモ型の変異の場合は、ターゲットのsort_uniqデータのみから検出され、ヘテロ型の変異の場合は、リファレンスと、ターゲットのsort_uniqデータの両方から検出されるはずである。このようにして、予想と一致する検索値を示した多型を抽出した。
(ジェノタイプの確認)
ジェノタイプの確認は、PCRで当該領域を増幅し、サンガー法で決定した。
各領域の増幅に用いたプライマー配列は以下のとおりであった。
Figure 0007166638000023
PCR反応の反応条件は以下のとおりであった。
Figure 0007166638000024
反応サイクルは、94℃ 0.5分、60℃ 0.5分、72℃ 1分を30サイクルとして行った。
増幅したDNA断片を1%アガロースゲル電気泳動で分離し、0.5μg/mlの濃度のエチジウムブロマイドで染色して、長波長紫外線ランプ(365nm)で蛍光を発するバンドを切り出し、Promega社のWizard(登録商標) SV Gel and PCR Clean-Up System(Cat.# A9282)で断片の精製を行った。
サンガー法でのSNPを含む塩基配列の確認
精製された断片をBigDye(登録商標) Terminator v3.1 Cycle Sequencing Kit(Thermo Fisher Scientific Cat.# 4337455)で反応を行い、DNAシーケンサー ABI PRISM 3130xlで塩基配列の確認を行った。
(結果)
(シーケンシング)
イネ個体(N1、N1S5、N1S6、N1S7、N1S10)を次世代シーケンサーで解析した結果の塩基配列データは、DDBJに送信されており、以下のアクセッション番号で登録されている。
Figure 0007166638000025
各サンプルについてのリード数(総データ数)は、以下:
Figure 0007166638000026

のとおりであった。
出力されたFastqファイルを処理したsort_uniqのデータ数は以下:
Figure 0007166638000027

のとおりであった。sort_uniqは、Nを含まないリードとその相補鎖のデータをsortしたのちuniqで同一配列を一つにまとめたものである。このデータは、k-merのデータと異なり配列のみのデータであり、頻度の数値データは含まない。
(多型の検出)
コントロール配列:N1、対象配列:N1S7について、k-mer配列をリファレンス配列に沿って整列させ、多型の検出を行った結果は、図7および8に示される。下線を付された塩基がコントロールと対象で異なっており、多型が検出されたことが示されている。図9においては、コントロール配列(N1)と対象配列(N1S5、N1S6、N1S7、N1S10)のk-mer配列の部分配列サブセットを整列させることによって、イネリファレンス配列の対応する位置から始まる配列と対応する配列の出現頻度を求めた結果が示される。染色体番号、染色体の位置に続いて、N1、N1S5、N1S6、N1S7、N1S10の20-merの頻度が示される。N1S7でヘテロ、N1S10でミュータントホモになり、野生型の20-merがゼロになっていることがわかる。すなわち、k-mer配列の出現頻度の変化から、N1S7においてヘテロ変異が生じ、N1S10においてホモ変異となったことを検出することができた。
さらに、最終的に上記手順によって、コントロール配列(N1)と対象配列(N1S5、N1S6、N1S7、N1S10)との間で検出された多型の一部を、図10および11に示した。これらの結果は、リファレンス配列を用いた確認と一致していた。Wが野生型、Hがヘテロ型、Mがミュータントホモ型を示す。
サンガー法により確認した各サンプルにおける多型を、図12に示す。N1、N1S1、N1S2、N1S3、N1S4、N1S5、N1S6、N1S7、N1S8、N1S9、N1S10と世代を重ねながら、ヘテロ変異が生じ、その後ホモ変異として定着する様子が観察でき、この結果は、本発明の方法によって検出された多型とよく一致していた(図12)。
したがって、本発明の方法によって、世代間に生じた多型を詳細に検出できることが示された。また、Polymorphic Edge Detectionによって検出された多型が、サンガー法によっても確認されていることから、リファレンス配列(ゲノムリファレンス配列)を必要とせずに、配列データ間での多型の検出を行うことができることが実証された。
(実施例2:ナイジェリアのヨルバ族男性(NA18507)の配列解析)
(材料および方法)
コントロール配列データとして、ヒトゲノムリファレンスhg38を用いた。配列は、ftp://hgdownload.cse.ucsc.edu/goldenPath/hg38/chromosomes/よりchr1~chr22とchrX、chrY、chrMの染色体データをダウンロードして用いた。alt、v1等のファイル名にコメントが付いているデータは除外した。
対象配列データとしては、ヒトゲノムの次世代シーケンス配列データNA18507をダウンロードして用いた。この配列データは、Illumina社の次世代シーケンサーで解析が行われたものであり、NCBIに登録・公開されており、当該データをダウンロードして用いた。塩基配列セットの実験IDのURLは、https://www.ncbi.nlm.nih.gov/sra/SRX016231であり、配列のアクセッション番号は、SRR034939~SRR034975の範囲であった。
情報処理は、実施例1の(多型の検出)と同様に行った。
(結果)
結果の一部を、図13に示す。相補鎖(r)で検出された野生型・変異型の塩基は順鎖に変換して表示している。ヘテロ型の場合はそれぞれの塩基を並べて示している。参照、対象の配列xのA、C、G、Tの数はk-1配列に続く配列xの各塩基の頻度を示している。P~Q列に関しては、野生型あるいは変異型の塩基をもつ参照あるいは対象配列の数を示している。
順鎖、逆鎖の双方で同じ変異が検出された場合は、ほぼ間違いなくSNPである。片側の鎖のみで検出される原因は、逆の鎖の配列がゲノム上でユニークでなく、一意で検出できなかったため、あるいは、挿入、欠失、転座、逆位の境界塩基を検出したためであると考えられる。
本発明の方法は、ヒトゲノムについても、正確に多型を検出できることが実証された。また、コントロール配列として、データベースから取得したリファレンス配列を用いることができることも示される。
(実施例3:同一個体の組織間での多型の検出)
(概要)
本発明の方法により、同一個体の組織間での多型の検出が可能であることを実証する。
(材料および方法)
NCBIのSRAよりfastq-dumpを用いて配列データを取得した。本データは、Texas Cancer Research Biobank Open Access Data Sharing: Genome Projectが登録したデータであり、詳細データについて、以下のURL:https://www.ncbi.nlm.nih.gov/Traces/study/?acc=SRP060654において提供されている(Becnel, L. et al. NCBI Sequence Read Archive PRJNA284598 (2015).)。本発明の方法により、前記配列データにおいて多型を検出し、同一個体の正常組織と腫瘍組織との間の多型を本発明の方法によって検出できるかを検証した。
配列データの起源のサンプルは2種類で、配列データ名とサンプルの内容は以下の通りであった。
SRR2096532 コントロール血液 (Normal)
SRR2096535 濾胞性リンパ腫 (9690/3: Follicular lymphoma)
リード数 (配列長101塩基)
SRR2096532 1300353764
SRR2096535 1339310760
sort_uniqの配列数
SRR2096532 2056683322
SRR2096535 2181081390
SRR2096532(正常組織)をコントロール配列データとして用い、SRR2096535(腫瘍組織)を対象配列データとして用いて解析を行った。
情報処理は、実施例2の(多型の検出)と同様に行った。
なお、本実施例においては、確認(verify)工程で、検出されたそれぞれの変異に対して、ターゲットでは、変異型が5リード以上、野生型が1リード以下、コントロールでは変異型が1リード以下、野生型が5リード以上の場合、ホモ型変異(M)とマークし、ターゲットのリード総数に対する変異型の割合が0.3より大きいか、0.7より小さく、かつ、コントロールで変異型が1リード以下、かつ、ターゲットで野生型リードが5以上の場合、ヘテロ型変異(H)とした。
[従来法]
この分野で一般的に広く用いられている、Samtoolsを用いて同じデータを処理した。
従来法による解析は、以下の工程によって行った。工程0は準備で1回のみ実施し、工程1から5はサンプル毎に実施した。
0.準備:リファレンス配列にインデックス付加
1.ショートリードデータのマッピング
2.SAM形式をBAMに変換(マッピング位置でソートも)
3.Samtoolsで多型部位の検出
リファレンス配列データとして、ヒトゲノムリファレンスhg38を用いた。配列は、ftp://hgdownload.cse.ucsc.edu/goldenPath/hg38/chromosomes/よりchr1~chr22とchrX、chrY、chrMの染色体データをダウンロードして用いた。alt、v1等のファイル名にコメントが付いているデータは除外した。
(結果)
本発明の方法(PED)および従来法(bwa+Samtools)で検出された多型の数は以下のとおりであった。
Figure 0007166638000028
PEDでVerify(確認工程)にかけた座位数は22601で、そのうち順鎖、相補鎖共に置換変異が検出されたのは、514座位であった。Verifyの結果片側鎖のみで検出されたものの、ポジティブと判断された座位を合わせて1042座位に置換変異があると予想された。PEDで検出された1042座位はすべてヘテロ型であった。このことは、同一人物の血液と癌組織からの配列の比較であるため、原理的にホモ型変異が起こる可能性は非常に低いと考えられることと整合する。PEDではマッピングする前に直接対象とコントロールの比較を行っているため、対象・コントロールのSNP数は出力されない。
PEDで検出されたが、bwa+Samtoolsで検出できなかったSNPは20であった。bwa+Samtoolsではホモ型変異も多く検出されており、結果に非常に多くのノイズを含んでいると考えられる。
(考察)
同一人物からの組織間の比較なので、ミュータントホモはあり得ないと予想されていたところ、予想通りにヘテロ型のみが検出された。この精度で多型を検出できる系は、他にはなく、本発明の方法は従来技術に対して顕著に有利なものであると考えられる。
上記の結果から、同一個体の正常組織と癌組織の配列を直接比較して多型を検出することも可能であることが示された。本発明の方法は、いったんリファレンスゲノムにマップしてから差を調べる方法に比べて、多型の捕捉率、捕捉精度ともに非常に高くなる。
PEDではマッピングの前に、リード同士の比較からSNPを検出するので確度の高いSNPが得られる。本実施例で得られた結果を、挿入・欠失を検出するプログラムの結果と併せて考慮すると、癌細胞では、置換変異よりむしろ、二重鎖切断の後の除去修復に起因すると考えられる挿入・欠失変異の頻度が高いと考えられる。
(実施例4:コピー数多型の検出)
(材料および方法)
実施例1におけるN1S6と同じ世代の種子1粒から細胞培養を行い、1ヶ月、3ヶ月、5ヶ月後に再分化してイネの個体にした葉よりDNAを抽出し、それぞれ1M1、3M1、5M1のサンプルとして用いた。実施例1におけるN1種子と同じ世代の種子を5ヶ月培養して再分化した個体を4世代自殖した個体からDNAを抽出し、TTM2とTTM5のサンプルとして用いた。抽出したDNAから次世代シーケンサーによって配列データを取得した。シーケンシングのプロトコルは実施例1におけるものと同様であった。これらの配列データおよび実施例1のN1S5、N1S6、N1S7、N1S10の配列データを対象配列データとして、N1を参照配列データとして用いた。これらの配列データについてのアクセッション番号、リード数、sort_uniqの配列数は以下の表に示される。なお、TTM5のデータは、SRR556174とSRR556175の2つのアクセッション番号に分割されている。TTM5のsort_uniqは2つのリードを合わせて一つのファイルとして作成した。
Figure 0007166638000029
情報処理は、実施例1の(多型の検出)と同様に行い、参照配列データと、対象配列データの間で出現頻度が大きく異なったk-merを検出した。コントロールのN1より2倍程度以上の頻度が連続して検出された場合にCNVと判定した。誤検出も出てくるので、2.5倍以上と厳しくすることもできる。
(結果)
結果の一部を図17に示した。第7染色体の26694795位置(図7中矢印で示される)に対応するk-mer配列から、それまでの位置に対応するk-mer配列の出現頻度の4倍程度の値が、N1~N1S10で現れていることが分かる。
この位置からレトロトランスポゾンTos17の配列が始まる。このトランスポゾンはゲノムに2コピー存在しており、それぞれのトランスポゾンの末端部分は同じ配列(Long Terminal Repeat、LTR)を有する。そのため、第7染色体の26694795より、それまでの4倍程度の値が、N1~N1S10で現れたと考えられる。このトランスポゾンの全長は4.1kbあるため、図17には最初のジャンクションの部分のみ示している。
1M1、3M1、5M1はN1S6と同じ世代の種子1粒から細胞培養を行い、それぞれ1ヶ月、3ヶ月、5ヶ月後に再分化してイネの個体にした葉よりDNAを抽出して解析したものであり、図17に示される結果から、培養時間に応じてカウントが増えていることがわかる。これは、培養時間に比例してトランスポゾンが転移してコピー数が増加し、コピー数多型(CNV)が生じたためと考えられる。TTM2とTTM5はN1種子と同じ世代の種子を5ヶ月培養して再分化した個体を4世代自殖した個体のDNAであり、こちらも同様にコピー数の増加が認められる。
(考察)
Tos17は培養時のみ活性化されて転移するイネのトランスポゾンとして知られている。Tos17はレトロトランスポゾンなので、オリジナルは切り出されて転移することなく、Tos17のコピーの転移によってゲノム上のTos17のコピー数が増加する。そのため、Tos17は培養で転移してコピー数が増えることが以前から知られている。
本発明の方法によって、コピー数の変異を検出した結果、理論とよく一致して、培養時間が長くなるにつれて、Tos17のコピー数が増加していることが図17に示されるカウントから理解される。
したがって、本発明の方法によって、コピー数多型を検出することが可能であることが実証される。
例えば、このようなコピー数多型の検出は、培養細胞(例えば、iPS細胞等)において、品質の管理に用いることが可能であると考えられる。上記イネの培養細胞と同様にトランスポゾン等によるコピー数変異が、例えば、ヒトの培養細胞(iPS等)で観察されている場合には、例えば、治療に用いるのは危険である可能性が高いという判断に用いることが可能である。
(関連出願)
本出願は、2017年7月24日に出願された特願2017-142781号の優先権の利益を主張し、当該出願は、全ての目的において、その開示全体が本明細書において参考として援用される。さらに、本明細書において、本出願人により本出願と同日に出願された「挿入・欠失・逆位・転座・置換検出法」との名称の出願(整理番号NG012PCT/F5-18PCT075)およびその基礎出願である2017年7月24日に出願された特願2017-142782号(整理番号J1-17369162)は、全ての目的において、その開示全体が本明細書において参考として援用される。
塩基配列解析で多型を検出するすべての分野で利用が可能で、DNA育種利用の他、臨床検査、iPS細胞の検査、メタゲノム解析、発現解析等、幅広い分野で利用することができる。
配列番号1~60:図7のk-mer参照配列
配列番号61~80:図7のk-mer対象配列(変異が存在する部分)
配列番号81~140:図8のk-mer参照配列
配列番号141~160:図8のk-mer対象配列(変異が存在する部分)
配列番号161~190:図13のk-1(k=20)配列
配列番号191~221:(具体的な例)で用いられた配列
配列番号222~232:(例示的計算フロー)で用いられた配列
配列番号233~266:実施例1で用いられたプライマーの配列
配列番号267~275:図11の配列

Claims (47)

  1. 対象配列データにおいてコントロール配列データに対する多型をコンピュータに検出させる方法であって、該コンピュータが、
    a)該対象配列データの長さkの部分配列のサブセットの各々の部分配列の出現頻度を提供する工程と、
    b)該コントロール配列データの長さkの部分配列のサブセットの各々の部分配列の出現頻度を提供する工程と、
    c)対象配列とコントロール配列とを比較し、該出現頻度の分布の比較に基づいて、多型を検出する工程と
    を包含し、ここで、kは該対象配列および該コントロール配列のいずれか短いほうの全長以下の整数である、方法。
  2. 前記部分配列中の長さk-xの配列部分が共通する配列ごとに、長さxの部分について出現頻度の分布を算出する工程をさらに含み、ここで、xはk未満の正の整数である、請求項1に記載の方法。
  3. 前記比較が、前記部分配列中の長さk-xの配列部分が共通する配列における、長さxの部分の出現頻度の分布の差異の比較を含む、請求項2に記載の方法。
  4. 前記部分配列中の長さk-xの配列部分を、ユニークな配列ごとにグルーピングする工程を含み、ここで、xはk未満の正の整数である、請求項1~3のいずれか1項に記載の方法。
  5. 前記長さk-xの配列部分をソートする工程を含む、請求項4に記載の方法。
  6. 前記長さk-xの配列部分を文字列としてソートする工程を含む、請求項5に記載の方法。
  7. 前記kが、前記対象配列における偶然同一を排除する長さである、請求項1~6のいずれか1項に記載の方法。
  8. 前記対象配列データおよび前記コントロール配列データが、生物のゲノムに由来する塩基配列データであり、前記kが、前記生物のゲノムにおいて、異なる箇所での偶然同一を排除する長さである、請求項1~7のいずれか1項に記載の方法。
  9. 長さxが1~2である、請求項2~6のいずれか1項、または請求項7および8が請求項2~6のいずれか1項に従属する場合の請求項7および8のいずれか1項に記載の方法。
  10. 長さxが1である、請求項9に記載の方法。
  11. 前記長さxの部分が、前記部分配列の末端に存在する、請求項2~10のいずれか1項に記載の方法。
  12. 前記対象配列データおよび前記コントロール配列データが塩基配列データであり、前記長さxの部分が、前記部分配列の3’末端である、請求項9に記載の方法。
  13. 前記コントロール配列データのサブセットと前記対象配列データのサブセットとの間で、前記長さxの部分の配列の出現頻度が異なる場合、該長さxの部分の配列を、対象配列データにおけるコントロール配列データに対する多型として検出する、請求項2~12のいずれか1項に記載の方法。
  14. 前記コントロール配列データのサブセットと前記対象配列データのサブセットとの間で、前記長さxの部分の配列で最も高頻度のものが異なっている長さk-xの配列部分が存在する場合、該長さxの部分の配列を、対象配列データにおける多型として検出する、請求項2~13のいずれか1項に記載の方法。
  15. 前記対象配列データおよび前記コントロール配列データが塩基配列データであり、前記対象配列データのサブセットにおける前記長さxの部分の配列で、前記コントロール配列データのサブセットにおける最も高頻度のものと同一の長さxの部分の配列がノイズ以下のカウントしか存在しない長さk-xの配列部分が存在する場合、該長さxの部分の配列を、対象配列データにおけるホモ多型として検出する、請求項2~14のいずれか1項に記載の方法。
  16. 前記対象配列データおよび前記コントロール配列データが塩基配列データであり、対象サブセットにおける前記長さxの部分の配列で、コントロール配列データのサブセットにおける最も高頻度のものと同一の長さxの部分の配列が存在し、かつ、コントロール配列データのサブセットにおける最も高頻度のものと異なる長さxの部分の配列が存在する長さk-xの配列部分が存在する場合、該長さxの部分の配列を、対象配列データにおけるヘテロ多型として検出する、請求項2~15のいずれか1項に記載の方法。
  17. 対象配列データ量から予測される出現頻度と比較して、前記出現頻度が少ない部分配列をノイズとする、請求項1~16のいずれか1項に記載の方法。
  18. 前記対象配列データおよび前記コントロール配列データが塩基配列データであり、[(対象配列データ量)×(1-精度)]/(対象ゲノムサイズ)+1未満の出現頻度の部分配列をノイズとする、請求項17に記載の方法。
  19. 前記対象配列データが、次世代シーケンシングによって得られた塩基配列データである、請求項1~18のいずれか1項に記載の方法。
  20. 前記対象配列データが、個体から得られた配列データであり、前記コントロール配列データが、該個体と同種の別の個体、またはデータベースから得られた配列データである、請求項1~19のいずれか1項に記載の方法。
  21. 前記対象配列データが、個体の組織試料から得られた配列データであり、前記コントロール配列データが、該個体の別の組織、またはデータベースから得られた配列データである、請求項1~20のいずれか1項に記載の方法。
  22. 前記対象配列データが、細胞試料から得られた配列データであり、前記コントロール配列データが、別の細胞、またはデータベースから得られた配列データである、請求項1~21のいずれか1項に記載の方法。
  23. 前記多型が、置換、挿入、欠失、コピー数多型(Copy Number Variation, CNV)、STRP(short tandem repeat polymorphism)、逆位または転座である、請求項1~22のいずれか1項に記載の方法。
  24. 前記多型が、置換である、請求項23に記載の方法。
  25. 前記対象配列に対するリファレンス配列における前記多型の位置を特定する工程をさらに含む、請求項1~24のいずれか1項に記載の方法。
  26. 前記対象配列データおよび前記コントロール配列データが、生物のゲノムに由来する塩基配列データであり、前記多型のゲノム上の位置を特定する工程をさらに含む、請求項1~25のいずれか1項に記載の方法。
  27. 検出された多型の部位について、リファレンス配列またはコントロール配列から作成したクエリ配列セットを用いて、対象配列データおよび/またはコントロール配列データとの比較を行い確認する工程をさらに含む、請求項25または26に記載の方法。
  28. 前記クエリ配列セットが、リファレンス配列またはコントロール配列において前記多型に該当する部位の文字を異なる文字に置換した変異型クエリ配列セットを含む、請求項27に記載の方法。
  29. 前記対象配列データおよび前記コントロール配列データが塩基配列データであり、検出された多型の部位について、対象配列データおよび/またはコントロール配列データの相補鎖の配列データに対して、リファレンス配列またはコントロール配列から作成したクエリ配列セットとの比較を行い確認する工程をさらに含む、請求項27または28に記載の方法。
  30. 前記対象配列データおよび前記コントロール配列データが塩基配列データであり、検出された多型の部位について、対象配列データおよび/またはコントロール配列データの変異型の塩基を有する配列データに対して、リファレンス配列またはコントロール配列から作成したクエリ配列セットとの比較を行い確認する工程をさらに含む、請求項27~29のいずれか1項に記載の方法。
  31. 前記対象配列データおよび前記コントロール配列データが、生物のゲノムに由来する塩基配列データであり、前記ゲノムの配列が不明である、請求項1~30のいずれか1項に記載の方法。
  32. 実験結果またはデータベースから対象配列データまたはコントロール配列データを取得する工程をさらに含む、請求項1~31のいずれか1項に記載の方法。
  33. 対象配列データにおけるコントロール配列データに対する多型を含む部分配列中の多型ではない部分の少なくとも一部を含む配列を、該多型の識別子として割り当てることをさらに含む、請求項1~32のいずれか1項に記載の方法。
  34. 前記多型の識別子をリファレンス配列にマッピングし、リファレンス上の該多型の位置を特定することを含む、請求項33に記載の方法。
  35. 対象配列データにおいてコントロール配列データに対する多型を検出する方法をコンピュータに実行させるためのプログラムであって、該方法は、
    a)該対象配列データの長さkの部分配列のサブセットをコンピュータに保存する工程であって、ここで、kは対象配列およびコントロール配列の全長以下の長さである、工程と、
    b)該対象配列データの長さkのサブセットの各々の部分配列の出現頻度を算出する工程と、
    c)該コントロール配列データの長さkの部分配列のサブセットにおける各々の部分配列の出現頻度をコンピュータに保存する工程と、
    d)対象配列とコントロール配列とを比較し、該出現頻度の分布の差異に基づいて、多型を検出する工程と
    を包含する、プログラム。
  36. 前記方法が、前記部分配列中の多型ではない部分の少なくとも一部を含む配列を、検出された前記多型の名称として表示する工程をさらに含む、請求項35に記載のプログラム。
  37. 対象配列データにおいてコントロール配列データに対する多型を検出する方法をコンピュータに実行させるためのプログラムを格納する記録媒体であって、該方法は、
    a)該対象配列データの長さkの部分配列のサブセットをコンピュータに保存する工程であって、ここで、kは対象配列およびコントロール配列の全長以下の長さである、工程と、
    b)該対象配列データの長さkのサブセットの各々の部分配列の出現頻度を算出する工程と、
    c)該コントロール配列データの長さkの部分配列のサブセットにおける各々の部分配列の出現頻度をコンピュータに保存する工程と、
    d)対象配列とコントロール配列とを比較し、該出現頻度の分布の差異に基づいて、多型を検出する工程と
    を包含する、記録媒体。
  38. 前記方法が、前記部分配列中の多型ではない部分の少なくとも一部を含む配列を、検出された前記多型の名称として表示する工程をさらに含む、請求項37に記載の記録媒体。
  39. 対象配列データにおいてコントロール配列データに対する多型を検出するためのシステムであって、該システムは、
    該対象配列データおよび該コントロール配列データの長さkの部分配列のサブセットの各々の部分配列の出現頻度を提供するように構成された配列データ処理部であって、ここで、kは対象配列およびコントロール配列の全長以下の長さである、配列データ処理部と、
    対象配列とコントロール配列とを比較し、該出現頻度の分布の差異に基づいて、多型を検出する工程とを行うように構成された、配列データ計算部と
    を備える、システム。
  40. 前記システムが、前記部分配列中の多型ではない部分の少なくとも一部を含む配列を、検出された前記多型の名称として表示する表示手段をさらに含む、請求項39に記載のシステム。
  41. 部分配列を多型の識別子としてコンピュータに利用させる方法であって、該コンピュータが、
    請求項1~34のいずれか一項に記載の方法によって検出された、前記対象配列データにおける前記コントロール配列データに対する多型を含む部分配列中の多型ではない部分の少なくとも一部を含む配列を、多型の識別子として割り当てる工程を含む、方法。
  42. 部分配列を多型の識別子としてコンピュータに利用させる方法であって、該コンピュータが、
    請求項1~34のいずれか一項に記載の方法によって検出された、前記対象配列データにおける前記コントロール配列データに対する多型を含む部分配列中の多型ではない部分を含む配列を、多型の識別子として割り当てる工程を含む、方法。
  43. 部分配列を多型の識別子としてコンピュータに利用させる方法であって、該コンピュータが、
    請求項1~34のいずれか一項に記載の方法によって検出された、前記対象配列データにおける前記コントロール配列データに対する多型を含む部分配列を、多型の識別子として割り当てる工程を含む、方法。
  44. 前記識別子は、前記多型の識別子のマッピング、連鎖解析、形質との関係を調べる多型マーカーとしての使用、診断、育種、鑑定、品質管理、分類および/または検査において使用される、請求項41~43のいずれか一項に記載の方法。
  45. 前記部分配列は、長さkの部分配列またはk-xの部分配列のいずれかである、請求項41~44のいずれか一項に記載の方法。
  46. 前記部分配列は、前記対象配列データおよび前記コントロール配列データにおいてユニークな配列である、請求項41~45のいずれか一項に記載の方法。
  47. 前記多型の識別子をリファレンス配列にマッピングし、リファレンス上の該多型の位置を特定することを含む、請求項41~46のいずれか一項に記載の方法。
JP2019532603A 2017-07-24 2018-07-23 多型検出法 Active JP7166638B2 (ja)

Applications Claiming Priority (3)

Application Number Priority Date Filing Date Title
JP2017142781 2017-07-24
JP2017142781 2017-07-24
PCT/JP2018/027535 WO2019022018A1 (ja) 2017-07-24 2018-07-23 多型検出法

Publications (2)

Publication Number Publication Date
JPWO2019022018A1 JPWO2019022018A1 (ja) 2020-05-28
JP7166638B2 true JP7166638B2 (ja) 2022-11-08

Family

ID=65039682

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2019532603A Active JP7166638B2 (ja) 2017-07-24 2018-07-23 多型検出法

Country Status (3)

Country Link
JP (1) JP7166638B2 (ja)
TW (1) TW201920682A (ja)
WO (1) WO2019022018A1 (ja)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115910197B (zh) * 2021-12-29 2024-03-22 上海智峪生物科技有限公司 基因序列处理方法、装置、存储介质及电子设备

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2001112486A (ja) 1999-08-05 2001-04-24 Takeda Chem Ind Ltd 遺伝子解析結果の記録方法
JP2001167123A (ja) 1999-12-13 2001-06-22 Iyaku Bunshi Sekkei Kenkyusho:Kk 遺伝子配列多型クラスの管理方法
JP2008165375A (ja) 2006-12-27 2008-07-17 Canon Inc 塩基配列を識別する変異セットの選別法
JP2016504667A (ja) 2012-11-26 2016-02-12 コーニンクレッカ フィリップス エヌ ヴェKoninklijke Philips N.V. 患患者固有の関連性評価を用いた変異と疾患の関連付けを使用する診断的遺伝子分析
JP2017045451A (ja) 2012-03-29 2017-03-02 三菱レイヨン株式会社 βグロビン遺伝子の変異を検出するためのマイクロアレイ及びその検出方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2001112486A (ja) 1999-08-05 2001-04-24 Takeda Chem Ind Ltd 遺伝子解析結果の記録方法
JP2001167123A (ja) 1999-12-13 2001-06-22 Iyaku Bunshi Sekkei Kenkyusho:Kk 遺伝子配列多型クラスの管理方法
JP2008165375A (ja) 2006-12-27 2008-07-17 Canon Inc 塩基配列を識別する変異セットの選別法
JP2017045451A (ja) 2012-03-29 2017-03-02 三菱レイヨン株式会社 βグロビン遺伝子の変異を検出するためのマイクロアレイ及びその検出方法
JP2016504667A (ja) 2012-11-26 2016-02-12 コーニンクレッカ フィリップス エヌ ヴェKoninklijke Philips N.V. 患患者固有の関連性評価を用いた変異と疾患の関連付けを使用する診断的遺伝子分析

Also Published As

Publication number Publication date
TW201920682A (zh) 2019-06-01
WO2019022018A1 (ja) 2019-01-31
JPWO2019022018A1 (ja) 2020-05-28

Similar Documents

Publication Publication Date Title
Minnoye et al. Chromatin accessibility profiling methods
KR102662186B1 (ko) 임신 중 긴 세포유리 단편을 사용한 분자 분석
US20060286566A1 (en) Detecting apparent mutations in nucleic acid sequences
Lange et al. Analysis pipelines for cancer genome sequencing in mice
Larson et al. A clinician’s guide to bioinformatics for next-generation sequencing
US20190139628A1 (en) Machine learning techniques for analysis of structural variants
Kockum et al. Overview of genotyping technologies and methods
Lindner et al. Performance of methods to detect genetic variants from bisulphite sequencing data in a non‐model species
Goswami et al. RNA-Seq for revealing the function of the transcriptome
CN111826429B (zh) 一种基于简化基因组测序和snp次等位基因频率的非杂交后代鉴定方法
JP7166638B2 (ja) 多型検出法
JP2022549823A (ja) キットおよびキットの使用方法
US20230332205A1 (en) Linked dual barcode insertion constructs
JP7122006B2 (ja) 挿入・欠失・逆位・転座・置換検出法
KR20170046315A (ko) 차세대 시퀀싱을 이용한 분자마커 개발 방법
Fletcher et al. AFLAP: Assembly-Free Linkage Analysis Pipeline using k-mers from whole genome sequencing data
AU780824B2 (en) DNA marker profile data analysis
JP7444488B2 (ja) 混入検出法
Orr Methods for detecting mutations in non-model organisms
WO2022168195A1 (ja) 遺伝情報解析システム、及び遺伝情報解析方法
Porter Mapping bisulfite-treated short DNA reads
Conceição Differential DNA Methylation in Aging: in Silico Exploration Using High-Throughput Datasets
Lübberstedt et al. Markers and Sequencing
Data Check for updates Chapter 7 Genetic Diversity Assessment in Plants from Reduced Representation Sequencing Data José P. Jiménez-Madrigal, Bradley J. Till®, and Andrés Gatica-Arias
Pala Sequence Variation Of Copy Number Variable Regions In The Human Genome

Legal Events

Date Code Title Description
A621 Written request for application examination

Free format text: JAPANESE INTERMEDIATE CODE: A621

Effective date: 20210304

A131 Notification of reasons for refusal

Free format text: JAPANESE INTERMEDIATE CODE: A131

Effective date: 20220502

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20220630

TRDD Decision of grant or rejection written
A01 Written decision to grant a patent or to grant a registration (utility model)

Free format text: JAPANESE INTERMEDIATE CODE: A01

Effective date: 20221014

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20221019

R150 Certificate of patent or registration of utility model

Ref document number: 7166638

Country of ref document: JP

Free format text: JAPANESE INTERMEDIATE CODE: R150

RD04 Notification of resignation of power of attorney

Free format text: JAPANESE INTERMEDIATE CODE: R3D04