CN111667883A

CN111667883A - 一种基于复合微单倍型焦磷酸测序图谱解析的法医学混合dna的分析方法

Info

Publication number: CN111667883A
Application number: CN202010493118.8A
Authority: CN
Inventors: 张; 王玉芳; 曹悦岩; 宋凤; 罗海玻; 陈晓刚
Original assignee: Sichuan University
Current assignee: Sichuan University
Priority date: 2020-06-03
Filing date: 2020-06-03
Publication date: 2020-09-15
Anticipated expiration: 2040-06-03
Also published as: CN111667883B

Abstract

本发明提供了一种基于复合微单倍型焦磷酸测序图谱解析的法医学混合DNA的分析方法，属于法医遗传学和分子生物学技术领域，包括：提取样本的DNA；对样本DNA进行荧光定量，得到定量样本DNA，将所述定量样本DNA进行PCR；对PCR产物进行焦磷酸测序，得到微单倍型的基因分型结果；使用AdvISER‑M‑PYRO算法对微单倍型的基因分型结果进行解析，得到贡献系数、相关系数和贡献者；使用支持向量机算法对微单倍型的基因分型结果进行分型，对贡献者的构成进行分类。本发明获得微单倍型等位基因分型，计算出各位点混合DNA检测能力，获得个体相应微单倍型构成，获得混合DNA贡献者构成，获得混合DNA贡献者的混合比。

Description

一种基于复合微单倍型焦磷酸测序图谱解析的法医学混合 DNA的分析方法

技术领域

本发明属于法医遗传学和分子生物学技术领域，尤其涉及一种基于复合微单倍型焦磷酸测序图谱解析的法医学混合DNA的分析方法。

背景技术

混合DNA是指待测样本中包含了两个及以上个体DNA的样本。混合DNA的检测和分析在法医学领域非常重要且极具挑战。随着当前STR扩增试剂盒以及毛细管电泳技术检测灵敏度逐渐提高，混合DNA中的次要成分检出率也得以提升，使得利用原有技术无法发现的混合DNA得以检出。

由于正常人类是二倍体，每个STR基因座上最多含有两个不同的等位基因。当出现以下现象时，提示该检测样本为混合DNA：在多个STR基因座上出现了两个以上的等位基因峰；多个STR基因座上虽然只有两个等位基因峰，但峰高比严重失衡；牙釉蛋白基因座上出现了X和Y染色体等位基因峰的峰高严重失衡；stutter峰的峰高异常等。使用毛细管电泳技术检测STR扩增产物的方法(STR-CE)在混合DNA的检测和分析中发挥了很大的作用，但同时也存在伪峰(artefact)，等位基因或位点丢失(drop-out)以及共享等位基因(allelicsharing)等现象，以上这些现象均会导致STR位点上峰的数量与真实等位基因峰的数量不符，从而增加混合DNA的检测和分析难度。

微单倍型(microhaplotypes，microhaps或MHs)是由Kidd等提出的一种新型遗传标记，可以作为检测混合DNA样本的一种新的思路。微单倍型位点的片段范围通常小于300个核苷酸，一般由两个或更多个紧密连接的单核苷酸多态性(single nucleotidepolymorphisms,SNPs)组成，能够形成三个和(或)以上的SNP组合。由于一个微单倍型内部的SNP之间物理距离很近，因此SNP之间的重组率极低，分子上非常接近(<10kb)的SNP之间重组率低于10^-4。由于微单倍型是由多个SNP构成的单倍型，它们可以在每个位点基础上提供比单独的SNP标记更丰富的多态性。一直以来，Sanger测序是DNA测序的主要方法，但是当两个或多个SNP位点是杂合状态时，传统的Sanger测序则无法从基因组DNA中确定单个SNP等位基因分别属于哪条染色体，也就是实际SNP单倍型测定(亦称为phase的过程)。随着新一代测序技术(next generation sequencing,NGS)的发展，使原本不能通过Sanger测序定义phase的高密度紧密排列的SNP位点，可以通过NGS的“单分子测序”进行排序，从而在一个特定单倍型位点区分其等位基因(实际单倍型)。受限于NGS测序读长，目前已开发的SNP微单倍型位点的扩增片段大小一般不超过400bp。微单倍型的杂合度水平取决于不同因素，包括目标区域内不同位置的等位基因突变的历史积累，罕见重组体的出现，随机遗传漂移发生和/或选择。微单倍型的多态性与该区域内包含的SNP位点数量和多态性有关。

相较于传统的STR，由于构成微单倍型的SNP没有短串联重复序列结构，在PCR过程中可以避免出现复制滑脱，从而不会出现stutter这样的干扰伪峰。Stutter片段往往会增加贡献者比例不均衡的混合DNA的分析难度，尤其是当stutter峰的峰高与次要成分贡献者的等位基因峰高相当和/或stutter峰的片段与贡献者真实等位基因峰发生重叠时。因此，微单倍型这一遗传标记不会出现stutter所带来的干扰。与STR相比，微单倍型的另一个优势在于同一个微单倍型位点的等位基因都具有相同的片段长度。而同一个STR位点，其等位基因的片段长度差异可达100nt。若DNA降解或DNA中存在PCR抑制物，那么片段较小的等位基因会发生优势扩增，扩增效能明显高于同一个贡献者的其他片段较大的等位基因，这会增加混合DNA的分析难度。

然而，微单倍型与STR相比也有其自身的弊端，由于微单倍型作为一种新的遗传标记，其数据库的构建还不够完善，目前ALFRED数据库(https://alfred.med.yale.edu/alfred/index.asp)收录了198个微单倍型在部分参考群体中的基因频率，因此微单倍型的群体数据有待进一步的扩充。此外，微单倍型在不同群体间的等位基因频率差异很大，因此在计算随机匹配概率或似然率时需要获得群体特异性的等位基因频率。微单倍型的另一个缺点是其需要使用新一代测序这一工作流程十分耗时耗力的分型方法。尽管每个碱基的成本相对较低，但是较之毛细管电泳平台仍然具有很高的成本，而且需要复杂的流程进行分析，这可能会限制微单倍型在全球法医DNA分析领域的快速推广。

发明内容

本发明的目的在于提供一种基于复合微单倍型焦磷酸测序图谱解析的法医学混合DNA的分析方法。

为了实现上述发明目的，本发明提供了以下技术方案：

本发明提供了一种基于复合微单倍型焦磷酸测序图谱解析的法医学混合DNA的分析方法，包括以下步骤：

1)提取样本的DNA，得到样本DNA；

