MonteCarlo方法在遗传育种中的应用

Monte Carlo 方法在遗传育种中的应用Monte Carlo 方法在遗传育种中有着十分广泛的应用，下面举出一些例子。例1 遗传漂变的模拟遗传漂变是指在一个没有选择、迁移、实施随机交配的小群体中，群体基因频率和基因型频率会偏离Hardy-Weinbger平衡。利用Monte Carlo 方法，我们可以模拟随机遗传漂变的过程，研究各种因素，如群体大小、公母比例、繁殖力等，对遗传漂变的影响。具体做法如下：1 设定系统参数N 基础群大小S 基础群公母比例no 每头母畜每胎所产后代数P1、P2、P3 所考查的基因座位在基础群中的基因型频率（假设只有两个等位基因）g 所考查的世代数nr 模拟重复次数2 假定1） 公、母随机交配，2） 每头公畜与相同数目的母畜交配3） 后代为雄性和雌性的概率各为1/2，4） 基因从上代到下代的传递遵从孟德尔分离定律5） 世代不重叠6） 群体规模保持不变7） 每代中随机选留种畜3 模拟试验1） 基础群中每个个体的性别和基因型的确定 设雄性个体的代码为1，雌性个体的代码为2。再设两个等位基因的代码分别为1和2，三种基因型的代码分别为11、12和22。基础群中每个个体的性别和基因型可以硬性规定，也可以随机确定，如硬性规定，即人为指定哪些个体是雄性，哪些个体是雌性，哪些个体的基因型为11，哪些为12，哪些为22。但要注意必须符合事先给定的公母比例和基因型比例。若随机确定，可用以下方法：对于i = 1, , Na. 产生随机数，b. 如 u1 s/(1+s)，则Si = 1，否则，Si = 2c. 如 u2 P1，则Gi = 11如P1 u2 P1+P2，则Gi = 12如 u2 P1+P2，则Gi = 222） 交配组合的确定设雄性个体数为nm，雌性个体数为nf（nm + nf = N）。对于i = 1, ，nm ,对于j = 1, , nf/ nm a. 产生在区间1，nf上的离散均匀分布随机数rb. 若第r头母畜尚未被使用过，则设定第i头公畜的第j个配偶为第r头母畜，否则返回a3） 后代性别和基因型的确定对于每一个交配组合，对于k = 1, , noa. 产生u1 U(0, 1)，若u1 0.5，Sk = 1, 否则，Sk = 2b. 产生u2 U(0, 1)，若u2 0.5，则后代k从父亲处获得第一个等位基因，否则获得第二个等位基因c. 产生u3 U(0, 1)，若u3 0.5，则后代k从母亲处获得第一个等位基因，否则获得第二个等位基因d. 将从父亲和母亲处获得的等位基因合并，即得后代k的基因型4） 种畜的选留a. 在雄性后代中随机选留nm个公畜b. 在雌性后代中随机选留nf 个母畜5） 计算群体中的等位基因频率和基因型频率6） 重复2） 5）直至完成所有世代7） 重复1）6）直至完成所有重复，并计算各世代基因频率和基因型频率的平均数和标准差，进行与Hardy-Weinbger平衡的差异显著性检验。通过改变1 中的有关参数，可考查各种因素对遗传漂变的影响。例2 人工选择的模拟人工选择是家畜育种的中心工作，其最终目的是使家畜群体在重要经济性状获得理想的遗传进展。影响遗传进展的因素总体来说有：性状在群体中的遗传标准差，选择的准确性，选择强度，世代间隔。而其中的每一个因素又要受一些其他因素的影响，例如选择的准确性要受群体结构、估计育种值的方法等因素的影响。对于复杂的实际情况，我们很难从理论上定量地分析各种因素及它们的不同组合对遗传进展的影响，分析遗传进展及群体遗传结构在长期选择过程中的变化，而这些都是我们在制定最优化育种方案时必须要考虑的。这时就可借助Monte Carlo方法，通过模拟选择过程来对这些问题进行探讨。1）设置系统参数N 基础群大小ns 基础群中公畜数nd 基础群中母畜数（是公畜数的倍数）no U (5, 10)随机数，为每头母畜每胎所产后代数ng 选择的世代数ps 每世代选留公畜比例pd - 每世代选留母畜比例- 基础群中性状的遗传方差- 群体中性状的环境方差nr 模拟重复次数2）作系统假设- 基础群中个体间无亲缘关系- 性状表型值用以下模型模拟：其中m 是总平均，gi是第i个世代的固定环境效应（i = 0, 1, , ng ），sj是第j个性别的固定效应（j = 1, 2），ak是第k个个体的随机加性遗传效应（育种值），eijk是随机环境效应，- 所有个体都有一个性状观察值- 选择的依据是动物模型BLUP育种值，所用模型与模拟模型相同- 世代间不重叠- 所有世代中都实施避免同胞交配的随机交配- 环境方差在各世代中都相同3）模拟试验 基础群（0时代）的模拟 对于k = 1, L, Na. 如果k ns，个体k为公畜，否则为母畜 (即规定编号为1 ns的个体为公畜，编号为ns +1 N的个体为公畜)。b. 产生c. 产生d. 令，作为个体k的性状表型值，其中g0是0世代的效应值，s是该个体所在性别的性别效应值。 第j（j =1，2，ng）世代的模拟A. 确定在第j-1世代中选留的公母畜间的交配组合。可按例1中的方法确定交配组合，但要避开同胞交配。B. 对于每一交配组合：产生随机数n0 DU（5，15），作为该交配组合所产生的后代数对于 k =1, ，n0a. 产生u U(0, 1)，如果u 0.5，则个体k为公畜，否则为母畜b. 产生，其中 , fs和fd分别是父亲和母亲的近交系数（从2世代开始，将会有近交个体产生）c. 令 d. 产生e. 令，作为个体k的性状表型值，其中gj是世代j的效应值f. 计算个体k的近交系数C. 计算第j世代中的加性遗传效应的平均值和方差以及平均近交系数， ，其中n是该世代中的总个体数，ai是第i个个体的加性遗传效应，fi是第i个个体的近交系数。D. 估计每个个体的育种值E. 按公畜和母畜的留种比例根据估计育种值选留种畜 重复 nr次，我们可考察在当前的群体结构和系统假设下，在ng个世代中群体遗传进展（即平均加性遗传效应在各世代中的变化情况）、加性遗传方差、近交程度的平均变化情况。改变系统参数和系统假设，我们可考察不同群体结构和系统假设（如群体大小、群体加性遗传方差、留种率、育种值估计方法等）对遗传进展、遗传方差和近交程度的影响。 例3方差组分估计方法MIVQUE和REML的模拟比较在家畜育种中，估计群体遗传参数（方差组分）是一项十分重要的基础工作。由于家畜生产记录资料的复杂性（多为不平衡资料），给方差组分估计带来了很大难度，因此方差组分估计多年来一直是各国统计学家和遗传育种学家的重要研究内容。在70年代后，出现了多种估计方法，其中主要的有最小方差二次无偏估计（MIVQUE）和约束最大似然（REML）法。MIVQUE可以通过解析法求得估计值，其估计值具有无偏性，且当先验值等于真值时，估计值的方差最小。REML只能通过迭代近似求解，因而不能确切求得其估计值的理论偏差和方差。但当样本足够大时，从理论上可以证明其估计值具有渐近无偏性和有效性，因而一般认为对于大样本来说，REML要优于MIVQUE。但是究竟多大样本才能保证REML的渐近无偏性和有效性? 当样本不是足够大时MIVQUE和REML哪一个更好? 对这些问题并不能从理论上给出答案，但用Monte Carlo方法却可以给出至少是经验性的答案。张勤等（1995）用Monte Carlo模拟比较了在不同样本大小和数据结构下MIVQUE和REML估计值的经验抽样特性（无偏性和方差），其基本思想是，在数学模型和估计方法给定的情况下，方差组分估计值的性质只受相应的真值的大小和数据样本大小及其结构的影响，对某一特定的方差组分真值和数据结构，可根据给定的数学模型模拟观测值，再用模拟的数据来估计方差组分，经重复多次模拟（即从同一总体中重复多次抽样）可得方差组分的多次估计值，只要重复的次数足够多，这些估计值的平均数和方差(即经验期望和经验方差)就可近似作为估计值的理论数学期望和方差来看待。该研究考虑了4个与奶牛生产有关的数据集，它们的大小及结构如下表所示。第一个数据集取自北京市当时的荷斯坦牛头胎产奶量实际记录资料，以这个数据集为基础，从中抽取部分数据构成了数据集2、3和4。每个数据集中各头公牛的女儿数、每个场年季中的记录数、公牛之间的亲缘关系以及公牛的分组都是固定的。 3个数据集的样本大小及其结构数据集记录数公牛数场年季数女儿数/公牛记录数/场年季1128474777827321216.516.525097206922552077.47.73135116418841973.22.842002014810+01.40.7对这4个数据集都用以下模型进行模拟和估计方差组分：其中y为所有观察值的向量，1是单位向量，m 是总平均，h是所有场年季效应(固定效应)的向量，X是场年季效应的结构矩阵，g是所有公牛组效应(固定效应)的向量，B是公牛组效应的结构矩阵，s是所有公牛效应(随机效应)的向量，Z是公牛效应的结构矩阵，e是随机残差效应向量，且，其中为公牛方差组分，为残差方差组分，A为公牛间的分子血缘相关矩阵。在模拟中考虑了三组方差组分真值：第一组：第二组：第三组：对每一数据结构和每一组方差组分真值，都用以上模型重复模拟1000次。每次模拟都相当于从总体中进行一次抽样。具体地，可按以下步骤进行模拟：对每一数据结构：1） 假定一个m 值。2）模拟产生h，例如可假定h服从正态分布，且彼此独立。因而可用前节介绍的产生正态分布随机数的方法产生hi (i = 1, 2, L，nh)（ nh是在该数据集中的场年季水平数）。3）模拟产生g，可假定g也服从正态分布，且彼此独立，因而也可用产生正态分布随机数的方法产生gj (j = 1, 2, L，ng)（ ng是在该数据集中的公牛组水平数）。8） 对每组方差组分真值A. 产生公牛效应向量值s由于，可按上节介绍的产生多维正态分布随机数的方法产生s 。B. 对于每个场年季、公牛组和公牛的水平组合 a. 令 f = m + h + g + s b. 对于i = 1, L , n （= 该水平组合中的观察值个数） - 产生随机数e，作为残差效应抽样值 - 令 y = f + e, 作为表型值的抽样值C. 对模拟的表型值分别用REML和MIVQUE法估计方差组分和D. 重复A C 1000次，得到1000个和的估计值，计算这1000个估计值的平均数和方差，作为在该数据结构和参数真值下和估计值的经验期望和经验方差，通过以上模拟，我们可比较REML和MIVQUE这两种方法在不同数据结构和参数真值情况下所得方差组分估计值的偏差（估计值的期望值与真值的差）、精确性（由估计值的方差来度量）和均方误（=偏差2 + 方差）。例4 用于标记-QTL连锁分析的资源群体的模拟Monte Carlo方法在QTL检测和定位以及标记辅助选择的研究中发挥了重要作用，迄今为止所发表的这个领域的研究报道大部分与Monte Carlo方法有关，这一方面是由于Monte Carlo方法在这类研究中有其他方法不可比拟的优越性，另一方面也由于目前还很少有遗传标记的实际资料可供研究。下面介绍几种目前常用的用于QTL检测和定位的试验设计的模拟方法。回交群体设有两个近交系，其中一个在所有标记和QTL位点上都是某一等位基因的纯合体，而另一个在所有标记和QTL位点上都是另一等位基因的纯合体，它们交配后得到F1代个体，F1代个体再与两个近交系中的一个回交，所得的后代即为回交群体。模拟方法如下：1） 设置系统参数N 回交群体的个体数L 染色体的长度（以Morgen为单位）Nm 染色体（连锁群）上的标记数Nq 染色体上的QTL数P 标记和QTL在染色体上的位置（Nm +Nq维向量）sP 性状的表型标准差H2 性状的广义遗传力 遗传标准差 ，环境标准差v2 总的QTL效应对遗传方差的贡献：，剩余多基因效应标准差 si 第i个QTL的效应，表示为两个纯合体的差与表型标准差的比值ti 第i个QTL的显性度，即杂合基因型值与有利纯合基因型值之比2） 作系统假设- 每个QTL位点上有两个等位基因，分别用1和2表示，基因型值可用下图表示对于第i个QTL，2ai = sisP，ai = 0.5 sisP ，di = ti ai 。 在回交群体中，在一个QTL上只有11（或22）和12两种基因型，概率各为0.5，因而在第i个QTL上基因型值的平均值为0.5（ai + di），方差为所有QTL的总方差为因而si和ti的取值必须满足如各QTL的效应相等且无显性效应，则- 各个QTL之间无互作效应（上位效应）- 各QTL的效应相等或服从正态分布或伽吗分布- 性状表型值可用以下模型描述其中qi是第i个QTL的基因型值，u是剩余多基因效应，e是随机环境效应，- 随机环境效应方差在各群体中都相等3） 模拟试验A. 亲代标记-QTL基因型的模拟对于所有标记来说，将每个位点上的两个等位基因用1和2表示。可以将一个亲本定为等位基因1的纯合体，另一亲本定为等位基因2的纯合体。对于QTL来说，由于两个等位基因有不同的效应，这就需要确定每个亲本在每个QTL上是等位基因1的纯合体，还是等位基因2的纯合体。这可以人为规定，也可以随机确定。如果随机确定，可假设某一亲本在某个QTL上为等位基因1的纯合体概率为50%（或其他值）。于是，对于第一个亲本的第i（i = 1, L, Nq）个QTL，先产生U(0,1)随机数u，该QTL的基因型为而另一亲本在该QTL上的等位基因则