基于芯片數據的長白豬繁殖性狀基因組選擇研究

摘 要： 旨在比較不同基因組預測模型預測準確性與運行效率，以探究支持向量機（SVM）回歸與隨機森林（RandomForest）回歸在基因組預測中的應用價值與應用前景。本研究使用博瑞迪豬50K液相芯片，采用GBLUP、BayesB、BayesLASSO、SVM回歸和RandomForest回歸等基因組預測模型，對1 001頭長白豬繁殖性狀進行基因組預測評估。研究發現，在總產仔數、產活仔數、窩重等繁殖性狀中，使用SVM回歸徑向基函數核的評估準確性均最高；產活仔數、窩重在參數C值為1時評估準確性達到最大值，總產仔數在參數C值為2時評估準確性達到最大值。在總產仔數、產活仔數、窩重等繁殖性狀中，使用RandomForest回歸評估ntree、mtry、nodesize等參數時發現，基因組預測準確性隨著參數的變化展現一定的隨機性。RandomForest回歸模型在總產仔數、產活仔數、窩重中的評估準確性均最高，其次為SVM回歸，GBLUP、BayesB、BayesLASSO等模型遺傳評估準確性較差且保持一致。交叉驗證相關性顯示，不同模型遺傳評估結果存在較強的相關性，為0.806～0.995。SVM回歸與RandomForest回歸等非參數機器學習模型在豬繁殖性狀基因組選擇中具有一定的優勢，但運行時間在一定程度上限制了這些算法的使用。隨著算法的研究優化，SVM回歸與RandomForest回歸等非參數機器學習模型將具有良好的應用前景。

關鍵詞： 長白豬；繁殖性狀；基因芯片；基因組選擇；機器學習

中圖分類號：S828.3"""" 文獻標志碼：A"""" 文章編號： 0366-6964（2025）01-0213-09

收稿日期：2024-05-31

基金項目：福建省種業企業培優項目（2120814-農業生產發展支出）

作者簡介：陽文攀（1994-），男，湖北孝感人，碩士，主要從事動物遺傳育種與繁殖研究，E-mail：945226087@qq.com

*通信作者：丁能水，主要從事種豬遺傳及育種新技術、生豬養殖新技術方向的研究與開發工作，E-mail：13631698@qq.com

YANG" Wenpan 大量SNPs標記被廣泛應用于畜禽分子育種中^［1-3］。基因組選擇技術通過覆蓋全基因組的高密度SNPs標記信息對育種值進行估計，不同統計模型與計算方法是影響基因組選擇準確性與效率的主要因素之一^［4-7］。GBLUP（genomic best linear unbiased prediction）與Bayes回歸模型等參數方法是應用于基因組選擇的常用方法^［8］，非加性效應參數的增加將導致巨大的計算量^［9］，因此這些線性模型通常只考慮加性效應，而忽略標記與表型間復雜的非加性關聯^［10］。近年來，有較多將機器學習方法作為復雜性狀基因組遺傳評估參數模型替代方法的研究^［11-14］，機器學習能以適應性的方式捕捉基因型和表型之間的隱藏關系，且較少對性狀潛在的遺傳結構進行假設^{［10，15］}。支持向量機（support vector machines，SVM）是一種著名的機器學習算法^［16-18］，它使用多個特征向量通過在兩個類之間創建決策邊界來完成預測^［19］，基于不同的核方法，SVM還可以在一定程度上處理表型和基因組之間的非線性關系。隨機森林（RandomForest）采用引導聚合抽樣^［20］，可捕獲復雜的交互，并且對于過度擬合數據具有一定的穩健性^［21］。

本研究擬通過GBLUP、BayesB、BayesLASSO等參數模型與SVM回歸、RandomForest回歸等非參數機器學習模型來評估長白豬繁殖性狀基因組選擇的準確性與評估效率，以探討機器學習在豬遺傳評估中的應用價值和前景。

1 材料與方法

1.1 試驗群體

本試驗所用數據來自福建傲農集團2019年4月至2024年3月1 914頭長白種豬的6 195條繁殖記錄，包括總產仔數、產活仔數、窩重等繁殖性狀。使用磁珠法組織基因組提取試劑盒從1 001頭豬組織中提取基因組DNA，檢測DNA濃度及瓊脂糖凝膠電泳符合分型標準后，在石家莊博瑞迪生物技術有限公司進行豬50K液相芯片分型。

1.2 基因型填充及質控

為避免缺失基因型對基因組選擇計算產生影響，本試驗使用Beagle v5.4^［22］對原始芯片數據進行基因型填充，使用SNPRelate v1.36.1^［23］對芯片填充數據進行質控，選擇次等位基因頻率（minor allele frequency，MAF）≥ 0.05、哈迪-溫伯格平衡（Hardy-Weinberg equilibrium，HWE）檢驗Pgt;1×10^-6的SNPs進行基因組選擇計算。

1.3 表型數據處理

為保證數據質量，本試驗剔除了總產仔數小于3的繁殖數據。考慮后續分析的便利性，使用BLUPF90通過單性狀重復力模型計算固定效應值，去除個體相應固定效應值后取個體表型平均值作為全基因組選擇表型值。該單性狀重復力模型如下所示：

y=Xb+Za+Wp+e

其中：y為表型值向量；b為固定效應向量，包含場-年-季、胎次；a為個體加性遺傳效應向量，且a～N（0，Hσ2 a），H為基因型與系譜親緣關系矩陣^［6］；p為永久環境效應向量，且p～N（0，Iσ2 p），I為單位矩陣；e為隨機殘差向量，且e～N（0，Iσ2 e）；X、Z、W為b、a、p的關聯矩陣。

1.4 基因組預測模型

1.4.1 參數預測模型" 對于GBLUP，模型描述如下所示：

y^*=μ+Zg+e

其中，y^*為校正后表型，即y^*=mean（y-Xb）；μ為表型均值向量；g為個體加性遺傳效應向量，且a～N（0，Gσ2 a），G為基因型親緣關系矩陣^［5］；Z為g的關聯矩陣；e為隨機殘差向量，且e～N（0，Iσ2 e）。

對于BayesB^［4］與BayesLASSO^［24］，模型描述如下所示：

y* i=μ+∑sj=1m ijα j+e i，

其中，y* i為個體i校正后表型；μ為表型均值向量；m ij為第i個個體在第j個SNP處觀察到的基因型協變量；α j是與第j個SNP相關的等位基因替代效應；e i為第i個個體隨機殘差，且e～N（0，Iσ2 e）。對于BayesB，SNP α j服從正態分布，不同α j方差不同為σ2 α j，p（σ2 α j）=x^-2（ν，S），即SNP效應方差σ2 α j服從自由度為ν，參數尺度為S的逆卡方分布。σ2 α j先驗分布假設為σ2 α j=0的概率為π，p（σ2 α j）=x^-2（ν，S）的概率為1-π。通過最小化殘差平方和及約束回歸系數絕對值的和獲得BayesLASSO回歸系數估計值，通常采用的先驗分布為雙指數分布（又稱拉普拉斯分布）：p（g）=q j=1λ2exp（-λα j），它是一個兩水平的層級模型分布，由p（α jσ2 α j）=N（0，σ2 α j）與p（σ2 α j）=gamma（ λ2/2）混合組成。GBLUP、BayesB、BayesLASSO等模型使用BGLR v1.1.1^［25］計算。

1.4.2 SVM回歸

本試驗使用ε-支持向量回歸，為了執行非線性回歸，數據被核函數映射到更高維的空間中，模型為：

y=β 0+f x（X|β）+e=β 0+Kx，xT+e

其中Kx，xT是n×n核矩陣，β是n×1向量（未知）。有許多不同的核，它們被定義為徑向基函數核（Radial Basis Function Kernel，或高斯核）：

K ij（x i，xT i）=exp［-γ（x i-x j）（x i-x j）T］

多項式核（polynomial kernel）：

