單細(xì)胞轉(zhuǎn)錄組測序分析結(jié)腸癌微環(huán)境中免疫細(xì)胞亞群與癌癥進(jìn)程的關(guān)系

周彥彤，李春曉，王勁松，孫芳洲，許東奎，張雪燕，袁芃，錢海利*

1國家癌癥中心/國家腫瘤臨床醫(yī)學(xué)研究中心/中國醫(yī)學(xué)科學(xué)院北京協(xié)和醫(yī)學(xué)院腫瘤醫(yī)院分子腫瘤學(xué)國家重點(diǎn)實(shí)驗(yàn)室，北京 100021；2國家癌癥中心/國家腫瘤臨床醫(yī)學(xué)研究中心/中國醫(yī)學(xué)科學(xué)院北京協(xié)和醫(yī)學(xué)院腫瘤醫(yī)院腫瘤結(jié)直腸外科，北京 100021；3國家癌癥中心/國家腫瘤臨床醫(yī)學(xué)研究中心/中國醫(yī)學(xué)科學(xué)院北京協(xié)和醫(yī)學(xué)院腫瘤醫(yī)院腫瘤特需醫(yī)療部，北京 100021

根據(jù)2020年世界衛(wèi)生組織(WHO)發(fā)布的《2020年世界癌癥進(jìn)展報(bào)告》，結(jié)直腸癌發(fā)病率在世界范圍內(nèi)居第3位，病死率居第2位，是嚴(yán)重威脅人類健康的疾病[1]。為了攻克癌癥，新的治療技術(shù)層出不窮，其中免疫治療展現(xiàn)出不俗的治療效果，于2013年被Science雜志評為該年度最重要的科學(xué)突破。為了使免疫治療這一手段讓更多的患者受益，針對腫瘤微環(huán)境中的免疫細(xì)胞與腫瘤細(xì)胞的相互調(diào)控及其影響腫瘤細(xì)胞行為的研究成為熱點(diǎn)[2]。有研究顯示，腫瘤微環(huán)境中的免疫細(xì)胞具有多種協(xié)助腫瘤發(fā)生發(fā)展的生物學(xué)行為，包括參與腫瘤局部血管生成、協(xié)助腫瘤細(xì)胞發(fā)生轉(zhuǎn)移及免疫逃逸等[3-6]。但目前在探究免疫細(xì)胞與腫瘤細(xì)胞相關(guān)調(diào)控行為的研究中，對免疫細(xì)胞表型的劃分較為籠統(tǒng)，雖然基于單細(xì)胞測序及分析技術(shù)的提升，已經(jīng)能在結(jié)直腸癌免疫微環(huán)境中識別更細(xì)分的免疫亞群，但這些免疫亞群的分布及其在大樣本臨床隊(duì)列中的潛在價值尚不清楚。CIBERSORTx算法是一種機(jī)器學(xué)習(xí)方法，可用于估計(jì)細(xì)胞比例及推斷細(xì)胞類型特異性基因表達(dá)譜[7-8]。本研究以單細(xì)胞轉(zhuǎn)錄組測序(scRNAseq)結(jié)果為參考，利用CIBERSORTx算法構(gòu)建特征表達(dá)譜，測算癌癥基因組圖譜(TCGA)結(jié)腸癌隊(duì)列中腫瘤的免疫細(xì)胞浸潤情況，并評價其對腫瘤患者預(yù)后的預(yù)測價值。

1 資料與方法

1.1資料數(shù)據(jù)來源 通過基因表達(dá)數(shù)據(jù)庫(GEO)獲得GSE146771[9]中43 817個單細(xì)胞轉(zhuǎn)錄組，分析CD45+細(xì)胞樣本的每千堿基記錄本(transcripts per kilobase million，TPM)數(shù)據(jù)，調(diào)整樣本特征(nFeatures)數(shù)值，進(jìn)一步篩選獲得更為勻質(zhì)的13 193個細(xì)胞。通過10×測序平臺，獲得這些細(xì)胞的轉(zhuǎn)錄組數(shù)據(jù)。結(jié)腸癌單細(xì)胞轉(zhuǎn)錄組測序數(shù)據(jù)包含來自9例結(jié)腸癌患者的外周血、術(shù)中切除的正常腸黏膜組織、結(jié)腸腺瘤及腺癌組織的單細(xì)胞測序基因表達(dá)譜；從加州大學(xué)圣克魯茲分校癌癥功能基因組數(shù)據(jù)庫(UCSC XENA)獲取TCGA結(jié)直腸癌隊(duì)列，該隊(duì)列包含臨床結(jié)直腸癌樣本513例，收集其臨床資料如年齡、性別、TNM分期、隨訪時間及生存狀態(tài)等，并將每千個堿基的轉(zhuǎn)錄每百萬映射讀取片段(fragments per kilobase per million，F(xiàn)PKM)值的轉(zhuǎn)錄組測序(RNAseq)數(shù)據(jù)轉(zhuǎn)換為TPM數(shù)據(jù)，以便與scRNA-seq組數(shù)據(jù)保持一致。

