融合多組學數據的乳腺癌生存期預測模型研究

方秋蓮，周儀璇，伍 幸

(中南大學 數學與統計學院，湖南 長沙 410083)

0 引言

乳腺癌是一種發生于女性身上的惡性腫瘤[1]，是目前女性中發病率最高的惡性腫瘤[2].研究表明，大多數已確診乳腺癌的患者可以通過及時有效的治療被治愈或者延長生存期[3].當前乳腺癌診療研究的發展趨勢是根據患者的遺傳和臨床信息對其進行精準預后，并基于預測結果為其制定個性化的醫療方案.生存期預測是乳腺癌預后的重要研究內容之一[4].乳腺癌的發生和發展是一個復雜、多階段和連續的過程，不僅受遺傳因素(基因突變、表觀遺傳變化、細胞生物環境、個體特有特征)影響，也受患者生活環境因素影響，涉及多種生物分子水平的變化[5].近年來，得益于大規模組學數據的收集和計算方法的發展，從多個生物學層面解讀乳腺癌疾病機制的研究取得了重大進展.由高通量實驗技術產生的海量數據統稱為組學數據，組學數據為細胞的分子機制提供了詳細描述，不同層次組學數據在全基因組的分布狀態與乳腺癌的發生、發展、預后之間均存在關聯關系[6].目前的研究趨勢是綜合運用多個組學的信息對乳腺癌患者進行預后.與單一組學數據相比，多組學數據可以提供更全面的基因信息，利于更好地解釋分類結果[7]、改善預測效果[8]、理解復雜的分子機制[9].然而，多組學數據的融合面臨許多阻礙.一方面，組學數據自身具有缺失、類不平衡、噪聲大、復雜及維數災難等特點，這使計算變得困難[10].另一方面，不同來源的數據集具有異質性，具體表現為數據分布、類型、規模的差異.組學數據的這些特點為數據的融合過程帶來了挑戰.

現有的數據融合方法可按融合的時間點分為5類：早期、混合、中期、晚期和分層融合[11].早期融合即在建模前將所有數據集拼接為一個大型矩陣，這種方法思想簡單且易解釋，但會得到一個更加復雜高維的矩陣，增加了后續建模難度[11].混合融合是先將原始數據轉換為低維、少噪聲的簡單表達，然后融合數據的簡單表達[11-12].這類方法在進行融合前不僅降低了數據的維度，還減弱了數據集間的異質性，如數據規模.中期融合即直接融合未經轉換或拼接的數據集并輸出可用于后續分析的新的數據表達，這類融合方法假定所有數據集可投影到一個公共的隱空間[11，13].晚期融合也可稱為模型融合，即先分別對所有數據集建模，然后組合模型預測結果[11，14].分層融合即在神經網絡的連接結構中引入分子通路關系的先驗信息，可以在模型中體現分子層面的模塊結構[11，15].目前盡管已經有一些利用多組學數據預測乳腺癌患者生存期的研究，但是使用基于模型的晚期融合方法來提升預測精度的研究較少.本文采用基于模型的晚期融合方法，其優勢在于可以根據每一個數據集的特點分別選擇最合適的模型，充分考慮了數據集的異質性，且無須直接整合異質的數據.在融合環節，本文采用加權平均法組合模型預測結果，其中權重的取值由模型效果和預測結果間的互信息決定，這樣的組合方式可以改善預測效果，增強預測結果的可解釋性.

1 數據介紹及預處理

1.1 組學數據介紹

本文使用的數據來自 METABRIC數據庫[16]，包括1 981位乳腺癌患者的臨床信息數據、1 904位乳腺癌患者的基因表達信息數據以及2 173位乳腺癌患者的拷貝數變異信息數據，這3個數據集可通過共同字段——“PATIENT_ID”和“Hugo_Symbol”(均代表患者編號)進行連接，患者的生存時間變量包含在臨床信息數據集中.定義患者的生存時間是否大于5年為因變量，剔除隨訪期小于5年且生存狀態為存活的樣本，最終得到的基因表達數據集樣本量為1 844，臨床信息數據集和拷貝數變異數據集樣本量均為1 917，所有患者的最短生存期為0.1個月，最長生存期為29.60年，中位生存期為9.96年.

組學數據具有缺失、高維等特征，因此，需要先進行數據清洗、缺失值插補、特征選擇等數據預處理工作.

1.2 缺失值處理