2)对所述步骤1)得到的样本DNA进行荧光定量，得到定量样本DNA，将所述定量样本DNA进行PCR，得到PCR产物；

3)对所述步骤2)得到的PCR产物进行焦磷酸测序，得到微单倍型的基因分型结果；

4)使用AdvISER-M-PYRO算法对步骤3)得到的微单倍型的基因分型结果进行解析，得到微单倍型的基因分型结果、贡献者构成、贡献系数和相关系数；

5)使用支持向量机算法对步骤3)得到的微单倍型的基因分型结果进行分类，得到每个微单倍型的基因分型构成，得到每个混合DNA的贡献者构成；

所述步骤2)PCR使用的引物包括SEQ ID No.1～8所示的核苷酸序列；所述SEQ IDNo.1～2扩增微单倍型mh03ZJ-001；所述SEQ ID No.3～4扩增微单倍型mh06ZJ-001；所述SEQ ID No.5～6扩增微单倍型mh07ZJ-001；所述SEQ ID No.7～8扩增微单倍型mh19ZJ-001；

所述步骤3)焦磷酸测序使用的引物包括SEQ ID No.9～12所示的核苷酸序列；所述SEQ ID No.9测序微单倍型mh03ZJ-001；所述SEQ ID No.10测序微单倍型mh06ZJ-001；所述SEQ ID No.11测序微单倍型mh07ZJ-001；所述SEQ ID No.12测序微单倍型mh19ZJ-001。

优选的，所述步骤3)焦磷酸测序，当所述PCR产物中仅含有微单倍型mh03ZJ-001时，所述焦磷酸测序的体系每50μl包括：2×Go Taq Colorless Master Mix 25μl、浓度为10μM的SEQ IDNo.9 1.2μl、无核酸酶水21.8μl和浓度为5ng/μL的PCR产物2μl。

优选的，所述步骤3)焦磷酸测序，当所述PCR产物中仅含有微单倍型mh06ZJ-001时，所述焦磷酸测序的体系每50μl包括：2×Go Taq Colorless Master Mix 25μl、浓度为10μM的SEQ IDNo.10 1μl、无核酸酶水22μl和浓度为5ng/μL的PCR产物2μl。

优选的，所述步骤3)焦磷酸测序，当所述PCR产物中仅含有微单倍型mh07ZJ-001时，所述焦磷酸测序的体系每50μl包括：2×Go Taq Colorless Master Mix 25μl、浓度为10μM的SEQ IDNo.11 1.4μl、无核酸酶水21.6μl和浓度为5ng/μL的PCR产物2μl。

优选的，所述步骤3)焦磷酸测序，当所述PCR产物中仅含有微单倍型mh19ZJ-001时，所述焦磷酸测序的体系每50μl包括：2×Go Taq Colorless Master Mix 25μl、浓度为10μM的SEQ IDNo.12 0.8μl、无核酸酶水22.2μl和浓度为5ng/μL的PCR产物2μl。

优选的，制备焦磷酸测序使用的模板的程序为：预变性95℃2min；Touchdown阶段95℃30s，71℃45s，72℃30s，20个循环，每一个循环71℃下降0.5℃；循环扩增阶段95℃30s，61℃45s，72℃30s，25个循环；终延伸72℃5min。

优选的，当所述PCR产物中同时含有微单倍型mh03ZJ-001、mh06ZJ-001、mh07ZJ-001和mh19ZJ-001时，所述焦磷酸测序的体系每50μl包括：2×Go Taq Colorless MasterMix 25μl、浓度为0.24μM的SEQ IDNo.1～2各1.2μl、浓度为0.2μM的SEQ IDNo.3～4各1μl、浓度为0.28μM的SEQ IDNo.5～6各1.4μl、浓度为0.16μM的SEQ IDNo.7～8各0.8μl、无核酸酶水14.2μl和浓度为5ng/μL的PCR产物2μl。

优选的，所述焦磷酸测序的程序为预变性95℃2min；Touchdown阶段95℃30s，71℃45s，72℃30s，20个循环，每一个循环71℃下降0.5℃；循环扩增阶段95℃30s，61℃45s，72℃30s，25个循环；终延伸72℃5min。

本发明提供了一种基于复合微单倍型焦磷酸测序图谱解析的法医学混合DNA的分析方法，包括以下步骤：1)提取样本的DNA，得到样本DNA；

2)对所述步骤1)得到的样本DNA进行荧光定量，得到定量样本DNA，将所述定量样本DNA进行PCR，得到PCR产物；3)对所述步骤2)得到的PCR产物进行焦磷酸测序，分别得到每个微单倍型的基因分型结果；4)使用AdvISER-M-PYRO算法对步骤3)得到的每个微单倍型的基因分型结果进行解析，得到微单倍型的基因分型结果、贡献者构成、贡献系数和相关系数；5)使用支持向量机算法对步骤3)得到的微单倍型的基因分型结果进行分类，得到每个微单倍型的基因分型构成，得到每个混合DNA的贡献者构成；所述步骤2)PCR使用的引物包括SEQ ID No.1～8所示的核苷酸序列；所述SEQ ID No.1～2扩增微单倍型mh03ZJ-001；所述SEQ ID No.3～4扩增微单倍型mh06ZJ-001；所述SEQ ID No.5～6扩增微单倍型mh07ZJ-001；所述SEQ ID No.7～8扩增微单倍型mh19ZJ-001；所述步骤3)焦磷酸测序使用的引物包括SEQ ID No.9～12所示的核苷酸序列；所述SEQ ID No.9测序微单倍型mh03ZJ-001；所述SEQ ID No.10测序微单倍型mh06ZJ-001；所述SEQ ID No.11测序微单倍型mh07ZJ-001；所述SEQ ID No.12测序微单倍型mh19ZJ-001。

