基于多重PCR引物组的中华穿山甲地理溯源方法及其装置

本技术涉及罚没穿山甲溯源鉴定，特别是一种基于多重pcr引物组的中华穿山甲地理溯源方法及其装置。

背景技术：

1、中华穿山甲（

manis pentadactyla）是鳞甲目，穿山甲科的哺乳动物，其地理分布范围较广，西至喜马拉雅山脉南部的尼泊尔，东至中国东部的台湾岛，并从我国长江以南延展到中南半岛北部。目前中华穿山甲已经被列为我国国家一级重点保护野生动物，被 iucn濒危物种红色名录列为极度濒危物种。

2、从源头精准打击非法贸易是保护中华穿山甲最直接有效的手段，衍生出开发溯源方法的需求。现存的中华穿山甲地理溯源方法包括：贸易记录分析、分子遗传学分析方法等，分子遗传学分析方法通过建立穿山甲遗传信息和地理区域或地理坐标之间的关系，从而根据样品特定分子标记信息预测待测样本来源地。

3、现有种群研究发现中华穿山甲存在明显群体结构，且各种群拥有相对独立的地理分布区域，为其遗传溯源研究提供核心基础。同时，筛选少量群体间高度分化的snp 标记进行溯源的研究方法已经成功运用到树穿山甲非法贸易溯源、鸟类迁徙追踪等领域，为中华穿山甲遗传溯源研究提供了重要借鉴。

4、中华穿山甲非法贸易样品鉴定由于走私案件中待鉴定样品数量庞大且贸易记录缺失或多种穿山甲甲片混杂，现有鉴定方法用于此类甲片鉴定存在以下问题：

5、（1）贸易记录分析方法局限性：该方法高度依赖贸易记录的完整性和准确性，而非法贸易中的甲片常常因为运输路线长，转运环节多，贸易信息隐秘等，难以保证贸易记录的完整性和准确性，因此贸易记录分析方法所得结果不可靠；同时贸易记录分析方法难以准确地定位非法贸易样品的来源地，往往只能指向某一大致的区域，影响甲片样本品种精准确定；贸易记录分析方法虽然简单直接，但难以保证结果的可靠性和精准。

6、（2）分子遗传学分析方法局限性：该方法包括多种分析方法，一类方法根据样品的遗传信息确定待测样品的群体，根据该样品对应群体的分布区域确定其地理来源。但该方法在检测时所需位点数量较多，而样本提取所得遗传信息位点无法做到与种群位点完全一一对应，导致地理溯源结果的精度较低，实际应用中性价比较低；

7、另一类方法，在全基因组水平筛选群体间高分化位点，需要提前进行群体结构分析，为每个样品提供群体标签，群体结构分析则又需要多种分析方法相互印证才能确定待测样品对应群体，该方法分析流程复杂，耗时长，结果准确解读过程复杂需要大量专业知识；

8、还有一类方法，使用全基因组变异位点信息和深度学习框架进行地理溯源，这种方法需要全基因组测序数据及大量计算资源，所需时间较长，在实际应用中成本较高，耗时较长。

9、分子遗传学分析方法中的多类方法，虽然可以用多种方式衡量其结果的精确度和可靠性，但在实际应用场景下效率低、成本高，两者均无法同时满足高准确率、低耗时、低成本等执法核心需求，导致基层执法部门在鉴定大规模走私甲片时，面临人力和经济的双重负担，难以实现对非法贸易的高效打击。

10、公开于背景技术部分的信息仅仅旨在增加对本发明的总体背景的理解，而不应当被视为承认或以任何形式暗示该信息构成已为本领域普通技术人员所公知的现有技术。

技术实现思路

1、本技术针对上述技术问题提供了一种基于多重pcr引物组的中华穿山甲地理溯源方法及其装置。

2、本技术提供了一种基于多重pcr引物组的中华穿山甲地理溯源方法，包括以下步骤：