1.2建立scRNA-seq特征矩陣 使用CIBERSORTx在線分析平臺(https://cibersortx.stanford.edu/)，在不進(jìn)行物理分離的情況下，分析細(xì)胞類型特異性基因表達(dá)譜。按照指令上傳單細(xì)胞表達(dá)矩陣，在默認(rèn)參數(shù)保持不變的條件下，用CIBERSORTx計(jì)算獲得38個免疫細(xì)胞類型的簽名矩陣。

1.3CIBERSORTx計(jì)算細(xì)胞組分 依據(jù)CIBERSORTx的提示，上傳TCGA及GEO，獲得結(jié)直腸癌表達(dá)譜數(shù)據(jù)集，選取之前以scRNA-seq為參考數(shù)據(jù)得到的簽名矩陣，選擇S模式進(jìn)行批量校正。將排列數(shù)設(shè)置為500，其他參數(shù)保持默認(rèn)值。運(yùn)行CIBERSORTx后，獲得了每個樣本中38個腫瘤浸潤免疫細(xì)胞亞群的相對比例，將P＜0.05的樣本納入進(jìn)一步研究。

1.4免疫細(xì)胞亞群與轉(zhuǎn)移的相關(guān)性分析 為了探究免疫細(xì)胞亞群與腫瘤轉(zhuǎn)移的關(guān)系，將513例TCGA結(jié)直腸癌隊(duì)列分為轉(zhuǎn)移組(216例)與未轉(zhuǎn)移組(297例)，比較兩組免疫細(xì)胞亞群分?jǐn)?shù)的差異，并分析免疫細(xì)胞亞群與淋巴結(jié)轉(zhuǎn)移數(shù)的相關(guān)性。

1.5基于免疫細(xì)胞亞群的生存分析 為探究不同免疫細(xì)胞亞群對結(jié)腸癌預(yù)后的影響，使用Kaplan-Meier法繪制患者的生存曲線，用log-rankχ2檢驗(yàn)分析患者總生存期與免疫細(xì)胞亞群細(xì)胞分?jǐn)?shù)的關(guān)系。

1.6免疫細(xì)胞亞群功能分析 采用功能富集分析數(shù)據(jù)庫(the database for annotation, visualization and integrated discovery，DAVID)以及單細(xì)胞轉(zhuǎn)錄組分析和注釋工具(single-cell transcriptome analysis of pathway activation signatures，scTPA)進(jìn)行功能聚類。將CIBERSORTx分析獲得的CD4+T亞群的評分與其對應(yīng)的表達(dá)譜的基因列表進(jìn)行相關(guān)性分析，獲得的相關(guān)性列表剔除在先前評估細(xì)胞分?jǐn)?shù)中使用的基因列表，在剔除后的列表中取正負(fù)相關(guān)性最高的前500個基因進(jìn)行功能聚類。對非免疫反應(yīng)通路進(jìn)行比較，觀察不同CD4+T亞群對結(jié)腸癌細(xì)胞功能調(diào)控的傾向。在聚類分析結(jié)果中，除由于免疫細(xì)胞亞群差異導(dǎo)致的免疫功能相關(guān)通路的差異，還著重關(guān)注免疫細(xì)胞亞群差異對結(jié)腸癌腫瘤細(xì)胞生物功能的影響，關(guān)注其對腫瘤細(xì)胞周期、細(xì)胞生長、細(xì)胞增殖、能量代謝、氧化磷酸化、缺氧應(yīng)激等方面的調(diào)控差異。

1.7基于免疫細(xì)胞亞群的預(yù)后模型的構(gòu)建 采用LASSO算法篩選免疫細(xì)胞亞群，并使用Cox回歸分析構(gòu)建風(fēng)險預(yù)后模型。

1.8統(tǒng)計(jì)學(xué)處理 采用SPSS 23.0軟件進(jìn)行統(tǒng)計(jì)學(xué)分析。正態(tài)分布的計(jì)量資料組間比較采用獨(dú)立樣本t檢驗(yàn)。不符合正態(tài)分布的計(jì)量資料比較采用Wilcoxon秩和檢驗(yàn)。P＜0.05為差異有統(tǒng)計(jì)學(xué)意義。

