GLMM：广义线性混合模型（遗传参数评估）

0. 飞哥感言

这篇文章，主要是介绍了抗性数据，如何利用GLMM模型进行的分析，文中，他将9级分类性状变为了二分类性状，进行分析。

分析中用到了加性效应（A矩阵），空间分析（行列信息）。

对比了SAS和ASReml，结果基本一致。

其实，9分级性状，可以直接使用ASReml进行有序多分类性状分析，用累计Logistic模型分析，也可以考虑系谱数据和空间位置信息。这样效果应该更好。

回头找下数据，测试一下。

---

1. 文献

Genetic analysis of resistance to Pseudomonas syringae pv. actinidiae (Psa) in a kiwifruit progeny test: an application of generalised linear mixed models (GLMMs)

❝De Silva N H , Gea L , Lowe R . Genetic analysis of resistance to Pseudomonas syringae pv. actinidiae (Psa) in a kiwifruit progeny test: an application of generalised linear mixed models (GLMMs)[J]. Springerplus, 2014, 3(1):1-12. ❞

2. 摘要

「要点：」

• LMM模型可以结合遗传（G矩阵）和空间分析（误差R矩阵），估算BLUP值
• SAS中的GLIMMIX可以处理GLMM模型，但是门槛较高
• ASReml可以处理GLMM模型

❝线性混合模型（LMMs）结合了遗传和空间协方差结构，在动植物育种中用于估计遗传参数和预测育种值。虽然将LMM扩展到广义线性混合模型（GLMMs）的理论方面已经有一段时间了，但是直到最近十几年才开发出合适的软件。SAS®中的GLIMMIX程序越来越流行，可以在各种学科中安装GLMMs。GLMMs在遗传分析中的应用受到限制，可能是因为所用模型的复杂性。Proc-GLIMMIX尤其如此，因为与ASReml软件不同，它不是专门为分析育种数据而量身定做的，而且需要一些预处理程序编码。符合GLMM框架的二进制数据通常在育种实验中遇到，例如通过观察疾病的存在与否来评估个体的抗性。丁香假单胞菌引起的细菌性溃疡病。猕猴桃是新西兰和其他猕猴桃生产国的一种严重的猕猴桃病害。来自后代试验试验的数据可用于鉴定具有高抗性育种价值的亲本。为此，我们成功地应用了GLIMMIX程序。抗性遗传力中等，经鉴定，两个亲本及其家系具有较高的Psa抗性育种潜力。在使用二进制数据的GLMMs时，有几个潜在的陷阱，并简要讨论这些陷阱。 ❞

3. 介绍

「要点：」

• 介绍抗病亲本的利用
• 如何更有效的评估和利用