本发明的有益效果是：

(1)获得微单倍型等位基因(实际单倍型)分型，无需使用大规模平行测序或者单倍型分型(haplotypephasing)软件计算其理论单倍型分型；

(2)能够根据已有的微单倍型等位基因频率，计算出各位点混合DNA检测能力(Probabilities of Detecting aMixture,PDM)，并且可以评估位点联合使用时的累计概率；

(3)可以根据单一个体复合微单倍型焦磷酸测序结果，获得该个体相应微单倍型构成；

(4)该方法对单一DNA复合微单倍型进行焦磷酸测序的灵敏度可达40pg/μL；

(5)识别出由2,3,4个个体组成的混合DNA贡献者构成。以此为基础，可对混合DNA中贡献者数目进行推断。

(6)识别两个体不均衡混合DNA中贡献者构成，最大混合比达1:30，且混合比与贡献者的贡献系数之比存在良好的相关性，回归直线对观测值具有较高的拟合程度(见图5，R²＝0.9864)。以此为基础，本方法可对混合DNA中贡献者DNA之间的比例进行推断。

附图说明

图1为四个微单倍型的所有单倍型分型所得焦磷酸测序图谱，其中(a)-(e)分别表示微单倍型mh03ZJ-001测得的5种微单倍型基因型，(f)-(k)分别表示微单倍型mh06ZJ-001测得的6种微单倍型基因型，(l)-(n)分别表示微单倍型mh07ZJ-001测得3种微单倍型基因型，(o)-(q)分别表示微单倍型mh19ZJ-001测得的3种微单倍型基因型；

图2为复合微单倍型的单一DNA模板焦磷酸测序图谱，(a)为61号DNA复合微单倍型焦磷酸测序图谱，(b)为76号DNA复合微单倍型焦磷酸测序图谱，(c)为79号DNA复合微单倍型焦磷酸测序图谱；

图3为两个体等比混合DNA复合微单倍型测序图谱，(a)图为33和75号构成的等比混合样本，标记为M33-75；(b)图为61和79号构成的等比混合样本，标记为M61-79；(c)图为61和65号构成的等比混合样本，标记为M61-65；

图4为两个体不等比混合DNA复合微单倍型测序图谱，(a)-(d)均表示33和75号个体混合的焦磷酸测序结果，混合比分别为1:5,1:10,1:20,1:30(33号个体为次要成分，75号个体为主要成分),依次标记为M1:5，M1:10，M1:20和M1:30；

图5为两贡献者贡献系数比与混合比之间的回归分析，X轴为两个体初始DNA模板量的混合比，Y轴为解析后各贡献者的贡献系数比值；

图6为三个体等比混合DNA复合微单倍型测序图谱，(a)为33,75和61号个体等比混合的焦磷酸测序图谱，标记为M33-75-61；(b)为76,61和65号个体等比混合的焦磷酸测序图谱，标记为M76-61-65；

图7为四个体等比混合DNA复合微单倍型测序图谱，(a)为33,61,65和75号DNA等比混合的焦磷酸测序图谱，标记为M33-61-65-75；(b)为33,61,65和76号DNA等比混合的焦磷酸测序图谱，标记为M33-61-65-76。

具体实施方式

1)提取样本的DNA，得到样本DNA；

3)对所述步骤2)得到的PCR产物进行焦磷酸测序，分别得到每个微单倍型的基因分型结果；

4)使用AdvISER-M-PYRO算法对步骤3)得到的每个微单倍型的基因分型结果进行解析，得到微单倍型的基因分型结果、贡献者构成、贡献系数和相关系数；

本发明对提取样本的DNA的方法没有特殊限定，采用常规方法即可，在本发明具体实施方式中优选采用天根生化科技(北京)有限公司销售的血液/细胞/组织基因组提取试剂盒，按照试剂盒的说明书提取即可。

本发明将得到的样本DNA进行荧光定量，得到定量样本DNA，将所述定量样本DNA进行PCR，得到PCR产物。所述PCR使用的引物包括SEQ ID No.1～8所示的核苷酸序列；所述SEQID No.1～2扩增微单倍型mh03ZJ-001；所述SEQ ID No.3～4扩增微单倍型mh06ZJ-001；所述SEQ ID No.5～6扩增微单倍型mh07ZJ-001；所述SEQ ID No.7～8扩增微单倍型mh19ZJ-001。在本发明中，所述焦磷酸测序使用的引物包括SEQ ID No.9～12所示的核苷酸序列；所述SEQ ID No.9测序微单倍型mh03ZJ-001；所述SEQ ID No.10测序微单倍型mh06ZJ-001；所述SEQ ID No.11测序微单倍型mh07ZJ-001；所述SEQ ID No.12测序微单倍型mh19ZJ-001。所述PCR扩增使用的引物和测序使用的引物的具体序列见表1。在本发明中，所述PCR使用的引物优选以生物素进行标记。

表1 PCR引物与测序引物的序列信息

在本发明中，所述焦磷酸测序，当所述PCR产物中仅含有微单倍型mh03ZJ-001时，所述焦磷酸测序的体系每50μl优选包括：2×Go Taq Colorless Master Mix 25μl、浓度为10μM的SEQ IDNo.91.2μl、无核酸酶水21.8μl和浓度为5ng/μL的PCR产物2μl。在本发明中，所述焦磷酸测序的程序优选为：预变性95℃2min；Touchdown阶段95℃30s，71℃45s，72℃30s，20个循环，每一个循环71℃下降0.5℃；循环扩增阶段95℃30s，61℃45s，72℃30s，25个循环；终延伸72℃5min。