3、步骤s1：获取参考样品的全基因组测序数据和来源地地理坐标信息，并将各样品的全基因组测序数据与来源地地理坐标信息一一对应后存储为参考样品数据集；

4、步骤s2：将参考样品数据集中的各参考样品测序数据比对到已发表的组装到染色体水平的中华穿山甲参考基因组上，使用picard去除测序过程中由于pcr产生的冗余序列，利用genome analysis toolkit v4.0获取snp数据集，利用haplotypecaller模块鉴定每个样品对应个体的snp，将每个个体的snp合并到一个文件中形 成 vcf 文 件，用variantfiltration模块中的vcftools去除次要等位基因频率小于0.05 的 snps，去除群体水平上同一位置缺失位点覆盖大于20% 样本的snps，保证群体中等位基因数为2，去除测序深度分布在所有位点中低于2.5%或高于97.5%的 snps，完成对获得的snp的过滤，去除假阳性位点和多等位基因位点，得到参考样品snps数据集；

5、步骤s3：基于king分析对参考样品snps数据集中各样品进行亲缘关系分析，得到多组相同个体组和单个体数据集，仅保留相同个体组中有具体地理来源或测序质量较好的个体信息，去除相同个体组中冗余数据后，得到清洗数据集，过滤清洗数据集，得到的高质量snps数据集；

6、步骤s4：用eigensoft v.4.2 软件包中的 smartpca 程序对高质量snps数据集进行pca分析，获得每个样品的主成分值、各主成分值对应的方差值，以及各snp对主成分的贡献度值，通过tracy-widom 检验，确定对高质量snps数据集中数据解释度大且真实的主成分，根据贡献度的大小降序排列高质量snps数据集中各snp位点，选取snp位点，合并所选snp位点作为地理溯源snp分子标记数据集选取snp位点为：根据主成分方差比值确定每个主成分里选择位点的比例，然后根据要选择的位点总数乘对应的比例得到每个主成分里要选前多少个位点，根据所得待选位点数量在每个主成分下选取对应数量的贡献前列snp位点；

7、步骤s5：取数个梯度的小规模位点组合，分别将参考样品中这些位点组合的基因型信息和各参考样品的地理来源地经纬度坐标输入地理溯源软件origen，使用留一法交叉验证评估每一组位点用于地理溯源的误差距离情况，即每一轮将一个参考样品的地理来源地经纬度坐标设为未知，经过地理溯源软件origen预测得到预测坐标，使用r包geodist计算预测坐标与真实坐标之间的距离，即地理溯源的误差距离，统计每组所有样品的地理溯源误差距离，比较每组的地理溯源误差距离平均数以及误差距离大于500km的个体数，将以上2个指标中至少一个为最小对应的地理溯源snp分子标记数据集中各组位点，作为最佳snp分子标记集；

8、步骤s6：使用primer3分别生成目标位点的候选引物对，再使用mfeprimer-3.0作为核心质控模块，经多轮迭代筛选与全局检查，筛选获得在同一反应体系中无相互干扰、可高效特异性扩增的多重 pcr 引物池；

9、步骤s7：提取未知来源待鉴定样品的总dna，使用定制的多重pcr引物组对该样品进行靶向捕获并构建二代测序文库，使用高通量二代测序平台对文库进行测序得到二代测序数据集，对二代测序数据集进行基因分型数据处理得到该待测样品最佳snp分子标记集的基因型信息；

10、步骤s8：采用r包 origen进行溯源分析，将s7中获得的最佳snp分子标记集的基因型信息输入origen软件中获得待鉴定样品来源地理位置的预测经度和预测纬度信息，以及每一个样品的预测后验概率，作为预测结果输出。

11、优选地，参考样品的全基因组测序数据和来源地地理坐标信息来源：已发表文献、新收集的样品。

12、优选地，相同个体组按以下步骤筛选并过滤：

