CN109207606A

CN109207606A - 用于亲权鉴定的ssr位点的筛选方法和应用

Info

Publication number: CN109207606A
Application number: CN201811125248.5A
Authority: CN
Inventors: 兰道亮; 陈通; 吉文汇; 李键; 朱育星; 蔡雯祎; 熊显荣; 张大伟
Original assignee: Southwest Minzu University
Current assignee: Southwest Minzu University
Priority date: 2018-09-26
Filing date: 2018-09-26
Publication date: 2019-01-15
Anticipated expiration: 2038-09-26
Also published as: CN109207606B

Abstract

本发明提供了一种用于亲权鉴定的SSR位点的筛选方法和应用，涉及生物技术领域。该筛选方法首先从待筛选SSR位点的物种的全基因组序列中找出所有的潜在SSR位点，然后经过依次的筛选和过滤筛选出候选的SSR位点；然后使用模拟筛选的方法，依据待筛选SSR位点的物种的遗传特征，构建模拟亲代并杂交产生子代，以验证能够用于亲权鉴定的SSR位点的特征；然后将候选SSR位点按照PE值排序，然后筛选出符合如下条件的SSR位点：(b1)MAF值≤M；(b2)缺失值比例＜10％；(b3)任意两个SSR位点的距离＞LD衰减距离；(b4)重复单元≥N；满足上述条件的SSR位点可用于亲权鉴定。

Description

用于亲权鉴定的SSR位点的筛选方法和应用

技术领域

本发明涉生物技术领域，尤其是涉及一种用于亲权鉴定的SSR位点的筛选方法和应用。

背景技术

亲权鉴定亦称亲子鉴定，是将生物学、分子遗传学、医学方法结合起来，根据亲本与后代的形态结构和遗传物质进行遗传相似性分析，确定亲代与子代的关系。亲权鉴定的基本原则是应用遗传规则，即后代的遗传标记必须是父母各半。亲权鉴定的遗传学基础是孟德尔分离法则和自由组合法则，法医鉴定最先进行了应用，现在逐渐被应用于遗传育种及保护濒临灭绝的动物。

在现代畜牧业发展过程中，遗传谱系已经得到了广泛的关注。错误的系谱信息会导致育种偏差进而导致巨大的经济损失。实际生产中会因各种原因造成系谱记录不详，导致育种选择无法确定优良种公畜，动物亲权鉴定就显得十分重要。亲权鉴定不但在家畜育种具有重要意义，还适用于野生动物繁殖和濒危动物的保护。亲权鉴定是研究动物生活史和估计个体近亲繁殖水平的有力工具。在保护管理中，统计有效的群体、减少近亲繁殖的水平也很重要。在动物具有未知谱系的情况下，还可以使用分子标记来构建群体结构。

分子标记在动物育种中变得越来越重要，并且已经成功用于物种的识别，亲权鉴定以及建立群体之间的关系等。分子标记能够直接反映物种之间核苷酸的不同，多态性高、准确度高、具有显性或共显性标记。其主要分为DNA指纹标记、微卫星DNA和单核苷酸多态性(SNP)。其中，SSR 标记具有高度多态性、遗传信息丰富、在整个基因组中分布广泛、灵敏度高、分型结果可靠、片段大小合适(一般小于400bp)，适于PCR扩增等优点。因此，一种用于亲权鉴定的SSR位点的筛选方法是目前需要的。

有鉴于此，特提出本发明。

发明内容

本发明的第一目的在于提供一种用于亲权鉴定的SSR位点的筛选方法，该方法可以有效的筛选出用于亲权鉴定的SSR位点。

本发明的第二目的在于提供上述用于亲权鉴定的SSR位点的筛选方法的应用。

为解决上述技术问题，本发明特采用如下技术方案：

一种用于亲权鉴定的SSR位点的筛选方法，包括：

(S1)筛选出待筛选SSR位点的物种的全基因组序列所有潜在的SSR 位点，再筛选出位于基因区间内的SSR位点，再过滤掉单核苷酸重复类型的SSR位点；然后筛选出具有多态性的SSR位点；

(S2)然后过滤掉满足如下条件的SSR位点，得到候选SSR位点：(a1) 缺失值比例＞50％；(a2)哈温平衡检验显著性＜0.01；

(S3)计算候选SSR位点的LD衰减距离；

(S4)将候选SSR位点按照PE值排序，然后筛选出符合如下条件的 SSR位点：(b1)MAF值≤M；(b2)缺失值比例＜10％；(b3)任意两个 SSR位点的距离＞LD衰减距离；(b4)重复单元≥N；

其中，所述M和所述N按照如下方法得到：

模拟产生具有SSR位点的亲代个体，所述亲代个体满足如下条件：(c1) 亲代个体的SSR位点距离≥LD衰减距离；(c2)亲代个体的基因型满足哈温平衡；

然后将亲代个体作为父本，按照如下条件产生子代个体：(d1)杂交过程完全随机杂交；(d2)重组过程无连锁现象；(d3)个体不产生随机突变；