在本发明中，所述焦磷酸测序，当所述PCR产物中仅含有微单倍型mh06ZJ-001时，所述焦磷酸测序的体系每50μl优选包括：2×Go Taq Colorless Master Mix 25μl、浓度为10μM的SEQ IDNo.101μl、无核酸酶水22μl和浓度为5ng/μL的PCR产物2μl。在本发明中，所述焦磷酸测序的程序优选为：预变性95℃2min；Touchdown阶段95℃30s，71℃45s，72℃30s，20个循环，每一个循环71℃下降0.5℃；循环扩增阶段95℃30s，61℃45s，72℃30s，25个循环；终延伸72℃5min。

在本发明中，所述焦磷酸测序，当所述PCR产物中仅含有微单倍型mh07ZJ-001时，所述焦磷酸测序的体系每50μl优选包括：2×Go Taq Colorless Master Mix 25μl、浓度为10μM的SEQ IDNo.11 1.4μl、无核酸酶水21.6μl和浓度为5ng/μL的PCR产物2μl。在本发明中，所述焦磷酸测序的程序优选为：预变性95℃2min；Touchdown阶段95℃30s，71℃45s，72℃30s，20个循环，每一个循环71℃下降0.5℃；循环扩增阶段95℃30s，61℃45s，72℃30s，25个循环；终延伸72℃5min。

在本发明中，所述焦磷酸测序，当所述PCR产物中仅含有微单倍型mh19ZJ-001时，所述焦磷酸测序的体系每50μl优选包括：2×Go Taq Colorless Master Mix 25μl、浓度为10μM的SEQ IDNo.12 0.8μl、无核酸酶水22.2μl和浓度为5ng/μL的PCR产物2μl。在本发明中，所述焦磷酸测序的程序优选为：预变性95℃2min；Touchdown阶段95℃30s，71℃45s，72℃30s，20个循环，每一个循环71℃下降0.5℃；循环扩增阶段95℃30s，61℃45s，72℃30s，25个循环；终延伸72℃5min。

在本发明中，当所述PCR产物中同时含有微单倍型mh03ZJ-001、mh06ZJ-001、mh07ZJ-001和mh19ZJ-001时，所述焦磷酸测序的体系每50μl优选包括：2×Go TaqColorless Master Mix 25μl、浓度为0.24μM的SEQ IDNo.1～2各1.2μl、浓度为0.2μM的SEQIDNo.3～4各1μl、浓度为0.28μM的SEQ IDNo.5～6各1.4μl、浓度为0.16μM的SEQ IDNo.7～8各0.8μl、无核酸酶水14.2μl和浓度为5ng/μL的PCR产物2μl。在本发明中，所述焦磷酸测序的程序优选为预变性95℃2min；Touchdown阶段95℃30s，71℃45s，72℃30s，20个循环，每一个循环71℃下降0.5℃；循环扩增阶段95℃30s，61℃45s，72℃30s，25个循环；终延伸72℃5min。

在本发明中，所述PCR产物优选按照说明书(

Q96 ID，QIAGEN，德国)进行处理，处理后再进行焦磷酸测序。

本发明使用AdvISER-M-PYRO算法对步骤3)得到的每个微单倍型的基因分型结果进行解析，得到微单倍型的基因分型结果、贡献者构成、贡献系数和相关系数。拥有贡献系数的项为可能的微单倍型。使用该算法对混合DNA贡献者构成进行解析，贡献系数最高者为混合DNA中最有可能的贡献者来源。

本发明使用支持向量机算法对步骤3)得到的微单倍型的基因分型结果进行分类，得到每个微单倍型的基因分型构成，得到每个混合DNA的贡献者构成。

下面结合实施例对本发明提供的技术方案进行详细的说明，但是不能把它们理解为对本发明保护范围的限定。

以下实施例用到的仪器、试剂和软件：

a.实验所用仪器：

Eppendorf Centrifuge 5427 R离心机(Eppendorf，德国)

Applied Biosystems^TM 7500(Thermo Fisher Scientific，美国)

Veriti Thermal Cycler PCR仪(Applied Biosystems，美国)

Pyrosequencing Q96ID焦磷酸测序仪(QIAGEN，德国)

Vacuum Prep Tool焦磷酸测序样本处理制备平台(QIAGEN，德国)

Q96 Plate Low焦磷酸反应板(QIAGEN，德国)

恒温混匀仪(Eppendorf，德国)

单模块金属浴(Thermo Fisher Scientific，美国)

b.试剂：

血液/细胞/组织基因组提取试剂盒(天根生化科技(北京)有限公司)

Investigator Quantiplex Pro Kit(QIAGEN，德国)

链霉亲和素包被的琼脂糖珠(GE医疗集团，美国)

Tris-base、NaCl、EDTA、Tween-20、HCl、四水合乙酸镁、冰乙酸、NaOH、70％乙醇。

PyroMark GoldQ96 SQA Reagents(1×96)试剂盒(QIAGEN，德国)，包括：

Enzyme Mixture，

Substrate Mixture，

dNTP：dATP、dCTP、dGTP、dTTP

c.软件：

Q96 software(1.0.6版本，QIAGEN，德国)

PyroMark Assay Design(2.0版本，QIAGEN，德国)

AutoDimer(1.0版本，STRBase，美国)

Rstudio R软件包：AdvlSER-M-PYRO(Ambroise等人开发)

PyroMaker(1.1版本，Alan O’Neill等人开发)

Python(3.6版本，美国)

PyCharm(2019.3.3版本，JetBrains，捷克)

Modified-Powerstates(赵方等人开发)

实施例1