13、在plink v.2.0中使用“make-king”命令，遍历计算参考样品snps数据集任意两个体间的king-robust coefficients，将king-robust coefficients大于 0.354的任意两个体全基因序列构建为一个相同个体组，每组相同个体组中取测序覆盖度最高的一个个体进行后续分析，其余相同个体从数据集中过滤。

14、优选地，步骤s5还包括：位点组合精简操作；

15、位点组合精简操作包括以下步骤：

16、步骤s511：使用软件hitsnp的祖先估计功能，基于纯种样本模拟杂交群体构建训练数据，选用 classifier i 中的 lr-i、svm-i、rf-i 输出的预测概率训练 classifierii；

17、步骤s512：通过 classifier ii 区分最佳snp分子标记集中样本为纯种或杂交种并推断祖先数量，再依据 lr-i 预测概率大于 0.05 的阈值锁定杂交样本的具体祖先品种；

18、步骤s513：使用最小特征 snp 子集筛选功能，以特征snp 为初始集，采用序列前向选择算法结合 bagnb-iii 分类器，按数量梯度迭代筛选，以分类准确率为核心评估标准，获取可精准区分纯种群体的最小特征 snp 子集作为最佳snp分子标记集用于步骤s6中。

19、优选地，步骤s7中构建二代测序文库过程包括：两轮pcr扩增与纯化，之后进行文库片段长度分布的检测，得到符合所有质控标准的二代测序文库。

20、优选地，步骤s7中的二代测序数据集进行数据处理包括以下步骤：

21、步骤s721：通过 fastqc 对获取的原始测序数据进行数据质量评估，使用trimmomatic去除接头污染序列、低质量碱基及长度过短的读段，得到高质量测序数据；

22、步骤s722：利用 burrows-wheeler aligner将高质量测序数据比对至中华穿山甲参考基因组，通过 samtools进行比对结果的排序、索引与格式转换，利用genome analysistoolkit v4.0获得snp数据，最后使用vcftools获得该待测样品最佳snp分子标记集的基因型信息。

23、本技术的另一方面还提供了一种用于运行如上述的方法的装置，包括：

24、样本信息收集模块，用于获取参考样品的全基因组测序数据和来源地地理坐标信息，并将各参考样品的全基因组测序数据与来源地地理坐标信息一一对应后存储为参考样品数据集；

25、变异位点检测过滤模块，用于将参考样品数据集中的各参考样品测序数据比对到已发表的组装到染色体水平的中华穿山甲参考基因组上，使用picard去除测序过程中由于pcr产生的冗余序列，利用genome analysis toolkit v4.0获取snp数据集，利用haplotypecaller模块鉴定每个样品对应个体的snp，将每个个体的snp合并到一个文件中形 成 vcf 文 件，用 variantfiltration模块中的vcftools去除次要等位基因频率小于0.05 的 snps，去除群体水平上同一位置缺失位点覆盖大于20% 样本的snps，保证群体中等位基因数为2，去除测序深度分布在所有位点中低于2.5%或高于97.5%的 snps，完成对获得的snp的过滤，去除假阳性位点和多等位基因位点，得到参考样品snps数据集；

26、亲缘关系鉴定模块，用于基于king分析对参考样品snps数据集中各样品进行亲缘关系分析，得到多组相同个体组和单个体数据集，仅保留相同个体组中有具体地理来源或测序覆盖度较好的个体信息，去除相同个体组中冗余数据后，得到清洗数据集，过滤清洗数据集，得到的高质量snps数据集；