K ijx i，xT i=γx ixT j+rd

線性核（linear kernel）：

K ij（x i，xT i）=x ixT j

Sigmoid核（Sigmoid kernel）：

K ij（x i，xT i）=tanh（γx ixT j+r）

在求解SVM的過程中，最終會轉化為一個最佳化問題：

min12ω+C∑ni=1ε i

約束為：y i（ωTx i+b）≥1-ε i；ε i≥0，i= …，n。

其中ω是要求解的超平面，ε i是第i個樣本點的回歸損失，C是懲罰系數，即誤差的容忍度。γ是RBF核函數的一個參數。本試驗使用e1071 v1.7-14^［26］R包svm函數進行評估。

為篩選SVM最佳參數，在e1071包svm函數默認參數值附近進行參數調試，即C值為1、γ值為1/N snp（N snp為用于分析的SNP數）。本試驗使用svm函數默認C值與γ值分別對徑向基函數核、多項式核、線性核、Sigmoid核進行1×10重交叉驗證，篩選準確性最高的核函數；使用最優核函數與默認γ值，其次將C值設置為0.1、0.5、1、2、3、4、5、6、7、8、9、10進行1×10重交叉驗證，篩選準確性最高的C值；使用最優核函數與默認C值，其次將RBF核函數的γ值設置為10^-1、10^-2、10^-3、10^-4、10^-5、10^-6、10^-7、10^-8、10^-9進行1×10重交叉驗證，篩選準確性最高的γ值。

1.4.3 RandomForest回歸

RandomForest回歸模型通過隨機抽取樣本和特征，建立多棵相互不關聯的決策樹，通過并行的方式獲得預測結果。每棵決策樹都能通過抽取的樣本和特征得出一個預測結果，通過綜合所有樹的結果取平均值，得到整個森林的回歸預測結果。

每個樹節點中，從s個預測變量（其中s代表SNPs的數量）中隨機抽取mtry變量，使用給定的損失函數作為標準選擇，選擇風險最低的預測變量和分裂閾值。在完成所有獨立樹之后，RandomForest回歸聚合各樹信息以計算最終預測，其計算公式如下：

y^*=1ntree∑ntreeb=1ψ b（y^*，M）

其中，ψ b為獨立RandomForest回歸數，包含抽樣樣本、每個樹節點上的預測變量（標記）、分裂閾值及終端節點值。對于未觀察到的值，通過在每棵樹的流程圖中傳遞預測變量來獲得預測，并將終端節點的相應估計分配為預測值。將RandomForest回歸中每棵樹的預測值平均以計算未觀察數據的最終預測。本試驗使用randomForest v4.7-1.1包randomForest函數進行評估^［27］。

為篩選RandomForest最佳參數，將GIANOLA^［20］使用參數作為默認參數，在默認參數附近進行參數調試，即ntree為1 000、mtry為N snp、nodesize為5。使用推薦mtry值與nodesize值，分別將生長樹數目ntree設置為800、1 000、1 200、1 400、1 600、1 800、2 000進行1×10重交叉驗證，篩選最佳ntree值；使用推薦ntree值與nodesize值，分別將每個樹節點數mtry設置為150、200、250、300、350、400進行1×10重交叉驗證，篩選最佳mtry值；使用推薦ntree值與mtry值，其次將每個樹終端最小變量數nodesize設置為4、5、6、7、8、9、10進行1×10重交叉驗證，篩選最佳nodesize值。

1.5 基因組預測

使用GBLUP、BayesB、BayesLASSO、SVM、RandomForest等模型進行基因組預測，通過5×10交叉驗證獲取各模型預測準確性與運行時間，并比較各模型準確性相關性值。預測準確性為預測值與真實值的相關系數。

2 結 果

2.1 基因型數據處理與遺傳力估計