(1)使用shell脚本(https://github.com/YuancunZhao/vcfscripts/blob/master/bcftools_genotype_output.sh)通过循环迭代调用vcftools实现快速输出指定群体(1000Genomes中全球26个群体中的任意群体)和位点的基因型，获得CHB群体的SNPs位点的基因型。步骤(1)是从1000Genomes中根据基因频率对SNP进行初步筛选。

(2)根据SNP最小等位基因频率(Minor Allele Frequency,MAF)>0.4，以及SNPs位点间的物理距离<30bp获得MAF>0.4这两个标准，通过PYTHON编写的脚本针对每一条染色体进一步筛选参与构建微单倍型的SNPs。步骤(2)是将初筛的SNP根据其物理距离进行进一步筛选。

(3)通过PYTHON编写的脚本，使用上一步骤中获得的SNPs构建微单倍型。计算微单倍型组合在1000Genome数据库的CHB群体中频率分布情况，计算Ae值及个人识别概率，以此作为筛选微单倍型的标准，从而获得候选微单倍型。步骤(3)是获得微单倍型并计算其相关参数。

(4)根据已有的微单倍型等位基因频率，计算出各位点混合DNA检测能力(Probabilities of Detecting a Mixture,PDM)，并且评估位点联合使用时的累计概率。步骤(4)是计算各个微单倍型的混合DNA检测能力，并且评估位点联合使用时的累计概率。

(5)使用PyroMark Assay Design软件(2.0版本，Qiagen)设计PCR引物及焦磷酸测序引物。使用AutoDimer软件对以上PCR引物及测序引物进行测试，检验复合扩增和复合焦磷酸测序体系中各引物间是否存在相互作用以及引物自身是否存在发夹结构。所有引物的详细序列信息均列于表1中。其中，F/R表示正上游/下游扩增引物，S表示测序引物；在扩增引物的上游或下游引物上标记了生物素并以“*”标识。见表1。步骤(5)设计各微单倍型的扩增引物和测序引物。

(6)通过调用AdvISER-M-PYRO程序包中的“senator”函数，将每个微单倍型所有可能的等位基因分型列在csv格式的文件中，传入上述“senator”函数中，产生焦磷酸测序反应所需的碱基分配顺序：

SEQ ID No.36:gtgcatcagctagctcagatacgacgtcgcgtgtatgcgcagtcagca。为最大限度地区分不同微单倍型等位基因(实际单倍型)分型结果，使得不同微单倍型等位基因(实际单倍型)所生成的焦磷酸测序图谱各不相同，因此需要综合考虑不同微单倍型等位基因(实际单倍型)分型的特有序列(unique nucleotide sequence，UNS)以生成一段碱基分配顺序使得不同微单倍型等位基因(实际单倍型)所生成的焦磷酸测序图谱的差异性最大化。

(7)使用血液/细胞/组织基因组提取试剂盒(天根生化科技(北京)有限公司，中国)，按照试剂盒配套的标准操作流程对样本(健康个体的外周静脉血)进行DNA模板提取。步骤(6)是产生焦磷酸测序反应所需的碱基分配顺序。

(8)利用荧光定量试剂盒Investigator Quantiplex Kit(QIAGEN，德国)按照试剂盒配套的标准操作流程对所提取的DNA进行实时荧光定量PCR。

(9)对单一DNA的单一微单倍型以及单一DNA和混合DNA的复合微单倍型进行实时荧光定量PCR。

(10)根据说明书(

Q96 ID，QIAGEN，德国)对PCR产物进行测序前模板的处理以及焦磷酸测序。

(11)使用基于稀疏表示理论的AdvISER-M-PYRO算法对单一DNA复合微单倍型进行解析，软件输出的分析结果以贡献系数(contribution coefficients)及相关系数(correlation coefficients,r)呈现。拥有贡献系数的项为可能的微单倍型。使用该算法对混合DNA贡献者构成进行解析，贡献系数最高者为混合DNA中最有可能的贡献者来源。

(12)使用支持向量机(SVM)算法对单一DNA模板的复合微单倍型进行分型，并对混合DNA模板的贡献者构成进行分类。

(1)对单一微单倍型进行焦磷酸测序以获得其基因型分型结果

单一微单倍型测序模板制备的扩增体系组成如表2所示，其PCR的反应条件如表3所示。

表2单一微单倍型焦磷酸测序模板制备的反应体系

表3焦磷酸测序模板制备的反应条件

四个微单倍型基因型分型的焦磷酸测序图谱见图1。四个微单倍型的实际单倍型(等位基因)分型结果，相应微单倍型在CHB群体中的A_e值，杂合度以及个人识别概率见表4。

表4各微单倍型序列分型结果

(2)计算微单倍型的混合DNA检测能力(Probabilities of Detecting aMixture,PDM)

本方法则能够根据已有的微单倍型等位基因频率，准确地计算出各位点混合DNA检测能力(PDM)，并且可以评估位点联合使用时的累计概率。

本方法的基本原理如下：

微单倍型基因座上有t个等位基因(x1，x2，...，xt)，基因频率依次为(p1，p2...，pt)，则p1+p2+...Pt＝1。

混合斑中有来自m个独立个体的DNA，混合斑样本在该基因座上存在n个等位基因，对于常染色体遗传标记来说n＝2m。

根据孟德尔遗传定律，混合斑中的等位基因组合是m个随机无关个体的Hardy-Weinberg基因型阵列的乘积，函数为

(x1+x2+...xt)×(x1+x2+...xt)×...×(x1+x2+...xt)＝(x1+x2+...+xt)ⁿ

根据多项式定理证明，

展开多项式后，可以计算混合斑中该基因座等位基因组合的构成。展开式中共有tⁿ项，合并同类项之后，

的系数为

各等位基因组合出现的概率为各等位基因频率累乘与系数之积。

例如：该基因座在群体中的等位基因(频率)为p(P1)，q(P2)，r(P3)，2个无关个体混合斑中同时检出上述三种等位基因组合的系数为

等位基因组合概率的计算见下表15：

表15等位基因组合概率的计算

m个独立个体的遗传物质混合后，一基因座上检出n个或n-1个等位基因时才能判定为m个独立个体遗传物质混合。故含有n个和n-1个等位基因组合的概率累加得到该位点对m个无关个体混合的识别概率(P)。

彼此满足自由组合定律的y个基因座联合时，其累积的识别概率为

经过上述方法的计算，四个微单倍型对于由2，3和4个贡献者构成的混合DNA的累计检出能力见表。

表5微单倍型的混合DNA检测能力

(3)单一DNA和混合DNA的复合微单倍型测序模板制备

单一DNA和混合DNA的复合微单倍型测序模板制备的扩增体系组成如表6所示，反应条件如表3所示。混合DNA的各贡献者成分如表7。根据实验结果，单一DNA的复合微单倍型焦磷酸测序体系灵敏度可达40pg/μL，因此混合DNA中次要贡献者成分在体系中的最低浓度为40pg/μL。

表6复合微单倍型焦磷酸测序模板制备的反应体系

表7混合DNA中各贡献者成分构成

(4)测序前模板预处理及焦磷酸测序

根据说明书(

Q96 ID，QIAGEN，德国)进行前期测序模板的处理，以获得焦磷酸测序所需的单链DNA模板，并进行焦磷酸测序。

(5)使用AdvISER-M-PYRO算法解析单一DNA复合微单倍型分型以及解析混合DNA贡献者构成。

将四个微单倍型位点所测得的17个微单倍型基因型作为字典，以61，76，79号DNA样本的四个微单倍型复合微单倍型焦磷酸测序结果为测试信号。这三个个体在四个微单倍型位点上的等位基因(实际单倍型)见表8，其中每个微单倍型等位基因(实际单倍型)序列后标记有表4中对应基因型的编号。AdvISER-M-PYRO算法的解析结果如表9所示。图2显示三个DNA的复合微单倍型焦磷酸测序图谱，其中(a)为61号DNA，(b)为76号DNA，(c)为79号DNA。表9列出了三个DNA样本复合焦磷酸测序图谱的解析结果，由于微单倍型的等位基因(实际单倍型)序列较长，于是用表4中的相应字母编号表示。其中：61号样本正确解析出mh06ZJ-001和mh19ZJ-001两种微单倍型基因型，分别对应表4中相应编号；76号样本正确解析出mh03ZJ-001、mh06ZJ-001和mh07ZJ-001三种微单倍型基因型，分别对应表4中相应编号；79号样本正确解析出mh03ZJ-001和mh06ZJ-001两种微单倍型基因型，分别对应表4中相应编号。以上微单倍型中存在解析错误的基因型用下划线字体表示。

表8 61，76，79号DNA样本的四个微单倍型等位基因(实际单倍型)

表9对单一DNA复合微单倍型测序结果的贡献系数与相关系数

构建了三组由两个体等比组成的DNA模拟混合样本，其中三组混合DNA的测序结果均能够解析正确(图3，表10)。

表10对两个体等比混合DNA复合微单倍型测序结果解析的贡献系数与相关系数

使用33和75号个体构建了两个体不等比混合的样本，混合比分别为1:5,1:10,1:20,1:30(33号个体为次要成分，75号个体为主要成分)，以上每个比例的混合样本所对应的编号分别为：M1:5，M1:10，M1:20和M1:30。这四组不同混合比例的样本在由6个个体样本作为字典的条件下，可以正确解析出各自的贡献个体(图4，表11)。经过对表11中的数据进行回归分析，发现两贡献者的贡献系数比与混合比之间存在良好的相关关系,回归直线对观测值具有较高的拟合程度(见图5，R²＝0.9864)。

表11对两个体不等比混合DNA复合微单倍型测序结果解析的贡献系数与相关系数

为了检测该体系对三个个体组成的混合DNA样本的检测能力，我们构建了两组由三个个体等比混合的DNA混合样本，分别编号为M33-75-61和M76-61-65。解析结果见表12，测序图谱见图6。

表12对三个体等比混合DNA复合微单倍型测序结果解析的贡献系数与相关系数

为了进一步检测该体系对于四个个体组成混合DNA样本焦磷酸测序图谱的解析情况，我们构建了两组由四个个体等比混合的DNA混合样本，分别编号为M33-61-65-75和M33-61-65-76。字典构成与之前一样，还是相同的6个个体。解析结果见表13，其中编号为M33-61-65-75的混合样本包含33，61，65和75号个体，测序图谱见图7(a)，该混合样本会出现75号样本无法检出的错误解析结果(贡献系数为0，用红色字体表示)；编号为M33-61-65-76的混合样本包含33,61,65和76号个体，测序图谱见图7(b)，该混合样本会出现76号样本无法检出的错误解析结果(贡献系数为0，用下划线字体表示)。

表13对四个体等比混合DNA复合微单倍型测序结果解析的贡献系数与相关系数

(5)使用支持向量机(SVM)算法解析单一DNA复合微单倍型分型以及解析混合DNA贡献者构成。

在开始使用支持向量机模型进行数据训练之前，需要根据已经获得的焦磷酸测序图谱的峰高值来模拟混合DNA和单一DNA的复合微单倍型焦磷酸测序结果。

首先对单一及混合DNA的复合微单倍型进行焦磷酸测序(设置二至四个平行复孔)，将碱基分配顺序相应位置上的峰高进行均值和标准差的计算，然后导入Python3.6中的一个扩展程序库-Numpy(Numerical Python)。Numpy支持大量的维度数组与矩阵运算，此外也针对数组运算提供大量的数学函数库。调用Numpy中的random.normal函数生成高斯分布的概率密度随机数，以得到指定数量的模拟峰高。在获得充足数量的模拟焦磷酸测序结果后，将其作为数据集，导入Python 3.6中的Scikit learn(sklearn)模块，调用其中的sklearn.model_selection.train_test_split模块对数据集进行划分，其中参数设置为test_size＝0.3，即数据集中有70％的数据用于训练，其余30％用于测试模型性能。

在模型训练阶段，首先调用LinearSVC函数构建二分类器，其中参数设置为random_state＝0,max_iter＝1000000。然后调用OneVsRestClassifier函数根据二分类器构建多分类器。最后调用clf.fit函数和clf.predict函数分别进行模型训练和样本预测。最终模型的性能以F1值(F1 score)，精准率(Precision)和召回率(Recall)这三个指标进行评估。

对33、61、65、75、76和79号DNA分别进行了单一DNA的复合微单倍型焦磷酸测序。其中33、61和75号DNA有四个平行复孔的测序结果，65和76号DNA有三个平行复孔的测序结果，79号DNA有两个平行复孔的测序结果。将每个DNA的测序结果分别在相应碱基分配顺序上对其峰高进行均值和方差的计算，通过上述方法各获得1000个模拟测序峰高结果，这些数据作为x导入训练模型。每一个DNA的测序结果对应一个标签(y)，如33号DNA对应1，61号DNA对应2，以此类推。这些标签作为y导入训练模型。

在对单一DNA的复合微单倍型组合进行分类时，测试集中出现了一个分类错误(表14)，即将本该分型组合为76号个体的样本错误地分类为65号个体。总的测试样本量为1800，因此精准率为1799/1800≈0.999444。

设置了三组混合DNA样本，第一组是33和75号等比混合；第二组是61和79号等比混合；第三组是33，61和75号等比混合。对这三组混合DNA分别进行了复合微单倍型焦磷酸测序。其中第一组混合DNA有两个平行复孔的测序结果，第二组混合DNA有四个平行复孔的测序结果，第三组混合DNA有三个平行复孔的测序结果。将每组混合DNA的测序结果分别在相应碱基分配顺序上对其峰高进行均值和标准差的计算，通过上述方法获得1000个模拟测序峰高结果，这些数据作为x导入训练模型。每一组混合DNA的测序结果对应一个标签(y)，如第一组混合DNA对应1，第二组混合DNA对应2，以此类推。这些标签作为y导入训练模型。

在对混合DNA的贡献者构成进行分类时，测试集全部得到正确分类(表14)。

表14支持向量机算法对于单一/混合DNA焦磷酸测序结果的解析

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

序列表

<110> 四川大学

<120> 一种基于复合微单倍型焦磷酸测序图谱解析的法医学混合DNA的分析方法

<160> 36

<170> SIPOSequenceListing 1.0

<210> 1

<211> 23

<212> DNA

<213> 人工序列(Artificial Sequence)

<400> 1

cctgataaac gacggtttct tga 23

<210> 2

<211> 19

<212> DNA

<213> 人工序列(Artificial Sequence)

<400> 2

tgggtggggt catcctgat 19

<210> 3

<211> 21

<212> DNA

<213> 人工序列(Artificial Sequence)

<400> 3

gctgcttctg gtcaaaactg g 21

<210> 4

<211> 24

<212> DNA

<213> 人工序列(Artificial Sequence)

<400> 4

accaagatct taagcctccc aaag 24

<210> 5

<211> 24

<212> DNA

<213> 人工序列(Artificial Sequence)

<400> 5

aattattcca ttcccagggt agtg 24

<210> 6

<211> 22

<212> DNA

<213> 人工序列(Artificial Sequence)

<400> 6

tcacattcca ttctcccaac tg 22

<210> 7

<211> 23

<212> DNA

<213> 人工序列(Artificial Sequence)

<400> 7

ttttggaagg gatgggtgga tag 23

<210> 8

<211> 20

<212> DNA

<213> 人工序列(Artificial Sequence)

<400> 8

agtctcaggg tcccagtcat 20

<210> 9

<211> 15

<212> DNA

<213> 人工序列(Artificial Sequence)

<400> 9

tgctttgcgc cctcc 15

<210> 10

<211> 19

<212> DNA

<213> 人工序列(Artificial Sequence)

<400> 10

tcaaaactgg ctttaacaa 19

<210> 11

<211> 18

<212> DNA

<213> 人工序列(Artificial Sequence)

<400> 11

agagactcct ttaaagcg 18

<210> 12

<211> 15

<212> DNA

<213> 人工序列(Artificial Sequence)

<400> 12

cctgcaacag ccctg 15

<210> 13

<211> 28

<212> DNA

<213> 人工序列(Artificial Sequence)

<400> 13

ctgcacagaa agggggctgt gcacgccc 28

<210> 14

<211> 28

<212> DNA

<213> 人工序列(Artificial Sequence)

<400> 14

ttgccaacgc aggctgtgca aaccctag 28

<210> 15

<211> 28

<212> DNA

<213> 人工序列(Artificial Sequence)

<400> 15

ctgcacagaa agggggctgt gcaaaccc 28

<210> 16

<211> 28

<212> DNA

<213> 人工序列(Artificial Sequence)

<400> 16

ttgccaacgc aggctgtgca cgccctag 28

<210> 17

<211> 29

<212> DNA

<213> 人工序列(Artificial Sequence)

<400> 17

ggctttacat cccttagact ctatggatg 29

<210> 18

<211> 29

<212> DNA

<213> 人工序列(Artificial Sequence)

<400> 18

gcctttacgt cccttagact ctatggatg 29

<210> 19

<211> 25

<212> DNA

<213> 人工序列(Artificial Sequence)

<400> 19

ggctttactg tcccttgact ctaga 25

<210> 20

<211> 25

<212> DNA

<213> 人工序列(Artificial Sequence)

<400> 20

gcttttactg tcccttgact ctaga 25

<210> 21

<211> 27

<212> DNA

<213> 人工序列(Artificial Sequence)

<400> 21

ggctttatgt ccctttgact ctacaga 27

<210> 22

<211> 27

<212> DNA

<213> 人工序列(Artificial Sequence)

<400> 22

ggctttacat cccttacact ctataac 27

<210> 23

<211> 21

<212> DNA

<213> 人工序列(Artificial Sequence)

<400> 23

gcctttacgt ctgttagggt g 21

<210> 24

<211> 28

<212> DNA

<213> 人工序列(Artificial Sequence)

<400> 24

ggccttacat cccttagact ctatgaga 28

<210> 25

<211> 28

<212> DNA

<213> 人工序列(Artificial Sequence)

<400> 25

ggctttatgt ccctttgact ctatgatg 28

<210> 26

<211> 29

<212> DNA

<213> 人工序列(Artificial Sequence)

<400> 26

ggccttacat cccttagact ctatggatg 29

<210> 27

<211> 19

<212> DNA

<213> 人工序列(Artificial Sequence)

<400> 27

ggccttacat ctcttagac 19

<210> 28

<211> 20

<212> DNA

<213> 人工序列(Artificial Sequence)

<400> 28

cgtgcgtgcg tgcacacgca 20

<210> 29

<211> 21

<212> DNA

<213> 人工序列(Artificial Sequence)

<400> 29

cgtgcgtgcg tgcgcgcgca c 21

<210> 30

<211> 20

<212> DNA

<213> 人工序列(Artificial Sequence)

<400> 30

cgtgcgtgcg tgcacacaca 20

<210> 31

<211> 21

<212> DNA

<213> 人工序列(Artificial Sequence)

<400> 31

cgtgcgtgcg tgcgcgcgcg c 21

<210> 32

<211> 19

<212> DNA

<213> 人工序列(Artificial Sequence)

<400> 32

ggtcagacac tccacagtc 19

<210> 33

<211> 19

<212> DNA

<213> 人工序列(Artificial Sequence)

<400> 33

agtcagagcc gccacggtc 19

<210> 34

<211> 19

<212> DNA

<213> 人工序列(Artificial Sequence)

<400> 34

agtcagagcc gccacagtc 19

<210> 35

<211> 19

<212> DNA

<213> 人工序列(Artificial Sequence)

<400> 35

ggtcagacac tccacggtc 19

<210> 36

<211> 48

<212> DNA

<213> 人工序列(Artificial Sequence)

<400> 36

gtgcatcagc tagctcagat acgacgtcgc gtgtatgcgc agtcagca 48

Claims

1.一种基于复合微单倍型焦磷酸测序图谱解析的法医学混合DNA的分析方法，其特征在于，包括以下步骤：

1)提取样本的DNA，得到样本DNA；