然后模拟计算如下指标：

(e1)不同MAF的条件下亲子关系准确率，得到亲子关系准确率至少为99.9％时的MAF值，记为M；

(e2)不同SSR分型类型、不同位点规模下CPE的变化，以确定SSR 位点的重复单元数，记为N。

优选地，以待筛选SSR位点的物种的高通量测序数据作为全基因组序列的数据来源。

优选地，从NCBI SRA数据库中筛选所述物种的相同测序平台下，并且数据质量整齐度一致的测序样本，然后进行数据校准和整理以及重复序列标记和合并；

优选地，所述数据校准和整理包括：收集到的SRA数据首先用sra-tool kit的“fastq-dump-split-3”命令将其转换为FASTQ数据格式，用cutadapt 程序去除序列分数＜20的低质量测序片段/残基，得到过滤后测序数据。再用Burrows-Wheeler Alignment(BWA)tool对比工具的“bwa mem-T 30-h 5”程序将过滤后测序数据比对到牦牛基因组上，同时检查比对质量，把质量分值低于30的过滤掉以减少假阳性。然后用Samtools软件将比对结果进行格式转换及排序，产生对比后BAM数据文件格式；

优选地，所述重复序列标记和合并包括：用Picard Mark Duplicates软件来鉴定和标记每个比对后产生的BAM数据文件，将过度重复测到的冗余BAM数据文件进行合并，以保证每个样本产生一个对应的单独BAM数据文件。

优选地，使用scikit-learn package模拟产生具有SSR位点的亲代个体。

优选地，使用VCFtools对候选SSR位点数据进行LD分析，计算位点间的LD水平，并绘制LD衰减图。

优选地，所述筛选方法还包括分析潜在的SSR位点的特征，以用于作为模拟筛选的SSR位点库；

所述特征包括采用PSR Modify version统计具有多态性的SSR位点，及未检测到多态性的位点；统计SSR重复类型及其重复性占比；统计群体中每个个体纯合位点及杂合位点的数量；和，统计位点分型及MAF值。

优选地，所述筛选方法还包括分析待筛选SSR位点的物种和所述物种的候选SSR位点的群体遗传特征，以用于作为模拟筛选的SSR位点库；

所述群体遗传特征包括：群体结构与亲缘关系、SSRs所在区间的多态性水平和统计SSRs在基因各区域的数目。

优选地，所述待筛选SSR位点的物种为牦牛。

优选地，将牦牛的候选SSR位点按照PE值排序，然后筛选出符合如下条件的SSR位点：(b1)MAF值≤0.3；(b2)缺失值比例＜10％；(b3) 任意两个SSR位点的距离＞2kb；(b4)重复单元≥3。

本发明还提供了上述SSR位点的筛选方法在如下(x1)-(x4)中的应用：(x1)个体识别；(x2)家系管理；(x3)种质资源鉴定；(x4)遗传多态性位点分析。

与现有技术相比，本发明具有如下有益效果：

本发明提供的用于亲权鉴定的SSR位点的筛选方法，可以直接从待筛选SSR位点的物种的全基因组序列中筛选出用于亲权鉴定SSR位点，即使待筛选SSR位点的物种在公共数据库或己发表的相关文献中都没有相关的 SSR位点的报道，也可以高效的筛选出能够用于亲权鉴定的SSR位点。因此该方法具有适用范围广泛，筛选效率高，筛选出的SSR位点准确率高的优点。

本发明提供的上述用于亲权鉴定的SSR位点的筛选方法适用范围广泛，例如可以应用于个体识别，家系管理，种质资源鉴定，和遗传多态性位点分析。

附图说明

为了更清楚地说明本发明具体实施方式或现有技术中的技术方案，下面将对具体实施方式或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施方式，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例提供的具有多态性的SSR位点数及未检测到多态性的位点数；

图2为本发明实施例提供的SSR重复类型及各重复类型的占比；

图3A为本发明实施例提供的样本中每个个体纯合位点数量和杂合位点数量；

图3B为本发明实施例提供的样本中每个个体纯合位点数量与杂合位点数量的比值；

图4为本发明实施例提供的SSR位点在不同分型下的数目及其在不同 MAF值下的占比；

图5为本发明实施例提供的牦牛基因LD衰减图；

图6为本发明实施例提供的牦牛75个样本群体结构见图；

图7为本发明实施例提供的牦牛A组和B组群内和群间遗传距离；

图8为本发明实施例提供的牦牛SSR位点在基因Intergenic、CDS、 Intron、Downstream、Upstream区间上的多态性水平；

图9为本发明实施例提供的牦牛SSR位点在不同区间类型的数目；

图10为本发明实施例提供的牦牛SSR位点在区间长度矫正后不同区间类型的数目；

图11为本发明实施例提供的牦牛SSR位点模拟法筛选目标位点的预实验中相似度水平的模拟计算结果；

图12A为本发明实施例提供的牦牛SSR位点模拟法筛选目标位点中 MAF为0.1时的相似度水平的模拟计算结果；

