陳鵬麗 孔祥勇 林勇

上海理工大學(xué)醫(yī)療器械與食品學(xué)院，上海 200093

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

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

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

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

1 材料和方法

1.1 實(shí)驗(yàn)數(shù)據(jù)

本研究中選取的實(shí)驗(yàn)樣本為2 263例白種人，包括555名男性和1 708名女性。該研究得到相關(guān)機(jī)構(gòu)審查委員會批準(zhǔn)，所有研究參與者在進(jìn)入項(xiàng)目前都簽署了知情同意書[9]。

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

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

1.2 引入遺傳因素的機(jī)器學(xué)習(xí)骨密度回歸模型

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

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

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

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

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

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

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

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

(1)

定義統(tǒng)計值MIC：

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

(2)

其中B是關(guān)于樣本大小的函數(shù)，通常設(shè)B=n0.6，n為訓(xùn)練集中的樣本數(shù)。

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

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

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

1.2.3預(yù)測模型：本研究采用隨機(jī)森林[13]作為主要的預(yù)測模型，構(gòu)建機(jī)器學(xué)習(xí)預(yù)測模型的一個關(guān)鍵步驟是優(yōu)化超參數(shù)以獲得最佳模型性能，選擇優(yōu)化以下兩個隨機(jī)森林模型的超參數(shù)：①n_estimators，森林中決策樹的數(shù)量。②max_features，建立決策樹時選擇的最大特征數(shù)。首先使用網(wǎng)格搜索確定超參數(shù)n_estimators的最優(yōu)值，然后選擇max_features。

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

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

(3)

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

2 結(jié)果

本研究使用2 263例白種人樣本的臨床風(fēng)險因素和遺傳因素數(shù)據(jù)集，對上文提出的基于機(jī)器學(xué)習(xí)的骨質(zhì)疏松癥致病因素分析方法進(jìn)行驗(yàn)證。

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

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

在訓(xùn)練模型之前對隨機(jī)森林回歸模型的超參數(shù)進(jìn)行選擇，表2為訓(xùn)練和優(yōu)化隨機(jī)森林算法所測試的超參數(shù)及其對應(yīng)的測試值和最終確定的取值。

表2 隨機(jī)森林模型超參數(shù)取值Table 2 Hyperparameters in random forest model

網(wǎng)格搜索發(fā)現(xiàn)分類器個數(shù)的最優(yōu)值為170、建立決策樹時選擇的最大特征數(shù)的最優(yōu)值為25，此時模型的均方根誤差最低。

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

在模型輸入為臨床因子特征變量、遺傳因子加臨床因子特征變量兩種情況下，以隨機(jī)森林為回歸模型，進(jìn)行3組實(shí)驗(yàn)。測試集的 RMSE值如圖4所示，每組實(shí)驗(yàn)結(jié)果的均值以數(shù)字標(biāo)注在條形圖之上，標(biāo)準(zhǔn)差以誤差線形式表示。

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

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

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

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

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

如圖5所示，僅使用最大互信息系數(shù)進(jìn)行特征選擇，3組十折交叉驗(yàn)證實(shí)驗(yàn)數(shù)據(jù)的RMSE分別為(9.47±0.53)×10-2g/cm3、(9.48±0.52)×10-2g/cm3、(9.46±0.55)×10-2g/cm3，30次實(shí)驗(yàn)的平均RMSE為(9.47±0.54)×10-2g/cm3；第一階段使用最大互信息系數(shù)，第二階段換用遞歸特征消除(RFE)選擇最優(yōu)特征子集，3組十折交叉驗(yàn)證實(shí)驗(yàn)數(shù)據(jù)的RMSE分別為(9.49±0.55)×10-2g/cm3、(9.49±0.53)×10-2g/cm3、(9.48±0.52)×10-2g/cm3，30次實(shí)驗(yàn)的平均RMSE為(9.48±0.54)×10-2g/cm3；本研究所提混合最大互信息系數(shù)與序列浮動前向選擇的特征選擇算法30次實(shí)驗(yàn)的平均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 使用不同回歸模型的對比分析

為證明所提兩階段特征選擇方法以及遴選出來的特征子集的穩(wěn)定性，將最優(yōu)特征子集輸入不同的回歸模型(RF, SVR, LR, XGBoost)進(jìn)行十折交叉驗(yàn)證。其中，SVR模型采用的核函數(shù)為徑向基核函數(shù)(RBF)，已經(jīng)對輸入特征和標(biāo)簽值進(jìn)行數(shù)據(jù)標(biāo)準(zhǔn)化處理。LR及XGBoost采用標(biāo)準(zhǔn)參數(shù)。對比結(jié)果如圖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都有明顯降低，其中隨機(jī)森林回歸模型在臨床因素混合遺傳因素特征集上RMSE最低。

3 討論

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

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

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

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

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

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

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