CN113593636B - 测序结果分析方法、系统及计算机可读存储介质和电子设备 - Google Patents
测序结果分析方法、系统及计算机可读存储介质和电子设备 Download PDFInfo
- Publication number
- CN113593636B CN113593636B CN202010865293.5A CN202010865293A CN113593636B CN 113593636 B CN113593636 B CN 113593636B CN 202010865293 A CN202010865293 A CN 202010865293A CN 113593636 B CN113593636 B CN 113593636B
- Authority
- CN
- China
- Prior art keywords
- reads
- sequencing
- read
- quality
- sequencing data
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000012163 sequencing technique Methods 0.000 title claims abstract description 468
- 238000004458 analytical method Methods 0.000 title claims abstract description 48
- 238000000034 method Methods 0.000 claims abstract description 86
- 238000012937 correction Methods 0.000 claims abstract description 28
- 230000000903 blocking effect Effects 0.000 claims description 49
- 230000008569 process Effects 0.000 claims description 35
- 238000006243 chemical reaction Methods 0.000 claims description 27
- 230000037430 deletion Effects 0.000 claims description 16
- 238000012217 deletion Methods 0.000 claims description 16
- 238000001914 filtration Methods 0.000 claims description 15
- 230000037431 insertion Effects 0.000 claims description 15
- 238000003780 insertion Methods 0.000 claims description 15
- 108010014303 DNA-directed DNA polymerase Proteins 0.000 claims description 14
- 102000016928 DNA-directed DNA polymerase Human genes 0.000 claims description 14
- 230000000295 complement effect Effects 0.000 claims description 14
- 108010008286 DNA nucleotidylexotransferase Proteins 0.000 claims description 11
- 102100029764 DNA-directed DNA/RNA polymerase mu Human genes 0.000 claims description 11
- 238000010276 construction Methods 0.000 claims description 11
- 125000002887 hydroxy group Chemical group [H]O* 0.000 claims description 4
- 238000004590 computer program Methods 0.000 claims description 3
- 150000007523 nucleic acids Chemical class 0.000 claims description 3
- 108020004707 nucleic acids Proteins 0.000 claims description 2
- 102000039446 nucleic acids Human genes 0.000 claims description 2
- 230000002194 synthesizing effect Effects 0.000 claims description 2
- 238000012549 training Methods 0.000 claims description 2
- 238000012545 processing Methods 0.000 claims 2
- 238000009396 hybridization Methods 0.000 description 37
- 238000004140 cleaning Methods 0.000 description 35
- 239000007788 liquid Substances 0.000 description 31
- 239000003153 chemical reaction reagent Substances 0.000 description 22
- FAPWRFPIFSIZLT-UHFFFAOYSA-M Sodium chloride Chemical compound [Na+].[Cl-] FAPWRFPIFSIZLT-UHFFFAOYSA-M 0.000 description 20
- 239000000243 solution Substances 0.000 description 17
- 108020004414 DNA Proteins 0.000 description 15
- 239000000203 mixture Substances 0.000 description 11
- 230000002829 reductive effect Effects 0.000 description 11
- 239000011780 sodium chloride Substances 0.000 description 10
- 239000002773 nucleotide Substances 0.000 description 9
- 125000003729 nucleotide group Chemical group 0.000 description 9
- 238000007789 sealing Methods 0.000 description 9
- TWRXJAOTZQYOKJ-UHFFFAOYSA-L Magnesium chloride Chemical compound [Mg+2].[Cl-].[Cl-] TWRXJAOTZQYOKJ-UHFFFAOYSA-L 0.000 description 8
- WCUXLLCKKVVCTQ-UHFFFAOYSA-M Potassium chloride Chemical compound [Cl-].[K+] WCUXLLCKKVVCTQ-UHFFFAOYSA-M 0.000 description 8
- 239000000047 product Substances 0.000 description 8
- 238000001816 cooling Methods 0.000 description 7
- 238000001514 detection method Methods 0.000 description 7
- 238000003786 synthesis reaction Methods 0.000 description 7
- 239000011324 bead Substances 0.000 description 6
- 230000015572 biosynthetic process Effects 0.000 description 6
- 238000010586 diagram Methods 0.000 description 6
- 230000035772 mutation Effects 0.000 description 6
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 6
- 238000004925 denaturation Methods 0.000 description 5
- 230000036425 denaturation Effects 0.000 description 5
- 239000012634 fragment Substances 0.000 description 5
- 238000013507 mapping Methods 0.000 description 5
- 238000002156 mixing Methods 0.000 description 5
- 239000000523 sample Substances 0.000 description 5
- 102100031051 Cysteine and glycine-rich protein 1 Human genes 0.000 description 4
- 101710185487 Cysteine and glycine-rich protein 1 Proteins 0.000 description 4
- LFQSCWFLJHTTHZ-UHFFFAOYSA-N Ethanol Chemical compound CCO LFQSCWFLJHTTHZ-UHFFFAOYSA-N 0.000 description 4
- ZHNUHDYFZUAESO-UHFFFAOYSA-N Formamide Chemical compound NC=O ZHNUHDYFZUAESO-UHFFFAOYSA-N 0.000 description 4
- 108091092584 GDNA Proteins 0.000 description 4
- 229920004890 Triton X-100 Polymers 0.000 description 4
- 239000013504 Triton X-100 Substances 0.000 description 4
- BFNBIHQBYMNNAN-UHFFFAOYSA-N ammonium sulfate Chemical compound N.N.OS(O)(=O)=O BFNBIHQBYMNNAN-UHFFFAOYSA-N 0.000 description 4
- 229910052921 ammonium sulfate Inorganic materials 0.000 description 4
- 235000011130 ammonium sulphate Nutrition 0.000 description 4
- 239000000872 buffer Substances 0.000 description 4
- DBLXOVFQHHSKRC-UHFFFAOYSA-N ethanesulfonic acid;2-piperazin-1-ylethanol Chemical compound CCS(O)(=O)=O.OCCN1CCNCC1 DBLXOVFQHHSKRC-UHFFFAOYSA-N 0.000 description 4
- 229910001629 magnesium chloride Inorganic materials 0.000 description 4
- 239000001103 potassium chloride Substances 0.000 description 4
- 235000011164 potassium chloride Nutrition 0.000 description 4
- 239000001509 sodium citrate Substances 0.000 description 4
- NLJMYIDDQXHKNR-UHFFFAOYSA-K sodium citrate Chemical compound O.O.[Na+].[Na+].[Na+].[O-]C(=O)CC(O)(CC([O-])=O)C([O-])=O NLJMYIDDQXHKNR-UHFFFAOYSA-K 0.000 description 4
- 108010064535 CCAAT-Enhancer-Binding Protein-beta Proteins 0.000 description 3
- 102100031621 Cysteine and glycine-rich protein 2 Human genes 0.000 description 3
- 102000053602 DNA Human genes 0.000 description 3
- HEMHJVSKTPXQMS-UHFFFAOYSA-M Sodium hydroxide Chemical compound [OH-].[Na+] HEMHJVSKTPXQMS-UHFFFAOYSA-M 0.000 description 3
- 239000007983 Tris buffer Substances 0.000 description 3
- 239000008367 deionised water Substances 0.000 description 3
- 229910021641 deionized water Inorganic materials 0.000 description 3
- 238000007865 diluting Methods 0.000 description 3
- 239000000463 material Substances 0.000 description 3
- 238000012986 modification Methods 0.000 description 3
- 230000004048 modification Effects 0.000 description 3
- 238000000746 purification Methods 0.000 description 3
- 230000008439 repair process Effects 0.000 description 3
- 230000000717 retained effect Effects 0.000 description 3
- 239000006228 supernatant Substances 0.000 description 3
- OAKPWEUQDVLTCN-NKWVEPMBSA-N 2',3'-Dideoxyadenosine-5-triphosphate Chemical compound C1=NC=2C(N)=NC=NC=2N1[C@H]1CC[C@@H](CO[P@@](O)(=O)O[P@](O)(=O)OP(O)(O)=O)O1 OAKPWEUQDVLTCN-NKWVEPMBSA-N 0.000 description 2
- KWIUHFFTVRNATP-UHFFFAOYSA-N Betaine Natural products C[N+](C)(C)CC([O-])=O KWIUHFFTVRNATP-UHFFFAOYSA-N 0.000 description 2
- 108091035707 Consensus sequence Proteins 0.000 description 2
- 238000001712 DNA sequencing Methods 0.000 description 2
- 229910021380 Manganese Chloride Inorganic materials 0.000 description 2
- GLFNIEUTAYBVOC-UHFFFAOYSA-L Manganese chloride Chemical compound Cl[Mn]Cl GLFNIEUTAYBVOC-UHFFFAOYSA-L 0.000 description 2
- KWIUHFFTVRNATP-UHFFFAOYSA-O N,N,N-trimethylglycinium Chemical compound C[N+](C)(C)CC(O)=O KWIUHFFTVRNATP-UHFFFAOYSA-O 0.000 description 2
- 108091034117 Oligonucleotide Proteins 0.000 description 2
- 238000012408 PCR amplification Methods 0.000 description 2
- DBMJMQXJHONAFJ-UHFFFAOYSA-M Sodium laurylsulphate Chemical compound [Na+].CCCCCCCCCCCCOS([O-])(=O)=O DBMJMQXJHONAFJ-UHFFFAOYSA-M 0.000 description 2
- HDRRAMINWIWTNU-NTSWFWBYSA-N [[(2s,5r)-5-(2-amino-6-oxo-3h-purin-9-yl)oxolan-2-yl]methoxy-hydroxyphosphoryl] phosphono hydrogen phosphate Chemical compound C1=2NC(N)=NC(=O)C=2N=CN1[C@H]1CC[C@@H](COP(O)(=O)OP(O)(=O)OP(O)(O)=O)O1 HDRRAMINWIWTNU-NTSWFWBYSA-N 0.000 description 2
- ARLKCWCREKRROD-POYBYMJQSA-N [[(2s,5r)-5-(4-amino-2-oxopyrimidin-1-yl)oxolan-2-yl]methoxy-hydroxyphosphoryl] phosphono hydrogen phosphate Chemical compound O=C1N=C(N)C=CN1[C@@H]1O[C@H](COP(O)(=O)OP(O)(=O)OP(O)(O)=O)CC1 ARLKCWCREKRROD-POYBYMJQSA-N 0.000 description 2
- 229960003237 betaine Drugs 0.000 description 2
- 238000007664 blowing Methods 0.000 description 2
- 238000004364 calculation method Methods 0.000 description 2
- GVPFVAHMJGGAJG-UHFFFAOYSA-L cobalt dichloride Chemical compound [Cl-].[Cl-].[Co+2] GVPFVAHMJGGAJG-UHFFFAOYSA-L 0.000 description 2
- URGJWIFLBWJRMF-JGVFFNPUSA-N ddTTP Chemical compound O=C1NC(=O)C(C)=CN1[C@@H]1O[C@H](COP(O)(=O)OP(O)(=O)OP(O)(O)=O)CC1 URGJWIFLBWJRMF-JGVFFNPUSA-N 0.000 description 2
- 230000001934 delay Effects 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 238000002474 experimental method Methods 0.000 description 2
- 238000010438 heat treatment Methods 0.000 description 2
- 239000005457 ice water Substances 0.000 description 2
- 230000000670 limiting effect Effects 0.000 description 2
- 239000011565 manganese chloride Substances 0.000 description 2
- 235000002867 manganese chloride Nutrition 0.000 description 2
- 229940099607 manganese chloride Drugs 0.000 description 2
- 239000012452 mother liquor Substances 0.000 description 2
- 238000005086 pumping Methods 0.000 description 2
- 239000002096 quantum dot Substances 0.000 description 2
- 239000004593 Epoxy Substances 0.000 description 1
- 206010064571 Gene mutation Diseases 0.000 description 1
- 108700019961 Neoplasm Genes Proteins 0.000 description 1
- 102000048850 Neoplasm Genes Human genes 0.000 description 1
- 101100526762 Neurospora crassa (strain ATCC 24698 / 74-OR23-1A / CBS 708.71 / DSM 1257 / FGSC 987) rpl-28 gene Proteins 0.000 description 1
- 101100092791 Neurospora crassa (strain ATCC 24698 / 74-OR23-1A / CBS 708.71 / DSM 1257 / FGSC 987) rps-14 gene Proteins 0.000 description 1
- 102100030569 Nuclear receptor corepressor 2 Human genes 0.000 description 1
- 101710153660 Nuclear receptor corepressor 2 Proteins 0.000 description 1
- 108091028043 Nucleic acid sequence Proteins 0.000 description 1
- 229910019142 PO4 Inorganic materials 0.000 description 1
- 125000003277 amino group Chemical group 0.000 description 1
- 238000000137 annealing Methods 0.000 description 1
- 239000002981 blocking agent Substances 0.000 description 1
- 238000004422 calculation algorithm Methods 0.000 description 1
- 239000003795 chemical substances by application Substances 0.000 description 1
- 230000001427 coherent effect Effects 0.000 description 1
- 238000007796 conventional method Methods 0.000 description 1
- SUYVUBYJARFZHO-RRKCRQDMSA-N dATP Chemical compound C1=NC=2C(N)=NC=NC=2N1[C@H]1C[C@H](O)[C@@H](COP(O)(=O)OP(O)(=O)OP(O)(O)=O)O1 SUYVUBYJARFZHO-RRKCRQDMSA-N 0.000 description 1
- SUYVUBYJARFZHO-UHFFFAOYSA-N dATP Natural products C1=NC=2C(N)=NC=NC=2N1C1CC(O)C(COP(O)(=O)OP(O)(=O)OP(O)(O)=O)O1 SUYVUBYJARFZHO-UHFFFAOYSA-N 0.000 description 1
- RGWHQCVHVJXOKC-SHYZEUOFSA-J dCTP(4-) Chemical compound O=C1N=C(N)C=CN1[C@@H]1O[C@H](COP([O-])(=O)OP([O-])(=O)OP([O-])([O-])=O)[C@@H](O)C1 RGWHQCVHVJXOKC-SHYZEUOFSA-J 0.000 description 1
- HAAZLUGHYHWQIW-KVQBGUIXSA-N dGTP Chemical compound C1=NC=2C(=O)NC(N)=NC=2N1[C@H]1C[C@H](O)[C@@H](COP(O)(=O)OP(O)(=O)OP(O)(O)=O)O1 HAAZLUGHYHWQIW-KVQBGUIXSA-N 0.000 description 1
- NHVNXKFIZYSCEB-XLPZGREQSA-N dTTP Chemical compound O=C1NC(=O)C(C)=CN1[C@@H]1O[C@H](COP(O)(=O)OP(O)(=O)OP(O)(O)=O)[C@@H](O)C1 NHVNXKFIZYSCEB-XLPZGREQSA-N 0.000 description 1
- 238000013506 data mapping Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000018109 developmental process Effects 0.000 description 1
- 238000001035 drying Methods 0.000 description 1
- 230000009977 dual effect Effects 0.000 description 1
- 238000001976 enzyme digestion Methods 0.000 description 1
- 238000013467 fragmentation Methods 0.000 description 1
- 238000006062 fragmentation reaction Methods 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 230000002452 interceptive effect Effects 0.000 description 1
- 230000036961 partial effect Effects 0.000 description 1
- NBIIXXVUZAFLBC-UHFFFAOYSA-K phosphate Chemical compound [O-]P([O-])([O-])=O NBIIXXVUZAFLBC-UHFFFAOYSA-K 0.000 description 1
- 239000010452 phosphate Substances 0.000 description 1
- 125000002467 phosphate group Chemical group [H]OP(=O)(O[H])O[*] 0.000 description 1
- 238000011002 quantification Methods 0.000 description 1
- 230000011218 segmentation Effects 0.000 description 1
- 241000894007 species Species 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- LENZDBCJOHFCAS-UHFFFAOYSA-N tris Chemical compound OCC(N)(CO)CO LENZDBCJOHFCAS-UHFFFAOYSA-N 0.000 description 1
- 230000003612 virological effect Effects 0.000 description 1
- 238000005406 washing Methods 0.000 description 1
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
- G16B20/00—ICT specially adapted for functional genomics or proteomics, e.g. genotype-phenotype associations
- G16B20/30—Detection of binding sites or motifs
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
- G06N7/00—Computing arrangements based on specific mathematical models
- G06N7/01—Probabilistic graphical models, e.g. probabilistic networks
-
- 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
- G16B20/00—ICT specially adapted for functional genomics or proteomics, e.g. genotype-phenotype associations
- G16B20/50—Mutagenesis
-
- 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
Landscapes
- Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Life Sciences & Earth Sciences (AREA)
- Theoretical Computer Science (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Health & Medical Sciences (AREA)
- Bioinformatics & Computational Biology (AREA)
- Biotechnology (AREA)
- Spectroscopy & Molecular Physics (AREA)
- Medical Informatics (AREA)
- General Health & Medical Sciences (AREA)
- Evolutionary Biology (AREA)
- Proteomics, Peptides & Aminoacids (AREA)
- Biophysics (AREA)
- Analytical Chemistry (AREA)
- Chemical & Material Sciences (AREA)
- General Physics & Mathematics (AREA)
- Molecular Biology (AREA)
- Genetics & Genomics (AREA)
- Probability & Statistics with Applications (AREA)
- Pure & Applied Mathematics (AREA)
- Computing Systems (AREA)
- Artificial Intelligence (AREA)
- Algebra (AREA)
- Data Mining & Analysis (AREA)
- Evolutionary Computation (AREA)
- General Engineering & Computer Science (AREA)
- Mathematical Optimization (AREA)
- Mathematical Analysis (AREA)
- Mathematical Physics (AREA)
- Software Systems (AREA)
- Computational Mathematics (AREA)
- Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)
Abstract
本发明提供了一种测序结果分析方法。该方法包括：所述测序结果包括第一测序数据和第二测序数据，其中，所述第一测序数据和所述第二测序数据均由多个读段构成，所述第一测序数据中的至少一部分所述读段在所述第二测序数据中存在对应读段，所述测序结果分析方法包括：(a)基于所述第一测序数据和所述第二测序数据各自的至少一部分进行相互校正，以便获得最终序列信息。
Description
技术领域
本发明涉及生物信息学领域，具体地，涉及测序结果分析方法、测序结果分析系统、计算机可读存储介质和电子设备。
背景技术
在上世纪八十年代人们提出了单分子测序概念，2003年斯坦福大学生物工程系的教授Stephen Quake博士成功演示了第一个单分子DNA测序实验，2008年Helicos公司的第一台单分子测序仪(HeliScope)上市，2009年Korlach与Turner在《科学》杂志上发表文章，介绍了PacBio单分子测序技术原理，随后，2010年PacBio公司推出了PacBio RS测序系统，并于2011年正式商用，2014年Oxford Nanopore公司在AGBT(基因组生物学技术进展年会)上展示了其MinION测序系统。然而据报道，无论是Helicos、PacBio还是MinION测序平台，其单重测序(Single-Pass)的测序错误率均很高，最高可达30％。许多研究显示，上述测序平台的错误类型主要是InDel，且是随机发生的，可以通过重复读取的方法降低其测序错误率。
有文献报道PacBio可以采用CCS(环形一致序列)克服其SMRT测序技术的高错误率问题。另外，MinION通过2D和1D2的测序方法可以将测序准确率大幅提升，最高可达97％的准确率。
有文献报道，Helicos通过双重测序法(Two-Pass)进行测序可以将其测序中的缺失类型的错误率降低至1％以内，但是其所使用的文库是5’端为特异性接头，3’端为polyA接头，其操作过程比较繁琐复杂；同时其变性洗脱DNA链的步骤使用的是热水，可能会洗脱不完全，从而干扰后续的测序引物杂交以及测序；并且芯片表面的DNA链未经处理也是随后的测序引入错误的来源之一。
由此可见，现有的单分子测序技术需要进一步改进。
发明内容
发明人在使用GenoCare测序平台对人类基因组进行测序研究时发现其测序输出的数据存在较多的噪声，且测序准确率较低，与参考基因组进行比对，比对率(MappedRate)为53.59％±9.14％，唯一比对率(Unique Mapped Rate)为36.82％±8.71％，错误率(Error Rate)为6.65％±1.04％。
GenoCare测序平台技术原理与Helicos的类似，发明人利用合成序列使用two-pass的测序方法进行测序，按照文献(Harris T D,Buzby P R,Babcock H,et al.Single-Molecule DNA Sequencing of a Viral Genome[J].Science,2008,320(5872):106-109.)中所描述的方法进行分析，发现双重测序碱基(Two-Pass base)的错误率与单重测序碱基(Single-Pass base)测序相比仅下降了约30％。
同时，发明人还发现对于GenoCare单分子测序平台，其测序结果中一些特定的碱基组合间或者一些特定的序列后更容易出现缺失，例如连续的G碱基反应后出现缺失的概率更高等。
本发明旨在至少在一定程度上解决相关技术中的技术问题之一。为此，本发明提供了一种有效的测序结果分析方法。
在本发明的第一方面，本发明提出了一种测序结果分析方法。根据本发明的实施例，所述测序结果包括第一测序数据和第二测序数据，其中，所述第一测序数据和所述第二测序数据均由多个读段构成，所述第一测序数据中的至少一部分所述读段在所述第二测序数据中存在对应读段，所述测序结果分析方法包括：(a)基于所述第一测序数据和所述第二测序数据各自的至少一部分进行相互校正，以便获得最终序列信息。根据本发明实施例的方法通过单分子测序平台的双重测序法(two-pass)获得第一测序数据及第二测序数据，并利用第一测序数据和第二测序数据构建校正模型，以便获得在核酸序列中不同的前后碱基组合下，中间碱基发生缺失、插入、突变的概率，并将第一测序数据和第二测序数据进行相互校正，针对两个测序数据中的差异位点，利用校正模型确定该位点的碱基是插入、缺失或突变，以便判断该位点的正确碱基。需要注意的是，对于两个测序数据中的差异位点，可以为碱基的插入、缺失、突变等，所构建的校正模型也可以预测碱基的插入、缺失或者突变，上述突变可以为任意碱基之间的突变。
在本发明的第二方面，本发明提出了一种获得在本发明第一方面所提到的测序结果的方法。根据本发明的实施例，所述测序结果是通过下列步骤获得的：(1)对芯片表面的测序模板进行第一测序，以便通过形成第一新生测序链获得第一测序数据，所述测序模板通过测序接头连接在所述芯片表面上；(2)对至少一部分所述第一新生测序链的3’-末端进行第一封闭处理；和(3)对所述测序模板进行第二测序，以便通过形成第二新生测序链获得第二测序数据。根据本发明实施例的方法可以实现对测序模板的两轮测序，针对同一模板，获得两组测序数据，在第一测序与第二测序之间的封闭处理，可以防止残余的第一新生测序链在第二测序时继续延伸，可以有效保证第二测序的准确性。
在本发明的第三方面，本发明提出了一种测序结果分析系统。根据本发明的实施例，所述测序结果分析系统包括：测序设备，所述测序设备适于通过双重测序法获得测序结果，所述测序结果包括第一测序数据和第二测序数据，所述第一测序数据和所述第二测序数据均由多个读段构成，所述第一测序数据中的至少一部分所述读段在所述第二测序数据中存在对应读段；分析设备，所述分析设备包括校正模块，适用于基于所述第一测序数据和所述第二测序数据各自的至少一部分进行相互校正，以便获得最终序列信息。根据本发明实施例的测序结果分析系统能够有效地实施在本发明第一方面所提出的测序结果分析方法，通过对两轮测序的结果进行相互校正，进而提高测序结果的准确率。另外，如前所述，针对同一模板，进行两轮测序，获得两组测序数据，在第一测序与第二测序之间的封闭处理，可以防止残余的第一新生测序链在第二测序时继续延伸，能够有效地避免在第二轮测序即第二测序过程中产生干扰信号，有效保证第二测序的准确性。
另外，本发明还提供了一种计算机可读存储介质，其上存储有计算机程序，根据本发明的实施例，该程序被处理器执行时实现前面所述方法的步骤。
本发明还提供了一种电子设备，其包括：前面所述的计算机可读存储介质；以及一个或者多个处理器，用于执行所述计算机可读存储介质中的程序。
本发明的附加方面和优点将在下面的描述中部分给出，部分将从下面的描述中变得明显，或通过本发明的实践了解到。
附图说明
图1是根据本发明一个实施例的测序结果分析方法的流程示意图；
图2是根据本发明又一个实施例的测序结果分析方法的流程示意图；
图3是根据本发明又一个实施例的测序结果分析方法的流程示意图；
图4是根据本发明一个实施例的获得Consensus Reads(一致序列/共同序列)的分析方法流程示意图
图5是根据本发明一个实施例的测序方法的流程示意图；
图6是根据本发明又一个实施例的测序方法的流程示意图；
图7是根据本发明又一个实施例的测序方法的流程示意图；
图8是根据本发明一个实施例的获得Reads1和Reads2的测序方法示意图；
图9是根据本发明一个实施例的测序文库构建示意图；
图10是根据本发明又一个实施例的测序结果分析系统的结构示意图；
图11是根据本发明又一个实施例的测序结果分析系统的结构示意图；
图12是根据本发明又一个实施例的测序结果分析系统的结构示意图。
具体实施方式
下面详细描述本发明的实施例，所述实施例的示例在附图中示出。下面通过参考附图描述的实施例是示例性的，旨在用于解释本发明，而不能理解为对本发明的限制。
在本发明的第一方面，本发明提出了一种测序结果分析方法。根据本发明实施例的方法可以对使用单分子测序平台的双重测序法所获得的两轮测序结果进行分析，该分析方法充分考虑单分子测序平台的特点，针对单分子测序平台测序结果的缺陷对测序数据进行优化，充分利用双重测序法(Two-Pass)所获得的坐标一一对应的两轮测序数据，对容易发生突变、插入和缺失的位点进行预测并且建立校正模型，大大提高了测序分析的准确性，避免测序误差。需要注意的是，获得第一测序数据与第二测序数据的方法不限于双重测序法，只要能够针对同一模板获得两则测序数据，并且该两组测序数据可以一一对应即可。
参考图1～4，所述测序结果包括第一测序数据和第二测序数据，所述第一测序数据和第二测序数据是采用双重测序法获得的，其中，所述第一测序数据和所述第二测序数据均由多个读段构成，所述第一测序数据中的至少一部分所述读段在所述第二测序数据中存在对应读段，所述测序结果分析方法包括：基于所述第一测序数据和所述第二测序数据各自的至少一部分进行相互校正，以便获得最终序列信息。
根据本发明的实施例，所述相互校正包括下列步骤：在所述第一测序数据和所述第二测序数据中选择高质量读段和所述高质量读段的对应读段，所述读段的长度不低于预定长度，所述读段具有不低于预定质量阈值的测序质量；和将所述高质量读段与所述高质量读段的对应读段进行比对，并基于所述比对结果进行序列信息校正。根据本发明的一个实施例，所述预定长度可以根据常规测序中读段长度的阈值来确定，在本发明的实施例中，所述预定长度一般为25bp左右，所述预定长度用于过滤噪声序列，提高测序数据比对的准确性。
根据本发明的实施例，通过对两轮测序的结果进行相互校正能够提高测序结果的准确率。另外，在进行第一轮测序，即第一测序之后，通过对残留在芯片表面的新生测序链的3’-末端进行封闭处理，能够有效地避免在第二轮测序即第二测序过程中产生干扰信号。由此，能够进一步提高测序结果的准确性。
参考图2，所述相互校正包括下列步骤：
S100构建第一读段集合
在该步骤中，根据所述读段的长度，基于所述第一测序数据，构建第一读段集合，所述第一读段集合中的每一个读段长度均不低于第一预定长度。
S200构建第二读段集合和第三读段集合
在该步骤中，根据所述对应读段的长度，基于所述第一读段集合，构建第二读段集合和第三读段集合，所述第二读段集合中每一个读段的所述对应读段的长度均不低于第二预定长度，所述第三读段集合中每一个读段的所述对应读段的长度均处于预定长度范围内。
S300构建第四读段集合和第五读段集合
在该步骤中，根据所述第二读段集合中所述读段及其所述对应读段的测序质量，基于所述第二读段集合及其所述对应读段，构建第四读段集合和第五读段集合。
根据本发明的实施例，所述第四读段集合和所述第五读段集合分别是按照下列原则确定的：
将所述第二读段集合中的所述读段与其所述对应读段进行测序质量比较；
选择测序质量高的一方作为所述第四读段集合的元素，选择测序质量低的一方作为所述第五读段集合的元素；
对于测序质量相同的情形，则选择来自所述第二读段集合的所述读段作为所述第四读段集合的元素，则选择所述对应读段作为所述第五读段集合的元素。
S400构建第六读段集合
在该步骤中，利用测序质量，对所述第四读段集合进行过滤处理，以便构建第六读段集合，所述第六读段集合中的所述读段的测序质量均不低于第一预定质量阈值。
S500构建第七读段集合
在该步骤中，利用所述第六读段集合，从所述第五读段集合中选择与所述第六读段集合中的所述读段对应的所述读段，以便构建第七读段集合。
S600将第六读段集合和第七读段集合进行比较，确定第一差异位点
在该步骤中，将所述第六读段集合与所述第七读段集合进行读段比对，并在所述第六读段集合的所述读段上确定第一差异位点
S700对第一差异位点进行校正
在该步骤中，利用预先确定的测序误差预测模型，对所述第一差异位点进行校正，以便确定第一序列信息，所述测序误差预测模型用于确定在测序过程中，差异位点发生插入或者缺失的概率。
参考图3，在得到第一测序信息后，还可以进一步包括：
S400a构建第八读段集合
在该步骤中，利用测序质量，对所述第三读段集合进行过滤处理，以便构建第八读段集合，其中，所述第八读段集合中的所述读段的测序质量均不低于第二预定质量阈值。
S500a构建第九读段集合
在该步骤中，利用所述第八读段集合，从所述第二测序数据中选择与所述第七读段集合中的所述读段对应的所述读段，以便构建第九读段集合。
S600a将第八读段集合和第九读段集合进行比对，确定第二差异位点
在该步骤中，将所述第八读段集合与所述第九读段集合进行读段比对，并在所述第八读段集合的所述读段上确定第二差异位点。
S700a对第二差异位点进行校正确定第二序列信息
在该步骤中，利用所述测序误差预测模型，对所述第二差异位点进行校正，以便确定第二序列信息。
根据本发明的实施例，所述测序误差预测模型是基于所述第一测序数据和所述第二测序数据与参考基因组的比对结果，对朴素贝叶斯模型进行训练获得的。
根据本发明的实施例，针对所述第一差异位点和所述第二差异位点：
如果来自所述第六读段集合的读段在所述差异位点存在碱基，来自所述第七读段集合的对应读段在所述差异位点不存在碱基，并且在所述差异位点发生缺失的概率为50％以上，则保留所述第六读段集合的读段在所述差异位点的碱基作为最终测序结果。
如果来自所述第六读段集合的读段在所述差异位点不存在碱基，来自所述第七读段的读段集合的定读段在所述差异位点存在碱基，并且在所述差异位点发生插入的概率为50％以上，则保留所述第六读段集合的读段在所述差异位点的碱基作为最终测序结果；和
如果来自所述第六读段集合的读段在所述差异位点存在碱基，来自所述第七读段的读段集合的定读段在所述差异位点也存在碱基，则选择所述第六读段集合的读段在所述差异位点的碱基作为最终测序结果。
根据本发明的实施例，所述第一预定长度和所述第二预定长度分别独立地不低于20bp，优选不低于25bp，所述预定长度范围为10～25bp；所述第一预定质量阈值和所述第二预定质量阈值分别独立地不低于50，优选不低于60。
测序方法
在本发明的第二方面，本发明提出了一种降低单分子测序平台(例如GenoCare单分子测序平台)测序输出序列噪声及错误率的测序，参考图4～9，对根据本发明实施例的测序方法进行描述。
根据本发明的实施例，该方法包括：
S10第一测序获得第一测序数据
在该步骤中，对芯片表面的测序模板进行第一测序，以便通过形成第一新生测序链获得第一测序数据，所述测序模板通过测序接头连接在所述芯片表面上。
需要说明的是，在本文中所使用的术语“芯片”是指在测序平台中所使用的测序芯片，只要是通过边合成边测序的原理进行测序，就可以采用本发明的方法进行处理，其中，优选采用单分子测序平台，例如GenoCare单分子测序平台。当然，本领域技术人员能够理解的是，还可以适用于其他的单分子测序平台，在此不再赘述。
参考图6，在步骤S10之前，还可以通过下列步骤获得可以用于单分子测序平台的芯片：
S10a：使测序文库中的文库分子与芯片表面的测序接头进行杂交；
S10b：利用所述文库分子作为初始模板，通过合成互补链形成所述测序模板；和
S10c除去所述初始模板，并对所述芯片表面的核酸分子的3’-末端进行第二封闭处理。
由此，通过第二封闭处理，可以进一步去除残留的活性3’-末端对后续反应的影响。
参考图7，在进行S10c步骤之前，还可以进一步包括S11b：对步骤S10b中延伸不完全的所述互补链的3’-末端进行第三封闭处理。由此，可以进一步提高测序的准确性，降低不期望的测序噪音。
S20对至少一部分所述第一新生测序链的3’-末端进行第一封闭处理在该步骤中，对至少一部分所述第一新生测序链的3’-末端进行第一封闭处理，通过第一封闭处理，可以有效增加有效数据量，降低无效数据对信息分析的干扰。
根据本发明的一个实施例，步骤S20包括去除所述芯片表面的所述第一新生测序链，并对残余在所述芯片表面的所述第一新生测序链的3’-末端进行第一封闭处理。
根据本发明的一个实施例，步骤S20包括对所述第一新生测序链的3’-末端进行第一封闭处理，并去除封闭后的第一新生测序链。
S30第二测序获得第二测序数据
在该步骤中，对所述测序模板进行第二测序，以便通过形成第二新生测序链获得第二测序数据。
根据本发明的实施例，通过进行两轮测序，在进行第一轮测序，即第一测序之后，通过对残留在芯片表面的新生测序链的3’-末端进行封闭处理，能够有效地避免在第二轮测序即第二测序过程中产生干扰信号。由此，能够提高测序结果的准确性。
根据本发明的实施例，所述第一封闭处理、所述第二封闭处理和所述第三封闭处理可以分别独立地通过使3’-末端羟基与延伸反应阻断剂相连而进行的。由此，可以进一步提高封闭效果，从而进一步提高测序的准确性，降低不期望的测序噪音。
根据本发明的实施例，所述延伸反应阻断剂为ddNTP或其衍生物。由此，可以进一步提高封闭效果，从而进一步提高测序的准确性，降低不期望的测序噪音。
根据本发明的实施例，所述第一封闭处理、所述第二封闭处理和所述第三封闭处理分别独立地采用DNA聚合酶和末端转移酶的至少之一进行。由此，可以进一步提高封闭效果，从而进一步提高测序的准确性，降低不期望的测序噪音。
根据本发明的实施例，所述第一封闭处理和所述第三封闭处理分别独立地通过聚合酶连接所述ddNTP或其衍生物，所述第二封闭处理通过所述末端转移酶连接所述ddNTP或其衍生物。由此，可以进一步提高封闭效果，从而进一步提高测序的准确性，降低不期望的测序噪音。
根据本发明具体的实施例，用于GenoCare单分子Two-Pass测序文库构建的接头以及测序引物，所述接头由寡核苷酸链D7-S1-T和带有5’磷酸基团修饰的D9-S2退火获得，所述测序引物为D7S1T-R2P。其中，所述D7-S1-T的序列为SEQ ID NO：1，所述D9-S2的序列为SEQ ID NO：2，所述D7S1T-R2P的序列为SEQ ID NO：3。
首先，使用上述接头和测序引物在GenoCare单分子测序平台进行Two-Pass测序获得Reads1和Reads2：
步骤一：构建Two-Pass测序文库，使用文库制备试剂盒(Universal DNALibrary Prep Kit for IllμMina V2(ND606-01))将退火后的D7-S1-T/D9-S2接头与准备好的片段化的人类gDNA进行连接，连接后无需进行PCR扩增，直接使用纯化试剂盒(VAHTSDNA Clean Beads(N411-01))进行纯化获得目的文库。
步骤二：将步骤一获得的文库与测序芯片表面接头进行杂交。
步骤三：对步骤二中杂交于芯片表面的初始模板进行互补链合成。
步骤四(可选地)：对步骤三中延伸不完全的新生链的3’OH进行封闭，减少其对测序过程的干扰；
步骤五：变性去除步骤二中杂交于芯片表面的初始模板。
步骤六：对芯片表面残余接头的3’OH进行封闭，减少其对测序过程的干扰。
步骤七：以步骤三中合成的互补链为模板杂交所述测序引物D7S1T-R2P。
步骤八：以步骤三中合成的互补链为模板，以步骤七中杂交的测序引物D7S1T-R2P为引物进行Read1测序。
步骤九：变性去除步骤八中的新生测序链。
步骤十：对经过步骤九处理后可能残余的步骤八中的新生测序链的3’OH进行封闭，防止其在Read2测序时继续延伸。
步骤十一：以步骤三中合成的互补链为模板杂交所述测序引物D7S1T-R2P。
步骤十二：以步骤三中合成的互补链为模板，以步骤十一中杂交的测序引物D7S1T-R2P为引物进行Read2测序。
步骤十三：对步骤八和步骤十二获得的测序数据进行拆分获得坐标一一对应的Reads1和Reads2两部分序列。
进一步地，对上述获得的Reads1和Reads2进行分析获得Consensus Reads的分析方法，包括：
步骤十四：校正模型构建，提取步骤十三获得的Reads1和Reads2序列中同一坐标两次测序读长均≥25bp的Reads，分别输出为T1(Read1)和T2(Read2)两个文件，将T1和T2中的Reads分别与参考基因组进行比对，通过朴素贝叶斯方法计算在不同的前后碱基组合下，中间碱基发生Deletion或Insertion的概率。在预测过程中，对于不同的前后碱基组合下的中间碱基，根据模型中发生Deletion或Insertion的概率决定是否保留中间碱基。若Deletion的概率大于50％，则保留中间碱基，反之舍弃该中间碱基。
步骤十五：对步骤十三获得的Reads1数据按读长进行过滤，将读长≥25bp的Reads1序列集合名为Fa1。使用25bp长度对短读序列进行过滤可以去除一部分噪声序列，提高测序数据mapping的准确性。
步骤十六：根据步骤十三获得的Read2的序列读长将与其对应的步骤十五获得的Fa1拆分为两个集合，其中Read2≥25bp时对应的Fa1中的Reads的集合命名为Fa2，10bp≤Read2＜25bp时对应的Fa1中的Reads的集合命名为Fa3。这里将Fa1拆分为两部分进行分析的目的是在提高Consensus Reads的准确率的同时，减少因长度过滤而导致的数据通量的损失。
步骤十七：将步骤十六获得的Fa2中的Reads和与其坐标相对应的步骤十三获得的Reads2中的Reads的Q值进行比较，取Q值较高的Reads(Q值相等时取Fa2中的Reads)的集合命名为Fa4，取Q值较低的Reads(Q值相等时取Reads2中的Reads)的集合命名为Fa5。该步骤的目的是将Reads1和Reads2中的序列划分为测序质量相对较高和较低的两个集合，保证最终输出的Consensus Reads中的序列是两次测序中测序质量较高相对更准确的Reads。
步骤十八：对步骤十七获得的Fa4和Fa5进行进一步的过滤，将Fa4中Q值≥60的Reads的集合命名为Fa6，将与Fa6中Reads坐标一一对应的Fa5中的Reads的集合命名为Fa7。
步骤十九：将步骤十八获得的Fa6和Fa7中的Reads一一对齐，根据其序列的相似性进行分级，并以Fa7为参考序列对Fa6进行校正，在Fa6序列中标记出与Fa7不同的位置，根据步骤十四构建的校正模型逐个判断不同位置碱基是否为Deletion或者Insertion，从而得到校正后的用于输出的Consensus Reads Part1。
该步骤中的不同的位置表述仅针对某一位置上只在Fa6或者Fa7的一条Reads上测出碱基。这时步骤十四构建的矫正模型将判断该碱基是否应该被保留。对于某一位置上Fa6和Fa7上均测出碱基，但是碱基种类不一致，则以Fa6的碱基为准，本模型不对上述情况进行矫正。
步骤二十：对所述步骤十六获得的Fa3中的Reads进行进一步过滤，将所述Fa3中Q值≥60的Reads的集合命名为Fa8。
步骤二十一：提取所述步骤十三获得的Reads2中与所述Fa8中坐标一一对应的Reads的集合，命名为Fa9。
步骤二十二：将所述Fa8与Fa9中的Reads一一对齐，根据其序列的相似性进行分级，并以Fa9为参考序列对Fa8进行校正，在Fa8中标记出与Fa8不同的位置，根据步骤十四构建的校正模型逐个判断不同位置碱基是否为Deletion或者Insertion，从而得到校正后的用于输出的Consensus Reads Part2。
步骤二十三：根据不同的应用需求将不同的相似性级别的所述Consensus ReadsPart1与所述Consensus Reads Part2进行合并，获得输出的Consensus Reads。
根据本发明的实施例，对于步骤二所述文库与测序芯片表面接头杂交过程包括(试剂常规)：
1)将用于杂交的文库在90-100℃预变性2-5分钟；
2)将从所述步骤1)获得的产物迅速在冰水混合物上冷却2分钟以上，获得变性的杂交文库母液；
3)将从所述步骤2)获得的变性的杂交文库母液使用80％的GenoCare杂交液稀释至合适的浓度，优选地为0.1～2nM，得到稀释的杂交文库；
4)将从所述步骤3)获得的30～50μL稀释的杂交文库通入使用复溶试剂预处理好的测序芯片通道，在40～60℃杂交10～30分钟；
5)向芯片通道通入200～1000μL体积的清洗液1，去除步骤4)杂交后剩余的稀释的杂交文库；
6)向芯片通道通入200～1000μL体积的清洗液2，去除步骤5)中的清洗液1，完成文库与测序芯片表面接头杂交。
根据本发明的实施例，所述GenoCare杂交液为3*SSC溶液。
根据本发明的实施例，所述复溶试剂的组分包括：清洗液1，组分包括：150mM的氯化钠，15mM的柠檬酸钠，150mM的4-羟乙基哌嗪乙磺酸，0.1％的十二烷基硫酸钠。
根据本发明的实施例，清洗液3，组分包括：450mM的氯化钠，45mM的柠檬酸钠。
根据本发明的实施例，所述清洗液2的组分包括：150mM的氯化钠，150mM的4-羟乙基哌嗪乙磺酸。
根据本发明的实施例，对于步骤三所述初始模板的互补链合成的过程包括：
1)向芯片通道中通入200～1000μL体积的延伸试剂，在50～70℃条件下反应5～10分钟；
2)向芯片道中通入200～1000μL体积的清洗液1，去除步骤1)中反应后的延伸试剂；
3)向芯片道中通入200～1000μL体积的清洗液2，去除步骤2)中的清洗液1，完成初始模板的互补链的合成。
根据本发明的实施例，所述延伸试剂的组分包括：10～100U/mL的DNA聚合酶，优选地，可以是Bst DNA聚合酶、Bsu DNA聚合酶、Klenow DNA聚合酶等，0.2～2mM的dNTP，0.5～2M的甜菜碱，20mM的三羟甲基氨基甲烷，10mM的氯化钠，10mM的氯化钾，10mM的硫酸铵，3mM的氯化镁，0.1％的Triton X-100，pH值为8.3。
根据本发明的实施例，对于步骤四所述对延伸不完全的链的3’OH进行封闭的过程包括：
1)向芯片通道中通入200～1000μL体积的封闭试剂1，在30～60℃条件下反应5～30分钟；
2)向芯片通道中通入200～1000μL体积的清洗液1，去除步骤1)中反应后的封闭试剂1，完成对延伸不完全的链的3’OH的封闭。
根据本发明的实施例，所述封闭试剂1的组分包括：10～100U/mL的DNA聚合酶，优选地，可以是Klenow DNA聚合酶、Bsu DNA聚合酶、N9 DNA聚合酶等，10～100μM的ddNTP，5mM的氯化锰，20mM的三羟甲基氨基甲烷，10mM的氯化钠，10mM的氯化钾，10mM的硫酸铵，3mM的氯化镁，0.1％的Triton X-100，pH值为8.3。
根据本发明的实施例，对于步骤五所述去除初始模板的过程包括：
1)向芯片通道中通入200～1000μL体积的变性试剂，优选地，所述变性试剂可以是甲酰胺、0.1M的NaOH等，于50～60℃反应2～5分钟；
2)向芯片通道中通入200～1000μL体积的清洗液1，去除步骤1)中反应后的变性试剂和从芯片上变性分离的初始模板；
3)重复步骤1)和步骤2)一次，完成初始模板的去除。
根据本发明的实施例，对于步骤六所述封闭芯片表面残余接头3’-OH的过程包括：
1)向芯片通道中通入200～1000μL体积的清洗液2；
2)向芯片通道中通入200～1000μL体积的封闭试剂2，在30～60℃条件下反应5～30分钟；
3)向芯片通道中通入200～1000μL体积的清洗液1，去除步骤2)中反应后的封闭试剂2，完成对芯片表面残余接头3’-OH的封闭。
根据本发明的实施例，所述封闭试剂2的组分包括：100U/mL的TerminalTransferase(NEB，M0315L)，1×Terminal Transferase Buffer，0.25mM的氯化钴，10～100μM的ddNTP。
根据本发明的实施例，对于步骤七所述杂交测序引物D7S1T-R2P的过程包括：
1)将所述测序引物D7S1T-R2的母液使用清洗液3稀释至合适的浓度，优选地为0.1～1μM，获得稀释的测序引物杂交液；
2)将从步骤1)获得的200～1000μL稀释的测序引物杂交液通入芯片通道，在50～60℃条件下杂交10～30分钟；
3)向芯片通道中通入200～1000μL体积的清洗液1，去除步骤2)杂交后残余的测序引物；
4)向芯片通道中通入200～1000μL体积的清洗液2，去除步骤3)中的清洗液1，完成测序引物的杂交。
根据本发明的实施例，对于步骤八所述Read1的测序过程参照GenoCare单分子双色测序通用试剂盒(备案号：粤深械备20190887号)中的描述进行。
根据本发明的实施例，对于步骤九所述去除新生测序链的过程参照步骤五进行。
根据本发明的实施例，对于步骤十所述残余新生链的3’OH的封闭过程参照步骤四进行。
根据本发明的实施例，对于步骤十一所述杂交测序引物D7S1T-R2P的过程参照步骤七进行。
根据本发明的实施例，对于步骤十二所述Read2的测序过程参照步骤八进行。
根据本发明的实施例，对于步骤十三所述将测序数据进行拆分获得坐标一一对应的Reads1和Reads2两部分序列的过程包括：
使用python语言按测序循环数将BaseCalling输出的“.fa_”文件中的每条Read从中间平均分为两个部分，分别输出为两份序列坐标一致的“.fa_”文件“Reads1.fa_”、“Reads2.fa_”；
使用python语言将从步骤1)获得的“Reads1.fa_”、“Reads2.fa_”文件中所用Reads中的字符“_”移除，输出“Reads1.fa”、“Reads2.fa”文件，完成将测序数据进行拆分获得坐标一一对应的Reads1和Reads2两部分序列。
根据本发明的实施例，对于步骤十四所述校正模型的构建过程包括：
1)提取步骤十三获得的Reads1和Reads2序列中同一坐标两次测序读长均≥25bp的Reads，分别输出为T1(Read1)和T2(Read2)两个fast文件；
2)将从步骤1)获得的T1和T2文件中对应的两条reads进行滑动对齐，标记对齐结果中两条reads相同和不同的Base，得到Common Reads；
3)将从步骤1)获得的T1和T2文件分别和参考序列做mapping得到Sam1和Sam2文件；
4)根据从步骤3)获得的Sam1和Sam2中对应的且mapping到同一位置的reads，找到参考序列中的最长公共子串Ref Reads；
5)将从步骤2)获得的Common Reads中两次测序不同的base与从步骤4)获得的RefReads进行比较，使用朴素贝叶斯方法计算在不同的前后碱基组合下中间碱基发生Deletion或Insertion的概率，完成校正模型的构建。
根据本发明实施例的方法，结合GenoCare单分子测序平台的特点，提供了一套联合使用接头D7-S1-T/D9-S2和测序引物D7S1T-R2P，利用GenoCare单分子测序平台进行Two-Pass测序获得Reads1和Reads2的测序方法。另一方面，本发明提供了一种使用所述Two-Pass测序方法获得的Reads1和Reads2进行分析获得Consensus Reads的分析方法。该分析方法可以显著降低输出的Consensus Reads(一致序列/共同序列)中的噪声序列以及碱基错误率。
测序结果分析系统
在本发明的第三方面，本发明还提出了一种能够实施上述测序结果分析方法的测序结果分析系统。参考图10～12，根据本发明的实施例，该系统包括：测序设备，所述测序设备适于通过双重测序法获得测序结果，所述测序结果包括第一测序数据和第二测序数据，所述第一测序数据和所述第二测序数据均由多个读段构成，所述第一测序数据中的至少一部分所述读段在所述第二测序数据中存在对应读段；分析设备，所述分析设备包括校正模块，适用于基于所述第一测序数据和所述第二测序数据各自的至少一部分进行相互校正，以便获得最终序列信息。
利用该系统能够有效地实施前面所述的测序结果分析方法及测序方法，从而通过对两轮测序的结果进行相互校正能够提高测序结果的准确率。另外，如前所述，在进行第一轮测序，即第一测序之后，通过对残留在芯片表面的新生测序链的3’-末端进行封闭处理，能够有效地避免在第二轮测序即第二测序过程中产生干扰信号。由此，能够进一步提高测序结果的准确性。
另外，本发明还提供了一种计算机可读存储介质，其上存储有计算机程序，其特征在于，该程序被处理器执行时实现前面所述方法的步骤。
本发明还提供了一种电子设备，其包括：前面所述的计算机可读存储介质；以及一个或者多个处理器，用于执行所述计算机可读存储介质中的程序。
下面将结合具体实施例对本发明进行进一步解释说明。下述实施例中所使用的实验方法如无特殊说明，均为常规方法。下述实施例中所用的材料、试剂等，如无特殊说明，均可从商业途径得到。
实施例
本实施例提供了一套降低GenoCare单分子测序平台测序输出序列噪声及错误率的测序与分析方法。其中Genocare单分子测序平台是使用TIRF成像系统检测掺入核苷酸种类的平台。Genocare测序过程有多种方式，第一种方式：四种核苷酸带有同种荧光信号，每轮反应加入一种核苷酸进行信号检测；第二种方式：四种核苷酸带有两种不同的荧光信号，每轮反应加入两种核苷酸进行信号检测；第三种方式：四种核苷酸带有四种不同的荧光信号，每轮反应加入四种核苷酸进行信号检测。具体测序过程可参看文章Single molecμLetargeted sequencing for cancer gene mutation detection，Scientific RepoRts|6:26110|DOI:10.1038/srep26110、专利CN201680047468.3、CN201910907555.7、CN201880077576.4或CN201911331502.1中测序过程的描述。
进一步地，本实施例提供的所述测序与分析方法包括：
一套联合联合使用接头D7-S1-T/D9-S2和测序引物D7S1T-R2P，利用GenoCare单分子测序平台进行Two-pass测序获得Reads1和Reads2的测序方法。其中，所述接头D7-S1-T/D9-S2由寡核苷酸链D7-S1-T和带有5’磷酸基团修饰的D9-S2组成。所述D7-S1-T的序列为SEQ ID NO：1，所述D9-S2的序列为SEQ ID NO：2，所述测序引物D7S1T-R2P的序列为SEQ IDNO：3。具体地，本发明中所涉及的引物序列和名称如表1所示。
表1：引物序列和名称
2)一套对上述Two-Pass测序方法获得的Reads1和Reads2进行分析获得ConsensusReads的分析方法。该分析方法可以显著降低输出的Consensus Reads中的噪声序列以及碱基错误率。
进一步地，本实施例提供的所述一套联合联合使用接头D7-S1-T/D9-S2和测序引物D7S1T-R2P，利用GenoCare单分子测序平台进行Two-Pass测序获得Reads1和Reads2的测序方法包括：
步骤一：构建Two-Pass测序文库。使用Universal DNA Library Prep Kitfor IllμMina V2(ND606-01)将退火后的D7-S1-T/D9-S2接头与准备好的片段化的人类gDNA进行连接，连接后无需进行PCR扩增，直接使用VAHTS DNA Clean Beads(N411-01)进行纯化获得目的文库。
具体地，本实施例中构建Two-Pass测序文库的步骤包括：
1)人类gDNA片段化：使用Covaris设置参数Peak Power，75；Duty Factor，25；Cycles/Burst，50；Time(s)，250对0.1～1ug人类gDNA进行超声打断获得100～300bp的DNA片段。可选地，本步骤也可以使用酶切的方法实现。
2)DNA片段进行末端修复和加A尾，反应体系如表2所示。
表2：反应体系
反应条件为：20℃反应15分钟，接着在65℃条件下反应10分钟。
3)末端修复加A产物与接头进行连接，反应体系如表3所示。
表3：反应体系
| 末端修复加A产物 | 20μL |
| D7-S1-T/D9-S2接头(20μM) | 5μL |
| Ligation Mix | 25μL |
| Total | 50μL |
反应条件为，混匀后室温放置15min。
4)连接产物纯化
纯化使用VAHTS DNA Clean Beads(N411-01)的试剂和说明书所示步骤进行纯化，回收产物10μL，完成测序文库的构建。具体步骤如下：
a)将连接后的PCR体系转移至1.5mL EP管中，加入0.8×(40μL)磁珠，吹打混匀10次，室温放置3分钟；
b)将1.5mL EP管放置在磁力架上，静置2-3分钟，移去上清；
c)用200μL体积80％乙醇洗涤，漂洗磁珠，室温孵育30sec，小心移除上清；
d)开盖干燥磁珠约5-10分钟至残余乙醇完全挥发；
e)加入22μL体积的去离子水从磁力架上去取进行洗脱，充分混匀后室温静置3分钟，置于磁力架上3分钟，待液体澄清后，回收产物20μL，再加入1.2x(24μL)磁珠，吹打混匀10次，室温放置3分钟；
f)将1.5mL EP管放置在磁力架上，静置2-3分钟，移去上清；
g)重复步骤c)—d)一次；
h)加入11μL体积的去离子水从磁力架上取下进行洗脱，充分混匀后室温静置3分钟，置于磁力架上3分钟，待液体澄清后，回收产物10μL，完成测序文库构建。
5)定量及检测
使用Qubit 3.0仪器和Qubit dsDNA HS检测试剂盒对所述构建的文库进行浓度检测。
使用Labchip DNA HS检测试剂盒和LabChip仪器对所述构建的文库进行片段分布检测。
步骤二：将步骤一获得的文库与测序芯片表面探针进行杂交。
芯片选择：
1)芯片选择：所用的芯片为环氧基修饰的芯片，通过探针上的氨基和芯片，序列为5’-TTTTTTTTTTTCCTTGATACCTGCGACCATCCAGTTCCACTCAGATGTGTATAAGAGACAGT-3’(SEQ IDNO:4)。
文库与芯片上探针杂交过程如下：
1)取3μL体积20nM浓度的步骤一构建的测序文库，加入3μL的去离子水，混合均匀，于95℃热变性5分钟；
2)将从步骤1)获得的变性文库迅速置于冰水混合物冷却2分钟以上；
3)向步骤2)的产物中加入24μL体积的GenoCare杂交液，将文库稀释至2nM的工作浓度。其中杂交液3xSSC缓冲液，3xSSC溶液是将20*SSC缓冲液((西格玛，#S6639-1L))用无核酸酶水(Rnase-free水)稀释而成。
4)将从步骤3)获得的30μL体积稀释的杂交文库通入从芯片的一条通道中，于42℃杂交反应30分钟，然后冷却至室温；
5)向步骤4)获得的杂交通道中通入200μL体积的清洗液1，去除未杂交至芯片表面的文库；
向芯片杂交通道通入200μL体积的清洗液2，替换通道内的清洗液1，完成文库与测序芯片表面接头的杂交。
清洗液1组分包括：150mM的氯化钠，15mM的柠檬酸钠，150mM的4-羟乙基哌嗪乙磺酸，0.1％的十二烷基硫酸钠。
清洗液2的组分包括：150mM的氯化钠，150mM的4-羟乙基哌嗪乙磺酸。
步骤三：初始模板进行互补链合成。
初始模板为步骤二中与探针进行杂家的文库，初始模板互补链合成的具体步骤如下：
1)将步骤二中完成文库杂交的芯片置于GenoCare测序仪；
2)向芯片杂交通道泵入750μL体积的延伸试剂，其中，延伸试剂组分为：120U/mLBst DNA聚合酶(NEB，#M0275M)，0.2mM dNTP(dATP、dTTP、dCTP、dGTP各0.2μM的混合物)，1M甜菜碱，20mM的三羟甲基氨基甲烷，10mM的氯化钠，10mM的氯化钾，10mM的硫酸铵，3mM的氯化镁，0.1％的Triton X-100，pH值为8.3；
3)将芯片升温至60±0.5℃，反应10分钟；
4)向芯片杂交通道泵入220μL体积的清洗液1，去除延伸试剂；
5)向芯片杂交通道泵入440μL体积的清洗液2，去除步骤4)中的清洗液1，完成初始模板互补链的合成。
步骤四(可选地)：对步骤三中延伸不完全的新生链的3’OH进行封闭,封闭的具体步骤如下：
1)将芯片降温至37±0.5℃，维持90秒；
2)向步骤三所述延伸后的通道中泵入750μL体积的封闭试剂1，反应10分钟。所述封闭试剂1的组分为：100U/mLKlenow DNA聚合酶大片段(3′→5′exo-,NEB，#M0212M)12.5μMddNTP mix(ddATP、ddTTP、ddCTP、ddGTP各12.5μM的混合物)，5mM的氯化锰，20mM的三羟甲基氨基甲烷，10mM的氯化钠，10mM的氯化钾，10mM的硫酸铵，3mM的氯化镁，0.1％的Triton X-100，pH值为8.3；
3)向步骤2)所述封闭后的通道中通入220μL体积的清洗液1，去除封闭反应后剩余的封闭液，完成对延伸不完全的新生链的3’OH的封闭。
步骤五：变性去除初始模板，去除初始模板的过程如下：
1)将芯片降温至55±0.5℃
2)向步骤四所述封闭后的通道中通入800μL体积的甲酰胺，变性2分钟；
3)向步骤2)所述变性后通道中通入220μL体积的清洗液1，去除变性后的初始模板；
4)重复步骤2)和步骤3)一次，完成对初始模板的去除。
步骤六：对芯片表面残余接头的3’OH进行封闭，封闭芯片表面残余接头3’OH的过程包括：
1)将芯片降温至37±0.5℃；
2)向步骤五所述封闭后的通道中通入440μL体积的清洗液2，替换通道内剩余的清洗液1；
3)向所述步骤2)处理后的通道中通入750μL体积的封闭试剂2，反应15分钟。其中，封闭试剂2的组分为：100U/mL末端转移酶(Terminal Transferase(NEB，M0315L))，1×Terminal Transferase Buffer，0.25mM氯化钴，100μM ddNTP mix(ddATP、ddTTP、ddCTP、ddGTP各100μM的混合物)；
4)向步骤3)所述封闭后的通道中通入220μL体积的清洗液1，完成对芯片表面残余接头3’OH的封闭。
步骤七：杂交测序引物D7S1T-R2P，杂交测序引物D7S1T-R2P的过程如下：
1)将芯片升温至55±0.5℃，保持1分钟；
2)向步骤六所述封闭后的通道中通入800μL体积的稀释的测序引物杂交液，杂交反应30分钟。所述稀释的测序引物杂交液为含有0.1μM引物D7S1T-R2P的清洗液3，清洗液3组分包括：450mM的氯化钠，45mM的柠檬酸钠；
3)将芯片降温至37±0.5℃，保持90秒；
4)向步骤2)所述杂交通道中通入220μL体积的清洗液1，去除通道中未被杂交的测序引物；
5)向所述步骤4)处理后的通道中通入440μL体积的清洗液2，替换通道中剩余的清洗液1，完成测序引物的杂交。
步骤八：进行Read1测序，Read1测序的过程如下：
利用Genocare单分子测序平台进行80个循环的测序，测序过程中采用四种核苷酸带有两种不同的荧光信号，每轮反应加入两种标记不同荧光信号的核苷酸进行信号检测的方式进行测序。
步骤九：去除新生测序链。
去除新生测序链的过程按照步骤五中的步骤进行。
步骤十：封闭残余新生链的3’OH。
封闭残余新生链的3’OH的过程按照步骤四中的步骤进行。
步骤十一：杂交所述测序引物D7S1T-R2P。
杂交测序引物D7S1T-R2P的过程按照步骤七中的步骤进行。
步骤十二：进行Read2测序。
Read2测序的过程按照步骤八中的步骤进行。
步骤十三：将测序数据进行拆分获得坐标一一对应的Reads1和Reads2两部分序列。
具体地，本实施例中所述将测序数据进行拆分获得坐标一一对应的Reads1和Reads2两部分序列的过程包括：
使用python语言将160循环测序BaseCalling输出的“.fa_”文件中的每条Read拆分为前80循环和后80循环两个部分，并将所有Reads中的字符“_”移除，分别输出为两份序列坐标一致的“.fa”文件“Reads1.fa”、“Reads2.fa”，完成将测序数据进行拆分获得坐标一一对应的Reads1和Reads2两部分序列。
进一步地，本实施例中提供的所述一套对上述Two-pass测序方法获得的Reads1和Reads2进行分析获得Consensus Reads的分析方法包括：
步骤十四：构建校正模型。
具体地，本实施例中所述构建校正模型的过程包括：
1)使用python语言，提取步骤十三获得的Reads1和Reads2序列中同一坐标两次测序读长均≥25bp的Reads，分别输出为T1(Read1)和T2(Read2)两个文件。其中同一坐标的对应方法是在生成Reads文件时将同一坐标Reads在不同文件中的Reads ID设置为一致；
2)将T1和T2中位置对应的Reads相互间做Align，在Align结果中标记两条Reads一致和不一致的Base，得到Common Reads。其中位置对应是通过比较两条Reads将的Reads ID是否一致实现；
3)分别将文件T1和T2和Reference做Mapping，得到Sam1和Sam2文件。将Sam1和Sam2中位置对应且mapping到同一位置的Reads，找到Reference中最长公共子串RefReads。公共子串指两条对应的Reads mapping后均覆盖的区域；
4)比较步骤2)中的Common Reads和步骤3)中的Ref Reads。对于Common Reads中不一致的Base，标记其是否真实存在于Reference中。若存在，对于没有测到的Reads则为Deletion。若不存在，对于测到的Reads则为Insertion；
5)统计步骤4)中的Deletion和Insertion情况，同时统计该不一致位置上前后Base的种类。因此得到在不同Base类型前或后引起Insertion或Deletion的概率。
具体地，本实例中运用的朴素贝叶斯模型如下：
其中：P(D|XY)表示对于某碱基在前后分别为X和Y碱基时发生Deletion的概率，X,Y∈[A,C,G,T]。P(D)表示对于某碱基发生Deletion的概率；P(I)表示对于某碱基发生Insertion的概率。
通过统计不同碱基下发生Deletion或Insertion时，前后碱基出现频率即可得到P(XY|D)和P(XY|I)，从而可以计算得到P(D|XY)和P(I|XY)。
步骤十五：过滤读长得到Fa1。
具体地，本实施例中所述读长过滤的过程包括：
使用Python语言逐行读取Reads1文件中所有reads，若Reads长度大于等于25bp，则输出的文本文件Fa1中。
步骤十六：根据Reads2读长，将Fa1中Reads进行分类。
具体地，本实例中所述根据Reads2读长分类Fa1中Reads的过程包括：
将Fa1中所有Reads对应在Reads2中的Reads读出，根据Reads2中Reads的长度，若Read2≥25bp，则将对应的Fa1中Reads保存于Fa2文件中；若10bp≤Read2＜25bp，则将对应的Fa1中Reads保存于Fa3文件中。
步骤十七：根据Q值输出置信Reads。
具体地，本实例中所述根据Q值重新输出置信Reads的过程包括：
1)将步骤十六得到的Fa2中所有Reads取出，并同时取出其对应的Reads2中的Reads。从Reads ID中分割得到该Reads的Quality Score值(简称Q值)。
2)比较两条对应的Reads的Q值，将Q值较大的Reads输出到文件Fa4中，Q值较小的Reads输出到文件Fa5中。若两者Q值相等，则默认将Reads1中的Reads输出至Fa4中，Reads2中的Reads输出至Fa5中。
步骤十八：根据Q值过滤Fa4和Fa5中Reads。
具体地，本实例中所述根据Q值过滤Fa4和Fa5中Reads的过程包括：
取Fa4中Reads，根据其Q值，若大于等于60，则将输出到文件Fa6中，同时将该Reads对应的Fa5中的Reads输出到文件Fa7中。
步骤十九：使用Fa7中Reads矫正Fa6中Reads，得到Consensus Reads Parts1(简称CRP1)。
具体地，本实例中所述使用Fa7中的Reads矫正Fa6中Reads的过程包括：
1)取Fa6中Reads和其对应的在Fa7中的Reads。将两条对应Reads相互配准，得到共同的一致性序列部分。其中两条序列配准使用Smith-Waterman算法，一致性序列指配准后通过在序列中增加、删除或修改部分Base，得到的局部最佳匹配序列。
2)得到一致性序列后，根据步骤十四构建的矫正模型，逐个判断一致性序列中不一致的Base位置。根据该Base位置前后的碱基类型计算该位置出现Deletion或Insertion的概率。若Deletion的概率大于50％，则认为该位置所测Base不应该出现，从而删除该位置Base。反之，保留该位置上的Base。
3)矫正所有不一致Base后，输出矫正后的Reads，即为CRP1。这里的不一致Base特指两条对应Reads中没有同时被测出的Base。若两次均测出该Base，但Base类型不一致，不在本实例矫正的候选范围内，该情况下，最终Base类型以Fa6中Reads的Base类型为准。
步骤二十：根据Q值过滤Fa3中Reads。
具体地，本实例中所述根据Q值过滤Fa3中Reads的过程包括：
取Fa3中所有Reads，分割Fa3中Reads的Reads ID，得到每条Reads的Q值。将Q值≥60的Reads输出到文件Fa8中。
步骤二十一：输出Fa8文件中对应的Reads2中的Reads。
具体地，本实例中所述根据Fa8中Reads输出Reads2中Reads的过程包括：
取Fa8文件中所有Reads，取出其对应的Reads2中的Reads，将其输出到文件Fa9中。
步骤二十二：使用Fa9中Reads矫正Fa8中Reads，得到Consensus Reads Parts2(简称CRP2)。
具体地，本实例中所述的使用Fa9中Reads矫正Fa8中Reads的过程参照所述步骤十九进行。
步骤二十三：根据不同应用对应测序数据准确率的需求，将符合相似度阈值的CRP1和CRP2中的Reads合并输出，得到Consensus Reads。
具体地，本实例中所述不同应用对应测序数据准确率的需求，过滤ConsensusReads Part中Reads并输出的过程包括：
1)根据不同应用对应测序数据准确率的需求，设定对应的相似度阈值。其中对Part1和Part2的相似度阈值可以不同；
2)分别计算CRP1和CRP2中的Reads相似度，相似度是指某Reads在Reads1和Reads2中对应的Reads的相似度。相似度计算步骤是先将两条对应Reads相互配准。再计算配准得到的一致性序列中一致的Base数占总Base数的比值。其中配准方法、一致性序列和不一致Base定义参照步骤十九。
3)根据不同应用对应测序数据准确率的需求，分别将CRP1和CRP2中符合相似度阈值要求的Reads输出到最终的文件中，得到Consensus Reads，参考表4。
表4:不同相似度阈值过滤输出序列与参考基因组mapping分析比较
注：数据损失主要发生在读长过滤步骤，由于Read1和Read2测序是相互独立事件，所以必然存在部分读长不一致的序列。
此外，术语“第一”、“第二”仅用于描述目的，而不能理解为指示或暗示相对重要性或者隐含指明所指示的技术特征的数量。由此，限定有“第一”、“第二”的特征可以明示或者隐含地包括至少一个该特征。在本发明的描述中，“多个”的含义是至少两个，例如两个，三个等，除非另有明确具体的限定。
在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不必须针对的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任一个或多个实施例或示例中以合适的方式结合。此外，在不相互矛盾的情况下，本领域的技术人员可以将本说明书中描述的不同实施例或示例以及不同实施例或示例的特征进行结合和组合。
尽管上面已经示出和描述了本发明的实施例，可以理解的是，上述实施例是示例性的，不能理解为对本发明的限制，本领域的普通技术人员在本发明的范围内可以对上述实施例进行变化、修改、替换和变型。
序列表
<110> 深圳市真迈生物科技有限公司
<120> 测序结果分析方法、系统及计算机可读存储介质和电子设备
<130> PIDC4200100
<160> 4
<170> PatentIn version 3.5
<210> 1
<211> 51
<212> DNA
<213> Artificial Sequence
<220>
<223> 引物
<400> 1
ctcagatcct acaacgacgc tctaccgatg aagatgtgta taagagacag t 51
<210> 2
<211> 51
<212> DNA
<213> Artificial Sequence
<220>
<223> 引物
<400> 2
ctgtctctta tacacatctg agtggaactg gatggtcgca ggtatcaagg a 51
<210> 3
<211> 43
<212> DNA
<213> Artificial Sequence
<220>
<223> 引物
<400> 3
ctacaacgac gctctaccga tgaagatgtg tataagagac agt 43
<210> 4
<211> 62
<212> DNA
<213> Artificial Sequence
<220>
<223> 探针序列
<400> 4
tttttttttt tccttgatac ctgcgaccat ccagttccac tcagatgtgt ataagagaca 60
gt 62
Claims (39)
1.一种测序结果分析方法，其特征在于，
所述测序结果包括第一测序数据和第二测序数据，其中，所述第一测序数据和所述第二测序数据均由多个读段构成，所述第一测序数据中的至少一部分所述读段在所述第二测序数据中存在对应读段，
所述测序结果分析方法包括：
(a)基于所述第一测序数据和所述第二测序数据各自的至少一部分进行相互校正，以便获得最终序列信息；
所述相互校正包括下列步骤：
在所述第一测序数据和所述第二测序数据中选择高质量读段和所述高质量读段的对应读段，所述读段的长度不低于预定长度，所述读段具有不低于预定质量阈值的测序质量；和
将所述高质量读段与所述高质量读段的对应读段进行比对，并基于所述比对结果进行序列信息校正；
其中，对芯片表面的测序模板进行第一测序，以便通过形成第一新生测序链获得第一测序数据，所述测序模板通过测序接头连接在所述芯片表面上；
对所述测序模板进行第二测序，以便通过形成第二新生测序链获得第二测序数据。
2.根据权利要求1所述的方法，其特征在于，步骤(a)进一步包括：
(a-1)根据所述读段的长度，基于所述第一测序数据，构建第一读段集合，所述第一读段集合中的每一个读段长度均不低于第一预定长度；
(a-2)根据所述对应读段的长度，基于所述第一读段集合，构建第二读段集合和第三读段集合，所述第二读段集合中每一个读段的所述对应读段的长度均不低于第二预定长度，所述第三读段集合中每一个读段的所述对应读段的长度均处于预定长度范围内；
(a-3)根据所述第二读段集合中所述读段及其所述对应读段的测序质量，基于所述第二读段集合及其所述对应读段，构建第四读段集合和第五读段集合，其中，所述第四读段集合和所述第五读段集合分别是按照下列原则确定的：
将所述第二读段集合中的所述读段与其所述对应读段进行测序质量比较，
选择测序质量高的一方作为所述第四读段集合的元素，选择测序质量低的一方作为所述第五读段集合的元素，
对于测序质量相同的情形，则选择来自所述第二读段集合的所述读段作为所述第四读段集合的元素，则选择所述对应读段作为所述第五读段集合的元素；
(a-4)利用测序质量，对所述第四读段集合进行过滤处理，以便构建第六读段集合，所述第六读段集合中的所述读段的测序质量均不低于第一预定质量阈值；
(a-5)利用所述第六读段集合，从所述第五读段集合中选择与所述第六读段集合中的所述读段对应的所述读段，以便构建第七读段集合；
(a-6)将所述第六读段集合与所述第七读段集合进行读段比对，并在所述第六读段集合的所述读段上确定第一差异位点；和
(a-7)利用预先确定的测序误差预测模型，对所述第一差异位点进行校正，以便确定第一序列信息，所述测序误差预测模型用于确定在测序过程中，差异位点发生插入或者缺失的概率。
3.根据权利要求2所述的方法，其特征在于，进一步包括：
(a-4a)利用测序质量，对所述第三读段集合进行过滤处理，以便构建第八读段集合，其中，所述第八读段集合中的所述读段的测序质量均不低于第二预定质量阈值；
(a-5a)利用所述第八读段集合，从所述第二测序数据中选择与所述第七读段集合中的所述读段对应的所述读段，以便构建第九读段集合；
(a-6a)将所述第八读段集合与所述第九读段集合进行读段比对，并在所述第八读段集合的所述读段上确定第二差异位点；
(a-7a)利用所述测序误差预测模型，对所述第二差异位点进行校正，以便确定第二序列信息。
4.根据权利要求2或3所述的方法，其特征在于，所述测序误差预测模型是基于所述第一测序数据和所述第二测序数据与参考基因组的比对结果，对朴素贝叶斯模型进行训练获得的。
5.根据权利要求3所述的方法，其特征在于，针对所述第一差异位点和所述第二差异位点：
如果来自所述第六读段集合的读段在所述差异位点存在碱基，来自所述第七读段集合的对应读段在所述差异位点不存在碱基，并且在所述差异位点发生缺失的概率为50％以上，则保留所述第六读段集合的读段在所述差异位点的碱基作为最终测序结果；
如果来自所述第六读段集合的读段在所述差异位点不存在碱基，来自所述第七读段的读段集合的定读段在所述差异位点存在碱基，并且在所述差异位点发生插入的概率为50％以上，则保留所述第六读段集合的读段在所述差异位点的碱基作为最终测序结果；和
如果来自所述第六读段集合的读段在所述差异位点存在碱基，来自所述第七读段的读段集合的定读段在所述差异位点也存在碱基，则选择所述第六读段集合的读段在所述差异位点的碱基作为最终测序结果。
6.根据权利要求2所述的方法，其特征在于，所述预定长度包括第一预定长度和第二预定长度；所述预定长度范围为10～25bp。
7.根据权利要求6所述的方法，其特征在于，所述第一预定长度和所述第二预定长度分别独立地不低于20bp。
8.根据权利要求7所述的方法，其特征在于，所述第一预定长度和所述第二预定长度分别独立地不低于25bp。
9.根据权利要求1所述的方法，其特征在于，所述预定质量阈值包括第一预定质量阈值和第二预定质量阈值，所述第一预定质量阈值和所述第二预定质量阈值分别独立地不低于50。
10.根据权利要求9所述的方法，其特征在于，所述第一预定质量阈值和所述第二预定质量阈值分别独立地不低于60。
11.根据权利要求1-3、5-10任一项所述的方法，其特征在于，在获得第一测序数据后，在进行所述第二测序前，包括：去除所述芯片表面的所述第一新生测序链，并对残余在所述芯片表面的所述第一新生测序链的3’-末端进行第一封闭处理。
12.根据权利要求11所述的方法，其特征在于，所述第一封闭处理是通过使3’-末端羟基与延伸反应阻断剂相连而进行的。
13.根据权利要求12所述的方法，其特征在于，所述延伸反应阻断剂为ddNTP或其衍生物。
14.根据权利要求13所述的方法，其特征在于，所述第一封闭处理是采用DNA聚合酶和末端转移酶的至少之一进行的。
15.根据权利要求14所述的方法，其特征在于，所述第一封闭处理通过聚合酶连接所述ddNTP或其衍生物。
16.根据权利要求1-3、5-10任一项所述的方法，其特征在于，在获得第一测序数据前，包括：
(1-a)使测序文库中的文库分子与芯片表面的测序接头进行杂交；
(1-b)利用所述文库分子作为初始模板，通过合成互补链形成所述测序模板；
(1-c)除去所述初始模板，并对所述芯片表面的核酸分子的3’-末端进行第二封闭处理。
17.根据权利要求16所述的方法，其特征在于，所述第二封闭处理是通过使3’-末端羟基与延伸反应阻断剂相连而进行的。
18.根据权利要求17所述的方法，其特征在于，所述延伸反应阻断剂为ddNTP或其衍生物。
19.根据权利要求18所述的方法，其特征在于，所述第二封闭处理是采用DNA聚合酶和末端转移酶的至少之一进行的。
20.根据权利要求19所述的方法，其特征在于，所述第二封闭处理通过所述末端转移酶连接所述ddNTP或其衍生物。
21.根据权利要求16所述的方法，其特征在于，在(1-c)之前，进一步包括：
(1-b-1)对步骤(1-b)中延伸不完全的所述互补链的3’-末端进行第三封闭处理。
22.根据权利要求21所述的方法，其特征在于，所述第三封闭处理是通过使3’-末端羟基与延伸反应阻断剂相连而进行的。
23.根据权利要求22所述的方法，其特征在于，所述延伸反应阻断剂为ddNTP或其衍生物。
24.根据权利要求23所述的方法，其特征在于，所述第三封闭处理是采用DNA聚合酶和末端转移酶的至少之一进行的。
25.根据权利要求24所述的方法，其特征在于，所述第三封闭处理通过聚合酶连接所述ddNTP或其衍生物。
26.一种测序结果分析系统，其特征在于，包括:
测序设备，所述测序设备适于通过双重测序法获得测序结果，所述测序结果包括第一测序数据和第二测序数据，所述第一测序数据和所述第二测序数据均由多个读段构成，所述第一测序数据中的至少一部分所述读段在所述第二测序数据中存在对应读段；
分析设备，所述分析设备包括校正模块，适用于基于所述第一测序数据和所述第二测序数据各自的至少一部分进行相互校正，以便获得最终序列信息；
所述校正模块适用于以下步骤：
在所述第一测序数据和所述第二测序数据中选择高质量读段和所述高质量读段的对应读段，所述读段的长度不低于预定长度，所述读段具有不低于预定质量阈值的测序质量；和
将所述高质量读段与所述高质量读段的对应读段进行比对，并基于所述比对结果进行序列信息校正；
其中，对芯片表面的测序模板进行第一测序，以便通过形成第一新生测序链获得第一测序数据，所述测序模板通过测序接头连接在所述芯片表面上；
对所述测序模板进行第二测序，以便通过形成第二新生测序链获得第二测序数据。
27.根据权利要求26所述的系统，其特征在于，所述校正模块进一步包括：
第一读段集合确定单元，根据所述读段的长度，基于所述第一测序数据，构建第一读段集合，所述第一读段集合中的每一个读段长度均不低于第一预定长度；
第二读段集合和第三读段集合确定单元，根据所述对应读段的长度，基于所述第一读段集合，构建第二读段集合和第三读段集合，所述第二读段集合中每一个读段的所述对应读段的长度均不低于第二预定长度，所述第三读段集合中每一个读段的所述对应读段的长度均处于预定长度范围内；
第四读段集合和第五读段集合确定单元，根据所述第二读段集合中所述读段及其所述对应读段的测序质量，基于所述第二读段集合及其所述对应读段，构建第四读段集合和第五读段集合，其中，所述第四读段集合和所述第五读段集合分别是按照下列原则确定的：
将所述第二读段集合中的所述读段与其所述对应读段进行测序质量比较，
选择测序质量高的一方作为所述第四读段集合的元素，选择测序质量低的一方作为所述第五读段集合的元素，
对于测序质量相同的情形，则选择来自所述第二读段集合的所述读段作为所述第四读段集合的元素，则选择所述对应读段作为所述第五读段集合的元素；
第六读段集合确定单元，利用测序质量，对所述第四读段集合进行过滤处理，以便构建第六读段集合，所述第六读段集合中的所述读段的测序质量均不低于第一预定质量阈值；
第七读段集合确定单元，利用所述第六读段集合，从所述第五读段集合中选择与所述第六读段集合中的所述读段对应的所述读段，以便构建第七读段集合；
第一差异位点确定单元，将所述第六读段集合与所述第七读段集合进行读段比对，并在所述第六读段集合的所述读段上确定第一差异位点；和
第一序列信息确定单元，利用预先确定的测序误差预测模型，对所述第一差异位点进行校正，以便确定第一序列信息，所述测序误差预测模型用于确定在测序过程中，差异位点发生插入或者缺失的概率。
28.根据权利要求27所述的系统，其特征在于，所述校正模块进一步包括：
第八读段集合确定单元，利用测序质量，对所述第三读段集合进行过滤处理，以便构建第八读段集合，其中，所述第八读段集合中的所述读段的测序质量均不低于第二预定质量阈值；
第九读段集合确定单元，利用所述第八读段集合，从所述第二测序数据中选择与所述第七读段集合中的所述读段对应的所述读段，以便构建第九读段集合；
第二差异位点确定单元，将所述第八读段集合与所述第九读段集合进行读段比对，并在所述第八读段集合的所述读段上确定第二差异位点；
第二序列信息确定单元，利用所述测序误差预测模型，对所述第二差异位点进行校正，以便确定第二序列信息。
29.根据权利要求27所述的系统，其特征在于，所述第一预定长度和所述第二预定长度分别独立地不低于20bp。
30.根据权利要求29所述的系统，其特征在于，所述第一预定长度和所述第二预定长度分别独立地不低于25bp。
31.根据权利要求27、28或29任一项所述的系统，其特征在于，所述预定长度范围为10～25bp。
32.根据权利要求27所述的系统，其特征在于，所述第一预定质量阈值不低于50。
33.根据权利要求28所述的系统，其特征在于，所述第二预定质量阈值不低于50。
34.根据权利要求32所述的系统，其特征在于，所述第一预定质量阈值不低于60。
35.根据权利要求33所述的系统，其特征在于，所述第二预定质量阈值不低于60。
36.根据权利要求26所述的系统，其特征在于，所述测序结果分析系统进一步包括：测序误差预测模型构建模块，所述测序误差预测模型构建模块适用于基于所述第一测序数据和所述第二测序数据与参考基因组的比对结果，对朴素贝叶斯模型进行训练，以便获得所述测序误差预测模型。
37.根据权利要求28所述的系统，其特征在于，针对所述第一差异位点和所述第二差异位点：
如果来自所述第六读段集合的读段在所述差异位点存在碱基，来自所述第七读段集合的对应读段在所述差异位点不存在碱基，并且在所述差异位点发生缺失的概率为50％以上，则保留所述第六读段集合的读段在所述差异位点的碱基作为最终测序结果；
如果来自所述第六读段集合的读段在所述差异位点不存在碱基，来自所述第七读段的读段集合的定读段在所述差异位点存在碱基，并且在所述差异位点发生插入的概率为50％以上，则保留所述第六读段集合的读段在所述差异位点的碱基作为最终测序结果；和
如果来自所述第六读段集合的读段在所述差异位点存在碱基，来自所述第七读段的读段集合的定读段在所述差异位点也存在碱基，则选择所述第六读段集合的读段在所述差异位点的碱基作为最终测序结果。
38.一种计算机可读存储介质，其上存储有计算机程序，其特征在于，该程序被处理器执行时实现权利要求1～25中任一项所述方法的步骤。
39.一种电子设备，其特征在于，包括：
权利要求38中所述的计算机可读存储介质；以及
一个或者多个处理器，用于执行所述计算机可读存储介质中的程序。
Priority Applications (3)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| EP21797265.2A EP4144745A4 (en) | 2020-04-30 | 2021-04-30 | SEQUENCING METHOD, ANALYSIS METHOD AND ANALYSIS SYSTEM, COMPUTER READABLE STORAGE MEDIUM AND ELECTRONIC DEVICE |
| PCT/CN2021/091279 WO2021219114A1 (zh) | 2020-04-30 | 2021-04-30 | 测序方法及其分析方法和系统、计算机可读存储介质和电子设备 |
| US17/922,340 US20230178183A1 (en) | 2020-04-30 | 2021-04-30 | Sequencing method, analysis method therefor and analysis system thereof, computer-readable storage medium, and electronic device |
Applications Claiming Priority (2)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN2020103625876 | 2020-04-30 | ||
| CN202010362587 | 2020-04-30 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113593636A CN113593636A (zh) | 2021-11-02 |
| CN113593636B true CN113593636B (zh) | 2024-05-03 |
Family
ID=77467752
Family Applications (2)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010865293.5A Active CN113593636B (zh) | 2020-04-30 | 2020-08-25 | 测序结果分析方法、系统及计算机可读存储介质和电子设备 |
| CN202110264389.0A Pending CN113337576A (zh) | 2020-04-30 | 2021-03-11 | 文库制备方法、试剂盒及测序方法 |
Family Applications After (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110264389.0A Pending CN113337576A (zh) | 2020-04-30 | 2021-03-11 | 文库制备方法、试剂盒及测序方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (2) | CN113593636B (zh) |
Families Citing this family (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2023066255A1 (zh) * | 2021-10-18 | 2023-04-27 | 深圳市真迈生物科技有限公司 | 测序方法、测序数据处理方法、设备和计算机设备 |
| CN114134204B (zh) * | 2021-12-07 | 2023-03-21 | 上海交通大学 | 一种用于微量标记核酸样本检测和测序的样品制备方法 |
| CN115747301B (zh) * | 2022-08-01 | 2023-12-22 | 深圳赛陆医疗科技有限公司 | 一种测序文库的构建方法、构建测序文库的试剂盒及基因测序方法 |
| WO2024093421A1 (zh) * | 2022-10-31 | 2024-05-10 | 深圳市真迈生物科技有限公司 | 表面处理方法、测序方法和试剂盒 |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107077538A (zh) * | 2014-12-10 | 2017-08-18 | 深圳华大基因研究院 | 测序数据处理装置和方法 |
| CN110021351A (zh) * | 2018-07-19 | 2019-07-16 | 深圳华大生命科学研究院 | 分析碱基连锁强度以及基因分型方法和系统 |
| CN110520542A (zh) * | 2017-03-23 | 2019-11-29 | 华盛顿大学 | 用于靶向核酸序列富集的方法及在错误纠正的核酸测序中的应用 |
Family Cites Families (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| EP2483428A1 (en) * | 2009-09-29 | 2012-08-08 | Agency For Science, Technology And Research | Methods and arrays for dna sequencing |
| EP3434776A1 (en) * | 2012-12-12 | 2019-01-30 | The Broad Institute, Inc. | Methods, models, systems, and apparatus for identifying target sequences for cas enzymes or crispr-cas systems for target sequences and conveying results thereof |
| US10102337B2 (en) * | 2014-08-06 | 2018-10-16 | Nugen Technologies, Inc. | Digital measurements from targeted sequencing |
| CN108165618B (zh) * | 2017-12-08 | 2021-06-08 | 东南大学 | 一种包含核苷酸和3’端可逆封闭核苷酸的dna测序方法 |
| CN108611346A (zh) * | 2018-05-06 | 2018-10-02 | 湖南大地同年生物科技有限公司 | 一种单细胞基因表达量检测文库的构建方法 |
-
2020
- 2020-08-25 CN CN202010865293.5A patent/CN113593636B/zh active Active
-
2021
- 2021-03-11 CN CN202110264389.0A patent/CN113337576A/zh active Pending
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107077538A (zh) * | 2014-12-10 | 2017-08-18 | 深圳华大基因研究院 | 测序数据处理装置和方法 |
| CN110520542A (zh) * | 2017-03-23 | 2019-11-29 | 华盛顿大学 | 用于靶向核酸序列富集的方法及在错误纠正的核酸测序中的应用 |
| CN110021351A (zh) * | 2018-07-19 | 2019-07-16 | 深圳华大生命科学研究院 | 分析碱基连锁强度以及基因分型方法和系统 |
Non-Patent Citations (2)
| Title |
|---|
| "肝细胞癌基因变异分析及循环肿瘤DNA检测对肝细胞癌的诊断价值";熊宇;《中国博士学位论文全文数据库医药卫生科技辑》;全文 * |
| Higgins Jacob E 等."Duplex Sequencing Accurately Detects Variants below 1/100,000 in Genes Recurrently Mutated in Acute Myeloid Leukemia".《Blood》.2018,全文. * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113593636A (zh) | 2021-11-02 |
| CN113337576A (zh) | 2021-09-03 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN113593636B (zh) | 测序结果分析方法、系统及计算机可读存储介质和电子设备 | |
| US20200354773A1 (en) | High multiplex pcr with molecular barcoding | |
| CN106795514B (zh) | 泡状接头及其在核酸文库构建及测序中的应用 | |
| US20220251549A1 (en) | Method for constructing capture library and kit | |
| CN106086162A (zh) | 一种用于检测肿瘤突变的双标签接头序列及检测方法 | |
| JP2022510723A (ja) | 遺伝子標的エリアの富化方法及びキット | |
| CN107604046B (zh) | 用于微量dna超低频突变检测的双分子自校验文库制备及杂交捕获的二代测序方法 | |
| CN109593757B (zh) | 一种探针及其适用于高通量测序的对目标区域进行富集的方法 | |
| CN111041563B (zh) | 一种靶向序列捕获及pcr建库方法 | |
| CN110396516B (zh) | 一种基于特有识别序列的绝对定量转录组文库构建方法 | |
| CN112251821A (zh) | 一种快速高效的构建二代测序文库的试剂盒 | |
| EP3198063A1 (en) | Rna stitch sequencing: an assay for direct mapping of rna : rna interactions in cells | |
| CN110468211B (zh) | 膀胱癌肿瘤突变基因特异性引物、试剂盒和文库构建方法 | |
| CN113593637B (zh) | 测序方法及其分析方法和系统、计算机可读存储介质和电子设备 | |
| CN112941635A (zh) | 一种提高文库转化率的二代测序建库试剂盒及其方法 | |
| CN112301430B (zh) | 一种建库方法及应用 | |
| CN114807300A (zh) | 单引物多重扩增技术在检测片段化稀有特征核酸分子中的应用及试剂盒 | |
| US20180100180A1 (en) | Methods of single dna/rna molecule counting | |
| WO2023202030A1 (zh) | 一种小分子rna的高通量测序文库构建方法 | |
| WO2023066255A1 (zh) | 测序方法、测序数据处理方法、设备和计算机设备 | |
| CN114277114B (zh) | 一种扩增子测序添加唯一性标识符的方法及应用 | |
| CN116536308A (zh) | 测序封闭剂及其应用 | |
| CN113584135B (zh) | 一种混样检测rna修饰并实现精准定量的方法 | |
| CN114807324A (zh) | 单引物扩增建库技术在检测片段化稀有dna分子突变中的应用及试剂盒 | |
| WO2020135650A1 (zh) | 一种基因测序文库的构建方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |