基于生物信息學的急性髓系白血病患者預后風險模型的構建與驗證

張斯娜，馬書杰

1.淄博市疾病預防控制中心性病艾滋病防制所（山東淄博 255020）

2.河南理工大學醫學院（河南焦作 454003）

急性髓系白血病（acute myeloid leukemia,AML）是一種未成熟髓系造血干細胞異常克隆性惡性腫瘤，具有高度異質性，其主要特征是克隆性造血干細胞或祖細胞（hematopoietic stem and progenitor cells, HSPCs）分化和凋亡障礙。AML可導致克隆性HSPCs 在骨髓中惡性增殖和聚集，從而影響骨髓正常造血功能[1-2]。AML 的高度異質性與多種因素有關，如基因組突變、表觀遺傳學異常和基因表達異常等[3]。全球范圍內，AML發病人數從1990 年的63 840 例增加到2017 年的119 570 例，發病率在過去的28 年間增加了87.3%[4]。化療、生物制劑和干細胞移植是AML的主要治療方式[5-7]。盡管AML 的治療取得了重大進展，但化療藥物毒性可能導致急性和危及生命的并發癥。AML 通常預后不良，患者中位生存時間僅為8.5 個月，2 年和5 年總生存率（overall survival, OS）均低于35%，僅35%～40%的年輕患者（＜60 歲）和5%～15%的老年患者（≥60 歲）在化療后可存活5 年以上[8-9]。隨著分子生物學的迅速發展，基因組學、蛋白質組學和生物信息學分析方法已被用于開發新的個性化治療策略，研究相關生物分子的功能，收集臨床數據基因組匹配的新趨勢信息是改善AML 患者預后的有效方法[10-11]。盡管許多研究分析了AML 的基因組變異，但基因組變異與AML 分子機制間的關聯仍不明確。因此，需要進一步探索更有效的AML早期診斷和預后方法。AML 的基因表達譜已被證明在診斷不同的細胞遺傳學亞型、發現新的AML亞類和預測預后方面具有重要價值[12-13]。基因組分析是篩選潛在的疾病診斷和預后生物標志物的一種有效方法，本研究基于公共數據庫的基因表達譜數據集構建并驗證AML 患者預后風險模型，識別AML 患者的預后生物標志物。

1 資料與方法

1.1 數據來源與差異表達基因篩選

基于癌癥基因組圖譜（The Cancer Genome Atla, TCGA）數據庫（https://gdc.xenahubs.net/download/TCGA-LAML/Xena_Matrices/TCGALAML.htseq_fpkm.tsv.gz），獲取AML 基因數據（151 例）作為病例組。基于基因型-組織表達（Genotype-Tissue Expression, GTEx） 數據庫（https://toil.xenahubs.net/download/TCGatargettex_rsem_gene_fpkm.gz），獲取正常骨髓樣本（70 例）作為對照組。利用Wilcoxon 秩和檢驗分析AML樣本基因數據，篩選|logFC|＞2 且錯誤發現率（false discovery rate, FDR）＜0.05 的差異表達基因（differentially expressed genes, DEGs）。采用單因素Cox 回歸分析，篩選DEGs 中與AML 生存相關基因，以P＜0.01 為差異有統計學意義。

1.2 風險評分模型構建與特征基因驗證

151 例AML 樣本中，剔除臨床特征信息缺失和隨訪時間為0 的樣本后，最終納入132 例，將132 例AML 樣本按照1 ∶ 1 的比例隨機分成兩組，分別作為訓練集和測試集。在訓練集中，通過單變量Cox 比例風險回歸分析mRNA 的表達與患者OS 之間的關系。在構建風險預測模型前，采用隨機生存森林算法對AML 預后相關的差異基因進行篩選。在隨機生存森林監督分類算法中，采用迭代過程縮小基因集。隨機森林的樹數（ntree）為1 000，將每一步輸入節點數的平方根設為單個分類樹的每個節點隨機選取的預后相關的差異mRNA 的大小，對袋外數據樣本進行泛化誤差估計。為了選擇極少最重要的基因變量，可依據變量的重要性（per-mutation variable important,VIMP）對變量進行篩選，VIMP 值越大表明其預測能力越強。風險預測評分模型是基于AML 預后相關的差異基因而構建，通過其估計回歸系數加權如下：

