CN111354417B - 一种基于admixture-mcp模型估计水产动物基因组品种组成的新方法 - Google Patents

Abstract

本发明提供一种基于ADMIXTURE‑MCP模型估计水产动物基因组品种组成的新方法，属于水产动物品种组成领域，方法包括步骤：获得未知品种水产动物基因组SNP标记数据；对参考水产动物进行筛选；估计水产动物基因组品种组成。针对以往的ADMIXTURE模型对不同种扇贝的相似基因序列产生小分量的问题，本发明使用稀疏优化方法产生稀疏解来降低传统模型由基因相似性而产生的估计GBC噪音。我们发现传统L1范数约束模型产生的稀疏解不稳定，虽然可以减少ADMIXTURE小分量的问题，但估计GBC组成与真实GBC组成往往会发较大偏离。为了克服传统算法使用L1范数约束导致严重GBC估计偏差的缺陷，本发明提出使用非凸函数MCP约束来替代L1范数约束，其能够在减少GBC噪声同时，更能产生更接近无偏估计GBC结果。

Description

一种基于ADMIXTURE-MCP模型估计水产动物基因组品种组成 的新方法

技术领域

本发明属于水产动物品种组成领域，具体涉及一种基于ADMIXTURE-MCP模型估计水产动物基因组品种组成的新方法。

背景技术

利用基础科学研究推动畜牧业高质量发展并为消费者提供高质、低价和安全的畜禽产品是动物基因组学研究的出发点与落脚点。动物个体的基因组品种构成(genomicbreed composition,GBC)是指动物祖先对于动物个体基因组的遗传贡献比例。对动物个体GBC的准确估计对动物遗传学研究和育种工作具有重大意义，如可以从基因组水平实现动物亲缘关系鉴定、评估某动物品种的育成历史和品种纯度、保护地方特色品种、预测杂种优势(即评估杂交品种对特定生产环境的适应性)从而实现杂交计划和杂交育种方案的精准设计。通过系谱信息确定品种组成的传统方法具有一定的局限性，其评估的可靠性可能会受到系谱缺失、不准确或不完整记录的影响。例如，基于系谱的预估函数对于那些父母具有相同品种组成的个体会产生相同的估计值。众所周知,对于不同于F1代的杂交个体而言,由于在减数分裂过程中发生的染色体重组与交叉互换，导致在同一品种群体中品种的实际组成各不相同。因此,基于系谱的品种组成估计可能与实际组成大不相同。利用基因组数据来估计基因组品种组成(GBC)，能够在基因组水平上衡量父母的实际贡献，并有助于纠正谱系错误，这种方法比使用谱系记录更有效。GBC的祖先和品种起源基本上是用微卫星标记、单核苷酸多态性(SNP)和序列数据估计的。

目前分析GBC的方法都使用遗传数据分析不同的种群:包括基于模型的聚类方法，如STRUCTURE和ADMIXTURE，以及主成分分析(PCA)。大多数基于模型的聚类方法对多个个体的祖先进行无监督分析，并一同估计祖先群体的等位基因频率。除此之外，ADMIXTURE的监督分析对于准确估计已知祖先历史个体的ADMIXTURE分数是很有价值的，其采用适当的线性回归方法(最小二乘法)估计SNP基因型的品种组成。

但是，由于同一水产动物的不同品种在进化上的同源性(来源于同一祖先)，使得不同品种水产动物的基因序列具有一定的相似度，因此在用现有的ADMIXTURE模型估计水产动物品种组成时，模型会对这些相似的基因序列产生一些小的分量，换言之，并不是所有纯种动物都有100％的GBC，如果按照100％的GBC作为判断纯种扇贝的标准，那么这些小分量的存在会导致鉴定结果具有很高的假阴性，若直接使用较低的估计值作为判断标准，就会大大降低估计结果的准确性。

发明内容

本发明提出了一种基于ADMIXTURE-MCP模型估计水产动物基因组品种组成的新方法。针对以往的ADMIXTURE模型对不同种扇贝的相似基因序列产生小分量的问题，本发明使用稀疏优化方法产生稀疏解来降低传统模型由基因相似性而产生的估计GBC噪音。我们发现传统L1范数约束模型产生的稀疏解不稳定，虽然可以减少ADMIXTURE小分量的问题，但估计GBC组成与真实GBC组成往往会发较大偏离。为了克服传统算法使用L1范数约束导致严重GBC估计偏差的缺陷，本发明提出使用非凸函数MCP约束来替代L1范数约束，其能够在减少GBC噪声同时，更能产生更接近无偏估计GBC结果。

本发明提供一种基于ADMIXTURE-MCP模型估计水产动物基因组品种组成的新方法，所述方法包括步骤：

1)获得某一水产动物未知品种的基因组SNP标记数据；