图12B为本发明实施例提供的牦牛SSR位点模拟法筛选目标位点中 MAF为0.2时的相似度水平的模拟计算结果；

图12C为本发明实施例提供的牦牛SSR位点模拟法筛选目标位点中 MAF为0.3时的相似度水平的模拟计算结果；

图12D为本发明实施例提供的牦牛SSR位点模拟法筛选目标位点中 MAF为0.4时的相似度水平的模拟计算结果；

图12E为本发明实施例提供的牦牛SSR位点模拟法筛选目标位点中 MAF为0.5时的相似度水平的模拟计算结果；

图13为本发明实施例提供的牦牛不同分型SSR位点的数量与CPE的关系。

具体实施方式

下面将结合附图对本发明的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。实施例中未注明具体条件者，按照常规条件或制造商建议的条件进行。所用试剂或仪器未注明生产厂商者，均为可以通过市售购买获得的常规产品。

本发明提供的用于亲权鉴定的SSR位点的筛选方法，首先从待筛选 SSR位点的物种的全基因组序列中找出所有的潜在SSR位点，然后经过依次的筛选和过滤筛选出候选的SSR位点，包括：筛选出位于基因区间内的 SSR位点；过滤掉单核苷酸重复类型的SSR位点；筛选出具有多态性的SSR 位点；过滤掉缺失值比例＞50％和哈温平衡检验显著性＜0.01的SSR位点。进一步的，计算出上述候选SSR位点的LD衰减距离以进一步应用于模拟筛选。

进一步的，将候选SSR位点按照PE值排序，然后筛选出符合如下条件的SSR位点：(b1)MAF值≤M；(b2)缺失值比例＜10％；(b3)任意两个SSR位点的距离＞LD衰减距离；(b4)重复单元≥N；满足上述条件的SSR位点可用于亲权鉴定。

本发明提供的用于亲权鉴定的SSR位点的筛选方法使用模拟筛选的方法，依据待筛选SSR位点的物种的遗传特征，构建模拟亲代并杂交产生子代，以验证能够用于亲权鉴定的SSR位点的特征，包括：

然后模拟计算如下指标：

[术语解释]

微卫星(Microsatellite)，又称短串联重复序列(short Tandem Repeat， SSR)或简单重复序列(Simple Sequences Repeat，SSR)，由1-6个碱基对为核心单元进行串联重复构成。同一类微卫星DNA可分布在整个基因组不同位置，由于重复次数不用，或重复程度不完全，而形成每个座位的多态性。

最小等位基因频率(MAF)，指在给定群体中的最不常见的等位基因发生频率。

多态信息含量(Polymorphism information content，PIC)表示后代所获等位基因标记来自它的母亲或父亲的同一个等位标记的可能性，反映微卫星多态性高低的一个重要指标。多态信息含量公式如下：

I和j分别表示为第i和第j个等位基因；pi和pj分别表示第i和第j 个等位基因频率；n表示某一位点的等位基因数；N表示群体中的个体数； ii表示第i个等位基因纯合的个体数；jn表示与i共显的第n个等位基因。

杂合度(Heterozygosity，H)表示微卫星座位在群体中为杂合子的比率，主要分为期望杂合度(Expectedheterozygosity，He)和观察杂合度 (Observedheterozygosity，Ho)。Ho是指在一个群体中观察到的杂合个体总数与观察到的个体总数之比。He是在哈温平衡假设下杂合度的期望值，公式如下：

非父排除率(Probability of paternity exclution，PE)经过遗传标记检测后可以将不是生父的个体排除掉的概率叫作非父排除率，可以衡量每个遗传标记在亲权鉴定中的价值，PE的大小与被检测对象无关，与等位基因频率、等位基因数和系统遗传方式有关。

单个位点的非父排除率PE(只测父母的其中一个与子代)，计算公式如下：

pi为等位基因在群体中的频率，n为等位基因的数目。

累计非父排除率(Comulate PE，CPE)在亲权鉴定中大多使用多个遗传标记，如果各个遗传标记之间没有遗传连锁不平衡现象，累计非父排除率的公式如下。

m个位点累计非父排除率CPE为：

亲权指数(Paternity index，PI)即假设生物学父亲提供生父基因成为子代生父的可能性和随机雄性个体提供生父基因成为子代生父的可能性的比值，用来判定是否是亲生关系。

NCBI SRA数据库：为美国国立生物技术信息中心(National Center forBiotechnology Information，NCBI)推出的SRA(Sequence Read Archive) 数据库，用于存储、显示、提取和分析高通量测序数据。

哈迪-温伯格平衡定律(Hardy-Weinberg Equilibrium,HWE)：简称哈温平衡定律，由英国数学家D.H.ardy和德国医生W.Weinberg于1908年分别独立证明，是群体遗传中最重要的原理，常被用作某些遗传检验的基础。哈温平衡定律的主要内容为：在理想状况下，当基因代代传递时，群体的基因频率及基因型频率将保持不变，两者的关系也不会改变，且前者确定后者。

