基于單細(xì)胞RNA測序的結(jié)直腸癌預(yù)后預(yù)測模型的建立和驗證

馬燕如，季林華，童天穎，嚴(yán)宇青，沈超琴，張昕雨，曹穎穎，洪 潔，陳豪燕

1.上海交通大學(xué)醫(yī)學(xué)院附屬仁濟(jì)醫(yī)院消化科，上海200001；2.上海交通大學(xué)醫(yī)學(xué)院附屬仁濟(jì)醫(yī)院胃腸外科，上海200001

目前，結(jié)直腸癌（colorectal cancer，CRC）是世界上第三大常見的癌癥，高居癌癥相關(guān)死因第二位［1］。約25%的CRC 患者存在轉(zhuǎn)移，而約20%的患者在患病期間就出現(xiàn)遠(yuǎn)處轉(zhuǎn)移［2-3］。轉(zhuǎn)移性CRC的治療對臨床來說仍是一個挑戰(zhàn)，實際治愈率較低［4］。因此，越來越多的研究開始嘗試建立各種模型來完善CRC 的臨床分級和預(yù)后判斷，為臨床決策提供參考［5-7］。

RNA 測序（RNA sequence，RNA-seq）技術(shù)已廣泛應(yīng)用于轉(zhuǎn)錄組分析，用以研究轉(zhuǎn)錄結(jié)構(gòu)、剪接模式、基因表達(dá)差異等［8-9］。然而，普通RNA-seq 方法忽視了組織內(nèi)存在的異質(zhì)性，無法識別不同細(xì)胞群的內(nèi)部特征。近年來，單細(xì)胞RNA 測序（single cell RNA sequence，scRNA-seq）技術(shù)的發(fā)展為根據(jù)腫瘤內(nèi)基因表達(dá)活性分析細(xì)胞類型和狀態(tài)提供了新思路［10-11］，有助于腫瘤細(xì)胞的進(jìn)一步分類。該技術(shù)已經(jīng)成功應(yīng)用于胰腺導(dǎo)管腺癌、乳腺癌和肺癌等［12-13］。此外，scRNA-seq在難治性腫瘤的研究中具有一定的應(yīng)用價值，在監(jiān)測循環(huán)腫瘤細(xì)胞、分析瘤內(nèi)異質(zhì)性和利用其敏感性評估復(fù)發(fā)性腫瘤預(yù)后等方面具有一定的優(yōu)勢［14］。Zhang 等［15］利用scRNA-seq 對患者來源的CRC 樣本進(jìn)行了分析，發(fā)現(xiàn)發(fā)生和未發(fā)生肝轉(zhuǎn)移的腫瘤組織在腫瘤微環(huán)境的信號轉(zhuǎn)導(dǎo)、代謝、免疫調(diào)節(jié)等多個方面的基因表達(dá)上存在差異。

本研究利用GEO 數(shù)據(jù)庫的CRC 患者scRNA-seq 數(shù)據(jù)篩選在轉(zhuǎn)移組織中差異表達(dá)的候選基因，并對候選基因進(jìn)行大樣本分析，以驗證其在預(yù)測CRC 預(yù)后中的意義，為進(jìn)一步的個體化治療提供參考。

1 材料與方法

1.1 數(shù)據(jù)獲取

從美國國家生物技術(shù)信息中心的Gene Expression Omnnibus（GEO）數(shù)據(jù)庫（https：//www.ncbi.nlm.nih.gov/geo/）下載CRC scRNA-seq數(shù)據(jù)集GSE97693。原數(shù)據(jù)集對12 例CRC 患者的腫瘤原發(fā)灶、淋巴結(jié)轉(zhuǎn)移灶和遠(yuǎn)處轉(zhuǎn)移灶樣本進(jìn)行了單細(xì)胞測序，共190個單細(xì)胞數(shù)據(jù)通過了質(zhì)量控制。本研究提取了其中淋巴結(jié)轉(zhuǎn)移灶和相應(yīng)的腫瘤原發(fā)灶數(shù)據(jù)以供后續(xù)分析。從癌癥基因組圖譜（The Cancer Genome Atlas，TCGA）數(shù)據(jù)庫下載588例CRC樣本轉(zhuǎn)錄組表達(dá)譜以及對應(yīng)的臨床信息，包括年齡、性別、腫瘤原發(fā)灶分級（T stage）、是否有淋巴結(jié)轉(zhuǎn)移（N stage）、是否發(fā)生遠(yuǎn)處轉(zhuǎn)移（M stage）、無進(jìn)展生存期（progression free survival，PFS）和無進(jìn)展生存狀態(tài)共7個指標(biāo)。

1.2 數(shù)據(jù)處理