2)对步骤1)中所述的SNP数据进行筛选，筛选出的SNP要有较高的信息含量，可以准确鉴定和区分现有的水产动物品种；

3)对步骤2)中所筛选的SNP使用均态分布的SNP，减少或删除处于高度连锁不平衡的SNP；

4)对参考水产动物进行筛选；

5)对于步骤4)中筛选的水产动物应为纯种或者是每个品种中代表性的动物个体，认为这些动物个体SNP的基因型具有独立的多项分布，用以计算某一水产动物的某品种属于某些特定品种的可能性。

6)估计水产动物基因组品种组成；

7)在对步骤6)所述计算品种组成时，用G_i代表个体基因型在SNP等位基因的数量(0、1或2)，q_ij表示j种群中第i个SNP处等位基因的频率。已知k参考值或已知等位基因频率的祖先群体，让x_i代表第j个种群的ADMIXTURE比例，并且χ＝[x₁,x₂,…,x_k]是ADMIXTURE系数的向量。在给定等位基因频率和ADMIXTURE比例的情况下，将

定义为SNPi的权重等位基因频率。

8)在步骤7)中所述群体遵循Hardy-Weinberg平衡法则(HWE)，观察基因型Gi在i位点的概率为：

对于给定ADMIXTURE比例的向量，个体所观察到的基因型g的对数似然可定义为：

其中g_i是位点i观察到的基因型。上述似然函数(1)也可写成各位点基因型的函数：

其中C是常数。对于k群体给定的等位基因频率矩阵q_ij(1≤i≤nand1≤j≤k)，目标为确定当ADMIXTURE系数受L(χ)最大约束(x_j≥1和∑_jx_j)时，向量χ＝[x₁,x₂,…,x_k]的估计。

进一步的，对步骤1)中使用LD的r2>0.30作为删除SNP的尺度。

进一步的，对步骤7)中传统模型ADMIXTURE对某一水产动物不同品种的相似基因序列产生小分量的问题进行了优化，使用MCP稀疏优化方法产生稀疏解来降低估计GBC噪音。

进一步的，对步骤8)中通过最大化数据稀疏性对参数的非凸惩罚MCP的似然对数，得到模型(2)的稀疏解的估计:χ{x_j}(j＝1,…,k)

其中λ(λ＞0)并且

MCP的参数为凹性参数。

附图说明

通过以下附图对本发明进行说明

图1是本发明的传统ADMIXTURE模型估计GBC结果图。

图2是本发明的稀疏正则化ADMIXTURE-MCP模型估计GBC结果图。

图3是本发明方法实施流程图。

具体实施方式

以下结合附图和实施例对本发明进一步的说明。

本发明提供一种基于ADMIXTURE-MCP模型估计水产动物基因组品种组成的新方法，方法包括步骤：

1)获得某一水产动物未知品种的基因组SNP标记数据；

2)对步骤1)中的SNP数据进行筛选，筛选出的SNP要有较高的信息含量，可以准确鉴定和区分现有的水产动物品种；

3)对步骤2)中所筛选的SNP使用均态分布的SNP，减少或删除处于高度连锁不平衡的SNP；

4)对参考水产动物进行筛选；

5)对于步骤4)中筛选的水产动物应为纯种或者是每个品种中代表性的动物个体，认为这些动物个体SNP的基因型具有独立的多项分布，用以计算某一水产动物的某品种属于某些特定品种的可能性。

6)估计水产动物基因组品种组成；

7)在对步骤6)计算品种组成时，用G_i代表个体基因型在SNP等位基因的数量(0、1或2)，q_ij表示j种群中第i个SNP处等位基因的频率。已知k参考值或已知等位基因频率的祖先群体，让x_i代表第j个种群的ADMIXTURE比例，并且χ＝[x₁,x₂,…,x_k]是ADMIXTURE系数的向量。在给定等位基因频率和ADMIXTURE比例的情况下，将

定义为SNPi的权重等位基因频率。

8)在步骤7)中群体遵循Hardy-Weinberg平衡法则(HWE)，观察基因型Gi在i位点的概率为：

对于给定ADMIXTURE比例的向量，个体所观察到的基因型g的对数似然可定义为：

其中g_i是位点i观察到的基因型。上述似然函数(1)也可写成各位点基因型的函数：

其中C是常数。对于k群体给定的等位基因频率矩阵q_ij(1≤i≤nand1≤j≤k)，目标为确定当ADMIXTURE系数受L(χ)最大约束(x_j≥1和∑_j x_j)时，向量χ＝[x₁,x₂,…,x_k]的估计。

作为发明的进一步方案，对步骤1)中使用LD的r2>0.30作为删除SNP的尺度。

作为发明的进一步方案，对步骤7)中传统模型ADMIXTURE对某一水产动物不同品种的相似基因序列产生小分量的问题进行了优化，使用MCP稀疏优化方法产生稀疏解来降低估计GBC噪音。