所述步骤2)PCR使用的引物包括SEQ ID No.1～8所示的核苷酸序列；所述SEQ ID No.1～2扩增微单倍型mh03ZJ-001；所述SEQ ID No.3～4扩增微单倍型mh06ZJ-001；所述SEQ IDNo.5～6扩增微单倍型mh07ZJ-001；所述SEQ ID No.7～8扩增微单倍型mh19ZJ-001；

所述步骤3)焦磷酸测序使用的引物包括SEQ ID No.9～12所示的核苷酸序列；所述SEQID No.9测序微单倍型mh03ZJ-001；所述SEQ ID No.10测序微单倍型mh06ZJ-001；所述SEQID No.11测序微单倍型mh07ZJ-001；所述SEQ ID No.12测序微单倍型mh19ZJ-001。

2.根据权利要求1所述的分析方法，其特征在于，所述步骤3)焦磷酸测序，当所述PCR产物中仅含有微单倍型mh03ZJ-001时，所述焦磷酸测序的体系每50μl包括：2×GoTaqColorless Master Mix 25μl、浓度为10μM的SEQ ID No.91.2μl、无核酸酶水21.8μl和浓度为5ng/μL的PCR产物2μl。

3.根据权利要求1所述的分析方法，其特征在于，所述步骤3)焦磷酸测序，当所述PCR产物中仅含有微单倍型mh06ZJ-001时，所述焦磷酸测序的体系每50μl包括：2×GoTaqColorless Master Mix 25μl、浓度为10μM的SEQ ID No.101μl、无核酸酶水22μl和浓度为5ng/μL的PCR产物2μl。