提取下載的scRNA-seq 數(shù)據(jù)，通過Scanpy 包進(jìn)行數(shù)據(jù)分析［16］，參考基因組為GRCh38。利用統(tǒng)一流形逼近與投影（uniform manifold approximation and projection，UMAP） 算法對scRNA-seq 數(shù)據(jù)進(jìn)行降維分析，采用Louvain 算法進(jìn)行聚類分析。以標(biāo)準(zhǔn)化P=0.05 為截止標(biāo)準(zhǔn)，再根據(jù)基因表達(dá)比值對數(shù)的絕對值（|log fold change|，|log FC|）由高到低排序，分別篩選出轉(zhuǎn)移組和未轉(zhuǎn)移組前15 個差異表達(dá)基因，組成候選基因集，結(jié)果以熱圖的形式展現(xiàn)。

1.3 CRC患者中標(biāo)記基因的篩選和鑒定

將來自TCGA 數(shù)據(jù)庫的588 例CRC 患者表達(dá)譜隨機(jī)分為訓(xùn)練集和驗證集（49.8%∶50.2%）。提取訓(xùn)練集中293 例患者標(biāo)記基因的轉(zhuǎn)錄組圖譜。在scRNA-seq 獲得的候選基因集中，使用glmnet包選取其中9個與預(yù)后相關(guān)的關(guān)鍵基因，以是否復(fù)發(fā)為因變量，通過Logistic 回歸分析的方法建立套索回歸算法（the least absolute shrinkage and selection operator，LASSO）回歸模型。根據(jù)建立的回歸模型對每例患者的關(guān)鍵基因表達(dá)進(jìn)行評分，評分=-log10［-?（βi×Expi）］；其中βi為系數(shù)，代表賦予的每個關(guān)鍵基因表達(dá)的加權(quán)。隨后，分別在訓(xùn)練集和驗證集中對復(fù)發(fā)和未復(fù)發(fā)患者的關(guān)鍵基因評分利用Wilcoxon檢驗進(jìn)行比較；并通過Kaplan-Meier 分析評估訓(xùn)練集、驗證集和總TCGA集預(yù)后結(jié)局的差異。

1.4 預(yù)后預(yù)測模型的建立和評估

合并總TCGA 集的評分和其他臨床變量。采用Logistic 回歸分析對重要的臨床變量進(jìn)行評估。排除無統(tǒng)計學(xué)意義的變量后，利用廣義線性模型（generalized linear model，GLM）建立完整的模型。采用具有曲線下面積（area under curve，AUC） 的受試者工作特征（receiver operating characteristic，ROC）曲線評估列線圖的實際預(yù)測價值。為了評估模型的臨床使用價值，通過對總TCGA 集中不同閾值概率下患者的凈獲益進(jìn)行量化，采用決策曲線分析（decision curve analysis，DCA）確定列線圖的臨床有效性。

1.5 生物信息學(xué)分析和統(tǒng)計學(xué)分析

采用glmnet 和survival 包進(jìn)行LASSO 回歸、Logistic回歸和Kaplan-Meier 分析。GLM 和列線圖用rms 包建立。使用Wilcoxon 秩和檢驗判斷評分在復(fù)發(fā)和未復(fù)發(fā)患者中的差異是否具有統(tǒng)計學(xué)意義。所有統(tǒng)計分析均在R studio（3.6.1版）軟件中進(jìn)行，P<0.05表示差異有統(tǒng)計學(xué)意義。

2 結(jié)果

2.1 在單細(xì)胞水平篩選CRC預(yù)后相關(guān)基因

為從單細(xì)胞水平篩選CRC 預(yù)后相關(guān)基因，從GEO 數(shù)據(jù)庫獲取CRC樣本scRNA-seq數(shù)據(jù)集，通過UMAP算法，對scRNA-seq 數(shù)據(jù)進(jìn)行降維分析，結(jié)果顯示淋巴結(jié)轉(zhuǎn)移灶與腫瘤原發(fā)病灶來源的細(xì)胞群分布具有較大差異（圖1A）。進(jìn)一步采用Louvain算法進(jìn)行聚類分析，將scRNAseq數(shù)據(jù)分為2個亞群，其中腫瘤轉(zhuǎn)移相關(guān)基因鈣黏蛋白1（cadherin 1，CDH1）的表達(dá)在細(xì)胞群0和細(xì)胞群1中具有顯著差異（圖1B、C）。研究［17］表明，CDH1突變與細(xì)胞黏附減少而增殖、侵襲、轉(zhuǎn)移增加相關(guān)；因此，本研究將細(xì)胞群0 和1 分別定義為未轉(zhuǎn)移組和轉(zhuǎn)移組，以進(jìn)一步分析CRC 轉(zhuǎn)移相關(guān)基因（圖1D）。分別在轉(zhuǎn)移組和未轉(zhuǎn)移組選擇具有統(tǒng)計學(xué)意義（P<0.05）的前15 個差異表達(dá)基因作為后續(xù)研究的候選基因（圖1E）。

圖1 單細(xì)胞水平篩選CRC預(yù)后相關(guān)基因Fig 1 Screening genes associated with prognosis of CRC from single-cell level

2.2 對關(guān)鍵基因的評分和驗證

在訓(xùn)練集人群中，采用LASSO 回歸從候選基因集選擇了9個與預(yù)后相關(guān)的關(guān)鍵基因（圖2），分別是膜聯(lián)蛋白A13（annexin A13，ANXA13）、凋亡蛋白酶1（caspase 1，CASP1）、卵黃樣羧肽酶（carboxypeptidase vitellogenic like，CPVL）、胱抑素7（cystatin-7，CST7）、長鏈非編碼RNA H19、嗅覺受體家族7 亞家族C 成員2（olfactory receptor family 7 subfamily C member 2，OR7C2）、信號序列受體3（signal sequence receptor subunit 3，SSR3）、WD重 復(fù) 結(jié) 構(gòu) 域26 （WD repeat domain 26，WDR26） 和WDR62。并根據(jù)建立的回歸模型對每例患者的關(guān)鍵基因表達(dá)情況進(jìn)行評分。CRC患者的年齡、性別、腫瘤原發(fā)灶分級、是否有淋巴結(jié)轉(zhuǎn)移、是否發(fā)生遠(yuǎn)處轉(zhuǎn)移、關(guān)鍵基因評分（score）在訓(xùn)練集和驗證集中的差異均無統(tǒng)計學(xué)意義（表1）。

圖2 采用LASSO回歸選擇關(guān)鍵基因Fig 2 Texture feature selection using LASSO binary Logistic regression model

分別在訓(xùn)練集和驗證集中對復(fù)發(fā)和未復(fù)發(fā)患者的評分進(jìn)行比較，差異均有統(tǒng)計學(xué)意義（圖3A、B，P值分別為0.000 和0.002）。根據(jù)評分的中位數(shù)，分別在訓(xùn)練集、驗證集和總TCGA 集中將患者分為高分組和低分組，再進(jìn)行Kaplan-Meier 分析。結(jié)果表明，訓(xùn)練集中評分高的患者PFS 結(jié)局明顯較差（圖3C，P=0.002），并在隨訪時間>1 個月的驗證集和總TCGA 集中得到了一致驗證（圖3D、E，P值分別為0.017和0.000）。

表1 CRC患者的臨床特征Tab 1 Clinical characteristics of patients with CRC

2.3 構(gòu)建預(yù)后預(yù)測模型

將評分與其他臨床變量相結(jié)合，構(gòu)建完整的預(yù)測CRC 患者預(yù)后的模型。排除單因素Logistic回歸分析中無統(tǒng)計學(xué)意義的年齡和性別變量（表2，P 值分別為0.406和0.779）以及在多變量Logistic 回歸模型中無統(tǒng)計學(xué)意義的淋巴結(jié)轉(zhuǎn)移分級變量（N stage）。選擇是否發(fā)生轉(zhuǎn)移、腫瘤原發(fā)灶分級和關(guān)鍵基因評分3個獨(dú)立變量作為最終的模型變量。利用GLM 回歸算法建立了包含這3 個變量的列線圖（圖4A）。列線圖預(yù)測預(yù)后結(jié)局的AUC 分別達(dá)到0.775和0.705（圖4B、C）。

表2 訓(xùn)練集中評分和臨床候選預(yù)測變量的單因素Logistic回歸分析Tab 2 Univariate Logistic regression analysis of scores and clinical candidate predictors in training set

2.4 模型的臨床應(yīng)用價值

評分和評分-臨床變量整合模型的DCA 曲線如圖5 所示。該曲線表明，如果復(fù)發(fā)風(fēng)險閾值在4%～60%，使用評分-臨床變量整合模型預(yù)測預(yù)后比治療所有患者的方案或無治療方案更有益；且在此范圍內(nèi)，使用整合模型預(yù)測的效果優(yōu)于單純使用評分。

圖3 關(guān)鍵基因評分驗證Fig 3 Hub gene score validation

圖4 基于評分和部分臨床數(shù)據(jù)預(yù)測復(fù)發(fā)的列線圖Fig 4 Nomogram for predicting prognosis based on the score and some clinical variables

圖5 列線圖和關(guān)鍵基因評分的DCAFig 5 DCA of the nomogram and the hub gene scores

3 討論

