基于R语言的亚麻种质资源农艺性状相关及聚类分析

张丽丽 耿立格 孙娟 王玉祥 曲志华 李世芳 乔海明

摘要：本研究旨在通过对河北省种质资源库部分亚麻资源进行进一步的表型鉴定，以期探索亚麻种质资源的农艺性状相关性和聚类研究，用于亚麻新品种选育研究工作。本试验以499份亚麻种质资源为材料，应用R语言对亚麻主要农艺性状进行相关性分析，得出各性状间的相关系数，并进行聚类分析。研究结果显示亚麻单株粒重与分茎数、主茎分枝数、单株有效果数存在显著正相关，株高与工艺长显著正相关。聚类分析将499份资源分成了3类，种群1更有利于亚麻高产育种，种群2有利于抗倒伏亚麻新品种的选育。

关键词：亚麻;农艺性状;相关性;聚类分析;R语言

中图分类号：S563.2

文献标志码：A

论文编号：cjas20190700124

0引言

亚麻（Linum usitatissimum L.）是亚麻科亚麻属一年或多年生草本植物[1]，按用途可分为油用、纤用和油纤兼用3种类型[2]。河北省是全国六大油用亚麻产区之一，年种植面积在3.7万hm2左右。亚麻籽富含α一亚麻酸、膳食纤维及木酚素等保健成分[3]，能降低三高、抗肿瘤、抗衰老、预防老年痴呆、增加智力和保护视力等[4]。种质资源是农业科学研究尤其是育种研究不可缺少的重要物质基础，通过对亚麻种质资源进行鉴定评价，可挖掘种质资源潜力，拓宽遗传多样性，有利于突破亚麻育种瓶颈[5]。

“河北省农业生物资源保存中心”始建于1983年，2005年初步建立河北省农作物种质资源特性评价鉴定信息系统[6]。拥有种子低温保存长期库、中期库、短期库以及试管苗库、超低温库等保存设施，共收集、保存了57种作物45962份种质资源，开展了花生口[7]、大豆[8]、黑豆[9]、玉米[10]、小麦[11]等作物资源的相关研究。张家口市农业科学院在“八五”期间为该中心登记入库亚麻种质资源1800多份[12]。这些资源在种子库己保存了30年之久，当时国内的亚麻调查记载标准并不统一，本研究旨在将其中部分资源进行繁种更新，并根据《亚麻种质资源描述规范和数据标准》[13]对其进行更规范化的表型鉴定，以期通过探索亚麻种质资源的农艺性状相关性和聚类分析研究，为亚麻新品种选育的亲本选择提供数据支持，用于辅助亚麻新品种选育研究工作。

随着数据量级的不断增大，大数据挖掘算法提出了新的需求与挑战。本研究应用数据聚类的K-Means算法，结合R语言的实现，通过迭代方法实现了基于Map-Reduce函数的K-means优化算法，并用R程序得以实现[14-15];最后给出了程序的测试与应用，验证算法的可行性。

1材料与方法

1.1试验材料

“河北省农业生物资源保存中心”提供的500份亚麻种质资源（见表1），2018年在河北省张家口市农业科学院的坝上试验基地繁种更新，每行长6.67m，行距0.25m，4行区种植。其中1份由于蒴果开裂，果实脱落未能及时收获，共收获499份材料。

1.2试验方法

试验于2018年在河北省张家口市农业科学院张北基地进行。每小区1mx6.67m，4行区种植，田间管理同一般大田。亚麻成熟收获后，每小区随机收获10株用于考种，并测定小区产量，记载按照《亚麻种质资源描述规范和数据标准》[13]进行。

1.3数据分析

将499份亚麻资源的8个农艺性状指标（株高、工艺长、分茎数、主茎分枝数、单株有效果数、单株无效果数、每果粒数和单株粒重）作为细分变量，基于该数据，采用k-means聚类分析方法，将具有相似属性的亚麻品种聚为一类，使得同一类品种具有高度的相似性。采用excel计算各农艺性状平均数，用R语言进行各农艺性状的相关性分析（Pearson相关系数、Kendall相关系数、Spearman相关系数）和聚类分析。

2结果与分析

2.1亚麻种质资源各主要农艺性状的相关性分析

图1为亚麻种质资源各主要农艺性状指标的相关系数矩阵图，主对角线为各指标直方图;主对角线下方为各指标间的散点图;主对角线上方从上到下依次为各指标间的Spearman相关系数、Kendall相关系数、Pearson相关系数以及Pearson相关系数的显著性检验的P值[16]，背景颜色越接近蓝色，则相关系数越接近于1，背景颜色越接近红色，则相关系数越接近于-1。由图可看出，亚麻种质资源的主茎分枝数、单株粒重均与其他7个农艺性状存在正相关关系，其他各农艺性状间存在正相关也存在负相关。其中，株高与工艺长显著正相关，Spearman相关系数最大，达到0.87，与主茎分枝数和单株无效果数存在弱相关关系;单株有效果数与单株粒重、主茎分枝数、分茎数显著正相关，Spearman相关系数分别为0.73、0.5、0.52;主茎分枝数与单株粒重、单株有效果数均为显著正相关，相关系数分别为0.54、0.50;单株有效果数与分茎数显著正相关，相关系数0.52;主茎分枝数与各性状均呈正相关，其中与分茎数和工艺长的相关性没有达到显著水平，与其他5个性状均达到显著水平，与单株粒重、单株有效果数的Spearman相关系数分别为0.54、0.50;单株无效果数除与单株粒重相关系数为0外，与其他各性状相关显著，其中与株高、工艺长和每果粒数显著负相关，与分茎数、主茎分枝数和单株有效果数显著正相关;每果粒数与单株粒重存在弱相关关系。Pearson相关系数的显著性检验的P值可用于比较不同类别数值时的聚类结果，从而找出最优聚类结果，该值越大表明组内差距越小，组间差距越大，聚类效果越好。

2.2亚麻种质资源的K-means聚类分析

应用R语言通过计算不同K值下簇集中各对象的轮廓系数确定最优聚类数[17];然后通过凝聚层次聚类的方法获得数据集的分布，确定不同类别的中心坐标点;最后利用k-means方法完成聚類，将这499份资源聚成3类。

2.2.1轮廓系数轮廓系数是聚类效果好坏的一种评价方式，可以更好实现对于聚类效果的判断[18]。将各指标数据进行标准化处理后，根据轮廓系数图（图2）可以看到，在聚类数为3时轮廓系数达到了峰值，所以最佳聚类数为3。

2.2.2各类别的中心点坐标中心坐标点是用于比较不同类别数值时的聚类结果，从而找出最优聚类结果，该值越大表明组内差距越小，组间差距越大，聚类效果越好。由表2可见，3个分群所包含的样本量分别为141、215、143;各分群的组内平方和分别为15353.26、20657.56、15937.71，分群3最高;另外组间平方和占总平方和的53.2%。

2.2.3各分群的概率函数由各亚麻种质资源分群的概率密度函数图进行分群特点分析如下，详情见图3，图4，图5。

分群1的特点：株高主要分布在50-60cm之间，工艺长主要分布在25-35cm之间，分茎数主要集中在0.5-1.5个之间，主茎分枝数主要集中在4-6个之间，单株有效果数主要在25-35个之间，单株无效果数在0-2个之间，每果粒数在6-9个之间，单株粒重在0.8-1.3g之间。

分群2的特点：株高主要分布在35-50cm之间，工艺长主要分布在18-30cm之间，分茎数主要集中在0.5-1.5个之间，主茎分枝数主要集中在2.5-5个之间，单株有效果数在10-25个之间，单株无效果数在0-2个之间，每果粒数在5-9个之间，单株粒重在0.3-0.8g之间。

分群3的特点：株高主要分布在55-65cm之间，工艺长主要分布在30-45cm之间，分茎数主要集中在0-1.5个之间，主茎分枝数主要集中在2.5-5个之间，单株有效果数在10-25个之间，单株无效果数在0-1个之间，每果粒数在5-9个之间，单株粒重在0.3-0.8g之间。

3结论与讨论

人们对农作物种质资源尤其亚麻种质资源的研究结果分析采用的都是SPSS、DPS或者SAS等软件[2，19]。目前，R语言分析方法被越来越多的农业科研人员接受。温岚等[20]应用R语言对长蒴黄麻5个产量性状进行回归与相关分析;肖海霞等[21]采用R语言对吐鲁番驴、疆岳驴及和田青驴的体重和体尺性状进行了相关和回归分析;张祯勇等[22]使用R语言对”3414”肥料效应试验结果拟合了二元二次、三元二次肥料效应方程;盛坤等[23]用R语言计算冬小麦品种品质性状的安全指数;郭敏杰等[24]用R语言对花生区试进行品种的适应性、丰产性和稳产性，试点环境的相关性、区分力和代表性分析。

本研究通过对499份亚麻种质资源的8个主要农艺性状指标按照标准进行了规范化鉴定，并应用R语言分析方法进行相关性分析，得出亚麻种质资源的主茎分枝数、单株粒重均与其他7个农艺性状存在正相关关系，其他6个农艺性状间存在正相关也存在负相关，其中单株粒重和分茎数、主茎分枝数、单株有效果数显著正相关，株高和工艺长显著正相关。从主要农艺性状的相关性及相关系数分析结果看，分茎数和主茎分枝数指标会对单株有效果数指标产生显著影响，而单株有效果数和主茎分枝数指标会对单株粒重指标产生显著影响。因此，在亚麻育种亲本选择时，注重对单株有效果数和单株粒重的选择，从而能提高分茎数、主茎分枝数，提高单株生产力及种子产量，从而达到高产育种的目的，这为亚麻新品种选育工作的亲本选择提供了新的理论依据。

应用R语言将这499份资源聚成了3类，对比3个种群聚类分析结果发现，种群1在单株有效果数、单株粒重等性状上表现更为优良，株高、工艺长等性状表现较差，容易发生倒伏。种群2在株高、工艺长等性状上表现较好，不易发生倒伏，其他性状表现一般。种群3在株高、工艺长等性状上表现较差，容易发生倒伏，其他性状表现和种群2相似。整体上，种群1各农艺性状表现更为优良，更有利于亚麻高产育种，种群2有利于抗倒伏亚麻新品种的选育。此结论为亚麻新品种选育提供了新的目标亲本材料。

参考文献

[1]米君.亚麻（胡麻）高产栽培技术[M].北京：金盾出版社，2006.

[2]崔翠，周清元，王利鹃，等.亚麻种质主要农艺性状主成分分析与综合评价[J].西南大学学报：自然科学版，2016，38（12）：10-18.

[3]党占海，赵玮.胡麻产业技术体系[M].兰州：兰州大学出版社，2015

[4]郭永利，范丽娟.亚麻籽的保健功效和药用价值[J].中国麻业科学，2007，29（3）：147-149.

[5]党占海，赵玮中国现代农业产业可持续发展战略研究胡麻分册[M].北京：中国农业出版社，2016.

[6]耿立格，李灵芝，王丽娜，等.河北省农作物种质资源特性评价鉴定信息系统的建立[J].河北农业科学，2005，9（02）：70-72.

[7]刘立峰，耿立格，王静华，等.河北省花生地方品種农艺性状和品质性状的遗传分化[J].植物遗传资源学报，2008，9（02）：190-194.

[8]耿立格，宋春风，王丽娜，等.近红外光谱无损测定大豆种子生活力方法研究[J].植物遗传资源学报，2013，14（06）：1208-1212.

[9]耿立格，王丽娜，张磊，等.河北省绿子叶黑豆种质资源表现型和ISSR标记遗传多样性分析[J].植物遗传资源学报，2010.11（03）：266-270.

[10]张磊，耿立格，王丽娜，等.不同玉米自交系萌芽期的抗旱性研究[J].玉米科学，2010，18（04）：77-81.

[11]许红星，许云峰，耿立格，等.我国小麦农家品种和近缘种对白粉病的苗期抗性[J].中国生态农业学报，2011，19（05）：1210-1214.

[12]米君.河北省胡麻生产调研报告[J].现代农村科技，2009（20）：49-50.

[13]王玉富，粟建光.亚麻种质资源描述规范和数据标准[M].中国农业出版社.2006.

[14]郭显娥.K-Means优化算法的R语言实现[J].山西大同大学学报：自然科学版，2018，34（2）：27-29，33.

[15]李晓瑜，俞丽颖，雷航，等.一种K-means改进算法的并行化实现与应用[J].电子科技大学学报，2017，43（1）：61-68.

[16]金林，李研.几种相关系数辨析及其在R语言中的实现[J].统计与信息论坛，2019，34（4）：3-11.

[17]夏士雄，李文超，周勇，等.一种改进的k-means聚类算法（英文）[J].Joumal of Southeast University （English Edition），2007（03）：435-438.

[18]朱连江，马炳先，赵学泉.基于轮廓系数的聚类有效性分析[J].计算机应用，2010，30（S2）：139-141，198.

[19]张丽丽，刘晶晶，乔海明，等.从俄罗斯引进亚麻种质资源的农艺性状评价[J].中国油料作物学报，2017，39（05）：698-703.

[20]温岚，陈基权，戴志刚，等.长蒴黄麻产叶量的多元回归与偏相关的R语言分析[J].作物杂志，2013（01）：49-53.

[21]肖海霞，托乎提·阿及德，石国庆，等.基于R语言的吐鲁番驴体尺和体质量相关分析[J].河南农业科学，2012，41（10）：153-157.

[22]张祯勇，高明文，肖启银，等.基于R语言的“3414”肥效试验的统计分析[J].中国农学通报，2011，27（27）： 127-134.

[23]盛坤，李晓航，王映红，等.用R语言计算冬小麦品种品质性状的安伞指数[J].中国农学通报，2017，33（25）：8-12.

[24]郭敏杰，邓丽，任丽，等基于R语言的AMMI和GGE双标图在花生区试中的应用[J].花生学报，2017，46（02）：24-31.

[25]张丽丽，米君，李世芳.胡麻种间杂交种主要农艺性状与产量的关系研究[J]。河北农业科学，2014，18（03）： 76-78，88.

基金项目：现代农业产业技术体系建设专项资金资助项目“胡麻抗逆育种岗位”（CARS-14-1-08）;国家科技资源共享服务平台“国家农作物种质资源共享服务平台”（NICGR2018-23）;河北省现代农业产业技术体系油料产业创新团队“特色油料岗位”（HBCT2018090204）。

第一作者简介：张丽丽，女，1983年出生，河北保定人，副研究员，硕士研究生，主要从事胡麻新品种选育及栽培技术研究。通信地址：075000張家口市经开区惠通街张家口市农业科学院，Tel：0313-7155779，E-mail：zhanglili57@126.com。

通讯作者：乔海明，男，1965年出生，张北人，研究员，本科，主要从事胡麻新品种选育及栽培技术研究。通信地址：075000张家口市经开区惠通街张家口市农业科学院，Tel： 0313-7155774，E-mail：qhm1965@163.com。

收稿日期：2019-07-16，修回日期：2019-08-30。