建立作物遗传资源核心样品的统计方法

建立作物遗传资源核心样品的统计方法

王磊魏兴华

(中国水稻研究所，浙江杭州 310006)

提要：本文主要介绍了在建立作物遗传资源核心样品中的统计方法：（1）在系统聚类分析中如何同时利用定量和定性变量；（2）非系统聚类方法中的正态混合模型聚类方法，以及（3）利用多环境试验数据进行聚类分析。

关键词：核心样品；聚类方法；正态混合模型；三元分析；模式分析

所谓作物遗传资源核心样品是指以最少量的种质材料代表一个物种及其近缘野生种最大限度的遗传多样性。核心样品的建立，不仅可以提高种质库的利用效率，而且有利于种质库管理、新种质收集、种质创新及种质资源的深层次的研究。核心样品概念在1984年提出后，立即引起了全球广泛的兴趣，先后在西非秋葵、澳州大豆、野生稻等作物得以实践，而我国的核心样品研究始于1994年，近年来有了快速的发展。

建立核心样品中用到的统计方法主要是聚类分析和排序。本文首先介绍同时利用定性和定量变量进行系统聚类分析的Gower方法；然后介绍非系统聚类方法中的一种最常见的正态混合模型聚类方法，最后介绍利用多环境试验数据进行聚类分析的模式分析（Pattern Analysis）和三元分析（Three-Mode Analysis）。

1 同时利用定性和定量变量进行聚类

聚类分析涉及如何结合不同类型的变量，选择合适的距离度量，应用合理的聚类策略，决定最佳的类的数目，并鉴别出有显著区别力度（discriminatory power）的变量。标准的聚类分析使用的或者是定量变量，或者是分类变量，而不是同时在这两类变量上进行聚类。Gower（1971）提出了可同时利用连续和定性变量的两样品之间的距离度量，Wishart（1986）利用该距离度量提出了系统聚类策略。

聚类分析中的样品之间关系有二种：相似或相异系数和距离。而它们的定义与样品的变量的类型关系极大，通常变量按照测量它们的尺度来进行分类：

(1)间隔尺度：指标用连续的量来表示如产量、株高等。如在间隔尺度存在绝对零点，又称比例尺度。

(2)有序尺度：指标度量时没有明确的数量表示，只有次序关系，如米质分成好、中、次三等，三个等次有次序关系，但没有数量表示。

(3)名义尺度：指标度量时既没有数量表示，也没有次序关系。如品种的外观颜色。名义尺度中只取两种状态的变量（binomial）是最简单的一种。

一般说来，通常所说的连续变量或定量变量是指按照间隔尺度度量的变量，而按照有序尺度或名义尺度度量的变量称为定性变量。Cole-Rodgers 等（1997）提出了能同时利用定性和定量变量计算样品之间欧氏距离的统计方法：标准两元度量方法（normed binary scale），进而对样品进行聚类分析。对定性（名义尺度）变量，譬如外观颜色变量C有四种颜色：绿色（G），棕色（B），红（R），和黄（Y），则生成四个子变量CG，CB，CR和CY，如某一样品的颜色为棕色和

红色，则，而CG=0，CY=0，其中4为颜色数目。如某一样品

A1的颜色为B和R，而另一样品A2的颜色为B和Y，则由颜色变量而决定的该两样品之间的相异系数为：

。

事实上该相异系数为。而对于定量变量，譬如产量Yield，在

计算欧氏距离之前先进行变量转化。如样品中的最低和最高产量分别为

min(yield)和max(yield)，变幅R=max(yield)－min(yield)，则变换为

，其实在计算欧氏距离时，变换只需。

p个两元性状在两样品i和j出现或不出现的记录组成的向量记为 (x1，…，x p) 和(y1，…，y p)，其中x k=1或0表明性状k在样品i出现或没有出现。

定义：

分别表示 (x i, y i) =(1,1)， (0,1)，(1,0) 和 (0,0) 的频率。两样品的相似系数可定义为：s1= 或s2 = 。

当同时有定性和定量的变量时，Gower（1971）提出了两样品i和j之间的相似系数：

如果第k个变量为定性变量，则w k=1，不然，R k是第k个变量的变幅。而距离。对更一般的情形，

，其中d ijk是第k个性状对样品i和j之间距离的贡献，

w ijk是d ijk对应的权重，当变量值缺失，w ijk=0。。对于名义尺变变量，如果两样品有相同的值，d ijk=0，不然d ijk=1。对于有序尺度和间隔尺度变量，

。

进一步，对于Gower距离，Wishart（1986）给出了如何计算类之间距离。第k 个变量在第s类的均值和权重为:

类i和j之间的距离（对第k个变量，非名义尺度变量）

Manbattar距离：

欧氏距离平方：

对于名义尺度变量（假设有t个状态）X ikl：为在第i类的第k个名义尺度变量的第l个状态出现的比例。注意变量对总的距离的贡献是

上面两等式中的分母2是为了保证u ijk的值在0和1之间。

2正态混合模型聚类方法

