引入遺傳因子的骨密度機器學習回歸模型研究

陳鵬麗 孔祥勇 林勇

上海理工大學醫療器械與食品學院，上海 200093

骨質疏松癥是一種常見的復雜疾病，受到多種因素的影響，包括臨床危險因素(個人的年齡、性別、體重、身高、飲食、既往骨折史、長期使用糖皮質激素等)以及遺傳因素[1]。有多項研究[2-4]證明，遺傳因素對骨密度(bone mineral density，BMD)具有顯著影響，骨密度差異的遺傳率在50 %～85 %之間。全基因組關聯研究(genome-wide association study, GWAS)及薈萃分析(Meta-analysis)發現，相較只使用臨床危險因素數據，加入骨密度相關的遺傳變異因素可以顯著提高骨折預測的準確度[5]。

目前，全基因組關聯研究和薈萃分析已經發現了許多與骨密度、骨質疏松癥和骨質疏松性骨折相關的位點[6-7]。然而復雜疾病的多因素性質使得傳統的基于統計學分析方法的GWAS效果有限[8]。近年來研究者們[8]提出多種用于檢測骨質疏松癥致病因素的機器學習模型和方法，在對多種復雜疾病建模以及易感位點識別中取得了突出的成果。

但相關研究大多數只以單一因素對骨質疏松癥進行分析，很少考慮遺傳因素和臨床因素特征之間的相互作用，且缺乏識別在生物學上具有可解釋性的易感位點的能力。

本文提出了一種基于機器學習的骨質疏松癥致病因素分析方法，引入遺傳因素及臨床風險因素，識別影響個體對骨質疏松易感性的最優特征組合。研究采用兩階段的特征選擇方法，首先通過最大互信息系數(maximal information coefficient, MIC)篩選出與預測變量高度相關的SNP子集，然后將所選子集混合臨床風險因素作為輸入進行序列特征選擇(sequential feature selection, SFS)。下面將對數據預處理、兩階段特征選擇和預測的過程進行詳細描述，并以模型均方根誤差(root mean square error, RMSE)為指標，在不同的回歸模型中進行對比分析，衡量這一方法的性能。

1 材料和方法

1.1 實驗數據

本研究中選取的實驗樣本為2 263例白種人，包括555名男性和1 708名女性。該研究得到相關機構審查委員會批準，所有研究參與者在進入項目前都簽署了知情同意書[9]。

研究使用全身骨密度值作為因變量，將對骨密度的影響因素劃分為臨床風險因素和遺傳因素。以數據較全的年齡、性別(男性：1，女性：2)、體重、身高、體質量指數(BMI)作為臨床風險因素；遺傳因素為單核苷酸多態性(SNP)基因分型數據。受試者臨床風險因素基本信息見表1。

表1 受試者臨床風險因素基本信息Table 1 Basic information of clinical risk factors of subjects

1.2 引入遺傳因素的機器學習骨密度回歸模型

本文提出的對骨密度回歸模型的分析流程如圖1所示，包括4個階段：①數據預處理；②特征選擇；③建立預測模型；④模型評估。

圖1 骨密度回歸模型分析流程Fig.1 Flow chart of BMD regression model

1.2.1數據預處理：①數據轉換：使用Plink全基因組關聯分析工具通過加性顯性編碼將SNP基因分型數據格式(AA, AC, CC, NC)轉換成數值形式(0, 1, 2, NC)并存放在文本文件中。②數據補缺：SNP基因分型數據集中缺失的值通常標記為“NC”。數據補缺遵循以下標準：分型后“NC”值多于10 %的SNP位點將被丟棄，其余位點的“NC”值將被該位點的眾數(SNP數據集中每個位點出現頻率最多的編碼)替換。

預處理完成后，SNP數據集包含2 263例樣本，35 780個特征。由此，特征組合分為臨床風險因素特征組、臨床風險因素結合遺傳因素特征組。臨床風險因素特征組由年齡、性別(男性：1，女性：2)、體重、身高、BMI及其對應項平方組成，共10維特征。在此基礎上加入SNP位點特征即為臨床風險因素結合遺傳因素特征組。