2 結(jié) 果

2.138個免疫細(xì)胞亞群的簽名矩陣 取GSE146771檢測基因數(shù)目nFeatures_RNA值在800～1200的勻質(zhì)細(xì)胞的測序數(shù)據(jù)作為后續(xù)分析的參考集，基于CIBERSORTx算法創(chuàng)建包括38個免疫細(xì)胞集群共3801個基因的特征矩陣。如圖1所示。

2.2不同TNM分期的免疫細(xì)胞亞群的比例差異

與TCGA結(jié)腸癌未轉(zhuǎn)移患者比較，轉(zhuǎn)移患者中hT04_CD4-TCF7分?jǐn)?shù)上調(diào)(P＜0.01)，hT11_CD4-CTLA4、hT12_CD8-LEF1、hB01_PlasmaB-IgG、hM10_Macro-IL1B分?jǐn)?shù)下調(diào)(P＜0.05，圖2A)，且上述hT04_CD4-TCF7、hT12_CD8-LEF1、hB01_PlasmaB-IgG、hM10_Macro-IL1B免疫細(xì)胞亞群與淋巴結(jié)轉(zhuǎn)移數(shù)存在明顯相關(guān)關(guān)系(P＜0.05，圖2B)。

2.3免疫細(xì)胞亞群與總生存期的相關(guān)性 如圖3所示，hT04_CD4-TCF7、hT03_CD4-GNLY、hT09_CD4-CXCL13、hT16_CD8-CD6、hT17_CD8-CD160、hT14_CD8-CX3CR1、hI02_NK-GZMK、hB03_FollicularB-IgD、hB01_PlasmaB-IgG、hB02_GALTB-IgA、hM01_Mast-TPSAB1、hM12_TAMC1QC、hM04_cDC1-BATF3、hM03_cDC2-CD1C、hM05_Mono-CD14與患者的總生存期明顯相關(guān)(P＜0.05)。

2.410個CD4+T細(xì)胞亞群對結(jié)腸癌細(xì)胞調(diào)控功能的異質(zhì)性 CIBERSORTx分析共獲得10個CD4+T亞群的評分。其中CD4-GNLY、CD4-ANXA1、CD4-CXCR6對結(jié)腸癌腫瘤的調(diào)控作用優(yōu)先聚集在血管生成相關(guān)通路，CD4-IL23R優(yōu)先聚集在能量代謝相關(guān)通路，CD4-GZMK優(yōu)先聚集在缺氧應(yīng)激相關(guān)通路，CD4-ANXA1、CD4-TCF7、CD4-CTLA4優(yōu)先聚集在細(xì)胞增殖相關(guān)通路，CD4-ANXA1、CD4-TCF7、CD4-CTLA4、CD4-CXCL13優(yōu)先聚集在細(xì)胞黏附及遷移相關(guān)通路(圖4A)。此外，CCR7、CXCR5兩個亞群富集結(jié)果顯示，這兩個亞群在促癌生物學(xué)過程中的作用較其他亞群不顯著。

比較轉(zhuǎn)移組與非轉(zhuǎn)移組間存在明顯差異的免疫細(xì)胞亞群，CD4+T細(xì)胞亞群中T04_CD4-TCF7與轉(zhuǎn)移呈正相關(guān)，T11_CD4-CTLA4與轉(zhuǎn)移呈負(fù)相關(guān)。使用scTPA線上工具比較兩個亞群細(xì)胞內(nèi)通路激活的差異，發(fā)現(xiàn)兩者在細(xì)胞黏附及免疫激活或抑制相關(guān)通路中的激活狀態(tài)存在明顯差異(P＜0.05，圖4B)。生存相關(guān)性分析結(jié)果顯示，T04_CD4-TCF7與患者生存呈正相關(guān)關(guān)系，T09_CD4-CXCL13與患者生存呈負(fù)相關(guān)關(guān)系(P＜0.05)。與CD4-TCF7相比，CD4-CXCL13在抗腫瘤及TCR相關(guān)通路中的激活上明顯上調(diào)(P＜0.05，圖4C)。生存分析結(jié)果顯示，T17_CD8-CD160對TCR相關(guān)通路、Fc受體及免疫刺激細(xì)胞因子相關(guān)的通路激活明顯高于T16_CD8-CD6細(xì)胞亞群(P＜0.05，圖4D)。