其中，N 為預后相關的差異基因數量，Expi為預后相關的差異基因的表達值，Coei為單變量Cox 回歸分析中預后相關差異基因的估計回歸系數。模型構建后，模型中涉及的特征基因通過以下方式進一步驗證：①對模型高低風險組的生存分析，最佳的截斷表達值由“survminer”R 包的“surv_cutpoint”函數確定。使用Kaplan-Meier和log-rank 方法在訓練集中進行生存分析，再使用AML 數據和臨床信息在測試集中驗證該模型。②在訓練集和測試集中使用“survivalROC”R 軟件包進行時間依賴的受試者工作特征（receiver operating characteristic, ROC）曲線分析，并計算1 年、2 年和3 年AML 患者生存率的ROC 曲線下面積（area under the curve, AUC）。通過決策曲線分析量化不同時間概率下的凈收益，從而確定風險評分的臨床實用性，決策曲線分析通過“rmda”R 軟件包進行。基于所有樣本的基因表達，使用“Rtsne”R 軟件包進行t-SNE 分析，探索不同組的分布。

1.3 分層分析

在全集中，對AML 患者總生存率顯著相關的臨床特征數據進行分層分析，其中年齡以60歲為截斷值，≥60 歲為老年組，＜60 歲為年輕組。

1.4 功能富集分析

為了闡明與風險評分相關的生物學功能和途徑，高風險組和低風險組之間的DEGs 被用于進行基因本體（Gene Ontology, GO）富集分析。利用“clusterProfiler”R 包進行分析，以校正后P值小于0.05 為差異有統計學意義。

2 結果

2.1 AML患者臨床特征

訓練集、測試集和全集中AML 患者的臨床特征見表1。

表1 AML患者的臨床特征Table 1.Clinical characteristics of patients with AML

2.2 差異表達基因鑒別

基于AML 樣本和正常骨髓樣本數據，共獲取19 022 個基因的mRNA 測序數據，通過基因差異顯著性分析篩選出2 270 個DEGs（圖1）。通過單變量Cox 回歸分析，得到335 個AML 預后相關特征基因（P＜0.01）。

圖1 差異基因表達熱圖Figure 1.Heatmap of differential genes expression

2.3 特征基因篩選

為了進一步縮小預后基因的數量，通過隨機森林算法對335 個與AML 患者OS 顯著相關的DEGs 進行分析，在每一步中，根據通過置換測試使用袋外樣品對每個基因的重要分數進行估計，丟棄不重要的基因，并根據每一步中的排列重要評分，得到與AML 患者OS 顯著相關的5 個基因標志物 （ACOT7、SFXN3、SLC29A2、HINT2和ZMAT1），見圖2-A。

圖2 特征基因風險評分分析Figure 2.Risk scoring analysis of characteristic genes

風險評分模型構建如下：風險評分=（-0.959 24×ZMAT1表達值）+（0.525 21×SFXN3表達值）+（1.060 22×HINT2表達值）+（0.669 68×SLC29A2表達值）+（0.079 19×ACOT7表達值）。基于風險評分模型計算訓練集每個AML 患者的風險評分，根據評分將患者分為高風險組（n=24）和低風險組（n=42），使用cutoff=19.686 作為風險臨界值。圖2-B、圖2-C、圖2-D 展示了66 例AML 患者的預后評分、生存狀態和腫瘤mRNA 表達的分布，根據5 個特征基因的預后評分值進行排序。在5 個基因中，4 個基因（ACOT7、SFXN3、SLC29A2和HINT2）與預后高風險相關 [風險比（hazard ratio, HR）＞1]，1 個基因（ZMAT1）為保護性因素（HR ＜1）。預后評分高的基因傾向于表達高風險的mRNA，而預后評分低的基因傾向于表達保護性mRNA，見圖2-C 和圖2-D。Kaplan-Meier 生存分析和對數秩檢驗表明，高風險組和低風險組的總體生存率有顯著差異，高風險組患者總生存期明顯短于低風險組（中位生存期3.79 年 vs.0.57 年，P＜0.001），見圖2-E。5 個基因生存預測的時間依賴ROC 曲線分析顯示，1 年、2 年、3 年的AUC 值分別為0.778、0.743、0.677，曲線接近左上角，具有良好的敏感性和特異性，表明該預測模型具有較高的準確率，見圖2-F。