连锁不平衡(linkagedisequilibrium，LD)是指相邻基因座上等位基因的非随机相关，当位于某一基因座上的特定等位与同一条染色体另一基因座位上的某等位同时出现的几率高于或低于人群中的随机分布，就称这两个位点处于连锁不平衡状态。连锁不平衡检验就是测试两个观测位点的配子频率是否与基于单个等位基因频率的期望值相符。

实施例1

下面以用于牦牛亲权鉴定的SSR分子标记的筛选为实施例，进一步解释本发明。

(一)公共数据选择

公共数据选择：从NCBI SRA数据库中，筛选以牦牛genomic DNA测序为数据源的全基因组测序数据，包括全基因组测序项目和重测序项目，共计114个样本，126个Run(测序过程中的一次)。为了避免不同平台间的短序列片段识别的特异性，随机选择其中的相同平台-Hiseq2000测序平台，且测序数据质量整齐度一致的测序样本共75个。

(二)数据校准和数据整理

数据校准和整理：收集到的SRA数据首先用sra-tool kit(v2.5.2)的“fastq-dump-split-3”命令将其转换为FASTQ数据格式，用cutadapt(v.1.8.1) 程序(默认参数)去除序列分数＜20的低质量测序片段/残基，得到过滤后测序数据。再用Burrows-WheelerAlignment(BWA)tool对比工具的“bwa mem-T 30-h 5”程序将过滤后测序数据比对到牦牛基因组上，同时检查比对质量，把质量分值低于30的过滤掉以减少假阳性。然后用Samtools(V0.1.16) 将比对结果进行格式转换及排序，产生对比后BAM数据文件格式。

优选地，所述重复序列标记和合并包括：使用Picard Mark Duplicates (v1.46)(参数默认)，运行SORT_ORDER＝coordinate MERGE_SEQVENCE_DICTIONARIES＝ture命令，来鉴定每个比对后产生的BAM数据文件，并对过度重复测到的冗余BAM数据文件进行合并，以保证每个样本产生一个对应的单独BAM数据文件。

SNP、InDel的命名和过滤：为了随后的群体分析，用GATK(v3.5)中 GATK-HaplotypeCaller程序以执行GATK的最优实践来命名SNP、InDels。用vcftools(v0.1.15)软件根据以下条件进行筛选：

(1)GATK PASS(通过VCF文件中过滤柱评价为PASS标记并移除所有其他变体来维持变体)来自内部GATK；

(2)变体质量＞60；

(3)错误等级片段＞样本的50％；

(4)哈温平衡测试显著性＞0.01。

(三)筛选牦牛基因中的SSR候选位点

用PSR算法来辨别有2个残基支持的(SSRs不在残基末端)所有样本的SSR多态性位点和基因分型。牦牛SSRs的总结和模式概要由以下3个软件进行分析和绘制：pythonscripts with numpy(v1.11.1)、scipy(v0.17.0) 和matplotlib(v1.4.3)packages。在全基因组范围内公筛选出候选的SSR位点共876700个，位于基因区间内的SSR位点共194967个，对基因区间内的SSR位点，过滤单核苷酸重复类型后有69112个。

(四)收集文献报道中牦牛近缘物种(牛、羊)的亲权鉴定SSR位点，根据引物序列，筛选能够比对到牦牛基因组上的SSR位点，并与牦牛初步筛选位点进行比较。

收集文献报道中牦牛近缘物种(牛、羊)的亲权鉴定SSR位点，根据引物序列，筛选能够比对到牦牛基因组上的SSR位点，并与牦牛初步筛选位点进行比较。文献报道中总共收集了202个潜在的SSR位点，其中根据引物序列，能够比对到牦牛基因组上的序列共148个。其中147个位点能在步骤(三)筛选出的候选的SSR位点鉴定结果中找到。

(五)PSR Modify version软件分析候选SSR位点

采用PSR Modify version软件对筛选出的牦牛微卫星位点，进行如下分析：统计具有多态性的SSR位点和未检测到多态性的位点，统计SSR重复类型及各重复类型的占比，统计群体中每个个体纯合位点及杂合位点的数量，统计SSR位点分型及MAF值。

5.1多态性位点统计

用PSR Modify version识别群体的基因型，鉴定全基因组范围SSR位点的多态性，具有多态性的SSR位点及未检测到多态性的位点统计结果见图1。non-polymorphism为未检测到多态性的位点，polymorphism为具有多态性的位点，纵坐标为位点数量。从图中可以看到，具有多态性的潜在位点有56750个，未检测到多态性的位点有5965个，大部分的潜在的SSR位点具有多态性。

5.2 SSR重复类型及各重复类型的占比

SSR重复类型及其重复性占比统计结果见图2。横坐标为潜在SSR位点的单位数，纵坐标为位点数量，不同颜色代表不同SSR重复数的比例。从图中可以看到，SSR重复次数5～8次的占比最大，SSR单位数为2和3 的位点数最多，且占比远高于重复单位3和4的位点数。

