(57)【要約】塩基配列における塩基種を推定する方法、装置、プログラムを提供する。前記方法では、各核酸分子由来のリード配列を取得し（Ｓ１１）、取得したリード配列を参照配列上にマッピングし（Ｓ１２）、マッピング結果に基づきリード配列を分類し（Ｓ１３）、分類結果に基づきコンセンサスリードとコンセンサスリード中の各塩基に対応する信頼性スコアを求め（Ｓ１４）、コンセンサスリードより特定位置における塩基種を任意に選択し、コンセンサスリード総数に対する、選択した特定位置の各塩基種を含むコンセンサスリード数を取得し、選択した特定位置の塩基種に対応する信頼性スコアを用いたパラメトリックなエラーモデルにより（Ｓ１５）、設定された信頼性で特定位置における塩基種の解析エラーによる出現率を推定し、出現率に基づき、特定位置における特定の塩基種が、設定された信頼性でクローナル核酸混合物中に存在するか否かを推定する（Ｓ１６）。
2019132010.tif 000002
【特許請求の範囲】
【請求項１】
コンピュータを用いて、クローナル核酸混合物を構成する核酸分子の塩基配列における特定位置に存在する特定塩基種の存在を推定する方法であって、
ａ）各核酸分子に由来するリード配列を取得する工程、
ｂ）工程ａで取得したリード配列を参照配列上にマッピングしマッピング結果を取得する工程、
ｃ）工程ｂで得られたマッピング結果を基にリード配列のクラスタリング結果を取得する工程、
ｄ）クラスタリング結果を基にコンセンサスリードとコンセンサスリード中の各塩基に対応する信頼性スコアからなるコンセンサスリード情報を取得する工程、
ｅ）クローナル核酸混合物を構成する各核酸分子から工程ａ〜ｄを経て得られるコンセンサスリードの集合に含まれる各コンセンサスリードの配列情報より、参照配列上の特定位置に対応する位置に存在する塩基種とその信頼性スコアを選抜する工程、
ｆ）工程ｅで使用した特定位置を含むコンセンサスリードの総数に対する、工程ｅで選抜した特定位置に特定の各塩基種を含有するコンセンサスリード数の比率を取得する工程、
ｇ）工程ｄで取得した信頼性スコアを用いたパラメトリックなエラーモデルにより、設定された閾値または有意水準での工程ｅで選抜した特定位置において解析エラーによって特定の塩基が出現する出現率を推定する工程、
ｈ）工程ｆで得られた比率が工程ｇにより得られた出現率より有意に高い場合に、工程ｅで選抜した特定位置に特定の塩基種を有する核酸分子が、工程ｇで設定された閾値または有意水準でクローナル核酸混合物中に存在すると推定する工程、
を含む方法。
【請求項２】
コンピュータを用いて、クローナル核酸混合物を構成する核酸分子の塩基配列における特定位置に存在する特定塩基種がクローナル核酸混合物中で存在する割合を推定する方法であって、
ａ）各核酸分子に由来するリード配列を取得する工程、
ｂ）工程ａで取得したリード配列を参照配列上にマッピングしマッピング結果を取得する工程、
ｃ）工程ｂで得られたマッピング結果を基にリード配列のクラスタリング結果を取得する工程、
ｄ）クラスタリング結果を基にコンセンサスリードとコンセンサスリード中の各塩基に対応する信頼性スコアからなるコンセンサスリード情報を取得する工程、
ｅ）クローナル核酸混合物を構成する各核酸分子から工程ａ〜ｄを経て得られるコンセンサスリードの集合に含まれる各コンセンサスリードの配列情報より、参照配列上の特定位置に対応する位置に存在する塩基種とその信頼性スコアを選抜する工程、
ｆ）工程ｅで使用した特定位置を含むコンセンサスリードの総数に対する、工程ｅで選抜した特定位置に特有の各塩基種を含有するコンセンサスリード数の比率を取得する工程、
ｇ）工程ｄで取得した信頼性スコアを用いたパラメトリックなエラーモデルにより、設定された閾値または有意水準での工程ｅで選抜した特定位置において解析エラーによって特有の塩基が出現する出現率を推定する工程、
ｈ）工程ｆで得られた比率が工程ｇにより得られた出現率より有意に高い場合に、工程ｅで選抜した特定位置に特定の塩基種を有する核酸分子が、工程ｇで設定された閾値または有意水準でクローナル核酸混合物中に存在すると推定する工程、及び
ｉ）工程ｈにおいて工程ｅで選抜した塩基種が存在すると判断された場合、工程ｆで得られた特定位置において塩基種を含有するコンセンサスリードの数から、工程ｇで設定された閾値または有意水準における工程ｅで選抜した特定位置における特定塩基種の比率を算出し、クローナル核酸混合物中で存在する割合として採用する工程、
を含む方法。
【請求項３】
コンピュータを用いて、クローナル核酸混合物を構成する各核酸分子の特定領域の塩基配列推定方法であって、
ａ）各核酸分子に由来するリード配列を取得する工程、
ｂ）工程ａで取得したリード配列を参照配列上にマッピングしマッピング結果を取得する工程、
ｃ）工程ｂで得られたマッピング結果を基にリード配列のクラスタリング結果を取得する工程、
ｄ）クラスタリング結果を基にコンセンサスリードとコンセンサスリード中の各塩基に対応する信頼性スコアからなるコンセンサスリード情報を取得する工程、
ｅ）クローナル核酸混合物を構成する各核酸分子から工程ａ〜ｄを経て得られるコンセンサスリードの集合に含まれる各コンセンサスリードの配列情報より、参照配列上の特定位置に対応する位置に存在する塩基種とその信頼性スコアを選抜する工程、
ｆ）工程ｅで使用した特定位置を含むコンセンサスリードの総数に対する、工程ｅで選抜した特定位置に特定の各塩基種を含有するコンセンサスリード数の比率を取得する工程、
ｇ）工程ｅで選抜した特定位置の塩基種に対応する工程ｄで取得した信頼性スコアを用いたパラメトリックなエラーモデルにより、設定された閾値または有意水準での工程ｅで確認した特定位置における塩基種の解析エラーによる出現率を推定する工程、
ｈ）工程ｆで得られた比率が工程ｇにより得られた出現率より有意に高い場合に、工程ｅで選抜した特定位置に特定の塩基種を有する核酸分子が、工程ｇで設定された閾値または有意水準でクローナル核酸混合物中に存在する、あるいは該出現率より低い場合は特定位置の塩基種は不明、と判断する工程、
ｉ）工程ｈで得られた結果を該塩基の各コンセンサスリードに適用し二次配列情報を作成する工程、及び
ｊ）工程ｅにおける塩基位置を変更し、工程ｅ〜ｉを繰り返すことにより得られた特定領域の全二次配列情報を統合し、クローナル核酸混合物を構成する各核酸分子の特定領域の塩基配列として採用する工程、
を含む方法。
【請求項４】
工程ｂで採用した参照配列とコンセンサスリード上の塩基配列が同じである場合は工程ｅ〜工程ｈを経ずにそのまま工程ｉで作成する二次配列情報の塩基種として採用することを特徴とする請求項３記載の塩基配列決定方法。
【請求項５】
工程ｅで選抜する特定の塩基位置が、工程ｂで採用した参照配列とコンセンサスリード上の塩基種が異なる位置であることを特徴とする請求項１〜３のいずれかに記載の方法。
【請求項６】
核酸分子が少なくとも一方の末端に分子バーコードを有することを特徴とする請求項１〜３のいずれかに記載の方法。
【請求項７】
各核酸分子が異なる分子バーコードを有することを特徴とする請求項６に記載の方法。
【請求項８】
各核酸分子が一つ以上の異なる分子バーコードを有することを特徴とする請求項７記載の方法。
【請求項９】
工程ｃにおけるクラスタリングをリード配列における分子バーコード領域の配列を指標に実施することを特徴とする請求項１〜８のいずれかに記載の方法。
【請求項１０】
分子バーコードがランダムな並びの塩基配列であることを特徴とする請求項９記載の方法。
【請求項１１】
分子バーコード領域のリード配列におけるミスマッチが、該分子バーコード領域のリード配列間の配列類似度が許容範囲内であることを特徴とする請求項１０記載の方法。
【請求項１２】
信頼性スコアが、クラスタリングされたリード配列の情報から、ベイズアプローチによって取得されることを特徴とする請求項１〜１１のいずれかに記載の方法。
【請求項１３】
リード配列が核酸分子の増幅産物のリード配列であることを特徴とする請求項１２に記載の方法
【請求項１４】
増幅物がポリメラーゼ連鎖反応により得られることを特徴とする請求項１３記載の方法。
【請求項１５】
ベイズアプローチの尤度関数がシーケンサーエラーと増幅エラーを考慮した関数であることを特徴とする請求項１３または１４に記載の方法。
【請求項１６】
工程ｂにおけるマッピング工程において参照配列に対しリード配列で異なる塩基が挿入または欠失に相当すると判断された場合、シーケンサーエラーに関し更に補正を加えた関数を採用することを特徴とする請求項１５記載の方法。
【請求項１７】
工程ｇのパラメトリックなエラーモデルがポアソン分布に基づくことを特徴とする請求項１〜３いずれかに記載の方法。
【請求項１８】
特定位置の特定塩基種が低頻度な変異の対象塩基であることを特徴とする請求項１または２記載の方法。
【請求項１９】
工程ｃから工程ｄへ直接進むことを特徴とする請求項１〜３のいずれかに記載の方法。
【請求項２０】
工程ｃの次に、
ｋ）工程ｃで得られた各クラスタリングに含まれるリード配列を用いて、ローカルアセンブリを実行し、アセンブル配列のクラスタリング結果を取得する工程、
を行うことを特徴とする請求項１〜３のいずれかに記載の方法。
【請求項２１】
請求項１９及び請求項２０に記載の方法によって得られた結果を基に、クローナル核酸混合物を構成する拡散分子の塩基配列における特定位置に存在する特定塩基種を決定する工程を含む方法。
【請求項２２】
請求項１〜２１のいずれかに記載の方法をコンピュータで実行させるためのコンピュータ読み取り可能なプログラム。
【請求項２３】
クローナル核酸混合物を構成する核酸分子の塩基配列における特定位置の特定塩基種の存在を推定するための装置であって、
各核酸分子由来のリード配列を取得するリード配列取得部と、
取得したリード配列をマッピング装置に供給し、前記マッピング装置による参照配列上へのマッピング結果を取得するマッピング情報取得部と、
マッピングしたリード配列のクラスタリング結果を得るクラスタリング結果取得部と、
クラスタリング結果を基にコンセンサスリードとコンセンサスリード中の各塩基に対応する信頼性スコアからなるコンセンサスリード情報取得部と、
を備えた
推定装置。
【請求項２４】
クローナル核酸混合物を構成する核酸分子の塩基配列における特定位置の特定塩基種の存在を推定するための装置であって、
各核酸分子由来のリード配列を取得するリード配列取得部と、
取得したリード配列をマッピング装置に供給し、前記マッピング装置による参照配列上へのマッピング結果を取得するマッピング情報取得部と、
マッピングしたリード配列のクラスタリング結果を得る第１のクラスタリング結果取得部と、
各クラスタリングに含まれるリード配列を用いて、ローカルアセンブリを実行し、アセンブル配列のクラスタリング結果を得る第２のクラスタリング結果取得部と、
第２のクラスタリング結果取得部のクラスタリング結果を基にコンセンサスリードとコンセンサスリード中の各塩基に対応する信頼性スコアからなるコンセンサスリード情報取得部と、
を備えた
推定装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、塩基配列における塩基種を推定する技術に関する。
【背景技術】
【０００２】
がん細胞は、正常な細胞中のがん遺伝子及びがん抑制遺伝子（がん関連遺伝子）に塩基配列の突然変異が生じることにより発生すると言われている。近年、これらのがん関連遺伝子の遺伝子異常を応用して、正常な細胞とがん細胞を見分ける手法の開発が進められている。例えば、非特許文献１および非特許文献２では、血中や便中に存在する微量ながん細胞由来の異常なＤＮＡを検出できたことが開示されている。さらに、がん関連遺伝子の変異によって翻訳された異常タンパク質を標的とする分子標的薬の登場により、塩基配列を高精度に決定することが求められる。【０００３】
次世代シーケンサーによる腫瘍検体の変異の包括的なプロファイリングの結果、一つの腫瘍のなかには、正常な細胞とは異なる配列のゲノムを持つ、がん細胞由来の複数のゲノムクローンが存在することが報告され、この現象は腫瘍内不均一性と呼ばれている。この現象により、腫瘍自体がポリクローナルな混合状態となっていることから、低頻度な変異も含めて網羅的に検出する必要性がある。【０００４】
現在の次世代シーケンス技術では、核酸増幅時及びシーケンシング時にある頻度でエラーが生じる。したがって、シーケンス解析を実施したゲノムＤＮＡの塩基配列には技術的なエラーに起因する偽陽性の変異（ｆａｌｓｅ  ｐｏｓｉｔｉｖｅ  ｖａｒｉａｎｔ）が含まれるため、１〜５％以下の低頻度な変異の検出は難しい。しかし近年、分子バーコード技術を応用することによって、低頻度な変異検出時のエラーが抑制されることが報告されている。【０００５】
例えば、非特許文献３には、配列解析装置より得られた塩基配列情報の統計処理により検出精度を高める技術としてｓｍＣｏｕｎｔｅｒが開示されているが、リキッドバイオプシーなど、より高い検出感度と精度を求められる場合には精度が足らず、不十分である。
【先行技術文献】
【非特許文献】
【０００６】
【非特許文献１】Ａｒｍｅｎｇｏｌ  Ｇ等、Ｊ  Ｍｏｌ  Ｄｉａｇｎ．  ２０１６  Ｊｕｌ；１８（４）：４７１−９．  ｐｍｉｄ：  ２７１５５０４８
【非特許文献２】Ｄｉｅｈｌ  Ｆ等、Ｎａｔ  Ｍｅｄ．  ２００８；１４（９）：  ９８５-９０．  ｐｍｉｄ：１８６７０４２２
【非特許文献３】Ｃｈａｎｇ  Ｘｕ等、ＢＭＣ  Ｇｅｎｏｍｉｃｓ、第１８巻、ｅ５ページ、２０１７年（ｄｏｉ：１０．１１８６／ｓ１２８６４−０１６−３４２５−４））
【発明の概要】
【発明が解決しようとする課題】
【０００７】
本発明は、塩基配列における塩基種を推定する方法、装置、及びその方法を実施するためのコンピュータープログラムを提供することを目的とする。
【課題を解決するための手段】
【０００８】
本発明者は、テンプレートＤＮＡに分子バーコードを付加させて調製したＤＮＡ断片から得られた次世代シーケンサーの多量の塩基配列データ（リード情報）から、分子バーコード毎に信頼性スコアを決定し、パラメトリック解析することで、従来法よりも高い感度で変異を検出できることを見出した。【０００９】
また、本発明者は、シーケンスデータから低頻度に存在する、特定位置の塩基種を推定するために、核酸増幅又はシーケンシングで発生するエラーを効率的に取り除く指標を算出し、真の塩基種を高確度に判定する統計手法を開発した。【００１０】
本発明の検出精度を高める技術は、例えば、より疾患早期ステージでの変異体核酸の検出を期待できる。【００１１】
本開示の第一の態様において、コンピュータを用いて、クローナル核酸混合物を構成する核酸分子の塩基配列における特定位置に存在する特定塩基種の存在を推定する第一の推定方法が提供される。
第一の推定方法は、
ａ）各核酸分子に由来するリード配列を取得する工程、
ｂ）工程ａで取得したリード配列を参照配列上にマッピングしマッピング結果を取得する工程、
ｃ）工程ｂで得られたマッピング結果を基にリード配列のクラスタリング結果を取得する工程、
ｄ）クラスタリング結果を基にコンセンサスリードとコンセンサスリード中の各塩基に対応する信頼性スコアからなるコンセンサスリード情報を取得する工程、
ｅ）クローナル核酸混合物を構成する各核酸分子から工程ａ〜ｄを経て得られるコンセンサスリードの集合に含まれる各コンセンサスリードの配列情報より、参照配列上の特定位置に対応する位置に存在する塩基種とその信頼性スコアを選抜する工程、
ｆ）工程ｅで使用した特定位置を含むコンセンサスリードの総数に対する、工程ｅで選抜した特定位置に特定の各塩基種を含有するコンセンサスリード数の比率を取得する工程、
ｇ）工程ｄで取得した信頼性スコアを用いたパラメトリックなエラーモデルにより、設定された閾値または有意水準での工程ｅで選抜した特定位置において解析エラーによって特定の塩基が出現する出現率を推定する工程、
ｈ）工程ｆで得られた比率が工程ｇにより得られた出現率より有意に高い場合に、工程ｅで選抜した特定位置に特定の塩基種を有する核酸分子が、工程ｇで設定された閾値または有意水準でクローナル核酸混合物中に存在すると推定する工程、
を含む。【００１２】
本開示の第二の態様において、コンピュータを用いて、クローナル核酸混合物を構成する核酸分子の塩基配列における特定位置に存在する特定塩基種がクローナル核酸混合物中で存在する割合を推定する第二の推定方法が提供される。
第二の推定方法は、
ａ）各核酸分子に由来するリード配列を取得する工程、
ｂ）工程ａで取得したリード配列を参照配列上にマッピングしマッピング結果を取得する工程、
ｃ）工程ｂで得られたマッピング結果を基にリード配列のクラスタリング結果を取得する工程、
ｄ）クラスタリング結果を基にコンセンサスリードとコンセンサスリード中の各塩基に対応する信頼性スコアからなるコンセンサスリード情報を取得する工程、
ｅ）クローナル核酸混合物を構成する各核酸分子から工程ａ〜ｄを経て得られるコンセンサスリードの集合に含まれる各コンセンサスリードの配列情報より、参照配列上の特定位置に対応する位置に存在する塩基種とその信頼性スコアを選抜する工程、
ｆ）工程ｅで使用した特定位置を含むコンセンサスリードの総数に対する、工程ｅで選抜した特定位置に特有の各塩基種を含有するコンセンサスリード数の比率を取得する工程、
ｇ）工程ｄで取得した信頼性スコアを用いたパラメトリックなエラーモデルにより、設定された閾値または有意水準での工程ｅで選抜した特定位置において解析エラーによって特有の塩基が出現する出現率を推定する工程、
ｈ）工程ｆで得られた比率が工程ｇにより得られた出現率より有意に高い場合に、工程ｅで選抜した特定位置に特定の塩基種を有する核酸分子が、工程ｇで設定された閾値または有意水準でクローナル核酸混合物中に存在すると推定する工程、及び
ｉ）工程ｈにおいて工程ｅで選抜した塩基種が存在すると判断された場合、工程ｆで得られた特定位置において塩基種を含有するコンセンサスリードの数から、工程ｇで設定された閾値または有意水準における工程ｅで選抜した特定位置における特定塩基種の比率を算出し、クローナル核酸混合物中で存在する割合として採用する工程、
を含む。【００１３】
本開示の第三の態様において、コンピュータを用いて、クローナル核酸混合物を構成する各核酸分子の特定領域の塩基配列を推定する第三の推定方法が提供される。
第三の推定方法は、
ａ）各核酸分子に由来するリード配列を取得する工程、
ｂ）工程ａで取得したリード配列を参照配列上にマッピングしマッピング結果を取得する工程、
ｃ）工程ｂで得られたマッピング結果を基にリード配列のクラスタリング結果を取得する工程、
ｄ）クラスタリング結果を基にコンセンサスリードとコンセンサスリード中の各塩基に対応する信頼性スコアからなるコンセンサスリード情報を取得する工程、
ｅ）クローナル核酸混合物を構成する各核酸分子から工程ａ〜ｄを経て得られるコンセンサスリードの集合に含まれる各コンセンサスリードの配列情報より、参照配列上の特定位置に対応する位置に存在する塩基種とその信頼性スコアを選抜する工程、
ｆ）工程ｅで使用した特定位置を含むコンセンサスリードの総数に対する、工程ｅで選抜した特定位置に特定の各塩基種を含有するコンセンサスリード数の比率を取得する工程、
ｇ）工程ｅで選抜した特定位置の塩基種に対応する工程ｄで取得した信頼性スコアを用いたパラメトリックなエラーモデルにより、設定された閾値または有意水準での工程ｅで確認した特定位置における塩基種の解析エラーによる出現率を推定する工程、
ｈ）工程ｆで得られた比率が工程ｇにより得られた出現率より有意に高い場合に、工程ｅで選抜した特定位置に特定の塩基種を有する核酸分子が、工程ｇで設定された閾値または有意水準でクローナル核酸混合物中に存在する、あるいは該出現率より低い場合は特定位置の塩基種は不明、と判断する工程、
ｉ）工程ｈで得られた結果を該塩基の各コンセンサスリードに適用し二次配列情報を作成する工程、及び
ｊ）工程ｅにおける塩基位置を変更し、工程ｅ〜ｉを繰り返すことにより得られた特定領域の全二次配列情報を統合し、クローナル核酸混合物を構成する各核酸分子の特定領域の塩基配列として採用する工程、
を含む。【００１４】
本開示の第四の態様において、上記第一ないし第三の推定方法のいずれかをコンピュータで実行させるためのコンピュータ読み取り可能なプログラムが提供される。【００１５】
本開示の第五の態様において、クローナル核酸混合物を構成する核酸分子の塩基配列における特定位置の特定塩基種の存在を推定するための推定装置が提供される。
推定装置は、
各核酸分子由来のリード配列を取得するリード配列取得部と、
取得したリード配列をマッピング装置に供給し、前記マッピング装置による参照配列上へのマッピング結果を取得するマッピング情報取得部と、
マッピングしたリード配列のクラスタリング結果を得るクラスタリング結果取得部と、
クラスタリング結果を基にコンセンサスリードとコンセンサスリード中の各塩基に対応する信頼性スコアからなるコンセンサスリード情報取得部と、
を備える。【００１６】
本開示の第六の態様において、クローナル核酸混合物を構成する核酸分子の塩基配列における特定位置の特定塩基種の存在を推定するための推定装置が提供される。
推定装置は、
各核酸分子由来のリード配列を取得するリード配列取得部と、
取得したリード配列をマッピング装置に供給し、前記マッピング装置による参照配列上へのマッピング結果を取得するマッピング情報取得部と、
マッピングしたリード配列のクラスタリング結果を得る第１のクラスタリング結果取得部と、
各クラスタリングに含まれるリード配列を用いて、ローカルアセンブリを実行し、アセンブル配列のクラスタリング結果を得る第２のクラスタリング結果取得部と、
第２のクラスタリング結果取得部のクラスタリング結果を基にコンセンサスリードとコンセンサスリード中の各塩基に対応する信頼性スコアからなるコンセンサスリード情報取得部と、
を備える。
【発明の効果】
【００１７】
本発明によれば、被験体の身体サンプル由来の一部の核酸上に発生する低頻度な変異を検出でき、塩基配列における塩基種を推定できる。
【図面の簡単な説明】
【００１８】
【図１】本開示の実施の形態に係るＤＮＡ配列の変異を検出するシステムの構成を示す図
【図２】ＤＮＡ配列における新規な変異検出方法の手順を示した図
【図３】塩基配列データのファイル形式であるＦＡＳＴＱ形式を説明した図
【図４】変異解析の具体的な手順を示した図
【図５】マッピングデータのファイル形式であるＳＡＭ／ＢＡＭ形式を説明した図
【図６】ＵＭＩ分類作業の概要を説明した図
【図７】ローカルアセンブリに使用するデブルーイングラフ（ｄｅ  Ｂｒｕｉｊｎ  Ｇｒａｐｈ）の作成を説明した図
【図８】変異検出するための統計モデルを説明するための図
【図９】Ｌａｒｇｅ  Ｉｎｄｅｌ変異を検出するための統計モデルを説明するための図
【図１０】実施例１における、Ｐｈｒｅｄクオリティスコアから算出されたエラー率とＵＱから算出されたエラー率の比較を示す図
【図１１】実施例１における、核酸混合物中の８種類の変異検出感度の評価結果を示す図
【図１２】実施例２における、核酸混合物中のリード比率より推定された変異分子の存在比率比較を示す図
【図１３】変異検出装置の内部構成を示したブロック図
【図１４】変異検出装置の制御部により実行される処理を示すフローチャート
【図１５】変異検出装置の制御部により実行される処理を示すフローチャート
【発明を実施するための形態】
【００１９】
以下、添付の図面を参照しながら本開示の実施の形態を説明する。
【００２０】
＜用語の定義＞
本開示における「クローナル核酸混合物」とは、核酸一分子から生物学的な複製を繰り返すことにより発生した複数の核酸の混合物を指す。このようなクローナル核酸が得られる試料としては、特に生物種による限定も無いが、例えばヒトの生体試料が挙げられ、組織、血液、血漿、血清、尿、唾液、粘膜排出物、痰、便および涙が例示される。また組織のＦＦＰＥ（ホルマリン固定パラフィン包埋）試料の他、ｉｎ  ｖｉｔｒｏでの培養による細胞集団であってもよい。ＦＦＰＥ試料由来ＤＮＡは、短く断片化されている、ニックやギャップが入っている、シトシンの脱アミノ化（ウラシル化）が生じているなどダメージを受けていることが多い。そのままでは解析不良の一因となるため、ＮＥＢＮｅｘｔ  ＦＦＰＥ  ＤＮＡ  Ｒｅｐａｉｒ  Ｍｉｘ（ＮＥＢ社製）など市販の試薬を用いて該ＤＮＡを修復することが望ましい。
【００２１】
本開示における「リード配列」とはシーケンサーより出力された塩基配列データである。
【００２２】
本開示における「低頻度に存在する、特定位置の塩基種」とは、特に限定するものではないが、例えば「低頻度な変異」が挙げられる。ここで、「低頻度な変異」とは、生体内で生じた核酸中の塩基の（突然）変異であって、以下の二つの条件ａ、ｂを満たす変異を意図する。
ａ）核酸分子において、１×１０−３／塩基以下の頻度（すなわち、１，０００塩基に一つ以下の頻度）で出現する。
ｂ）核酸分子を含む試料において、特定の位置の塩基に異なる配列の塩基を有する核酸分子の割合が試料中の全核酸分子数の１０％以下となる。
【００２３】
以下、「低頻度に存在する、特定位置の塩基種」として「低頻度な変異」を例に挙げて本開示を説明するが、本開示は低頻度な変異の検出に何ら限定されるものではない。
【００２４】
本開示における（突然）変異は、置換、挿入及び欠失のいずれでも良い。また、該（突然）変異は、既知であってもよく、新規であってもよい。
【００２５】
（実施の形態）
１．変異検出のための構成
図１は、本発明の実施の形態に係る塩基配列の推定方法を実施する変異検出システム（推定システムの一例）の構成を示した図である。変異検出システム１００はＤＮＡ配列における変異を検出する。変異検出システム１００は、断片化されたＤＮＡをシーケンサー５０で解析し、その解析の結果得られたリード配列データを変異検出装置１０（推定装置の一例）に入力する。変異検出装置１０はシーケンサー５０からリード配列データを取得し、解析してＤＮＡ配列における変異を検出する。
【００２６】
２．変異検出方法
図２は、本実施の形態のＤＮＡ配列における新規な変異検出方法の手順を示した図である。図２を参照して、ＤＮＡ配列における変異検出方法の手順を説明する。
【００２７】
まず、生体試料からゲノムＤＮＡやＲＮＡ（クローナル核酸）を抽出し、シーケンサーに適した長さのＤＮＡ（例えば１０００ｂａｓｅ以下）が含まれるように断片化し、断片化クローナル核酸を生成する（Ｓ１）。生体試料が、すでに断片化した状態のＤＮＡ（ｃｅｌｌ−ｆｒｅｅ  ＤＮＡ、ＦＦＰＥサンプルなど）や短いＲＮＡ（ｎｏｎ−ｃｏｄｉｎｇ  ＲＮＡなど）の場合、断片化工程は行わなくともよい。
【００２８】
本発明は特に、低頻度な変異を含む可能性のある核酸の解析に好適に使用される。解析したいクローナル核酸は、ＤＮＡまたはＲＮＡであり、好ましくはＤＮＡであり、さらに好ましくはゲノムＤＮＡである。
【００２９】
次世代シーケンサーの機種によっては解析可能な配列長が数百ｂｐと短いものもある。このため、例えば、クローナル核酸がゲノムＤＮＡの場合、シーケンサーで解析可能な長さに断片化する必要がある。断片化の方法に限定は無く、公知の方法で断片化すればよい。本発明を何ら限定するものではないが、例えば、超音波を用いた物理的な切断を行うＤＮＡ  Ｓｈｅａｒｉｎｇ  Ｓｙｓｔｅｍ  Ｍ２２０／ＭＥ２２０（Ｃｏｖａｒｉｓ社製）、酵素を用いた切断を行うＱＩＡｓｅｑ  ＦＸ  ＤＮＡ  Ｌｉｂｒａｒｙ  Ｋｉｔ（Ｑｉａｇｅｎ社製）など市販の機器や試薬が使用できる。好適には、超音波を用いた物理的な切断が用いられる。
【００３０】
次に各断片化クローナル核酸に「分子バーコード」と呼ばれる核酸を付加する（Ｓ２）。これにより得られたテンプレート核酸はのちにコピーされるが、分子バーコードを付加することにより、コピーがどのクローナル核酸由来かを識別することが可能となる。通常、クローナル核酸を構成する各核酸分子のコピーを作製し、当該コピーのリード配列を得る場合、リード配列の塩基情報のみから本来の核酸分子単位の由来を特定することは不可能である。これに対処するため、該分子バーコードを用いる。本開示において、クローナル核酸に付与される分子バーコードは、一つ以上でユニークであれば良い。また、その配列長や配列パターンは、特に限定されないが、１０００種類のクローナル核酸を識別するには最低５塩基長（４５＝１０２４通り）が必要であり、通常６〜１０塩基長の配列で使用される。なお分子バーコードは様々な呼称があり、“ｕｎｉｑｕｅ  ｍｏｌｅｃｕｌａｒ  ｉｄｅｎｔｉｆｉｅｒ（ＵＭＩ）”や“ｕｎｉｑｕｅ  ｍｏｌｅｃｕｌａｒ  ｔａｇ（ＵＭＴ）”、または単に“分子タグ（ｍｏｌｅｃｕｌａｒ  ｔａｇ）”と呼ばれることもある。クローナル核酸への分子バーコードの付加方法に特に限定は無いが、リガーゼ等の酵素により結合させることで調製することができる。
【００３１】
以下、分子バーコードを付加したクローナル核酸を、テンプレート核酸または解析対象核酸と称することがある。
【００３２】
テンプレート核酸に分子バーコードの他、ＰＣＲ（ｐｏｌｙｍｅｒａｓｅ  ｃｈａｉｎ  ｒｅａｃｔｉｏｎ）増幅用プライマーがアリーニング可能な配列を加えておくと後述の核酸増幅操作が行いやすくなる。このようなテンプレート核酸の調製には、何ら限定するものではないが、例えば、市販製品であるＴｈｒｕＰＬＥＸ（登録商標）Ｔａｇ−Ｓｅｑ  Ｋｉｔ（タカラバイオ社製）を利用することができる。また解析したい核酸領域があらかじめ決まっている場合、当該領域を増幅可能なＰＣＲプライマーと分子バーコードを連結させた核酸を用いることにより、クローナル核酸への分子バーコードの付加と増幅を同時に行うことができる。あるいは、いわゆるプローブキャプチャー法などにより、目的領域を包含するテンプレート核酸を濃縮しても良い。該濃縮操作はコピーを得た後に行っても良い。
【００３３】
該分子バーコートはさらに、「ステム（ｓｔｅｍ）」などと呼ばれる配列を含んでいてもよい。当該ステム配列は、分子バーコードの開始点の目印となる配列である。その配列長や配列パターンは特に限定されないが、例えば、８〜１０塩基長の配列で使用される。当該ステム配列は、クローナル核酸と分子バーコードに挟まれる位置に配置される。
【００３４】
クローナル核酸への分子バーコードの付加後、得られたテンプレート核酸を増幅してライブラリーを作製する（Ｓ３）。本実施の形態では、シーケンサー５０が次世代シーケンサーの場合、配列解析を行うためにはある程度の核酸量を必要とするため、低頻度な変異を含む核酸分子も核酸量を増やす、すなわちコピーを得る必要がある。
【００３５】
テンプレート核酸のコピーの製造法は特に限定しないが、核酸増幅法が挙げられ、操作の簡便さなどからＰＣＲに基づく方法によって行うことが好ましい。このようにして得られたシーケンサーによる解読対象となるサンプルはライブラリーと呼ばれる。
【００３６】
このライブラリーのシーケンシングにより得られた塩基配列データを用いてＤＮＡ配列の変異が解析される。本実施の形態では、テンプレート核酸の増幅は、ＰＣＲに基づく方法によって行うが、これに限定されない。ライブラリーの作製方法は特に限定されない。ライブラリー作製過程において、ＰＣＲ法やプローブキャプチャー法により解析対象とする領域を濃縮する場合もあるが、実施において濃縮手法に関する制限はない。
【００３７】
ライブラリーをシーケンサー５０へ入力し、塩基配列を解析する（Ｓ４）。すなわち、リード配列は、公知のシーケンシング法による塩基配列解読法により取得することができる。シーケンサーのプラットホームは特に限定されないが、解析対象がゲノムＤＮＡのように情報量が膨大になる場合には、次世代シーケンサーが望ましい。次世代シーケンサーとして、例えば、ＨｉＳｅｑ（登録商標）及びＭｉＳｅｑ（登録商標）（ｉｌｌｕｍｉｎａ社製）、Ｉｏｎ  ＰｒｏｔｏｎＴＭおよびＩｏｎ  ＰＧＭＴＭ（Ｔｈｅｒｍｏ  Ｆｉｓｈｅｒ  Ｓｃｉｅｎｔｉｆｉｃ社製）、ＰａｃＢｉｏ（登録商標）  ＲＳ  ＩＩおよびＰａｃＢｉｏ（登録商標）  ＳｅｑｕｅｌＴＭ（Ｐａｃｉｆｉｃ  Ｂｉｏｓｃｉｅｎｃｅｓ社製）、ＭｉｎＩＯＮ、ＧｒｉｄＩＯＮ  Ｘ５、ＰｒｏｍｅｔｈＩＯＮ、およびＳｍｉｄｇＩＯＮ（Ｏｘｆｏｒｄ  Ｎａｎｏｐｏｒｅ  Ｔｅｃｈｎｏｌｏｇｉｅｓ社製）などが挙げられ、数千万から数億のＤＮＡ断片を同時並行的に処理して塩基配列を決定する装置が好ましい。
【００３８】
シーケンサー５０での解析が終了すると、シーケンサー５０からの解析結果を用いてＤＮＡ配列における変異が検出される（Ｓ５）。この工程（Ｓ５）は変異検出装置１０により実施される。
【００３９】
シーケンサー５０から得られる塩基配列データ（リード配列）のファイル形式はＦＡＳＴＱ形式であるが、他の形式でもよい。ＦＡＳＴＱ形式は、当該技術において公知の形式であり、図３（Ａ）に示すような形式となる。データの各列の意味は図３（Ｂ）に示すとおりである。ＦＡＳＴＱ形式のシーケンスデータでは、各塩基のベースコール（塩基の指定）に対する正確さをあらわす指標として、Ｐｈｒｅｄクオリティスコア（Ｐｈｒｅｄ  ｑｕａｌｉｔｙ  ｓｃｏｒｅ）が出力される。本実施の形態では、一塩基あたりのシーケンスによるエラー率をエラー頻度（Ｐｅｒｒｏｒ）で示す。Ｐｈｒｅｄクオリティスコア（ａ）とエラー頻度（Ｐｅｒｒｏｒ）の関係は次式で表される。
Ｐｅｒｒｏｒ  ＝  １０−ａ／１０  （／塩基）
ａ  ＝  ＡＳＣＩＩ−３３
ＡＳＣＩＩコードの値が、ＦＡＳＴＱ形式のデータにおいて出力される。
【００４０】
Ｐｈｒｅｄクオリティスコアは、次世代シーケンサーにより自動的に算出される。エラー頻度への変換式はそれぞれのプラットホームによって異なるが、本実施の形態においては限定される要因はない。
【００４１】
各核酸分子由来のリード配列を取り込む工程において、低品質なリード配列、すなわち、ベースコールの精度（Ｐｈｒｅｄクオリティスコア）が低い塩基を含む配列データ自体を取り除く、いわゆるリードフィルタリング（リードクリーニングともいう）処理を行っても良い。この処理には公知の技術を利用することができ、例えば公知のソフトウエアであるＣｕｔａｄａｐｔ（ドルトムント工科大学）やＴｒｉｍｍｏｍａｔｉｃ（アーヘン工科大学）などが挙げられる。
【００４２】
また、リード配列の一部分の塩基が低品質の場合、該当する一部分の塩基配列のみを取り除く、いわゆるトリミング（ｔｒｉｍｍｉｎｇ）処理を行ってもよい。トリミング処理は、マスキング（ｍａｓｋｉｎｇ）処理と呼ばれることもある。該処理には公知の技術を利用することができ、例えば公知のソフトウエアであるＴｒｉｍｍｏｍａｔｉｃやｓｉｃｋｌｅ（ｈｔｔｐｓ：／／ｇｉｔｈｕｂ．ｃｏｍ／ｎａｊｏｓｈｉ／ｓｉｃｋｌｅ）などが挙げられる。
【００４３】
図４（ａ）（ｂ）は、図２の変異解析処理（Ｓ５）のより具体的な手順を示した図である。後で詳しく説明するように、これらのうち図４（ａ）は、ＰＣＲエラーやシーケンスのエラーを考慮したデータ解析による手順を示したものであり、図４（ｂ）は、リードのアライメントを考慮したデータ解析による手順を示したものである。両方の図から明らかなように、図４（ｂ）は、「ローカルアセンブリ」（Ｓ１７）の工程を含むことのみ、図４（ａ）と異なる。以下、図４（ａ）と図４（ｂ）の両方を参照しつつ、変異解析処理（Ｓ５）の各工程について具体的に説明する。
【００４４】
１）リード配列の取得（Ｓ１１）
まず、シーケンサー５０から、リード配列データを取得する（Ｓ１１）。
【００４５】
２）参照配列へのマッピング（Ｓ１２）  
リード配列はクローナル核酸を構成する各核酸分子由来の様々な断片由来のコピー分子の配列情報の集合体であるため、先ず、参照配列との比較によりリード配列を分類するマッピング処理を行う。ここで、参照配列とは、基準とする配列であり、任意の配列を用いることができる。本発明を何ら限定するものではないが、例えば、ＮＣＢＩ  ＧｅｎＢａｎｋ、ＤＤＢＪ、ＵＣＳＣ  Ｇｅｎｏｍｅ  Ｂｒｏｗｓｅｒなどに登録されている配列を用いることができる。例えば、解析対象がヒトゲノムＤＮＡの場合、Ｇｅｎｏｍｅ  Ｒｅｆｅｒｅｎｃｅ  Ｃｏｎｓｏｒｔｉｕｍ  ｈｕｍａｎ  ｂｕｉｌｄ  ３８（ＧＲＣｈ３８）やＵＣＳＣ  ｈｕｍａｎ  ｇｅｎｏｍｅ  １９（ｈｇ１９）を参照配列とすることができる（これらのバージョンは随時更新されており、必要に応じて適切なバージョンを選択すればよい）。また、参照配列とする領域は、目的に応じて適宜選択することができる。レコードに収録されている配列全長でもよいし、目的とする任意の領域のみであってもよい。該参照配列と後述のコンセンサスリードを比較し、参照配列とは異なる塩基がコンセンサスリードにあった場合、当該塩基は変異の候補となる。
【００４６】
マッピング技術において、ゲノム配列がすでに解読されているヒトを含む生物種については、そのゲノム配列を参照配列として一般にＦＡＳＴＡ形式にて取得可能である。本開示では、テンプレート核酸から解析したリード配列を参照配列にマッピングした後に、比較することで、変異の有無を判定することが可能であり、次世代シーケンスによる解析では、多量のリードシーケンスから変異の有無を検出できる。
【００４７】
本実施の形態において、上記で取得した参照配列に対してリード配列をマッピングすることになる。特に限定するものではないが、たとえば、ＢＷＡ（Ｗｅｌｌｃｏｍｅ  Ｔｒｕｓｔ  Ｓａｎｇｅｒ  Ｉｎｓｔｉｔｕｔｅ）、Ｂｏｗｔｉｅ（ジョンズ・ホプキンズ大学）、ＤＮＡＳＴＡＲ（登録商標）  ＳｅｑＭａｎ  ＮＧｅｎ（ＤＮＡＳＴＡＲ社製）、Ｈｉｓａｔ２（ジョンズ・ホプキンズ大学）、Ｎｏｖｏａｌｉｇｎ（ＮＯＶＡＣＲＡＦＴ社製）等のソフトがマッピング作業に使用できる。
【００４８】
本実施の形態において、マッピングの結果は、図５に示すようなＳＡＭ／ＢＡＭ形式ファイルと呼ばれる公知のファイル形式で出力される。ＢＡＭ形式ファイルはＳＡＭ形式ファイルのバイナリーファイルとなる。マッピングの結果を出力するファイル形式はＢＡＭ形式に限定されるものではない。図５（Ａ）は、参照配列へのマッピング（Ｓ１２）の結果が出力されたＳＡＭ／ＢＡＭ形式ファイルであり、図５（Ｂ）は、分子バーコードのクラスタリング（Ｓ１３）の結果が追加されたＳＡＭ／ＢＡＭ形式ファイルである。ＳＡＭ／ＢＡＭ形式データの各列の意味は図５（Ｃ）に示すとおりである。
【００４９】
このようにマッピングされたリード配列に対し、次に個々のテンプレート核酸に由来するリード配列を分類するクラスタリング処理を行う。テンプレート核酸に前述の分子バーコードを付加した場合は、リード配列中に存在する分子バーコード領域の配列情報を基に分類することが可能である。なお、以下では、分子バーコードを”ＵＭＩ”と称することがある。
【００５０】
ベイズ推定は解析に使用するリード配列の数が多いほど解析精度が高くなるため、得られたリード配列をいかに解析に使用するか、は重要である。
既存の解析手法において、Ｐｈｒｅｄクオリティスコア（Ｐｈｒｅｄ  ｑｕａｌｉｔｙ  ｓｃｏｒｅ）が低いリード配列は解析不良の一因となるため、以降の解析に使用されていない。そのため、大量のシーケンスデータを得ながら、解析に使用できるリード配列が少なくなってしまうという問題があった。
【００５１】
一方、本開示の方法では、Ｐｈｒｅｄクオリティスコア（Ｐｈｒｅｄ  ｑｕａｌｉｔｙ  ｓｃｏｒｅ）が低い箇所の塩基を「Ｎ」（塩基種の判別不可な塩基）とし、参照配列と比較する際には「完全一致した塩基」とみなすことによって以降の解析に使用可能とした。
【００５２】
マッピングによって１つのグループとなったとしても、異なる分子バーコードのリード配列が混在していることもある。その場合、分子バーコード配列に基づいたグループ分けを行うことにより、正しいグループへ再編することができる。
【００５３】
既存の解析手法において、マッピングされなかったリード配列は解析不良の一因となるため、以降の解析に使用されていない。そのため、大量のシーケンスデータを得ながら、解析に使用できるリード配列が少なくなってしまうという問題があった。
【００５４】
しかしながら、本発明において、上記マッピングから外れた（マッピングされなかった、ummapped）リード配列は判断保留とし、分子バーコード配列に基づいたグループ分け工程において、作成されたファミリーと配列比較して解析に使用できるか判断する。このようにして、解析に使用可能となるリード配列を増やしている（解析不可とするリード配列を減らしている）。
【００５５】
３）分子バーコード（ＵＭＩ）のクラスタリング（Ｓ１３）  
本実施の形態では、ＳＡＭ／ＢＡＭ形式ファイルを介して、リード配列のマッピング情報（ｅ．ｇ．リード配列、ＵＭＩ配列、参照配列に対するマッピング位置）を含む一連の情報を取得した後に、各リード配列に付与されたＵＭＩ配列を決定して、同じクローナル核酸に由来するリード配列群を一つに取りまとめる、ＵＭＩの分類工程（分子バーコードのクラスタリング）を行う。図６（Ａ）は、リード配列を、参照配列上でのマッピング開始点と終了点によりグループ分けした後の状態の一例を示している。図６（Ａ）は、参照配列上でのマッピング位置に基づき分類された３つのグループ（Ｐｏｓｉｔｉｏｎ  ｆａｍｉｌｙ）を例示している。図６（Ｂ）は、その後さらに、ＵＭＩ（分子バーコード）の配列によりグループ分けした状態の一例を示している。図６（Ｂ）は、ＵＭＩ（分子バーコード）の配列による分類された２つのグループ（ＵＭＩ  ｆａｍｉｌｙ）を例示している。
【００５６】
ＵＭＩの分類の工程（分子バーコードの配列によりグループ分け）では、以下の工程ａ〜ｄを行う。
【００５７】
工程ａ：参照配列に対するマッピング位置からなる各リード配列の位置に基づき、ＵＭＩ配列に非依存的にリード配列群を分類する。分類されたグループはＰｏｓｉｔｉｏｎ  ｆａｍｉｌｙとして管理する。すなわち、ＳＡＭ／ＢＡＭ形式ファイルに記載されている参照配列上におけるマッピング開始点と終了点が同一のリード配列群をＰｏｓｉｔｉｏｎ  ｆａｍｉｌｙとして管理する。
工程ｂ：各Ｐｏｓｉｔｉｏｎ  ｆａｍｉｌｙ内のリード配列数を集計する。集計値が閾値より少ないＰｏｓｉｔｉｏｎ  ｆａｍｉｌｙにはフラグを立てる。
工程ｃ：各Ｐｏｓｉｔｉｏｎ  ｆａｍｉｌｙ内でＵＭＩ配列の類似度に基づきリード配列を分類する。ＵＭＩ配列の類似度に基づき分類したグループをＵＭＩ  ｆａｍｉｌｙと称する。また、二種類以上のＵＭＩが同一ライブラリー上に存在する場合には、全てのＵＭＩ配列を結合させる。具体的には、図６の例では、Ｌ−ＵＭＩとＲ−ＵＭＩを結合させる。ＵＭＩ配列の類似度に基づきＰｏｓｉｔｉｏｎ  ｆａｍｉｌｙ内のリード配列群をさらに分類する。類似度の閾値を設定して、ある程度のミスマッチを許容することで、閾値内のミスマッチのＵＭＩ同士は同じ配列とする。Ｐｏｓｉｔｉｏｎ  ｆａｍｉｌｙ内のＵＭＩ配列の出現数によって、全リード配列数に対して出現数の少ないＵＭＩ配列を持つリード配列にはＳＡＭ／ＢＡＭファイル上でフラグを立てる。
工程ｄ：Ｐｏｓｉｔｉｏｎ  ｆａｍｉｌｙ内のリード配列数の多寡に基づき、近傍に位置するＰｏｓｉｔｉｏｎ  ｆａｍｉｌｙに限定して、最終ＵＭＩ  ｆａｍｉｌｙの再分類を実施する。すなわち、フラグを立てた、リード数が少ないＰｏｓｉｔｉｏｎ  ｆａｍｉｌｙについて、近傍のＰｏｓｉｔｉｏｎ  ｆａｍｉｌｙに統合できるかをＵＭＩ配列の類似度に基づき再判定する。統合できると判定した場合、フラグを立てた、リード数が少ないＰｏｓｉｔｉｏｎ  ｆａｍｉｌｙを近傍のＰｏｓｉｔｉｏｎ  ｆａｍｉｌｙに統合する。
【００５８】
本実施形態において、ＵＭＩ配列が同一かどうかの配列類似度の判定はハミング距離などの編集距離を使用する。しかし、配列類似度の判定はこれに限定するものではなく、例えば、モチーフ出現頻度に基づいた（類似度）ネットワーク解析を想定したＺｓｃｏｒｅベースのベクトルの類似度や分子系統解析で用いる系統・進化距離といった指標も選択可能である。さらに類似度の判定基準は、ＵＭＩ配列の長さやテンプレート核酸のコピー数によって変更してもよい。
【００５９】
次世代シーケンスの作業では、ＵＭＩ配列自体にもＰＣＲ増幅やシーケンスによるエラーが生じる可能性がある。そのため本実施形態では、このエラー発生の可能性に対して、下記の工程が実施可能である。
工程ａ：シーケンスデータのＰｈｒｅｄクオリティスコアが低い塩基については、Ｎ（任意の塩基を示す記号）でマスクして、配列類似度の比較対象から除外する。
工程ｂ：ＵＭＩ配列間のミスマッチを許容するように編集距離の閾値を調整する。すなわち、ＵＭＩ配列の類似度（ハミング距離等）に基づき、Ｐｏｓｉｔｉｏｎ  ｆａｍｉｌｙ内のリード配列群をさらに分類する。類似度の閾値を設定して、ある程度のミスマッチを許容することで、閾値の条件を満たすミスマッチのＵＭＩ同士は同じ配列とする。
【００６０】
さらに国際公開第２０１６／１７６０９１号では、ＰＣＲ増幅及びシーケンス過程において、最初に付与されたＵＭＩ配列とは異なる配列に置き換わるＵＭＩ置換現象が示唆されている（ｅ．ｇ．  ＵＭＩ  Ｊｕｍｐｉｎｇ、ＰＣＲ  Ｊｕｍｐｉｎｇ  ［ＵＭＩ−ｔｏｏｌｓ：  ｍｏｄｅｌｉｎｇ  ｓｅｑｕｅｎｃｉｎｇ  ｅｒｒｏｒｓ  ｉｎ  Ｕｎｉｑｕｅ  Ｍｏｌｅｃｕｌａｒ  Ｉｄｅｎｔｉｆｉｅｒｓ  ｔｏ  ｉｍｐｒｏｖｅ  ｑｕａｎｔｉｆｉｃａｔｉｏｎ  ａｃｃｕｒａｃｙ，  Ｔｏｍ  Ｓｍｉｔｈ  ｅｔ  ａｌ．  ｈｔｔｐ：／／ｄｘ．ｄｏｉ．ｏｒｇ／  １０．１１０１／ｇｒ．２０９６０１．１１６．］、ｉｎｄｅｘ  ｓｗａｐｐｉｎｇ  ［Ｃｈａｒａｃｔｅｒｉｚａｔｉｏｎ  ａｎｄ  ｒｅｍｅｄｉａｔｉｏｎ  ｏｆ  ｓａｍｐｌｅ  ｉｎｄｅｘ  ｓｗａｐｓ  ｂｙ  ｎｏｎ−ｒｅｄｕｎｄａｎｔ  ｄｕａｌ  ｉｎｄｅｘｉｎｇ  ｏｎ  ｍａｓｓｉｖｅｌｙ  ｐａｒａｌｌｅｌ  ｓｅｑｕｅｎｃｉｎｇ  ｐｌａｔｆｏｒｍｓ，  Ｍａｕｒａ  Ｃｏｓｔｅｌｌｏ  ｅｔ  ａｌ．（ｈｔｔｐ：／／ｄｘ．ｄｏｉ．ｏｒｇ／１０．１１０１／２００７９０）］、Ｃｈａｒａｃｔｅｒｉｚａｔｉｏｎ  ａｎｄ  ｒｅｍｅｄｉａｔｉｏｎ  ｏｆ  ｓａｍｐｌｅ  ｉｎｄｅｘ  ｓｗａｐｓ  ｂｙ  ｎｏｎ−ｒｅｄｕｎｄａｎｔ  ｄｕａｌ  ｉｎｄｅｘｉｎｇ  ｏｎ  ｍａｓｓｉｖｅｌｙ  ｐａｒａｌｌｅｌ  ｓｅｑｕｅｎｃｉｎｇ  ｐｌａｔｆｏｒｍｓ，  Ｍａｕｒａ  Ｃｏｓｔｅｌｌｏ  ｅｔ  ａｌ．（ｈｔｔｐ：／／ｄｘ．ｄｏｉ．ｏｒｇ／１０．１１０１／２００７９０）にて報告／示唆済み）。この可能性に対処するため、ＵＭＩ配列が複数あるようなライブラリー構成の場合には、Ｐｏｓｉｔｉｏｎ  ｆａｍｉｌｙ内の範囲において、少なくとも一つのＵＭＩが完全一致すれば同一のＵＭＩを持つリード配列として分類する工程を含めている。すなわち、複数のＵＭＩが存在するライブラリー構成の場合には、Ｐｏｓｉｔｉｏｎ  ｆａｍｉｌｙ内において、少なくとも一つのＵＭＩが完全に一致すれば同一のＵＭＩを持つリード配列としてまとめる。
【００６１】
本開示を何ら限定するものではないが、例えば、分子バーコードが６塩基からなる場合、一方または両方の分子バーコードにおいて２塩基以下のミスマッチ、または両末端いずれかの分子バーコード領域のリード配列６塩基が一致すればクラスタリングの対象となるリード配列とみなす。これにより、ＵＭＩ配列自体に存在するエラーの可能性を考慮したクラスタリングが可能となる。
【００６２】
ここで、該ＵＭＩ配列自体に存在するエラーとは、ミスマッチ塩基、ベースコール不良塩基、マッピングエラーなどを指す。既存の解析手法ではミスマッチ塩基のみ考慮されていた。一方、本発明ではベースコール不良塩基やマッピングエラーなども考慮することにより、解析精度を向上することができた。
【００６３】
本発明を何ら限定するものではないが、分子バーコード配列自体に存在するエラーの存在率は、分子バーコード塩基長の５０％以下、好ましくは４０％以下、さらに好ましくは３５％以下、特に好ましくは３０％以下であれば、該エラーの可能性を考慮したクラスタリングが可能である。
【００６４】
分子バーコードがクローナル核酸の両末端に付加される場合、該ＵＭＩ配列自体に存在するエラーが一方の分子バーコードのみに存在しても、両方の分子バーコードに存在しても、本発明の方法では、当該エラーの可能性を考慮したクラスタリングが可能である。
【００６５】
本実施の形態において、ＵＭＩの分類結果は前述のようにＳＡＭ／ＢＡＭ形式ファイルのデータに追記する。すなわち、図５（Ｂ）のＳＡＭ／ＢＡＭ形式（出力形式）に示すように、一列目のリード配列の名称に［染色体番号］−［マッピング開始点］−［マッピング終了点］−［アライメントパターン］−［分類後のＵＭＩ配列］−［分類前のＵＭＩ配列］という形式としている。さらにＢＣ：Ｚ：ｓｅｑｕｅｎｃｅ：ｓｅｑｕｅｎｃｅ，ＸＤ：ｉ：ｃｏｕｎｔ，ＸＰ：ｉ：ｓｔａｔｕｓに必要な情報を追加したものをＳＡＭ／ＢＡＭ形式ファイルのオプション情報として出力する（以下の表１参照)。追記の形式はこの形式に限定されるものではない。
【表１】
2019132010.tif 000003
本実施の形態では、リードのアライメントパターンのステータス（ＸＰ：ｓｔａｔｕｓ）情報は、上述のように、０、１、２、３の値を採り、このうちＸＰ＝０の場合が「正常」のステータスを表すものとされている。また、ＸＰ＝１、２、３のときのデータを取り扱うことにより、別染色体にマッピングされるリードの場合や、片方のリードが参照配列上にマッピングされない場合なども、ＵＭＩの分類結果として出力することができる。
【００６６】
本実施の形態において、ＵＭＩ分類結果を追加したＳＡＭ／ＢＡＭ形式ファイルに基づき変異検出を実施する。
【００６７】
４）ローカルアセンブリ（Ｓ１７）
このローカルアセンブリの工程は、図４（ｂ）に示す変異解析処理に含まれる工程である。ローカルアセンブリの工程では、参照配列にはない検体特有のリード配列長以上の配列領域を決定する処理が行われる。この配列と参照配列を比較することで、リード配列長以上の配列サイズを有する塩基置換や挿入や欠損が確認できる。なお本実施の形態のローカルアセンブリは、Ｒｉｍｍｅｒ  ｅｔ  ａｌ． (Ｎａｔｕｒｅ  Ｇｅｎｅｔｉｃｓ，  ２０１４) に掲載されたＰｌａｔｙｐｕｓのアルゴリズムを参考にコーディングされており、参照配列とリード配列よりｔｗｏ−ｃｏｌｏｒｅｄ ｄｅ  Ｂｒｕｉｊｎ  Ｇｒａｐｈを構築してアセンブル配列を構築する。なおローカルアセンブリにて使用可能なオプションおよびその例を、以下の表２に示す。
【００６８】
【表２】
2019132010.tif 000004
【００６９】
ｄｅ  Ｂｒｕｉｊｎ  ｇｒａｐｈでは、まず塩基配列を長さｋの塩基配列（ｋ−ｍｅｒ）に切断して、Ｎｏｄｅ（頂点）とする。さらに各Ｎｏｄｅ（頂点）をＥｄｇｅ（辺）で結ぶ。複数のリード配列もしくは参照配列にてｋ−ｍｅｒのＮｏｄｅ−Ｅｄｇｅを発生させて、同じｋ−ｍｅｒを有するＮｏｄｅや同じＮｏｄｅを合体させたグラフがｄｅ  Ｂｒｕｉｊｎ  ｇｒａｐｈとなる（図７参照）。本実施の形態では、ｄｅ  Ｂｒｕｉｊｎ  ｇｒａｐｈのＮｏｄｅとＥｄｇｅに対して、塩基配列の由来（参照配列もしくはリード配列から由来、分子バーコードの有無およびその配列情報）や信頼性スコア（シーケンスのＰｈｒｅｄクオリティスコア）から算出される重みを付けることができる。
【００７０】
（４−１）ｄｅ  Ｂｒｕｉｊｎ  ｇｒａｐｈ（デブルーイングラフ）の作成
本実施の形態において、ｄｅ  Ｂｒｕｉｊｎ  ｇｒａｐｈは解析遺伝子の一部領域を対象に、配列長(表２のｗｉｎｄｏｗサイズオプションにより変更可能)を指定した対象領域に相当する参照配列と同領域にマッピングされたリード配列を使用して作成する。さらに使用するリード配列には、同領域にマッピングされたリード配列自体とペアリード(つまり、参照配列にマッピングされなかったｕｎｍａｐｐｅｄリードまたは対象領域外にマップされたリード)も含まれる。また、本実施の形態では、グラフのＥｄｇｅ（辺）のｃｏｌｏｒを、参照配列側を０、リード配列側を１として、両方が存在する部分は０．５と設定してＥｄｇｅ（辺）の由来を区別することが可能である。
【００７１】
（４−２）グラフの再構築
本実施の形態において、（例えば、図７の場合、）構築したグラフの５’末端から探索して、環状化していた場合には、ｋ−ｍｅｒの長さを伸長して、再度グラフを作成することができる。
【００７２】
（４−３）最短経路の検出
本実施の形態において、グラフ上の経路探索に、辺の重みが非負数の場合の単一始点最短経路問題を解くための最良優先探索によるアルゴリズムであるダイクストラ法(優先度付きキュー)を採用する (https://ja.wikipedia.org/wiki/ダイクストラ法  または  https://en.wikipedia.org/wiki/Dijkstra%27s_algorithm)。つまり、個々の探索経路がアセンブル配列となり、さらに経路上のＥｄｇｅ(辺)のｃｏｌｏｒの違いにより、参照配列特有もしくはリード配列特有のＥｄｇｅが判別可能であるため、変異候補を同定することができる。しかしながら、経路探索の方法やアルゴリズムは、上記に限定されるものではない。
【００７３】
５）コンセンサスリード（Ｓ１４）
上述の、特に「３）分子バーコード（ＵＭＩ）のクラスタリング（Ｓ１３）」までで述べたように、ＵＭＩ配列に基づき、ＰＣＲエラーやシーケンスのエラーを考慮して、同じＵＭＩ配列（同じＵＭＩ  ｆａｍｉｌｙ）をもつリード配列から共通のリード配列すなわちコンセンサスリードを出力する。当該コンセンサスリードは、単にコンセンサス、コンセンサス配列、またはコンセンサスリード配列と呼ぶこともある。図８（Ａ）は、コンセンサスリードを求めるコンセンサスリード処理の処理工程を説明した図である。（以下、「（ａ）ＰＣＲエラーやシーケンスエラーを考慮したデータ解析」参照）。
一方、ローカルアセンブリを用いた解析では、シーケンスのエラーとリードのマッピング距離を考慮して、上述の、特に「４）ローカルアセンブリ（Ｓ１７）」にて述べたように、同じＵＭＩ配列(同じＵＭＩ  ｆａｍｉｌｙ)をもつリード配列からアセンブル配列を構築し、上述の「（４−３）最短経路の検出」に基づいて、その中から信頼性が高いアセンブル配列をコンセンサスリードとする。図９（Ａ） は、コンセンサスリードを求める処理工程を説明した図である。（後の「（ｂ）リードのアライメントを考慮したデータ解析（Ｌａｒｇｅ  Ｉｎｄｅｌの検出）」参照）。
更に、「（ａ）ＰＣＲエラーやシーケンスエラーを考慮したデータ解析」には、別手法もある（後の「（ａ’）ＰＣＲエラーやシーケンスエラーを考慮したデータ解析の別手法」参照）。
【００７４】
（ａ）ＰＣＲエラーやシーケンスエラーを考慮したデータ解析
本実施の形態では、まず、ＵＭＩ配列に基づき、ＰＣＲエラーやシーケンスエラーを考慮して、同じＵＭＩ配列を持つ（すなわち、同じＵＭＩ  ｆａｍｉｌｙ内の）リード配列から共通のリード配列（コンセンサスリード）を求める。
【００７５】
なお、特定の塩基種が配列解析時のエラーによるものかクローナル核酸中に存在する塩基に由来するものか判断するために、コンセンサスリード及び当該配列中の各塩基に対する信頼性スコア（ＵＱ）を使用する。当該信頼性スコア（ＵＱ）は各ＵＭＩ  ｆａｍｉｌｙ内で各塩基に対して算出され、信頼性スコア（ＵＱ）は次式で表される。
【数１】
2019132010.tif 000005
【００７６】
信頼性スコア（ＵＱ）を算出する際、単純な多数決ではなく、統計処理を行う。コンセンサスリードは、市販製品、例えばＣｏｎｎｏｒ（ｈｔｔｐｓ：／／ｇｉｔｈｕｂ．ｃｏｍ／ｕｍｉｃｈ−ｂｒｃｆ−ｂｉｏｉｎｆ／Ｃｏｎｎｏｒ）を利用して得ることができる。
【００７７】
コンセンサスリード処理の開始前に、以下の工程ａ、ｂを行う。
工程ａ）各ｐｏｓｉｔｉｏｎ（Ｎｏ  ＝  １，２，３…ｋ）にマッピングされたＵＭＩ  ｆａｍｉｌｙ（Ｎｏ  ＝  １，２，３…ｊ）を抽出、
工程ｂ）各ＵＭＩ  ｆａｍｉｌｙ内の全リード配列に対して、番号（Ｎｏ  ＝  １，２，３…ｉ）を割り振り、該当ｐｏｓｉｔｉｏｎ上の塩基種類、各塩基のリード配列数、及びシーケンスのエラー頻度（Ｐｅｒｒｏｒ）のデータを取得する。
【００７８】
これらの処理は、ＵＭＩ分類後のＳＡＭ／ＢＡＭファイルのデータに対して実施する。図８（Ｂ）に示す例の場合、第３染色体の３８，９３０番目の塩基対（ｂｐ）にマッピングされた特定の、すなわちｊ番目のＵＭＩ  ｆａｍｉｌｙ  （Ｎｏ＝ｊ）では、リード配列が合計７本（リード配列数Ｎｏ＝１…７）存在する。第３染色体の３８，９３０番目の位置の塩基がＡとなるリード配列は４本であり、Ｃとなるリード配列は３本である。それぞれのエラー頻度（Ｐｅｒｒｏｒ）は、ＳＡＭ／ＢＡＭ形式ファイルの１１列目のＡＳＣＩＩ情報から入手できる（図５（Ｂ），（Ｃ）参照）。
【００７９】
各ＵＭＩ  ｆａｍｉｌｙ内で各塩基に対する信頼性スコア（ＵＱ）は、特に限定されるものではないが、例えば、ベイズアプローチでは事後分布として推定することができる。具体的には各ＵＭＩ  ｆａｍｉｌｙにおいて、該当する場所での対象塩基（ｅａｃｈ  ｂａｓｅ）毎の確度を推定する。対象塩基は全体集合θの部分集合であり、全４種類の塩基対（Ａ，Ｔ，Ｇ，Ｃ）と、同参照配列上にマッピングされる全ＵＭＩ  ｆａｍｉｌｙで検出されるすべての挿入または欠損塩基（ｏｔｈｅｒ  ｏｂｓ．  ａｌｌｅｌｅｓ）とから構成される。図８（Ｂ）の例では、ＵＭＩ  ｆａｍｉｌｙ（Ｎｏ＝ｊ）の該当箇所（ｃｈｒ３：３８，９３０ｂｐ）において、Ａ，  Ｃの二種類の塩基が検出されているが、集合θの構成はＡ，Ｔ，Ｇ，Ｃの４種類となる（挿入と欠損塩基は検出されないため含まれない）。
【００８０】
公知のソフトウエアＳｍＣｏｕｎｔｅｒの統計モデルは、挿入または欠損塩基に対して検出感度が低い仕様となっている。一方、本発明の統計モデルでは、θにおいて挿入または欠損塩基も考慮することにより、塩基置換に加え、塩基挿入および塩基欠損についても検出感度が高くなった。
【００８１】
信頼性スコア（ＵＱ）の式において、Ｐ（ＵＭＩｊ｜ｅａｃｈ  ｂａｓｅ）は尤度関数であり、Ｐ（ｅａｃｈ  ｂａｓｅ）は対象塩基の出現率を示す事前分布である。事前分布は１／θのような一様分布を採用することができる。しかし、事前データ等からθもしくはＰ（ｅａｃｈ  ｂａｓｅ）に関する何らかの情報があれば変更することが可能である。
【００８２】
尤度関数Ｐ（ＵＭＩｊ｜ｅａｃｈ  ｂａｓｅ）は次式で表される。
【数２】
2019132010.tif 000006
【００８３】
尤度関数Ｐ（ＵＭＩｊ｜ｅａｃｈ  ｂａｓｅ）は、シーケンスエラーとＰＣＲ増幅エラーを考慮した上で、対象塩基の事前分布から各ＵＭＩ  ｆａｍｉｌｙの構成するリード配列データが生じる確率を表現する。例えば、図８（Ｂ）に示す例で検討すると、対象塩基がＡの場合（ｅａｃｈ  ｂａｓｅ  ＝  Ａ）の尤度関数は、ＵＭＩ  ｆａｍｉｌｙでは本来Ａが正しいとした場合に、対象塩基がＡとなるリード配列が４本、Ｃとなるリード配列が３本となる７本のリード配列が生じる尤度Ｐ（ＵＭＩｊ｜Ａ）を計算する。
【００８４】
尤度関数では、シーケンスエラー（Ｓｅｒｒｏｒ）とＰＣＲ増幅エラーを表現している。尤度関数の式において、第一項がシーケンスエラー（Ｓｅｒｒｏｒ）の頻度を示し、第二項がＰＣＲ増幅エラーの発生頻度を示す。図８（Ｂ）の例において、対象塩基がＡの場合（ｅａｃｈ  ｂａｓｅ  ＝  Ａ）、対象塩基がＣとなる３本のリード配列では、シーケンスエラー（Ａ→Ｃへのエラー）がＳｅｒｒｏｒの頻度で生じる。一方、対象塩基がＡとなる４本のリード配列では（１−Ｓｅｒｒｏｒ）の頻度で正しくシーケンスされたことになり、この条件での尤度を求めることになる。
【００８５】
なお、第一項では、ＰＣＲ増幅エラーが生じない場合を想定しているため、ｃｏｒＣｐ（ＰＣＲ増幅エラーが無しとなる確率）を乗算している。さらに、ｎｙは集合θの各構成塩基（Ｙ）のリード配列の数（冗長度）を示している。図８（Ｂ）の例では、Ｙ＝Ａの場合にはｎｙは４本となる。一方、集合θの内で実際に検出できなかったＴ及びＧについては、ｎｙは０本となり、この第一項の計算は実施されないこととなる。すなわち、ＰＣＲ増幅エラーによる尤度（第二項）のみが計算される。
【００８６】
本開示の方法において、シーケンスエラー（Ｓｅｒｒｏｒ）は上述のＰｈｒｅｄクオリティスコアから算出されるＰｅｒｒｏｒに基づき算出している。しかし、挿入・欠損塩基がシーケンスエラーとして生じる可能性は一塩基置換（ｓｉｎｇｌｅ  ｎｕｃｌｅｏｔｉｄｅ  ｖａｒｉａｎｔ、ＳＮＶ）より低くなるという報告もある（Ｍａｒｉｎｉｅｒ  Ｅ等、ＢＭＣ  Ｂｉｏｉｎｆｏｒｍａｔｉｃｓ．  ２０１５  Ｊａｎ  １６；１６：１０．  ｄｏｉ：  １０．１１８６／ｓ１２８５９−０１４−０４３５−６；  Ｓｃｈｉｒｍｅｒ  Ｍ等、Ｎｕｃｌｅｉｃ  Ａｃｉｄｓ  Ｒｅｓ．  ２０１５  Ｍａｒ  ３１；４３（６）：ｅ３７．  ｄｏｉ：  １０．１０９３／ｎａｒ／ｇｋｕ１３４１参照）。そのため本実施の形態では、対象塩基（ｅａｃｈ  ｂａｓｅ）が挿入・欠損塩基の場合には、Ｐｅｒｒｏｒに補正値ｍ（１／１００）を乗算してエラー率を低く抑えるよう処理を加えて、検出性能を向上させている。
【００８７】
尤度関数の第二項は、ＰＣＲ増幅エラーの頻度（ｍｉｓＣｐ）に基づいたコードとなる。ｍｉｓＣｐは、ＰＣＲｅｒｒｏｒ（ｘ）から算出される。
【数３】
2019132010.tif 000007
【００８８】
ＰＣＲｅｒｒｏｒ（ｘ）は、ｘサイクルのＰＣＲサイクルで一つの変異が生じる確率と、ｘサイクルのＰＣＲ産物が最終サイクル（ｎ）後のＰＣＲ産物に含まれる確率との積により求められる。
【数４】
2019132010.tif 000008
前者の変異出現の確率はポアソン分布に従う。ポアソン分布のパラメーターはＰＣＲ酵素のエラー出現率（μ）と、ＰＣＲ酵素の増幅効率（λ）とにより期待値として求めることができる。後者の確率は二項分布の期待値（最大尤度）として算出する。よってＰＣＲｅｒｒｏｒ（ｘ）は下記式で求められる。
【数５】
2019132010.tif 000009
【００８９】
ＰＣＲｅｒｒｏｒ（ｘ）の算出において、変異が出現するＰＣＲサイクル数ｘを推定することが重要であるが、本実施の形態では、この推定をＵＭＩ  ｆａｍｉｌｙ内の各塩基のリード配列数に基づき実施する。図８（Ｂ）の例では、Ａが対象塩基の場合ではＣがＰＣＲ増幅エラーによって出現したケースに該当することから、その割合は３／７（４３．８％）となる。このリード配列数の割合は、全ＰＣＲ産物中に含まれるＰＣＲ増幅エラーに由来するＰＣＲ産物の割合の下記数式で近似的に示される。
Ｓｍ（１＋λｍ）ｎ−ｘ／Ｓｏ（１＋λｗ）ｎ
ここで、Ｓｍは、ＰＣＲ増幅エラーの初発核酸コピー数であり、１である。Ｓｏは、最初の初発核酸コピー数であり、１である。λｍはＰＣＲエラーの核酸産物に対するＰＣＲ増幅効率である。λｗはＰＣＲエラー無しの核酸産物に対するＰＣＲ増幅効率である。
【００９０】
ＰＣＲエラー産物と正常産物のＰＣＲ増幅効率が同じ（λｍ＝λｗ）であると仮定した場合、ＰＣＲサイクル数ｎは既知であることから、エラーが発生したＰＣＲサイクル数ｘを推定可能である。
【数６】
2019132010.tif 000010
【００９１】
尤度関数Ｐ（ＵＭＩｊ｜ｅａｃｈ  ｂａｓｅ）の第二項では、ＰＣＲ酵素のエラー出現率（μ）とＰＣＲ酵素の増幅効率（λ）が重要な設定パラメーターとなる。
【００９２】
本実施の形態では、尤度関数では、ＰＣＲ増幅エラーとシーケンスエラーという二種類のみを考慮し、互いに排反なイベントとして扱っているが、相互作用の効果（ＰＣＲ増幅エラーとシーケンスエラーの両方が生じるケースが該当）やその他の要因も式に追加することができる。その他の要因には、マッピングエラー（異なる箇所と違う場所にリード配列がマッピングされるケース）や酸化反応による塩基置換などが含まれる。
【００９３】
上記の工程を経て、各ＵＭＩ  ｆａｍｉｌｙ内で、各位置の塩基に対して信頼性スコア（ＵＱ）が算出される。図８（Ｂ）のＮｏ．ｊのＵＭＩ  ｆａｍｉｌｙを例にすると、３８，９３０番目の塩基に対して、Ａ，Ｔ，Ｇ，Ｃのそれぞれに対する信頼性スコアが算出されることとなり、その中で最も信頼性スコアが高い塩基をコンセンサスリードの配列とする。図８（Ｂ）の例では、Ａに対する信頼性スコアが最も高くなるため、Ａをコンセンサスリードの配列とする。本手法はコンセンサスリードに対して信頼性スコアを算出する点が従来法と異なる。
【００９４】
（ｂ）リードのアライメントを考慮したデータ解析（Ｌａｒｇｅ  Ｉｎｄｅｌの検出）
本実施の形態では、ローカルアセンブリを用いる解析においては、シーケンスのエラーとリードのマッピング距離を考慮して、同じＵＭＩ配列（同じＵＭＩ  ｆａｍｉｌｙ）を持つリード配列から共通のリード配列、即ちコンセンサスリードを求める。図９（Ａ）は、コンセンサスリードを求める処理工程を示した図である。
【００９５】
つまり、ローカルアセンブリの最短経路(アセンブル配列、ハプロタイプ候補, H)とそれぞれに由来する変異リストから、10base以上のＩｎｓｅｒｔおよび／またはＤｅｌｅｔｉｏｎを持つ経路を対象に下記の手順に準じて、最大事後確率Ｐ{ｅａｃｈ ｈａｐｌｏｔｙｐｅ|ＵＭＩj}を決定し、（後で説明する）変異検出に繋げる。本実施の形態における最大事後確率Ｐは、前述（ａ）の信頼性スコア（ＵＱ）に相当する。
【数７】
2019132010.tif 000011
【００９６】
上述のそれぞれの変数の説明を、下記に示す。図９（Ｂ）を参照されたい。
【数８】
2019132010.tif 000012
［数７］の（２）において、前半部分はリードのマッピング（アライメント）距離の考慮部分を、後半部分はシーケンスのエラーの考慮部分を示す。ここで、前半部分の「Ｐ{distk|μ，σ}：Ｉｎｓｅｒｔサイズ(正規分布μ，σ)に対するリード間の距離の確率」は、同じＵＭＩ配列（同じＵＭＩ Ｆａｍｉｌｙ、ＵＭＩｊ）を有する全リード配列（ｋ∈ＲｅａｄＰａｉｒ）をハプロタイプの配列に再アラインメントして、リード間の距離（ｄｉｓｔｋ）に対する確率を表すものである。なお、本実施の形態における母集団は、全リードの参照配列上のリード間の距離の平均値（μ）と分散（σ２）によって定義される正規分布を、Ｉｎｓｅｒｔサイズの分布(正規分布)とする。
【００９７】
（ｂ−１）リード配列の抽出
本実施の形態に使用されるリード配列は、参照配列に対するマッピング情報によって選抜することが可能である。例えば、本実施の形態では、ＳＡＭ／ＢＡＭ形式のマッピング結果から取得できるリード配列のＣｉｇａｒ（図５（Ｃ）参照）にてＩ／Ｄ／Ｓ／Ｈを持つリード(参照配列にマッピングされないｕｎｍａｐｐｅｄリードも含む)を抽出する。さらに同じＵＭＩを持つリードも抽出する(上記条件に合致しない場合も抽出対象となる)。
【００９８】
（ｂ−２）再アラインメント
特定ＵＭＩ（ＵＭＩｊ）の各リードをハプロタイプ(ｈｉ)に再アラインメントして、リード間の距離とマッチ/ミスマッチ塩基の情報を入手する。本実施の形態の再アラインメントでは、Ｓｍｉｔｈ−Ｗａｔｅｒｍａｎ法を採用している。なお、再アラインメントのアルゴリズムは、Ｓｍｉｔｈ−Ｗａｔｅｒｍａｎ法とは異なるものを採用してもよい。
【００９９】
（ｂ−３）最大事後確率Ｐ{ｅａｃｈ ｈａｐｌｏｔｙｐｅ｜ＵＭＩｊ}の決定 (ＭＡＰ推定によるハプロタイプ推定)
特定ＵＭＩ（ＵＭＩｊ）に対して各ハプロタイプの事後確率を求めて、(全ハプロタイプの中で)最大事後確率を示すハプロタイプ(ｈｉ)から特定ＵＭＩ（ＵＭＩｊ）が生み出されたと仮定する。つまり、Ｐ{ｅａｃｈ ｈａｐｌｏｔｙｐｅ｜ＵＭＩｊ}は特定ＵＭＩ（ＵＭＩｊ）が特定ハプロタイプ(hi)に由来する確率を示す。後続の処理では、最大事後確率に基づき、ポアソン分布によって変異検出を実施する。
【０１００】
（ａ’）ＰＣＲエラーやシーケンスエラーを考慮したデータ解析の別手法
更に、上述の「（ａ）ＰＣＲエラーやシーケンスエラーを考慮したデータ解析」の別手法について説明する。
【０１０１】
（ａ’−１）ＳＮＶ，ＭＮＰ，Ｓｍａｌｌ  Ｉｎｄｅｌｓの検出
上述のＬａｒｇｅ  Ｉｎｄｅｌ（上記「（ｂ）リードのアライメントを考慮したデータ解析」参照）と同じく、ベイズ推定によって実施する。Ｐ{Ａｉ｜ＵＭＩｋ}は、特定ＵＭＩがＡｉというアリルに由来する事後確率を求めることになる（数９）。その際の尤度は、1)ＰＣＲエラーが無くベースコールエラーのみ生じる場合、２)ＰＣＲエラーが生じてベールコールエラーが無い場合、の二種類の条件から構成されている(数１０)。
【数９】
2019132010.tif 000013
【数１０】
2019132010.tif 000014
この数９は、上述の「（ａ）ＰＣＲエラーやシーケンスエラーを考慮したデータ解析」における［数２］と対応するが、ＰＣＲエラーのＣｅの部分が異なる。
【０１０２】
上述のそれぞれの変数の説明は、以下のようになる。
【数１１】
2019132010.tif 000015
例えば、「Ｓｅｒｒｏｒ（ｙ）ｊ」は、数２と同じ定義である。
【０１０３】
（ａ’−２）リードの抽出
該当Ｐｏｓｉｔｉｏｎにマッピングされるリード配列において、該当箇所に変異をもつリードをCigarとNMタグに基づいて選抜する。更に同じＵＭＩのリードも変異の有無にかかわらず全て抽出する。該当箇所の全変異は集合θで表され、Ａ、Ｔ、Ｇ、Ｃと該当箇所の全変異（ｏｂｓ Ａｌｌｅｌｅ）から構成される要素Ａｉが属する。
【０１０４】
（ａ’−３）ＰＣＲエラー(Ｃe)の推定
ある特定のＵＭＩをもつ全１０リードについて、７本が“Ａ”で３本が“Ｔ”となり、“Ｔ”がＰＣＲエラーで生じたと仮定した場合、「Ｃe」は下記のように算出する。
【数１２】
2019132010.tif 000016
【０１０５】
また、それぞれの値の定義と算出方法は以下のものとした。
ＰＣＲエラーを持たないＰＣＲ産物の割合（noFq）：
【数１３】
2019132010.tif 000017
これはポアソン分布を使った割合である。個々の変数の意味は以下の通りである。
n:ＰＣＲサイクル数
λ:ＰＣＲ増幅効率∈[0,1]
μ:ＰＣＲエラー率（ＰＣＲの１サイクルで1塩基に対して生じるエラー頻度）
G: テンプレートの長さ（本実施の形態では１ｂｐである）
【０１０６】
ＰＣＲエラーを持たないＰＣＲ産物数（noPCR）：
【数１４】
2019132010.tif 000018
【０１０７】
ＰＣＲエラーを持つＰＣＲ産物数（ePCR）：
【数１５】
2019132010.tif 000019
【０１０８】
ＰＣＲエラーを持つＰＣＲ産物の割合（eFq）：
【数１６】
2019132010.tif 000020
θ:集合θの要素数 (A,T,C,G,obsAllele)
【０１０９】
（ａ’−４）最大事後確率Ｐ{ｅａｃｈ  ｂａｓｅ｜ＵＭＩｊ}の決定
特定ＵＭＩ（ＵＭＩｊ）に対して各アリルの事後確率を求めて、(全アリルの中で)最大事後確率を示すアリル(ｅａｃｈ  ｂａｓｅ)から特定ＵＭＩ（ＵＭＩｊ）が生み出されたと仮定する。つまり、Ｐ{ｅａｃｈ ｂａｓｅ｜ＵＭＩｊ}は特定ＵＭＩ（ＵＭＩｊ）が特定アリル(ｅａｃｈ ｂａｓｅ)に由来する確率を示す。後続の処理では、最大事後確率に基づき、ポアソン分布によって変異検出を実施する。
【０１１０】
６）エラーモデル構築（Ｓ１５）  
続いて本実施の形態では、全ＵＭＩ  ｆａｍｉｌｙから得られたコンセンサスリード及び各塩基に対する信頼性スコア（ＵＱ）を使用して特定位置の塩基種を推定する（図８（Ａ）、図９（Ａ）参照）。
【０１１１】
該解析を行うために、まずコンセンサスリードの集合より、特定の位置に存在する特定の塩基種を任意に選抜する。ここで「任意」とは「あらゆる場合」、「すべての場合」という意味である。「特定の位置」とは、参照配列とコンセンサスリードの塩基種が異なる位置という意味である。「特定の塩基種」とは、参照配列の塩基と異なるコンセンサスリードの塩基という意味である。
【０１１２】
なお、クローナル核酸において確認したい位置と塩基種については、任意に設定（選択）すればよい。公知の変異塩基の存在を確認したい場合は、その変異に関する部位および塩基種をあらかじめ選択すればよい。また、未知の変異部位を確認したい場合は、参照配列とコンセンサスリードとの比較で参照配列とは異なる塩基種が存在する位置を抽出しても良い。また、塩基種も塩基置換タイプ（変異としてはＳＮＰ）の他、挿入部位に相当する塩基種や欠失相当の塩基種（この場合は該等塩基種無との処理）であっても良い。また、参照配列と異なるコンセンサスリードが存在しない位置については、以降の推定工程を行わなくともよい。
【０１１３】
次いで、塩基種の選抜に用いたコンセンサスリードの総数に対する、特定の位置に特定の各塩基種を含有するコンセンサスリード数の比率を取得する。
【０１１４】
一方、信頼性スコアを用いたパラメトリックなエラーモデルにより、設定された閾値（または有意水準）で特定位置において解析エラーによって特定の塩基が出現する出現率を推定する。
【０１１５】
次いで、上記のコンセンサスリード数を基にした比率が解析エラーによる出現率より有意に高い場合に、上記閾値（または有意水準）で特定塩基位置における特定塩基種がクローナル核酸混合物中に存在すると推定する。
【０１１６】
塩基種推定の基本モデルはポアソン分布によるパラメトリックなエラーモデルを構築して、任意の箇所で塩基種推定を実施する。実施形態によっては、エラーモデルを変更することが可能である。ポアソン分布の他、例えば、多値変数（事前分布はディリクレ分布）や二値変数でのベータ二項分布での検定、教師付きの学習モデルなども、本発明のパラメトリックなエラーモデルとして使用できる。
【０１１７】
塩基種推定は、以下の工程ａ〜ｄを含む。
工程ａ）対象となる箇所を選択する（図８（Ｂ）ではｃｈｒ３：３８，９３０）。
工程ｂ）該当する箇所にマッピングされた全ＵＭＩ  ｆａｍｉｌｙを抽出する（図８（Ｂ）では、ｊ個のＵＭＩ  ｆａｍｉｌｙを抽出）。
工程ｃ）各ＵＭＩ  ｆａｍｉｌｙの該当箇所におけるコンセンサスとＵＱのスコアを取得する。
工程ｄ）確認したい塩基種を有する参照配列（Ａ）とは異なるコンセンサスリードとＵＱのスコアのみを選択する。
工程ｅ）工程ａ〜ｄを経て得られた情報（選択した塩基種を有するコンセンサスリードのパターンとＵＱスコア）からエラーモデル（ポアソン分布）を構築する。
【０１１８】
ポアソン分布のパラメーターλについては、好ましくは以下のように表現できる。
【０１１９】
ある位置において、ｍ個のＵＭＩ  ｆａｍｉｌｙが得られ、そのうちｎ個のＵＭＩ  ｆａｍｉｌｙで塩基Ｘ（＝Ａ，Ｔ，Ｇ  ｏｒ  Ｃ）のコンセンサスリードが得られた場合、下記式によりポアソン分布のパラメーターλの期待値Ｅ［λ］が算出できる。
【数１７】
2019132010.tif 000021
ここで、ａｒｇｍｉｎＰ｛Ｘ│ＵＭＩ｝は、塩基Ｘとなったｎ個のＵＭＩ  ｆａｍｉｌｙのＵＱ（Ｘ）の中の最小値を示す。
【０１２０】
例えば、ｃｈｒ３：３８，９３０の箇所で１０００個のＵＭＩ  ｆａｍｉｌｙが得られ、そのうち５個のＵＭＩ  ｆａｍｉｌｙでＣのコンセンサスリードが得られた場合、下記式によりポアソン分布のパラメーターλが算出できる。
λ＝１０００＊（１−ａｒｇｍｉｎＰ｛Ｃ｜ＵＭＩ｝）
ここで、ａｒｇｍｉｎＰ｛Ｃ｜ＵＭＩ｝は、塩基がＣとなった５個のコンセンサスリードに対する５個のＵＱスコアＵＱ（ｃ）（Ｐ｛Ｃ｜ＵＭＩ｝）の中の最小値を示す。
【０１２１】
上記のようにしてλを求めることで、ポアソン分布を用いたエラーモデルを構築できる。
【０１２２】
７）変異検出（Ｓ１６）
上記のようにして構築したエラーモデル（ポアソン分布）を利用して変異を検出する。すなわち、エラーモデルから、ｍ個のＵＭＩを取得した場合に、エラーによってｎ個のコンセンサスリードが得られる確率が推定できる。
【０１２３】
例えば、１，０００個のＵＭＩを取得した場合に、エラーによって５個のコンセンサスリードが得られる確率が推定できる。ここで、前述の閾値は有意水準とも呼ばれ、一般的には０．０５（５％）に設定されるが、必要に応じて０．０１（１％）や０．００１（０．１％）に設定される。そして、エラーモデルから推定した確率を所定の閾値（１％水準又は５％水準等）と比較する。推定した確率が所定の閾値より低い場合には、エラーによる仮説を棄却することになるため、選択した塩基種であると決定する。これにより、１０００分子の中の１分子のみ存在する塩基種の存在の推定が可能となる。したがって、本開示は低頻度な変異の検出に好適に使用できる。
【０１２４】
実施形態によっては、エラーモデルは適宜変更することが可能である。
【０１２５】
前述のとおり推定された特定塩基位置における特定塩基種がクローナル核酸混合物中に存在する比率は、下記計算式によって算出することができる。また、該比率は、変異頻度とも呼ばれる。
比率（変異頻度）＝（変異箇所にマッピングされた変異を有するコンセンサスリードのＵＭＩ数）／（変異箇所にマッピングされたコンセンサスリードの全ＵＭＩ数）
【０１２６】
算出された該比率は、クローナル核酸混合物中に存在する低頻度な変異の割合（存在率）と見なすことができる。
【０１２７】
各ＵＭＩ  ｆａｍｉｌｙ内において、解析対象とする特定領域の全塩基に対して、特定塩基位置における特定塩基種がクローナル核酸混合物中に存在すると推定できるか、あるいは特定塩基位置における塩基種が不明であると推定できるか、検証し、当該特定塩基種が存在する情報を集約することにより、当該集約された配列情報はクローナル核酸混合物中に存在する低頻度な変異体の配列とすることができる。即ち、特定塩基位置における特定塩基種がクローナル核酸混合物中に存在すると推定できるか、あるいは特定塩基位置における塩基種が不明であると推定できるか、検証することにより得られる結果を、当該塩基のコンセンサスリードに適用して二次配列情報を作成し、更に塩基位置を変更しつつ上述の処理を繰り返すことにより得られる特定領域の全二次配列情報を統合して、クローナル核酸混合物を構成する各拡散分子の特定領域の塩基配列として採用することができる。なお、参照配列と異なるコンセンサスリードが存在しない（参照配列とすべてのコンセンサスリードが同じ）位置については、特定塩基種は存在しないと判断し、クローナル核酸の当該位置の塩基は参照配列と同一であるとする。
【０１２８】
必要に応じて、さらに変異フィルタリングを行うことができる。例えば、Ｐ．  Ｆｌａｈｅｒｔｙ等、Ｎｕｃｌｅｉｃ  Ａｃｉｄｓ  Ｒｅｓ．、第４０巻、第１号、ｅ２ページ、２０１２年（ＤＯＩ：  １０．１０９３／ｎａｒ／ｇｋｒ８６１）に記載の方法を利用し、正常検体で検出された変異を腫瘍検体で検出された変異から除外してもよい。これにより、腫瘍検体に特異的な変異を抽出（フィルタリング）することができる。
【０１２９】
３．実施例
以下、本実施の形態の方法を適用した実施例を説明する。本開示の適用は以下の実施例に限定されるものではない。
【０１３０】
（１）実施例１：断片化された核酸混合物中の低頻度な変異の検出
核酸混合物中に存在する稀な変異の検出感度を調査するために、ＥＧＦＲ、ＫＲＡＳ、ＮＲＡＳ、ＰＩＫ３ＣＡの４種のヒト遺伝子群において、野生型遺伝子のコピー、さらに同じ遺伝子に稀な変異が生じている８種類の変異型遺伝子のコピーを異なる混合比で有するヒト由来のクローナル核酸として、ＭｕｌｔｉｐｌｅｘＩ  ｃｆＤＮＡ  Ｒｅｆｅｒｅｎｃｅ  Ｓｔａｎｄａｒｄ  Ｓｅｔ（ｈｏｒｉｚｏｎ社製）を用いた。本実施例では混合比として、変異型遺伝子のコピー数が１％（Ｐａｒｔ  Ｎｏ．：  ＨＤ７７８）、０．１％（Ｐａｒｔ  Ｎｏ．：  ＨＤ７７９）の二種類のクローナル核酸を使用した。
【０１３１】
ＭｕｌｔｉｐｌｅｘＩ  ｃｆＤＮＡ  Ｒｅｆｅｒｅｎｃｅ  Ｓｔａｎｄａｒｄ  Ｓｅｔは、血漿中の無細胞核酸（ｃｅｌｌ−ｆｒｅｅ  ＤＮＡ）を模して、人工的に１６０ｂｐ程度に断片化したクローナル核酸の標準品である。
【０１３２】
それぞれのクローナル核酸から５０ｎｇを用いて、後述の操作を実施した。なお、ヒトゲノムＤＮＡは、約３０億塩基対から構成され、１塩基の平均分子量は３３０となるため、１０ｎｇのクローナル核酸には約３０００コピーのゲノムＤＮＡが存在する。変異型遺伝子の混合比が１％の場合には、３０コピーの変異型遺伝子が存在し、混合比が０．１％の場合には３コピーの変異遺伝子が存在することになる。
【０１３３】
ＴｈｒｕＰＬＥＸ  Ｔａｇ−ｓｅｑ  ６Ｓ  （１２）  Ｋｉｔ（タカラバイオ社製）のプロトコールを使用して、５０ｎｇのクローナル核酸の両末端にＳｔｅｍ−Ｌｏｏｐアダプターをライゲーション法により付加し、テンプレート核酸とした。次いで、得られたテンプレート核酸を鋳型とし、アダプター配列を有するユニバーサルプライマーによりＰＣＲ増幅を実施してイルミナ社のシーケンサー用の増幅核酸物（アダプター＋テンプレート核酸）を調製した。ＴｈｒｕＰＬＥＸ技術では、６塩基から成るランダムな配列のＵＭＩがＳｔｅｍ−Ｌｏｏｐアダプター中に含まれているため、増幅前に個々のテンプレート核酸にＵＭＩを付加することになる（合計２個のＵＭＩが含まれる）。
【０１３４】
対象領域となる８個の遺伝子領域に由来するテンプレート核酸のみをシーケンスするために、カスタムのＳｕｒｅＳｅｌｅｃｔのキャプチャプローブ（アジレント社製）を用いて、増幅核酸物から８個の遺伝子領域を含む特定ターゲット領域を濃縮した。具体的な操作は、ＴｈｒｕＰＬＥＸ  Ｔａｇ−ｓｅｑ  Ｋｉｔにて推奨されているマニュアルに従って行った。
【０１３５】
イルミナ社の次世代シーケンサーＨｉＳｅｑを使用する超並列シーケンシングアプローチを用いて、各検体に対する両末端１００ｂｐのシーケンスを実施して塩基配列データ（リード配列データ）を取得した。
【０１３６】
取得したリード配列データに対して、Ｔｒｉｍｍｏｍａｔｉｃ（ｖｅｒｓｉｏｎ：０．３２）を使用して低品質なリード配列データを除去した。
【０１３７】
リード配列データを、ＢＷＡ−ＭＥＭ（ｖｅｒｓｉｏｎ：０．７．１５）を用いてヒトの参照配列（ｈｇ１９）にマッピングを実施した。得られたＳＡＭ／ＢＡＭ形式ファイルより、各リード配列について参照配列に対するマッピング情報とＵＭＩの配列情報を入手した。
【０１３８】
上記の手法に従って、マッピングの位置情報とＵＭＩの配列情報に基づいて、特定ターゲット領域に対して分子バーコードのクラスタリングを実施した。分類の条件は下記のとおりである。
条件１：  ミスマッチ２塩基以内のものは同一ＵＭＩとする。
条件２：  条件１を満たさない場合でも、少なくとも一つのＵＭＩ配列が同一であれば、同じＵＭＩとする。
【０１３９】
上記の手法に従って、コンセンサスリードの作成を実施した。その後、第３染色体に座乗するＰＩＫ３ＣＡ遺伝子のエキソン１０に該当する１７８，９３５，９９７−１７８，９３６，１２１ｂｐの各塩基について、１）イルミナ社シーケンサーのリード配列に付与されるＰｈｒｅｄクオリティスコアから算出されるエラー率Ｐｅｒｒｏｒと、２）ＵＱから算出されたエラー率１−ＵＱ（ｂａｓｅ）を算出して、二種類のエラー率を比較した。その結果を図１０に示す。図１０では、分子バーコードのクラスタリング後のＳＡＭ／ＢＡＭファイルを用いて、各箇所でＵＭＩ数が最も多い塩基を対象にして上記二種類のエラー率を算出し、Ｐｈｒｅｄスケール（−１０ｌｏｇ１０（エラー率））でプロットした。
【０１４０】
図１０に示すように、Ｐｈｒｅｄクオリティスコアから算出されたエラー率（図１０中の×）は、Ｐｈｒｅｄスケールで３０〜４０の間となっているが、ＵＱに基づいたエラー率（図１０中の◆）は５０を超える値を示しており、全領域でエラー率が約１／１０倍に抑えられた。
【０１４１】
本開示の方法にて、８種類の変異型遺伝子が検出できるかどうかを検証した。さらに従来法との比較として、１）  ＬｏＦｒｅｑ（ｖｅｒｓｉｏｎ：  ２．１．２  ｈｔｔｐ：／／ｃｓｂ５．ｇｉｔｈｕｂ．ｉｏ／ｌｏｆｒｅｑ／）と、２）  ＳｍＣｏｕｎｔｅｒ  （ｖｅｒｓｉｏｎ：  ｈｔｔｐｓ：／／ｇｉｔｈｕｂ．ｃｏｍ／ｘｕｃｈａｎｇ１１６／ｓｍＣｏｕｎｔｅｒ  １７０８２３取得）とを用いて変異検出を実施した。その結果を図１１に示す。図１１において、グレーのハッチング部の変異は、検出された変異を示し、ＲＥＦは参照配列のリード配列数を示し、ＡＬＴは変異配列を有するリード配列数を示し、ＡＦはリード配列比から算出されるアリル頻度を示す。
【０１４２】
ここで、ＬｏＦｒｅｑは、論文（Ｗｉｌｍ  Ａ等、Ｎｕｃｌｅｉｃ  Ａｃｉｄｓ  Ｒｅｓ．  ２０１２  ４０（２２）：  １１１８９−１１２０１）で報告された変異検出ソフトであり、Ｐｈｒｅｄ  ｑｕａｌｉｔｙ  ｓｃｏｒｅから算出されたシーケンスエラー率よりベルヌーイ試行を実施して変異を検出するソフトである。図１１では、ＢＷＡでマッピングしたＳＡＭ／ＢＡＭファイルに対して、Ｃｏｎｎｏｒ（ｖｅｒｓｉｏｎ：  ０．５．１  ｈｔｔｐｓ：／／ｇｉｔｈｕｂ．ｃｏｍ／ｕｍｉｃｈ−ｂｒｃｆ−ｂｉｏｉｎｆ／Ｃｏｎｎｏｒ）を用いてコンセンサスリードを取得した後に変異検出を実施した。
【０１４３】
ここで、ＳｍＣｏｕｎｔｅｒは、ベイズアプローチにより分子バーコードに基づいてエラー補正を実施して変異を検出するソフトである（非特許文献３）。図１１では、分子バーコードのクラスタリング（本開示手法）を実施したＳＡＭ／ＢＡＭ形式ファイルを入力ファイルとして、ＳｍＣｏｕｎｔｅｒを用いて変異検出を行った結果を示した。
【０１４４】
図１１に示すとおり、ＬｏＦｒｅｑでは混合比１．０％の８種の変異中、ハッチングで表示した、ＥＧＦＲ遺伝子のΔ７４６−７５０の欠損変異（Ｄｅｌｅｔｉｏｎ）、Ｔ７９０Ｍの一塩基置換（ＳＮＶ）、ＮＲＡＳ遺伝子のＱ６１ＫおよびＡ５９ＴのＳＮＶの４種類が検出され、検出感度は５０％（４種類／８種類）であった。一方、ＳｍＣｏｕｎｔｅｒでは混合比１．０％の８種の変異中、ハッチングで表示した、ＥＧＦＲ遺伝子のΔ７４６−７５０の欠損変異とＶ７６９−Ｄ７７０の挿入変異以外の６種類の変異が検出され、検出感度は７５％（６種類／８種類）であった。しかしながら、混合比０．１％の変異にて調査した際には、両ソフトとも８種類のいずれの変異遺伝子も検出できなかった。
【０１４５】
図１１に記載の通り、本発明の解析コード（図中：「独自パイプライン」）は、ハッチングで表示した、混合比１．０％の８種類の変異型遺伝子を全て検出し、検出感度は１００％（８種類／８種類）であった。さらに混合比０．１％の変異中、ハッチングで表示した、ＥＧＦＲ遺伝子のＬ８５８ＲのＳＮＶおよびΔ７４６−７５０のＤｅｌｅｔｉｏｎ、ＮＲＡＳ遺伝子のＧ１２ＤおよびＡ５９ＴのＳＮＶ、ならびにＰＩＫ３ＣＡのＥ５４５ＫのＳＮＶの５種類の変異を検出し、検出感度は６２．５％（５種類／８種類）であった。検出できた該変異のｐ値はいずれもＥ−０６  （１０−６）オーダーであり、閾値（１％水準）より低かった。
【０１４６】
以上、実施例１の結果より、本実施形態の検出手法は従来法より優れた検出感度を示し、さらには従来法では検出できなかった極めて低い変異においても、変異を検出することができた。
【０１４７】
（２）実施例２：核酸混合物中における変異頻度の定量
核酸混合物中に存在する稀な変異の混合比率（変異頻度）を調査するため、１９種の遺伝子群において、野生型遺伝子のコピー、さらに同じ遺伝子に稀な変異が生じている３５種類の変異型遺伝子のコピーを１．０％−３０％の様々な混合比で有するＴｒｕ−Ｑ７（１．３％  Ｔｉｅｒ）Ｒｅｆｅｒｅｎｃｅ  Ｓｔａｎｄａｒｄ（Ｐａｒｔ  Ｎｏ．：  ＨＤ７３４）（ｈｏｒｉｚｏｎ社製）をテンプレート核酸に用いた。
【０１４８】
実施例１と同じく、ＴｈｒｕＰＬＥＸ  Ｔａｇ−ｓｅｑ  ６Ｓ  （１２）  Ｋｉｔを用いてイルミナ社シーケンサー用のライブラリーを作製して、ＨｉＳｅｑによりシーケンスリードを取得した。
【０１４９】
実施例１と同様に変異検出を実施した結果、３５種類の変異中３４種類の変異を検出できた（検出感度  ９７．１％）。結果を図１２に示した。これらの検出された変異の頻度を分子バーコードのクラスタリング後のＳＡＭ／ＢＡＭファイルより算出した。算出方法は、下記のとおりである。
【０１５０】
変異頻度は以下の式によって算出される。
変異頻度＝（変異箇所にマッピングされた変異を有するコンセンサスリードのＵＭＩ数）／（変異箇所にマッピングされたコンセンサスリードの全ＵＭＩ数）
【０１５１】
上記計算式により算出された変異頻度（Ｏｂｓｅｒｖｅｄ  Ａｌｌｅｌｅ  Ｆｒｅｑｕｅｎｃｙ）とＤｒｏｐｌｅｔ  Ｄｉｇｉｔａｌ  ＰＣＲＴＭ（Ｂｉｏ−Ｒａｄ社製）によって求められた変異頻度（Ｅｘｐｅｃｔｅｄ  Ａｌｌｅｌｅ  Ｆｒｅｑｕｅｎｃｙ、製品情報に明記された値）を比較して算出した。結果を図１２に示した。その結果、本開示の方法により推定された頻度は、Ｄｉｇｉｔａｌ  ＰＣＲで求められた値と極めて近く相関係数Ｒ２は０．９９４５であった。
【０１５２】
以上、実施例２の結果より、本開示の方法での分子バーコードのクラスタリングにより精度よくＵＭＩ分類が実現できた。また、分解が進んだ（断片化された）低品質のゲノムＤＮＡであっても低頻度な変異を高精度に検出することができた。
【０１５３】
４．変異検出装置
図１３は、上述した変異検出方法を実施する変異検出装置１０の具体的な内部構成を示したブロック図である。変異検出装置１０はパーソナルコンピュータのような情報処理装置で構成される。変異検出装置１０は、その全体動作を制御する制御部１１と、画面表示を行う表示部１３と、ユーザが操作を行う操作部１５と、データやプログラムを記憶するデータ記憶部１７と、外部機器やネットワークと接続し、通信を行なうインタフェース部１８と、制御部１１が処理を行うためのメインメモリであるＲＡＭ（Ｒａｎｄｏｍ  Ａｃｃｅｓｓ  Ｍｅｍｏｒｙ）１６とを備える。
【０１５４】
表示部１３は例えば、液晶ディスプレイや有機ＥＬディスプレイで構成される。操作部１５はキーボード、マウス、タッチパネル等を含む。
【０１５５】
インタフェース部１８は、ＵＳＢ、ＨＤＭＩ（登録商標）、ＳＣＳＩ等のインタフェースに準拠した種々の機器（プリンタ、通信装置、入力装置等）を接続可能であり、接続した機器と変異検出装置１０間のデータや制御コマンドの通信を可能とする。特に、変異検出装置１０はインタフェース部１８を介して次世代シーケンサー５０からリード配列データを取得する。
【０１５６】
制御部１１は変異検出装置１０全体の動作を制御するものであり、プログラムを実行することで所定の機能を実現するＣＰＵやＭＰＵで構成される。制御部１１で実行されるプログラムは通信回線や、ＣＤ、ＤＶＤ、メモリカード等の記録媒体を介して提供されることができる。または、制御部１１は、所定の機能を実現するように設計された専用のハードウェア回路（ＦＰＧＡ，ＡＳＩＣ等）で構成されてもよい。
【０１５７】
該プログラムは、例えば、ＪＡＶＡ（登録商標）、Ｃ、Ｃ＋＋、Ｃ＃、ＯＢＪＥＣＴＩＶＥ−Ｃ、ＳＷＩＦＴなどのコンピュータ言語、またはＰＥＲＬ、ＰＹＴＨＯＮ、ＢＩＯ−ＲＵＢＹなどのスクリプト言語など任意の適切な言語を用いて、例えば、従来型又はオブジェクト指向技術を用いて、ソフトウェアコードとして実行することができる。
【０１５８】
該ソフトウェアコードは、ＲＡＭ１６など記憶及び／又は送信のためのコンピュータ読み取り可能な媒体において、一連の命令として記憶することができる。あるいは、コード化され、インターネットを含む種々のプロトコールに適合する有線、光、及び／又は無線通信回線を介して送信するために適合されたキャリア信号を用いて送信することができる。
【０１５９】
データ記憶部１７はデータやプログラムを記憶する装置であり、例えばハードディスク（ＨＤＤ）、ＳＳＤ（ＳＯＬＩＤ  ＳＴＡＴＥ  ＤＲＩＶＥ）、半導体メモリ素子、光ディスクで構成することができる。データ記憶部２１は、制御部１１で実行される制御プログラムや、リード配列データ等を格納する。
【０１６０】
図１４（ａ）（ｂ）は、図１３の変異検出装置１０の制御部１１により実行される処理を示すフローチャートである。このうち図１４（ａ）のフローチャートは、ＰＣＲエラーやシーケンスのエラーを考慮したデータ解析に係る処理を示している。図１４（ｂ）のフローチャートは、リードのアライメントを考慮したデータ解析に係る処理を示している。
【０１６１】
まず、図１４（ａ）に示すＰＣＲエラーやシーケンスエラーを考慮したデータ解析に係る処理を説明する。
【０１６２】
変異検出装置１０はインタフェース部１８を介してシーケンサー５０からリード配列データを取得する（Ｓ５１）。変異検出装置１０の制御部１１は、各リード配列データについて参照配列へのマッピング処理を行う（Ｓ５２）。制御部１１は、上述の方法にしたがい参照配列上の位置及びＵＭＩに基づきリード配列データの分類を行う（Ｓ５３）。その後、制御部１１は、コンセンサスリード処理を行い、各ＵＭＩ  ｆａｍｉｌｙについてコンセンサスリードを作製する（Ｓ５４）。
【０１６３】
制御部１１はエラーモデルを構築する（Ｓ５５）。具体的には、制御部１１は、信頼性スコアＵＱを計算し、コンセンサスリードのパターンと信頼性スコアＵＱからエラーモデル（ポアソン分布）のパラメーターλを求めて、エラーモデルを構築する。
【０１６４】
そして、制御部１１は、エラーモデルを参照して変異を検出する（Ｓ５６）。具体的には、制御部１１は、変異を確認したい位置と塩基種を選択し、エラーモデルを参照してその変異の発生確率を推定する。そして、制御部１１は、推定値が閾値より小さいときは、変異が発生していないと判断し、選択した塩基種が正しいと判断する。図１４（ａ）に示すフローチャートでは、小さいサイズの変異が検出され得る。
【０１６５】
任意の工程として、必要であれば、制御部１１は変異フィルタリングを行い（Ｓ５７）、所定の条件を満たす変異を変異検出結果から除外する。
【０１６６】
次に、図１４（ｂ）に示す、リードのアライメントを考慮したデータ解析に係る処理を説明する。変異検出装置１０はインタフェース部１８を介してシーケンサー５０からリード配列データを取得する（Ｓ５１）。変異検出装置１０の制御部１１は、各リード配列データについて参照配列へのマッピング処理を行う（Ｓ５２）。制御部１１は、上述の方法にしたがい参照配列上の位置及びＵＭＩに基づきリード配列データの分類を行う（Ｓ５３）。
【０１６７】
図１４（ｂ）に示すフローチャート（処理）では、「分子バーコードのクラスタリング（Ｓ５３）」と（後続の）「コンセンサスリード処理（Ｓ５４）」との間に、制御部１１は、参照配列にはない検体特有のリード配列長以上の配列領域を決定するローカルアセンブリ処理（Ｓ５８）を行う。
【０１６８】
その後、制御部１１は、コンセンサスリード処理を行い、各ＵＭＩ  ｆａｍｉｌｙについてコンセンサスリードを作製する（Ｓ５４）。次に、制御部１１はエラーモデルを構築する（Ｓ５５）。具体的には、制御部１１は、信頼性スコアＵＱを計算し、コンセンサスリードのパターンと信頼性スコアＵＱからエラーモデル（ポアソン分布）のパラメーターλを求めて、エラーモデルを構築する。更に、制御部１１は、エラーモデルを参照して変異を検出する（Ｓ５６）。具体的には、制御部１１は、変異を確認したい位置と塩基種を選択し、エラーモデルを参照してその変異の発生確率を推定する。そして、制御部１１は、推定値が閾値より小さいときは、変異が発生していないと判断し、選択した塩基種が正しいと判断する。更に必要であれば、制御部１１は変異フィルタリングを行い（Ｓ５７）、所定の条件を満たす変異を変異検出結果から除外する。
【０１６９】
図１４（ｂ）に示すフローチャート（処理）では、「分子バーコードのクラスタリング（Ｓ５３）」と（後続の）「コンセンサスリード処理（Ｓ５４）」との間に、「ローカルアセンブリ処理（Ｓ５８）」が存在すること以外は、図１４（ａ）の処理と同様である。但し、図１４（ｂ）に示すフローチャートでは、大きいサイズの変異が検出され得る。
【０１７０】
なお、図１５のフローチャートに示すように、
・リード配列処理（Ｓ５１）、
・参照配列のマッピング（Ｓ５２）、及び、
・分子バーコードのクラスタリング（ＵＭＩの分類）（Ｓ５３）
の処理以降に、
・コンセンサスリード処理（Ｓ５４ａ）、
・エラーモデル構築（Ｓ５５ａ）、及び、
・変異検出（小さい変異）（Ｓ５６ａ）
の、小さい変異の検出のための処理フローと、
・ローカルアセンブリ処理（Ｓ５８）
・コンセンサスリード処理（Ｓ５４ｂ）、
・エラーモデル構築（Ｓ５５ｂ）、及び、
・変異検出（大きい変異）（Ｓ５６ｂ）
の、大きい変異の検出のための処理フローとを並行して行い、その後、「小さい変異」の検出と「大きい変異」の検出との和集合の出力データを作成する（「和集合処理（Ｓ５６ｃ）」、という処理の流れであってもよい。
【０１７１】
以上のようにして変異検出装置１０は、シーケンサー５０からのリード配列のデータに基づきＤＮＡ配列の変異を検出することができる。
【０１７２】
（本開示）
上記の実施の形態は以下の技術思想を開示している。
【０１７３】
（１）コンピュータを用いて、クローナル核酸混合物を構成する核酸分子の塩基配列における特定位置に存在する特定塩基種の存在を推定する方法であって、
ａ）各核酸分子に由来するリード配列を取得する工程、
ｂ）工程ａで取得したリード配列を参照配列上にマッピングしマッピング結果を取得する工程、
ｃ）工程ｂで得られたマッピング結果を基にリード配列のクラスタリング結果を取得する工程、
ｄ）クラスタリング結果を基にコンセンサスリードとコンセンサスリード中の各塩基に対応する信頼性スコアからなるコンセンサスリード情報を取得する工程、
ｅ）クローナル核酸混合物を構成する各核酸分子から工程ａ〜ｄを経て得られるコンセンサスリードの集合に含まれる各コンセンサスリードの配列情報より、参照配列上の特定位置に対応する位置に存在する塩基種とその信頼性スコアを選抜する工程、
ｆ）工程ｅで使用した特定位置を含むコンセンサスリードの総数に対する、工程ｅで選抜した特定位置に特定の各塩基種を含有するコンセンサスリード数の比率を取得する工程、
ｇ）工程ｄで取得した信頼性スコアを用いたパラメトリックなエラーモデルにより、設定された閾値または有意水準での工程ｅで選抜した特定位置において解析エラーによって特定の塩基が出現する出現率を推定する工程、
ｈ）工程ｆで得られた比率が工程ｇにより得られた出現率より有意に高い場合に、工程ｅで選抜した特定位置に特定の塩基種を有する核酸分子が、工程ｇで設定された閾値または有意水準でクローナル核酸混合物中に存在すると推定する工程、
を含む方法。
【０１７４】
（２）コンピュータを用いて、クローナル核酸混合物を構成する核酸分子の塩基配列における特定位置に存在する特定塩基種がクローナル核酸混合物中で存在する割合を推定する方法であって、
ａ）各核酸分子に由来するリード配列を取得する工程、
ｂ）工程ａで取得したリード配列を参照配列上にマッピングしマッピング結果を取得する工程、
ｃ）工程ｂで得られたマッピング結果を基にリード配列のクラスタリング結果を取得する工程、
ｄ）クラスタリング結果を基にコンセンサスリードとコンセンサスリード中の各塩基に対応する信頼性スコアからなるコンセンサスリード情報を取得する工程、
ｅ）クローナル核酸混合物を構成する各核酸分子から工程ａ〜ｄを経て得られるコンセンサスリードの集合に含まれる各コンセンサスリードの配列情報より、参照配列上の特定位置に対応する位置に存在する塩基種とその信頼性スコアを選抜する工程、
ｆ）工程ｅで使用した特定位置を含むコンセンサスリードの総数に対する、工程ｅで選抜した特定位置に特有の各塩基種を含有するコンセンサスリード数の比率を取得する工程、
ｇ）工程ｄで取得した信頼性スコアを用いたパラメトリックなエラーモデルにより、設定された閾値または有意水準での工程ｅで選抜した特定位置において解析エラーによって特有の塩基が出現する出現率を推定する工程、
ｈ）工程ｆで得られた比率が工程ｇにより得られた出現率より有意に高い場合に、工程ｅで選抜した特定位置に特定の塩基種を有する核酸分子が、工程ｇで設定された閾値または有意水準でクローナル核酸混合物中に存在すると推定する工程、及び
ｉ）工程ｈにおいて工程ｅで選抜した塩基種が存在すると判断された場合、工程ｆで得られた特定位置において塩基種を含有するコンセンサスリードの数から、工程ｇで設定された閾値または有意水準における工程ｅで選抜した特定位置における特定塩基種の比率を算出し、クローナル核酸混合物中で存在する割合として採用する工程、
を含む方法。
【０１７５】
（３）コンピュータを用いて、クローナル核酸混合物を構成する各核酸分子の特定領域の塩基配列推定方法であって、
ａ）各核酸分子に由来するリード配列を取得する工程、
ｂ）工程ａで取得したリード配列を参照配列上にマッピングしマッピング結果を取得する工程、
ｃ）工程ｂで得られたマッピング結果を基にリード配列のクラスタリング結果を取得する工程、
ｄ）クラスタリング結果を基にコンセンサスリードとコンセンサスリード中の各塩基に対応する信頼性スコアからなるコンセンサスリード情報を取得する工程、
ｅ）クローナル核酸混合物を構成する各核酸分子から工程ａ〜ｄを経て得られるコンセンサスリードの集合に含まれる各コンセンサスリードの配列情報より、参照配列上の特定位置に対応する位置に存在する塩基種とその信頼性スコアを選抜する工程、
ｆ）工程ｅで使用した特定位置を含むコンセンサスリードの総数に対する、工程ｅで選抜した特定位置に特定の各塩基種を含有するコンセンサスリード数の比率を取得する工程、
ｇ）工程ｅで選抜した特定位置の塩基種に対応する工程ｄで取得した信頼性スコアを用いたパラメトリックなエラーモデルにより、設定された閾値または有意水準での工程ｅで確認した特定位置における塩基種の解析エラーによる出現率を推定する工程、
ｈ）工程ｆで得られた比率が工程ｇにより得られた出現率より有意に高い場合に、工程ｅで選抜した特定位置に特定の塩基種を有する核酸分子が、工程ｇで設定された閾値または有意水準でクローナル核酸混合物中に存在する、あるいは該出現率より低い場合は特定位置の塩基種は不明、と判断する工程、
ｉ）工程ｈで得られた結果を該塩基の各コンセンサスリードに適用し二次配列情報を作成する工程、及び
ｊ）工程ｅにおける塩基位置を変更し、工程ｅ〜ｉを繰り返すことにより得られた特定領域の全二次配列情報を統合し、クローナル核酸混合物を構成する各核酸分子の特定領域の塩基配列として採用する工程、
を含む方法。
【０１７６】
（４）（３）の方法において、工程ｂで採用した参照配列とコンセンサスリード上の塩基配列が同じである場合は工程ｅ〜工程ｈを経ずにそのまま工程ｉで作成する二次情報の塩基種として採用してもよい。
【０１７７】
（５）（１）〜（３）のいずれかの方法において、工程ｅで選抜する特定の塩基位置が、工程ｂで採用した参照配列とコンセンサスリード上の塩基種が異なる位置であってもよい。
【０１７８】
（６）（１）〜（３）のいずれかの方法において、核酸分子が少なくとも一方の末端に分子バーコードを有してもよい。
【０１７９】
（７）（１）〜（６）のいずれかの方法において、各核酸分子が異なる分子バーコードを有してもよい。
【０１８０】
（８）（７）の方法において、各核酸分子が一つ以上の異なる分子バーコードを有してもよい。
【０１８１】
（９）（１）〜（８）のいずれかの方法において、工程ｃにおけるクラスタリングをリード配列における分子バーコード領域の配列を指標に実施してもよい。
【０１８２】
（１０）（９）の方法において、分子バーコードがランダムな並びの塩基配列であってもよい。
【０１８３】
（１１）（１０）の方法において、分子バーコード領域のリード配列におけるミスマッチが、該分子バーコード領域のリード配列間の配列類似度が許容範囲内であってもよい。
【０１８４】
（１２）（１）〜（１１）のいずれかの方法において、信頼性スコアが、クラスタリングされたリード配列の情報から、ベイズアプローチによって取得されてもよい。
【０１８５】
（１３）（１２）の方法において、リード配列が核酸分子の増幅産物のリード配列であってもよい。
【０１８６】
（１４）（１３）の方法において、増幅物がポリメラーゼ連鎖反応により得られてもよい。
【０１８７】
（１５）（１３）または（１４）の方法において、ベイズアプローチの尤度関数はシーケンサーエラーと増幅エラーを考慮した関数であってもよい。
【０１８８】
（１６）（１５）の方法において、工程ｂにおけるマッピング工程において参照配列に対しリード配列で異なる塩基が挿入または欠失に相当すると判断された場合、シーケンサーエラーに関し更に補正を加えた関数を採用してもよい。
【０１８９】
（１７）（１）〜（３）のいずれかの方法において、工程ｇのパラメトリックなエラーモデルがポアソン分布に基づいてもよい。
【０１９０】
（１８）（１）または（２）の方法において、特定位置の特定塩基種が低頻度な変異の対象塩基であってもよい。
【０１９１】
（１９）（１）〜（３）のいずれかの方法において、工程ｃから工程ｄへ直接進んでもよい。
【０１９２】
（２０）（１）〜（３）のいずれかの方法において、工程ｃの次に、
ｋ）工程ｃで得られた各クラスタリングに含まれるリード配列を用いて、ローカルアセンブリを実行し、アセンブル配列のクラスタリング結果を取得する工程、
を行ってもよい。
【０１９３】
（２１）（１９）及び（２０）の方法によって得られた結果を基に、クローナル核酸混合物を構成する拡散分子の塩基配列における特定位置に存在する特定塩基種を決定する工程を含む方法であってもよい。
【０１９４】
（２２）（１）ないし（２１）のいずれかの方法をコンピュータで実行させるためのコンピュータ読み取り可能なプログラム。
【０１９５】
（２３）クローナル核酸混合物を構成する核酸分子の塩基配列における特定位置の特定塩基種の存在を推定するための装置であって、
各核酸分子由来のリード配列を取得するリード配列取得部と、
取得したリード配列をマッピング装置に供給し、前記マッピング装置による参照配列上へのマッピング結果を取得するマッピング情報取得部と、
マッピングしたリード配列のクラスタリング結果を得るクラスタリング結果取得部と、
クラスタリング結果を基にコンセンサスリードとコンセンサスリード中の各塩基に対応する信頼性スコアからなるコンセンサスリード情報取得部と、
を備えた、推定装置。
【０１９６】
（２４）クローナル核酸混合物を構成する核酸分子の塩基配列における特定位置の特定塩基種の存在を推定するための装置であって、
各核酸分子由来のリード配列を取得するリード配列取得部と、
取得したリード配列をマッピング装置に供給し、前記マッピング装置による参照配列上へのマッピング結果を取得するマッピング情報取得部と、
マッピングしたリード配列のクラスタリング結果を得る第１のクラスタリング結果取得部と、
各クラスタリングに含まれるリード配列を用いて、ローカルアセンブリを実行し、アセンブル配列のクラスタリング結果を得る第２のクラスタリング結果取得部と、
第２のクラスタリング結果取得部のクラスタリング結果を基にコンセンサスリードとコンセンサスリード中の各塩基に対応する信頼性スコアからなるコンセンサスリード情報取得部と、
を備えた、推定装置。
【符号の説明】
【０１９７】
１０  変異検出装置
１１  制御部
１３  表示部
１５  操作部
１６  ＲＡＭ
１７  データ記憶部
１８  インタフェース部
５０  シーケンサー
１００  変異検出システム
【図１】
2019132010.tif 000022
【図２】
2019132010.tif 000023
【図３】
2019132010.tif 000024
【図４】
2019132010.tif 000025
【図５】
2019132010.tif 000026
【図６】
2019132010.tif 000027
【図７】
2019132010.tif 000028
【図８】
2019132010.tif 000029
【図９】
2019132010.tif 000030
【図１０】
2019132010.tif 000031
【図１１】
2019132010.tif 000032
【図１２】
2019132010.tif 000033
【図１３】
2019132010.tif 000034
【図１４】
2019132010.tif 000035
【図１５】
2019132010.tif 000036
【国際調査報告】
2019132010.tif 000037
2019132010.tif 000038
2019132010.tif 000039
2019132010.tif 000040