作为发明的进一步方案，对步骤8)中通过最大化数据稀疏性对参数的非凸惩罚MCP的似然对数，得到模型(2)的稀疏解的估计:χ{x_j}(j＝1，…，k)

其中λ(λ＞0)并且

MCP的参数为凹性参数。

实施例

下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述。以下用水产动物中双壳贝类扇贝为实例，用于说明本发明，但不用来限制本发明的范围。

海湾扇贝是我国最重要的养殖贝类之一，但是近年来出现的种质退化问题成为制约我国海湾扇贝养殖业发展的严重障碍。紫扇贝和海湾扇贝同属于海湾扇贝属的优良品种，存在互补的性状，通过种间杂交有望培育出生长速度快、个体大且温度适应范围广的扇贝，可以改善海湾扇贝的种质，提高扇贝养殖业的经济产量。

紫扇贝与海湾扇贝的成功杂交为改善海湾扇贝的种质提供了新的途径。杂交F₁代中存在着少量雌性可育个体，可通过回交的方法培育生长优势极其显著的杂交扇贝后代，但紫扇贝和海湾扇贝均为雌雄同体的动物,在育种实践中存在着精子或卵子污染的问题,为确保育种进程的顺利进行并确证后代的谱系,需要确定某一后代是否为真正的杂交后代及其亲本来源。

步骤1：建立杂交扇贝家系。本发明实施例用控温促熟的方法培育海湾扇贝，使其与紫扇贝亲贝同时成熟，用温度刺激的方法使紫扇贝和海湾扇贝同时排放精卵，分别获得紫扇贝的精卵和海湾扇贝的精卵，建立以下家系组合：(1)紫海杂交扇贝(紫扇贝♀×海湾扇贝♂)；(2)海紫杂交扇贝(海湾扇贝♀×紫扇贝♂)。在回交试验中，同样用控温促熟的方法分别获得海紫杂交扇贝和紫海杂交扇贝的卵子及紫扇贝和海湾扇贝的精子，建立以下家系组合：(1)海紫×紫杂交扇贝(海紫杂交扇贝♀×紫扇贝♂)；(2)紫海×海杂交扇贝(紫海杂交扇贝♀×海湾扇♂)；(3)紫海×紫杂交扇贝(紫海杂交扇贝♀×紫扇贝♂)；(4)海紫×海杂交扇贝(海紫杂交扇贝♀×海湾扇贝♂)。用传统的海湾扇贝育苗方法对这几个家系组合进行幼虫培育，中间保苗和养成，4个月后对各组随机取贝进行试验。

步骤2：扇贝全基因组DNA的提取。取待鉴定扇贝闭壳肌，保存在75％的乙醇中于-20℃备用。用天根DNA提取试剂盒(天根生化科技有限公司，北京)分别提取全基因组DNA，用TE缓冲液稀释至终质量浓度100ng/μL保存于-20℃备用，提取的全基因组DNA经1.2％琼脂糖凝胶电泳，溴化乙锭染色，凝胶成像系统照相观察。

步骤3：获得扇贝基因组SNP数据。本发明实例应用2b-RAD高通量简化基因分型技术获得不同家系扇贝个体基因组SNP标记数据，并使用LD的r²>0.30作为删除SNP的尺度，删除高度连锁不平衡的SNP，筛选出有较高的信息含量的SNP数据，用以准确鉴定和区分扇贝品种。

步骤4：对参考扇贝进行筛选。筛选的扇贝应为每个家系中代表性的扇贝个体，认为这些扇贝个体SNP的基因型具有独立的多项分布，用以计算某种扇贝属于某一特定品种的可能性。

步骤5：估计扇贝基因组品种组成。

步骤5.1：本发明实例用G_i代表紫扇贝、海湾扇贝以及各杂交家系扇贝个体基因型在SNP等位基因的数量(0、1或2)，q_ij表示j种群中第i个SNP处等位基因的频率。已知k参考值或已知等位基因频率的祖先群体，让x_i代表第j个种群的ADMIXTURE比例，并且χ＝[x₁,x₂,…,x_k]是ADMIXTURE系数的向量。在给定等位基因频率和ADMIXTURE比例的情况下，将

定义为SNPi的权重等位基因频率。群体遵循Hardy-Weinberg平衡法则(HWE)，观察基因型Gi在i位点的概率为：

对于给定ADMIXTURE比例的向量，个体所观察到的基因型g的对数似然可定义为：

其中g_i是位点i观察到的基因型。上述似然函数(1)也可写成各位点基因型的函数：

其中C是常数。对于k群体给定的等位基因频率矩阵q_ij(1≤i≤nand1≤j≤k)，最终计算结果为当ADMIXTURE系数受L(χ)最大约束(x_j≥1和∑_jx_j)时，向量χ＝[x₁,x₂,…,x_k]的估计。