非系统聚类方法，比如像正态混合模型（Normal Mixture Model，可简写为GM模型或Nomix模型）和密度插索方法，则用于描述种质资源中的表现型多样性的模式，并进行聚类，进而建立核心样品。下面主要介绍极大似然方法和Nomix 模型聚类。

2.1极大似然方法

假定样本数为n，样本x1…，x n相互独立，假设来自g个可能的密度函数为f（x；θk）（k=1，2，…，g）的子群体中的某一个。记r=(r1，…，r n)'为一组鉴别标签（identifying labels），r i=k表明x j来自第k个子群体，i=1，…，n；k=1，…，g 。如果C k={x i：r i=k}，k=1，…，g，似然函数：

L(r；θ1，…，θg)=

假设极大似然估计（m.l.e.）为，而对应的一个划分为。如将中的一个样本点移入，似然函数值将减少：

所以

当密度函数为正态分布时，即的密度函数。很容易解得 k

和的m.l.e.：

，

和S k分别是中的n k个样本的均值和协方差。

所以r的m.l.e.是使极小的一种聚类。

如果假定（未知）这时r的m.l.e.为使|w|极小的聚类，其中

在上面的讨论中假设g是已知的。所以我们需要先确定g。当方差一致时，Marriott（1971）建议g取使为最小的值。经验公式为。

2.2混合模型

假设x k以概率p k来自第k个子群体；k=1，2，…，g那么x1，…，x n是来自f（x）的一个样本

这里的密度函数。

当

其中

如果，则x i属于第k个子群体。

那么类的数目g如何确定呢？这可以通过似然比检验确定：

通过Monte Carlo研究，Wolfe推荐近似关系：

，f=2p(g’-g)

注意这一混合模型其实是和前面的极大似然方法几乎是等价的，只是在极大似然方法中r是一个不可观测的随机变量。

3 多环境试验数据分析

保存和利用遗传资源的一个重要步骤是对其进行评价。当样品在不同环境进行评价，通常用聚类分析对样品进行分类，进而建立核心样品，而更常见的是同时利用一些多元统计方法，比如主成分分析，用于描述种质资源中的表现型多样性的模式，并进行聚类，进而建立核心样品。

3.1模式分析

所谓模式分析是指在单一分析变量的基础上，利用由样品和环境（G?E）组成的矩阵对样品和环境同时进行聚类和排序的一种方法。在对品种和环境聚类的基础上，将转化后的G?E矩阵简化为新的行数和列数都已变小的G?E矩阵，再通过排序，对品种和环境，以及样品和环境的互作进行更深入的分析。

一般对G? E矩阵需进行转化分析，一些常用的转化为：

(1)环境中心化

2.(2)环境标准化

(3)均值磨和

3.(4)环境标准残差

下面将转化用公式表达，并更进一步解释转化的具体意义。假设有n个参试品种p个地点，则在第j个地点的第i个品种的产量可用下列模型表述（为解释方便计算公式中不包括随机误差项）：

y ij = μ + g i + e j + ge ij

其中ge ij为第i个品种与第j个地点的互作。

1.环境中心化：y = y ij -

= g i + ge ij

这一转化剔除了环境效应，但仍保留品种的适应性（品种主效和互作之和）。相应的环境聚类则是根据环境对品种的适应性模式的鉴别影响力，而品种的聚类则是根据品种的适应性模式，从变序G? E互作的角度，聚合在一起的品种不是与性状的平均值接近的品种，而是彼此地点之间变序较少的品种。但环境之间的标准差通常不一致，采用下一转化更合适。

(2) 环境标准化：y= (y ij - ) / s j = (g i + ge ij) / s j = g i / s j + ge ij / s j

其中s j是环境j的品种表现的标准差：s j= ，其中n j为环境j 的参试品种。

(3) 均值磨和：y = y ij - = ge ij

每一环境的品种表现的变异是互作的重要方面，所以没有再进一步对每一环境标准化。

(4) 环境标准残差：y =（ y ij - ）/ s

其中s =，其中n j为环境j的参试品种。这一转化的解释与上一转化基本类似，只是这一转化的分析可能更可靠些。

3.2 聚类和排序

对品种和环境的聚类的原则是聚合在一起的类的品种或环境有类似的模式。品种的聚类是按照品种对环境的（转化后的）反映，而环境的聚类是按照环境对品种的鉴别影响。目前通常推荐使用的聚类方法为Ward方法，这一方法也叫做最小离差平方和方法，每次聚类是将类内平方和增加为最少的对象聚合。

那么，最终的品种或环境的类数怎么定呢? 因为每一聚类都损失一部分信息，所以我们希望通过聚类原始矩阵得到简化，同时损失的信息尽可能的小。一般的原则是损失的信息控制在15%和25%之间。

为研究品种在各环境中的表现模式，需对在转化矩阵和聚类后简化的矩阵进行主成分分析，便可对品种或环境进行排序。常用的图形工具是在同一图上同时给出地点和品种的图标的双标图（biplot）。

3.3 三元分析（Three-Mode Analysis）