其中，Mono表示SSR重复单元中包含1个碱基，Di表示SSR重复单元中包含2个碱基，Tri表示SSR重复单元中包含3个碱基，Quad表示SSR 重复单元中包含4个碱基，Penta表示SSR重复单元中包含5个碱基，Hexa 表示SSR重复单元中包含6个碱基；纵坐标为位点数量。不同颜色代表不同SSR重复单元的重复次数。

5.3统计群体中每个个体纯合位点及杂合位点的数量

75个样本中，每个个体纯合位点数量和杂合位点数量，和每个个体纯合位点数量与杂合位点数量的比值如图3A和图3B所示，从图中可以看出，纯合位点数远高于杂合位点数，纯合位点数是杂合位点数6倍。

5.4 SSR位点分型及MAF值

SSR位点分型及MAF值的统计见图4。横坐标为SSR位点具有的不同的等位等位基因数，例如SSR位点SSR_n1，具有一个等位基因，即具有的不同的等位等位基因数为1，该SSR_n1位点为二分型；例如SSR位点SSR_n2，具有两个等位基因，即具有的不同的等位等位基因数为2，该SSR_n2位点为三分型；以此类推。

纵坐标为候选SSR位点数的个数，不同颜色表示不同MAF值的占比。从图中可以看出，不同等位基因数为1和2，即二分型和三分型的位点数量最大，且不同等位基因数越大，位点数量越少。MAF值小于5％的位点数占比最大且远高于其它MAF值的位点数。

(六)分析群体遗传特征

首先，对全基因组SSR位点进行过滤，过滤掉如下2种位点：

a、在一半的样本中，都未能观察到基因型的位点，即Miss value(缺失值比例)>50％

b、不满足哈温平衡的位点，HWE(哈温平衡)检验显著性<0.01

6.1 LD衰减

用VCFtools对获得的位点数据进行LD分析，计算位点间的LD水平，并绘制LD衰减图，如图5所示，横坐标表示距离，单位Kb，纵坐标为R 平方值。从图中可以看出，牦牛在2Kb水平即从连锁不平衡演变到连锁平衡。即，目标SSR位点需要满足，任意目标SSR位点间距离>2Kb。

6.2群体结构与亲缘关系

用admixture软件，预先定义基因进化簇K的数量从2到10，并经过5 倍交叉实验确定了当K＝7时最优。然后使用admixture软件，构建分别假设存在2～7祖先型下，75个样本的群体结构。由亲缘关系较近的个体组成不同集合，将亲缘关系距离最远的两个集合标记为A、B两组，考察其群内和群间遗传距离上的差异水平。

75个样本群体结构见图6。其中，最左侧，聚类为绿色的标为A组，其亲缘关系较近的集合为28个个体，样本编号分别为：'SRR2059895' 'SRR3112434' 'SRR2059930' 'SRR3112443' 'SRR3112426' 'SRR2059927' 'SRR3112422' 'SRR3112418' 'SRR3112439''SRR3112417' 'SRR3112421' 'SRR3112433' 'SRR2059932' 'SRR3112430' 'SRR2059926''SRR3112437' 'SRR3112436' 'SRR2059929' 'SRR3112425' 'SRR3112415' 'SRR3112440''SRR3112431' 'SRR2059935' 'SRR2059939' 'SRR3112432' 'SRR962827' 'SRR3112441''SRR3112428'

最右侧的2个个体标记为B组，样本编号分别为：SRR2059942和 SRR2059947。

A组和B组群内和群间遗传距离见图7。第一个集合为群内遗传距离，第二个集合为群间遗传距离。群间距离显著高于群内距离。

6.3分析SSR所在区间的多态性水平

分析SSRs在基因Intergenic、CDS、Intron、Downstream、Upstream区间上的多态性水平，绘图并结合F检验，统计检验结果。5个基因区间SSR 位点的多态性水平如图8所示，从图中可以看出CDS区域的多态性水平显著小于其他区域，更为保守。Intergenic和Intron区域的多态性水平最高。

6.4统计SSRs在基因各区域的数目

统计SSRs在基因Intergenic、CDS、Intron、Downstream、Upstream区间内的数目，因为CDS，intron等区域的大小远小于基因间区，所以要进行区间长度矫正。矫正公式：

SSR位点在不同区间类型的数目如图9所示。从数量上看，基因间区显著高于其它区，但是CDS，intron等区间的大小远小于基因间区，所以只能经过长度矫正后再分析，区间长度矫正结果见图10。从图中可以看出，各区域位点数量从大到小依次为Intron、Upstream、Downstream、Intergenic、 CDS，CDS区显著小于其它区间。

通过上述对牦牛所有潜在的SSR位点的分析，得到了潜在的SSR位点在各区间的多态性、各SSR的MAF的分布，在进一步的筛选出可靠的SSR 位点以用于鉴定亲子关系的步骤中，既不宜选择几乎无多态性的SSRs，也不能选择过分多态的SSRs位点。

(七)模拟法筛选目标位点

7.1预实验