本研究使用博瑞迪液相芯片對1 001頭長白母豬進行基因分型，共獲得52 000個SNPs，使用Beagle v5.4進行基因型填充后共獲得50 254個SNPs。使用SNPRelate v1.36.1對芯片填充數據進行質控，6 354個SNPs因MAFlt;0.05而剔除，815個SNPs因HWE檢驗Plt;1×10^-6而剔除，最終43 085個SNPs被用于后續分析。

通過BLUPF90對長白豬繁殖性狀固定效應值進行計算，同時獲取長白豬繁殖性狀遺傳參數估計值如表1所示。其中總產仔數遺傳力與重復力最高，窩重遺傳力最低，產活仔數重復力最低。

2.2 SVM參數對基因組預測準確性的影響

為篩選SVM最佳參數，使用e1071包對核函數、C值分別進行1×10重交叉驗證，其評估準確性與計算時間如圖1所示。使用SVM函數默認參數，在總產仔數、產活仔數、窩重等繁殖性狀中，使用高斯函數核評估準確性均最高，線性核評估準確性均最低。而在運行時間上，線性核運行時間較長。在產活仔數、窩重等繁殖性狀中，評估準確性在C值為1時達到最大值，在C值為5后基本不變?？偖a仔數評估準確性在C值為2時達到最大值，在C值為6后基本不變。運行時間受C值影響較小。在總產仔數、產活仔數、窩重等繁殖性狀中，將高斯核函數的γ值設置為10^-1、10^-2、10^-3、10^-4、10^-5、10^-6、10^-7、10^-8、10^-9進行1×10重交叉驗證，發現評估準確性均保持不變。

2.3 RandomForest參數對基因組預測準確性的影響

為篩選RandomForest最佳參數，使用RandomForest包對ntree、mtry、nodesize等參數進行1×10重交叉驗證，其評估準確性與計算時間如圖2所示。在總產仔數、產活仔數、窩重等繁殖性狀中對ntree、mtry、nodesize等參數進行評估時，發現基因組預測準確性隨著參數的變化展現一定的隨機性。在運行時間上，隨著nodesize值的增加，運行時間逐漸減少，而隨著ntree、mtry值的增加，運行時間逐漸增加。

2.4 各模型計算準確性與運行時間

在對SVM與RandomForest參數進行評估后，使用SVM高斯函數核與軟件默認參數與RandomForest軟件默認參數設置模型參數進行基因組預測。GBLUP、BayesB、BayesLASSO、SVM、RandomForest等模型計算準確性與運行時間如圖3所示。RandomForest模型在總產仔數、產活仔數、窩重上均表現出最高的評估準確性，其次為SVM，而GBLUP、BayesB、BayesLASSO等遺傳評估準確性較差且保持一致。在運行時間上，GBLUP運行時間保持最優，其次為RandomForest模型，BayesB、BayesLASSO、SVM等模型評估時間均較久。

2.5 各模型相關性分析

為驗證各模型在交叉驗證時的評估一致性，將各模型評估準確性進行Person相關性分析，具體結果如圖4所示。其中GBLUP、BayesB、BayesLASSO間評估相關性最高，為0.928～0.995；SVM與GBLUP、BayesB、BayesLASSO評估相關性較高，為0.891～0.922；RandomForest與GBLUP、BayesB、BayesLASSO、SVM評估相關性較低，為0.806～0.871。

3 討 論

本研究使用BLUPF90通過一步法對長白豬繁殖性狀遺傳參數進行估計，得到的總產仔數、產活仔數、窩重的遺傳力分別為0.102±0.019、0.084±0.017、0.071±0.016，總產仔數、產活仔數、窩重的重復力分別為0.171±0.015、0.138±0.015、0.147±0.014，這與前人的研究類似^［28-29］。

在篩選SVM最佳參數時，在總產仔數、產活仔數、窩重等繁殖性狀中，使用高斯函數核評估準確性均最高，線性核評估準確性均最低，與Zhao等^［30］對豬育種數據的研究結果一致，這可能與所用的SNP數目相近有關。Kasnavi等^［31］通過模擬數據發現高斯函數核有較好的應用效果。由于挑選合適的C值與γ值需要更多的計算量^{［30，32］}，也有較多的研究直接選用推薦的參數值^［33-34］。在本研究中，SVM最佳參數與軟件推薦的默認值差別較小，在允許基因組預測準確性損失較小的情況下，可直接使用默認值進行計算。

