CN116064755B - 一种基于连锁基因突变检测mrd标志物的装置 - Google Patents
一种基于连锁基因突变检测mrd标志物的装置 Download PDFInfo
- Publication number
- CN116064755B CN116064755B CN202310063877.4A CN202310063877A CN116064755B CN 116064755 B CN116064755 B CN 116064755B CN 202310063877 A CN202310063877 A CN 202310063877A CN 116064755 B CN116064755 B CN 116064755B
- Authority
- CN
- China
- Prior art keywords
- mutation
- snv
- linkage
- mutations
- sequencing
- 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
- 206010064571 Gene mutation Diseases 0.000 title claims abstract description 11
- 239000003550 marker Substances 0.000 title abstract description 6
- 230000035772 mutation Effects 0.000 claims abstract description 494
- 238000012163 sequencing technique Methods 0.000 claims abstract description 110
- 238000012544 monitoring process Methods 0.000 claims abstract description 80
- 238000000034 method Methods 0.000 claims abstract description 69
- 239000012634 fragment Substances 0.000 claims abstract description 66
- 238000001514 detection method Methods 0.000 claims abstract description 57
- 238000010276 construction Methods 0.000 claims abstract description 21
- 210000004369 blood Anatomy 0.000 claims abstract description 12
- 239000008280 blood Substances 0.000 claims abstract description 12
- 239000000523 sample Substances 0.000 claims description 141
- 206010028980 Neoplasm Diseases 0.000 claims description 66
- 238000004364 calculation method Methods 0.000 claims description 47
- VYZAMTAEIAYCRO-UHFFFAOYSA-N Chromium Chemical compound [Cr] VYZAMTAEIAYCRO-UHFFFAOYSA-N 0.000 claims description 46
- 239000011159 matrix material Substances 0.000 claims description 44
- 230000036438 mutation frequency Effects 0.000 claims description 44
- 238000003745 diagnosis Methods 0.000 claims description 36
- 238000001914 filtration Methods 0.000 claims description 29
- 210000000349 chromosome Anatomy 0.000 claims description 24
- 238000012216 screening Methods 0.000 claims description 21
- 238000011282 treatment Methods 0.000 claims description 17
- 210000004602 germ cell Anatomy 0.000 claims description 14
- 230000001717 pathogenic effect Effects 0.000 claims description 12
- 238000012360 testing method Methods 0.000 claims description 12
- 238000002372 labelling Methods 0.000 claims description 11
- 230000002159 abnormal effect Effects 0.000 claims description 10
- 238000012545 processing Methods 0.000 claims description 10
- 230000008859 change Effects 0.000 claims description 9
- 238000010998 test method Methods 0.000 claims description 9
- 238000006243 chemical reaction Methods 0.000 claims description 7
- 238000005070 sampling Methods 0.000 claims description 7
- 238000003780 insertion Methods 0.000 claims description 6
- 230000037431 insertion Effects 0.000 claims description 6
- FZNCGRZWXLXZSZ-CIQUZCHMSA-N Voglibose Chemical compound OCC(CO)N[C@H]1C[C@](O)(CO)[C@@H](O)[C@H](O)[C@H]1O FZNCGRZWXLXZSZ-CIQUZCHMSA-N 0.000 claims description 4
- 230000002759 chromosomal effect Effects 0.000 claims description 4
- 230000003902 lesion Effects 0.000 claims description 4
- 210000001161 mammalian embryo Anatomy 0.000 claims description 4
- 239000002773 nucleotide Substances 0.000 claims description 2
- 125000003729 nucleotide group Chemical group 0.000 claims 1
- 238000005516 engineering process Methods 0.000 abstract description 16
- 108090000623 proteins and genes Proteins 0.000 abstract description 10
- 238000007796 conventional method Methods 0.000 abstract description 6
- 238000002474 experimental method Methods 0.000 abstract description 6
- 238000000684 flow cytometry Methods 0.000 abstract description 6
- 238000012806 monitoring device Methods 0.000 abstract description 4
- 208000027418 Wounds and injury Diseases 0.000 abstract description 3
- 230000006378 damage Effects 0.000 abstract description 3
- 238000011156 evaluation Methods 0.000 abstract description 3
- 208000014674 injury Diseases 0.000 abstract description 3
- 230000002980 postoperative effect Effects 0.000 abstract description 3
- 238000004445 quantitative analysis Methods 0.000 abstract description 3
- 230000002285 radioactive effect Effects 0.000 abstract description 3
- 210000004027 cell Anatomy 0.000 description 16
- 210000001519 tissue Anatomy 0.000 description 14
- 238000004422 calculation algorithm Methods 0.000 description 13
- 210000002381 plasma Anatomy 0.000 description 10
- 210000005259 peripheral blood Anatomy 0.000 description 9
- 239000011886 peripheral blood Substances 0.000 description 9
- 230000008569 process Effects 0.000 description 9
- 108091081062 Repeated sequence (DNA) Proteins 0.000 description 8
- 238000012070 whole genome sequencing analysis Methods 0.000 description 8
- 108020004414 DNA Proteins 0.000 description 7
- 238000004891 communication Methods 0.000 description 5
- 238000002156 mixing Methods 0.000 description 5
- 238000007481 next generation sequencing Methods 0.000 description 5
- 238000012408 PCR amplification Methods 0.000 description 4
- 201000011510 cancer Diseases 0.000 description 4
- 238000012937 correction Methods 0.000 description 4
- 238000013461 design Methods 0.000 description 4
- 230000000694 effects Effects 0.000 description 4
- 238000000605 extraction Methods 0.000 description 4
- 238000003908 quality control method Methods 0.000 description 4
- 230000009471 action Effects 0.000 description 3
- 238000012417 linear regression Methods 0.000 description 3
- 238000011160 research Methods 0.000 description 3
- 230000035945 sensitivity Effects 0.000 description 3
- 238000004088 simulation Methods 0.000 description 3
- 108091092584 GDNA Proteins 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 2
- 238000001574 biopsy Methods 0.000 description 2
- 238000007847 digital PCR Methods 0.000 description 2
- 238000007689 inspection Methods 0.000 description 2
- 238000000053 physical method Methods 0.000 description 2
- 238000000611 regression analysis Methods 0.000 description 2
- 230000000630 rising effect Effects 0.000 description 2
- 238000001356 surgical procedure Methods 0.000 description 2
- 206010065163 Clonal evolution Diseases 0.000 description 1
- WQZGKKKJIJFFOK-GASJEMHNSA-N Glucose Natural products OC[C@H]1OC(O)[C@H](O)[C@@H](O)[C@@H]1O WQZGKKKJIJFFOK-GASJEMHNSA-N 0.000 description 1
- 206010025323 Lymphomas Diseases 0.000 description 1
- 238000012614 Monte-Carlo sampling Methods 0.000 description 1
- 208000007660 Residual Neoplasm Diseases 0.000 description 1
- 230000032683 aging Effects 0.000 description 1
- 210000000601 blood cell Anatomy 0.000 description 1
- 238000002591 computed tomography Methods 0.000 description 1
- 238000004590 computer program Methods 0.000 description 1
- 239000000470 constituent Substances 0.000 description 1
- 238000013500 data storage Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 201000010099 disease Diseases 0.000 description 1
- 208000037265 diseases, disorders, signs and symptoms Diseases 0.000 description 1
- 239000008103 glucose Substances 0.000 description 1
- 230000036541 health Effects 0.000 description 1
- 201000001421 hyperglycemia Diseases 0.000 description 1
- 208000015181 infectious disease Diseases 0.000 description 1
- 230000002458 infectious effect Effects 0.000 description 1
- 238000007726 management method Methods 0.000 description 1
- 238000013507 mapping Methods 0.000 description 1
- 201000001441 melanoma Diseases 0.000 description 1
- 239000002207 metabolite Substances 0.000 description 1
- 239000000203 mixture Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 239000003147 molecular marker Substances 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 230000008520 organization Effects 0.000 description 1
- 230000003647 oxidation Effects 0.000 description 1
- 238000007254 oxidation reaction Methods 0.000 description 1
- 238000011084 recovery Methods 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 229920006395 saturated elastomer Polymers 0.000 description 1
- 238000000926 separation method Methods 0.000 description 1
- 230000000392 somatic effect Effects 0.000 description 1
- 210000004881 tumor cell Anatomy 0.000 description 1
- 238000002604 ultrasonography 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/50—Mutagenesis
-
- 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/6869—Methods for sequencing
-
- 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/6876—Nucleic acid products used in the analysis of nucleic acids, e.g. primers or probes
- C12Q1/6883—Nucleic acid products used in the analysis of nucleic acids, e.g. primers or probes for diseases caused by alterations of genetic material
- C12Q1/6886—Nucleic acid products used in the analysis of nucleic acids, e.g. primers or probes for diseases caused by alterations of genetic material for cancer
-
- 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
- G16B50/00—ICT programming tools or database systems specially adapted for bioinformatics
- G16B50/30—Data warehousing; Computing architectures
-
- 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
- C12Q2600/00—Oligonucleotides characterized by their use
- C12Q2600/118—Prognosis of disease development
-
- 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
- C12Q2600/00—Oligonucleotides characterized by their use
- C12Q2600/156—Polymorphic or mutational markers
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02D—CLIMATE CHANGE MITIGATION TECHNOLOGIES IN INFORMATION AND COMMUNICATION TECHNOLOGIES [ICT], I.E. INFORMATION AND COMMUNICATION TECHNOLOGIES AIMING AT THE REDUCTION OF THEIR OWN ENERGY USE
- Y02D10/00—Energy efficient computing, e.g. low power processors, power management or thermal management
Landscapes
- Life Sciences & Earth Sciences (AREA)
- Health & Medical Sciences (AREA)
- Chemical & Material Sciences (AREA)
- Engineering & Computer Science (AREA)
- Proteomics, Peptides & Aminoacids (AREA)
- Organic Chemistry (AREA)
- Physics & Mathematics (AREA)
- Bioinformatics & Cheminformatics (AREA)
- General Health & Medical Sciences (AREA)
- Zoology (AREA)
- Genetics & Genomics (AREA)
- Wood Science & Technology (AREA)
- Biotechnology (AREA)
- Analytical Chemistry (AREA)
- Biophysics (AREA)
- Theoretical Computer Science (AREA)
- Immunology (AREA)
- Molecular Biology (AREA)
- General Engineering & Computer Science (AREA)
- Spectroscopy & Molecular Physics (AREA)
- Pathology (AREA)
- Medical Informatics (AREA)
- Microbiology (AREA)
- Evolutionary Biology (AREA)
- Bioinformatics & Computational Biology (AREA)
- Biochemistry (AREA)
- Databases & Information Systems (AREA)
- Hospice & Palliative Care (AREA)
- Oncology (AREA)
- Bioethics (AREA)
- Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)
Abstract
本申请涉及生物基因技术领域，尤其涉及一种基于连锁基因突变检测MRD标志物的装置。本发明建立了一种基于连锁突变的MRD监测装置。相比于其他流式细胞方法和PET‑CET方法，本发明方法样本获取简单，仅需要患者血液，操作方便，无放射性伤害；相比于PCR定量方法，本发明能一次检测多个位点，并且检测下限可达10―6级别，可有效地进行MRD术后评价和复发监控，随后根据患者个体情况开展精准指导；本发明不依赖于UMI分子标签技术，当仅使用常规方法进行cfDNA建库时，本发明检测下限也可达到10―5级别，相对于双重测序技术，本发明节约实验建库成本，能有效提高原始片段的利用率。
Description
技术领域
本申请涉及生物基因技术领域，尤其涉及一种基于连锁基因突变检测MRD标志物的装置。
背景技术
根据国家癌症中心发布的最新统计数据，2016年我国癌症新发病例为406.4万例，死亡病例为241.4万例。而根据最新世界卫生组织国际癌症研究机构，2020年中国新发癌症457万人，死亡病例300多万。对比已有统计数据，我国整体癌症粗发病率仍处于持续上升趋势，粗死亡率也是仍然处于上升趋势。并且随着中国人口老龄化加剧，中国癌症负担将逐渐加大，加大对癌症预防和治疗的研究投入，将有助于降低未来的癌症负担。死亡癌症患者的绝大部分是晚期患者，晚期肿瘤患者死亡率较高的原因，大多数并非是手术或者治疗不成功，而是经过各种综合治疗之后，体内仍残余少部分肿瘤细胞，这部分细胞定位为微小残留病灶(MRD，minimal residual disease)，并且随着时间推移，这部分细胞残余会发展成新的病灶。
通过技术手段对MRD进行监测，能够监测肿瘤复发可能性，根据病人的风险程度进行分层管理，定制精准化治疗方案。MRD的技术监测手段包含正电子发射断层成像-X线计算机断层成像仪技术(PET-CT)、流式细胞技术、数字PCR技术和NGS技术等等。PET-CT技术通过监测细胞代谢物葡萄糖进行判断肿瘤组织判断，但是该技术具有放射性，并且针对一些感染性病变、血糖过高等情况存在假阳性。流式细胞技术运用不同荧光标记物的多种抗体组合对细胞表面表达情况进行检测，进而区分细胞来源，然而活细胞的提取限制了流式细胞在MRD的应用。数字PCR技术对异常的DNA或者RNA片段进行扩增，进而定量。由于成本限制，该技术仅能对1个或少数几个突变进行监测，故存在假阴性可能。通过对外周血游离DNA(circulating free DNA，cfDNA)进行NGS检测，对患者来说无需通过手术进行组织提取，能够一次检测多个位点，大大提高检测的灵敏度，并且对于肿瘤的克隆演化也有较好的监控作用。
通过cfDNA测序进行MRD监测，主要有两个关键制约因素：1.经过综合治疗之后，血浆中循环肿瘤DNA(circulating tumor DNA，ctDNA)剂量极少；2.存在较多的测序实验过程中引入的背景噪音。目前大部分通过分子标记测序技术(unique molecular identifiers，UMI)来降低背景噪音，但是该方法需要测序较为饱和来实现校正算法，并且该方法受限于DNA双链分子回收率，仅能有效利用到20至25％的原始片段。此外，目前该技术的主要应用方法，均为基于多个突变位点中每一个突变位点的真假进行独立判断，该策略将会大大影响检测的灵敏度，目前基于NGS测序进行单个突变位点的阳性判定的检测下限基本在0.1％至0.5％之间。
综上所述，现有技术中阳性判定的检测下限偏低，影响了检测结果容易造成误判。
发明内容
为了解决背景技术中存在的问题，本发明建立了一种利用连锁突变计算肿瘤MRD的方法及对应的数据存储计算装置，通过患者初诊时的任何检测手段(如靶向测序、全外显子测序或者全基因组测序等)所获得的初诊连锁突变作为基线，进行超高灵敏度的MRD监测。具体内容如下：
第一方面，本发明提供了一种连锁基因突变数据库构建方法，包括：
提取若干个正常人的cfDNA进行双端测序，获取测序数据作为样本数据，使用hg19参考基因组进行对比，获得比对信息，并标记重复片段；
针对每个样本数据，并行计算该样本每个位点的每类碱基突变频率；
针对每个样本数据，并行计算该样本所有可能的LV信息；
收集每个样本数据的SNV及LV信息，构建SNV和LV数据库；
其中SNV以SNV突变位置和突变类型为索引，以样本编号为列名，存储每个样本对应的SNV突变支持数、测序深度以及SNV突变频率；LV以LV突变位置和突变类型为索引，以样本编号为列名，存储每个样本对应的LV突变支持数、测序深度以及LV突变频率。
进一步的，所述针对每个样本数据，并行计算该样本每个位点的每类突变频率，具体包括：
将探针目标区域延伸340bp，使用本地计算装置按照目标延伸区域多线程并行计算每个位点的测序深度(depth，DP)以及4类碱基(ATCG)对应的数目(allele count，AD)，按照read(测序读段)比对质量(mapping quality，MQ)和碱基质量(base quantity，BQ)进行过滤，按照MQ≥20，BQ≥20设置，然后得出每个样本在该位置对应4类碱基的突变频率AF＝AD/DP。
进一步的，所述针对每个样本数据，并行计算该样本所有可能的LV信息，具体包括：
A.过滤soft-clip(软切除)的read，按照MQ≥20筛选；
B.通过read上MD标签，获得每个read上的错配碱基数目，如果read上错配碱基数量≥7个，则过滤该read；
C.将双端read合并为一个片段，如果Read1和Read2的重叠部分碱基不同，则修改为N碱基，如果两个read没有重叠，则将空白区域标记为N碱基；
D.标记合并片段上的SNV突变碱基，并排除N碱基，如果一个片段上≥2个alt，则认为该片段包含一个可能LV突变；
E.将所有存在LV突变的片段按照首个SNV的突变位置进行排序，迭代所有片段，获得以LV突变为索引，LV突变片段数目为列的LV突变矩阵；
F.合并所有LV突变，构成连锁突变区域，如果某一LV突变i，其染色体编号为chrom0，起始的SNV突变位置是start0,，终止的SNV突变位置是end0，另一LV突变j，其染色体编号是chrom1，起始SNV突变位置是start1，终止SNV突变位置是end1，突变j满足如下a)或者b)条件：
a)chrom0＝chrom1并且start1≥start0并且start1≤end0
b)chrom0＝chrom1并且start0≥start1并且start0≤end1
则将连锁突变i和连锁突变j合并到一个连锁突变区域；
G.以一个连锁突变区域为计算单元，针对该区域完成如下计算：
a)假设该连锁区域包含n个LV连锁突变，每个LV连锁突变，包含多个SNV单元，即为集合Ai，对应的片段支持数为ADi；
b)针对SNV单元数目大于等于3的LV突变，假设为该SNV单元数目为x，从中随机抽取2个SNV单元构成2连锁LV，则包含个组合，从中随机抽取3个构成3连锁LV，则包含个组合，最终将1个x连锁的LV突变拆成个3连锁或者2连锁的LV突变；
c)将SNV单元数目等于2的LV突变和上述拆分所有3连锁或者2连锁的LV突变合并成一个包含m个2连锁或者3连锁的LV突变集合M，集合中每个LV突变是包含2到3个SNV突变单元(即为集合Bj)，则对应的片段支持数ADj按照如下公式进行计算：
d)上述每个LV突变Bj，假设SNV的坐标位置为{POS1,POS2,…,POSj}(即为集合Pk)，具有相同的突变位置LV合并，然后构成k个连锁突变位置集合K；
e)假设该连锁突变区域，起始位置为s1，终止位置为sq，则将集合K，转变为k×q的01矩阵Mat：
f)同时构建一个k的0向量VecDP，用以存储集合K中每个元素对应的深度信息DPk；
g)迭代该连锁区域的所有read，过滤掉标签为重复、补充比对、次比对、比对异常、MQ≤20和错配碱基数目大于7个的read，合并双端read；
h)针对每个双端read，得到双端read覆盖的所有参考碱基位置{PEP1,PEP,…,PEPp}(即为集合R)，将该集合转为长度为q的向量Vecpe：
i)将矩阵Mat和向量Vecpe按照如下公式计算得到向量VecS：
j)按照下述方法，计算Vecdp数值变化，假设加入VecS前，Vecdp第k位数值是DPk,j―1，加入VecS计算之后，数值是DPk,j，则DPj,j和DPk,j―1的关系如下：
k)得到每个连锁位置k对应的DPk，因而每个LV突变Bj，突变支持数为ADj，深度为DPk，突变频率为ADj/DPk；
合并每个计算单元，得到该样本的所有2连锁或者3连锁的突变支持数、深度以及突变频率信息。
第二方面，本发明提供了一种黑名单存储装置，具体包括：
本发明涉及黑名单存储装置建立，该装置用于存储人群单核苷酸突变(single-nucleotide variant,SNV)变异信息和连锁基因突变(linkage variant,LV)信息，主要应用于定制化靶向捕获探针方法。这些信息可能来源于胚系突变和测序噪音等。在进行二代测序过程，由于实验过程造成的氧化损失、PCR错误和测序荧光信号读取误差容易引入低频噪音突变，在一些实施案例中，会需要到该装置对测序结果进行矫正。
1.提取20个正常人的cfDNA进行双端测序，获取测序数据，使用hg19参考基因组进行对比，获得比对信息，并标记重复片段；
2.针对每个样本数据，并行计算该样本每个位点的每类突变频率，具体方法如下：将探针目标区域延伸340bp(假设极端情况，3连锁SNV突变，一个SNV点位于探针目标区域，根据cfDNA片段峰值170bp，3个连锁SNV突变，最大SNV位点相距170乘以2，相距340bp)，使用本地计算装置按照目标延伸区域多线程并行计算每个位点的测序深度(depth，DP)以及4类碱基(ATCG)对应的数目(allele count，AD)，需要按照read(测序读段)比对质量(mappingquality，MQ)和碱基质量(base quantity，BQ)进行过滤，按照MQ≥20，BQ≥20设置，然后得出每个样本在该位置对应4类碱基的突变频率AF(AD/DP)；
3.针对每个样本，并行计算该样本所有可能的LV信息，具体算法如下：
A.过滤soft-clip(软切除)的read，按照MQ≥20；
B.通过read上MD标签，获得每个read上的错配碱基数目，如果read上错配碱基(排除unknown碱基，N碱基)≥7个，则过滤该read；
C.将双端read合并为一个片段，如果Read1和Read2的重叠部分碱基不同，则修改为N碱基，如果两个read没有重叠，则将空白区域标记为N碱基；
D.标记合并片段上的SNV突变碱基(alt，排除N碱基)，如果一个片段上≥2个alt，则认为该片段包含一个可能LV突变；
E.将所有存在LV突变的片段按照首个SNV的突变位置进行排序，迭代所有片段，获得以LV突变为索引，LV突变片段数目为列的LV突变矩阵；
F.合并所有LV突变，构成连锁突变区域，如果某一LV突变i，其染色体编号为chrom0，起始的SNV突变位置是start0,，终止的SNV突变位置是end0，另一LV突变j，其染色体编号是chrom1，起始SNV突变位置是start1，终止SNV突变位置是end1，突变j满足如下a)或者b)条件：
a)chrom0＝chrom1并且start1≥start0并且start1≤end0
b)chrom0＝chrom1并且start0≥start1并且start0≤end1
则将连锁突变i和连锁突变j合并到一个连锁突变区域；
G.以一个连锁突变区域为计算单元，针对该区域完成如下计算：
a)假设该连锁区域包含n个LV连锁突变，每个LV连锁突变，包含多个SNV单元，即为集合Ai，对应的片段支持数为ADi；
b)针对SNV单元数目大于等于3的LV突变，假设为该SNV单元数目为x，从中随机抽取2个构成2连锁LV，则包含个组合，从中随机抽取3个构成3连锁LV，则包含个组合，最终将1个x连锁的LV突变拆成个3连锁或者2连锁的LV突变；
c)将SNV单元数目等于2的LV突变和上述拆分所有3连锁或者2连锁的LV突变合并成一个包含m个2连锁或者3连锁的LV突变集合M，集合中每个LV突变是包含2到3个SNV突变单元(即为集合Bj)，则对应的片段支持数ADj按照如下公式进行计算：
d)上述每个LV突变Bj，假设SNV的坐标位置为{POS1,POS2,…,POSj}(即为集合Pk)，具有相同的突变位置LV合并，然后构成k个连锁突变位置集合K；
e)假设该连锁突变区域，起始位置为s1，终止位置为sq，则将集合K，转变为k×q的01矩阵Mat：
f)同时构建一个k的0向量VceDP，用以存储集合K中每个元素对应的深度信息DPk；
g)迭代该连锁区域的所有read，过滤掉标签为重复、补充比对、次比对、比对异常、MQ≤20和错配碱基数目大于7个的read，合并双端read；
h)针对每个双端read，得到双端read覆盖的所有参考碱基位置{PEP1,PEP,…,PEPp}(即为集合R)，将该集合转为长度为q的向量Vecpe：
i)将矩阵Mat和向量Vecpe按照如下公式计算得到向量VecS：
j)按照下述方法，计算Vecdp数值变化，假设加入VecS前，Vecdp第k位数值是DPk,j―1，加入VecS计算之后，数值是DPk,j，则DPk,j和DPk,j―1的关系如下：
k)这样得到每个连锁位置k对应的DPk，因而每个LV突变Bj，突变支持数是ADj，深度是DPk，突变频率是ADj/DPk；
H.合并每个计算单元，得到该样本的所有2连锁或者3连锁的突变支持数、深度以及突变频率信息；
4.收集每个样本的SNV及LV突变信息，构建SNV和LV突变数据库，SNV以SNV突变位置和突变类型为索引，以样本编号为列名，存储每个样本对应的SNV突变支持数、测序深度以及SNV突变频率，LV类似。
第三方面，本发明提供了一种电子设备，包括处理器、通信接口、存储器和通信总线，其中，处理器、通信接口和存储器通过通信总线完成相互间的通信；
存储器，用于存放计算机程序；
处理器，用于执行存储器上所存放的程序时，实现第一方面所述的连锁基因突变数据库构建方法。
第四方面，本发明提供了一种初诊监测位点检测装置，
本发明设计用于MRD监测的初诊样本监测位点检测装置建立。该装置主要利用组织切片或者cfDNA(tumor)，以及配对的血细胞样本(normal)，获得可能致病的LV突变信息，使用该信息进行MRD监测。初诊监测位点检测装置包括依次连接的比对信息单元、突变标注单元、插入片段长度计算单元、合并连锁突变单元、突变矩阵构建单元、连锁突变信息计算单元，以及筛选连锁突变单元；
比对信息单元，用于针对配对的tumor和normal样本，使用上述探针完成靶向捕获建库，并进行双端测序，获取测序数据，使用hg19参考基因组进行对比，获得比对信息，并标记重复片段；
突变标注单元，用于使用配对样本的比对信息进行SNV变异检测，使用fisher检验，标注胚系SNV突变和体系SNV突变；
插入片段长度计算单元，用于计算tumor样本的插入片段长度分布平均值μ和标准差σ，按照如下公式计算得到覆盖99.9％的插入片段最大长度L，L＝μ+3σ；
合并连锁突变单元，用于根据L将所有胚系和体系SNV突变合并成潜在的LV突变组，每组中从第一个体系SNV突变开始，从与之距离≤L的体系SNV集合依次抽出1个体系SNV，构成一对潜在2连锁SNV突变，迭代上述步骤最后得到一系列非重复的潜在2连锁SNV突变；从与之距离≤L的体系SNV集合依次抽出2个体系SNV，迭代上述步骤最后得到一系列非重复的潜在3连锁SNV突变；将潜在2连锁的SNV突变，从于与之距离≤L的胚系SNV集合抽出一个胚系SNV，构成潜在3连锁SNV突变；所有的潜在连锁突变，如果某一LV突变i，其染色体编号为chrom0，起始的SNV突变位置是start0,，终止的SNV突变位置是end0，另一LV突变j，其染色体编号是chrom1，起始SNV突变位置是start1，终止SNV突变位置是end1，chrom0＝chrom1并且start1≥start0并且end1≤end0，则将j纳入i突变区间内；
突变矩阵构建单元，用于针对所有合并的潜在LV突变区域，构建连锁突变矩阵；
连锁突变信息计算单元，用于按照上述区域并行计算tumor和normal对应LV的深度DP和突变支持数AD；
筛选连锁突变单元，用于根据第一方面方法得到的黑名单数据库，过滤掉噪音连锁突变。
进一步的，所述突变矩阵构建单元，具体用于：
针对所有合并的潜在LV突变区域，构建如下矩阵，假设该区域包含从编号为X1开始到编号为Xi的LV突变，假设该区域染色体编号为变量Y，从位置Z1开始，一直到位置Zj为止，通过如下碱基数字对应关系：
将每个潜在LV突变区域转为i×j的矩阵Mat：
该矩阵每一行代表一个潜在LV突变，Basenij代表该位置SNV碱基编号；
同时也构建一个包含i个元素的向量VecLV：
该向量的每个元素SNi代表潜在LV突变对应的SNV个数。
进一步的，所述连锁突变信息计算单元，具体用于：
按照上述区域并行计算tumor和normal对应LV的深度DP和突变支持数AD，主要方法如下：
1)构建长度为i的0元素向量Vec0，分别赋值给Vecdp和Vecad，Vecdp和Vecad分别用以存储对应LV的深度信息和突变支持数信息；
2)迭代比对到染色体编号为变量Y，起始位置Z1和Zj之间所有read，过滤掉标签为重复、补充比对、次比对、比对异常、MQ≤20和错配碱基数目大于7个的read，根据上述碱基数字转化公式以及read的参考基因组比对位置，将read序列转化为长度为j的向量Vecse，如果read对应的参考碱基是插入突变，则跳过该read上位置。
3)将单端read1的向量Vecse1和单端read2的向量Vecse2按照如下公式合并为双端read的向量Vecpe：
Vecse1＝[R1Basen1 … R1Basenj]
Vecse2＝[R2Basen1 … R2Basenj]
Vecpe＝[PRBasen1 … RRBasenj]
4)将矩阵Mat和Vecpe按照如下公式进行计算，得到Matpe:
5)根据Matpe，按照如下公式计算得到向量Vecdppe和向量Vecadpe：
6)按照下述方法，计算Vecdp和Vecad数值变化，假设加入Matpe前，Vecdp和Vecad第i位数值分别是DPi,，j―1和ADi,j―1，加入Matpe计算之后，数值分别是DPi,j和ADi,j，则DPi,j和DPi,j―1的关系如下：
ADi,j和ADi,j―1的关系如下：
7)迭代该区域对应的所有配对read，可得到每个潜在LV对应的深度和突变支持数矩阵，分别处理tumor和normal的样本比对数据，最后得到tumor和normal对应的所有LV深度、突变支持数和突变频率AF(AD/DP)。
进一步的，所述筛选连锁突变单元，具体用于：
通过对比tumor和normal的LV突变信息，使用fisher检验获得可能的致病LV突变，然后按照下述条件筛选监控致病突变：大于等于一定突变频率(当定制化panel时，组织样本可设置为2％，血浆样本时可设置为0.2％，当为全外显子探针时，组织样本设置为5％，血浆样本设置为1％，组织样本全基因测序时设置为10％)、测序深度(默认为100X，组织样本全基因组测序时设置为30X)和突变支持数(默认为2)，LV至少一个SNV突变位于探针目标区域前后延伸50bp内；
如果某N个SNV连锁突变，和黑名单数据库中SNV有交集，并且交集个数大于等于N-1，这N-1个SNV突变的黑名单人群频率大于等于5％并且突变频率大于等于1％，则将该连锁突变移除监控位点列表；
如果某个LV，存在黑名单数据库中，并且该LV人群频率大于等于5％，通过单边格拉布斯检验法判断该LV是否显著大于黑名单数据库中LV突变频率，如果不是，则将该LV移除监控位点列表。
第四方面，本发明提供了一种MRD监测装置，
本发明涉及MRD监测装置建立。该装置主要利用综合治疗之后的cfDNA，通过LV位点检测，来判断是非有MRD残留或者肿瘤复发可能性，主要通过随机抽样来判断噪音干扰。
MRD监测装置包括：依次连接的信息采集单元、数据构建单元和检验单元；
信息采集单元，用于首先提取患者综合治疗之后的血液中cfDNA，然后使用探针进行捕获建库及测序，在一些实施案例中，这里可是初诊建立基线所用到的基因组合(即Panel)，在另一些实施案例中，这里需要依据初诊结果进行Panel的定制；测序之后，通过比对标记重复获得比对信息；使用黑名单数据库建立的方法获得所有可能的连锁突变；
数据构建单元，用于提取初诊监控位点对应的噪音突变信息，采用蒙特卡洛抽样方法，假设为k个连锁突变，从剩余突变中筛选可能的噪音突变，每个监控位点突变筛选50个相应的噪音突变，噪音连锁突变的每个SNV突变类型需要和监控的连锁突变染色体编号一致，突变类型一致，并且每个SNV突变距离需要和监控连锁突变的SNV距离一致，最后得到所有监控位点对应的个噪音突变；
检验单元，用于统计监控位点对应的突变频率平均值μ0，完成10000次抽样，每次抽样依次从k个监测位点对应的噪音位点随机抽样，如从第i个监控位点的噪音突变{NSi,1,NSi,2…,NSi,50}随机抽一个，然后计算每次抽选的k个噪音突变的平均值μi，最后得到了{μ1,μ2…,μ10000}所有抽样结果，通过单边格拉布斯检验法判断μ0是否显著大于{μ1,μ2…,μ10000}，得到检验P值；当P<0.05时，判定MRD为阳性，定义MRD样本的监控位点对应的LV突变分子数目为x，这些监控位点总的测序深度为d，则样本的肿瘤占比为λ，λ等于x/d；当P≥0.05时，判定MRD结果阴性或低于检测下限。
本申请实施例提供的上述技术方案与现有技术相比具有如下优点：
本发明建立了一种基于连锁突变的MRD监测计算装置。相比于其他流式细胞方法和PET-CET方法，本发明方法样本获取简单，仅需要患者血液，操作方便，无放射性伤害；相比于PCR定量方法，本发明能一次检测多个位点，并且检测下限可达10―6级别，可有效地进行MRD术后评价和复发监控，根据患者个体情况开展精准指导；
本发明不依赖于UMI分子标签技术，当仅使用常规方法进行cfDNA建库时，本发明检测下限也可达到10―5级别，相对于双重测序duplex sequence，本发明节约实验建库成本，能有效提高原始片段的利用率；
当使用UMI分子标签技术时，本发明仅可使用该技术标记实验过程引物的重复序列，保留原始片段中的天然重复序列，进一步提高原始片段利用率，本发明无需测序过饱和来达到UMI分子校正所需数据量；对比表1和表2可看出，使用UMI标记重复技术之后，本发明的检测到ctDNA片段相对增多，即使在10―6级别本发明也能稳定检测肿瘤残余。
本发明技术适用领域广泛，可基于靶向panel测序或者全外显子、全基因组等测序技术获得的初诊连锁突变进行MRD监测；
利用本发明设计的探针以及对应的计算装置可检测肿瘤病人综合治疗之后外周血的肿瘤负荷，突变检测下限可达0.0001％。通过监测病人综合治疗之后外周血的肿瘤负荷，可评估治疗效果，监测病人复发可能性；
常规NGS检测方法在肿瘤占比较大时，如＞2％时，实际检测和理论检测比较符合；但是当肿瘤占比低于1％时，实际检测和理论检测相关性较差。通过模拟混样，从肿瘤突变占比0.25％依次降低，对实际检测和理论检测占比进行线性回归分析，R2相关性达到0.98以上，说明本发明在肿瘤占比较低的情况下，实际检测占比与理论值也有很高的符合度。
附图说明
此处的附图被并入说明书中并构成本说明书的一部分，示出了符合本发明的实施例，并与说明书一起用于解释本发明的原理。
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，对于本领域普通技术人员而言，在不付出创造性劳动性的前提下，还可根据这些附图获得其他的附图。
图1为本申请实施例装置示意图；
图2为本申请实施例预期突变和检测突变占比回归分析示意图。
具体实施方式
为使本申请实施例的目的、技术方案和优点更加清楚，下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本申请的一部分实施例，而不是全部的实施例。基于本申请中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本申请保护的范围。
连锁基因突变是指位于同一条染色体的基因突变。假设某个突变来自背景噪音的概率为千分之一，两个连锁突变都为背景噪音，则是两个独立事件同时发生，概率下降至百万分之一。连锁基因突变判断方法主要有三种：家系方法、群体遗传学方法和物理方法。家系法是使用目标样本和亲本的数据进行突变连锁判断的方法。群体遗传学方法是使用群体中大量无血缘关系的个体数据，根据连锁不平衡原理，构建模型，判断突变连锁的方法。物理方法就是简单地通过两个突变位于同一个片段来判断突变连锁。某些肿瘤，如淋巴瘤、黑色素瘤等，通常肿瘤组织中会存在大量的连锁突变，通过监测此类肿瘤基因组上连锁突变在治疗后的残留情况，理论上可达到超低检测下限的检测效果，本发明基于上述理论知识体系，开发了一套基于连锁突变的MRD标志物检测方法及其试剂盒，可达到远超目前市场上各类检测服务检测下限的效果。
如图1所示本发明提供了黑名单存储装置、初诊检测位点检测装置和MRD监测装置。通过本发明算法获得LV黑名单数据库和SNV黑名单数据库，然后对基于本发明算法检测的LV初诊位点进行过滤，最后将该位点集合用于MRD检测进行显著性计算和肿瘤负荷计算。
实施1基于靶向基因组合(Panel)的实施案例
在本实施案例中，靶向Panel可是某种疾病中常见突变基因的组合，优选容易发生体细胞高频突变的基因组合。
1.黑名单存储装置的建立
收集20例正常人的血浆样本，针对每个样本完成如下操作：
1)使用离心机进行血浆分离，使用商业方法提取cfDNA，使用常规方法或者分子标签(UMI)完成建库，使用定制化的panel获得目的基因，然后二次PCR扩增；
2)使用illumina或者其他平台完成双端测序，读长为150bp，数据量至少1000X；
3)上述样本数据下机之后，完成bcl2fastq数据拆分，获得每个样本的fastq数据；
4)利用fastp对原始测序数据进行过滤以及质控处理，获得过滤后的fastq数据；
5)将过滤后的数据利用fgbio将fastq格式文件转化为ubam格式文件以提取umi标签序列；
6)利用BWA mem将测序reads比对到hg19参考基因组上，使用picardMarkDuplicates根据UMI和比对位置标记重复序列；针对常规杂交捕获方法，可直接使用fastp过滤后的cfDNA的fastq文件进行比对，使用picard MarkDuplicates根据比对位置标记重复序列；
7)使用samtools view接口过滤比对质量小于20的read、标记为重复、不合理比对以及次级比对read，获得剩余read构成的比对文件；
8)针对每个样本数据，利用步骤7)的比对文件，并行计算该样本每个位点的每类突变频率，具体方法如下：将探针目标区域延伸340bp(假设极端情况，3连锁SNV突变，一个SNV点位于探针目标区域，根据cfDNA峰值170bp，3个连锁SNV突变，最大SNV位点相距170乘以2，相距340bp)，使用本地计算装置按照目标延伸区域多线程并行计算每个位点的测序深度(depth，DP)以及4类碱基(ATCG)对应的数目(allelecount，AD)，需要按照read(测序读段)比对质量(mapping quality，MQ)和碱基质量(basequantity，BQ)进行过滤，按照MQ≥20，BQ≥20设置，然后得出每个样本在该位置对应4类碱基的突变频率AF(AD/DP)；
9)针对每个样本，利用步骤6)的比对文件，并行计算该样本所有可能的LV信息，使用本地构建的LV检测算法程序，具体算法实现如下：
A.过滤soft-clip(软切除)的read，按照MQ≥20筛选read；
B.通过read上MD标签，获得每个read上的错配碱基数目，如果read上错配碱基(排除unknown碱基，N碱基)≥7个，则过滤该read；
C.将双端read合并为一个片段，如果Read1和Read2的重叠部分碱基不同，则修改为N碱基，如果两个read没有重叠，则将空白区域标记为N碱基；
D.标记合并片段上的SNV突变碱基(alt，排除N碱基)，如果一个片段上≥2个alt，则认为该片段包含一个可能LV突变；
E.将所有存在LV突变的片段按照首个SNV的突变位置进行排序，迭代所有片段，获得以LV突变为索引，LV突变片段数目为列的LV突变矩阵；
F.合并所有LV突变，构成连锁突变区域，如果某一LV突变i，其染色体编号为chrom0，起始的SNV突变位置是start0,，终止的SNV突变位置是end0，另一LV突变j，其染色体编号是chrom1，起始SNV突变位置是start1，终止SNV突变位置是end1，突变j满足如下a)或者b)条件：
a)chrom0＝chrom1并且start1≥start0并且start1≤end0
b)chrom0＝chrom1并且start0≥start1并且start0≤end1
则将连锁突变i和连锁突变j合并到一个连锁突变区域；
G.以一个连锁突变区域为计算单元，针对该区域完成如下计算：
a)假设该连锁区域包含n个LV连锁突变，每个LV连锁突变，包含多个SNV单元，即为集合Ai，对应的片段支持数为ADi；
b)针对SNV单元数目大于等于3的LV突变，假设为该SNV单元数目为x，从中随机抽取2个构成2连锁LV，则包含个组合，从中随机抽取3个构成3连锁LV，则包含个组合，最终将X连锁的LV突变拆成个3连锁或者2连锁的LV突变；
c)将SNV单元数目等于2的LV突变和上述拆分所有3连锁或者2连锁的LV突变合并成一个包含m个2连锁或者3连锁的LV突变集合M，集合中每个LV突变包含2到3个SNV突变单元，即为集合Bj，则对应的片段支持数ADj按照如下公式进行计算：
d)上述每个LV集合Bj，假设SNV的坐标位置为{POS1,POS2,…,POSj}，即为集合Pk，具有相同LV的突变位置合并，然后构成k个连锁突变位置集合K；
e)假设该连锁突变区域，起始位置为s1，终止位置为sq，则将集合K，转变为k×q的01矩阵Mat：
f)同时构建一个k的0向量VecDP，用以存储集合K中每个元素对应的深度信息DPk；
g)迭代该连锁区域的所有read，过滤掉标签为重复、补充比对、次比对、比对异常、MQ≤20和错配碱基数目大于7个的read，合并双端read；
h)针对每个双端read，得到双端read覆盖的所有参考碱基位置{PEP1,PEP,…,PEPp}，即为集合R，将该集合转为长度为q的向量Vecpe：
i)将矩阵Mat和向量Vecpe按照如下公式计算得到向量VecS：
j)按照下述方法，计算Vecdp数值变化，假设加入VecS前，Vecdp第k位数值是DPk,j―1，加入VecS计算之后，数值是DPk,j，则DPk,j和DPk,j―1的关系如下：
k)这样得到每个连锁位置k对应的DPk，因而每个LV突变Bj，突变支持数是ADj，深度是DPk，突变频率是ADj/DPk；
H.合并每个计算单元，得到该样本的所有2连锁或者3连锁的突变支持数、深度以及突变频率信息；
10)收集每个样本的SNV及LV突变信息，构建SNV和LV突变数据库，SNV以SNV突变位置和突变类型为索引，以样本编号为列名，存储每个样本对应的SNV突变支持数、深度以及SNV突变频率，LV类似；
2.获取初诊数据
患者初诊时，送检组织活检或者切片样本，进行本探针进行靶向捕获测序，获取初诊测序数据。此步骤实验过程如下：
1)使用商业方法提取DNA，使用超声打断，完成建库，使用定制化panel进行目标片段获取，然后进行二次PCR扩增；
2)将上述制备好的文库，进行测序反应，测序仪可是illumina任意品类测序仪(如:Nextseq550)，测序读长为双端150bp，数据量要求为500X以上；
3)上述数据测序下机之后，完成bcl2fastq数据拆分，获得该样本的fastq数据；
4)利用fastp对原始测序数据进行过滤以及质控处理，获得过滤后的fastq数据；
5)利用BWA mem将测序reads比对到hg19参考基因组上，获得比对数据，使用picard标记重复的插入片段，获得标记重复的比对数据。
送检初诊血液样本时，使用“1.黑名单存储装置的建立”的步骤1)-6)获得标记重复的比对文件；
3.获取MRD数据
患者经过治疗后，送检外周血样本，对外周血有核细胞及血浆ctDNA进行与初诊所用的Panel相同的Panel进行检测，获取MRD测序数据及配对数据。其具体实施方法如下：
1)使用离心机分离有核细胞和血浆，按照“1.黑名单存储装置的建立”的建立步骤1)获得cfDNA文库；
2)将上述制备好cfDNA的文库，进行测序反应，测序仪可是illumina任意品类测序仪(如:Nextseq550)，测序读长为双端150bp，深度要求为2000X；
3)将有核细胞按照上述“2.获得初诊数据”步骤1)操作，和制备好的cfDNA文库一起上机，有核细胞数据量要求为200X以上；
4)上述样本数据下机之后，完成bcl2fastq数据拆分，获得cfDNA样本和有核细胞的fastq数据；cfDNA测序数据按照“1.黑名单存储装置的建立”的步骤4)-6)进行操作获得标记重复的比对数据，有核细胞的测序数据按照“2.获得初诊数据”的4)-5)步骤操作，获得标记重复的比对数据；
4.获取MRD监控位点
基于“2.获得初诊数据”中初诊样本标记重复的比对和“3.获得MRD数据”中有核细胞标记重复的比对数据，使用本地算法获得监控位点，具体实现如下：
1)延伸探针目标区域340bp，获得SNV检测区域，使用VarDict根据有核细胞和cfDNA的标记重复比对数据，获得胚系SNV突变和体系SNV突变；
2)使用picard的CollectInsertSizeMetrics模块计算tumor样本的插入片段长度分布平均值μ和标准差σ，按照如下公式计算得到覆盖99.9％的插入片段最大长度L，L＝μ+3σ；
3)根据L将所有胚系和体系SNV突变合并成潜在的LV突变组，每组中从第一个体系SNV突变开始，从与之距离≤L的体系SNV集合依次抽出1个体系SNV，构成一对潜在2连锁SNV突变，迭代上述步骤最后得到一系列非重复的2连锁SNV突变；从与之距离≤L的体系SNV集合依次抽出2个体系SNV，迭代上述步骤最后得到一系列非重复3连锁SNV突变；
4)这些潜在2连锁的SNV突变，从于与之距离≤L的胚系SNV集合抽出一个胚系SNV，构成潜在3连锁SNV突变；
5)所有的潜在连锁突变，如果某一LV突变i，其染色体编号为chrom0，起始的SNV突变位置是start0,，终止的SNV突变位置是end0，另一LV突变j，其染色体编号是chrom1，起始SNV突变位置是start1，终止SNV突变位置是end1，chrom0＝chrom1并且start1≥start0并且end1≤end0，则将j纳入i突变区间内；
6)针对所有合并的潜在LV突变区域，构建如下矩阵，假设该区域包含从编号为X1开始到编号为Xi的LV突变，假设该区域染色体编号为变量Y，从位置Z1开始，一直到位置Zj为止，通过如下碱基数字对应关系：
将每个潜在LV突变区域转为i×j的矩阵Mat：
该矩阵每一行代表一个潜在LV突变，Basenij代表该位置SNV碱基编号；
同时也构建一个包含i个元素的向量VecLV：
该向量的每个元素SNi代表潜在LV突变对应的SNV个数。
7)按照上述区域并行计算tumor和normal对应LV的深度DP和突变支持数AD，主要方法如下：
A.构建长度为i的0元素向量Vec0，分别赋值给Vecdp和Vecad，Vecdp和Vecad分别用以存储对应LV的深度信息和突变支持数信息；
B.迭代比对到染色体编号为变量Y，起始位置Z1和Zj之间所有read，过滤掉标签为重复、补充比对、次比对、比对异常、MQ≤20和错配碱基数目大于7个的read，根据上述碱基数字转化公式以及read的参考基因组比对位置，将read序列转化为长度为j的向量Vecse，如果read对应的参考碱基是插入突变，则跳过该read上位置。
C.将单端read1的向量Vecse1和单端read2的向量Vecse2按照如下公式合并为双端read的向量Vecpe：
Vecse1＝[R1Basen1 … R1Basenj]
Vecse2＝[R2Basen1 … R2Basenj]
Vecpe＝[PRBasen1 … RRBasenj]
D.将矩阵Mat和Vecpe按照如下公式进行计算，得到Matpe:
E.根据Matpe，按照如下公式计算得到向量Vecdppe和向量Vecadpe：
F.按照下述方法，计算Vecdp和Vecad数值变化，假设加入Matpe前，Vecdp和Vecad第i位数值分别是DPi,j―1和ADi,j―1，加入Matpe计算之后，数值分别是DPi,j和ADi,j，则DPi,j和DPi,j―1的关系如下：
ADi,j和ADi,j―1的关系如下：
G.迭代该区域对应的所有配对read，可得到每个潜在LV对应的深度和突变支持数矩阵，分别处理tumor和normal的样本比对数据，最后得到tumor和normal对应的所有LV深度、突变支持数和突变频率AF(AD/DP)；
8)通过对比tumor和normal的LV突变信息，使用fisher检验获得可能的致病LV突变，然后按照下述条件筛选用于监控的致病突变：突变频率大于等于2％(当初诊样本为血液时设置为0.1％)、测序深度≥100X和突变支持数大于等于2，LV至少一个SNV突变位于探针目标区域前后延伸50bp内；
9)如果某N个SNV连锁突变，和黑名单数据库中SNV有交集，并且交集个数大于等于N-1，这N-1个SNV突变的黑名单人群频率大于等于5％并且突变频率大于等于1％，则将该连锁突变移除监控位点列表；
10)如果某个LV，存在黑名单数据库中，并且该LV人群频率大于等于5％，通过单边格拉布斯检验法判断该LV是否显著大于黑名单数据库中LV突变频率，如果不是，则将该LV移除监控位点列表；
5.MRD计算
基于上述获得的LV监控位点和“3.获得MRD数据”中cfDNA标记重复的比对数据，进行如下操作：
1)针对cfDNA标记重复的比对文件，按照“1.黑名单存储装置的建立”中9)LV检测算法获得所有LV真实突变和噪音突变；
2)使用本地计算装置按照如下操作计算MRD样本的显著性：
A.提取初诊监控位点对应的突变信息，假设为k个，从剩余突变中筛选可能的噪音突变，每个监控位点突变筛选50个相应的噪音突变，噪音连锁突变的每个SNV单位突变类型需要和监控的连锁突变的SNV单位突变的染色体编号致，突变类型一致，并且每个SNV突变距离需要和监控连锁突变的SNV距离一致，最后得到所有监控位点对应的个噪音突变；
B.统计监控位点对应的突变频率平均值μ0，完成10000次抽样，每次抽样依次从k个监测位点对应的噪音位点随机抽样，如从第i个监控位点的噪音突变{NSi,1,NSi,2…,NSi,50}随机抽一个，然后计算每次抽选的k个噪音突变的平均值μi，最后得到了{μ1,μ2…,μ10000}所有抽样结果，通过单边格拉布斯检验法判断μ0是否显著大于{μ1,μ2…,μ10000}，得到检验P值；
C.当MRD样本显著(P＜0.05)时，定义MRD样本的监控位点对应的LV突变分子数目为x，这些监控位点总的测序深度为d，则样本的肿瘤占比为λ，λ等于x/d；
实施2基于初诊为全外显子或者全基因组测序的MRD监测
本实施案例针对患者初诊选择全外显子测序或者全基因组测序时，通过个性化panel设计进行连锁突变MRD监测。
1.获取初诊数据
患者初诊时，送检组织活检或者切片样本，进行全外显子或者全基因组测序，获取初诊测序数据。此步骤实验过程如下：
1)使用商业方法提取DNA，使用超声打断，完成建库，使用全外显子探针进行目标片段获取(全基因测序无需这一步)，然后进行二次PCR扩增；
2)将上述制备好的文库，进行测序反应，测序仪可是illumina任意品类测序仪(如:Nextseq550)，测序读长为双端150bp，数据量要求为200X以上(全基因组测序至少50X)；
3)上述数据测序下机之后，完成bcl2fastq数据拆分，获得该样本的fastq数据；
4)利用fastp对原始测序数据进行过滤以及质控处理，获得过滤后的fastq数据；
5)利用BWA mem将测序reads比对到hg19参考基因组上，获得比对数据，使用picard标记重复的插入片段，获得标记重复的比对数据。
送检初诊血液样本时，进行全外显子测序，获得初诊样本测序数据，步骤如下：
A.使用离心机进行血浆分离，使用商业方法提取cfDNA，使用常规方法或者分子标签(UMI)完成建库，使用全外显子探针进行靶向捕获，然后二次PCR扩增；
B.使用illumina或者其他平台完成双端测序，读长为150bp，数据量至少200X；
C.上述样本数据下机之后，完成bcl2fastq数据拆分，获得每个样本的fastq数据；
D.利用fastp对原始测序数据进行过滤以及质控处理，获得过滤后的fastq数据；
E.将过滤后的数据利用fgbio将fastq格式文件转化为ubam格式文件以提取umi标签序列；
F.利用BWA mem将测序reads比对到hg19参考基因组上，使用picardMarkDuplicates根据UMI和比对位置标记重复序列；针对常规杂交捕获方法，可直接使用fastp过滤后的cfDNA的fastq文件进行比对，使用picard MarkDuplicates根据比对位置标记重复序列；
2.MRD监测位点获得
使用全外显子测序或者全基因组，利用初诊样本和治疗之后第一次MRD的血液样本提取的有核细胞获得连锁突变MRD监控位点：
1)使用离心机进行血浆分离，获得有核细胞，按照“1.获得初诊数据”组织样本处理步骤获得有核细胞标记重复的比对文件；
2)延伸探针目标区域340bp，获得SNV检测区域，使用VarDict根据有核细胞和cfDNA的标记重复比对数据，获得胚系SNV突变和体系SNV突变；
3)使用picard的CollectInsertSizeMetrics模块计算tumor样本的插入片段长度分布平均值μ和标准差σ，按照如下公式计算得到覆盖99.9％的插入片段最大长度L，L＝μ+3σ；
4)根据L将所有胚系和体系SNV突变合并成潜在的LV突变组，每组中从第一个体系SNV突变开始，从与之距离≤L的体系SNV集合依次抽出1个体系SNV，构成一对潜在2连锁SNV突变，迭代上述步骤最后得到一系列非重复的2连锁SNV突变；从与之距离≤L的体系SNV集合依次抽出2个体系SNV，迭代上述步骤最后得到一系列非重复3连锁SNV突变；
5)这些潜在2连锁的SNV突变，从于与之距离≤L的胚系SNV集合抽出一个胚系SNV，构成潜在3连锁SNV突变；
6)所有的潜在连锁突变，如果某一LV突变i，其染色体编号为chrom0，起始的SNV突变位置是start0,，终止的SNV突变位置是end0，另一LV突变j，其染色体编号是chrom1，起始SNV突变位置是start1，终止SNV突变位置是end1，chrom0＝chrom1并且start1≥start0并且end1≤end0，则将j纳入i突变区间内；
7)针对所有合并的潜在LV突变区域，构建如下矩阵，假设该区域包含从编号为X1开始到编号为Xi的LV突变，假设该区域染色体编号为变量Y，从位置Z1开始，一直到位置Zj为止，通过如下碱基数字对应关系：
将每个潜在LV突变区域转为i×j的矩阵Mat：
该矩阵每一行代表一个潜在LV突变，Basenij代表该位置SNV碱基编号；
同时也构建一个包含i个元素的向量VecLV：
该向量的每个元素SNi代表潜在LV突变对应的SNV个数。
8)按照上述区域并行计算tumor和normal对应LV的深度DP和突变支持数AD，主要方法如下：
A.构建长度为i的0元素向量Vec0，分别赋值给Vecdp和Vecad，Vecdp和Vecad分别用以存储对应LV的深度信息和突变支持数信息；
B.迭代比对到染色体编号为变量Y，起始位置Z1和Zj之间所有read，过滤掉标签为重复、补充比对、次比对、比对异常、MQ≤20和错配碱基数目大于7个的read，根据上述碱基数字转化公式以及read的参考基因组比对位置，将read序列转化为长度为j的向量Vecse，如果read对应的参考碱基是插入突变，则跳过该read上位置。
C.将单端read1的向量Vecse1和单端read2的向量Vecse2按照如下公式合并为双端read的向量Vecpe：
Vecse1＝[R1Basen1 … R1Basenj]
Vecse2＝[R2Basen1 … R2Basenj]
Vecpe＝[PRBasen1 … PRBasenj]
D.将矩阵Mat和Vecpe按照如下公式进行计算，得到Matpe:
E.根据Matpe，按照如下公式计算得到向量Vecdppe和向量Vecadpe：
F.按照下述方法，计算Vecdp和Vecad数值变化，假设加入Matpe前，Vecdp和Vecad第i位数值分别是DPi,j―1和ADi,j―1，加入Matpe计算之后，数值分别是DPi,j和ADi,j，则DPi,j和DPi,j―1的关系如下：
ADi,j和ADi,j―1的关系如下：
G.迭代该区域对应的所有配对read，可得到每个潜在LV对应的深度和突变支持数矩阵，分别处理tumor和normal的样本比对数据，最后得到tumor和normal对应的所有LV深度、突变支持数和突变频率AF(AD/DP)；
9)通过对比tumor和normal的LV突变信息，使用fisher检验获得可能的致病LV突变，然后按照下述条件筛选监控致病突变：突变频率大于等于一定阈值(当为全外显子探针时，组织样本设置为5％，血浆样本设置为1％，组织样本全基因测序时设置为10％)、测序深度(默认为100X，组织样本全基因组测序时设置为30X)、和突变支持数大于等于2，LV至少一个SNV突变位于探针目标区域前后延伸50bp内；
10)按照上述条件筛选最后得到MRD样本连锁突变监控位点。
3.MRD监测Panel定制
第一次MRD样本送检时，需要将外显子或者全基因初诊所获得的所有LV监控位点连锁区域合并，按照如下条件合并：
如果某一LV突变i，其染色体编号为chrom0，起始的SNV突变位置是start0,，终止的SNV突变位置是end0，另一LV突变j，其染色体编号是chrom1，起始SNV突变位置是start1，终止SNV突变位置是end1，突变j满足如下a)或者b)条件：
a)chrom0＝chrom1并且start1≥start0并且start1≤end0
b)chrom0＝chrom1并且start0≥start1并且start0≤end1
则将连锁突变i和连锁突变j合并到一个连锁突变区域；
上述方法合并之后获得连锁突变区域，然后进行靶向Panel设计，设计优选委托IDT公司进行设计服务，后续MRD检测无需该步骤。
4.MRD计算
患者经过治疗后，送检外周血样本，利用外周血样本进行MRD计算，具体方法如下：
1)按照按照“1.获得初诊数据”血液样本处理步骤获得cfDNA样本标记重复的比对数据，所不同的是，目标区域捕获使用“2.MRD监测Panel定制”中定制的Panel，数据至少10000X；
2)针对cfDNA标记重复的比对文件，过滤soft-clip(软切除)的read，按照MQ≥20；
3)通过read上MD标签，获得每个read上的错配碱基数目，如果read上错配碱基(排除unknown碱基，N碱基)≥7个，则过滤该read；
4)将双端read合并为一个片段，如果Read1和Read2的重叠部分碱基不同，则修改为N碱基，如果两个read没有重叠，则将空白区域标记为N碱基；
5)标记合并片段上的SNV突变碱基(alt，排除N碱基)，如果一个片段上≥2个alt，则认为该片段包含一个可能LV突变；
6)将所有存在LV突变的片段按照首个SNV的突变位置进行排序，迭代所有片段，获得以LV突变为索引，LV突变片段数目为列的LV突变矩阵；
7)按照“3.MRD监测Panel定制”的方法合并所有连锁突变，构建连锁突变区域；
8)以一个连锁突变区域为计算单元，针对该区域完成如下计算：
A.假设该连锁区域包含n个LV连锁突变，每个LV连锁突变，包含多个SNV单元，即为集合Ai，对应的片段支持数为ADi；
B.针对SNV单元数目大于等于3的LV突变，假设为该SNV单元数目为x，从中随机抽取2个构成2连锁LV，则包含个组合，从中随机抽取3个构成3连锁LV，则包含个组合，最终将X连锁的LV突变拆成个3连锁或者2连锁的LV突变；
C.将SNV单元数目等于2的LV突变和上述拆分所有3连锁或者2连锁的LV突变合并成一个包含m个2连锁或者3连锁的LV突变集合M，集合中每个LV突变包含2到3个SNV突变单元，即为集合Bj，则对应的片段支持数ADj按照如下公式进行计算：
D.上述每个LV突变Bj，假设SNV的坐标位置为{POS1,POS2,…,POSj}，即为集合Pk，具有相同LV的突变位置合并，然后构成k个连锁突变位置集合K；
E.假设该连锁突变区域，起始位置为s1，终止位置为sq，则将集合K，转变为k×q的01矩阵Mat：
F.同时构建一个k的0向量VceDP，用以存储集合K中每个元素对应的深度信息DPk；
G.迭代该连锁区域的所有read，过滤掉标签为重复、补充比对、次比对、比对异常、MQ≤20和错配碱基数目大于7个的read，合并双端read；
H.针对每个双端read，得到双端read覆盖的所有参考碱基位置{PEP1,PEP,…,PEPp}，即为集合R，将该集合转为长度为q的向量Vecpe：
I.将矩阵Mat和向量Vecpe按照如下公式计算得到向量VecS：
J.按照下述方法，计算Vecdp数值变化，假设加入VecS前，Vecdp第k位数值是DPk,j―1，加入VecS计算之后，数值是DPk,j，则DPj,j和DPk,j―1的关系如下：
K.这样得到每个连锁位置k对应的DPk，因而每个LV突变Bj，突变支持数是ADj，深度是DPk，突变频率是ADj/DPk；
3)合并每个计算单元，得到该样本的所有2连锁或者3连锁的突变支持数、深度以及突变频率信息；；
4)使用本地计算装置按照如下操作计算MRD样本的显著性：
A.提取初诊监控位点对应的突变信息，假设为k个，从剩余突变中筛选可能的噪音突变，每个监控位点突变筛选50个相应的噪音突变，噪音连锁突变的每个SNV突变类型需要和监控的连锁突变的染色体编号一致，突变类型一致，并且每个SNV突变距离需要和监控连锁突变的SNV距离一致，最后得到所有监控位点对应的个噪音突变；
B.统计监控位点对应的突变频率平均值μ0，完成10000次抽样，每次抽样依次从k个监测位点对应的噪音位点随机抽样，如从第i个监控位点的噪音突变{NSi,1,NSi,2…,NSi,50}随机抽一个，然后计算每次抽选的k个噪音突变的平均值μi，最后得到了{μ1,μ2…,μ10000}所有抽样结果，通过单边格拉布斯检验法判断μ0是否显著大于{μ1,μ2…,μ10000}，得到检验P值；
C.当MRD样本显著(P＜0.05)时，定义MRD样本的监控位点对应的LV突变分子数目为x，这些监控位点总的测序深度为d，则样本的肿瘤占比为λ，λ等于x/d。
本发明的技术效果如下：
利用本发明设计的探针以及对应的计算装置可检测肿瘤病人综合治疗之后外周血的肿瘤负荷，突变检测下限可达0.0001％。通过监测病人综合治疗之后外周血的肿瘤负荷，可评估治疗效果，监测病人复发可能性。
1.使用3例经临床医院科研合作获取的的gDNA样本，按照一定比例混合正常人的超声打断的gDNA样本，模拟不同肿瘤负荷浓度的样本，按照本发明方法进行检测，得到如下表1所述结果：
表1 混样模拟测试统计
如表1中25个混样结果展示，混样比例从0.25％到低至0.0001％时，本发明装置仍然能够显著检出肿瘤残余。
2.常规NGS检测方法在肿瘤占比较大时，如＞2％时，实际检测和理论检测比较符合；但是当肿瘤占比低于1％时，实际检测和理论检测相关性较差。通过模拟混样，从肿瘤突变占比0.25％依次降低，对实际检测和理论检测占比进行线性回归分析，R2相关性达到0.98以上，说明本发明在肿瘤占比较低的情况下，实际检测占比与理论值也有很高的符合度。
预期突变和检测突变占比回归分析结果如图2所示，其中X轴代表预期突变占比的数值，Y轴代表检测突变占比，虚线是线性回归趋势线。图中预期突变占比和检测突变占比相关性R2相关性达到0.98以上。
本发明至少具备以下有益效果：
本发明建立了一种基于连锁突变的MRD监测计算装置。相比于其他流式细胞方法和PET-CET方法，本发明方法样本获取简单，仅需要患者血液，操作方便，无放射性伤害；相比于PCR定量方法，本发明能一次检测多个位点，并且检测下限可达10―6级别，可有效地进行MRD术后评价和复发监控，根据患者个体情况开展精准指导；
本发明不依赖于UMI分子标签技术，当仅使用常规方法进行cfDNA建库时，本发明检测下限也可达到10―5级别，相对于duplex sequence，本发明节约实验建库成本，能有效提高原始提取序列的利用率；
表2cfDNA常规方法建库混样测试统计
当使用UMI分子标签技术时，本发明仅可使用该技术标记实验过程引物的重复序列，保留原始片段中的天然重复序列，进一步提高原始片段利用率，本发明无需测序过饱和来达到UMI分子校正所需数据量；对比表1和表2可看出，使用UMI标记重复技术之后，本发明的检测到ctDNA片段相对增多，即使在10―6级别本发明也能稳定检测肿瘤残余。
本发明技术适用领域广泛，可基于靶向panel测序或者全外显子、全基因组等测序技术获得的初诊突变进行MRD监测；
本发明提出了一种使用连锁突变检测MRD技术方案和计算装置，包含连锁突变黑名单存储数据库构建装置、初诊样本LV致病突变检测算法和MRD样本监测算法；
本发明提出了一种高效、可靠的连锁突变信息检测算法，该算法通过已知胚系和体系SNV突变，构建LV潜在致病突变连锁区域，使用多进程针对每个连锁区域进行计算，将碱基信号转为质数编号，通过矩阵和向量运算，计算连锁突变信息；
本发明提出了一种针对固定panel而不是全外显子或者全基因组检测的连锁突变二次降噪系统，该系统包含一种高效检测所有LV突变以及噪音突变的算法，该系统首先通过此算法构建黑名单存储装置，利用此装置针对初诊样本中可能的非致病突变以及系统噪音进行过滤，通过单边格拉布斯检验法过滤初诊监测位点噪音；其次，通过此算法检测MRD样本连锁真实和噪音突变，通过随机抽样和单边格拉布斯检验法，判断MRD样本中噪音突变对MRD监测突变的影响。
需要说明的是，在本文中，诸如“第一”和“第二”等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来，而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。而且，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者设备所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括所述要素的过程、方法、物品或者设备中还存在另外的相同要素。
以上所述仅是本发明的具体实施方式，使本领域技术人员能够理解或实现本发明。对这些实施例的多种修改对本领域的技术人员来说将是显而易见的，本文中所定义的一般原理可在不脱离本发明的精神或范围的情况下，在其它实施例中实现。因此，本发明将不会被限制于本文所示的这些实施例，而是要符合与本文所申请的原理和新颖特点相一致的最宽的范围。
Claims (3)
1.一种连锁基因突变数据库构建方法，其特征在于，包括：
提取若干个正常人的cfDNA进行双端测序，获取测序数据作为样本数据，使用hg19参考基因组进行对比，获得比对信息，并标记重复片段；
针对每个样本数据，并行计算该样本每个位点的每类碱基突变频率；
针对每个样本数据，并行计算该样本所有可能的LV信息；
收集每个样本数据的SNV及LV信息，构建SNV和LV数据库；
其中SNV数据库以SNV突变位置和突变类型为索引，以样本编号为列名，存储每个样本对应的SNV突变支持数、测序深度以及SNV突变频率；LV数据库以LV突变位置和突变类型为索引，以样本编号为列名，存储每个样本对应的LV突变支持数、测序深度以及LV突变频率；
所述SNV表示单核苷酸突变；所述LV表示连锁基因突变；
所述针对每个样本数据，并行计算该样本每个位点的每类碱基突变频率，具体包括：
将测序使用的捕获探针针对的目标区域延伸340bp，使用本地计算装置按照目标延伸区域多线程并行计算每个位点的测序深度DP以及ATCG 4类碱基对应的数目AD，按照测序读段比对质量MQ和碱基质量BQ进行过滤，按照MQ≥20，BQ≥20设置，然后得出每个样本在该位置对应4类碱基的突变频率AF＝AD/DP；
所述针对每个样本数据，并行计算该样本所有可能的LV信息，具体包括：
A.过滤软切除的测序读段，按照MQ≥20筛选；
B.通过测序读段上MD标签，获得每个测序读段上的错配碱基数目，如果测序读段上错配碱基数量≥7个，则过滤该测序读段；
C.将双端测序读段合并为一个片段，如果测序读段1和测序读段2的重叠部分碱基不同，则修改为N碱基，如果两个测序读段没有重叠，则将空白区域标记为N碱基；
D.标记合并片段上的SNV突变碱基，并排除N碱基，如果一个片段上≥2个alt，则认为该片段包含一个可能LV突变；
E.将所有存在LV突变的片段按照首个SNV的突变位置进行排序，迭代所有片段，获得以LV突变为索引，LV突变片段数目为列的LV突变矩阵；
F.合并所有LV突变，构成连锁突变区域，某一LV突变i，其染色体编号为chrom0，起始的SNV突变位置是start0，终止的SNV突变位置是end0，另一LV突变j，其染色体编号是chrom1，起始SNV突变位置是start1，终止SNV突变位置是end1，如果突变j满足如下a)或者b)条件：
a)chrom0＝chrom1并且start1≥start0并且start1≤end0；
b)chrom0＝chrom1并且start0≥start1并且start0≤end1;
则将连锁突变i和连锁突变j合并到一个连锁突变区域；
G.以一个连锁突变区域为计算单元，针对该区域完成如下计算：
a)假设该连锁突变区域包含n个LV连锁突变，每个LV连锁突变，包含多个SNV单元，即为集合Ai，对应的片段支持数为ADi；
b)针对SNV单元数目大于等于3的LV突变，假设为该SNV单元数目为x，从中随机抽取2个SNV单元构成2连锁LV，则包含个组合，从中随机抽取3个构成3连锁LV，则包含个组合，最终将1个x连锁的LV突变拆分成 个3连锁或者2连锁的LV突变；
c)将SNV单元数目等于2的LV突变和上述拆分所有3连锁或者2连锁的LV突变合并成一个包含m个2连锁或者3连锁的LV突变集合M，集合Bj是集合M中每个LV突变包含的2到3个SNV突变单元组成的SNV突变单元的集合，则对应的片段支持数ADj按照如下公式进行计算：
d)假设上述每个LV突变j的每个SNV突变单元的坐标位置为{POS1,POS2,…,POSj}，具有相同的突变位置合并，构成k个连锁突变位置集合Pk；
e)假设该连锁突变区域，起始位置为s1，终止位置为sq，则将集合Pk转变为k×q的01矩阵Mat：
f)同时构建一个k的0向量VecDP，用以存储集合K中每个元素对应的深度信息DPk；
g)迭代该连锁区域的所有测序读段，过滤掉标签为重复、补充比对、次比对、比对异常、MQ≤20和错配碱基数目大于7个的测序读段，合并双端测序读段；
h)针对每个双端测序读段，得到双端测序读段覆盖的所有参考碱基位置{PEP1,PEP,…,PEPp}设为集合R，将集合R转为长度为q的向量Vecpe：
i)将矩阵Mat和向量Vecpe按照如下公式计算得到向量VecS：
j)按照下述方法，计算Vecdp数值变化，假设加入VecS前，Vecdp第k位数值是DPk,j―1，加入VecS计算之后，数值是DPk,j，则DPk,j和DPk,j―1的关系如下：
k)得到每个连锁位置k对应的深度信息DPk，每个LV突变Bj的突变支持数为ADj，可计算得到突变频率为ADj/DPk；
H.合并每个计算单元，得到针对连锁突变数据库中每个正常人样本的所有2连锁或者3连锁的突变支持数、深度以及突变频率信息。
2.一种初诊监测位点检测装置，其特征在于，包括依次连接的比对信息单元、突变标注单元、插入片段长度计算单元、合并连锁突变单元、突变矩阵构建单元、连锁突变信息计算单元，以及筛选连锁突变单元；
比对信息单元，用于针对肿瘤和正常配对样本，使用探针完成靶向捕获建库，并进行双端测序，获取测序数据，使用hg19参考基因组进行对比，获得比对信息，并标记重复片段；
突变标注单元，用于使用配对样本比对信息进行SNV变异检测，使用fisher检验，标注胚系SNV突变和体系SNV突变；
插入片段长度计算单元，用于计算肿瘤样本的插入片段长度分布平均值μ和标准差σ，按照如下公式计算得到覆盖99.9％的插入片段最大长度L，L＝μ+3σ；
合并连锁突变单元，用于根据L将所有胚系和体系SNV突变合并成潜在的LV突变组，每组中从第一个体系SNV突变开始，从与之距离≤L的体系SNV集合依次抽出1个体系SNV，构成一对潜在2连锁SNV突变，迭代上述步骤最后得到一系列非重复的潜在2连锁SNV突变；从与之距离≤L的体系SNV集合依次抽出2个体系SNV，迭代上述步骤最后得到一系列非重复的潜在3连锁SNV突变；将潜在2连锁的SNV突变，从于与之距离≤L的胚系SNV集合抽出一个胚系SNV，构成潜在3连锁SNV突变；所有的潜在连锁突变，如果某一LV突变i，其染色体编号为chrom0，起始的SNV突变位置是start0，终止的SNV突变位置是end0，另一LV突变j，其染色体编号是chrom1，起始SNV突变位置是start1，终止SNV突变位置是end1，chrom0＝chrom1并且start1≥start0并且end1≤end0，则将j纳入i突变区间内；
突变矩阵构建单元，用于针对所有合并的潜在LV突变区域，构建连锁突变矩阵；
连锁突变信息计算单元，用于按照上述区域并行计算肿瘤和正常对应LV的深度DP和突变支持数AD；
筛选连锁突变单元，用于根据权利要求1所述的方法得到的数据库，过滤掉噪音连锁突变；
所述突变矩阵构建单元，具体用于：
针对所有合并的潜在LV突变区域，构建如下矩阵，假设该区域包含从编号为X1开始到编号为Xi的LV突变，假设该区域染色体编号为变量Y，从位置Z1开始，一直到位置Zj为止，通过如下碱基数字对应关系：
将每个潜在LV突变区域转为i×j的矩阵Mat：
该矩阵每一行代表一个潜在LV突变，Basenij代表该位置SNV碱基编号；
同时也构建一个包含i个元素的向量VecLV：
该向量的每个元素SNi代表潜在LV突变对应的SNV个数；
所述连锁突变信息计算单元，具体用于：
按照突变区域并行计算肿瘤和正常对应LV的深度DP和突变支持数AD，主要方法如下：
1)构建长度为i的0元素向量Vec0，分别赋值给Vecdp和Vecad，Vecdp和Vecad分别用以存储对应LV的深度信息和突变支持数信息；
2)迭代比对到染色体编号为变量Y，起始位置Z1和Zj之间所有测序读段，过滤掉标签为重复、补充比对、次比对、比对异常、MQ≤20和错配碱基数目大于7个的测序读段，根据上述碱基数字转化公式以及测序读段的参考基因组比对位置，将测序读段序列转化为长度为j的向量Vecse，如果测序读段对应的参考碱基是插入突变，则跳过该测序读段上位置；
3)将单端测序读段1的向量Vecse1和单端测序读段2的向量Vecse2按照如下公式合并为双端测序读段的向量Vecpe：
Vecse1＝[R1Basen1 … R1Basenj]
Vecse2＝[R2Basen1 … R2Basenj]
Vecpe＝[PRBasen1 … RRBasenj]
4)将矩阵Mat和Vecpe按照如下公式进行计算，得到Matpe:
5)根据Matpe，按照如下公式计算得到向量Vecdppe和向量Vecadpe：
6)按照下述方法，计算Vecdp和Vecad数值变化，假设加入Matpe前，Vecdp和Vecad第i位数值分别是DPi,j―1和ADi,j―1，加入Matpe计算之后，数值分别是DPi,j和ADi,j，则DPi,j和DPi,j―1的关系如下：
ADi,j和ADi,j―1的关系如下：
7)迭代该区域对应的所有配对测序读段，可得到每个潜在LV对应的深度和突变支持数矩阵，分别处理肿瘤和正常的样本比对数据，最后得到肿瘤和正常对应的所有LV深度、突变支持数和突变频率AF＝AD/DP；
所述筛选连锁突变单元，具体用于：
通过对比肿瘤和正常的LV突变信息，使用fisher检验获得可能的致病LV突变，然后按照下述条件筛选监控致病突变：大于等于一定突变频率、测序深度和突变支持数，LV至少一个SNV突变位于探针目标区域前后延伸50bp内；
如果某N个SNV连锁突变，和黑名单数据库中SNV有交集，并且交集个数大于等于N-1，所述N-1个SNV突变的黑名单人群频率大于等于5％并且突变频率大于等于1％，则将该连锁突变移除监控位点列表；
如果某个LV，存在黑名单数据库中，并且该LV人群频率大于等于5％，通过单边格拉布斯检验法判断该LV是否显著大于黑名单数据库中LV突变频率，如果不是，则将该LV移除监控位点列表。
3.一种MRD监测装置，其特征在于，包括：依次连接的信息采集单元、数据构建单元和检验单元；
信息采集单元，用于首先提取患者综合治疗之后的血液中cfDNA，然后使用探针进行捕获建库及测序；测序之后，通过比对标记重复获得比对信息；使用权利要求1所述的方法获得所有可能的连锁突变信息；
数据构建单元，用于提取初诊监控位点对应的突变信息，假设为k个，从剩余突变中筛选可能的噪音突变，每个监控位点突变筛选50个相应的噪音突变，噪音连锁突变的每个SNV突变类型需要和监控的连锁突变染色体编号一致，突变类型一致，并且每个SNV突变距离需要和监控连锁突变的SNV距离一致，最后得到所有监控位对应的个噪音突变；
检验单元，用于统计监控位点对应的突变频率平均值μ0，完成10000次抽样，每次抽样依次从k个监测位点对应的噪音位点随机抽样，如从第i个监控位点的噪音突变{NSi,1,NSi,2…,NSi,50}随机抽一个，然后计算每次抽选的k个噪音突变的平均值μi，最后得到了{μ1,μ2…,μ10000}所有抽样结果，通过单边格拉布斯检验法判断μ0是否显著大于{μ1,μ2…,μ10000}，得到检验P值；当P<0.05时，判定MRD为阳性，定义MRD样本的监控位点对应的LV突变分子数目为x，监控位点总的测序深度为d，则样本的肿瘤占比为λ，λ等于x/d；当P≥0.05时，判定MRD结果阴性或低于检测下限；
所述MRD表示微小残留病灶。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310063877.4A CN116064755B (zh) | 2023-01-12 | 2023-01-12 | 一种基于连锁基因突变检测mrd标志物的装置 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310063877.4A CN116064755B (zh) | 2023-01-12 | 2023-01-12 | 一种基于连锁基因突变检测mrd标志物的装置 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN116064755A CN116064755A (zh) | 2023-05-05 |
| CN116064755B true CN116064755B (zh) | 2023-10-20 |
Family
ID=86169500
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310063877.4A Active CN116064755B (zh) | 2023-01-12 | 2023-01-12 | 一种基于连锁基因突变检测mrd标志物的装置 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116064755B (zh) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116959579B (zh) * | 2023-09-19 | 2023-12-22 | 北京求臻医学检验实验室有限公司 | 一种用于降低二代测序系统错误的系统 |
| CN117373527A (zh) * | 2023-12-07 | 2024-01-09 | 中国科学院微生物研究所 | Hiv序列质控方法、设备及存储介质 |
Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103221551A (zh) * | 2010-11-23 | 2013-07-24 | 深圳华大基因科技有限公司 | Hla基因型别-snp连锁数据库、其构建方法、以及hla分型方法 |
| WO2016078067A1 (zh) * | 2014-11-21 | 2016-05-26 | 深圳华大基因研究院 | 个体单核苷酸多态性位点分型方法及装置 |
| CN107002093A (zh) * | 2014-06-26 | 2017-08-01 | 瑞泽恩制药公司 | 用于靶向遗传修饰的方法和组合物，以及这些组合物的使用方法 |
| CN113174431A (zh) * | 2021-03-11 | 2021-07-27 | 厦门市妇幼保健院（厦门市计划生育服务中心） | Topovibl作为靶标在诊断和治疗非梗阻性无精子症中的应用 |
| CN113284554A (zh) * | 2021-04-28 | 2021-08-20 | 中山大学肿瘤防治中心（中山大学附属肿瘤医院、中山大学肿瘤研究所） | 一种筛查结直肠癌术后微小残留病灶及预测复发风险的循环肿瘤dna检测系统及应用 |
| CN114187964A (zh) * | 2021-12-13 | 2022-03-15 | 深圳市海普洛斯生物科技有限公司 | 一种肺癌围手术期分子残留病灶基因检测panel及检测模型的构建方法 |
| CN114292912A (zh) * | 2021-12-24 | 2022-04-08 | 广州燃石医学检验所有限公司 | 一种变体核酸的检测方法 |
| CN115491423A (zh) * | 2022-09-23 | 2022-12-20 | 珠海横琴铂华医学检验有限公司 | 一种用于b细胞淋巴瘤mrd监测的基因组合、试剂盒与应用 |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2019055835A1 (en) * | 2017-09-15 | 2019-03-21 | The Regents Of The University Of California | DETECTION OF SOMATIC MONONUCLEOTIDE VARIANTS FROM ACELLULAR NUCLEIC ACID WITH APPLICATION TO MINIMUM RESIDUAL DISEASE SURVEILLANCE |
| US11781189B2 (en) * | 2019-08-27 | 2023-10-10 | Fundación Para La Investigación Biomédica Del Hospital Universitario 12 De Octubre | Method for determining the presence or absence of minimal residual disease (MRD) in a subject who has been treated for a disease |
-
2023
- 2023-01-12 CN CN202310063877.4A patent/CN116064755B/zh active Active
Patent Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103221551A (zh) * | 2010-11-23 | 2013-07-24 | 深圳华大基因科技有限公司 | Hla基因型别-snp连锁数据库、其构建方法、以及hla分型方法 |
| CN107002093A (zh) * | 2014-06-26 | 2017-08-01 | 瑞泽恩制药公司 | 用于靶向遗传修饰的方法和组合物，以及这些组合物的使用方法 |
| WO2016078067A1 (zh) * | 2014-11-21 | 2016-05-26 | 深圳华大基因研究院 | 个体单核苷酸多态性位点分型方法及装置 |
| CN113174431A (zh) * | 2021-03-11 | 2021-07-27 | 厦门市妇幼保健院（厦门市计划生育服务中心） | Topovibl作为靶标在诊断和治疗非梗阻性无精子症中的应用 |
| CN113284554A (zh) * | 2021-04-28 | 2021-08-20 | 中山大学肿瘤防治中心（中山大学附属肿瘤医院、中山大学肿瘤研究所） | 一种筛查结直肠癌术后微小残留病灶及预测复发风险的循环肿瘤dna检测系统及应用 |
| CN114187964A (zh) * | 2021-12-13 | 2022-03-15 | 深圳市海普洛斯生物科技有限公司 | 一种肺癌围手术期分子残留病灶基因检测panel及检测模型的构建方法 |
| CN114292912A (zh) * | 2021-12-24 | 2022-04-08 | 广州燃石医学检验所有限公司 | 一种变体核酸的检测方法 |
| CN115491423A (zh) * | 2022-09-23 | 2022-12-20 | 珠海横琴铂华医学检验有限公司 | 一种用于b细胞淋巴瘤mrd监测的基因组合、试剂盒与应用 |
Non-Patent Citations (3)
| Title |
|---|
| Circulating Tumor DNA and Minimal Residual Disease (MRD) in Solid Tumors: Current Horizons and Future Perspectives;Yan Peng等;Front Oncol.;第1-14页 * |
| 白血病基因诊断及其信息系统建立;孙新六;中国实验诊断学;第1182-1184页 * |
| 过继细胞转输治疗白血病微小残留病灶治疗效果的Meta分析;张庆义;世界最新医学信息文摘;第88-89页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN116064755A (zh) | 2023-05-05 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| JP7051900B2 (ja) | 不均一分子長を有するユニーク分子インデックスセットの生成およびエラー補正のための方法およびシステム | |
| CN116064755B (zh) | 一种基于连锁基因突变检测mrd标志物的装置 | |
| CN110800063B (zh) | 使用无细胞dna片段大小检测肿瘤相关变体 | |
| CN104781421B (zh) | 检测稀有突变和拷贝数变异的系统和方法 | |
| EP2926288B1 (en) | Accurate and fast mapping of targeted sequencing reads | |
| CN112086129B (zh) | 预测肿瘤组织cfDNA的方法及系统 | |
| CN106156543B (zh) | 一种肿瘤ctDNA信息统计方法 | |
| CN114574581A (zh) | 检测稀有突变和拷贝数变异的系统和方法 | |
| CN113903401B (zh) | 基于ctDNA长度的分析方法和系统 | |
| CN107523563A (zh) | 一种用于循环肿瘤dna分析的生物信息处理方法 | |
| CN105132407B (zh) | 一种脱落细胞dna低频突变富集测序方法 | |
| CN105525033A (zh) | 检测血液中微生物的方法及装置 | |
| CN114187964A (zh) | 一种肺癌围手术期分子残留病灶基因检测panel及检测模型的构建方法 | |
| CN115083521B (zh) | 一种单细胞转录组测序数据中肿瘤细胞类群的鉴定方法及系统 | |
| CN114694750A (zh) | 一种基于ngs平台的单样本肿瘤体细胞突变判别及tmb检测方法 | |
| WO2021061473A1 (en) | Systems and methods for diagnosing a disease condition using on-target and off-target sequencing data | |
| CN116356001B (zh) | 一种基于血液循环肿瘤dna的双重背景噪声突变去除方法 | |
| CN113862351B (zh) | 体液样本中鉴定胞外rna生物标志物的试剂盒及方法 | |
| CN105779435A (zh) | 试剂盒及其用途 | |
| CN117275585A (zh) | 基于lp-wgs和dna甲基化的肺癌早筛模型构建方法及电子设备 | |
| CN113362893A (zh) | 肿瘤筛查模型的构建方法及应用 | |
| CN112687341B (zh) | 一种以断点为中心的染色体结构变异鉴定方法 | |
| CN115331737A (zh) | 一种分析肠道菌群中致病菌和量化菌群地域特征的方法 | |
| CN114613436B (zh) | 血样Motif特征提取方法及癌症早筛模型构建方法 | |
| RU2766198C9 (ru) | Способы и системы для получения наборов уникальных молекулярных индексов с гетерогенной длиной молекул и коррекции в них ошибок |
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 | ||
| CB03 | Change of inventor or designer information |
Inventor after: Xiao Min Inventor after: Zhu Zhoujie Inventor after: Zhang Wei Inventor after: Yang Xiaoxia Inventor after: Shen Kefeng Inventor after: Zhang Meilan Inventor after: Mu Wei Inventor before: Xiao Min Inventor before: Zhu Zhoujie Inventor before: Zhang Wei Inventor before: Yang Xiaolu Inventor before: Shen Kefeng Inventor before: Zhang Meilan Inventor before: Mu Wei |
|
| CB03 | Change of inventor or designer information | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |