肖麗娜 魏瑞研 譚芷茵 張衛(wèi)華 林元震

（1.廣東省森林培育與保護利用重點實驗室/廣東省林業(yè)科學研究院，廣東 廣州 510520；2.華南農(nóng)業(yè)大學/廣東省森林植物種質(zhì)創(chuàng)新與利用重點實驗室，廣東 廣州 510642）

混合線性模型(mixed linear model)被廣泛應(yīng)用于動植物子代測試數(shù)據(jù)的分析，包括樹木在內(nèi)[1-2]?；旌暇€性模型是指含有固定效應(yīng)和隨機效應(yīng)的線性模型，來解釋目標性狀的變異程度，例如奶牛產(chǎn)奶量或林分材積增長量。育種者的主要興趣是預(yù)測某種動物或作物特定基因型在環(huán)境中的未來表現(xiàn)，他們往往將影響目標性狀的潛在遺傳因素視為待預(yù)測的隨機效應(yīng)。相比之下，試驗地點或試驗地點內(nèi)的重復(fù)一般視為固定效應(yīng)[3]。因此，混合線性模型非常適合于遺傳數(shù)據(jù)分析，在預(yù)測隨機遺傳效應(yīng)時通過固定效應(yīng)來調(diào)整表型數(shù)據(jù)。現(xiàn)在多數(shù)遺傳分析軟件采用限制性最大似然(Restricted Maximum Likelihood, REML)法來估算隨機效應(yīng)的方差分量，進而估計固定效應(yīng)和預(yù)測隨機效應(yīng)。

幾乎所有的動植物遺傳改良項目都廣泛使用子代測定試驗。在子代試驗數(shù)據(jù)分析中，育種者主要關(guān)注加性效應(yīng)。在一個給定的群體中，加性效應(yīng)是親本可以穩(wěn)定遺傳給子代的效應(yīng)，因此，將加性方差與表型方差的比例定義為性狀的狹義遺傳力(h2)。加性方差和遺傳力的數(shù)值大小對遺傳增益有著較大影響[4-5]。現(xiàn)在的子代測定數(shù)據(jù)，往往具有譜系數(shù)據(jù)、表型數(shù)據(jù)以及標記數(shù)據(jù)，專業(yè)遺傳分析軟件可基于譜系數(shù)據(jù)形成A 矩陣或基于標記數(shù)據(jù)形成G 矩陣，用于方差分量和隨機效應(yīng)BLUP 的估計[6-9]。這類軟件包括ASReml[10]、Echidna[11]、BLUPf90[12]和R 程 序 包sommer[13]、breedR[14]等。上述軟件中的免費軟件，存在功能有限或者用法復(fù)雜。因此，本文的目的是提供一個免費的R 程序包AFEchidna，并以半同胞家系測定數(shù)據(jù)為例，演示如何使用該包基于混合線性模型產(chǎn)生方差分量、遺傳參數(shù)的估計，BLUP 值與準確性的估算，隨機效應(yīng)顯著性的檢驗，以及單性狀的批量分析。

1 材料與方法

1.1 開發(fā)背景

對于植物試驗來說，尤其是林木試驗，由于試驗環(huán)境因子比較復(fù)雜，因此需要專業(yè)遺傳分析軟件可以擬合復(fù)雜的方差結(jié)構(gòu)[15-18]。目前，商業(yè)軟件ASReml 被公認為植物遺傳評估的先鋒軟件，但其費用昂貴。breedR 包雖然是針對林木試驗開發(fā)的，但其版本偏舊且方差結(jié)構(gòu)類型偏少，難以滿足林木試驗的分析需求。Echidna 是ASReml 主要研發(fā)者Gilmour 教授于2018 年開發(fā)的免費軟件，其也是使用REML 法估算參數(shù)值，功能非常接近ASReml，是目前動植物遺傳評估免費軟件中功能最強大的軟件，但其用法稍復(fù)雜，一般用戶難以適應(yīng)。因此，我們使用R 語言基于Echidna 軟件開發(fā)AFEchidna 包，降低用戶進行育種數(shù)據(jù)遺傳分析的使用難度，另外，該包不僅保留了Echidna軟件的專業(yè)遺傳評估功能，還可擴展使用R 語言強大的數(shù)據(jù)管理、統(tǒng)計分析與圖形繪制功能。

1.2 AFEchidna 包開發(fā)策略

Echidna 單機版通過.es 程序進行混合線性模型分析，并生成一系列結(jié)果文件。但其語法稍復(fù)雜，結(jié)果文件數(shù)較多，且數(shù)據(jù)處理和可視化功能較弱，一般用戶使用不太方便。因此，我們使用R 語言開發(fā)一個程序包，通過該包來輔助Echidna軟件，讓代碼運行和結(jié)果提取更直觀和更便利。

開發(fā)思路分3 個層次：第一層次是數(shù)據(jù)和es0模版程序的處理，主要是各類數(shù)據(jù)讀取和es0 模版程序生成，其通過get.es0.file()函數(shù)來完成；第二層次是混合線性模型的指定與運行，主要是將指定的混合線性模型通過R 語言傳遞給Echidna 運行，然后將生成的結(jié)果傳回到R 語言中，其通過主函數(shù)echidna()來實施，傳回的結(jié)果存為R 對象object；第三層次是模型結(jié)果的提取與運用，通過各種函數(shù)從R 對象object 中提取相應(yīng)結(jié)果，例如方差分量、模型解等，或者基于R 對象object 做進一步的遺傳參數(shù)計算或模型比較等。

圖1 展示了AFEchidna 包的運行流程。AFEchidna 包通過主函數(shù)echidna()結(jié)合模板程序.es0、表型數(shù)據(jù)、譜系甚至標記關(guān)系矩陣等文件進行混合線性模型分析，將運行結(jié)果存為R 對象object，然后通過函數(shù)Var()提取方差分量、coef()獲取模型解、predict()獲取模型預(yù)測值、pin()計算遺傳參數(shù)值、update()運行新模型和model.comp()比較不同模型。表1 列出了AFEchidna 包隸屬函數(shù)所需的R 依賴包、功能和簡易用法。

圖1 AFEchidna 包的運行流程Figure 1 The running process of AFEchidna package

表1 AFEchidna 包隸屬函數(shù)的匯總Table 1 The summary of listed functions in AFEchidna package

1.3 AFEchidna 的在線安裝