27、pca模块，用于用eigensoft v.4.2 软件包中的 smartpca 程序对高质量snps数据集进行pca分析，获得每个样品的主成分值、各主成分值对应的方差值，以及各snp对主成分的贡献度值，通过tracy-widom 检验，确定对高质量snps数据集中数据解释度大且真实的主成分，根据贡献度的大小降序排列高质量snps数据集中各snp位点，选取snp位点，合并所选snp位点作为地理溯源snp分子标记数据集，逐一对比地理溯源snp分子标记数据集中各组位点，相对地理来源数据集中各位点的误差距离平均数以及误差距离大于500km的样品数，然后选取误差距离平均数以及误差距离大于500km的样品数最小的地理溯源snp分子标记数据集中各组位点作为最佳位点集；选取snp位点为：根据主成分方差比值确定每个主成分里选择位点的比例，然后根据要选择的位点总数乘对应的比例得到每个主成分里要选前多少个位点，根据所得要选位点数量选取对应数量的snp位点；

28、获取最佳snp分子标记集模块，取数个梯度的小规模位点组合，分别将参考样品中这些位点组合的基因型信息和各参考样品的地理来源地经纬度坐标输入地理溯源软件origen，使用留一法交叉验证评估每一组位点用于地理溯源的误差距离情况，即每一轮将一个参考样品的地理来源地经纬度坐标设为未知，经过地理溯源软件origen预测得到预测坐标，使用r包geodist计算预测坐标与真实坐标之间的距离，即地理溯源的误差距离，统计不同位点组合所有样品的地理溯源误差距离，比较每组的地理溯源误差距离平均数以及误差距离大于500km的个体数，以上2个指标中至少一个为最小对应的地理溯源snp分子标记数据集中各组位点，作为最佳snp分子标记集；

29、多重pcr引物池获取模块，用于使用primer3分别生成目标位点的候选引物对，再使用mfeprimer-3.0作为核心质控模块，经多轮迭代筛选与全局检查，筛选获得在同一反应体系中无相互干扰、可高效特异性扩增的多重 pcr 引物池；

30、高质量测序数据获取和基因分型模块，用于提取未知来源待鉴定样品的总dna的提取，使用定制的多重pcr引物组对质控的样品进行靶向捕获并构建未知待鉴定样品的二代测序文库，使用高通量测序平台对二代测序文库进行测序得到二代测序数据集，对二代测序数据集进行数据处理得到高质量测序数据，通过基因分型数据处理得到该待测样品最佳snp分子标记集的基因型信息；

31、预测结果输出模块，用于采用r包 origen进行溯源分析，将待鉴定样品的最佳snp分子标记集的基因型信息输入origen软件中获得待鉴定样品来源地理位置的预测经度和预测纬度信息，以及每一个样品的预测后验概率，作为预测结果输出；

32、参考样品的目标snp分子标记集的基因型信息获取，使用vcftools从参考样品的全基因组snp数据集中提取得到最佳snp分子标记集的基因型信息。

33、本技术能产生的有益效果包括：

34、1）实现了成本低、操作简单的中华穿山甲地理溯源：

35、该方法无需全基因组测序中华穿山甲样品的基因组信息，只需通过多重pcr扩增目标片段，即可获得所需位点的基因型信息，有效降低检测成本；同时，该方法由于所需位点少，在变异检测和地理溯源过程中，操作较为简单，鉴定时所需算力较小，消耗时间短；

36、2）位点筛选简单

37、目前中华穿山甲地理溯源方法中，筛选位点的方法一般较为复杂。比如，一些方法需要先通过多种方法解析中华穿山甲的群体结构，为每一个样品提供一个群体标签后才能进行后续的位点筛选。本方法不需要群体结构分析这样的前置分析，得到snp数据集后就可以直接进行目标位点的筛选。

38、3）溯源精度高

39、该方法得到的最佳位点集的地理溯源误差距离平均数以及误差距离大于500km个体数分别为163.92km，6，而全基因组snps使用深度学习框架locator的溯源方法的误差距离平均数为175.43 km，误差距离大于500km的个体数为9。综合三个指标，两者相比较，针对中华穿山甲地理溯源筛选到的小规模snps的地理溯源效果接近全基因组snps的地理溯源效果。