基因表達對膀胱癌生存結局的預測：稀疏與混合Cox模型的實證比較研究*

徐州醫科大學公共衛生學院流行病與衛生統計學系(221004) 陸皓杰 曾 平 黃水平

【提 要】 目的 研究稀疏Cox(coxlasso)與混合Cox模型(coxlmm)在全基因表達數據中對膀胱癌預后的預測表現。方法 通過計算一致性指數(C-index)評價兩種模型在膀胱癌全基因表達數據中(TCGA，GSE31684和GSE13507)的預測精度，同時在混合Cox模型中將膀胱癌的生存方差劃分為臨床(PCE)和轉錄組(PGE)兩部分。結果 當TCGA數據集為訓練集時，coxlmm預測能力(C-index=0.676)高于coxlasso(C-index=0.655)，兩者在外部驗證集的C-index分別為0.527和0.534。當三個合并數據集為訓練集時coxlmm(C-index=0.671)比coxlasso(C-index=0.650)的預測精度提高2.1%。當GSE31684為訓練集時，coxlmm(C-index=0.553)比coxlasso(C-index=0.550)的預測精度提高0.3%，兩個模型在外部驗證集上的C-index分別為0.632和0.633。生存方差劃分表明膀胱癌的臨床貢獻高于轉錄組的貢獻(PCE=14.95%，PGE=10.88%)。結論 成功構建了一種用于膀胱癌的預后預測的coxlmm模型，揭示了整合全轉錄組信息可在一定程度上提高膀胱癌預后預測能力。

膀胱癌(bladder cancer)是最常見的癌癥之一，全世界每年新發病例約54.9萬，死亡病例約20萬，其中男性的發病率和死亡率均高于女性，分別為9.6/10萬和3.2/10萬[1]。根據其腫瘤分期，膀胱癌可以分為肌肉浸潤性膀胱癌(MIBC)和非肌肉浸潤性膀胱癌(NMIBC)，其中MIBC約占膀胱癌早期診斷病例的25%，NMIBC進展為MIBC的比例高達10%～15%[2-3]。隨著醫療水平的不斷提高，膀胱癌的預后治療也取得了長足進步；然而，由于膀胱癌疾病高度異質的緣故[4]，MIBC患者的預后仍然較差。不同的遺傳標記為深入研究疾病的遺傳基礎、發展新的診斷技術和治療方法提供了全新的視角[5]，因此整合大規模組學數據同時聯合臨床信息為膀胱癌預后評價帶來了新的工具[6]。

在先前的研究中組學信息已經被廣泛應用于膀胱癌的預后預測。然而，先前的研究只是將臨床信息和個別生物標志物納入預測模型；例如，有研究者僅將4個基因(TMPRSS11E，SCEL，KRT78和TMEM185A)納入膀胱癌預后預測模型[7]。另一些研究者通過降低數據維度的方法(如Lasso方法)提取重要生物標志物進行預后預測[8]。從統計方法角度看，這些方法可被認為是稀疏模型(sparse models)，因為其明確假設只有少部分組學信息對預測有用。雖然稀疏模型已經被證實能夠提高患者預后預測的準確性，也利于臨床實踐應用；然而，單個或少數生物標志物的預測性能往往不穩定，而生物標志物的組合能提高預測性能[9]。最近的研究表明，將100、300和5000個mRNA整合到模型中時，預測精度可從0.58、0.62提高到0.64[10]，表明在一定程度上納入更多生物標志物可以進一步提高預測效果。

在遺傳預測方面，混合模型(mixed models)同樣也展現出較高的準確性[11-13]；與稀疏模型不同，混合模型假設所有基因都參與了疾病進展并且對疾病均有影響。很顯然，混合模型和稀疏模型的預測效果優劣取決于上述假設與真實情況是否吻合。然而，實際中基因與疾病的真實關系往往未知；因此，對混合模型和稀疏模型的預測性能進行比較和評價具有重要的意義。有研究者利用全基因組表達數據對比評價了九種遺傳預測模型方法(包括稀疏模型和多基因模型)[14]。本文通過整合公開的膀胱癌數據，使用C-index來比較混合Cox模型(coxlmm)與稀疏Cox模型(coxlasso)對膀胱癌的預后預測表現；同時在混合Cox模型中將膀胱癌的生存方差劃分為臨床和轉錄組兩部分。

材料和方法

1.研究人群

