單細胞轉錄組測序分析結腸癌微環(huán)境中免疫細胞亞群與癌癥進程的關系

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

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

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

1 資料與方法

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

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

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

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

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

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

1.7基于免疫細胞亞群的預后模型的構建 采用LASSO算法篩選免疫細胞亞群，并使用Cox回歸分析構建風險預后模型。

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

2 結 果

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

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

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

2.3免疫細胞亞群與總生存期的相關性 如圖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與患者的總生存期明顯相關(P＜0.05)。

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

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

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

3 討 論

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

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

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

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