2.4 特征基因驗證

為了測試5 個特征基因在預測AML 患者總生存期中的預后價值，在測試集中測試5 個特征基因。通過使用相同的風險評分模型，測試集的66 例患者被分為高風險組（n=33）和低風險組（n=33）。結果顯示，高風險組患者的總生存期明顯短于低風險組患者（中位生存期0.71 年vs.2.5 年，P＜0.001），見圖3-A。時間依賴的ROC 曲線分析對5 個特征基因的生存預測的AUC值在1 年時達到0.775，2 年時達到0.746，3 年時達到0.696，見圖3-B。進一步將5 個特征基因應用于整個數據集的患者，以驗證其預測價值，來自訓練集的相同風險評分模型和中位風險截止標準將全集的132 例AML 患者分為高風險組（n=66）和低風險組（n=66）。高風險組患者的總生存期明顯短于低風險組患者（中位生存期0.68 年 vs.2.58 年，P＜0.001），見圖3-C。在132 例AML患者全集中驗證5 個特征基因的ROC 曲線在1 年、2 年、3 年的AUC 值分別為0.836、0.777、0.694，見圖3-D。

2.5 特征基因的預后價值與臨床變量的獨立性

如表2 所示，多變量Cox 回歸分析結果表明，經年齡、性別和各分子分析異常測試結果調整后，5 個特征基因仍與AML 患者總生存率存在顯著相關性。當控制其他臨床變量時，在訓練集中，高風險組相對于低風險組總生存率的風險比為3.47[95%CI（1.80，6.70），P＜ 0.001]，在測試集和全集中，高風險組相對于低風險組總生存率的風險比分別為3.25 [95%CI（1.71, 6.19），P＜0.001]和3.39 [95%CI（2.12，5.43），P＜0.001]。年齡與三組數據集AML 患者的總生存率顯著相關，因此根據年齡進行分層分析，所有患者被分為年輕組（n=77）和老年組（n=55），再根據風險評分模型，將年輕組和老年組患者分為高風險組和低風險組。如圖4-A 和圖4-B 所示，生存分析表明，在年輕組和老年組中，高風險組的總體生存時間明顯短于低風險組（P＜0.001）。

圖4 兩組總體生存率的估計Figure 4.Estimation of overall survival rates for two groups

表2 單變量和多變量Cox回歸分析結果Table 2.Results of univariate and multivariate Cox regression analysis

2.6 特征基因的性能評估

t-SNE 分析表明，不同風險組的患者分布在兩個方向（圖5-A）。5 個特征基因組合風險評分的決策曲線分析表明，風險評分預測AML 患者1 年、2 年和3 年的生存情況表現出比極端曲線更高的凈獲益，而1 年的閾值范圍更大，認為1 年風險評分表現更優，見圖5-B。

圖5 特征基因的性能評估Figure 5.Performance evaluation of characteristic genes

2.7 功能富集分析

如圖6-A 所示，在整個數據集中，DEGs 富集幾種與代謝過程相關的途徑，如固醇代謝過程、膽固醇代謝過程、類固醇代謝過程、仲醇代謝過程、脂質代謝過程和酒精代謝過程（校正后P＜0.05）；在許多免疫相關的生物過程中也明顯富集，如白細胞遷移、T 細胞分化和巨噬細胞遷移（校正后P＜0.05）；DEGs 在凋亡相關的生物過程中也明顯富集，如凋亡信號通路的負調控、線粒體凋亡途徑和線粒體釋放細胞色素C途徑（校正后P＜0.05）。圖6-B 展示了各途徑所富集到的基因。

圖6 GO富集分析Figure 6.GO enrichment analysis

3 討論

在最近的研究中，許多AML 相關的癌基因和腫瘤抑制基因已被確定。例如，Zhou 等發現HSPG2過表達與AML 患者的不良預后相關，AML 患者中HSPG2的表達水平顯著上調，完全緩解后則下調，而在復發時再次升高[14]。Li 等研究表明，SLC38A1是AML 的重要預后和預測指標，SLC38A1高表達患者的NPM1基因突變率較低，不良風險核型發生率較高，SLC38A1表達水平高的患者總生存期顯著縮短[15]。Wróbel 等發現IL-17F多態性與AML 的誘導密切相關[16]。盡管發現了一些與AML 預后和診斷相關的基因，但其在AML 中的作用仍有待探究。AML 的預后部分由遺傳因素驅動，多種基因的結合有助于提高預后預測的準確性。本研究中，AML 患者的臨床數據和基因表達數據來自TCGA 數據庫，以鑒定與預后相關的基因。通過單變量Cox 分析和隨機生存森林分析選擇生存相關特征基因，最終通過多變量Cox 回歸分析構建了基于5 個預后相關特征基因的風險評分模型。

在本研究篩選出的5 個特征基因中，僅ZMAT1被認為是低風險且與AML預后相關基因。目前尚未發現ZMAT1在AML 中的研究報告，但在胰腺導管腺癌中有過報道，研究發現ZMAT1的下調與胰腺導管腺癌不良的臨床病理特征和低生存率相關，ZMAT1通過誘導SIRT3/p53信號通路在胰腺導管腺癌中起到抑癌因子的作用，其中ZMAT1促進SIRT3基因的轉錄，隨后上調p53的表達，并參與調節胰腺導管腺癌的細胞周期進程和凋亡[17]。Zhang 等研究發現高水平的ACOT7對AML 患者的預后有顯著影響，ACOT7高表達組的無事件生存期（event-free survival, EFS）和OS顯著降低[18]。ACOT7的表達水平是影響AML 患者預后的獨立因素，ACOT7的高表達水平在男性患者中更常見，且表達水平與FAB 亞型以及細胞遺傳學之間可能存在潛在相關性[18]。SLC29A2也是ENT2，有研究表明肝細胞癌（hepatocellular carcinoma, HCC）組織中SLC29A2的上調與晚期、血管浸潤與患者生存率不良顯著相關[19]。有研究推測SLC29A2的異常上調可能改變細胞核苷酸合成和核苷酸池平衡，從而導致癌癥基因組不穩定，進而為腫瘤提供生長優勢[20]。由于SLC29A2是一種已知的腺苷轉運體，而從缺氧組織中釋放的腺苷刺激血管內皮生長因子的釋放，再結合內皮細胞上的受體，刺激內皮細胞增殖、遷移和腫瘤血管生成[20]。SFXN3和HINT2已被證明分別在口腔鱗狀細胞癌和肝細胞癌中作為癌基因[21-22]。Murase 等認為口腔鱗狀細胞癌患者血清抗SFXN3自身抗體水平顯著升高，血清抗SFXN3-autoAb水平與原發腫瘤大小輕度相關，但在早期癌癥中水平稍高，SFXN3蛋白在癌巢之間的基質成纖維細胞中表達，也在鱗狀上皮的基底層中表達，治療后血清抗SFXN3自身抗體水平的變化與臨床腫瘤負荷相關[21]。Zhou 等研究表明HINT2在HCC組織中的表達明顯低于鄰近的非腫瘤組織，且HINT2低表達HCC 患者復發的可能性更高，可作為HCC 復發的一個預后指標[22]。鑒于在其他類型癌癥中的重要作用，這些特征基因可用于預測AML 患者的生存率。基于5 個特征基因的風險評分模型，低風險和高風險患者的存活率有顯著差異，高風險患者比低風險患者存活時間明顯更短。

基于不同風險組間的DEGs 進行了GO 富集分析，發現許多免疫相關的生物過程和途徑得到了富集，說明AML 可能與腫瘤免疫有密切聯系的推測具有一定合理性。基于5 個特征基因的風險評分模型初步用于AML 患者的生存預測，并進一步證明了這些特征基因在AML 生存中的重要作用，可作為更好地評估AML 患者生存和預后的新方法。此外，風險評分模型在測試集及整個數據集中得到進一步驗證，表明該模型在AML患者生存預測中具有良好的可重復性。由于流行病學的限制，本研究未能分析風險評分與AML患者生存率之間的具體關系，未來有待進一步深入探究。

綜上所述，本研究確定了5 個與AML 預后密切相關的特征基因，并構建和驗證了基于5 個特征基因的風險評分模型，該模型可用于預測AML 患者的生存期，為AML 患者的個性化風險評估、治療和預后提供參考。