2.5基于免疫細(xì)胞亞群反映患者免疫狀態(tài)的預(yù)后風(fēng)險模型的構(gòu)建 通過LASSO算法從38個免疫細(xì)胞亞群中選擇4類免疫細(xì)胞亞群CD160-CD8+T細(xì)胞、CX3CR1-CD8+T細(xì)胞、BATF3-cDC1細(xì)胞、IgG-PlasmaB細(xì)胞(圖5A)。其中CD160-CD8+T細(xì)胞(P=0.0064)、IgG-PlasmaB細(xì)胞與風(fēng)險評分呈負(fù)相關(guān)，CX3CR1-CD8+T細(xì)胞(P=0.0098)、BATF3-cDC1細(xì)胞(P=0.011)與風(fēng)險評分呈正相關(guān)(P＜0.05)。采用Cox分析建立預(yù)后風(fēng)險評分模型，根據(jù)整個TCGA結(jié)腸癌隊(duì)列獲得的風(fēng)險評分，患者可分為高風(fēng)險組與低風(fēng)險組(圖5B)。Kaplan-Meier曲線表明，在訓(xùn)練集、驗(yàn)證集及整組中，高風(fēng)險組患者的死亡風(fēng)險較高(圖5C–E)。在免疫細(xì)胞亞群功能差異的分析中，CD160-CD8+T細(xì)胞比CD8+T細(xì)胞表現(xiàn)出更高的Fc受體(Fc receptors，F(xiàn)cR)激活水平及T細(xì)胞抗原受體(T cell receptor，TCR)激活水平(P＜0.05)。高風(fēng)險組患者的T細(xì)胞功能失調(diào)評分高于低風(fēng)險組，差異有統(tǒng)計(jì)學(xué)意義(P＜0.05，圖5F)。

3 討 論

免疫治療在結(jié)腸癌中的應(yīng)用較少，尚處于較初級的探索階段。鑒別結(jié)腸癌的腫瘤微環(huán)境特征，特別是免疫細(xì)胞亞群及功能狀態(tài)的特點(diǎn)，以及可能為臨床精準(zhǔn)實(shí)施免疫治療提供參考的可靠預(yù)測標(biāo)志物，對于提高結(jié)腸癌患者的免疫治療效果，提出可行的輔助治療方案具有重要意義。近年來，單細(xì)胞測序技術(shù)的發(fā)展使研究者能在更高分辨水平研究腫瘤微環(huán)境中的各種細(xì)胞成分，其中包括免疫細(xì)胞亞群在基因表達(dá)及功能狀態(tài)上的異質(zhì)性。一系列已經(jīng)發(fā)表的單細(xì)胞測序分析結(jié)果顯示，不同癌種的免疫細(xì)胞亞群譜存在差異[10-12]，癌種間免疫治療效果的差異可能與其免疫細(xì)胞譜系的差異相關(guān)，目前相關(guān)研究主要通過不同的表面標(biāo)記分子，結(jié)合免疫組化空間定位分析免疫細(xì)胞亞群的功能及預(yù)后價值[13-14]。免疫細(xì)胞調(diào)控腫瘤細(xì)胞功能及腫瘤發(fā)生發(fā)展的研究更多著眼于傳統(tǒng)的免疫細(xì)胞分群，而探討免疫微環(huán)境中免疫細(xì)胞異質(zhì)性的單細(xì)胞測序分析多側(cè)重于免疫細(xì)胞本身的基因表達(dá)特性及功能特性。為進(jìn)一步揭示結(jié)腸癌腫瘤微環(huán)境中各種免疫細(xì)胞亞群對結(jié)腸癌細(xì)胞功能調(diào)控的差異以及對腫瘤進(jìn)展的影響，本研究利用CIBERSORTx算法，基于單細(xì)胞測序分析定義的高分辨結(jié)腸癌腫瘤免疫微環(huán)境免疫細(xì)胞譜對TCGA的結(jié)腸癌患者隊(duì)列進(jìn)行研究，以期探究結(jié)腸癌微環(huán)境中特征免疫細(xì)胞譜對結(jié)腸癌腫瘤細(xì)胞的調(diào)控作用及其與腫瘤進(jìn)程的關(guān)系，結(jié)果顯示，hT04_CD4-TCF7、hT11_CD4-CTLA4、hT12_CD8-LEF1、hB01_PlasmaB-IgG、hM10_Macro-IL1B細(xì)胞亞群與結(jié)腸癌淋巴結(jié)轉(zhuǎn)移相關(guān)，hT04_CD4-TCF7、hT03_CD4-GNLY、hT09_CD4-CXCL13、hT16_CD8-CD6、hT17_CD8-CD160、hT14_CD8-CX3CR1、hI02_NK-GZMK、hB03_FollicularBIgD、hB01_PlasmaB-IgG、hB02_GALTB-IgA、hM01_Mast-TPSAB1、hM12_TAM-C1QC、hM04_cDC1-BATF3、hM03_cDC2-CD1C、hM05_Mono-CD14細(xì)胞亞群與患者的總生存期相關(guān)。