本文首先剔除數據集中含缺失的樣本得到無缺失的數據集，然后在無缺失的數據集上人工隨機生成缺失，缺失的生成比例即原始數據集中該特征的缺失比例，然后基于此數據集從多個缺失值插補方法中選出效果最好的方法，并將該方法用于原始數據的缺失值插補.本文用到的缺失值插補方法有：隨機森林預測插補法、K最近鄰(KNN)預測插補法、簡單均值眾數插補法以及基于K-means聚類的均值眾數插補法.對于數值型特征，本文采用實際值與插補值的均方誤差來衡量插補效果，對于類別型特征，采用實際值與插補值的F1值來衡量插補效果.

1.3 特征選擇

本文的特征選擇過程均是在單個數據集上進行，分為3步：

1)采用假設檢驗法分別對3個數據集的特征進行初篩，選擇在長生存期和短生存期兩類樣本上的分布具有顯著差異的特征.對于數值型特征：首先分別對兩類樣本上的特征進行Shapiro正態性檢驗和Levene方差齊性檢驗.若兩類樣本上的特征取值均服從正態分布且具有等方差性，則對特征進行方差分析.若兩類樣本上的特征取值不全服從正態分布或不具有等方差性，則對兩類樣本上的特征進行Wilcoxon秩和檢驗，選擇方差分析或Wilcoxon秩和檢驗顯著的特征，將其納入特征子集[17].對于類別型特征：首先構建特征與因變量之間的列聯表，根據列聯表單元格中的計數來進行相關檢驗：若全部單元格計數T≥5且樣本總數n≥40，則用Pearson卡方檢驗；若存在單元格計數1≤T<5且樣本總數n≥40，則用連續性校正的卡方檢驗；若存在單元格計數T<1或n<40，則用Fisher精確檢驗[17].

2)對第一步選出的特征建立隨機森林模型，計算每個特征的基尼重要性和基于排列的重要性，剔除重要性為0的特征，綜合兩種特征重要性對剩下的特征進行排序.

3)將上述特征依序逐個加入Logistic回歸模型，若加入的特征能使模型的AUC提升10-6，則將加入的特征保留下來，否則將其剔除.

1.4 數據標準化

完成缺失值插補及特征選擇后，在訓練模型之前，需要對數據進行標準化.這是因為部分機器學習模型如支持向量機、KNN等是基于距離度量進行樣本分類和預測的，而距離對量綱和取值范圍非常敏感.數據標準化的過程如式(1)所示：

(1)

式中：X為原始特征；Xmean為X的平均值；Xstd為X的標準差；z為標準化后的特征.

2 正類預測概率加權融合模型構建

對數據進行預處理后，為了充分利用每一個數據集的信息，本文先為單一數據集選取最優分類模型，模型的融合是通過對3個數據集上的最優模型的正類預測概率進行加權求和實現.圖1展示的是本文構建多組學數據融合模型的流程.

圖1 融合模型構建流程Fig.1 Process of integrating model building

2.1 最優分類模型選擇

本文采用Logistic[18]、支持向量機[19]、隨機森林[20]、Xgboost[21]，以及BP神經網絡[22]5種機器學習分類算法分別對基因表達、拷貝數變異和臨床信息數據集建立5年生存期預測模型，并以模型在訓練集上的AUC值為評價指標，分別為每種數據集選擇一個最優預測模型.

2.2 數據融合策略

在確定了3種數據集上的最優預測模型后，本文通過加權最優模型的預測結果實現數據融合，并在融合時給效果更好且包含更多信息的模型賦予更高的權重.

假設ge_prob，cli_prob，cnv_prob分別表示基因表達、臨床信息和拷貝數變異數據集上最優模型的正類預測概率，即模型預測的樣本生存期大于5年的概率，定義融合模型的預測結果final_weight_prob為：

final_weight_prob=a·ge_prob+b·cli_prob+c·cnv_prob.

(2)

若final_weight_prob大于0.5，則表示預測生存期大于5年，若final_weight_prob小于0.5，則表示預測生存期小于5年.a，b，c為權重參數，當a=b=c=1/3時，融合結果為等權重的正類預測概率加權平均.然而，由于3個模型的效果及其所利用信息存在差異，其對最終正確預測的貢獻是不同的，故其對應的權重也應不同.

設G、C、V是基因表達、臨床信息、拷貝數變異數據訓練集上的預測結果，其取值g、c、v為0或1.令I(G，V)、I(G，C)、I(C，V)表示G、C、V間的互信息，其計算公式分別如式(3)～式(5)所示，由互信息的定義——兩個隨機變量共享的信息量，可知兩個預測結果的互信息值體現了其對應模型共享信息的多少，互信息值越高，說明這兩個模型共享的信息量越大，即其中任一模型所含的信息越少.

(3)

(4)

(5)