❝丁香假单胞菌。猕猴桃（actinidiae，Psa）是猕猴桃的一种病原菌。自从2010年11月在新西兰首次发现Psa-V病毒后，该病毒的毒力形式现已在新西兰建立（Everett et al。2011年）。从生物学角度来看，Psa-V是一系列疾病症状的致病因子，包括叶斑和坏死、花萎蔫、甘蔗枯萎、枝干溃烂，在易感品种的情况下，常常导致葡萄树死亡(http://www.kvh.org.nz/). 这种疾病对新西兰猕猴桃产业造成了极大的损害：预计五年内净现值的成本将在3.1亿至4.1亿新西兰元之间（Greer和Saunders，2012年）。造成损失较大的主要原因是二倍体黄肉猕猴桃品种‘Hort16A’的完全易感性。相比之下，绿色肉质的美味多年生品种‘海沃德’表现出一定程度的抗性，这使得这种疾病可以通过果园的做法来管理。新西兰的高抗Psa品种是新西兰猕猴桃育种的重点。从这个意义上说，评估亲本并选择那些表现出某种程度的Psa抗性的亲本是当前所有育种计划的一个组成部分。 ❞

「后代检验是衡量标准」

❝由于猕猴桃属植物通常雌雄同株，后代试验是育种家预测雄性果实性状育种值的唯一方法。因子杂交设计是一个标准，每个杂交种大约有25个全同胞母株进行表型评估。据报道，四倍体（4x）羊草表现出比二倍体（2x）群体更高的Psa抗性（Gea等人。2012年；Montefiori 2013年），目前已成为培育黄肉猕猴桃抗Psa的核心群体。2008年进行了一项后代试验，包括4对母本和19对四倍体父本杂交。在2010年Psa爆发之后，这些葡萄藤已经被评为疾病严重程度，这些数据是本研究的主要分析重点，以了解猕猴桃对Psa的抗性的遗传结构。 ❞

「要点」

• 对于抗性数据，在进行分析时，将其变为连续性状
• 作为连续性状是有信息损失的
• 可以用二分类性状，用GLMM模型进行遗传评估

❝不同的方法和测量尺度被用来记录Psa疾病的进展（Gea等人。2012年）在单个葡萄藤水平。对Psa的易感性没有明显的表型进展。因此，没有一种不同的疾病评分监测量表能够通过几个评分水平记录疾病的有序进展，从而使人们能够假设基础量表是连续的。定量遗传分析的目的是估计遗传参数，假设测量尺度是连续的。由于本研究中疾病评分量表的非序贯性，我们将Psa评分数据转换为一个二进制量表（0 = 无疾病1 = 疾病）。为了便于分析，二进制数据通常以样本比例表示。方差异质性和样本比例的非正态性问题传统上是通过数据转换（如角度变换）来处理的。 ❞

4. GLMM相对于LMM的优势

❝通过对家系数据进行线性混合模型（LMMs）拟合，建立了连续性状基因型随机效应方差分量估计和最佳线性无偏预测（BLUPs）。具有复杂方差结构的统计模型解释了系谱和田间布局内的空间趋势，已广泛应用于此类数据，以评估感兴趣的性状是否具有显著的遗传成分和可遗传性（Piepho等人。2008年）。然而，对于转换的比例数据，LMMs的使用可能是有限的，并且结果可能不可靠，特别是当样本大小可变且很小时。此外，在一些变换的情况下，例如角度，模型预测反变换到原始比例标度不一定在区间[0,1]内有界。经验logit和probit变换不受这个问题的影响。当估计遗传参数，如遗传力，二元性状，参数化是更好地处理在一个潜在的无限连续责任量表，其中它是最可解释的（Lee et al。2011年）。 ❞

「对于分类数据，GLMM模型评估遗传评估是标配」

❝McCullagh和Nelder（1989）的广义线性模型（GLMs）将线性模型（LM）扩展到服从非正态分布但仍属于指数分布族（如Poisson和二项式）的数据。GLM仅在数据不相关时适用。广义线性混合模型（GLMMs）扩展了这一点，通过在线性预测中加入随机效应和/或直接建模相关性来实现（Schabenberger 2005）。Piepho（1999）提供了一个很好的讨论，并举例说明了使用GLMM分析设计实验的疾病发病率数据。GLMMs最初是为指数族成员开发的，通过使用准似然估计方法，GLMMs已扩展到更广泛的应用范围（Littell等人。2006年）。因此，GLMMs是拟合二元家族数据方差分量的逻辑选择。 ❞

5. 常用软件

• R中的lme4
• SAS中的GLIMMIX
• ASReml
• R中的ASReml-R

❝本研究的主要目的是对中国四倍体中国四倍体亲本Psa后代的测试数据进行深入分析，并对Psa易感的附加遗传、环境方差成分和狭义遗传力进行可靠估计。基于上述原因，我们使用了适用于二进制/二项分布式数据的GLMM方法。GLMM的基本原理在一段时间前就已经开发出来，但它在广泛可用的统计软件中的实现却发生了很晚。R lme4包装（Bates等人。2014）首次上传于2003年，SAS®（SAS Institute Inc.2013）Proc Glimix在2008年成为V9.2中的标准程序，尽管2005年发布了生产版本。阿斯雷姆（吉尔摩等人）。2009）是一个专门的独立软件包，用于育种数据，它使用平均信息（AI）算法和稀疏矩阵方法拟合LMM。GenStat使用相同的算法进行REML估计。ASReml-R是R中ASReml的实现（Butler等人）。2009年）。虽然ASReml软件能够适应GLMMs，但ASReml-R中GLMMs的拟合显得有限。SAS®中的GLIMMIX程序因在各个学科中适合GLMMs而受到欢迎，并且有一些例子说明了它在植物和动物育种中的应用（Fikret 2011，Maxa等人）。2009年）。我们还没有遇到Proc glimix在后代测试中的应用，在这些应用中，父母的家谱信息被纳入分析。与ASReml不同，Proc GLIMMIX不是专门为分析育种数据而开发的；因此，根据可用数据和分析目标，需要进行一些调整。因此，我们的第二个目标是演示Proc-glimix在拟合模型中的应用，这些模型包含父母和后代之间的家族相似性，这是本研究中Psa子代测试数据集的例子。 ❞

6. 试验设计

• 因子设计
• 采集空间行列信息
• 收集抗性数据
• 考虑系谱信息

❝猕猴桃可以是二倍体（2x = 58）或四倍体（4x = 116）。本研究中的家系来自析因交配设计，将4对母本（分别标记为GU、GZ、GO&GT）与19对父本（编号为28:46）进行杂交。这是新西兰猕猴桃育种中使用的标准符号，雌性在一个十字架上用两个字母命名，雄性用数字命名，家族通过两个标签连接起来。有些杂交是未经试验或不成功的：雌性围棋，除了与雄性32和33杂交外；雄性32，除了与雌性围棋杂交外；以及雌性顾对雄性33。我们只为缺失的父本提供一个平衡的gla×33杂交数据集。在55个全同胞家系中，53个被保留在排除一个女性和一个男性父母的子集中。杂交所用亲本均为四倍体。田间试验在普伦蒂湾（37.8°S 176.3°E）的Te Puke研究中心进行，普伦蒂湾是新西兰主要的猕猴桃种植区。2008年春季（10月），平均每个全同胞家系约36株幼苗以随机区组设计进行了田间种植，共有三个重复；除最后一个复制品有额外的一行，可容纳8个家系的额外幼苗外，每个复制品包含三个连续行。 ❞

❝因此，实验布局并不平衡。根据种子的可用性，家系大小从24到50不等。实验区块中的10排是南北走向的，其中一排的柱子间距为4米，行间距为3米。通常在果园中种植密度高于商业试验。这里的实验单元（地块）是一排内两个柱子之间的间隔，在那里，任何一个科的12个未固定的幼苗作为一个双排种植，在网格线的两侧每隔60厘米放置6个位置（a:f）的幼苗交替。每对双排的幼苗之间有40厘米的距离。在离实验区约4米的西侧和南侧有一道防护树树篱。每排苗木下有19个间隔（编号为0:18），第一个和最后一个间隔仅包含两个种植位置（a和b）。种植的秧苗中只有一半是雌性的，但是雌性和雄性都可以评估Psa症状。在本试验中，以中国玉米螟品种‘Hort16A’和‘Zesy003’为对照地，作为行尾守卫。 ❞

「系谱信息如下：」

• 可以构建A矩阵，加入到GLMM模型中

7. 表型数据分类标准

「分类性状，分为有序的9级」

• 0，无症状
• 1~3，叶片斑点
• 4，渗出
• 5~8 伤害加重~死亡

「这种分级数据，可以很容易转化为二分类性状，比如0~3编码为0,4~8编码为1」

❝Psa症状于2010年11月首次在Te Puke的商业猕猴桃果园检测到。Psa是一种系统性疾病，随着病情的发展，葡萄藤的不同部位会出现不同的症状。这使得很难建立一个真正有序的单一植物病害严重程度量表。在最初的观察中，我们记录了葡萄藤的主要症状：叶斑的存在/缺失、藤茎末端的短暂枯萎、甘蔗枯萎、溃疡和渗出。在此基础上，我们定义了一个准序数视觉评估量表（0-8）：0 = 无症状，1，2，3 = 叶片斑点和伤害强度增加，4 = 渗出，并根据涉及的植物部分对分数进行限定：s-芽，d-主要领导和b-芽。5，6，7，8 = 藤蔓死亡并被移除。去除Psa评分高于4的植物的决定是为了减少接种量，并遵循了负责管理该病的行业机构猕猴桃藤健康组织的建议。在2010年和随后的季节，每两周使用所述的量表监测疾病通过街区的移动和进展。Gea等报道了疾病评分随时间变化趋势的初步分析。（2013年）。在接下来的遗传分析中，我们决定在某个时间点监测疾病在区块内的传播情况，但当评分仍有一些变化时，可以估计由于潜在环境因素导致的疾病评分的可变性。在此基础上，我们选择了2012年11月，即在研究果园首次发现后20个月的病害测量进行分析。测量量表的非序数性以及分数≥4触发vine移除的事实证明数据被转换为二进制：0 = 分数≤ 3，否则为1。阈值评分与我们的观察结果一致，如果症状仅限于叶斑，猕猴桃藤表现出一定程度的抗性，被认为足以避免系统性传播和藤蔓死亡。此外，这种抗性可能足以通过其他管理措施控制果园的病害。 ❞

8. 表型数据的田间分布图

• 红色的为感染
• 蓝色的为健康

9. GLMM模型

「要点：」

• Y变量为二分类性状
• 考虑了加性效应（A矩阵）
• 用了Logistic模型

10 遗传力的计算方法

「注意，残差的值为pi^2/3 = 3.218，这个值是固定的。」

11. 主要结果

❝表2列出了三对母本和18对父本的估计繁殖值，按Psa发病概率的比例。母本GZ通过向后代传递加性遗传效应而表现出最大的抗性。在抗性方面排名第二的是雄性亲本37和42（表2），它们的亲缘关系值为0.0625（通过共同祖先4，图2），但与排名靠前的雌性GZ无关。这表明，四倍体中国对虾基因型间可能存在独立的Psa抗性遗传源。在53个全同胞家系中，抗性最强的是GZ37，在所有重复试验的48株藤本中，只有3株病藤（6.25%）。这个家族对Psa抗性的预测随机效应最大，这说明了非加性基因对家族绩效的影响。基于双亲的平均eBV值和家庭随机效应的预测，该家系的Pr（Psa）估计值为0.31。

❞

12. SAS 主要代码

「遗传力计算代码：」