用scikit-learn package(v 0.18.0)进行模拟计算，首先进行预实验，考虑SSRs不同的MAF水平，在每个多重等位基因多态基因座上基于随机频率模拟亲本个体，并且基因型频率符合哈迪韦伯格平衡(p值>0.01)。然后基于上述对潜在的SSR位点的过滤和相关特征的分析(多态性、SSR重复类型及其重复性占比、纯合位点及杂合位点的数量、位点分型及MAF值和物种和所述物种的候选SSR位点的群体遗传特征)，以筛选出的SSR位点为选择库，模拟产生20个SSR位点，并根据以下条件模拟随机杂交产生子代群体。

由亲本的个体组成的父代群体满足如下条件：各位点满足LD衰减结果和基因型满足哈温平衡；

杂交满足如下条件：杂交过程完全随机杂交；重组过程无连锁现象；和，个体不产生随机突变。

然后使用定义为相同基因型比率的孟德尔遗传相似性来描述任意两个个体之间的遗传关系。然后，为每两个个体建立一个遗传相似性矩阵。基于相似性矩阵数据，使用ETEtoolkit软件重构系统发生树，然后使用Naive bayes算法预测亲子关系。

相似度水平的模拟计算结果见图11。从图中可以看出，同父同母的子代最接近，其次是父子关系的个体，再次是同父异母或同母异父的个体，完全无血缘关系的个体间，差异最大，亲子关系准确率＝99.97％。

7.2多态性水平对相似度的影响：

然后根据预实验结果，选择模拟300个位点，MAF值取0.2，100个亲本共随机杂交产生500个子代。关系如下：

(1)父代个体间比较：无亲缘关系。

(2)父代个体与子代个体间比较：亲子关系或无亲缘关系；

(3)子代个体比较：无亲缘关系；同父异母或同母异父；或，同父同母。

然后模拟计算相似度水平，并且采用不同的MAF值进行计算，考察多态性水平对相似度的影响。

不同的MAF值下，相似度水平的模拟计算结果见图12A-图12E。MAF 为0.1时，亲子关系准确率＝99.97％，MAF为0.2时，亲子关系准确率＝99.97％，MAF为0.3时，亲子关系准确率＝99.97％，MAF为0.4时，亲子关系准确率＝99.87％，MAF为0.5时，亲子关系准确率＝99.64％。

从上面5个MAF取值结果看，随MAF增大，6组关系的相似度均下降。从父子关系与无关系差异相比，MAF值小于等于0.3较为合适，用以识别亲子关系，准确率均为99.97％，即MAF≤0.3亲子关系准确率均能大于99.9％。可见，相对过于多态的位点，不利于亲权关系的鉴定。

7.3不同分型下CPE变化

在不同规模位点数量的分型标记中(2～6分型)，分别随机测试500次。然后随机模拟一组CPE大于99.99％的五分型SSR标记，进行模型测试。

随机产生500个个体，并产生300个子代。推断其亲子关系的准确率，然后分析不同规模和不同分型类型下的模型推断的结果(进行30次，取平均)。

在不同规模位点数量的分型标记中(2～6分型)，分别随机测试500次的结果见图13。从图中可以看到，无论是什么分型的SSR标记，CPE都与位点数呈正相关，四分型以上的SSR标记，在20个位点左右，CPE就达到了99.99％以上，三分型在30个位点左右CPE能达到99.99％，而二分型的 SSR标记需要50个位点以上才能使CPE达到99.99％。CPE在多分型SSR标记中，明显高于二分型的SSR标记。

7.4亲子关系的准确率测试

模拟随机产生500个父本，300个子代。推断其亲子关系的准确率，结果见表1。从上面的结果看，600个真实的亲子关系，错判的有2例，但是位点较少时，假阳性确实较高，而精确性较差。但是敏感性和特异性都是非常高的。

表1模拟样本亲子关系准确率

sensitivity,recall	specificity	precision	false discovery rate	Accuracy
					99.67％	99.87％	86.18％	13.82％	99.87％

7.5不同分型类型、不同位点规模下CPE的变化

分析不同分型类型和不同位点规模下的模型推断的结果(进行30次，取平均)见表2。从表中可以明显看出，同一分型下，位点数量增加，CPE 明显增加。随着单个位点上，等位基因数目的增多，CPE也明显增加。当CPE趋近于1.0时，sensitivity(敏感性)，specificity(特异性)，precision (精确性)，false discovery rate(伪发现率)和Accuracy(准确率)表现都能令人满意。因此，我们计算单个位点的PE值，以筛选一个能用于预测亲缘关系的SSR标记。

表2不同分型类型、不同位点规模下CPE的变化

(八)位点筛选

根据各项测试结果(PE水平、MAF水平、LD衰减等)筛选可用于牦牛亲权鉴定的SSR位点。

(1)按SSRs位点的PE水平排序。

(2)MAF≤0.3

从MAF取值结果看，随MAF增大，6组关系的相似度均下降。从父子关系与无关系差异相比，MAF值小于等于0.3较为合适，用以识别亲子关系，准确率均为99.97％，即MAF≤0.3亲子关系准确率均能大于99.9％。可见，相对过于多态的位点，不利于亲权关系的鉴定。