在对样品进行评价时，对性状单独进行分析是不够的，需要对多个性状同时进行分析。聚类分析所使用的方法为混合极大似然聚类法，而排序是三元主成分分析。这里所要处理的数据为样品?地点?性状的三向（Three-way）数据，而所谓的三元（Three-mode）对应于样品、地点和性状。

3.3.1混合极大似然聚类方法

这是非系统聚类方法，根据性状和地点的信息对样品进行聚类。假设样品来自一固定数目群体的各种比例的混合样本，而各群体（类）可以有不同的均值向量和协方差矩阵，这样同一类中各性状的相互关系的关系结构得到了考虑。

假设有e个环境，a个样品，来自c个群体，第i个样本的观察值X i1=(X i1，…X ie)，X i1,…，X ie为该样本在e个环境的所有性状的观察值向量。假设这些向量X1，X2，…，X a相互独立，并来自C个群体未知比例P m（m=1，…，c）的混合群体。如果X i来自类m，则X ij的密度函数为N（μmj，Σm），（i=1，…，a；j=1，…，e；m=1，…，c）。

上面所提到的类群参数P m，μmj和Σm可以通过应于EM算法求得极大似然估计，将这些极大似然估计值代入带有观察值X i的样本i属于某一群体（类）的后验概率公式，从而计算出相应的后验概率值，样本i属于最大后验概率值对应的类。

3.3.2 三元主成分分析

第一元的分量个数分别记为P（样本），Q（环境）和R（性状），相应的成分矩阵分别记为A，B和C，第i个样本在环境j的第k个性状值

其中a ip，b iq，和c kr分别为A、B和C的元素，而g p q r

为三向（Three-way）核心矩阵G的元素，e ijk为随机误差。通过提供由三元（modes）的成分的第一组合解释的变异比例，核心矩阵给出了三元（样品、地点和性状）的主成分（第一，第二，第三，等等）的相互关系的强度。

三元模型拟合数据的合适程序和不同成分（以及组合）的重要性可以通过计算总的及单独的元（mode）的拟合平方和进行评价。增加模型中的成分自然使拟合程度增加，但相应模型的参数增加使模型的解释变成复杂。为解释模型自然希望能在低维（一般是二维空间）表达成分的结构。最有效的图形是双标图，和类似于双标图的联合图形（Joint plot）。

4 Internet网上软件资源

4.1三元分析（The three-mode analysis）

http://ruls01.fsw.leidenuniv.nl/~kroonenb/：居然有人还专门成立了一个三元公司（The Three-Mode Company），专门致力于三元分析。既然是公司，则开发的软件自然是商业版本，但可免费下载一个能分析小样本的三元主成分分析的演示软件。

4.2混合正态模型聚类分析

https://www.360docs.net/doc/689698519.html,/software.htm/：澳大利亚的Queensland大学的生物统计系提供了可免费下载的DOS环境下的利用混合正态模型聚类的程序：

Mixclus2 –两向两元（two-way and two mode）的聚类分析和Mixclus3 –三向三元(three way and three-mode)的聚类分析。

4.3聚类分析软件Clustan

https://www.360docs.net/doc/689698519.html,/ClustanPackage.htm/: 这是一家专门致力于开发聚类分析的软件公司Clustan Ltd的网页，该公司开发的软件Clustan具备非常丰富的聚类分析方法。只可惜是一个商业软件。

参考文献

1 魏新华，颜启传，应存山，张丽华，章林平. 建立浙江地方籼型稻种资源样品的研究.中国水稻科学，

1999，13（2）：81- 85

2.Basigalup, D.H., D.K. Barnes, and R.E. Stucker. Development of a

core collection for perennial

Medicago plant introductions. Crop Sci., 1995,35:1163-1168.

3.Cole-Rodgers, P. D.W. Smith, and P.W. Bosland. A novel statistical

approach to analyze genetic

resource evaluations using Capsicum as an example. Crop Sci.,

1997,37:1000-1002.

4.Diwan, N., M.S. McIntosh, G.R. Bauchan. Methods of developing a core

collection of annual

Medicago species. Theor. Appl.Genet., 1995, 90: 755-761

5.Franco, J., J. Crossa, J. Diaz et al.. A sequential clustering

strategy for classifiying gene

bank accessions. Crop Sci. 1997, 37:1656-1662.

6.Franco,J., J.Crossa, J. Villasenor, et al.. Classifying Mexican

maize accessions using

hierarchical and density search methods. Crop Sci., 1997, 37:972 – 980.

7.Franco,J., J.Crossa, J. Villasenor, et al.. Classifying genetic

resources by categorical and

continuous variables. Crop Sci., 1998, 38:1688 – 1696.

8.Franco,J., J.Crossa, J. Villasenor, et al.. A two-stage, three-way

method for classifying

genetic resources in multiple environments. Crop Sci., 1999, 39: 259 – 267.

9.Taba, S., J. Diaz, J. Franco, and J. Crossa. Evaluation of Caribbean

maize accessions to develop

a core subset. Crop Sci. 1998, 38:1378-1386.