本研究還探索了不同CD4+T細(xì)胞亞群對結(jié)腸癌腫瘤細(xì)胞調(diào)控功能的差異，發(fā)現(xiàn)對于腫瘤血管生成、腫瘤細(xì)胞生長增殖、能量代謝、氧化磷酸化、缺氧應(yīng)激等通路的調(diào)節(jié)，不同CD4+T細(xì)胞亞群表現(xiàn)出不同的優(yōu)先調(diào)控方向，從而影響腫瘤的發(fā)展進(jìn)程，導(dǎo)致患者在轉(zhuǎn)移和預(yù)后結(jié)局方面的不同。

在轉(zhuǎn)移與預(yù)后分析中，本研究發(fā)現(xiàn)，CD4+和CD8+T細(xì)胞的不同亞群對轉(zhuǎn)移及預(yù)后具有相反的相關(guān)關(guān)系，通過比較這些突出亞群的細(xì)胞內(nèi)通路激活狀態(tài)，發(fā)現(xiàn)這些細(xì)胞亞群間通路激活狀態(tài)的差異很好地解釋了其對腫瘤發(fā)生發(fā)展方向的影響。先前的分析只是單獨(dú)地看待這些免疫細(xì)胞亞群的功能差異，而經(jīng)過構(gòu)建模型，可以看到這種顯著影響患者生存的腫瘤微環(huán)境差異主要來自于哪些免疫細(xì)胞亞群，從而可以從整體觀察免疫細(xì)胞亞群對患者疾病進(jìn)程的影響。腫瘤部位基于免疫細(xì)胞亞群的免疫細(xì)胞構(gòu)成的異質(zhì)性及免疫狀態(tài)，對于患者的生存具有較高的影響，該模型在訓(xùn)練集、驗(yàn)證集及整組中顯示出較好的基于免疫狀態(tài)的死亡風(fēng)險預(yù)測能力，相較單純基于差異基因簇，基于免疫狀態(tài)預(yù)測患者的死亡風(fēng)險具有更高的臨床參考價值，展現(xiàn)了患者檢測節(jié)點(diǎn)的免疫狀態(tài)。在這一基于免疫細(xì)胞亞群構(gòu)建的預(yù)后模型中，筆者發(fā)現(xiàn)，低風(fēng)險組與高風(fēng)險組中有4個免疫細(xì)胞亞群的構(gòu)成存在差異，這可能表明某些免疫細(xì)胞亞群之間存在免疫應(yīng)答過程中的相互作用，降低腫瘤微環(huán)境中BATF3-cDC1樹突細(xì)胞的浸潤可降低CX3CR1-CD8+T細(xì)胞的浸潤，或提升IgG-PlasmaB細(xì)胞的浸潤，進(jìn)而提升CD160-CD8+T細(xì)胞的浸潤，改善CD8+T細(xì)胞群的功能活性，從而改善患者生存，或輔助免疫檢測點(diǎn)治療。

綜上所述，結(jié)腸癌腫瘤微環(huán)境中的免疫細(xì)胞構(gòu)成及免疫微環(huán)境狀態(tài)能夠通過對結(jié)腸癌細(xì)胞功能的差異調(diào)控影響結(jié)腸癌的進(jìn)展，進(jìn)而影響患者的預(yù)后。本研究基于單細(xì)胞測序數(shù)據(jù)，使用CIBERSORTx算法在TCGA隊(duì)列中研究免疫細(xì)胞譜對結(jié)腸癌進(jìn)展及預(yù)后的影響，為進(jìn)一步明確各免疫細(xì)胞亞群對腫瘤微環(huán)境的影響提供了方向和線索。但本研究尚有一定的局限，如研究結(jié)果基于公共數(shù)據(jù)和算法分析，未能闡述各免疫細(xì)胞亞群之間的相互關(guān)系網(wǎng)絡(luò)，且需要進(jìn)一步的細(xì)胞生物學(xué)實(shí)驗(yàn)對結(jié)果進(jìn)行驗(yàn)證。