通过最大化数据稀疏性对参数的非凸惩罚MCP的似然对数，得到模型(2)的稀疏解的估计:χ{x_j}(j＝1,…,k)

其中λ(λ>0)并且

MCP的参数为凹性参数；

上述步骤参数需通过交叉验证得到模型的适当正则化参数值，其次，通过惩罚每个个体的小的ADMIXTURE系数，具有最佳正则化参数的稀疏正则化ADMIXTURE模型更为简约，并且对于具有大量祖先的扇贝种群能产生更可信的结果，最后，利用已知祖先的个体数据对群体进行更精确的祖先估计。

步骤5.2：按照传统ADMIXTURE模型将上述扇贝个体基因组品种组成进行估计。

结果：在估计不同家系扇贝个体基因组品种组成时，MCP稀疏正则化ADMIXTURE模型较传统ADMIXTURE模型有更精准的GBC估计，传统ADMIXTURE模型在估计GBC时会产生较大的噪音，如图1所示，海紫×紫杂交扇贝中传统方法估计GBC组成海湾扇贝占比接近30％，紫扇贝占比接近50％，而由于两种扇贝基因相似性而产生的噪音扇贝高达20％。MCP正则化产生稀疏解很好的降低了由不同种扇贝基因相似性而引起的噪音，将不同家系(1)海紫×紫杂交扇贝(海紫杂交扇贝♀×紫扇贝♂)；(2)紫海×海杂交扇贝(紫海杂交扇贝♀×海湾扇♂)；(3)紫海×紫杂交扇贝(紫海杂交扇贝♀×紫扇贝♂)；(4)海紫×海杂交扇贝(海紫杂交扇贝♀×海湾扇贝♂)的数据集随机分为四个动物子集。随机选取三个子集用作训练样本(训练集)来估计模型中的品种等位基因频率，另外一个子集(测试集)以测试模型的GBC预测能力。通过交叉验证得到合适的参数λ后，MCP稀疏正则化ADMIXTURE模型比传统ADMIXTURE模型给了纯种动物更高的比例，如图2所示，海紫×紫杂交扇贝中新方法估计GBC组成海湾扇贝占比接近35％，紫扇贝占比接近62％，而噪音降至3％以下，这使得扇贝GBC估计更为准确。

Claims (4)

1.一种基于ADMIXTURE-MCP模型估计水产动物基因组品种组成的新方法，所述方法包括步骤：

1)获得某一水产动物未知品种的基因组SNP标记数据；

2)对步骤1)中所述的SNP数据进行筛选，筛选出的SNP要有较高的信息含量，可以准确鉴定和区分现有的水产动物品种；

3)对步骤2)中所筛选的SNP使用均态分布的SNP，减少或删除处于高度连锁不平衡的SNP；

4)对参考水产动物进行筛选；

5)对于步骤4)中筛选的水产动物应为纯种或者是每个品种中代表性的动物个体，认为这些动物个体SNP的基因型具有独立的多项分布，用以计算某一水产动物的某品种属于某些特定品种的可能性；

6)估计水产动物基因组品种组成；

7)在对步骤6)计算品种组成时，用G_i代表个体基因型在SNP等位基因的数量，所述G_i为0、1或2，q_ij表示j种群中第i个SNP处等位基因的频率，已知k参考值或已知等位基因频率的祖先群体，让x_i代表第j个种群的ADMIXTURE比例，并且

是ADMIXTURE系数的向量，在给定等位基因频率和ADMIXTURE比例的情况下，将

定义为SNPi的权重等位基因频率；

8)在步骤7)中所述群体遵循Hardy-Weinberg平衡法则(HWE)，观察基因型Gi在i位点的概率为：

对于给定ADMIXTURE比例的向量，个体所观察到的基因型g的对数似然可定义为：

其中g_i是位点i观察到的基因型；

上述似然函数(1)也可写成各位点基因型的函数：

其中C是常数；

对于k群体给定的等位基因频率矩阵q_ij(1≤i≤nand1≤j≤k)，目标为确定当ADMIXTURE系数受

最大约束(x_j≥1和∑_jx_j)时，向量

的估计。

2.根据权利要求1所述的方法，对步骤1)中使用LD的r2>0.30作为删除SNP的尺度。

3.根据权利要求1所述的方法，对步骤7)中传统模型ADMIXTURE对某一水产动物不同品种的相似基因序列产生小分量的问题进行了优化，使用MCP稀疏优化方法产生稀疏解来降低估计GBC噪音。

4.根据权利要求1所述的方法，对步骤8)中通过最大化数据稀疏性对参数的非凸惩罚MCP的似然对数，得到模型(2)的稀疏解的估计:

其中

λ(λ>0)并且

MCP的参数为凹性参数。

Images