從癌癥基因組圖譜(TCGA)[15]下載了430名膀胱癌患者的信息，保留同時具有基因表達數據和臨床信息的345名患者，同時將年齡、性別、腫瘤分期和淋巴結狀態作為臨床協變量納入研究。其中男性254人(73.6%)，女性91人(26.4%)，平均年齡為68.2歲(34～89歲)。從美國國家生物技術信息中心(National Center of Biotechnology Information，NCBI)的基因表達綜合數據庫(gene expression omnibus，GEO)中下載基因表達譜膀胱癌樣本數據集GSE31684和GSE13507[16-17]。GSE31684數據集中男性有58人(75.3%)，女性有19人(24.7%)，平均年齡為67.0歲(42～85歲)；在最后隨訪中有32名患者死于膀胱癌，有22名患者死于其他原因，其余23名患者仍然存活。GSE13507數據集中男性有134人(81.7%)，女性有30人(18.3%)，平均年齡為65.2歲(24～88歲)；在最后隨訪中有68人死于膀胱癌。這三個數據集的各項臨床特征和生存信息詳見表1。

表1 膀胱癌患者的臨床和生存信息

2.基因表達的數據集

我們將原始基因表達的count數轉換為log2(count+1)并通過R包limma中的normalizeQuantiles函數進行正規化處理[18-19]，剔除沒有名稱和表達量為0高于50%的基因之后還剩下29668個基因。使用R包GEOquery[20]下載GSE31684基因表達，該數據已經過GCRMA算法歸一化處理[21]；從R包hgu133plus2.db中的hgu133plus2SYMBOL獲得探針與基因之間的對應關系。GSE13507基因表達數據經過分位數標準化和log2轉換[17]，從R包illuminaHumanv2.db中的illuminaHumanv2SYMBOL獲得探針與基因的對應關系。剔除兩個數據集中沒有對應基因名稱的探針，當同一基因對應多個探針時，將表達量平均值最高的探針作為基因對應的唯一探針，然后把這組中的其他探針剔除，最后刪除表達量為0超過50%的基因之后兩個數據集分別保留了20174和15685個基因。將三個數據集的基因表達數據按基因名合并后，保留12517個共同的基因進行后續分析。此外，使用R包ComBat[22]中的sva消除三個數據集的批次效應[23]。

3.一般Cox模型

設Xi為個體i的臨床協變量(如TCGA數據集中的腫瘤分期、年齡、性別和淋巴結狀態[15])，并對每個X進行標準化，即X服從均值為0和方差為1的標準正態分布，ti表示個體i的生存時間，樣本量設為n。建立如下的一般Cox模型[24]：

(1)

其中h0(t)是所有協變量取值為0時的基線風險函數，a=(a1，a2，…，ap)是臨床協變量的系數。

4.具有臨床協變量和基因表達水平的稀疏Cox模型

設Gi為m維向量，表示個體i的一組遺傳標志物(如TCGA數據中基因表達水平)，并對每個標志物進行標準化，即G服從均值為0和方差為1的標準正態分布，包含Xi和Gi的Cox模型為：

(2)

其中b=(b1，b2，…，bm)為遺傳標志物效應大小的m維向量。在高維的背景下，基因數量遠遠大于樣本大小(即m?n)，所以傳統偏似然估計已不適用。在過去的幾年里，出現了許多懲罰來估計未知參數的正則化方法[25-27]。本文使用lasso降低數據維度[27]，lasso利用系數的絕對值函數來完成對模型系數的壓縮，將一些基因表達的系數收縮為零以同時執行特征選擇[28]。具體來說，lasso符合以下Cox模型(用coxlasso表示)：

(3)

5.具有臨床協變量和基因表達水平的線性混合Cox模型

如前所述，正則化方法(如coxlasso)本質上是稀疏模型，它假設只有少數基因與生存相關。然而線性混合Cox模型(coxlmm)將所有基因納入模型中[12，31-32]，則可認為coxlmm是多基因模型，因為其明確假設所有基因都參與了疾病的進展，并且對生存均具有影響[33-34]。在coxlmm模型中假設基因效應服從正態分布[35]

(4)

(5)

基于上述關系構造了一個新的等價的Cox混合模型：

(6)

其中Zi是K1/2和K=GGT的第i行向量，在遺傳預測中通常K被稱為遺傳相關矩陣[12-13]；δ是Zi效應大小的n維向量。在變換之后，隨機效應的維數從m減少到n。因為在TCGA和GEO數據集中n通常比m小得多，這使得coxlmm的計算效率大大提高。最后b的估計值為：

