JPWO2019132010A1 - 塩基配列における塩基種を推定する方法、装置及びプログラム - Google Patents
塩基配列における塩基種を推定する方法、装置及びプログラム Download PDFInfo
- Publication number
- JPWO2019132010A1 JPWO2019132010A1 JP2019562510A JP2019562510A JPWO2019132010A1 JP WO2019132010 A1 JPWO2019132010 A1 JP WO2019132010A1 JP 2019562510 A JP2019562510 A JP 2019562510A JP 2019562510 A JP2019562510 A JP 2019562510A JP WO2019132010 A1 JPWO2019132010 A1 JP WO2019132010A1
- Authority
- JP
- Japan
- Prior art keywords
- sequence
- base
- nucleic acid
- read
- consensus
- 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.)
- Pending
Links
- 238000000034 method Methods 0.000 title claims abstract description 178
- 108020004707 nucleic acids Proteins 0.000 claims abstract description 197
- 102000039446 nucleic acids Human genes 0.000 claims abstract description 197
- 150000007523 nucleic acids Chemical class 0.000 claims abstract description 197
- 238000013507 mapping Methods 0.000 claims abstract description 87
- 239000000203 mixture Substances 0.000 claims abstract description 62
- 238000004458 analytical method Methods 0.000 claims abstract description 45
- 230000035772 mutation Effects 0.000 claims description 149
- 241000894007 species Species 0.000 claims description 74
- 230000008569 process Effects 0.000 claims description 68
- 235000014676 Phragmites communis Nutrition 0.000 claims description 59
- 238000003752 polymerase chain reaction Methods 0.000 claims description 56
- 238000009826 distribution Methods 0.000 claims description 29
- GNFTZDOKVXKIBK-UHFFFAOYSA-N 3-(2-methoxyethoxy)benzohydrazide Chemical compound COCCOC1=CC=CC(C(=O)NN)=C1 GNFTZDOKVXKIBK-UHFFFAOYSA-N 0.000 claims description 23
- 244000273256 Phragmites communis Species 0.000 claims description 16
- 238000003199 nucleic acid amplification method Methods 0.000 claims description 14
- 230000003321 amplification Effects 0.000 claims description 13
- FGUUSXIOTUKUDN-IBGZPJMESA-N C1(=CC=CC=C1)N1C2=C(NC([C@H](C1)NC=1OC(=NN=1)C1=CC=CC=C1)=O)C=CC=C2 Chemical compound C1(=CC=CC=C1)N1C2=C(NC([C@H](C1)NC=1OC(=NN=1)C1=CC=CC=C1)=O)C=CC=C2 FGUUSXIOTUKUDN-IBGZPJMESA-N 0.000 claims description 12
- 238000012217 deletion Methods 0.000 claims description 9
- 230000037430 deletion Effects 0.000 claims description 9
- 238000003780 insertion Methods 0.000 claims description 9
- 230000037431 insertion Effects 0.000 claims description 9
- 238000013476 bayesian approach Methods 0.000 claims description 6
- 238000009792 diffusion process Methods 0.000 claims description 3
- 238000001514 detection method Methods 0.000 description 68
- 108020004414 DNA Proteins 0.000 description 22
- 238000012545 processing Methods 0.000 description 22
- 238000012408 PCR amplification Methods 0.000 description 19
- 108090000623 proteins and genes Proteins 0.000 description 19
- 238000007405 data analysis Methods 0.000 description 16
- 239000000523 sample Substances 0.000 description 16
- 230000006870 function Effects 0.000 description 14
- 238000010586 diagram Methods 0.000 description 13
- 206010028980 Neoplasm Diseases 0.000 description 12
- 230000035945 sensitivity Effects 0.000 description 12
- 102000054766 genetic haplotypes Human genes 0.000 description 11
- 238000002156 mixing Methods 0.000 description 10
- 238000012163 sequencing technique Methods 0.000 description 10
- 108700028369 Alleles Proteins 0.000 description 9
- 210000004027 cell Anatomy 0.000 description 9
- 108091028043 Nucleic acid sequence Proteins 0.000 description 8
- 230000036438 mutation frequency Effects 0.000 description 8
- 201000011510 cancer Diseases 0.000 description 7
- 108090000790 Enzymes Proteins 0.000 description 6
- 102000004190 Enzymes Human genes 0.000 description 6
- 238000004422 calculation algorithm Methods 0.000 description 5
- 210000000349 chromosome Anatomy 0.000 description 5
- 238000007796 conventional method Methods 0.000 description 5
- 230000037429 base substitution Effects 0.000 description 4
- 238000004364 calculation method Methods 0.000 description 4
- 238000004891 communication Methods 0.000 description 4
- 238000013500 data storage Methods 0.000 description 4
- 238000005516 engineering process Methods 0.000 description 4
- 238000001914 filtration Methods 0.000 description 4
- 238000013467 fragmentation Methods 0.000 description 4
- 238000006062 fragmentation reaction Methods 0.000 description 4
- 230000012447 hatching Effects 0.000 description 4
- 239000002773 nucleotide Substances 0.000 description 4
- 125000003729 nucleotide group Chemical group 0.000 description 4
- 238000013179 statistical model Methods 0.000 description 4
- 108091032973 (ribonucleotides)n+m Proteins 0.000 description 3
- 101150039808 Egfr gene Proteins 0.000 description 3
- 101100495925 Schizosaccharomyces pombe (strain 972 / ATCC 24843) chr3 gene Proteins 0.000 description 3
- 238000012300 Sequence Analysis Methods 0.000 description 3
- 239000012472 biological sample Substances 0.000 description 3
- 238000010276 construction Methods 0.000 description 3
- 238000012937 correction Methods 0.000 description 3
- 108700021358 erbB-1 Genes Proteins 0.000 description 3
- 239000012634 fragment Substances 0.000 description 3
- 238000006467 substitution reaction Methods 0.000 description 3
- 101150044182 8 gene Proteins 0.000 description 2
- 101000605639 Homo sapiens Phosphatidylinositol 4,5-bisphosphate 3-kinase catalytic subunit alpha isoform Proteins 0.000 description 2
- 101150073096 NRAS gene Proteins 0.000 description 2
- 102100038332 Phosphatidylinositol 4,5-bisphosphate 3-kinase catalytic subunit alpha isoform Human genes 0.000 description 2
- 230000002159 abnormal effect Effects 0.000 description 2
- 230000005540 biological transmission Effects 0.000 description 2
- 210000004369 blood Anatomy 0.000 description 2
- 239000008280 blood Substances 0.000 description 2
- 238000012512 characterization method Methods 0.000 description 2
- 239000003153 chemical reaction reagent Substances 0.000 description 2
- 235000019506 cigar Nutrition 0.000 description 2
- 238000003776 cleavage reaction Methods 0.000 description 2
- OPTASPLRGRRNAP-UHFFFAOYSA-N cytosine Chemical compound NC=1C=CNC(=O)N=1 OPTASPLRGRRNAP-UHFFFAOYSA-N 0.000 description 2
- 230000002950 deficient Effects 0.000 description 2
- 230000009977 dual effect Effects 0.000 description 2
- 238000000605 extraction Methods 0.000 description 2
- 230000009191 jumping Effects 0.000 description 2
- 238000007481 next generation sequencing Methods 0.000 description 2
- 230000003287 optical effect Effects 0.000 description 2
- 230000037361 pathway Effects 0.000 description 2
- 238000011002 quantification Methods 0.000 description 2
- 238000005067 remediation Methods 0.000 description 2
- 102200007380 rs121913528 Human genes 0.000 description 2
- 230000007017 scission Effects 0.000 description 2
- 210000001519 tissue Anatomy 0.000 description 2
- 238000009966 trimming Methods 0.000 description 2
- 101150040471 19 gene Proteins 0.000 description 1
- 125000003903 2-propenyl group Chemical group [H]C([*])([H])C([H])=C([H])[H] 0.000 description 1
- 108091035707 Consensus sequence Proteins 0.000 description 1
- 230000033616 DNA repair Effects 0.000 description 1
- 102100030708 GTPase KRas Human genes 0.000 description 1
- 102100039788 GTPase NRas Human genes 0.000 description 1
- 241000282412 Homo Species 0.000 description 1
- 101000584612 Homo sapiens GTPase KRas Proteins 0.000 description 1
- 101000744505 Homo sapiens GTPase NRas Proteins 0.000 description 1
- 102000003960 Ligases Human genes 0.000 description 1
- 108090000364 Ligases Proteins 0.000 description 1
- 238000007476 Maximum Likelihood Methods 0.000 description 1
- 108700020796 Oncogene Proteins 0.000 description 1
- 102000043276 Oncogene Human genes 0.000 description 1
- 101150063858 Pik3ca gene Proteins 0.000 description 1
- 241000532838 Platypus Species 0.000 description 1
- 206010036790 Productive cough Diseases 0.000 description 1
- 102000044209 Tumor Suppressor Genes Human genes 0.000 description 1
- 108700025716 Tumor Suppressor Genes Proteins 0.000 description 1
- 239000002253 acid Substances 0.000 description 1
- 150000007513 acids Chemical class 0.000 description 1
- 230000004931 aggregating effect Effects 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 230000000712 assembly Effects 0.000 description 1
- 238000000429 assembly Methods 0.000 description 1
- 230000015556 catabolic process Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 238000004140 cleaning Methods 0.000 description 1
- 238000004590 computer program Methods 0.000 description 1
- 239000000470 constituent Substances 0.000 description 1
- 101150091306 corC gene Proteins 0.000 description 1
- 238000005520 cutting process Methods 0.000 description 1
- 229940104302 cytosine Drugs 0.000 description 1
- 230000009615 deamination Effects 0.000 description 1
- 238000006481 deamination reaction Methods 0.000 description 1
- 230000007812 deficiency Effects 0.000 description 1
- 238000006731 degradation reaction Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000018109 developmental process Effects 0.000 description 1
- 238000007847 digital PCR Methods 0.000 description 1
- 201000010099 disease Diseases 0.000 description 1
- 208000037265 diseases, disorders, signs and symptoms Diseases 0.000 description 1
- 238000011304 droplet digital PCR Methods 0.000 description 1
- 239000003814 drug Substances 0.000 description 1
- 229940079593 drug Drugs 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 102000052116 epidermal growth factor receptor activity proteins Human genes 0.000 description 1
- 108700015053 epidermal growth factor receptor activity proteins Proteins 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 210000003608 fece Anatomy 0.000 description 1
- 230000005861 gene abnormality Effects 0.000 description 1
- 238000000338 in vitro Methods 0.000 description 1
- 230000010365 information processing Effects 0.000 description 1
- 230000010354 integration Effects 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 230000002601 intratumoral effect Effects 0.000 description 1
- 238000011528 liquid biopsy Methods 0.000 description 1
- 239000004973 liquid crystal related substance Substances 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 239000003550 marker Substances 0.000 description 1
- 230000000873 masking effect Effects 0.000 description 1
- YOHYSYJDKVYCJI-UHFFFAOYSA-N n-[3-[[6-[3-(trifluoromethyl)anilino]pyrimidin-4-yl]amino]phenyl]cyclopropanecarboxamide Chemical compound FC(F)(F)C1=CC=CC(NC=2N=CN=C(NC=3C=C(NC(=O)C4CC4)C=CC=3)C=2)=C1 YOHYSYJDKVYCJI-UHFFFAOYSA-N 0.000 description 1
- 238000003012 network analysis Methods 0.000 description 1
- 108091027963 non-coding RNA Proteins 0.000 description 1
- 102000042567 non-coding RNA Human genes 0.000 description 1
- 238000007254 oxidation reaction Methods 0.000 description 1
- 238000013081 phylogenetic analysis Methods 0.000 description 1
- 210000002381 plasma Anatomy 0.000 description 1
- 238000002360 preparation method Methods 0.000 description 1
- 102000004169 proteins and genes Human genes 0.000 description 1
- 239000013643 reference control Substances 0.000 description 1
- 230000008439 repair process Effects 0.000 description 1
- 230000010076 replication Effects 0.000 description 1
- 102200124923 rs121913254 Human genes 0.000 description 1
- 102200006539 rs121913529 Human genes 0.000 description 1
- 239000010979 ruby Substances 0.000 description 1
- 229910001750 ruby Inorganic materials 0.000 description 1
- 210000003296 saliva Anatomy 0.000 description 1
- 239000004065 semiconductor Substances 0.000 description 1
- 210000002966 serum Anatomy 0.000 description 1
- 239000007787 solid Substances 0.000 description 1
- 210000003802 sputum Anatomy 0.000 description 1
- 208000024794 sputum Diseases 0.000 description 1
- 238000007619 statistical method Methods 0.000 description 1
- 238000003860 storage Methods 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- -1 tissues Substances 0.000 description 1
- 238000009827 uniform distribution Methods 0.000 description 1
- 210000002700 urine Anatomy 0.000 description 1
- 239000013598 vector Substances 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
- G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
- G16B30/00—ICT specially adapted for sequence analysis involving nucleotides or amino acids
- G16B30/20—Sequence assembly
-
- C—CHEMISTRY; METALLURGY
- C12—BIOCHEMISTRY; BEER; SPIRITS; WINE; VINEGAR; MICROBIOLOGY; ENZYMOLOGY; MUTATION OR GENETIC ENGINEERING
- C12N—MICROORGANISMS OR ENZYMES; COMPOSITIONS THEREOF; PROPAGATING, PRESERVING, OR MAINTAINING MICROORGANISMS; MUTATION OR GENETIC ENGINEERING; CULTURE MEDIA
- C12N15/00—Mutation or genetic engineering; DNA or RNA concerning genetic engineering, vectors, e.g. plasmids, or their isolation, preparation or purification; Use of hosts therefor
-
- C—CHEMISTRY; METALLURGY
- C12—BIOCHEMISTRY; BEER; SPIRITS; WINE; VINEGAR; MICROBIOLOGY; ENZYMOLOGY; MUTATION OR GENETIC ENGINEERING
- C12Q—MEASURING OR TESTING PROCESSES INVOLVING ENZYMES, NUCLEIC ACIDS OR MICROORGANISMS; COMPOSITIONS OR TEST PAPERS THEREFOR; PROCESSES OF PREPARING SUCH COMPOSITIONS; CONDITION-RESPONSIVE CONTROL IN MICROBIOLOGICAL OR ENZYMOLOGICAL PROCESSES
- C12Q1/00—Measuring or testing processes involving enzymes, nucleic acids or microorganisms; Compositions therefor; Processes of preparing such compositions
- C12Q1/68—Measuring or testing processes involving enzymes, nucleic acids or microorganisms; Compositions therefor; Processes of preparing such compositions involving nucleic acids
- C12Q1/6806—Preparing nucleic acids for analysis, e.g. for polymerase chain reaction [PCR] assay
Landscapes
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Chemical & Material Sciences (AREA)
- Analytical Chemistry (AREA)
- Biophysics (AREA)
- Proteomics, Peptides & Aminoacids (AREA)
- Health & Medical Sciences (AREA)
- Engineering & Computer Science (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Bioinformatics & Computational Biology (AREA)
- Biotechnology (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)
- Apparatus Associated With Microorganisms And Enzymes (AREA)
Abstract
塩基配列における塩基種を推定する方法、装置、プログラムを提供する。前記方法では、各核酸分子由来のリード配列を取得し(S11)、取得したリード配列を参照配列上にマッピングし(S12)、マッピング結果に基づきリード配列を分類し(S13)、分類結果に基づきコンセンサスリードとコンセンサスリード中の各塩基に対応する信頼性スコアを求め(S14)、コンセンサスリードより特定位置における塩基種を任意に選択し、コンセンサスリード総数に対する、選択した特定位置の各塩基種を含むコンセンサスリード数を取得し、選択した特定位置の塩基種に対応する信頼性スコアを用いたパラメトリックなエラーモデルにより(S15)、設定された信頼性で特定位置における塩基種の解析エラーによる出現率を推定し、出現率に基づき、特定位置における特定の塩基種が、設定された信頼性でクローナル核酸混合物中に存在するか否かを推定する(S16)。
Description
本発明は、塩基配列における塩基種を推定する技術に関する。
がん細胞は、正常な細胞中のがん遺伝子及びがん抑制遺伝子(がん関連遺伝子)に塩基配列の突然変異が生じることにより発生すると言われている。近年、これらのがん関連遺伝子の遺伝子異常を応用して、正常な細胞とがん細胞を見分ける手法の開発が進められている。例えば、非特許文献1および非特許文献2では、血中や便中に存在する微量ながん細胞由来の異常なDNAを検出できたことが開示されている。さらに、がん関連遺伝子の変異によって翻訳された異常タンパク質を標的とする分子標的薬の登場により、塩基配列を高精度に決定することが求められる。
次世代シーケンサーによる腫瘍検体の変異の包括的なプロファイリングの結果、一つの腫瘍のなかには、正常な細胞とは異なる配列のゲノムを持つ、がん細胞由来の複数のゲノムクローンが存在することが報告され、この現象は腫瘍内不均一性と呼ばれている。この現象により、腫瘍自体がポリクローナルな混合状態となっていることから、低頻度な変異も含めて網羅的に検出する必要性がある。
現在の次世代シーケンス技術では、核酸増幅時及びシーケンシング時にある頻度でエラーが生じる。したがって、シーケンス解析を実施したゲノムDNAの塩基配列には技術的なエラーに起因する偽陽性の変異(false positive variant)が含まれるため、1〜5%以下の低頻度な変異の検出は難しい。しかし近年、分子バーコード技術を応用することによって、低頻度な変異検出時のエラーが抑制されることが報告されている。
例えば、非特許文献3には、配列解析装置より得られた塩基配列情報の統計処理により検出精度を高める技術としてsmCounterが開示されているが、リキッドバイオプシーなど、より高い検出感度と精度を求められる場合には精度が足らず、不十分である。
Armengol G等、J Mol Diagn. 2016 Jul;18(4):471−9. pmid: 27155048
Diehl F等、Nat Med. 2008;14(9): 985-90. pmid:18670422
Chang Xu等、BMC Genomics、第18巻、e5ページ、2017年(doi:10.1186/s12864−016−3425−4))
本発明は、塩基配列における塩基種を推定する方法、装置、及びその方法を実施するためのコンピュータープログラムを提供することを目的とする。
本発明者は、テンプレートDNAに分子バーコードを付加させて調製したDNA断片から得られた次世代シーケンサーの多量の塩基配列データ(リード情報)から、分子バーコード毎に信頼性スコアを決定し、パラメトリック解析することで、従来法よりも高い感度で変異を検出できることを見出した。
また、本発明者は、シーケンスデータから低頻度に存在する、特定位置の塩基種を推定するために、核酸増幅又はシーケンシングで発生するエラーを効率的に取り除く指標を算出し、真の塩基種を高確度に判定する統計手法を開発した。
本発明の検出精度を高める技術は、例えば、より疾患早期ステージでの変異体核酸の検出を期待できる。
本開示の第一の態様において、コンピュータを用いて、クローナル核酸混合物を構成する核酸分子の塩基配列における特定位置に存在する特定塩基種の存在を推定する第一の推定方法が提供される。
第一の推定方法は、
a)各核酸分子に由来するリード配列を取得する工程、
b)工程aで取得したリード配列を参照配列上にマッピングしマッピング結果を取得する工程、
c)工程bで得られたマッピング結果を基にリード配列のクラスタリング結果を取得する工程、
d)クラスタリング結果を基にコンセンサスリードとコンセンサスリード中の各塩基に対応する信頼性スコアからなるコンセンサスリード情報を取得する工程、
e)クローナル核酸混合物を構成する各核酸分子から工程a〜dを経て得られるコンセンサスリードの集合に含まれる各コンセンサスリードの配列情報より、参照配列上の特定位置に対応する位置に存在する塩基種とその信頼性スコアを選抜する工程、
f)工程eで使用した特定位置を含むコンセンサスリードの総数に対する、工程eで選抜した特定位置に特定の各塩基種を含有するコンセンサスリード数の比率を取得する工程、
g)工程dで取得した信頼性スコアを用いたパラメトリックなエラーモデルにより、設定された閾値または有意水準での工程eで選抜した特定位置において解析エラーによって特定の塩基が出現する出現率を推定する工程、
h)工程fで得られた比率が工程gにより得られた出現率より有意に高い場合に、工程eで選抜した特定位置に特定の塩基種を有する核酸分子が、工程gで設定された閾値または有意水準でクローナル核酸混合物中に存在すると推定する工程、
を含む。
第一の推定方法は、
a)各核酸分子に由来するリード配列を取得する工程、
b)工程aで取得したリード配列を参照配列上にマッピングしマッピング結果を取得する工程、
c)工程bで得られたマッピング結果を基にリード配列のクラスタリング結果を取得する工程、
d)クラスタリング結果を基にコンセンサスリードとコンセンサスリード中の各塩基に対応する信頼性スコアからなるコンセンサスリード情報を取得する工程、
e)クローナル核酸混合物を構成する各核酸分子から工程a〜dを経て得られるコンセンサスリードの集合に含まれる各コンセンサスリードの配列情報より、参照配列上の特定位置に対応する位置に存在する塩基種とその信頼性スコアを選抜する工程、
f)工程eで使用した特定位置を含むコンセンサスリードの総数に対する、工程eで選抜した特定位置に特定の各塩基種を含有するコンセンサスリード数の比率を取得する工程、
g)工程dで取得した信頼性スコアを用いたパラメトリックなエラーモデルにより、設定された閾値または有意水準での工程eで選抜した特定位置において解析エラーによって特定の塩基が出現する出現率を推定する工程、
h)工程fで得られた比率が工程gにより得られた出現率より有意に高い場合に、工程eで選抜した特定位置に特定の塩基種を有する核酸分子が、工程gで設定された閾値または有意水準でクローナル核酸混合物中に存在すると推定する工程、
を含む。
本開示の第二の態様において、コンピュータを用いて、クローナル核酸混合物を構成する核酸分子の塩基配列における特定位置に存在する特定塩基種がクローナル核酸混合物中で存在する割合を推定する第二の推定方法が提供される。
第二の推定方法は、
a)各核酸分子に由来するリード配列を取得する工程、
b)工程aで取得したリード配列を参照配列上にマッピングしマッピング結果を取得する工程、
c)工程bで得られたマッピング結果を基にリード配列のクラスタリング結果を取得する工程、
d)クラスタリング結果を基にコンセンサスリードとコンセンサスリード中の各塩基に対応する信頼性スコアからなるコンセンサスリード情報を取得する工程、
e)クローナル核酸混合物を構成する各核酸分子から工程a〜dを経て得られるコンセンサスリードの集合に含まれる各コンセンサスリードの配列情報より、参照配列上の特定位置に対応する位置に存在する塩基種とその信頼性スコアを選抜する工程、
f)工程eで使用した特定位置を含むコンセンサスリードの総数に対する、工程eで選抜した特定位置に特有の各塩基種を含有するコンセンサスリード数の比率を取得する工程、
g)工程dで取得した信頼性スコアを用いたパラメトリックなエラーモデルにより、設定された閾値または有意水準での工程eで選抜した特定位置において解析エラーによって特有の塩基が出現する出現率を推定する工程、
h)工程fで得られた比率が工程gにより得られた出現率より有意に高い場合に、工程eで選抜した特定位置に特定の塩基種を有する核酸分子が、工程gで設定された閾値または有意水準でクローナル核酸混合物中に存在すると推定する工程、及び
i)工程hにおいて工程eで選抜した塩基種が存在すると判断された場合、工程fで得られた特定位置において塩基種を含有するコンセンサスリードの数から、工程gで設定された閾値または有意水準における工程eで選抜した特定位置における特定塩基種の比率を算出し、クローナル核酸混合物中で存在する割合として採用する工程、
を含む。
第二の推定方法は、
a)各核酸分子に由来するリード配列を取得する工程、
b)工程aで取得したリード配列を参照配列上にマッピングしマッピング結果を取得する工程、
c)工程bで得られたマッピング結果を基にリード配列のクラスタリング結果を取得する工程、
d)クラスタリング結果を基にコンセンサスリードとコンセンサスリード中の各塩基に対応する信頼性スコアからなるコンセンサスリード情報を取得する工程、
e)クローナル核酸混合物を構成する各核酸分子から工程a〜dを経て得られるコンセンサスリードの集合に含まれる各コンセンサスリードの配列情報より、参照配列上の特定位置に対応する位置に存在する塩基種とその信頼性スコアを選抜する工程、
f)工程eで使用した特定位置を含むコンセンサスリードの総数に対する、工程eで選抜した特定位置に特有の各塩基種を含有するコンセンサスリード数の比率を取得する工程、
g)工程dで取得した信頼性スコアを用いたパラメトリックなエラーモデルにより、設定された閾値または有意水準での工程eで選抜した特定位置において解析エラーによって特有の塩基が出現する出現率を推定する工程、
h)工程fで得られた比率が工程gにより得られた出現率より有意に高い場合に、工程eで選抜した特定位置に特定の塩基種を有する核酸分子が、工程gで設定された閾値または有意水準でクローナル核酸混合物中に存在すると推定する工程、及び
i)工程hにおいて工程eで選抜した塩基種が存在すると判断された場合、工程fで得られた特定位置において塩基種を含有するコンセンサスリードの数から、工程gで設定された閾値または有意水準における工程eで選抜した特定位置における特定塩基種の比率を算出し、クローナル核酸混合物中で存在する割合として採用する工程、
を含む。
本開示の第三の態様において、コンピュータを用いて、クローナル核酸混合物を構成する各核酸分子の特定領域の塩基配列を推定する第三の推定方法が提供される。
第三の推定方法は、
a)各核酸分子に由来するリード配列を取得する工程、
b)工程aで取得したリード配列を参照配列上にマッピングしマッピング結果を取得する工程、
c)工程bで得られたマッピング結果を基にリード配列のクラスタリング結果を取得する工程、
d)クラスタリング結果を基にコンセンサスリードとコンセンサスリード中の各塩基に対応する信頼性スコアからなるコンセンサスリード情報を取得する工程、
e)クローナル核酸混合物を構成する各核酸分子から工程a〜dを経て得られるコンセンサスリードの集合に含まれる各コンセンサスリードの配列情報より、参照配列上の特定位置に対応する位置に存在する塩基種とその信頼性スコアを選抜する工程、
f)工程eで使用した特定位置を含むコンセンサスリードの総数に対する、工程eで選抜した特定位置に特定の各塩基種を含有するコンセンサスリード数の比率を取得する工程、
g)工程eで選抜した特定位置の塩基種に対応する工程dで取得した信頼性スコアを用いたパラメトリックなエラーモデルにより、設定された閾値または有意水準での工程eで確認した特定位置における塩基種の解析エラーによる出現率を推定する工程、
h)工程fで得られた比率が工程gにより得られた出現率より有意に高い場合に、工程eで選抜した特定位置に特定の塩基種を有する核酸分子が、工程gで設定された閾値または有意水準でクローナル核酸混合物中に存在する、あるいは該出現率より低い場合は特定位置の塩基種は不明、と判断する工程、
i)工程hで得られた結果を該塩基の各コンセンサスリードに適用し二次配列情報を作成する工程、及び
j)工程eにおける塩基位置を変更し、工程e〜iを繰り返すことにより得られた特定領域の全二次配列情報を統合し、クローナル核酸混合物を構成する各核酸分子の特定領域の塩基配列として採用する工程、
を含む。
第三の推定方法は、
a)各核酸分子に由来するリード配列を取得する工程、
b)工程aで取得したリード配列を参照配列上にマッピングしマッピング結果を取得する工程、
c)工程bで得られたマッピング結果を基にリード配列のクラスタリング結果を取得する工程、
d)クラスタリング結果を基にコンセンサスリードとコンセンサスリード中の各塩基に対応する信頼性スコアからなるコンセンサスリード情報を取得する工程、
e)クローナル核酸混合物を構成する各核酸分子から工程a〜dを経て得られるコンセンサスリードの集合に含まれる各コンセンサスリードの配列情報より、参照配列上の特定位置に対応する位置に存在する塩基種とその信頼性スコアを選抜する工程、
f)工程eで使用した特定位置を含むコンセンサスリードの総数に対する、工程eで選抜した特定位置に特定の各塩基種を含有するコンセンサスリード数の比率を取得する工程、
g)工程eで選抜した特定位置の塩基種に対応する工程dで取得した信頼性スコアを用いたパラメトリックなエラーモデルにより、設定された閾値または有意水準での工程eで確認した特定位置における塩基種の解析エラーによる出現率を推定する工程、
h)工程fで得られた比率が工程gにより得られた出現率より有意に高い場合に、工程eで選抜した特定位置に特定の塩基種を有する核酸分子が、工程gで設定された閾値または有意水準でクローナル核酸混合物中に存在する、あるいは該出現率より低い場合は特定位置の塩基種は不明、と判断する工程、
i)工程hで得られた結果を該塩基の各コンセンサスリードに適用し二次配列情報を作成する工程、及び
j)工程eにおける塩基位置を変更し、工程e〜iを繰り返すことにより得られた特定領域の全二次配列情報を統合し、クローナル核酸混合物を構成する各核酸分子の特定領域の塩基配列として採用する工程、
を含む。
本開示の第四の態様において、上記第一ないし第三の推定方法のいずれかをコンピュータで実行させるためのコンピュータ読み取り可能なプログラムが提供される。
本開示の第五の態様において、クローナル核酸混合物を構成する核酸分子の塩基配列における特定位置の特定塩基種の存在を推定するための推定装置が提供される。
推定装置は、
各核酸分子由来のリード配列を取得するリード配列取得部と、
取得したリード配列をマッピング装置に供給し、前記マッピング装置による参照配列上へのマッピング結果を取得するマッピング情報取得部と、
マッピングしたリード配列のクラスタリング結果を得るクラスタリング結果取得部と、
クラスタリング結果を基にコンセンサスリードとコンセンサスリード中の各塩基に対応する信頼性スコアからなるコンセンサスリード情報取得部と、
を備える。
推定装置は、
各核酸分子由来のリード配列を取得するリード配列取得部と、
取得したリード配列をマッピング装置に供給し、前記マッピング装置による参照配列上へのマッピング結果を取得するマッピング情報取得部と、
マッピングしたリード配列のクラスタリング結果を得るクラスタリング結果取得部と、
クラスタリング結果を基にコンセンサスリードとコンセンサスリード中の各塩基に対応する信頼性スコアからなるコンセンサスリード情報取得部と、
を備える。
本開示の第六の態様において、クローナル核酸混合物を構成する核酸分子の塩基配列における特定位置の特定塩基種の存在を推定するための推定装置が提供される。
推定装置は、
各核酸分子由来のリード配列を取得するリード配列取得部と、
取得したリード配列をマッピング装置に供給し、前記マッピング装置による参照配列上へのマッピング結果を取得するマッピング情報取得部と、
マッピングしたリード配列のクラスタリング結果を得る第1のクラスタリング結果取得部と、
各クラスタリングに含まれるリード配列を用いて、ローカルアセンブリを実行し、アセンブル配列のクラスタリング結果を得る第2のクラスタリング結果取得部と、
第2のクラスタリング結果取得部のクラスタリング結果を基にコンセンサスリードとコンセンサスリード中の各塩基に対応する信頼性スコアからなるコンセンサスリード情報取得部と、
を備える。
推定装置は、
各核酸分子由来のリード配列を取得するリード配列取得部と、
取得したリード配列をマッピング装置に供給し、前記マッピング装置による参照配列上へのマッピング結果を取得するマッピング情報取得部と、
マッピングしたリード配列のクラスタリング結果を得る第1のクラスタリング結果取得部と、
各クラスタリングに含まれるリード配列を用いて、ローカルアセンブリを実行し、アセンブル配列のクラスタリング結果を得る第2のクラスタリング結果取得部と、
第2のクラスタリング結果取得部のクラスタリング結果を基にコンセンサスリードとコンセンサスリード中の各塩基に対応する信頼性スコアからなるコンセンサスリード情報取得部と、
を備える。
本発明によれば、被験体の身体サンプル由来の一部の核酸上に発生する低頻度な変異を検出でき、塩基配列における塩基種を推定できる。
以下、添付の図面を参照しながら本開示の実施の形態を説明する。
<用語の定義>
本開示における「クローナル核酸混合物」とは、核酸一分子から生物学的な複製を繰り返すことにより発生した複数の核酸の混合物を指す。このようなクローナル核酸が得られる試料としては、特に生物種による限定も無いが、例えばヒトの生体試料が挙げられ、組織、血液、血漿、血清、尿、唾液、粘膜排出物、痰、便および涙が例示される。また組織のFFPE(ホルマリン固定パラフィン包埋)試料の他、in vitroでの培養による細胞集団であってもよい。FFPE試料由来DNAは、短く断片化されている、ニックやギャップが入っている、シトシンの脱アミノ化(ウラシル化)が生じているなどダメージを受けていることが多い。そのままでは解析不良の一因となるため、NEBNext FFPE DNA Repair Mix(NEB社製)など市販の試薬を用いて該DNAを修復することが望ましい。
本開示における「クローナル核酸混合物」とは、核酸一分子から生物学的な複製を繰り返すことにより発生した複数の核酸の混合物を指す。このようなクローナル核酸が得られる試料としては、特に生物種による限定も無いが、例えばヒトの生体試料が挙げられ、組織、血液、血漿、血清、尿、唾液、粘膜排出物、痰、便および涙が例示される。また組織のFFPE(ホルマリン固定パラフィン包埋)試料の他、in vitroでの培養による細胞集団であってもよい。FFPE試料由来DNAは、短く断片化されている、ニックやギャップが入っている、シトシンの脱アミノ化(ウラシル化)が生じているなどダメージを受けていることが多い。そのままでは解析不良の一因となるため、NEBNext FFPE DNA Repair Mix(NEB社製)など市販の試薬を用いて該DNAを修復することが望ましい。
本開示における「リード配列」とはシーケンサーより出力された塩基配列データである。
本開示における「低頻度に存在する、特定位置の塩基種」とは、特に限定するものではないが、例えば「低頻度な変異」が挙げられる。ここで、「低頻度な変異」とは、生体内で生じた核酸中の塩基の(突然)変異であって、以下の二つの条件a、bを満たす変異を意図する。
a)核酸分子において、1×10−3/塩基以下の頻度(すなわち、1,000塩基に一つ以下の頻度)で出現する。
b)核酸分子を含む試料において、特定の位置の塩基に異なる配列の塩基を有する核酸分子の割合が試料中の全核酸分子数の10%以下となる。
a)核酸分子において、1×10−3/塩基以下の頻度(すなわち、1,000塩基に一つ以下の頻度)で出現する。
b)核酸分子を含む試料において、特定の位置の塩基に異なる配列の塩基を有する核酸分子の割合が試料中の全核酸分子数の10%以下となる。
以下、「低頻度に存在する、特定位置の塩基種」として「低頻度な変異」を例に挙げて本開示を説明するが、本開示は低頻度な変異の検出に何ら限定されるものではない。
本開示における(突然)変異は、置換、挿入及び欠失のいずれでも良い。また、該(突然)変異は、既知であってもよく、新規であってもよい。
(実施の形態)
1.変異検出のための構成
図1は、本発明の実施の形態に係る塩基配列の推定方法を実施する変異検出システム(推定システムの一例)の構成を示した図である。変異検出システム100はDNA配列における変異を検出する。変異検出システム100は、断片化されたDNAをシーケンサー50で解析し、その解析の結果得られたリード配列データを変異検出装置10(推定装置の一例)に入力する。変異検出装置10はシーケンサー50からリード配列データを取得し、解析してDNA配列における変異を検出する。
1.変異検出のための構成
図1は、本発明の実施の形態に係る塩基配列の推定方法を実施する変異検出システム(推定システムの一例)の構成を示した図である。変異検出システム100はDNA配列における変異を検出する。変異検出システム100は、断片化されたDNAをシーケンサー50で解析し、その解析の結果得られたリード配列データを変異検出装置10(推定装置の一例)に入力する。変異検出装置10はシーケンサー50からリード配列データを取得し、解析してDNA配列における変異を検出する。
2.変異検出方法
図2は、本実施の形態のDNA配列における新規な変異検出方法の手順を示した図である。図2を参照して、DNA配列における変異検出方法の手順を説明する。
図2は、本実施の形態のDNA配列における新規な変異検出方法の手順を示した図である。図2を参照して、DNA配列における変異検出方法の手順を説明する。
まず、生体試料からゲノムDNAやRNA(クローナル核酸)を抽出し、シーケンサーに適した長さのDNA(例えば1000base以下)が含まれるように断片化し、断片化クローナル核酸を生成する(S1)。生体試料が、すでに断片化した状態のDNA(cell−free DNA、FFPEサンプルなど)や短いRNA(non−coding RNAなど)の場合、断片化工程は行わなくともよい。
本発明は特に、低頻度な変異を含む可能性のある核酸の解析に好適に使用される。解析したいクローナル核酸は、DNAまたはRNAであり、好ましくはDNAであり、さらに好ましくはゲノムDNAである。
次世代シーケンサーの機種によっては解析可能な配列長が数百bpと短いものもある。このため、例えば、クローナル核酸がゲノムDNAの場合、シーケンサーで解析可能な長さに断片化する必要がある。断片化の方法に限定は無く、公知の方法で断片化すればよい。本発明を何ら限定するものではないが、例えば、超音波を用いた物理的な切断を行うDNA Shearing System M220/ME220(Covaris社製)、酵素を用いた切断を行うQIAseq FX DNA Library Kit(Qiagen社製)など市販の機器や試薬が使用できる。好適には、超音波を用いた物理的な切断が用いられる。
次に各断片化クローナル核酸に「分子バーコード」と呼ばれる核酸を付加する(S2)。これにより得られたテンプレート核酸はのちにコピーされるが、分子バーコードを付加することにより、コピーがどのクローナル核酸由来かを識別することが可能となる。通常、クローナル核酸を構成する各核酸分子のコピーを作製し、当該コピーのリード配列を得る場合、リード配列の塩基情報のみから本来の核酸分子単位の由来を特定することは不可能である。これに対処するため、該分子バーコードを用いる。本開示において、クローナル核酸に付与される分子バーコードは、一つ以上でユニークであれば良い。また、その配列長や配列パターンは、特に限定されないが、1000種類のクローナル核酸を識別するには最低5塩基長(45=1024通り)が必要であり、通常6〜10塩基長の配列で使用される。なお分子バーコードは様々な呼称があり、“unique molecular identifier(UMI)”や“unique molecular tag(UMT)”、または単に“分子タグ(molecular tag)”と呼ばれることもある。クローナル核酸への分子バーコードの付加方法に特に限定は無いが、リガーゼ等の酵素により結合させることで調製することができる。
以下、分子バーコードを付加したクローナル核酸を、テンプレート核酸または解析対象核酸と称することがある。
テンプレート核酸に分子バーコードの他、PCR(polymerase chain reaction)増幅用プライマーがアリーニング可能な配列を加えておくと後述の核酸増幅操作が行いやすくなる。このようなテンプレート核酸の調製には、何ら限定するものではないが、例えば、市販製品であるThruPLEX(登録商標)Tag−Seq Kit(タカラバイオ社製)を利用することができる。また解析したい核酸領域があらかじめ決まっている場合、当該領域を増幅可能なPCRプライマーと分子バーコードを連結させた核酸を用いることにより、クローナル核酸への分子バーコードの付加と増幅を同時に行うことができる。あるいは、いわゆるプローブキャプチャー法などにより、目的領域を包含するテンプレート核酸を濃縮しても良い。該濃縮操作はコピーを得た後に行っても良い。
該分子バーコートはさらに、「ステム(stem)」などと呼ばれる配列を含んでいてもよい。当該ステム配列は、分子バーコードの開始点の目印となる配列である。その配列長や配列パターンは特に限定されないが、例えば、8〜10塩基長の配列で使用される。当該ステム配列は、クローナル核酸と分子バーコードに挟まれる位置に配置される。
クローナル核酸への分子バーコードの付加後、得られたテンプレート核酸を増幅してライブラリーを作製する(S3)。本実施の形態では、シーケンサー50が次世代シーケンサーの場合、配列解析を行うためにはある程度の核酸量を必要とするため、低頻度な変異を含む核酸分子も核酸量を増やす、すなわちコピーを得る必要がある。
テンプレート核酸のコピーの製造法は特に限定しないが、核酸増幅法が挙げられ、操作の簡便さなどからPCRに基づく方法によって行うことが好ましい。このようにして得られたシーケンサーによる解読対象となるサンプルはライブラリーと呼ばれる。
このライブラリーのシーケンシングにより得られた塩基配列データを用いてDNA配列の変異が解析される。本実施の形態では、テンプレート核酸の増幅は、PCRに基づく方法によって行うが、これに限定されない。ライブラリーの作製方法は特に限定されない。ライブラリー作製過程において、PCR法やプローブキャプチャー法により解析対象とする領域を濃縮する場合もあるが、実施において濃縮手法に関する制限はない。
ライブラリーをシーケンサー50へ入力し、塩基配列を解析する(S4)。すなわち、リード配列は、公知のシーケンシング法による塩基配列解読法により取得することができる。シーケンサーのプラットホームは特に限定されないが、解析対象がゲノムDNAのように情報量が膨大になる場合には、次世代シーケンサーが望ましい。次世代シーケンサーとして、例えば、HiSeq(登録商標)及びMiSeq(登録商標)(illumina社製)、Ion ProtonTMおよびIon PGMTM(Thermo Fisher Scientific社製)、PacBio(登録商標) RS IIおよびPacBio(登録商標) SequelTM(Pacific Biosciences社製)、MinION、GridION X5、PromethION、およびSmidgION(Oxford Nanopore Technologies社製)などが挙げられ、数千万から数億のDNA断片を同時並行的に処理して塩基配列を決定する装置が好ましい。
シーケンサー50での解析が終了すると、シーケンサー50からの解析結果を用いてDNA配列における変異が検出される(S5)。この工程(S5)は変異検出装置10により実施される。
シーケンサー50から得られる塩基配列データ(リード配列)のファイル形式はFASTQ形式であるが、他の形式でもよい。FASTQ形式は、当該技術において公知の形式であり、図3(A)に示すような形式となる。データの各列の意味は図3(B)に示すとおりである。FASTQ形式のシーケンスデータでは、各塩基のベースコール(塩基の指定)に対する正確さをあらわす指標として、Phredクオリティスコア(Phred quality score)が出力される。本実施の形態では、一塩基あたりのシーケンスによるエラー率をエラー頻度(Perror)で示す。Phredクオリティスコア(a)とエラー頻度(Perror)の関係は次式で表される。
Perror = 10−a/10 (/塩基)
a = ASCII−33
ASCIIコードの値が、FASTQ形式のデータにおいて出力される。
Perror = 10−a/10 (/塩基)
a = ASCII−33
ASCIIコードの値が、FASTQ形式のデータにおいて出力される。
Phredクオリティスコアは、次世代シーケンサーにより自動的に算出される。エラー頻度への変換式はそれぞれのプラットホームによって異なるが、本実施の形態においては限定される要因はない。
各核酸分子由来のリード配列を取り込む工程において、低品質なリード配列、すなわち、ベースコールの精度(Phredクオリティスコア)が低い塩基を含む配列データ自体を取り除く、いわゆるリードフィルタリング(リードクリーニングともいう)処理を行っても良い。この処理には公知の技術を利用することができ、例えば公知のソフトウエアであるCutadapt(ドルトムント工科大学)やTrimmomatic(アーヘン工科大学)などが挙げられる。
また、リード配列の一部分の塩基が低品質の場合、該当する一部分の塩基配列のみを取り除く、いわゆるトリミング(trimming)処理を行ってもよい。トリミング処理は、マスキング(masking)処理と呼ばれることもある。該処理には公知の技術を利用することができ、例えば公知のソフトウエアであるTrimmomaticやsickle(https://github.com/najoshi/sickle)などが挙げられる。
図4(a)(b)は、図2の変異解析処理(S5)のより具体的な手順を示した図である。後で詳しく説明するように、これらのうち図4(a)は、PCRエラーやシーケンスのエラーを考慮したデータ解析による手順を示したものであり、図4(b)は、リードのアライメントを考慮したデータ解析による手順を示したものである。両方の図から明らかなように、図4(b)は、「ローカルアセンブリ」(S17)の工程を含むことのみ、図4(a)と異なる。以下、図4(a)と図4(b)の両方を参照しつつ、変異解析処理(S5)の各工程について具体的に説明する。
1)リード配列の取得(S11)
まず、シーケンサー50から、リード配列データを取得する(S11)。
まず、シーケンサー50から、リード配列データを取得する(S11)。
2)参照配列へのマッピング(S12)
リード配列はクローナル核酸を構成する各核酸分子由来の様々な断片由来のコピー分子の配列情報の集合体であるため、先ず、参照配列との比較によりリード配列を分類するマッピング処理を行う。ここで、参照配列とは、基準とする配列であり、任意の配列を用いることができる。本発明を何ら限定するものではないが、例えば、NCBI GenBank、DDBJ、UCSC Genome Browserなどに登録されている配列を用いることができる。例えば、解析対象がヒトゲノムDNAの場合、Genome Reference Consortium human build 38(GRCh38)やUCSC human genome 19(hg19)を参照配列とすることができる(これらのバージョンは随時更新されており、必要に応じて適切なバージョンを選択すればよい)。また、参照配列とする領域は、目的に応じて適宜選択することができる。レコードに収録されている配列全長でもよいし、目的とする任意の領域のみであってもよい。該参照配列と後述のコンセンサスリードを比較し、参照配列とは異なる塩基がコンセンサスリードにあった場合、当該塩基は変異の候補となる。
リード配列はクローナル核酸を構成する各核酸分子由来の様々な断片由来のコピー分子の配列情報の集合体であるため、先ず、参照配列との比較によりリード配列を分類するマッピング処理を行う。ここで、参照配列とは、基準とする配列であり、任意の配列を用いることができる。本発明を何ら限定するものではないが、例えば、NCBI GenBank、DDBJ、UCSC Genome Browserなどに登録されている配列を用いることができる。例えば、解析対象がヒトゲノムDNAの場合、Genome Reference Consortium human build 38(GRCh38)やUCSC human genome 19(hg19)を参照配列とすることができる(これらのバージョンは随時更新されており、必要に応じて適切なバージョンを選択すればよい)。また、参照配列とする領域は、目的に応じて適宜選択することができる。レコードに収録されている配列全長でもよいし、目的とする任意の領域のみであってもよい。該参照配列と後述のコンセンサスリードを比較し、参照配列とは異なる塩基がコンセンサスリードにあった場合、当該塩基は変異の候補となる。
マッピング技術において、ゲノム配列がすでに解読されているヒトを含む生物種については、そのゲノム配列を参照配列として一般にFASTA形式にて取得可能である。本開示では、テンプレート核酸から解析したリード配列を参照配列にマッピングした後に、比較することで、変異の有無を判定することが可能であり、次世代シーケンスによる解析では、多量のリードシーケンスから変異の有無を検出できる。
本実施の形態において、上記で取得した参照配列に対してリード配列をマッピングすることになる。特に限定するものではないが、たとえば、BWA(Wellcome Trust Sanger Institute)、Bowtie(ジョンズ・ホプキンズ大学)、DNASTAR(登録商標) SeqMan NGen(DNASTAR社製)、Hisat2(ジョンズ・ホプキンズ大学)、Novoalign(NOVACRAFT社製)等のソフトがマッピング作業に使用できる。
本実施の形態において、マッピングの結果は、図5に示すようなSAM/BAM形式ファイルと呼ばれる公知のファイル形式で出力される。BAM形式ファイルはSAM形式ファイルのバイナリーファイルとなる。マッピングの結果を出力するファイル形式はBAM形式に限定されるものではない。図5(A)は、参照配列へのマッピング(S12)の結果が出力されたSAM/BAM形式ファイルであり、図5(B)は、分子バーコードのクラスタリング(S13)の結果が追加されたSAM/BAM形式ファイルである。SAM/BAM形式データの各列の意味は図5(C)に示すとおりである。
このようにマッピングされたリード配列に対し、次に個々のテンプレート核酸に由来するリード配列を分類するクラスタリング処理を行う。テンプレート核酸に前述の分子バーコードを付加した場合は、リード配列中に存在する分子バーコード領域の配列情報を基に分類することが可能である。なお、以下では、分子バーコードを”UMI”と称することがある。
ベイズ推定は解析に使用するリード配列の数が多いほど解析精度が高くなるため、得られたリード配列をいかに解析に使用するか、は重要である。
既存の解析手法において、Phredクオリティスコア(Phred quality score)が低いリード配列は解析不良の一因となるため、以降の解析に使用されていない。そのため、大量のシーケンスデータを得ながら、解析に使用できるリード配列が少なくなってしまうという問題があった。
既存の解析手法において、Phredクオリティスコア(Phred quality score)が低いリード配列は解析不良の一因となるため、以降の解析に使用されていない。そのため、大量のシーケンスデータを得ながら、解析に使用できるリード配列が少なくなってしまうという問題があった。
一方、本開示の方法では、Phredクオリティスコア(Phred quality score)が低い箇所の塩基を「N」(塩基種の判別不可な塩基)とし、参照配列と比較する際には「完全一致した塩基」とみなすことによって以降の解析に使用可能とした。
マッピングによって1つのグループとなったとしても、異なる分子バーコードのリード配列が混在していることもある。その場合、分子バーコード配列に基づいたグループ分けを行うことにより、正しいグループへ再編することができる。
既存の解析手法において、マッピングされなかったリード配列は解析不良の一因となるため、以降の解析に使用されていない。そのため、大量のシーケンスデータを得ながら、解析に使用できるリード配列が少なくなってしまうという問題があった。
しかしながら、本発明において、上記マッピングから外れた(マッピングされなかった、ummapped)リード配列は判断保留とし、分子バーコード配列に基づいたグループ分け工程において、作成されたファミリーと配列比較して解析に使用できるか判断する。このようにして、解析に使用可能となるリード配列を増やしている(解析不可とするリード配列を減らしている)。
3)分子バーコード(UMI)のクラスタリング(S13)
本実施の形態では、SAM/BAM形式ファイルを介して、リード配列のマッピング情報(e.g.リード配列、UMI配列、参照配列に対するマッピング位置)を含む一連の情報を取得した後に、各リード配列に付与されたUMI配列を決定して、同じクローナル核酸に由来するリード配列群を一つに取りまとめる、UMIの分類工程(分子バーコードのクラスタリング)を行う。図6(A)は、リード配列を、参照配列上でのマッピング開始点と終了点によりグループ分けした後の状態の一例を示している。図6(A)は、参照配列上でのマッピング位置に基づき分類された3つのグループ(Position family)を例示している。図6(B)は、その後さらに、UMI(分子バーコード)の配列によりグループ分けした状態の一例を示している。図6(B)は、UMI(分子バーコード)の配列による分類された2つのグループ(UMI family)を例示している。
本実施の形態では、SAM/BAM形式ファイルを介して、リード配列のマッピング情報(e.g.リード配列、UMI配列、参照配列に対するマッピング位置)を含む一連の情報を取得した後に、各リード配列に付与されたUMI配列を決定して、同じクローナル核酸に由来するリード配列群を一つに取りまとめる、UMIの分類工程(分子バーコードのクラスタリング)を行う。図6(A)は、リード配列を、参照配列上でのマッピング開始点と終了点によりグループ分けした後の状態の一例を示している。図6(A)は、参照配列上でのマッピング位置に基づき分類された3つのグループ(Position family)を例示している。図6(B)は、その後さらに、UMI(分子バーコード)の配列によりグループ分けした状態の一例を示している。図6(B)は、UMI(分子バーコード)の配列による分類された2つのグループ(UMI family)を例示している。
UMIの分類の工程(分子バーコードの配列によりグループ分け)では、以下の工程a〜dを行う。
工程a:参照配列に対するマッピング位置からなる各リード配列の位置に基づき、UMI配列に非依存的にリード配列群を分類する。分類されたグループはPosition familyとして管理する。すなわち、SAM/BAM形式ファイルに記載されている参照配列上におけるマッピング開始点と終了点が同一のリード配列群をPosition familyとして管理する。
工程b:各Position family内のリード配列数を集計する。集計値が閾値より少ないPosition familyにはフラグを立てる。
工程c:各Position family内でUMI配列の類似度に基づきリード配列を分類する。UMI配列の類似度に基づき分類したグループをUMI familyと称する。また、二種類以上のUMIが同一ライブラリー上に存在する場合には、全てのUMI配列を結合させる。具体的には、図6の例では、L−UMIとR−UMIを結合させる。UMI配列の類似度に基づきPosition family内のリード配列群をさらに分類する。類似度の閾値を設定して、ある程度のミスマッチを許容することで、閾値内のミスマッチのUMI同士は同じ配列とする。Position family内のUMI配列の出現数によって、全リード配列数に対して出現数の少ないUMI配列を持つリード配列にはSAM/BAMファイル上でフラグを立てる。
工程d:Position family内のリード配列数の多寡に基づき、近傍に位置するPosition familyに限定して、最終UMI familyの再分類を実施する。すなわち、フラグを立てた、リード数が少ないPosition familyについて、近傍のPosition familyに統合できるかをUMI配列の類似度に基づき再判定する。統合できると判定した場合、フラグを立てた、リード数が少ないPosition familyを近傍のPosition familyに統合する。
工程b:各Position family内のリード配列数を集計する。集計値が閾値より少ないPosition familyにはフラグを立てる。
工程c:各Position family内でUMI配列の類似度に基づきリード配列を分類する。UMI配列の類似度に基づき分類したグループをUMI familyと称する。また、二種類以上のUMIが同一ライブラリー上に存在する場合には、全てのUMI配列を結合させる。具体的には、図6の例では、L−UMIとR−UMIを結合させる。UMI配列の類似度に基づきPosition family内のリード配列群をさらに分類する。類似度の閾値を設定して、ある程度のミスマッチを許容することで、閾値内のミスマッチのUMI同士は同じ配列とする。Position family内のUMI配列の出現数によって、全リード配列数に対して出現数の少ないUMI配列を持つリード配列にはSAM/BAMファイル上でフラグを立てる。
工程d:Position family内のリード配列数の多寡に基づき、近傍に位置するPosition familyに限定して、最終UMI familyの再分類を実施する。すなわち、フラグを立てた、リード数が少ないPosition familyについて、近傍のPosition familyに統合できるかをUMI配列の類似度に基づき再判定する。統合できると判定した場合、フラグを立てた、リード数が少ないPosition familyを近傍のPosition familyに統合する。
本実施形態において、UMI配列が同一かどうかの配列類似度の判定はハミング距離などの編集距離を使用する。しかし、配列類似度の判定はこれに限定するものではなく、例えば、モチーフ出現頻度に基づいた(類似度)ネットワーク解析を想定したZscoreベースのベクトルの類似度や分子系統解析で用いる系統・進化距離といった指標も選択可能である。さらに類似度の判定基準は、UMI配列の長さやテンプレート核酸のコピー数によって変更してもよい。
次世代シーケンスの作業では、UMI配列自体にもPCR増幅やシーケンスによるエラーが生じる可能性がある。そのため本実施形態では、このエラー発生の可能性に対して、下記の工程が実施可能である。
工程a:シーケンスデータのPhredクオリティスコアが低い塩基については、N(任意の塩基を示す記号)でマスクして、配列類似度の比較対象から除外する。
工程b:UMI配列間のミスマッチを許容するように編集距離の閾値を調整する。すなわち、UMI配列の類似度(ハミング距離等)に基づき、Position family内のリード配列群をさらに分類する。類似度の閾値を設定して、ある程度のミスマッチを許容することで、閾値の条件を満たすミスマッチのUMI同士は同じ配列とする。
工程a:シーケンスデータのPhredクオリティスコアが低い塩基については、N(任意の塩基を示す記号)でマスクして、配列類似度の比較対象から除外する。
工程b:UMI配列間のミスマッチを許容するように編集距離の閾値を調整する。すなわち、UMI配列の類似度(ハミング距離等)に基づき、Position family内のリード配列群をさらに分類する。類似度の閾値を設定して、ある程度のミスマッチを許容することで、閾値の条件を満たすミスマッチのUMI同士は同じ配列とする。
さらに国際公開第2016/176091号では、PCR増幅及びシーケンス過程において、最初に付与されたUMI配列とは異なる配列に置き換わるUMI置換現象が示唆されている(e.g. UMI Jumping、PCR Jumping [UMI−tools: modeling sequencing errors in Unique Molecular Identifiers to improve quantification accuracy, Tom Smith et al. http://dx.doi.org/ 10.1101/gr.209601.116.]、index swapping [Characterization and remediation of sample index swaps by non−redundant dual indexing on massively parallel sequencing platforms, Maura Costello et al.(http://dx.doi.org/10.1101/200790)]、Characterization and remediation of sample index swaps by non−redundant dual indexing on massively parallel sequencing platforms, Maura Costello et al.(http://dx.doi.org/10.1101/200790)にて報告/示唆済み)。この可能性に対処するため、UMI配列が複数あるようなライブラリー構成の場合には、Position family内の範囲において、少なくとも一つのUMIが完全一致すれば同一のUMIを持つリード配列として分類する工程を含めている。すなわち、複数のUMIが存在するライブラリー構成の場合には、Position family内において、少なくとも一つのUMIが完全に一致すれば同一のUMIを持つリード配列としてまとめる。
本開示を何ら限定するものではないが、例えば、分子バーコードが6塩基からなる場合、一方または両方の分子バーコードにおいて2塩基以下のミスマッチ、または両末端いずれかの分子バーコード領域のリード配列6塩基が一致すればクラスタリングの対象となるリード配列とみなす。これにより、UMI配列自体に存在するエラーの可能性を考慮したクラスタリングが可能となる。
ここで、該UMI配列自体に存在するエラーとは、ミスマッチ塩基、ベースコール不良塩基、マッピングエラーなどを指す。既存の解析手法ではミスマッチ塩基のみ考慮されていた。一方、本発明ではベースコール不良塩基やマッピングエラーなども考慮することにより、解析精度を向上することができた。
本発明を何ら限定するものではないが、分子バーコード配列自体に存在するエラーの存在率は、分子バーコード塩基長の50%以下、好ましくは40%以下、さらに好ましくは35%以下、特に好ましくは30%以下であれば、該エラーの可能性を考慮したクラスタリングが可能である。
分子バーコードがクローナル核酸の両末端に付加される場合、該UMI配列自体に存在するエラーが一方の分子バーコードのみに存在しても、両方の分子バーコードに存在しても、本発明の方法では、当該エラーの可能性を考慮したクラスタリングが可能である。
本実施の形態において、UMIの分類結果は前述のようにSAM/BAM形式ファイルのデータに追記する。すなわち、図5(B)のSAM/BAM形式(出力形式)に示すように、一列目のリード配列の名称に[染色体番号]−[マッピング開始点]−[マッピング終了点]−[アライメントパターン]−[分類後のUMI配列]−[分類前のUMI配列]という形式としている。さらにBC:Z:sequence:sequence,XD:i:count,XP:i:statusに必要な情報を追加したものをSAM/BAM形式ファイルのオプション情報として出力する(以下の表1参照)。追記の形式はこの形式に限定されるものではない。
本実施の形態では、リードのアライメントパターンのステータス(XP:status)情報は、上述のように、0、1、2、3の値を採り、このうちXP=0の場合が「正常」のステータスを表すものとされている。また、XP=1、2、3のときのデータを取り扱うことにより、別染色体にマッピングされるリードの場合や、片方のリードが参照配列上にマッピングされない場合なども、UMIの分類結果として出力することができる。
本実施の形態において、UMI分類結果を追加したSAM/BAM形式ファイルに基づき変異検出を実施する。
4)ローカルアセンブリ(S17)
このローカルアセンブリの工程は、図4(b)に示す変異解析処理に含まれる工程である。ローカルアセンブリの工程では、参照配列にはない検体特有のリード配列長以上の配列領域を決定する処理が行われる。この配列と参照配列を比較することで、リード配列長以上の配列サイズを有する塩基置換や挿入や欠損が確認できる。なお本実施の形態のローカルアセンブリは、Rimmer et al. (Nature Genetics, 2014) に掲載されたPlatypusのアルゴリズムを参考にコーディングされており、参照配列とリード配列よりtwo−colored de Bruijn Graphを構築してアセンブル配列を構築する。なおローカルアセンブリにて使用可能なオプションおよびその例を、以下の表2に示す。
このローカルアセンブリの工程は、図4(b)に示す変異解析処理に含まれる工程である。ローカルアセンブリの工程では、参照配列にはない検体特有のリード配列長以上の配列領域を決定する処理が行われる。この配列と参照配列を比較することで、リード配列長以上の配列サイズを有する塩基置換や挿入や欠損が確認できる。なお本実施の形態のローカルアセンブリは、Rimmer et al. (Nature Genetics, 2014) に掲載されたPlatypusのアルゴリズムを参考にコーディングされており、参照配列とリード配列よりtwo−colored de Bruijn Graphを構築してアセンブル配列を構築する。なおローカルアセンブリにて使用可能なオプションおよびその例を、以下の表2に示す。
de Bruijn graphでは、まず塩基配列を長さkの塩基配列(k−mer)に切断して、Node(頂点)とする。さらに各Node(頂点)をEdge(辺)で結ぶ。複数のリード配列もしくは参照配列にてk−merのNode−Edgeを発生させて、同じk−merを有するNodeや同じNodeを合体させたグラフがde Bruijn graphとなる(図7参照)。本実施の形態では、de Bruijn graphのNodeとEdgeに対して、塩基配列の由来(参照配列もしくはリード配列から由来、分子バーコードの有無およびその配列情報)や信頼性スコア(シーケンスのPhredクオリティスコア)から算出される重みを付けることができる。
(4−1)de Bruijn graph(デブルーイングラフ)の作成
本実施の形態において、de Bruijn graphは解析遺伝子の一部領域を対象に、配列長(表2のwindowサイズオプションにより変更可能)を指定した対象領域に相当する参照配列と同領域にマッピングされたリード配列を使用して作成する。さらに使用するリード配列には、同領域にマッピングされたリード配列自体とペアリード(つまり、参照配列にマッピングされなかったunmappedリードまたは対象領域外にマップされたリード)も含まれる。また、本実施の形態では、グラフのEdge(辺)のcolorを、参照配列側を0、リード配列側を1として、両方が存在する部分は0.5と設定してEdge(辺)の由来を区別することが可能である。
本実施の形態において、de Bruijn graphは解析遺伝子の一部領域を対象に、配列長(表2のwindowサイズオプションにより変更可能)を指定した対象領域に相当する参照配列と同領域にマッピングされたリード配列を使用して作成する。さらに使用するリード配列には、同領域にマッピングされたリード配列自体とペアリード(つまり、参照配列にマッピングされなかったunmappedリードまたは対象領域外にマップされたリード)も含まれる。また、本実施の形態では、グラフのEdge(辺)のcolorを、参照配列側を0、リード配列側を1として、両方が存在する部分は0.5と設定してEdge(辺)の由来を区別することが可能である。
(4−2)グラフの再構築
本実施の形態において、(例えば、図7の場合、)構築したグラフの5’末端から探索して、環状化していた場合には、k−merの長さを伸長して、再度グラフを作成することができる。
本実施の形態において、(例えば、図7の場合、)構築したグラフの5’末端から探索して、環状化していた場合には、k−merの長さを伸長して、再度グラフを作成することができる。
(4−3)最短経路の検出
本実施の形態において、グラフ上の経路探索に、辺の重みが非負数の場合の単一始点最短経路問題を解くための最良優先探索によるアルゴリズムであるダイクストラ法(優先度付きキュー)を採用する (https://ja.wikipedia.org/wiki/ダイクストラ法 または https://en.wikipedia.org/wiki/Dijkstra%27s_algorithm)。つまり、個々の探索経路がアセンブル配列となり、さらに経路上のEdge(辺)のcolorの違いにより、参照配列特有もしくはリード配列特有のEdgeが判別可能であるため、変異候補を同定することができる。しかしながら、経路探索の方法やアルゴリズムは、上記に限定されるものではない。
本実施の形態において、グラフ上の経路探索に、辺の重みが非負数の場合の単一始点最短経路問題を解くための最良優先探索によるアルゴリズムであるダイクストラ法(優先度付きキュー)を採用する (https://ja.wikipedia.org/wiki/ダイクストラ法 または https://en.wikipedia.org/wiki/Dijkstra%27s_algorithm)。つまり、個々の探索経路がアセンブル配列となり、さらに経路上のEdge(辺)のcolorの違いにより、参照配列特有もしくはリード配列特有のEdgeが判別可能であるため、変異候補を同定することができる。しかしながら、経路探索の方法やアルゴリズムは、上記に限定されるものではない。
5)コンセンサスリード(S14)
上述の、特に「3)分子バーコード(UMI)のクラスタリング(S13)」までで述べたように、UMI配列に基づき、PCRエラーやシーケンスのエラーを考慮して、同じUMI配列(同じUMI family)をもつリード配列から共通のリード配列すなわちコンセンサスリードを出力する。当該コンセンサスリードは、単にコンセンサス、コンセンサス配列、またはコンセンサスリード配列と呼ぶこともある。図8(A)は、コンセンサスリードを求めるコンセンサスリード処理の処理工程を説明した図である。(以下、「(a)PCRエラーやシーケンスエラーを考慮したデータ解析」参照)。
一方、ローカルアセンブリを用いた解析では、シーケンスのエラーとリードのマッピング距離を考慮して、上述の、特に「4)ローカルアセンブリ(S17)」にて述べたように、同じUMI配列(同じUMI family)をもつリード配列からアセンブル配列を構築し、上述の「(4−3)最短経路の検出」に基づいて、その中から信頼性が高いアセンブル配列をコンセンサスリードとする。図9(A) は、コンセンサスリードを求める処理工程を説明した図である。(後の「(b)リードのアライメントを考慮したデータ解析(Large Indelの検出)」参照)。
更に、「(a)PCRエラーやシーケンスエラーを考慮したデータ解析」には、別手法もある(後の「(a’)PCRエラーやシーケンスエラーを考慮したデータ解析の別手法」参照)。
上述の、特に「3)分子バーコード(UMI)のクラスタリング(S13)」までで述べたように、UMI配列に基づき、PCRエラーやシーケンスのエラーを考慮して、同じUMI配列(同じUMI family)をもつリード配列から共通のリード配列すなわちコンセンサスリードを出力する。当該コンセンサスリードは、単にコンセンサス、コンセンサス配列、またはコンセンサスリード配列と呼ぶこともある。図8(A)は、コンセンサスリードを求めるコンセンサスリード処理の処理工程を説明した図である。(以下、「(a)PCRエラーやシーケンスエラーを考慮したデータ解析」参照)。
一方、ローカルアセンブリを用いた解析では、シーケンスのエラーとリードのマッピング距離を考慮して、上述の、特に「4)ローカルアセンブリ(S17)」にて述べたように、同じUMI配列(同じUMI family)をもつリード配列からアセンブル配列を構築し、上述の「(4−3)最短経路の検出」に基づいて、その中から信頼性が高いアセンブル配列をコンセンサスリードとする。図9(A) は、コンセンサスリードを求める処理工程を説明した図である。(後の「(b)リードのアライメントを考慮したデータ解析(Large Indelの検出)」参照)。
更に、「(a)PCRエラーやシーケンスエラーを考慮したデータ解析」には、別手法もある(後の「(a’)PCRエラーやシーケンスエラーを考慮したデータ解析の別手法」参照)。
(a)PCRエラーやシーケンスエラーを考慮したデータ解析
本実施の形態では、まず、UMI配列に基づき、PCRエラーやシーケンスエラーを考慮して、同じUMI配列を持つ(すなわち、同じUMI family内の)リード配列から共通のリード配列(コンセンサスリード)を求める。
本実施の形態では、まず、UMI配列に基づき、PCRエラーやシーケンスエラーを考慮して、同じUMI配列を持つ(すなわち、同じUMI family内の)リード配列から共通のリード配列(コンセンサスリード)を求める。
なお、特定の塩基種が配列解析時のエラーによるものかクローナル核酸中に存在する塩基に由来するものか判断するために、コンセンサスリード及び当該配列中の各塩基に対する信頼性スコア(UQ)を使用する。当該信頼性スコア(UQ)は各UMI family内で各塩基に対して算出され、信頼性スコア(UQ)は次式で表される。
信頼性スコア(UQ)を算出する際、単純な多数決ではなく、統計処理を行う。コンセンサスリードは、市販製品、例えばConnor(https://github.com/umich−brcf−bioinf/Connor)を利用して得ることができる。
コンセンサスリード処理の開始前に、以下の工程a、bを行う。
工程a)各position(No = 1,2,3…k)にマッピングされたUMI family(No = 1,2,3…j)を抽出、
工程b)各UMI family内の全リード配列に対して、番号(No = 1,2,3…i)を割り振り、該当position上の塩基種類、各塩基のリード配列数、及びシーケンスのエラー頻度(Perror)のデータを取得する。
工程a)各position(No = 1,2,3…k)にマッピングされたUMI family(No = 1,2,3…j)を抽出、
工程b)各UMI family内の全リード配列に対して、番号(No = 1,2,3…i)を割り振り、該当position上の塩基種類、各塩基のリード配列数、及びシーケンスのエラー頻度(Perror)のデータを取得する。
これらの処理は、UMI分類後のSAM/BAMファイルのデータに対して実施する。図8(B)に示す例の場合、第3染色体の38,930番目の塩基対(bp)にマッピングされた特定の、すなわちj番目のUMI family (No=j)では、リード配列が合計7本(リード配列数No=1…7)存在する。第3染色体の38,930番目の位置の塩基がAとなるリード配列は4本であり、Cとなるリード配列は3本である。それぞれのエラー頻度(Perror)は、SAM/BAM形式ファイルの11列目のASCII情報から入手できる(図5(B),(C)参照)。
各UMI family内で各塩基に対する信頼性スコア(UQ)は、特に限定されるものではないが、例えば、ベイズアプローチでは事後分布として推定することができる。具体的には各UMI familyにおいて、該当する場所での対象塩基(each base)毎の確度を推定する。対象塩基は全体集合θの部分集合であり、全4種類の塩基対(A,T,G,C)と、同参照配列上にマッピングされる全UMI familyで検出されるすべての挿入または欠損塩基(other obs. alleles)とから構成される。図8(B)の例では、UMI family(No=j)の該当箇所(chr3:38,930bp)において、A, Cの二種類の塩基が検出されているが、集合θの構成はA,T,G,Cの4種類となる(挿入と欠損塩基は検出されないため含まれない)。
公知のソフトウエアSmCounterの統計モデルは、挿入または欠損塩基に対して検出感度が低い仕様となっている。一方、本発明の統計モデルでは、θにおいて挿入または欠損塩基も考慮することにより、塩基置換に加え、塩基挿入および塩基欠損についても検出感度が高くなった。
信頼性スコア(UQ)の式において、P(UMIj|each base)は尤度関数であり、P(each base)は対象塩基の出現率を示す事前分布である。事前分布は1/θのような一様分布を採用することができる。しかし、事前データ等からθもしくはP(each base)に関する何らかの情報があれば変更することが可能である。
尤度関数P(UMIj|each base)は、シーケンスエラーとPCR増幅エラーを考慮した上で、対象塩基の事前分布から各UMI familyの構成するリード配列データが生じる確率を表現する。例えば、図8(B)に示す例で検討すると、対象塩基がAの場合(each base = A)の尤度関数は、UMI familyでは本来Aが正しいとした場合に、対象塩基がAとなるリード配列が4本、Cとなるリード配列が3本となる7本のリード配列が生じる尤度P(UMIj|A)を計算する。
尤度関数では、シーケンスエラー(Serror)とPCR増幅エラーを表現している。尤度関数の式において、第一項がシーケンスエラー(Serror)の頻度を示し、第二項がPCR増幅エラーの発生頻度を示す。図8(B)の例において、対象塩基がAの場合(each base = A)、対象塩基がCとなる3本のリード配列では、シーケンスエラー(A→Cへのエラー)がSerrorの頻度で生じる。一方、対象塩基がAとなる4本のリード配列では(1−Serror)の頻度で正しくシーケンスされたことになり、この条件での尤度を求めることになる。
なお、第一項では、PCR増幅エラーが生じない場合を想定しているため、corCp(PCR増幅エラーが無しとなる確率)を乗算している。さらに、nyは集合θの各構成塩基(Y)のリード配列の数(冗長度)を示している。図8(B)の例では、Y=Aの場合にはnyは4本となる。一方、集合θの内で実際に検出できなかったT及びGについては、nyは0本となり、この第一項の計算は実施されないこととなる。すなわち、PCR増幅エラーによる尤度(第二項)のみが計算される。
本開示の方法において、シーケンスエラー(Serror)は上述のPhredクオリティスコアから算出されるPerrorに基づき算出している。しかし、挿入・欠損塩基がシーケンスエラーとして生じる可能性は一塩基置換(single nucleotide variant、SNV)より低くなるという報告もある(Marinier E等、BMC Bioinformatics. 2015 Jan 16;16:10. doi: 10.1186/s12859−014−0435−6; Schirmer M等、Nucleic Acids Res. 2015 Mar 31;43(6):e37. doi: 10.1093/nar/gku1341参照)。そのため本実施の形態では、対象塩基(each base)が挿入・欠損塩基の場合には、Perrorに補正値m(1/100)を乗算してエラー率を低く抑えるよう処理を加えて、検出性能を向上させている。
PCRerror(x)は、xサイクルのPCRサイクルで一つの変異が生じる確率と、xサイクルのPCR産物が最終サイクル(n)後のPCR産物に含まれる確率との積により求められる。
前者の変異出現の確率はポアソン分布に従う。ポアソン分布のパラメーターはPCR酵素のエラー出現率(μ)と、PCR酵素の増幅効率(λ)とにより期待値として求めることができる。後者の確率は二項分布の期待値(最大尤度)として算出する。よってPCRerror(x)は下記式で求められる。
PCRerror(x)の算出において、変異が出現するPCRサイクル数xを推定することが重要であるが、本実施の形態では、この推定をUMI family内の各塩基のリード配列数に基づき実施する。図8(B)の例では、Aが対象塩基の場合ではCがPCR増幅エラーによって出現したケースに該当することから、その割合は3/7(43.8%)となる。このリード配列数の割合は、全PCR産物中に含まれるPCR増幅エラーに由来するPCR産物の割合の下記数式で近似的に示される。
Sm(1+λm)n−x/So(1+λw)n
ここで、Smは、PCR増幅エラーの初発核酸コピー数であり、1である。Soは、最初の初発核酸コピー数であり、1である。λmはPCRエラーの核酸産物に対するPCR増幅効率である。λwはPCRエラー無しの核酸産物に対するPCR増幅効率である。
Sm(1+λm)n−x/So(1+λw)n
ここで、Smは、PCR増幅エラーの初発核酸コピー数であり、1である。Soは、最初の初発核酸コピー数であり、1である。λmはPCRエラーの核酸産物に対するPCR増幅効率である。λwはPCRエラー無しの核酸産物に対するPCR増幅効率である。
尤度関数P(UMIj|each base)の第二項では、PCR酵素のエラー出現率(μ)とPCR酵素の増幅効率(λ)が重要な設定パラメーターとなる。
本実施の形態では、尤度関数では、PCR増幅エラーとシーケンスエラーという二種類のみを考慮し、互いに排反なイベントとして扱っているが、相互作用の効果(PCR増幅エラーとシーケンスエラーの両方が生じるケースが該当)やその他の要因も式に追加することができる。その他の要因には、マッピングエラー(異なる箇所と違う場所にリード配列がマッピングされるケース)や酸化反応による塩基置換などが含まれる。
上記の工程を経て、各UMI family内で、各位置の塩基に対して信頼性スコア(UQ)が算出される。図8(B)のNo.jのUMI familyを例にすると、38,930番目の塩基に対して、A,T,G,Cのそれぞれに対する信頼性スコアが算出されることとなり、その中で最も信頼性スコアが高い塩基をコンセンサスリードの配列とする。図8(B)の例では、Aに対する信頼性スコアが最も高くなるため、Aをコンセンサスリードの配列とする。本手法はコンセンサスリードに対して信頼性スコアを算出する点が従来法と異なる。
(b)リードのアライメントを考慮したデータ解析(Large Indelの検出)
本実施の形態では、ローカルアセンブリを用いる解析においては、シーケンスのエラーとリードのマッピング距離を考慮して、同じUMI配列(同じUMI family)を持つリード配列から共通のリード配列、即ちコンセンサスリードを求める。図9(A)は、コンセンサスリードを求める処理工程を示した図である。
本実施の形態では、ローカルアセンブリを用いる解析においては、シーケンスのエラーとリードのマッピング距離を考慮して、同じUMI配列(同じUMI family)を持つリード配列から共通のリード配列、即ちコンセンサスリードを求める。図9(A)は、コンセンサスリードを求める処理工程を示した図である。
つまり、ローカルアセンブリの最短経路(アセンブル配列、ハプロタイプ候補, H)とそれぞれに由来する変異リストから、10base以上のInsertおよび/またはDeletionを持つ経路を対象に下記の手順に準じて、最大事後確率P{each haplotype|UMIj}を決定し、(後で説明する)変異検出に繋げる。本実施の形態における最大事後確率Pは、前述(a)の信頼性スコア(UQ)に相当する。
上述のそれぞれの変数の説明を、下記に示す。図9(B)を参照されたい。
[数7]の(2)において、前半部分はリードのマッピング(アライメント)距離の考慮部分を、後半部分はシーケンスのエラーの考慮部分を示す。ここで、前半部分の「P{distk|μ,σ}:Insertサイズ(正規分布μ,σ)に対するリード間の距離の確率」は、同じUMI配列(同じUMI Family、UMIj)を有する全リード配列(k∈ReadPair)をハプロタイプの配列に再アラインメントして、リード間の距離(distk)に対する確率を表すものである。なお、本実施の形態における母集団は、全リードの参照配列上のリード間の距離の平均値(μ)と分散(σ2)によって定義される正規分布を、Insertサイズの分布(正規分布)とする。
[数7]の(2)において、前半部分はリードのマッピング(アライメント)距離の考慮部分を、後半部分はシーケンスのエラーの考慮部分を示す。ここで、前半部分の「P{distk|μ,σ}:Insertサイズ(正規分布μ,σ)に対するリード間の距離の確率」は、同じUMI配列(同じUMI Family、UMIj)を有する全リード配列(k∈ReadPair)をハプロタイプの配列に再アラインメントして、リード間の距離(distk)に対する確率を表すものである。なお、本実施の形態における母集団は、全リードの参照配列上のリード間の距離の平均値(μ)と分散(σ2)によって定義される正規分布を、Insertサイズの分布(正規分布)とする。
(b−1)リード配列の抽出
本実施の形態に使用されるリード配列は、参照配列に対するマッピング情報によって選抜することが可能である。例えば、本実施の形態では、SAM/BAM形式のマッピング結果から取得できるリード配列のCigar(図5(C)参照)にてI/D/S/Hを持つリード(参照配列にマッピングされないunmappedリードも含む)を抽出する。さらに同じUMIを持つリードも抽出する(上記条件に合致しない場合も抽出対象となる)。
本実施の形態に使用されるリード配列は、参照配列に対するマッピング情報によって選抜することが可能である。例えば、本実施の形態では、SAM/BAM形式のマッピング結果から取得できるリード配列のCigar(図5(C)参照)にてI/D/S/Hを持つリード(参照配列にマッピングされないunmappedリードも含む)を抽出する。さらに同じUMIを持つリードも抽出する(上記条件に合致しない場合も抽出対象となる)。
(b−2)再アラインメント
特定UMI(UMIj)の各リードをハプロタイプ(hi)に再アラインメントして、リード間の距離とマッチ/ミスマッチ塩基の情報を入手する。本実施の形態の再アラインメントでは、Smith−Waterman法を採用している。なお、再アラインメントのアルゴリズムは、Smith−Waterman法とは異なるものを採用してもよい。
特定UMI(UMIj)の各リードをハプロタイプ(hi)に再アラインメントして、リード間の距離とマッチ/ミスマッチ塩基の情報を入手する。本実施の形態の再アラインメントでは、Smith−Waterman法を採用している。なお、再アラインメントのアルゴリズムは、Smith−Waterman法とは異なるものを採用してもよい。
(b−3)最大事後確率P{each haplotype|UMIj}の決定 (MAP推定によるハプロタイプ推定)
特定UMI(UMIj)に対して各ハプロタイプの事後確率を求めて、(全ハプロタイプの中で)最大事後確率を示すハプロタイプ(hi)から特定UMI(UMIj)が生み出されたと仮定する。つまり、P{each haplotype|UMIj}は特定UMI(UMIj)が特定ハプロタイプ(hi)に由来する確率を示す。後続の処理では、最大事後確率に基づき、ポアソン分布によって変異検出を実施する。
特定UMI(UMIj)に対して各ハプロタイプの事後確率を求めて、(全ハプロタイプの中で)最大事後確率を示すハプロタイプ(hi)から特定UMI(UMIj)が生み出されたと仮定する。つまり、P{each haplotype|UMIj}は特定UMI(UMIj)が特定ハプロタイプ(hi)に由来する確率を示す。後続の処理では、最大事後確率に基づき、ポアソン分布によって変異検出を実施する。
(a’)PCRエラーやシーケンスエラーを考慮したデータ解析の別手法
更に、上述の「(a)PCRエラーやシーケンスエラーを考慮したデータ解析」の別手法について説明する。
更に、上述の「(a)PCRエラーやシーケンスエラーを考慮したデータ解析」の別手法について説明する。
(a’−1)SNV,MNP,Small Indelsの検出
上述のLarge Indel(上記「(b)リードのアライメントを考慮したデータ解析」参照)と同じく、ベイズ推定によって実施する。P{Ai|UMIk}は、特定UMIがAiというアリルに由来する事後確率を求めることになる(数9)。その際の尤度は、1)PCRエラーが無くベースコールエラーのみ生じる場合、2)PCRエラーが生じてベールコールエラーが無い場合、の二種類の条件から構成されている(数10)。
この数9は、上述の「(a)PCRエラーやシーケンスエラーを考慮したデータ解析」における[数2]と対応するが、PCRエラーのCeの部分が異なる。
上述のLarge Indel(上記「(b)リードのアライメントを考慮したデータ解析」参照)と同じく、ベイズ推定によって実施する。P{Ai|UMIk}は、特定UMIがAiというアリルに由来する事後確率を求めることになる(数9)。その際の尤度は、1)PCRエラーが無くベースコールエラーのみ生じる場合、2)PCRエラーが生じてベールコールエラーが無い場合、の二種類の条件から構成されている(数10)。
(a’−2)リードの抽出
該当Positionにマッピングされるリード配列において、該当箇所に変異をもつリードをCigarとNMタグに基づいて選抜する。更に同じUMIのリードも変異の有無にかかわらず全て抽出する。該当箇所の全変異は集合θで表され、A、T、G、Cと該当箇所の全変異(obs Allele)から構成される要素Aiが属する。
該当Positionにマッピングされるリード配列において、該当箇所に変異をもつリードをCigarとNMタグに基づいて選抜する。更に同じUMIのリードも変異の有無にかかわらず全て抽出する。該当箇所の全変異は集合θで表され、A、T、G、Cと該当箇所の全変異(obs Allele)から構成される要素Aiが属する。
また、それぞれの値の定義と算出方法は以下のものとした。
PCRエラーを持たないPCR産物の割合(noFq):
これはポアソン分布を使った割合である。個々の変数の意味は以下の通りである。
n:PCRサイクル数
λ:PCR増幅効率∈[0,1]
μ:PCRエラー率(PCRの1サイクルで1塩基に対して生じるエラー頻度)
G: テンプレートの長さ(本実施の形態では1bpである)
PCRエラーを持たないPCR産物の割合(noFq):
n:PCRサイクル数
λ:PCR増幅効率∈[0,1]
μ:PCRエラー率(PCRの1サイクルで1塩基に対して生じるエラー頻度)
G: テンプレートの長さ(本実施の形態では1bpである)
(a’−4)最大事後確率P{each base|UMIj}の決定
特定UMI(UMIj)に対して各アリルの事後確率を求めて、(全アリルの中で)最大事後確率を示すアリル(each base)から特定UMI(UMIj)が生み出されたと仮定する。つまり、P{each base|UMIj}は特定UMI(UMIj)が特定アリル(each base)に由来する確率を示す。後続の処理では、最大事後確率に基づき、ポアソン分布によって変異検出を実施する。
特定UMI(UMIj)に対して各アリルの事後確率を求めて、(全アリルの中で)最大事後確率を示すアリル(each base)から特定UMI(UMIj)が生み出されたと仮定する。つまり、P{each base|UMIj}は特定UMI(UMIj)が特定アリル(each base)に由来する確率を示す。後続の処理では、最大事後確率に基づき、ポアソン分布によって変異検出を実施する。
6)エラーモデル構築(S15)
続いて本実施の形態では、全UMI familyから得られたコンセンサスリード及び各塩基に対する信頼性スコア(UQ)を使用して特定位置の塩基種を推定する(図8(A)、図9(A)参照)。
続いて本実施の形態では、全UMI familyから得られたコンセンサスリード及び各塩基に対する信頼性スコア(UQ)を使用して特定位置の塩基種を推定する(図8(A)、図9(A)参照)。
該解析を行うために、まずコンセンサスリードの集合より、特定の位置に存在する特定の塩基種を任意に選抜する。ここで「任意」とは「あらゆる場合」、「すべての場合」という意味である。「特定の位置」とは、参照配列とコンセンサスリードの塩基種が異なる位置という意味である。「特定の塩基種」とは、参照配列の塩基と異なるコンセンサスリードの塩基という意味である。
なお、クローナル核酸において確認したい位置と塩基種については、任意に設定(選択)すればよい。公知の変異塩基の存在を確認したい場合は、その変異に関する部位および塩基種をあらかじめ選択すればよい。また、未知の変異部位を確認したい場合は、参照配列とコンセンサスリードとの比較で参照配列とは異なる塩基種が存在する位置を抽出しても良い。また、塩基種も塩基置換タイプ(変異としてはSNP)の他、挿入部位に相当する塩基種や欠失相当の塩基種(この場合は該等塩基種無との処理)であっても良い。また、参照配列と異なるコンセンサスリードが存在しない位置については、以降の推定工程を行わなくともよい。
次いで、塩基種の選抜に用いたコンセンサスリードの総数に対する、特定の位置に特定の各塩基種を含有するコンセンサスリード数の比率を取得する。
一方、信頼性スコアを用いたパラメトリックなエラーモデルにより、設定された閾値(または有意水準)で特定位置において解析エラーによって特定の塩基が出現する出現率を推定する。
次いで、上記のコンセンサスリード数を基にした比率が解析エラーによる出現率より有意に高い場合に、上記閾値(または有意水準)で特定塩基位置における特定塩基種がクローナル核酸混合物中に存在すると推定する。
塩基種推定の基本モデルはポアソン分布によるパラメトリックなエラーモデルを構築して、任意の箇所で塩基種推定を実施する。実施形態によっては、エラーモデルを変更することが可能である。ポアソン分布の他、例えば、多値変数(事前分布はディリクレ分布)や二値変数でのベータ二項分布での検定、教師付きの学習モデルなども、本発明のパラメトリックなエラーモデルとして使用できる。
塩基種推定は、以下の工程a〜dを含む。
工程a)対象となる箇所を選択する(図8(B)ではchr3:38,930)。
工程b)該当する箇所にマッピングされた全UMI familyを抽出する(図8(B)では、j個のUMI familyを抽出)。
工程c)各UMI familyの該当箇所におけるコンセンサスとUQのスコアを取得する。
工程d)確認したい塩基種を有する参照配列(A)とは異なるコンセンサスリードとUQのスコアのみを選択する。
工程e)工程a〜dを経て得られた情報(選択した塩基種を有するコンセンサスリードのパターンとUQスコア)からエラーモデル(ポアソン分布)を構築する。
工程a)対象となる箇所を選択する(図8(B)ではchr3:38,930)。
工程b)該当する箇所にマッピングされた全UMI familyを抽出する(図8(B)では、j個のUMI familyを抽出)。
工程c)各UMI familyの該当箇所におけるコンセンサスとUQのスコアを取得する。
工程d)確認したい塩基種を有する参照配列(A)とは異なるコンセンサスリードとUQのスコアのみを選択する。
工程e)工程a〜dを経て得られた情報(選択した塩基種を有するコンセンサスリードのパターンとUQスコア)からエラーモデル(ポアソン分布)を構築する。
ポアソン分布のパラメーターλについては、好ましくは以下のように表現できる。
ある位置において、m個のUMI familyが得られ、そのうちn個のUMI familyで塩基X(=A,T,G or C)のコンセンサスリードが得られた場合、下記式によりポアソン分布のパラメーターλの期待値E[λ]が算出できる。
ここで、argminP{X│UMI}は、塩基Xとなったn個のUMI familyのUQ(X)の中の最小値を示す。
例えば、chr3:38,930の箇所で1000個のUMI familyが得られ、そのうち5個のUMI familyでCのコンセンサスリードが得られた場合、下記式によりポアソン分布のパラメーターλが算出できる。
λ=1000*(1−argminP{C|UMI})
ここで、argminP{C|UMI}は、塩基がCとなった5個のコンセンサスリードに対する5個のUQスコアUQ(c)(P{C|UMI})の中の最小値を示す。
λ=1000*(1−argminP{C|UMI})
ここで、argminP{C|UMI}は、塩基がCとなった5個のコンセンサスリードに対する5個のUQスコアUQ(c)(P{C|UMI})の中の最小値を示す。
上記のようにしてλを求めることで、ポアソン分布を用いたエラーモデルを構築できる。
7)変異検出(S16)
上記のようにして構築したエラーモデル(ポアソン分布)を利用して変異を検出する。すなわち、エラーモデルから、m個のUMIを取得した場合に、エラーによってn個のコンセンサスリードが得られる確率が推定できる。
上記のようにして構築したエラーモデル(ポアソン分布)を利用して変異を検出する。すなわち、エラーモデルから、m個のUMIを取得した場合に、エラーによってn個のコンセンサスリードが得られる確率が推定できる。
例えば、1,000個のUMIを取得した場合に、エラーによって5個のコンセンサスリードが得られる確率が推定できる。ここで、前述の閾値は有意水準とも呼ばれ、一般的には0.05(5%)に設定されるが、必要に応じて0.01(1%)や0.001(0.1%)に設定される。そして、エラーモデルから推定した確率を所定の閾値(1%水準又は5%水準等)と比較する。推定した確率が所定の閾値より低い場合には、エラーによる仮説を棄却することになるため、選択した塩基種であると決定する。これにより、1000分子の中の1分子のみ存在する塩基種の存在の推定が可能となる。したがって、本開示は低頻度な変異の検出に好適に使用できる。
実施形態によっては、エラーモデルは適宜変更することが可能である。
前述のとおり推定された特定塩基位置における特定塩基種がクローナル核酸混合物中に存在する比率は、下記計算式によって算出することができる。また、該比率は、変異頻度とも呼ばれる。
比率(変異頻度)=(変異箇所にマッピングされた変異を有するコンセンサスリードのUMI数)/(変異箇所にマッピングされたコンセンサスリードの全UMI数)
比率(変異頻度)=(変異箇所にマッピングされた変異を有するコンセンサスリードのUMI数)/(変異箇所にマッピングされたコンセンサスリードの全UMI数)
算出された該比率は、クローナル核酸混合物中に存在する低頻度な変異の割合(存在率)と見なすことができる。
各UMI family内において、解析対象とする特定領域の全塩基に対して、特定塩基位置における特定塩基種がクローナル核酸混合物中に存在すると推定できるか、あるいは特定塩基位置における塩基種が不明であると推定できるか、検証し、当該特定塩基種が存在する情報を集約することにより、当該集約された配列情報はクローナル核酸混合物中に存在する低頻度な変異体の配列とすることができる。即ち、特定塩基位置における特定塩基種がクローナル核酸混合物中に存在すると推定できるか、あるいは特定塩基位置における塩基種が不明であると推定できるか、検証することにより得られる結果を、当該塩基のコンセンサスリードに適用して二次配列情報を作成し、更に塩基位置を変更しつつ上述の処理を繰り返すことにより得られる特定領域の全二次配列情報を統合して、クローナル核酸混合物を構成する各拡散分子の特定領域の塩基配列として採用することができる。なお、参照配列と異なるコンセンサスリードが存在しない(参照配列とすべてのコンセンサスリードが同じ)位置については、特定塩基種は存在しないと判断し、クローナル核酸の当該位置の塩基は参照配列と同一であるとする。
必要に応じて、さらに変異フィルタリングを行うことができる。例えば、P. Flaherty等、Nucleic Acids Res.、第40巻、第1号、e2ページ、2012年(DOI: 10.1093/nar/gkr861)に記載の方法を利用し、正常検体で検出された変異を腫瘍検体で検出された変異から除外してもよい。これにより、腫瘍検体に特異的な変異を抽出(フィルタリング)することができる。
3.実施例
以下、本実施の形態の方法を適用した実施例を説明する。本開示の適用は以下の実施例に限定されるものではない。
以下、本実施の形態の方法を適用した実施例を説明する。本開示の適用は以下の実施例に限定されるものではない。
(1)実施例1:断片化された核酸混合物中の低頻度な変異の検出
核酸混合物中に存在する稀な変異の検出感度を調査するために、EGFR、KRAS、NRAS、PIK3CAの4種のヒト遺伝子群において、野生型遺伝子のコピー、さらに同じ遺伝子に稀な変異が生じている8種類の変異型遺伝子のコピーを異なる混合比で有するヒト由来のクローナル核酸として、MultiplexI cfDNA Reference Standard Set(horizon社製)を用いた。本実施例では混合比として、変異型遺伝子のコピー数が1%(Part No.: HD778)、0.1%(Part No.: HD779)の二種類のクローナル核酸を使用した。
核酸混合物中に存在する稀な変異の検出感度を調査するために、EGFR、KRAS、NRAS、PIK3CAの4種のヒト遺伝子群において、野生型遺伝子のコピー、さらに同じ遺伝子に稀な変異が生じている8種類の変異型遺伝子のコピーを異なる混合比で有するヒト由来のクローナル核酸として、MultiplexI cfDNA Reference Standard Set(horizon社製)を用いた。本実施例では混合比として、変異型遺伝子のコピー数が1%(Part No.: HD778)、0.1%(Part No.: HD779)の二種類のクローナル核酸を使用した。
MultiplexI cfDNA Reference Standard Setは、血漿中の無細胞核酸(cell−free DNA)を模して、人工的に160bp程度に断片化したクローナル核酸の標準品である。
それぞれのクローナル核酸から50ngを用いて、後述の操作を実施した。なお、ヒトゲノムDNAは、約30億塩基対から構成され、1塩基の平均分子量は330となるため、10ngのクローナル核酸には約3000コピーのゲノムDNAが存在する。変異型遺伝子の混合比が1%の場合には、30コピーの変異型遺伝子が存在し、混合比が0.1%の場合には3コピーの変異遺伝子が存在することになる。
ThruPLEX Tag−seq 6S (12) Kit(タカラバイオ社製)のプロトコールを使用して、50ngのクローナル核酸の両末端にStem−Loopアダプターをライゲーション法により付加し、テンプレート核酸とした。次いで、得られたテンプレート核酸を鋳型とし、アダプター配列を有するユニバーサルプライマーによりPCR増幅を実施してイルミナ社のシーケンサー用の増幅核酸物(アダプター+テンプレート核酸)を調製した。ThruPLEX技術では、6塩基から成るランダムな配列のUMIがStem−Loopアダプター中に含まれているため、増幅前に個々のテンプレート核酸にUMIを付加することになる(合計2個のUMIが含まれる)。
対象領域となる8個の遺伝子領域に由来するテンプレート核酸のみをシーケンスするために、カスタムのSureSelectのキャプチャプローブ(アジレント社製)を用いて、増幅核酸物から8個の遺伝子領域を含む特定ターゲット領域を濃縮した。具体的な操作は、ThruPLEX Tag−seq Kitにて推奨されているマニュアルに従って行った。
イルミナ社の次世代シーケンサーHiSeqを使用する超並列シーケンシングアプローチを用いて、各検体に対する両末端100bpのシーケンスを実施して塩基配列データ(リード配列データ)を取得した。
取得したリード配列データに対して、Trimmomatic(version:0.32)を使用して低品質なリード配列データを除去した。
リード配列データを、BWA−MEM(version:0.7.15)を用いてヒトの参照配列(hg19)にマッピングを実施した。得られたSAM/BAM形式ファイルより、各リード配列について参照配列に対するマッピング情報とUMIの配列情報を入手した。
上記の手法に従って、マッピングの位置情報とUMIの配列情報に基づいて、特定ターゲット領域に対して分子バーコードのクラスタリングを実施した。分類の条件は下記のとおりである。
条件1: ミスマッチ2塩基以内のものは同一UMIとする。
条件2: 条件1を満たさない場合でも、少なくとも一つのUMI配列が同一であれば、同じUMIとする。
条件1: ミスマッチ2塩基以内のものは同一UMIとする。
条件2: 条件1を満たさない場合でも、少なくとも一つのUMI配列が同一であれば、同じUMIとする。
上記の手法に従って、コンセンサスリードの作成を実施した。その後、第3染色体に座乗するPIK3CA遺伝子のエキソン10に該当する178,935,997−178,936,121bpの各塩基について、1)イルミナ社シーケンサーのリード配列に付与されるPhredクオリティスコアから算出されるエラー率Perrorと、2)UQから算出されたエラー率1−UQ(base)を算出して、二種類のエラー率を比較した。その結果を図10に示す。図10では、分子バーコードのクラスタリング後のSAM/BAMファイルを用いて、各箇所でUMI数が最も多い塩基を対象にして上記二種類のエラー率を算出し、Phredスケール(−10log10(エラー率))でプロットした。
図10に示すように、Phredクオリティスコアから算出されたエラー率(図10中の×)は、Phredスケールで30〜40の間となっているが、UQに基づいたエラー率(図10中の◆)は50を超える値を示しており、全領域でエラー率が約1/10倍に抑えられた。
本開示の方法にて、8種類の変異型遺伝子が検出できるかどうかを検証した。さらに従来法との比較として、1) LoFreq(version: 2.1.2 http://csb5.github.io/lofreq/)と、2) SmCounter (version: https://github.com/xuchang116/smCounter 170823取得)とを用いて変異検出を実施した。その結果を図11に示す。図11において、グレーのハッチング部の変異は、検出された変異を示し、REFは参照配列のリード配列数を示し、ALTは変異配列を有するリード配列数を示し、AFはリード配列比から算出されるアリル頻度を示す。
ここで、LoFreqは、論文(Wilm A等、Nucleic Acids Res. 2012 40(22): 11189−11201)で報告された変異検出ソフトであり、Phred quality scoreから算出されたシーケンスエラー率よりベルヌーイ試行を実施して変異を検出するソフトである。図11では、BWAでマッピングしたSAM/BAMファイルに対して、Connor(version: 0.5.1 https://github.com/umich−brcf−bioinf/Connor)を用いてコンセンサスリードを取得した後に変異検出を実施した。
ここで、SmCounterは、ベイズアプローチにより分子バーコードに基づいてエラー補正を実施して変異を検出するソフトである(非特許文献3)。図11では、分子バーコードのクラスタリング(本開示手法)を実施したSAM/BAM形式ファイルを入力ファイルとして、SmCounterを用いて変異検出を行った結果を示した。
図11に示すとおり、LoFreqでは混合比1.0%の8種の変異中、ハッチングで表示した、EGFR遺伝子のΔ746−750の欠損変異(Deletion)、T790Mの一塩基置換(SNV)、NRAS遺伝子のQ61KおよびA59TのSNVの4種類が検出され、検出感度は50%(4種類/8種類)であった。一方、SmCounterでは混合比1.0%の8種の変異中、ハッチングで表示した、EGFR遺伝子のΔ746−750の欠損変異とV769−D770の挿入変異以外の6種類の変異が検出され、検出感度は75%(6種類/8種類)であった。しかしながら、混合比0.1%の変異にて調査した際には、両ソフトとも8種類のいずれの変異遺伝子も検出できなかった。
図11に記載の通り、本発明の解析コード(図中:「独自パイプライン」)は、ハッチングで表示した、混合比1.0%の8種類の変異型遺伝子を全て検出し、検出感度は100%(8種類/8種類)であった。さらに混合比0.1%の変異中、ハッチングで表示した、EGFR遺伝子のL858RのSNVおよびΔ746−750のDeletion、NRAS遺伝子のG12DおよびA59TのSNV、ならびにPIK3CAのE545KのSNVの5種類の変異を検出し、検出感度は62.5%(5種類/8種類)であった。検出できた該変異のp値はいずれもE−06 (10−6)オーダーであり、閾値(1%水準)より低かった。
以上、実施例1の結果より、本実施形態の検出手法は従来法より優れた検出感度を示し、さらには従来法では検出できなかった極めて低い変異においても、変異を検出することができた。
(2)実施例2:核酸混合物中における変異頻度の定量
核酸混合物中に存在する稀な変異の混合比率(変異頻度)を調査するため、19種の遺伝子群において、野生型遺伝子のコピー、さらに同じ遺伝子に稀な変異が生じている35種類の変異型遺伝子のコピーを1.0%−30%の様々な混合比で有するTru−Q7(1.3% Tier)Reference Standard(Part No.: HD734)(horizon社製)をテンプレート核酸に用いた。
核酸混合物中に存在する稀な変異の混合比率(変異頻度)を調査するため、19種の遺伝子群において、野生型遺伝子のコピー、さらに同じ遺伝子に稀な変異が生じている35種類の変異型遺伝子のコピーを1.0%−30%の様々な混合比で有するTru−Q7(1.3% Tier)Reference Standard(Part No.: HD734)(horizon社製)をテンプレート核酸に用いた。
実施例1と同じく、ThruPLEX Tag−seq 6S (12) Kitを用いてイルミナ社シーケンサー用のライブラリーを作製して、HiSeqによりシーケンスリードを取得した。
実施例1と同様に変異検出を実施した結果、35種類の変異中34種類の変異を検出できた(検出感度 97.1%)。結果を図12に示した。これらの検出された変異の頻度を分子バーコードのクラスタリング後のSAM/BAMファイルより算出した。算出方法は、下記のとおりである。
変異頻度は以下の式によって算出される。
変異頻度=(変異箇所にマッピングされた変異を有するコンセンサスリードのUMI数)/(変異箇所にマッピングされたコンセンサスリードの全UMI数)
変異頻度=(変異箇所にマッピングされた変異を有するコンセンサスリードのUMI数)/(変異箇所にマッピングされたコンセンサスリードの全UMI数)
上記計算式により算出された変異頻度(Observed Allele Frequency)とDroplet Digital PCRTM(Bio−Rad社製)によって求められた変異頻度(Expected Allele Frequency、製品情報に明記された値)を比較して算出した。結果を図12に示した。その結果、本開示の方法により推定された頻度は、Digital PCRで求められた値と極めて近く相関係数R2は0.9945であった。
以上、実施例2の結果より、本開示の方法での分子バーコードのクラスタリングにより精度よくUMI分類が実現できた。また、分解が進んだ(断片化された)低品質のゲノムDNAであっても低頻度な変異を高精度に検出することができた。
4.変異検出装置
図13は、上述した変異検出方法を実施する変異検出装置10の具体的な内部構成を示したブロック図である。変異検出装置10はパーソナルコンピュータのような情報処理装置で構成される。変異検出装置10は、その全体動作を制御する制御部11と、画面表示を行う表示部13と、ユーザが操作を行う操作部15と、データやプログラムを記憶するデータ記憶部17と、外部機器やネットワークと接続し、通信を行なうインタフェース部18と、制御部11が処理を行うためのメインメモリであるRAM(Random Access Memory)16とを備える。
図13は、上述した変異検出方法を実施する変異検出装置10の具体的な内部構成を示したブロック図である。変異検出装置10はパーソナルコンピュータのような情報処理装置で構成される。変異検出装置10は、その全体動作を制御する制御部11と、画面表示を行う表示部13と、ユーザが操作を行う操作部15と、データやプログラムを記憶するデータ記憶部17と、外部機器やネットワークと接続し、通信を行なうインタフェース部18と、制御部11が処理を行うためのメインメモリであるRAM(Random Access Memory)16とを備える。
表示部13は例えば、液晶ディスプレイや有機ELディスプレイで構成される。操作部15はキーボード、マウス、タッチパネル等を含む。
インタフェース部18は、USB、HDMI(登録商標)、SCSI等のインタフェースに準拠した種々の機器(プリンタ、通信装置、入力装置等)を接続可能であり、接続した機器と変異検出装置10間のデータや制御コマンドの通信を可能とする。特に、変異検出装置10はインタフェース部18を介して次世代シーケンサー50からリード配列データを取得する。
制御部11は変異検出装置10全体の動作を制御するものであり、プログラムを実行することで所定の機能を実現するCPUやMPUで構成される。制御部11で実行されるプログラムは通信回線や、CD、DVD、メモリカード等の記録媒体を介して提供されることができる。または、制御部11は、所定の機能を実現するように設計された専用のハードウェア回路(FPGA,ASIC等)で構成されてもよい。
該プログラムは、例えば、JAVA(登録商標)、C、C++、C#、OBJECTIVE−C、SWIFTなどのコンピュータ言語、またはPERL、PYTHON、BIO−RUBYなどのスクリプト言語など任意の適切な言語を用いて、例えば、従来型又はオブジェクト指向技術を用いて、ソフトウェアコードとして実行することができる。
該ソフトウェアコードは、RAM16など記憶及び/又は送信のためのコンピュータ読み取り可能な媒体において、一連の命令として記憶することができる。あるいは、コード化され、インターネットを含む種々のプロトコールに適合する有線、光、及び/又は無線通信回線を介して送信するために適合されたキャリア信号を用いて送信することができる。
データ記憶部17はデータやプログラムを記憶する装置であり、例えばハードディスク(HDD)、SSD(SOLID STATE DRIVE)、半導体メモリ素子、光ディスクで構成することができる。データ記憶部21は、制御部11で実行される制御プログラムや、リード配列データ等を格納する。
図14(a)(b)は、図13の変異検出装置10の制御部11により実行される処理を示すフローチャートである。このうち図14(a)のフローチャートは、PCRエラーやシーケンスのエラーを考慮したデータ解析に係る処理を示している。図14(b)のフローチャートは、リードのアライメントを考慮したデータ解析に係る処理を示している。
まず、図14(a)に示すPCRエラーやシーケンスエラーを考慮したデータ解析に係る処理を説明する。
変異検出装置10はインタフェース部18を介してシーケンサー50からリード配列データを取得する(S51)。変異検出装置10の制御部11は、各リード配列データについて参照配列へのマッピング処理を行う(S52)。制御部11は、上述の方法にしたがい参照配列上の位置及びUMIに基づきリード配列データの分類を行う(S53)。その後、制御部11は、コンセンサスリード処理を行い、各UMI familyについてコンセンサスリードを作製する(S54)。
制御部11はエラーモデルを構築する(S55)。具体的には、制御部11は、信頼性スコアUQを計算し、コンセンサスリードのパターンと信頼性スコアUQからエラーモデル(ポアソン分布)のパラメーターλを求めて、エラーモデルを構築する。
そして、制御部11は、エラーモデルを参照して変異を検出する(S56)。具体的には、制御部11は、変異を確認したい位置と塩基種を選択し、エラーモデルを参照してその変異の発生確率を推定する。そして、制御部11は、推定値が閾値より小さいときは、変異が発生していないと判断し、選択した塩基種が正しいと判断する。図14(a)に示すフローチャートでは、小さいサイズの変異が検出され得る。
任意の工程として、必要であれば、制御部11は変異フィルタリングを行い(S57)、所定の条件を満たす変異を変異検出結果から除外する。
次に、図14(b)に示す、リードのアライメントを考慮したデータ解析に係る処理を説明する。変異検出装置10はインタフェース部18を介してシーケンサー50からリード配列データを取得する(S51)。変異検出装置10の制御部11は、各リード配列データについて参照配列へのマッピング処理を行う(S52)。制御部11は、上述の方法にしたがい参照配列上の位置及びUMIに基づきリード配列データの分類を行う(S53)。
図14(b)に示すフローチャート(処理)では、「分子バーコードのクラスタリング(S53)」と(後続の)「コンセンサスリード処理(S54)」との間に、制御部11は、参照配列にはない検体特有のリード配列長以上の配列領域を決定するローカルアセンブリ処理(S58)を行う。
その後、制御部11は、コンセンサスリード処理を行い、各UMI familyについてコンセンサスリードを作製する(S54)。次に、制御部11はエラーモデルを構築する(S55)。具体的には、制御部11は、信頼性スコアUQを計算し、コンセンサスリードのパターンと信頼性スコアUQからエラーモデル(ポアソン分布)のパラメーターλを求めて、エラーモデルを構築する。更に、制御部11は、エラーモデルを参照して変異を検出する(S56)。具体的には、制御部11は、変異を確認したい位置と塩基種を選択し、エラーモデルを参照してその変異の発生確率を推定する。そして、制御部11は、推定値が閾値より小さいときは、変異が発生していないと判断し、選択した塩基種が正しいと判断する。更に必要であれば、制御部11は変異フィルタリングを行い(S57)、所定の条件を満たす変異を変異検出結果から除外する。
図14(b)に示すフローチャート(処理)では、「分子バーコードのクラスタリング(S53)」と(後続の)「コンセンサスリード処理(S54)」との間に、「ローカルアセンブリ処理(S58)」が存在すること以外は、図14(a)の処理と同様である。但し、図14(b)に示すフローチャートでは、大きいサイズの変異が検出され得る。
なお、図15のフローチャートに示すように、
・リード配列処理(S51)、
・参照配列のマッピング(S52)、及び、
・分子バーコードのクラスタリング(UMIの分類)(S53)
の処理以降に、
・コンセンサスリード処理(S54a)、
・エラーモデル構築(S55a)、及び、
・変異検出(小さい変異)(S56a)
の、小さい変異の検出のための処理フローと、
・ローカルアセンブリ処理(S58)
・コンセンサスリード処理(S54b)、
・エラーモデル構築(S55b)、及び、
・変異検出(大きい変異)(S56b)
の、大きい変異の検出のための処理フローとを並行して行い、その後、「小さい変異」の検出と「大きい変異」の検出との和集合の出力データを作成する(「和集合処理(S56c)」、という処理の流れであってもよい。
・リード配列処理(S51)、
・参照配列のマッピング(S52)、及び、
・分子バーコードのクラスタリング(UMIの分類)(S53)
の処理以降に、
・コンセンサスリード処理(S54a)、
・エラーモデル構築(S55a)、及び、
・変異検出(小さい変異)(S56a)
の、小さい変異の検出のための処理フローと、
・ローカルアセンブリ処理(S58)
・コンセンサスリード処理(S54b)、
・エラーモデル構築(S55b)、及び、
・変異検出(大きい変異)(S56b)
の、大きい変異の検出のための処理フローとを並行して行い、その後、「小さい変異」の検出と「大きい変異」の検出との和集合の出力データを作成する(「和集合処理(S56c)」、という処理の流れであってもよい。
以上のようにして変異検出装置10は、シーケンサー50からのリード配列のデータに基づきDNA配列の変異を検出することができる。
(本開示)
上記の実施の形態は以下の技術思想を開示している。
上記の実施の形態は以下の技術思想を開示している。
(1)コンピュータを用いて、クローナル核酸混合物を構成する核酸分子の塩基配列における特定位置に存在する特定塩基種の存在を推定する方法であって、
a)各核酸分子に由来するリード配列を取得する工程、
b)工程aで取得したリード配列を参照配列上にマッピングしマッピング結果を取得する工程、
c)工程bで得られたマッピング結果を基にリード配列のクラスタリング結果を取得する工程、
d)クラスタリング結果を基にコンセンサスリードとコンセンサスリード中の各塩基に対応する信頼性スコアからなるコンセンサスリード情報を取得する工程、
e)クローナル核酸混合物を構成する各核酸分子から工程a〜dを経て得られるコンセンサスリードの集合に含まれる各コンセンサスリードの配列情報より、参照配列上の特定位置に対応する位置に存在する塩基種とその信頼性スコアを選抜する工程、
f)工程eで使用した特定位置を含むコンセンサスリードの総数に対する、工程eで選抜した特定位置に特定の各塩基種を含有するコンセンサスリード数の比率を取得する工程、
g)工程dで取得した信頼性スコアを用いたパラメトリックなエラーモデルにより、設定された閾値または有意水準での工程eで選抜した特定位置において解析エラーによって特定の塩基が出現する出現率を推定する工程、
h)工程fで得られた比率が工程gにより得られた出現率より有意に高い場合に、工程eで選抜した特定位置に特定の塩基種を有する核酸分子が、工程gで設定された閾値または有意水準でクローナル核酸混合物中に存在すると推定する工程、
を含む方法。
a)各核酸分子に由来するリード配列を取得する工程、
b)工程aで取得したリード配列を参照配列上にマッピングしマッピング結果を取得する工程、
c)工程bで得られたマッピング結果を基にリード配列のクラスタリング結果を取得する工程、
d)クラスタリング結果を基にコンセンサスリードとコンセンサスリード中の各塩基に対応する信頼性スコアからなるコンセンサスリード情報を取得する工程、
e)クローナル核酸混合物を構成する各核酸分子から工程a〜dを経て得られるコンセンサスリードの集合に含まれる各コンセンサスリードの配列情報より、参照配列上の特定位置に対応する位置に存在する塩基種とその信頼性スコアを選抜する工程、
f)工程eで使用した特定位置を含むコンセンサスリードの総数に対する、工程eで選抜した特定位置に特定の各塩基種を含有するコンセンサスリード数の比率を取得する工程、
g)工程dで取得した信頼性スコアを用いたパラメトリックなエラーモデルにより、設定された閾値または有意水準での工程eで選抜した特定位置において解析エラーによって特定の塩基が出現する出現率を推定する工程、
h)工程fで得られた比率が工程gにより得られた出現率より有意に高い場合に、工程eで選抜した特定位置に特定の塩基種を有する核酸分子が、工程gで設定された閾値または有意水準でクローナル核酸混合物中に存在すると推定する工程、
を含む方法。
(2)コンピュータを用いて、クローナル核酸混合物を構成する核酸分子の塩基配列における特定位置に存在する特定塩基種がクローナル核酸混合物中で存在する割合を推定する方法であって、
a)各核酸分子に由来するリード配列を取得する工程、
b)工程aで取得したリード配列を参照配列上にマッピングしマッピング結果を取得する工程、
c)工程bで得られたマッピング結果を基にリード配列のクラスタリング結果を取得する工程、
d)クラスタリング結果を基にコンセンサスリードとコンセンサスリード中の各塩基に対応する信頼性スコアからなるコンセンサスリード情報を取得する工程、
e)クローナル核酸混合物を構成する各核酸分子から工程a〜dを経て得られるコンセンサスリードの集合に含まれる各コンセンサスリードの配列情報より、参照配列上の特定位置に対応する位置に存在する塩基種とその信頼性スコアを選抜する工程、
f)工程eで使用した特定位置を含むコンセンサスリードの総数に対する、工程eで選抜した特定位置に特有の各塩基種を含有するコンセンサスリード数の比率を取得する工程、
g)工程dで取得した信頼性スコアを用いたパラメトリックなエラーモデルにより、設定された閾値または有意水準での工程eで選抜した特定位置において解析エラーによって特有の塩基が出現する出現率を推定する工程、
h)工程fで得られた比率が工程gにより得られた出現率より有意に高い場合に、工程eで選抜した特定位置に特定の塩基種を有する核酸分子が、工程gで設定された閾値または有意水準でクローナル核酸混合物中に存在すると推定する工程、及び
i)工程hにおいて工程eで選抜した塩基種が存在すると判断された場合、工程fで得られた特定位置において塩基種を含有するコンセンサスリードの数から、工程gで設定された閾値または有意水準における工程eで選抜した特定位置における特定塩基種の比率を算出し、クローナル核酸混合物中で存在する割合として採用する工程、
を含む方法。
a)各核酸分子に由来するリード配列を取得する工程、
b)工程aで取得したリード配列を参照配列上にマッピングしマッピング結果を取得する工程、
c)工程bで得られたマッピング結果を基にリード配列のクラスタリング結果を取得する工程、
d)クラスタリング結果を基にコンセンサスリードとコンセンサスリード中の各塩基に対応する信頼性スコアからなるコンセンサスリード情報を取得する工程、
e)クローナル核酸混合物を構成する各核酸分子から工程a〜dを経て得られるコンセンサスリードの集合に含まれる各コンセンサスリードの配列情報より、参照配列上の特定位置に対応する位置に存在する塩基種とその信頼性スコアを選抜する工程、
f)工程eで使用した特定位置を含むコンセンサスリードの総数に対する、工程eで選抜した特定位置に特有の各塩基種を含有するコンセンサスリード数の比率を取得する工程、
g)工程dで取得した信頼性スコアを用いたパラメトリックなエラーモデルにより、設定された閾値または有意水準での工程eで選抜した特定位置において解析エラーによって特有の塩基が出現する出現率を推定する工程、
h)工程fで得られた比率が工程gにより得られた出現率より有意に高い場合に、工程eで選抜した特定位置に特定の塩基種を有する核酸分子が、工程gで設定された閾値または有意水準でクローナル核酸混合物中に存在すると推定する工程、及び
i)工程hにおいて工程eで選抜した塩基種が存在すると判断された場合、工程fで得られた特定位置において塩基種を含有するコンセンサスリードの数から、工程gで設定された閾値または有意水準における工程eで選抜した特定位置における特定塩基種の比率を算出し、クローナル核酸混合物中で存在する割合として採用する工程、
を含む方法。
(3)コンピュータを用いて、クローナル核酸混合物を構成する各核酸分子の特定領域の塩基配列推定方法であって、
a)各核酸分子に由来するリード配列を取得する工程、
b)工程aで取得したリード配列を参照配列上にマッピングしマッピング結果を取得する工程、
c)工程bで得られたマッピング結果を基にリード配列のクラスタリング結果を取得する工程、
d)クラスタリング結果を基にコンセンサスリードとコンセンサスリード中の各塩基に対応する信頼性スコアからなるコンセンサスリード情報を取得する工程、
e)クローナル核酸混合物を構成する各核酸分子から工程a〜dを経て得られるコンセンサスリードの集合に含まれる各コンセンサスリードの配列情報より、参照配列上の特定位置に対応する位置に存在する塩基種とその信頼性スコアを選抜する工程、
f)工程eで使用した特定位置を含むコンセンサスリードの総数に対する、工程eで選抜した特定位置に特定の各塩基種を含有するコンセンサスリード数の比率を取得する工程、
g)工程eで選抜した特定位置の塩基種に対応する工程dで取得した信頼性スコアを用いたパラメトリックなエラーモデルにより、設定された閾値または有意水準での工程eで確認した特定位置における塩基種の解析エラーによる出現率を推定する工程、
h)工程fで得られた比率が工程gにより得られた出現率より有意に高い場合に、工程eで選抜した特定位置に特定の塩基種を有する核酸分子が、工程gで設定された閾値または有意水準でクローナル核酸混合物中に存在する、あるいは該出現率より低い場合は特定位置の塩基種は不明、と判断する工程、
i)工程hで得られた結果を該塩基の各コンセンサスリードに適用し二次配列情報を作成する工程、及び
j)工程eにおける塩基位置を変更し、工程e〜iを繰り返すことにより得られた特定領域の全二次配列情報を統合し、クローナル核酸混合物を構成する各核酸分子の特定領域の塩基配列として採用する工程、
を含む方法。
a)各核酸分子に由来するリード配列を取得する工程、
b)工程aで取得したリード配列を参照配列上にマッピングしマッピング結果を取得する工程、
c)工程bで得られたマッピング結果を基にリード配列のクラスタリング結果を取得する工程、
d)クラスタリング結果を基にコンセンサスリードとコンセンサスリード中の各塩基に対応する信頼性スコアからなるコンセンサスリード情報を取得する工程、
e)クローナル核酸混合物を構成する各核酸分子から工程a〜dを経て得られるコンセンサスリードの集合に含まれる各コンセンサスリードの配列情報より、参照配列上の特定位置に対応する位置に存在する塩基種とその信頼性スコアを選抜する工程、
f)工程eで使用した特定位置を含むコンセンサスリードの総数に対する、工程eで選抜した特定位置に特定の各塩基種を含有するコンセンサスリード数の比率を取得する工程、
g)工程eで選抜した特定位置の塩基種に対応する工程dで取得した信頼性スコアを用いたパラメトリックなエラーモデルにより、設定された閾値または有意水準での工程eで確認した特定位置における塩基種の解析エラーによる出現率を推定する工程、
h)工程fで得られた比率が工程gにより得られた出現率より有意に高い場合に、工程eで選抜した特定位置に特定の塩基種を有する核酸分子が、工程gで設定された閾値または有意水準でクローナル核酸混合物中に存在する、あるいは該出現率より低い場合は特定位置の塩基種は不明、と判断する工程、
i)工程hで得られた結果を該塩基の各コンセンサスリードに適用し二次配列情報を作成する工程、及び
j)工程eにおける塩基位置を変更し、工程e〜iを繰り返すことにより得られた特定領域の全二次配列情報を統合し、クローナル核酸混合物を構成する各核酸分子の特定領域の塩基配列として採用する工程、
を含む方法。
(4)(3)の方法において、工程bで採用した参照配列とコンセンサスリード上の塩基配列が同じである場合は工程e〜工程hを経ずにそのまま工程iで作成する二次情報の塩基種として採用してもよい。
(5)(1)〜(3)のいずれかの方法において、工程eで選抜する特定の塩基位置が、工程bで採用した参照配列とコンセンサスリード上の塩基種が異なる位置であってもよい。
(6)(1)〜(3)のいずれかの方法において、核酸分子が少なくとも一方の末端に分子バーコードを有してもよい。
(7)(1)〜(6)のいずれかの方法において、各核酸分子が異なる分子バーコードを有してもよい。
(8)(7)の方法において、各核酸分子が一つ以上の異なる分子バーコードを有してもよい。
(9)(1)〜(8)のいずれかの方法において、工程cにおけるクラスタリングをリード配列における分子バーコード領域の配列を指標に実施してもよい。
(10)(9)の方法において、分子バーコードがランダムな並びの塩基配列であってもよい。
(11)(10)の方法において、分子バーコード領域のリード配列におけるミスマッチが、該分子バーコード領域のリード配列間の配列類似度が許容範囲内であってもよい。
(12)(1)〜(11)のいずれかの方法において、信頼性スコアが、クラスタリングされたリード配列の情報から、ベイズアプローチによって取得されてもよい。
(13)(12)の方法において、リード配列が核酸分子の増幅産物のリード配列であってもよい。
(14)(13)の方法において、増幅物がポリメラーゼ連鎖反応により得られてもよい。
(15)(13)または(14)の方法において、ベイズアプローチの尤度関数はシーケンサーエラーと増幅エラーを考慮した関数であってもよい。
(16)(15)の方法において、工程bにおけるマッピング工程において参照配列に対しリード配列で異なる塩基が挿入または欠失に相当すると判断された場合、シーケンサーエラーに関し更に補正を加えた関数を採用してもよい。
(17)(1)〜(3)のいずれかの方法において、工程gのパラメトリックなエラーモデルがポアソン分布に基づいてもよい。
(18)(1)または(2)の方法において、特定位置の特定塩基種が低頻度な変異の対象塩基であってもよい。
(19)(1)〜(3)のいずれかの方法において、工程cから工程dへ直接進んでもよい。
(20)(1)〜(3)のいずれかの方法において、工程cの次に、
k)工程cで得られた各クラスタリングに含まれるリード配列を用いて、ローカルアセンブリを実行し、アセンブル配列のクラスタリング結果を取得する工程、
を行ってもよい。
k)工程cで得られた各クラスタリングに含まれるリード配列を用いて、ローカルアセンブリを実行し、アセンブル配列のクラスタリング結果を取得する工程、
を行ってもよい。
(21)(19)及び(20)の方法によって得られた結果を基に、クローナル核酸混合物を構成する拡散分子の塩基配列における特定位置に存在する特定塩基種を決定する工程を含む方法であってもよい。
(22)(1)ないし(21)のいずれかの方法をコンピュータで実行させるためのコンピュータ読み取り可能なプログラム。
(23)クローナル核酸混合物を構成する核酸分子の塩基配列における特定位置の特定塩基種の存在を推定するための装置であって、
各核酸分子由来のリード配列を取得するリード配列取得部と、
取得したリード配列をマッピング装置に供給し、前記マッピング装置による参照配列上へのマッピング結果を取得するマッピング情報取得部と、
マッピングしたリード配列のクラスタリング結果を得るクラスタリング結果取得部と、
クラスタリング結果を基にコンセンサスリードとコンセンサスリード中の各塩基に対応する信頼性スコアからなるコンセンサスリード情報取得部と、
を備えた、推定装置。
各核酸分子由来のリード配列を取得するリード配列取得部と、
取得したリード配列をマッピング装置に供給し、前記マッピング装置による参照配列上へのマッピング結果を取得するマッピング情報取得部と、
マッピングしたリード配列のクラスタリング結果を得るクラスタリング結果取得部と、
クラスタリング結果を基にコンセンサスリードとコンセンサスリード中の各塩基に対応する信頼性スコアからなるコンセンサスリード情報取得部と、
を備えた、推定装置。
(24)クローナル核酸混合物を構成する核酸分子の塩基配列における特定位置の特定塩基種の存在を推定するための装置であって、
各核酸分子由来のリード配列を取得するリード配列取得部と、
取得したリード配列をマッピング装置に供給し、前記マッピング装置による参照配列上へのマッピング結果を取得するマッピング情報取得部と、
マッピングしたリード配列のクラスタリング結果を得る第1のクラスタリング結果取得部と、
各クラスタリングに含まれるリード配列を用いて、ローカルアセンブリを実行し、アセンブル配列のクラスタリング結果を得る第2のクラスタリング結果取得部と、
第2のクラスタリング結果取得部のクラスタリング結果を基にコンセンサスリードとコンセンサスリード中の各塩基に対応する信頼性スコアからなるコンセンサスリード情報取得部と、
を備えた、推定装置。
各核酸分子由来のリード配列を取得するリード配列取得部と、
取得したリード配列をマッピング装置に供給し、前記マッピング装置による参照配列上へのマッピング結果を取得するマッピング情報取得部と、
マッピングしたリード配列のクラスタリング結果を得る第1のクラスタリング結果取得部と、
各クラスタリングに含まれるリード配列を用いて、ローカルアセンブリを実行し、アセンブル配列のクラスタリング結果を得る第2のクラスタリング結果取得部と、
第2のクラスタリング結果取得部のクラスタリング結果を基にコンセンサスリードとコンセンサスリード中の各塩基に対応する信頼性スコアからなるコンセンサスリード情報取得部と、
を備えた、推定装置。
10 変異検出装置
11 制御部
13 表示部
15 操作部
16 RAM
17 データ記憶部
18 インタフェース部
50 シーケンサー
100 変異検出システム
11 制御部
13 表示部
15 操作部
16 RAM
17 データ記憶部
18 インタフェース部
50 シーケンサー
100 変異検出システム
Claims (24)
- コンピュータを用いて、クローナル核酸混合物を構成する核酸分子の塩基配列における特定位置に存在する特定塩基種の存在を推定する方法であって、
a)各核酸分子に由来するリード配列を取得する工程、
b)工程aで取得したリード配列を参照配列上にマッピングしマッピング結果を取得する工程、
c)工程bで得られたマッピング結果を基にリード配列のクラスタリング結果を取得する工程、
d)クラスタリング結果を基にコンセンサスリードとコンセンサスリード中の各塩基に対応する信頼性スコアからなるコンセンサスリード情報を取得する工程、
e)クローナル核酸混合物を構成する各核酸分子から工程a〜dを経て得られるコンセンサスリードの集合に含まれる各コンセンサスリードの配列情報より、参照配列上の特定位置に対応する位置に存在する塩基種とその信頼性スコアを選抜する工程、
f)工程eで使用した特定位置を含むコンセンサスリードの総数に対する、工程eで選抜した特定位置に特定の各塩基種を含有するコンセンサスリード数の比率を取得する工程、
g)工程dで取得した信頼性スコアを用いたパラメトリックなエラーモデルにより、設定された閾値または有意水準での工程eで選抜した特定位置において解析エラーによって特定の塩基が出現する出現率を推定する工程、
h)工程fで得られた比率が工程gにより得られた出現率より有意に高い場合に、工程eで選抜した特定位置に特定の塩基種を有する核酸分子が、工程gで設定された閾値または有意水準でクローナル核酸混合物中に存在すると推定する工程、
を含む方法。 - コンピュータを用いて、クローナル核酸混合物を構成する核酸分子の塩基配列における特定位置に存在する特定塩基種がクローナル核酸混合物中で存在する割合を推定する方法であって、
a)各核酸分子に由来するリード配列を取得する工程、
b)工程aで取得したリード配列を参照配列上にマッピングしマッピング結果を取得する工程、
c)工程bで得られたマッピング結果を基にリード配列のクラスタリング結果を取得する工程、
d)クラスタリング結果を基にコンセンサスリードとコンセンサスリード中の各塩基に対応する信頼性スコアからなるコンセンサスリード情報を取得する工程、
e)クローナル核酸混合物を構成する各核酸分子から工程a〜dを経て得られるコンセンサスリードの集合に含まれる各コンセンサスリードの配列情報より、参照配列上の特定位置に対応する位置に存在する塩基種とその信頼性スコアを選抜する工程、
f)工程eで使用した特定位置を含むコンセンサスリードの総数に対する、工程eで選抜した特定位置に特有の各塩基種を含有するコンセンサスリード数の比率を取得する工程、
g)工程dで取得した信頼性スコアを用いたパラメトリックなエラーモデルにより、設定された閾値または有意水準での工程eで選抜した特定位置において解析エラーによって特有の塩基が出現する出現率を推定する工程、
h)工程fで得られた比率が工程gにより得られた出現率より有意に高い場合に、工程eで選抜した特定位置に特定の塩基種を有する核酸分子が、工程gで設定された閾値または有意水準でクローナル核酸混合物中に存在すると推定する工程、及び
i)工程hにおいて工程eで選抜した塩基種が存在すると判断された場合、工程fで得られた特定位置において塩基種を含有するコンセンサスリードの数から、工程gで設定された閾値または有意水準における工程eで選抜した特定位置における特定塩基種の比率を算出し、クローナル核酸混合物中で存在する割合として採用する工程、
を含む方法。 - コンピュータを用いて、クローナル核酸混合物を構成する各核酸分子の特定領域の塩基配列推定方法であって、
a)各核酸分子に由来するリード配列を取得する工程、
b)工程aで取得したリード配列を参照配列上にマッピングしマッピング結果を取得する工程、
c)工程bで得られたマッピング結果を基にリード配列のクラスタリング結果を取得する工程、
d)クラスタリング結果を基にコンセンサスリードとコンセンサスリード中の各塩基に対応する信頼性スコアからなるコンセンサスリード情報を取得する工程、
e)クローナル核酸混合物を構成する各核酸分子から工程a〜dを経て得られるコンセンサスリードの集合に含まれる各コンセンサスリードの配列情報より、参照配列上の特定位置に対応する位置に存在する塩基種とその信頼性スコアを選抜する工程、
f)工程eで使用した特定位置を含むコンセンサスリードの総数に対する、工程eで選抜した特定位置に特定の各塩基種を含有するコンセンサスリード数の比率を取得する工程、
g)工程eで選抜した特定位置の塩基種に対応する工程dで取得した信頼性スコアを用いたパラメトリックなエラーモデルにより、設定された閾値または有意水準での工程eで確認した特定位置における塩基種の解析エラーによる出現率を推定する工程、
h)工程fで得られた比率が工程gにより得られた出現率より有意に高い場合に、工程eで選抜した特定位置に特定の塩基種を有する核酸分子が、工程gで設定された閾値または有意水準でクローナル核酸混合物中に存在する、あるいは該出現率より低い場合は特定位置の塩基種は不明、と判断する工程、
i)工程hで得られた結果を該塩基の各コンセンサスリードに適用し二次配列情報を作成する工程、及び
j)工程eにおける塩基位置を変更し、工程e〜iを繰り返すことにより得られた特定領域の全二次配列情報を統合し、クローナル核酸混合物を構成する各核酸分子の特定領域の塩基配列として採用する工程、
を含む方法。 - 工程bで採用した参照配列とコンセンサスリード上の塩基配列が同じである場合は工程e〜工程hを経ずにそのまま工程iで作成する二次配列情報の塩基種として採用することを特徴とする請求項3記載の塩基配列決定方法。
- 工程eで選抜する特定の塩基位置が、工程bで採用した参照配列とコンセンサスリード上の塩基種が異なる位置であることを特徴とする請求項1〜3のいずれかに記載の方法。
- 核酸分子が少なくとも一方の末端に分子バーコードを有することを特徴とする請求項1〜3のいずれかに記載の方法。
- 各核酸分子が異なる分子バーコードを有することを特徴とする請求項6に記載の方法。
- 各核酸分子が一つ以上の異なる分子バーコードを有することを特徴とする請求項7記載の方法。
- 工程cにおけるクラスタリングをリード配列における分子バーコード領域の配列を指標に実施することを特徴とする請求項1〜8のいずれかに記載の方法。
- 分子バーコードがランダムな並びの塩基配列であることを特徴とする請求項9記載の方法。
- 分子バーコード領域のリード配列におけるミスマッチが、該分子バーコード領域のリード配列間の配列類似度が許容範囲内であることを特徴とする請求項10記載の方法。
- 信頼性スコアが、クラスタリングされたリード配列の情報から、ベイズアプローチによって取得されることを特徴とする請求項1〜11のいずれかに記載の方法。
- リード配列が核酸分子の増幅産物のリード配列であることを特徴とする請求項12に記載の方法
- 増幅物がポリメラーゼ連鎖反応により得られることを特徴とする請求項13記載の方法。
- ベイズアプローチの尤度関数がシーケンサーエラーと増幅エラーを考慮した関数であることを特徴とする請求項13または14に記載の方法。
- 工程bにおけるマッピング工程において参照配列に対しリード配列で異なる塩基が挿入または欠失に相当すると判断された場合、シーケンサーエラーに関し更に補正を加えた関数を採用することを特徴とする請求項15記載の方法。
- 工程gのパラメトリックなエラーモデルがポアソン分布に基づくことを特徴とする請求項1〜3いずれかに記載の方法。
- 特定位置の特定塩基種が低頻度な変異の対象塩基であることを特徴とする請求項1または2記載の方法。
- 工程cから工程dへ直接進むことを特徴とする請求項1〜3のいずれかに記載の方法。
- 工程cの次に、
k)工程cで得られた各クラスタリングに含まれるリード配列を用いて、ローカルアセンブリを実行し、アセンブル配列のクラスタリング結果を取得する工程、
を行うことを特徴とする請求項1〜3のいずれかに記載の方法。 - 請求項19及び請求項20に記載の方法によって得られた結果を基に、クローナル核酸混合物を構成する拡散分子の塩基配列における特定位置に存在する特定塩基種を決定する工程を含む方法。
- 請求項1〜21のいずれかに記載の方法をコンピュータで実行させるためのコンピュータ読み取り可能なプログラム。
- クローナル核酸混合物を構成する核酸分子の塩基配列における特定位置の特定塩基種の存在を推定するための装置であって、
各核酸分子由来のリード配列を取得するリード配列取得部と、
取得したリード配列をマッピング装置に供給し、前記マッピング装置による参照配列上へのマッピング結果を取得するマッピング情報取得部と、
マッピングしたリード配列のクラスタリング結果を得るクラスタリング結果取得部と、
クラスタリング結果を基にコンセンサスリードとコンセンサスリード中の各塩基に対応する信頼性スコアからなるコンセンサスリード情報取得部と、
を備えた
推定装置。 - クローナル核酸混合物を構成する核酸分子の塩基配列における特定位置の特定塩基種の存在を推定するための装置であって、
各核酸分子由来のリード配列を取得するリード配列取得部と、
取得したリード配列をマッピング装置に供給し、前記マッピング装置による参照配列上へのマッピング結果を取得するマッピング情報取得部と、
マッピングしたリード配列のクラスタリング結果を得る第1のクラスタリング結果取得部と、
各クラスタリングに含まれるリード配列を用いて、ローカルアセンブリを実行し、アセンブル配列のクラスタリング結果を得る第2のクラスタリング結果取得部と、
第2のクラスタリング結果取得部のクラスタリング結果を基にコンセンサスリードとコンセンサスリード中の各塩基に対応する信頼性スコアからなるコンセンサスリード情報取得部と、
を備えた
推定装置。
Applications Claiming Priority (3)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
JP2017253510 | 2017-12-28 | ||
JP2017253510 | 2017-12-28 | ||
PCT/JP2018/048502 WO2019132010A1 (ja) | 2017-12-28 | 2018-12-28 | 塩基配列における塩基種を推定する方法、装置及びプログラム |
Publications (1)
Publication Number | Publication Date |
---|---|
JPWO2019132010A1 true JPWO2019132010A1 (ja) | 2021-01-21 |
Family
ID=67063856
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
JP2019562510A Pending JPWO2019132010A1 (ja) | 2017-12-28 | 2018-12-28 | 塩基配列における塩基種を推定する方法、装置及びプログラム |
Country Status (2)
Country | Link |
---|---|
JP (1) | JPWO2019132010A1 (ja) |
WO (1) | WO2019132010A1 (ja) |
Families Citing this family (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN113005188A (zh) * | 2020-12-29 | 2021-06-22 | 阅尔基因技术(苏州)有限公司 | 用一代测序评估样本dna中碱基损伤、错配和变异的方法 |
CN116741274B (zh) * | 2023-02-07 | 2024-07-26 | 杭州联川基因诊断技术有限公司 | 一种确定靶向测序数据中代表性序列的方法、设备和介质 |
Family Cites Families (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US8370079B2 (en) * | 2008-11-20 | 2013-02-05 | Pacific Biosciences Of California, Inc. | Algorithms for sequence determination |
WO2011037990A1 (en) * | 2009-09-22 | 2011-03-31 | President And Fellows Of Harvard College | Entangled mate sequencing |
EP3087204B1 (en) * | 2013-12-28 | 2018-02-14 | Guardant Health, Inc. | Methods and systems for detecting genetic variants |
GB201409282D0 (en) * | 2014-05-23 | 2014-07-09 | Univ Sydney Tech | Sequencing process |
JP6125731B2 (ja) * | 2014-07-02 | 2017-05-10 | 株式会社Dnaチップ研究所 | 核酸分子数計測法 |
JP2018513508A (ja) * | 2015-03-16 | 2018-05-24 | パーソナル ジノーム ダイアグノスティクス, インコーポレイテッド | 核酸を分析するためのシステムおよび方法 |
GB201607629D0 (en) * | 2016-05-01 | 2016-06-15 | Genome Res Ltd | Mutational signatures in cancer |
-
2018
- 2018-12-28 WO PCT/JP2018/048502 patent/WO2019132010A1/ja active Application Filing
- 2018-12-28 JP JP2019562510A patent/JPWO2019132010A1/ja active Pending
Also Published As
Publication number | Publication date |
---|---|
WO2019132010A1 (ja) | 2019-07-04 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
JP7119014B2 (ja) | まれな変異およびコピー数多型を検出するためのシステムおよび方法 | |
Sedlazeck et al. | Piercing the dark matter: bioinformatics of long-range sequencing and mapping | |
KR102638152B1 (ko) | 서열 변이체 호출을 위한 검증 방법 및 시스템 | |
EP3143537B1 (en) | Rare variant calls in ultra-deep sequencing | |
AU2022203114A1 (en) | Detecting mutations for cancer screening and fetal analysis | |
CN107849612B (zh) | 比对和变体测序分析管线 | |
US10468121B2 (en) | Phasing and linking processes to identify variations in a genome | |
US20210272649A1 (en) | Systems and methods for automating rna expression calls in a cancer prediction pipeline | |
Abnizova et al. | Computational errors and biases in short read next generation sequencing | |
JP6240210B2 (ja) | 標的シーケンシングリードの正確かつ迅速なマッピング | |
US20200327957A1 (en) | Detection of deletions and copy number variations in dna sequences | |
Mittempergher et al. | MammaPrint and BluePrint molecular diagnostics using targeted RNA next-generation sequencing technology | |
CN108292327A (zh) | 下一代测序中检测拷贝数变异的方法 | |
Brozynska et al. | Direct chloroplast sequencing: comparison of sequencing platforms and analysis tools for whole chloroplast barcoding | |
CN109461473B (zh) | 胎儿游离dna浓度获取方法和装置 | |
JPWO2019132010A1 (ja) | 塩基配列における塩基種を推定する方法、装置及びプログラム | |
Roy et al. | NGS-μsat: bioinformatics framework supporting high throughput microsatellite genotyping from next generation sequencing platforms | |
Heo | Improving quality of high-throughput sequencing reads | |
JP2019525308A (ja) | 合成wgsバイオインフォマティクスの検証 | |
Hu et al. | Processing UMI Datasets at High Accuracy and Efficiency with the Sentieon ctDNA Analysis Pipeline | |
Niehus et al. | PopDel identifies medium-size deletions jointly in tens of thousands of genomes | |
US11001880B2 (en) | Development of SNP islands and application of SNP islands in genomic analysis | |
US20220399079A1 (en) | Method and system for combined dna-rna sequencing analysis to enhance variant-calling performance and characterize variant expression status | |
Liu et al. | AsmMix: an efficient haplotype-resolved hybrid de novo genome assembling pipeline | |
WO2023031641A1 (en) | Methods and devices for non-invasive prenatal testing |