(3)缺失值比例(75个个体中缺值)<10％

由于牦牛重测序数据大多是30X，并不是每个位点都能被有效的覆盖到，对一个位点而言，其中70个牦牛个体的被测序read覆盖到了，而另外 5个完全没有被覆盖到，那么该位点的缺失值比率为5/75。当缺失值的比例大于10％时，会导致很多位点被遗漏，群体覆盖度不行，故要求该值，以达到良好的位点数及覆盖度。

(4)任意2个SSRs间的距离>2Kb，即大于LOD衰减距离。

(5)重复单元≥3。

从前述模拟筛选中可以看出，无论是什么分型的SSR标记，CPE都与位点数呈正相关，四分型以上的SSR标记，在20个位点左右，CPE就达到了99.99％以上，三分型在30个位点左右CPE能达到99.99％，而二分型的 SSR标记需要50个位点以上才能使CPE达到99.99％。CPE在多分型SSR 标记中，明显高于二分型的SSR标记。所以从准确性及临床应用的便捷性上，需要选三分型及其以上的。

最后共筛选出50个可用于牦牛亲权鉴定的SSR位点，50个位点的详细信息见表3和表4。从表中可以看出，筛选出的50个SSR分型都在4以上，且大部分在5以上，这符合不同分型下CPE变化的分析结果。50个位点的MAF值都小于0.1，且大部分位点都属于Intron区间。随着位点个数的减少，结合排除概率逐渐减小，前21个位点的CPE可以达到99.99％以上，当位点数取前10时，CPE达到99％，当位点数取前5时，CPE大于 90％。

表3 50个位点及位点信息

表4 50个位点及位点信息

实施例2

样本采集：颈静脉采取活体牦牛血液样本，添加血液抗凝剂EDTA后放入-80℃超低温冰箱冻存。样本为有明确记录的阳性样本，有6个家系，共18头(样本包括父、母、子代)，来自四川省龙日种蓄场，这些样本用来亲本分析，进行临床检验。

从筛选出的50个SSR位点中选择如下20个SSR位点，SSR geneID为 scaffold2072_1、scaffold341_1、scaffold1139_1、scaffold879_1、 scaffold2036_1、scaffold2058_1、scaffold1000_1、scaffold1649_1、 scaffold4112_1、scaffold506_1、scaffold1645_1、scaffold2687_1、 scaffold1210_1、scaffold547_1、scaffold1214_1、scaffold1343_2、scaffold1141_1、scaffold94_1、scaffold738_1和scaffold629_1。

微卫星引物设计：针对20个牦牛亲权鉴定微卫星位点，分别设计PCR 扩增引物，并在每个正向引物的5’端用6-FAM荧光基团进行修饰用于荧光PCR分析，引物由派森诺生物工程股份有限公司合成。

DNA提取：使用改良的缓冲液系统法提取牦牛血液样本的DNA，包括如下步骤：

(1)提取DNA所用的器材要经过高压灭菌，以防止杂质污染。

(2)将冻存血样放入37℃水浴锅中解冻。

(3)取3mL含抗凝剂的血液放入15mL离心管中，加入3mL细胞裂解液(Tris-Cl10mmol/L，pH8.0；EDTA 0.1mol/LSDS 0.5％；不含DNase 的胰RNase 20ug/ml)，充分混匀后，3600rpm离心2min，弃上清。注意混匀时尽量避免产生气泡。

(4)再次加入3mL细胞裂解液，充分混匀至无沉淀，3 600rpm离心 2min，弃上清。

(5)按照10:1的比例混合缓冲液(Tris-Cl 100mmol/L,pH8.0；EDTA 50mmol/L,pH8.0；Nacl 500mmol/L)与蛋白酶K(20mg/mL)。

(6)加入1mL缓冲液与蛋白酶K的混合液，涡旋震荡至无团块，65℃水浴30分钟。

(7)加入1mL异丙醇，颠倒充分混匀至出现丝状或簇状基因组DNA。

(8)3 600rpm离心8min，弃上清。将离心管倒置于干净滤纸上，确保沉淀存在。

(9)加入3mL 70％乙醇，涡旋震荡5s，3600rpm离心3min，弃上清。

(10)将离心管倒置于干净滤纸上5min，确保沉淀存在，然后空气干燥5min。

(11)加入300μL双蒸水，低速涡旋5s，65℃水浴加热1h溶解DNA。

荧光PCR：

使用荧光引物进行PCR，反应体系以及程序设定如表5和表6所示。

表5 PCR体系

表6 PCR程序

纯化：

(1)PCR结束以后，瞬离样品以除去管壁样品，随机挑取样品2μL 进行凝胶电泳。以确定样品浓度、片段大小范围等。

(2)取一新的96孔板标明板号。根据电泳情况调整加样量(如有必要需稀释后加样)，加入冷的70％乙醇至终体积50μL，震荡充分混匀。

(3)3 700rpm/min，4℃离心30min，以纯化样品。倒置瞬离，以除去乙醇。静置15min待乙醇挥发干净。