本研究通過GEO 數(shù)據(jù)庫獲取CRC 的scRNA-seq 數(shù)據(jù)，從單細(xì)胞水平探究CRC 預(yù)后相關(guān)的基因組特征。既往研究表明，CDH1 基因編碼產(chǎn)物鈣黏蛋白（E-cadherin）是一種鈣依賴的細(xì)胞間黏附分子和抑癌蛋白，通過細(xì)胞間黏附復(fù)合物在上皮細(xì)胞形成和分化中起關(guān)鍵作用［17-18］。CDH1的突變與腫瘤的浸潤和轉(zhuǎn)移能力相關(guān)，黏附復(fù)合物的破壞使細(xì)胞間黏附喪失，細(xì)胞活性增加［18］?；诖?，本研究根據(jù)CDH1 的表達(dá)將細(xì)胞群分為轉(zhuǎn)移和未轉(zhuǎn)移2個群。分別在轉(zhuǎn)移組和未轉(zhuǎn)移組選擇具有統(tǒng)計學(xué)意義（P<0.05）的前15個差異表達(dá)基因作為后續(xù)研究的候選基因，利用LASSO 回歸算法進(jìn)行篩選，確認(rèn)了9 個關(guān)鍵基因，其中部分基因在惡性腫瘤的進(jìn)展中起重要作用。與未發(fā)生轉(zhuǎn)移的CRC患者相比，SSR3在發(fā)生轉(zhuǎn)移的CRC患者中顯著升高，提示其可能與CRC 的疾病進(jìn)展相關(guān)［19］。在急性淋巴細(xì)胞白血病中，CASP1 的過表達(dá)導(dǎo)致糖皮質(zhì)激素受體分裂，糖皮質(zhì)激素誘導(dǎo)的轉(zhuǎn)錄反應(yīng)減弱，糖皮質(zhì)激素耐藥性增強(qiáng)［20］。H19 是一個長鏈非編碼RNA，其可以通過H19/S-腺苷高半胱氨酸水解酶/DNA 甲基轉(zhuǎn)移酶3B（H19/SAHH/DNMT3B）軸導(dǎo)致細(xì)胞自噬激活，從而引起乳腺癌患者對他莫昔芬治療產(chǎn)生耐藥［21］。CST7則在腫瘤逃逸和耐受、肝癌晚期復(fù)發(fā)、肝癌進(jìn)展和某些肝癌亞類中顯著富集［22］。這為后續(xù)的研究提供了一個方向。

根據(jù)建立的回歸模型對每例患者的關(guān)鍵基因表達(dá)進(jìn)行評分后，利用隨機(jī)分成的驗證集檢驗評分在臨床的應(yīng)用價值，發(fā)現(xiàn)復(fù)發(fā)組評分顯著高于未復(fù)發(fā)組，這與生存分析結(jié)果是一致的。至于其是否具有潛在的治療預(yù)測價值尚不清楚，這將是另一個有意義的研究方向。隨后采用Logistic 回歸分析，排除無統(tǒng)計學(xué)意義的年齡、性別和淋巴結(jié)轉(zhuǎn)移分級變量，選擇是否發(fā)生轉(zhuǎn)移、腫瘤原發(fā)灶分級和關(guān)鍵基因評分3個獨(dú)立的變量作為最終的整合模型變量。ROC 曲線顯示整合模型對CRC 患者的預(yù)后具有較高的預(yù)測價值，且臨床也可以據(jù)此對患者的個體化治療方案進(jìn)行及時調(diào)整和完善。

本研究將GEO 數(shù)據(jù)庫獲得的scRNA-seq 數(shù)據(jù)分析與TCGA 數(shù)據(jù)庫中的大樣本人群數(shù)據(jù)驗證相結(jié)合。與CRC傳統(tǒng)轉(zhuǎn)錄組測序分析相比［23-24］，傳統(tǒng)RNA-seq 可能會覆蓋潛在的重要標(biāo)記基因，而scRNA-seq 在這方面具有一定的優(yōu)勢。本研究仍存在以下局限性：首先，下載的TCGA數(shù)據(jù)大多來自美國或歐洲人群，而這是否適用于亞洲人種尚不清楚，因此，該研究結(jié)論需在中國人群中進(jìn)一步驗證；其次，雖然該整合模型在大量CRC 人群中得到了很好的驗證，但仍需補(bǔ)充基礎(chǔ)實驗來揭示其促進(jìn)腫瘤發(fā)展的具體機(jī)制。

綜上所述，本研究基于scRNA-seq 技術(shù)篩選標(biāo)記基因，并對標(biāo)記基因進(jìn)行大樣本分析，對預(yù)測臨床CRC 患者預(yù)后、優(yōu)化治療方案具有一定的指導(dǎo)意義。