在篩選RandomForest最佳參數時，在總產仔數、產活仔數、窩重等繁殖性狀中對ntree、mtry、nodesize等參數進行評估，發現基因組預測準確性隨著參數的變化展現一定的隨機性。Sarkar等^［35］在對nodesize與mtry最佳參數進行篩選時，發現其具有一定的規律性，這可能是由于使用參數較少且參數值差別較大所致。同時，RandomForest模型本身的隨機性（如節點特征值選擇與樣本分割）也可能導致基因組預測準確性隨著參數的變化展現一定的隨機性。

RandomForest模型在總產仔數、產活仔數、窩重中均表現出最高的評估準確性，其次為SVM，GBLUP、BayesB、BayesLasso等遺傳評估準確性較差且保持一致。與GBLUP、BayesB、BayesLasso等參數模型相比，RandomForest在繁殖性狀中的遺傳評估準確性可提升22.35%～34.25%，SVM在繁殖性狀中的遺傳評估準確性可提升6.15%～26.04%。在運行時間上，GBLUP運行時間保持最優，其次為RandomForest模型，BayesB、BayesLASSO、SVM等模型評估時間均較久。Banerjee等^［36］在使用SNP數據進行水稻遺傳評估時發現，RandomForest、SVM比傳統GS模型有更高的準確性。Liang等^［37］在中國西門塔爾肉牛研究中發現SVM、RandomForest較傳統GBLUP模型準確性更高，且RandomForest均有更好的穩定性。Merrick和Carter^［38］認為，隨著多年來對訓練群體的結合，在育種計劃中使用非參數機器學習模型來準確預測和實施復雜性狀的準確性更高。同時，在評估不同模型交叉驗證準確性時發現，其結果存在較強的一致性。然而，SVM與RandomForest在運行上需要更多的時間，且隨著數據量的增加，所需的時間呈指數式增加^［39-40］，這在一定程度上限制了機器學習方法的使用。

4 結 論

SVM與RandomForest非參數機器學習模型在豬繁殖性狀基因組選擇中具有一定的優勢，其評估結果與GBLUP、BayesB、BayesLasso等參數模型具有較高的相關性，但運行時間相對更長在一定程度上限制了該機器學習算法的使用。不過，隨著算法的研究優化和運行時長的逐步縮短，SVM與RandomForest仍具有較好的應用前景。

參考文獻（References）：

［1］ PILES M，BERGSMA R，GIANOLA D，et al.Feature selection stability and accuracy of prediction models for genomic prediction of residual feed intake in pigs using machine learning［J］.Front Genet，202 12：611506.

［2］ MUOZ M，BOZZI R，GARCA-CASCO J，et al.Genomic diversity，linkage disequilibrium and selection signatures in European local pig breeds assessed with a high density SNP chip［J］.Sci Rep，2019，9（1）：13546.

［3］ WELLMANN R，PREUSS S，THOLEN E，et al.Genomic selection using low density marker panels with application to a sire line in pigs［J］.Genet Sel Evol，2013，45（1）：28.

［4］ MEUWISSEN T H，HAYES B J，GODDARD M E.Prediction of total genetic value using genome-wide dense marker maps［J］.Genetics，200 157（4）：1819-1829.

［5］ VANRADEN P M.Efficient methods to compute genomic predictions［J］.J Dairy Sci，2008，91（11）：4414-4423.

［6］ LEGARRA A，AGUILAR I，MISZTAL I.A relationship matrix including full pedigree and genomic information［J］.J Dairy Sci，2009，92（9）：4656-4663.

［7］ MONTESINOS-L PEZ O A，MONTESINOS-L PEZ A，PREZ-RODRGUEZ P，et al.A review of deep learning applications for genomic selection［J］.BMC Genomics，202 56（1）：222-231