1.2.2特征選擇：本文提出一種兩階段特征選擇算法，首先以MIC作為過濾式方法剔除SNP數據集中大量的噪聲數據，再結合臨床特征變量使用序列浮動前向選擇算法(sequential floating forward selection, SFFS)這一封裝式算法，獲得信息最豐富的特征子集。

MIC由 Reshef等[10]提出，它基于互信息度量變量對之間的相關性，在數據量巨大的情況下，互信息能夠有效地表述變量間的非線性相關關系[11]。

特征Xi與預測變量Y的MIC定義如下：給定雙變量(Xi,Y)組成的數據集D，首先進行網格分區形成維度為(x,y)的網格G。對于給定的D，改變網格G的劃分方式，D∣G表示D在G上的概率分布，與落在每個子格內的散點的數量成正比，計算不同劃分方式下的最大互信息。定義特征矩陣M(D)x,y=(mx,y)，其中mx,y是任意x×y網格所獲得的最高歸一化互信息值，特征矩陣的第(x,y)項mx,y為：

M(D)x,y=[maxI(D∣G)]/log2min(x,y)

(1)

定義統計值MIC：

MIC(Xi,Y)=maxM(D)x,y，xy

(2)

其中B是關于樣本大小的函數，通常設B=n0.6，n為訓練集中的樣本數。

本文利用MINE(maximal information-based nonparametric exploration)算法[12]計算位點特征與預測變量BMD之間的MIC ，選擇MIC得分最高的前m個SNPs作為下一階段特征選擇方法的輸入。

SFS算法從一個空的特征子集開始，通過不斷添加(或移除)特征直到選擇出最優特征子集或達到預先指定的子集大小。SFFS是增加了回退機制的SFS算法。

具體的特征選擇流程見圖2。為了選擇出信息最豐富的特征子集，首先利用MINE算法計算SNP特征與BMD之間的MIC[10]，選擇MIC得分最高的前m個SNPs結合臨床數據作為下一階段SFFS的輸入。特征在SFFS算法中經過k次迭代，在每次迭代過程中，以準則函數最大化為目的，從特征空間中選擇一個最佳特征，并通過額外的排除或包含步驟，檢查若移除一個特征后，特征子集能否提高預測的性能。最后根據不同數量的特征組合，選擇達到最優預測效果的特征子集N。提出的特征選擇算法在python 3.7上實現。

1.2.3預測模型：本研究采用隨機森林[13]作為主要的預測模型，構建機器學習預測模型的一個關鍵步驟是優化超參數以獲得最佳模型性能，選擇優化以下兩個隨機森林模型的超參數：①n_estimators，森林中決策樹的數量。②max_features，建立決策樹時選擇的最大特征數。首先使用網格搜索確定超參數n_estimators的最優值，然后選擇max_features。

通過十折交叉驗證確定最終的最優特征子集。除隨機森林之外，另有幾種經典的回歸算法，即支持向量機回歸(SVR)、線性回歸(LR)、XGBoost用于測試我們提出的兩階段特征選擇方法以及最優特征子集的有效性和穩定性。

1.2.4評估方法：為衡量回歸模型的預測精度和泛化能力，使用均方根誤差為指標對模型進行評估。均方根誤差是評估回歸模型與數據集擬合程度的一種方法，其計算公式為：

(3)

其中，Yi是數據集中第i個樣本的預測值，f(xi)是數據集中第i個樣本的實際值，N是樣本容量。

2 結果

本研究使用2 263例白種人樣本的臨床風險因素和遺傳因素數據集，對上文提出的基于機器學習的骨質疏松癥致病因素分析方法進行驗證。

在對數據集進行預處理后，SNP數據集包含35 780個位點特征。首先以MIC作為過濾式方法剔除SNP數據集中大量的噪聲數據，最終保留MIC得分最高的前100個SNP作為遺傳特征，與臨床特征混合進行下一步的SFFS，選擇過程基于隨機森林回歸模型，在樣本數據集中進行十折交叉驗證，評估標準為RMSE。

如圖3所示，當特征數為57時，模型的RMSE達到最低為0.093 598 g/cm3，即這組特征(包含51個SNP位點，6個臨床特征)能夠最好地擬合實際的骨密度值。

在訓練模型之前對隨機森林回歸模型的超參數進行選擇，表2為訓練和優化隨機森林算法所測試的超參數及其對應的測試值和最終確定的取值。

表2 隨機森林模型超參數取值Table 2 Hyperparameters in random forest model

網格搜索發現分類器個數的最優值為170、建立決策樹時選擇的最大特征數的最優值為25，此時模型的均方根誤差最低。

2.1 加入遺傳因素對模型均方根誤差的影響

在模型輸入為臨床因子特征變量、遺傳因子加臨床因子特征變量兩種情況下，以隨機森林為回歸模型，進行3組實驗。測試集的 RMSE值如圖4所示，每組實驗結果的均值以數字標注在條形圖之上，標準差以誤差線形式表示。

圖4 加入遺傳因素前后模型的RMSE值Fig.4 The effects of introducing genetic factors on RMSE value of the model

分析加入遺傳因素前后模型的均方根誤差，加入遺傳因素前，3組十折交叉驗證實驗的RMSE分別為(9.88±0.51)×10-2g/cm3、(9.90±0.48)×10-2g/cm3、(9.88±0.52)×10-2g/cm3，30次實驗的平均RMSE為(9.89±0.51)×10-2g/cm3；加入遺傳因素后，3組十折交叉驗證實驗數據的RMSE分別為(9.35±0.50)×10-2g/cm3、(9.37±0.50)×10-2g/cm3、(9.36±0.49)×10-2g/cm3，30次實驗的平均RMSE為(9.36±0.50)×10-2g/cm3。

相比只以臨床危險因素為特征，在加入遺傳因素后，3組交叉驗證中模型的RMSE分別降低了5.36 %、5.35 %、5.26 %，30次實驗的平均RMSE降低了5.36 %。說明加入遺傳因素作為特征變量能夠降低模型的均方根誤差，使骨密度回歸模型更好地擬合數據集。

2.2 不同特征選擇策略對模型均方根誤差的影響

為證明如上提出的融合最大互信息系數和序列浮動前向選擇這一兩階段特征選擇方法的有效性，使用相同的數據集，混合遺傳因素與臨床因素特征，將本研究所提兩階段特征選擇方法與僅使用MINE算法、以及第一階段使用MINE算法，第二階段使用另一種經典的封裝式算法遞歸特征消除(recursive feature elimination, RFE)進行對比。

如圖5所示，僅使用最大互信息系數進行特征選擇，3組十折交叉驗證實驗數據的RMSE分別為(9.47±0.53)×10-2g/cm3、(9.48±0.52)×10-2g/cm3、(9.46±0.55)×10-2g/cm3，30次實驗的平均RMSE為(9.47±0.54)×10-2g/cm3；第一階段使用最大互信息系數，第二階段換用遞歸特征消除(RFE)選擇最優特征子集，3組十折交叉驗證實驗數據的RMSE分別為(9.49±0.55)×10-2g/cm3、(9.49±0.53)×10-2g/cm3、(9.48±0.52)×10-2g/cm3，30次實驗的平均RMSE為(9.48±0.54)×10-2g/cm3；本研究所提混合最大互信息系數與序列浮動前向選擇的特征選擇算法30次實驗的平均RMSE為(9.36±0.50)×10-2g/cm3，為3種方法中最低。

圖5 所提特征選擇方法與僅使用MINE算法、第一階段MINE算法+第二階段RFE選擇特征的效果對比Fig.5 Comparison of RMSE for MINE, MINE+RFE and the proposed MINE+SFS method

2.3 使用不同回歸模型的對比分析

為證明所提兩階段特征選擇方法以及遴選出來的特征子集的穩定性，將最優特征子集輸入不同的回歸模型(RF, SVR, LR, XGBoost)進行十折交叉驗證。其中，SVR模型采用的核函數為徑向基核函數(RBF)，已經對輸入特征和標簽值進行數據標準化處理。LR及XGBoost采用標準參數。對比結果如圖6所示。

圖6 特征子集在使用不同回歸模型時的RMSEFig.6 Comparison of RMSE for different regressors

加入遺傳因素作為特征變量前后，RF的RMSE分別為(9.90±0.48)×10-2g/cm3和(9.35±0.50)×10-2g/cm3；SVR的RMSE分別為(10.90±0.46)×10-2g/cm3和(9.47±0.49)×10-2g/cm3；LR的RMSE分別為(9.95±0.29)×10-2g/cm3和(9.53±0.52)×10-2g/cm3；XGBoost的RMSE分別為(10.41±0.51)×10-2g/cm3和(10.32±0.49)×10-2g/cm3。在使用本文提出的兩階段特征選擇算法所遴選出來的特征子集后，幾種回歸模型的RMSE都有明顯降低，其中隨機森林回歸模型在臨床因素混合遺傳因素特征集上RMSE最低。

3 討論

相較于既往的骨質疏松單一致病因素分析研究，本文引入遺傳因素，旨在基于SNP數據集結合臨床危險因素實現對這一復雜疾病更準確、更魯棒的預測。為提高算法性能并減少時間復雜度，同時保留特征的生物學含義和解釋性，采用最大互信息系數作為過濾式方法剔除SNP數據集中大量的噪聲數據，最終保留MIC得分最高的前100個SNP作為遺傳特征，與臨床特征混合，基于隨機森林回歸模型進行下一步的序列浮動前向選擇。初步構建出具有良好的預測準確度和穩定性的骨質疏松癥致病因素分析方法。

研究提出的兩階段特征選擇方法兼顧封裝式方法的精度和過濾式方法的效率，可以實現時間復雜度較低的非線性預測模型，降低預測誤差，為骨質疏松癥以及類似的復雜疾病的致病因素探明、預測模型的建立提供有價值的參考。

低骨密度可能是多種致病途徑的共同結果，這些途徑受到遺傳因素的影響，本研究基于2 263例白種人樣本，建立骨密度機器學習回歸模型，篩選出如表2的51個骨質疏松癥易感位點。

表3 方法篩選出的51個易感位點Table 3 Characteristics of 51 SNPs selected by the proposed method

這些位點位于13個基因，其多態性通過多種途徑影響骨密度。基因GNG12-AS1、WLS、MEF2C、CDKAL1和SFRP4通過激活Wnt/β-catenin信號通路調控骨形成，Wnt/β-catenin信號通路及相關蛋白在骨細胞分化、增殖和凋亡的過程中至關重要[14-17]。基因SUPT3H與PKDCC調控間充質干細胞與軟骨細胞的分化過程[18-19]。AKAP11、RPS6KA5與SMG6參與骨細胞生長發育過程，在此前的研究中被證實與骨質疏松特征具有顯著相關性[20-24]。除了基質蛋白和骨細胞的平衡外，骨骼的完整性還依賴于礦物質的穩態。GALNT3調控循環中的磷酸鹽水平[17]，KCNMA編碼細胞上Big K+(BK)大電導鈣和電壓激活的K+通道的成孔性α亞基[25]，參與骨吸收。基因ESR1編碼雌激素受體α，在調節骨量和骨質疏松癥的發生中發揮重要作用[26]。

值得注意的是，在我們識別出的易感位點中，已經有一部分在此前的GWAS研究被鑒定出與骨質疏松、骨折等性狀顯著相關，即rs726282[27]、rs932477[28]、rs2179922[29]、rs2941740[30]、rs4952590[31]、rs6721582[32]，這進一步證實了本文所提方法的有效性。本文創新性的將臨床風險因素和遺傳因素結合，通過機器學習方法識別骨質疏松癥易感位點。所提分析方法識別出的其他致病SNP揭示了遺傳因素之間、遺傳因素與臨床因素之間相互作用的存在，既有在臨床上對骨質疏松癥的預測意義，也為我們未來對骨質疏松癥致病因素的進一步探究提供了潛在的靶點。

研究存在一定的局限性：其一，數據樣本規模較小。后續研究將結合更大樣本的數據集對分析方法的性能進行進一步的驗證，并在全基因組范圍內探索骨質疏松癥致病因素，完善分析方法，提高模型的計算效率和識別能力，以期全方位且高效地識別骨質疏松癥的關聯位點和基因；其二，研究沒有對分析方法的運算時間進行定量分析，后續研究可以嘗試選擇不同數量的特征子集、優化機器學習算法等方法以提升模型的時間性能。