(7)

6.TCGA中膀胱癌臨床和轉錄組信息的相對重要性

為了揭示臨床和轉錄組信息在膀胱癌生存變異中的不同作用，需要量化其對生存表型的重要性[39-41]，我們首先對風險函數進行對數變換，然后定義兩個變量PCE和PGE[35]

(8)

前者代表有多少生存變異可以被臨床信息單獨解釋，后者揭示有多少生存變異可以用轉錄組的信息來解釋，其中var(x)是x的方差，通過E(x-E(x)2)/n計算，E(x)是x的期望，var(e)表示除X和G之外所有未解釋部分的方差。PCE和PGE的總和可以看作是模型中現有臨床和轉錄組信息共同解釋的生存方差(PVE)的比例。我們使用Jackknife方法來得到PCE或PGE的95%置信區間[42]。

7.模型評價

本研究使用一致性指數(C-index)來評估模型預測的準確性[43]。C-index為0.5表示完全不一致，說明模型預測能力不強；而C-index為1表示完全一致，表明模型的預測結果完全符合實際情況。參照先前的工作[13]，本研究使用蒙特卡羅交叉驗證(Monte Carlo cross validations，MCCV)評估預測性能：首先將訓練組患者(訓練集的80%)和試驗組患者(訓練集的20%)從訓練集中隨機抽取100次，然后測試組用于驗證之前使用訓練組估計的模型。在本文的訓練集預測分析中，與coxlasso模型相比，coxlmm模型通過整合全基因表達數據顯示出了更加穩健的預測性能。然而，coxlmm性能的提高是否與加入了有用的基因表達數據有關仍是未知的。因此，本研究在外部驗證數據集上計算C-index來進一步評估coxlmm的預測性能。

結 果

1.訓練集中Cox、coxlasso和coxlmm模型的預測能力

在TCGA訓練集中，Cox、coxlasso和coxlmm的C-index分別為0.652、0.655和0.676(圖1)。在三個合并數據訓練集中，Cox、coxlasso和coxlmm的C-index分別為0.651、0.650和0.671(圖2)。在GSE31684訓練集中，Cox、coxlasso和coxlmm的C-index分別為0.553、0.550和0.553(圖3)。在大多數情況下與另外兩個模型相比，coxlmm預測性能最好。例如，當TCGA數據集為訓練集時，coxlmm的預測精度比Cox平均提高2.4%；比coxlasso平均提高2.1%。當三個合并數據集為訓練集時，coxlmm的預測精度比Cox平均提高2.0%，比coxlasso的預測精度平均提高2.1%。然而當GSE31684數據集為訓練集時，coxlmm與Cox的預測精度相同，coxlmm的預測精度比coxlasso平均高出0.3%。

圖1 TCGA訓練集中的蒙特卡羅交叉驗證(MCCV)

圖2 三個合并數據訓練集中的蒙特卡羅交叉驗證(MCCV)

圖3 GSE31684訓練集中的蒙特卡羅交叉驗證(MCCV)

2.驗證集中Cox、coxlasso和coxlmm的預測能力

本文在使用外部數據集的情況下通過計算C-index來進一步評估coxlmm對預后的預測能力(表2)，結果顯示coxlmm比其他模型的預測精度更高。例如，當TCGA數據集作為訓練集同時合并GEO數據集作為驗證集時，coxlmm的C-index為0.534，coxlasso的C-index為0.527，coxlmm比coxlasso的預測精度提高0.7%，具有更強的預測能力；當GSE31684作為訓練集，同時TCGA和GSE13507合并的數據集作為驗證集時，Cox、coxlasso和coxlmm的C-index分別為0.624、0.632和0.633，coxlmm的預測精度比Cox和coxlasso分別提高0.9%和0.1%，同樣具有更強的預測能力。三種模型在使用兩種外部驗證集的情況下，C-index結果如表2所示。

表2 在驗證集中三種模型的一致性指數(C-index)

3.TCGA膀胱癌中估計PCE和PGE

本研究定義PCE和PGE以揭示臨床和轉錄組信息在膀胱癌生存變異中的不同作用，TCGA中膀胱癌數據集PCE和PGE的估計值分別為14.95%和10.88%，其中PCE的95%置信區間為8.76%～24.35%，PGE的95%置信區間為3.98%～6.42%，PCE和PGE的總和PVE為25.8%。本研究中膀胱癌PCE的估計值略高于PGE的估計值。PVE較大(例如大于10.0%)表明臨床信息和轉錄組信息在生存變異中均起重要作用；而本研究中膀胱癌的PVE為25.8%，結果表明臨床信息和轉錄組信息在膀胱癌生存變異中均起重要作用。