我們已將AFEchidna 包上傳至github 官網(wǎng)，并免費對外開放。感興趣的讀者可通過在R 語言 中 運 行remotes::install_github(‘yzhlinscau/AFEchidna’)來安裝程序包AFEchidna，然后運行AFEchidna::checkPack()檢測R 依賴包是否安裝，如R 依賴包未安裝，將自動進行安裝。首次使用Echidna 或AFEchidna，需要先進行軟件注冊，具體請見AFEchidna 包github 官網(wǎng)主頁(https://github.com/yzhlinscau/AFEchidna)。

1.4 示例數(shù)據(jù)

數(shù)據(jù)來源于文獻[2]中的例4.5，該數(shù)據(jù)集為來自4 個種源的某松樹36 個半同胞家系測定數(shù)據(jù)，試驗設(shè)計為完全隨機區(qū)組設(shè)計，5 個區(qū)組，每個小區(qū)中測量2～6 棵樹，測量性狀為樹高height、胸徑diameter 和材積volume。

1.5 統(tǒng)計方法

以樹高為示例性狀，假定混合線性模型如下：

狹義單株遺傳力公式如下：

其中，h2是狹義單株遺傳力，Va是加性方差，Vp是表型方差，Vf是母本方差，Vbf是區(qū)組與母本互作方差，Ve是誤差方差。

母本育種值及其準確性公式如下：

其中，bv是母本育種值，gca 是母本一般配合力，μ是總體均值；r是母本育種值準確性，SE 是母本一般配合力的標準誤差，Vf是母本方差。

上述統(tǒng)計分析將采用我們開發(fā)的AFEchidna，并與ASReml 軟件結(jié)果[2]進行比較。

2 結(jié)果與分析

2.1 方差分量與遺傳力的估計

AFEchidna 包和ASReml 軟件估計的方差分量如表2 所示，包括母本方差Vf等四個方差分量及其標準誤差的估計值均一致，具體地，區(qū)組方差Vb為0.107，母本方差Vf為0.190，區(qū)組與母本互作方差Vbf為0.198，誤差方差Ve為2.527，且各項方差值(除了區(qū)組外)均大于其標準誤差值的1.5倍，說明這兩個軟件估計的方差分量準確性高?；谏鲜龇讲罘至窟M一步估算的單株遺傳力等遺傳參數(shù)結(jié)果（表3）顯示AFEchidna 包和ASReml軟件的結(jié)果也完全一致，且估算的加性方差、表型方差與遺傳力值均大于其標準誤差值的1.5 倍。這表明AFEchidna 包估計混合線性模型的方差分量值等效于ASReml 軟件。

表2 AFEchidna 和ASReml 估計的方差分量Table 2 The estimation of variance components from AFEchidna and ASReml

表3 AFEchidna 和ASReml 估算的遺傳參數(shù)Table 3 The estimation of genetic parameters from AFEchidna and ASReml

2.2 固定效應(yīng)估計

本研究中，固定效應(yīng)包含總體均值和種源Prov，它們的無偏估計值如表4 所示，在AFEchidna 中，將總體均值固定為0，直接給出種源各個水平的效應(yīng)值；在ASReml 中，將種源首個水平的效應(yīng)值固定為0，然后求解出總體均值和種源其它水平值。由此可見，AFEchidna 和ASReml 對于固定效應(yīng)的求解方法稍有不同，結(jié)果也稍有差異。不過，ASReml 得到的種源效應(yīng)值是相對值，當它們加上總體均值后，與AFEchidna 的結(jié)果基本一致。

表4 固定效應(yīng)的估計Table 4 The estimate of fixed effects

2.3 母本育種值及其準確性

如表5 所示，羅列了前6 個母本的一般配合力、育種值及其準確性，可以看出AFEchidna 和ASReml 得到的母本一般配合力及其標準誤差結(jié)果是一致的，總體均值稍有差異，但其標準誤差也是相同的。由于總體均值的輕微不同，導(dǎo)致估計的母本育種值在兩軟件間也稍有不同，但它們準確性結(jié)果是一致的，不過兩軟件估計育種值之間的相關(guān)達到了0.999(P＜.001, 圖2)，因此，可以認為AFEchidna 和ASReml 估計的母本育種值也是等效的。

表5 估計的前6 個母本育種值及其準確性Table 5 The first six female breeding values and their accuracy

圖2 來自ASReml 和AFEchidna 估計的母本育種值的散點圖Figure 2 Scatter plot of female breeding values from ASReml and AFEchidna

2.4 隨機效應(yīng)的顯著性檢驗

將本研究的混合線性模型視為全階模型，其隨機效應(yīng)含有區(qū)組效應(yīng)Block、母本效應(yīng)Female以及區(qū)組與母本互作效應(yīng)BlockFemale，依次刪除它們中的一項作為降階模型，然后采用LRT 法來檢驗相應(yīng)隨機效應(yīng)的顯著性[2]。AFEchidna 包中的model.comp()函數(shù)可以進行LRT 檢驗。LRT 檢驗結(jié)果(表6)顯示，區(qū)組效應(yīng)Block、母本效應(yīng)Female 以及區(qū)組與母本互作效應(yīng)BlockFemale 對樹高的影響均為極顯著(P＜0.01)。

表6 隨機效應(yīng)的顯著性檢驗Table 6 Significant test for random effects

2.5 性狀批量分析

所謂批量分析即混合線性模型中，固定效應(yīng)、隨機效應(yīng)和殘差效應(yīng)均相同條件下，將待分析的多個性狀依次傳入程序中進行分析并輸出結(jié)果，其可是單性狀的批量運行，也可是多性狀的批量運行。尤其是待分析的性狀數(shù)較多時，通過批量分析可大大減少分析工作量。AFEchidna 包中的主函數(shù)echidna()可通過參數(shù)batch=TRUE 來實現(xiàn)性狀批量分析。本研究中，采用同樣的混合線性模型，對樹高height、胸徑diameter 和材積volume做了單性狀的批量分析，得到的方差分量結(jié)果如表7 所示。

表7 基于單性狀批量分析的方差分量估計值Table 7 The estimation of variance components from single trait batch analysis

3 討論與結(jié)論

本研究利用R 語言4.0.2 基于Echidna 單機版(V1.54)開發(fā)了R 程序包AFEchidna，以半同胞數(shù)據(jù)為例，使用混合線性模型進行了隨機效應(yīng)方差分量與遺傳力的估計、固定效應(yīng)和育種值及其準確性的估算，這些結(jié)果與商業(yè)軟件ASReml 的結(jié)果高度一致。此外，還演示了隨機效應(yīng)顯著性的檢驗和單性狀的批量分析。

王德源和童春發(fā)[19]開發(fā)了R 程序包HalfsibMS，用于分析林木多地點半同胞子代測定數(shù)據(jù)，采用他們指定的數(shù)據(jù)變量格式直接錄入數(shù)據(jù)進行遺傳分析。這種做法雖然看似方便，但實則應(yīng)用范圍很窄。由于植物(包含林木)的遺傳測定類型廣，不僅田間試驗設(shè)計類型繁多，而且遺傳材料來源復(fù)雜，簡單的一個固定混合模型，根本不足以囊括，更別說涉及不同的統(tǒng)計方法，比如空間分析模型[15]、因子分析模型[16,18]和基因組選擇模型[20]。因此，主流的遺傳評估軟件(例如ASReml、Echidna、SAS、breedR 等)均未采用固定的混合模型?；谶@種理念，我們開發(fā)AFEchidna包時，在主函數(shù)echidna()中通過參數(shù)fixed、random 和residual 保留了混合線性模型可以指定任意試驗因子及其方差結(jié)構(gòu)的靈活性，延續(xù)了Echidna分析混合模型的強大功能。

Isik 等[2]指出商業(yè)軟件ASReml 采用平均信息(average information, AI)和稀疏矩陣(sparse matrix)算法求解大量的混合模型方程組，其比SAS Proc Mixed 模版(采用牛頓-拉夫遜代數(shù)法，Newton-Raphson algorithm)速度更快、效率更高。ASReml 由于可擬合多樣化的方差結(jié)構(gòu)，其可輕松處理不同交配設(shè)計、不同田間設(shè)計、多變量模型等分析。ASReml 的不足之處是它缺乏數(shù)據(jù)管理、綜合性統(tǒng)計和圖形可視化功能。Echidna 作為ASReml的免費姊妹版軟件，同樣繼承了ASReml 的優(yōu)缺點，因此，我們借助R 語言開發(fā)AFEchidna 包，加入R 語言在數(shù)據(jù)管理、綜合性統(tǒng)計和圖形可視化方面的優(yōu)勢。例如，LRT 檢驗法，用于檢驗?zāi)Ｐ碗S機項的顯著性，或者比較不同模型，我們根據(jù)LRT 檢驗原理編寫了函數(shù)model.comp()。再比如，隨著育種進程的推進與目標性狀數(shù)量的累積，對性狀的批量分析需求是必然的，Echidna 軟件雖可以處理單性狀的批量分析(麻煩是結(jié)果文件偏多，分析結(jié)果批量提取難度大)，但不能處理多性狀模型的批量分析。于是，我們在主函數(shù)echidna()中通過參數(shù)batch 來運行單性狀或多性狀的批量分析，并可直接輸出結(jié)果。此外，主函數(shù)echidna()中還有參數(shù)batch.G 和batch.R，可用于多種G 結(jié)構(gòu)和多種R 結(jié)構(gòu)的批量分析，這是我們首次提出的方法，在其它遺傳評估軟件中均未涉及。

簡而言之，AFEchidna 包不僅繼承了Echidna軟件的優(yōu)點，可擬合不同交配設(shè)計(測交、巢式交配、雙列雜交等)模型、空間分析模型、多變量模型、多地點模型和基因組BLUP 模型，還延續(xù)了R 語言的優(yōu)勢，可自編函數(shù)批量分析性狀、檢驗隨機項的顯著性等，并可利用R 語言管理數(shù)據(jù)、綜合統(tǒng)計和圖形可視化的功能。未來我們會陸續(xù)介紹AFEchidna 包如何擬合空間分析模型、多變量模型、多地點模型和基因組BLUP 模型。