4.根据权利要求1所述的分析方法，其特征在于，所述步骤3)焦磷酸测序，当所述PCR产物中仅含有微单倍型mh07ZJ-001时，所述焦磷酸测序的体系每50μl包括：2×GoTaqColorless Master Mix 25μl、浓度为10μM的SEQ ID No.111.4μl、无核酸酶水21.6μl和浓度为5ng/μL的PCR产物2μl。

5.根据权利要求1所述的分析方法，其特征在于，所述步骤3)焦磷酸测序，当所述PCR产物中仅含有微单倍型mh19ZJ-001时，所述焦磷酸测序的体系每50μl包括：2×GoTaqColorless Master Mix 25μl、浓度为10μM的SEQ ID No.120.8μl、无核酸酶水22.2μl和浓度为5ng/μL的PCR产物2μl。

6.根据权利要求2～5任一项所述的分析方法，其特征在于，制备焦磷酸测序使用的模板的程序为：预变性95℃2min；Touchdown阶段95℃30s，71℃45s，72℃30s，20个循环，每一个循环71℃下降0.5℃；循环扩增阶段95℃30s，61℃45s，72℃30s，25个循环；终延伸72℃5min。

7.根据权利要求1所述的分析方法，其特征在于，当所述PCR产物中同时含有微单倍型mh03ZJ-001、mh06ZJ-001、mh07ZJ-001和mh19ZJ-001时，所述焦磷酸测序的体系每50μl包括：2×GoTaq Colorless MasterMix 25μl、浓度为0.24μM的SEQ ID No.1～2各1.2μl、浓度为0.2μM的SEQ ID No.3～4各1μl、浓度为0.28μM的SEQ ID No.5～6各1.4μl、浓度为0.16μM的SEQ ID No.7～8各0.8μl、无核酸酶水14.2μl和浓度为5ng/μL的PCR产物2μl。

8.根据权利要求7所述的分析方法，其特征在于，制备焦磷酸测序使用的模板的程序为预变性95℃2min；Touchdown阶段95℃30s，71℃45s，72℃30s，20个循环，每一个循环71℃下降0.5℃；循环扩增阶段95℃30s，61℃45s，72℃30s，25个循环；终延伸72℃5min。