令AUCge、AUCcli和AUCcnv表示模型在訓練集上的AUC值，定義A、B、C的計算公式如式(6)～式(8)所示：

(6)

(7)

(8)

權重參數a，b，c的計算公式分別如式(9)～式(11)所示，可見效果越好且跟其他兩種模型共享信息越少的模型，其預測概率將具有越高的權重.

(9)

(10)

(11)

3 基于正類預測概率加權數據融合的乳腺癌生存期預測

3.1 數據預處理

首先，采用1.2節所述的5種缺失值插補方法對3個數據集中的缺失值進行插補，選定的插補方法如表1所示.

表1 選定的缺失值插補方法

然后分3步進行特征選擇，剔除了方差小于0.1的特征后，采用假設檢驗、特征重要性排序和剔除冗余特征，得到3個數據集的最優特征數目如表2所示.

表2 3個數據集的最優特征數目

3.2 模型評價指標與調參策略

本文采用基于混淆矩陣的敏感度、特異性、準確率、ROC曲線和AUC值來評價模型分類效果.敏感度、特異性、準確率的計算公式分別如式(12)～式(14)所示.

(12)

(13)

(14)

式中，nTP、nTN、nFP、nFN分別為真正例、真負例、假正例、假負例個數.

本文采用分層抽樣的方式隨機抽取80%的樣本作為訓練集，20%的樣本作為測試集.在訓練集上采用分層五折交叉驗證、隨機網格搜索以及網格搜索的方式調整模型的參數.

3.3 最優模型選擇

本文先分別對所有數據按照式(1)進行標準化，然后在基因表達、臨床信息和拷貝數變異數據集上訓練分類模型.圖2中的3幅圖分別為5種分類模型在3種數據對應訓練集上的ROC曲線.在基因表達數據集上，BP神經網絡的AUC值最高，為0.826 8；在臨床信息數據集上，BP神經網絡的AUC值最高，為0.757 7；在拷貝數變異數據集上，支持向量機的AUC值最高，為0.686 0.以AUC值為主要評價指標，確定基因表達和臨床信息數據集的最優分類模型為BP神經網絡模型，拷貝數變異數據集的最優分類模型為支持向量機模型.

圖2 5種分類算法在3個訓練集上的ROC曲線及AUC值：(a)基因表達數據；(b)臨床信息數據；(c)拷貝數變異數據 Fig.2 ROC and AUC of 5 types of classifiers for the three train sets：(a)Gene expression dataset；(b)Clinical dataset；(c)Copy number variation dataset

3.4 融合模型應用及其性能比較

為驗證正類預測概率加權融合模型的有效性，本文綜合對比基于單一數據集的模型和融合模型在測試集上各評價指標的結果.表3展示了各模型在測試集上的敏感度、特異性、準確率和AUC值.其中BP_gene和BP_cli分別表示基于基因表達和臨床信息數據集預測的BP神經網絡模型，SVM_cnv表示基于拷貝數變異數據集預測的支持向量機模型，equal_fusion表示等權重的正類預測概率加權融合模型，weight_fusion表示正類預測概率加權融合模型，其中權重參數a=0.275 0、b=0.488 1、c=0.236 9.敏感度最高的是基于拷貝數變異數據集預測的支持向量機模型，其敏感度為0.536 6；特異性最高的是正類預測概率加權融合模型，其特異性為0.968 6；準確率最高的是等權重的正類預測概率加權融合模型，其準確率均為0.813 0；AUC值最高的是正類預測概率加權融合模型，其AUC值為0.815 6.與基于單一組學數據集的模型相比，基于多個數據集的模型具有更高的特異性、準確率和AUC，說明融合多組學數據可提升預測效果；與等權重的正類預測概率加權融合模型相比，加權融合模型的特異性和AUC更高，說明本文提出的權重是有效的.

表3 各模型在測試集上的評價指標結果

4 結論

通過高通量技術獲得的大量組學數據為個性化診療奠定了基礎.本文選取基因表達、臨床信息和拷貝數變異數據進行分析，分別對3個數據集進行數據清洗、缺失值插補、特征選擇、最優分類模型選取，并通過求3個最優模型的正類預測概率的加權平均得到融合模型預測結果.正類預測概率加權融合模型不僅同時利用了3個數據集的信息，還在確定權重的過程中綜合考慮了各個模型的準確度和它們之間的互信息.本文通過實例證明了正類預測概率加權融合模型具有良好的分類效果，其預測結果可為醫生制定高效、低副作用的治療方案提供參考.正類預測概率加權融合模型仍有提升空間，例如考慮非隨機缺失，應用更先進的降維、分類算法，融合更多種類的組學數據及提升預測效果等.