討 論

本文重點研究了三種方法(一般Cox模型、coxlasso模型和coxlmm模型)在TCGA和GEO膀胱癌數據集中對預后的預測效果，并系統地評估了轉錄組的預后價值[15-17]。研究結果顯示，coxlmm模型比其他兩種模型有更好的預測性能，表明全面整合臨床和轉錄組信息可以提高預測的準確性。換言之，在預測模型中納入所有基因比只納入少部分基因具有更好的預測精度。本文利用TCGA數據集進一步評估了coxlmm模型中臨床協變量和基因表達對膀胱癌生存變異的貢獻。先前的研究表明，將除了基因表達之外的其他基因組測量數據結合到臨床協變量中不能獲得顯著的功效增益[44]，所以PCE和PGE能夠涵蓋幾乎所有可用信息解釋的生存變異。當轉錄組信息在PGE量化的變異中占很大比例時，使用基因表達信息的預測模型的準確性將得到顯著提高。需要說明的是，本文的主要目的不是探討如何有效地將多個組學數據整合到膀胱癌預測模型中[6，45-48]。相反，而是利用單個組學數據(即基因表達)通過比較coxlasso和coxlmm預測精度的差異來實現更好的預后預測。此外，將本研究的分析方法用于膀胱癌其他組學數據也是可行的。

本研究存在以下不足，首先，雖然coxlmm模型比其他兩種模型的預測性能更好，但是C-index結果顯示提升的精度并非特別明顯。其次，本研究在分析時，只納入了轉錄組信息而忽略了其他組學數據集(如拷貝數變化和DNA甲基化)可能會影響預測精度。由于在TCGA和GEO數據庫中與膀胱癌有關的臨床信息(如吸煙)缺失嚴重，不能將有效的臨床信息納入模型，這也可能會導致模型的預測性能被低估。再次，膀胱癌的有效樣本量比較少，TCGA和GEO中的截尾數據的占比都比較高[15-17]，所以不可避免地降低了模型預測的準確性，進而會影響預測模型在實踐中應用的效果。另一個不可忽視的問題是，外部驗證的樣本量仍然較少，因此進一步使用更大的樣本量外部驗證模型是必要的。最后，TCGA和GSE31684數據集的種群都是歐洲人群，而GSE13507數據集的種群為亞洲人群，不同種群之間的差異可能也會影響模型預測的準確性。

本研究也有以下幾個優勢。首先，在coxlmm模型中整合全轉錄組數據同時聯合臨床信息對膀胱癌預后進行預測。其次，本研究中多個數據集的交叉驗證以及外部驗證保證了coxlmm模型預測性能的穩定性。第三，coxlmm是一種基于機器學習并且考慮了線性核的預測方法，它為進一步研究如何改善預后評價提供了可能。

我們注意到在其他現有的高維數據統計分析方法中，隨機森林也用于生存結局的預后預測[49]。隨機森林具有抗噪聲、防止過擬合、可處理非線性相關等優點，在基因表達數據判別分類研究中有著較好的分類性能，但大量的噪聲仍會影響隨機森林的分類效能，例如在基因表達差異很大的情況下，隨機森林中的每棵樹的節點很少，極有可能使其他差異基因不被抽到，從而漏掉相當一部分有差異的基因[50]。也即是，當基因表達差異很大時，coxlmm模型的預測效能可能優于隨機森林方法。然而，當生存數據刪失比例較大時，隨機森林和coxlmm模型可能均達不到預期的效果。因此，未來有必要進一步在高維數據背景下比較coxlmm模型與其他方法的優缺點。

最后，我們意識到混合模型(mixed models)和稀疏模型(sparse models)的預測效果優劣取決于對基因與疾病關系的假設與真實情況是否吻合。然而，實際中基因與疾病之間真實的關系往往未知。在今后的研究中，我們可以在遺傳預測的基礎上，結合這兩個模型構建一種混合稀疏預測方法[13，34]，混合稀疏Cox模型有望具有更好的預測性能。

綜上所述，在預測膀胱癌預后時，與僅納入少部分轉錄組學信息相比，納入全轉錄組信息可以提高預測能力。這些也為更深入地了解膀胱癌腫瘤進展機制提供了機會。