毛细管电泳及分型：

(1)在乙醇已经挥发完全的板中加入内标LIZ500和Hi-Ditm Formamide，震荡充分混匀，瞬离以除去管壁样品。

(2)放入PCR仪，95℃、4min变性。

(3)放到ABI 3730XL遗传分析仪中进行毛细管电泳及分型。

随机加入部分公牛数据与阳性公牛一起作为候选亲本，运用Cervus 3.0 软件对其进行亲本分析，结合纸质记录，分析所筛选微卫星位点进行亲权鉴定的准确率。

对有明确记录的6个家系18个样本进行亲本分析结果见表7，运用 Cervus 3.0软件中Analysis程序的Parentage Analysis模块进行数据分析， LOD值是亲子指数(Paternity index)的对数值，LOD值大于0的则表示与任意个体相比，候选亲本(Candidateparent)最有可能是真实的亲本；LOD 值小于0表示与任意个体相比，候选亲本不可能是真实的亲本。Cervus会显示最可能的候选亲本。最终分析结果与记录结果一致，且LOD值均为正数，说明了用这些位点进行亲权鉴定的准确性，结果如表7所示。

表7亲本分析

最后应说明的是：以上各实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述各实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围。

Claims

1.一种用于亲权鉴定的SSR位点的筛选方法，其特征在于，包括：

(S1)筛选出待筛选SSR位点的物种的全基因组序列所有潜在的SSR位点，再筛选出位于基因区间内的SSR位点，再过滤掉单核苷酸重复类型的SSR位点；然后筛选出具有多态性的SSR位点；

(S2)然后过滤掉满足如下条件的SSR位点，得到候选SSR位点：(a1)缺失值比例＞50％；(a2)哈温平衡检验显著性＜0.01；

(S3)计算候选SSR位点的LD衰减距离；

(S4)将候选SSR位点按照PE值排序，然后筛选出符合如下条件的SSR位点：(b1)MAF值≤M；(b2)缺失值比例＜10％；(b3)任意两个SSR位点的距离＞LD衰减距离；(b4)重复单元≥N；

其中，所述M和所述N按照如下方法得到：

模拟产生具有SSR位点的亲代个体，所述亲代个体满足如下条件：(c1)亲代个体的SSR位点距离≥LD衰减距离；(c2)亲代个体的基因型满足哈温平衡；

然后模拟计算如下指标：

(e2)不同SSR分型类型、不同位点规模下CPE的变化，以确定SSR位点的重复单元数，记为N。

2.根据权利要求1所述的SSR位点的筛选方法，其特征在于，以待筛选SSR位点的物种的高通量测序数据作为全基因组序列的数据来源。

3.根据权利要求2所述的SSR位点的筛选方法，其特征在于，从NCBISRA数据库中筛选所述物种的相同测序平台下，并且数据质量整齐度一致的测序样本，然后进行数据校准和整理以及重复序列标记和合并；

优选地，所述数据校准和整理包括：收集到的SRA数据首先用sra-tool kit的“fastq-dump-split-3”命令将其转换为FASTQ数据格式，用cutadapt程序去除序列分数＜20的低质量测序片段/残基，得到过滤后测序数据。再用Burrows-Wheeler Alignment(BWA)tool对比工具的“bwa mem-T 30-h 5”程序将过滤后测序数据比对到牦牛基因组上，同时检查比对质量，把质量分值低于30的过滤掉以减少假阳性。然后用Samtools软件将比对结果进行格式转换及排序，产生对比后BAM数据文件格式；

4.根据权利要求1所述的SSR位点的筛选方法，其特征在于，使用scikit-learnpackage模拟产生具有SSR位点的亲代个体。

5.根据权利要求1所述的SSR位点的筛选方法，其特征在于，使用VCFtools对候选SSR位点数据进行LD分析，计算位点间的LD水平，并绘制LD衰减图。

6.根据权利要求1所述的SSR位点的筛选方法，其特征在于，所述筛选方法还包括分析潜在的SSR位点的特征，以用于作为模拟筛选的SSR位点库；

7.根据权利要求1所述的SSR位点的筛选方法，其特征在于，所述筛选方法还包括分析待筛选SSR位点的物种和所述物种的候选SSR位点的群体遗传特征，以用于作为模拟筛选的SSR位点库；

8.根据权利要求1-7中任一项所述的SSR位点的筛选方法，其特征在于，所述待筛选SSR位点的物种为牦牛。

9.根据权利要求8所述的SSR位点的筛选方法，其特征在于，将牦牛的候选SSR位点按照PE值排序，然后筛选出符合如下条件的SSR位点：(b1)MAF值≤0.3；(b2)缺失值比例＜10％；(b3)任意两个SSR位点的距离＞2kb；(b4)重复单元≥3。

10.权利要求1-9中任一项所述的SSR位点的筛选方法在如下(x1)-(x4)中的应用：

(x1)个体识别；

(x2)家系管理；

(x3)种质资源鉴定；

(x4)遗传多态性位点分析。