基于多平臺數據識別胰腺導管腺癌預后相關腫瘤微環境基因

蒲垠全 馬雨凡 彭莉 湯小偉 彭燕

西南醫科大學附屬醫院消化內科，瀘州 646000

PDAC是最常見的惡性腫瘤之一，其預后不良，目前居癌癥死亡原因第7位，預計將在未來上升到第3位[1]。目前研究認為，胰腺周圍復雜的區域解剖結構、高侵襲性和獨特的腫瘤微環境是導致其預后不良的主要原因。盡管在過去的幾十年里，PDAC在外科手術、放療和化療方面都取得了長足的進展，但其5年總生存率仍不到5%。免疫治療的臨床試驗也未顯示出對大多數PDAC患者的臨床益處[2]，同時也缺乏可靠的生物學標志物來選擇有效的分子亞型。腫瘤微環境通過對T細胞浸潤形成物理和化學屏障、阻止藥物釋放等方式降低治療效果[3]，PDAC獨特的腫瘤微環境被認為是免疫治療抵抗的主要原因之一[4]。本研究旨在探索PDAC相關腫瘤微環境，以確定影響患者預后的生物學標志物和潛在的免疫治療靶點。

材料與方法

一、數據集的選擇

從腫瘤基因組圖譜(The Cancer Genome Atlas, TCGA)數據庫下載142例PDAC患者的基因表達譜和臨床隨訪數據，從基因表達綜合(Gene Expression Omnibus，GEO)數據庫的兩個微陣列數據集(GSE21501和GSE62452)下載168例具有完整臨床隨訪信息的PDAC患者的mRNA數據。

二、細胞類型富集分析及比較

xCell是一個用來對不同平臺的轉錄組數據進行腫瘤微環境數字剖析的工具[5]。它整合了反卷積基因富集分析的優點，能根據基因表達譜計算64種細胞類型(包括免疫細胞和基質細胞)富集分數，其富集分數范圍為-1～1。本研究運用xCell網絡工具對基因表達數據進行細胞類型富集分析。TCGA數據庫中的數據分析結果可以從xCell門戶網站上直接下載。對于GEO數據則需要使用R/Bioconductor軟件中的sva包來移除數據集中的批間差，然后再將規范化后的數據傳輸到xCell網站進行分析。以細胞富集評分中位數將TCGA的142例患者分為高分、低分2組，使用R/Bioconductor軟件中的survival包進行單因素生存分析確定有預后價值的細胞類型，再使用GEO收集的數據集進行驗證，P<0.05為差異有統計學意義。

三、基因表達差異分析及比較

四、差異表達基因的功能富集分析

為了進一步了解與預后相關DEGs的功能，使用R/Bioconductor軟件中的clusterprofiler包對預后相關DEGs進行功能富集分析，包括基因本體(gene ontology, GO)分析和京都基因與基因組百科全書(Kyoto Encyclopedia of Genes and Genomes, KEGG) 通路富集分析，其中GO包括分子功能(molecular function, MF)、生物學過程(biological process, BP)和細胞成分(cellular component, CC)。FDR<0.05為差異有統計學意義。

五、蛋白質-蛋白質交互作用網絡的建立

運用STRING(Search Tool for the Retrieval of Interacting proteins)數據庫建立預后相關DEGs的蛋白質-蛋白質相互作用(protein-protein interaction, PPI)網絡，然后通過Cytoscape軟件進行PPI網絡的重建，并使用Cytoscape MCODE插件查找集群，以定位密集連接的區域，同時計算網絡各節點的連接度。

六、TISIDB數據庫分析共同預后相關DEGs

TISIDB數據庫是一個用于腫瘤和免疫系統相互作用的集成存儲門戶網站。本研究先確定TCGA數據集及GEO數據集兩個平臺數據的共同預后相關DEGs，再使用TISIDB數據庫來評估共同預后相關DEGs在腫瘤與免疫系統的交互作用。

結 果

一、細胞類型富集評分及其與PDAC患者生存率的相關性

通過xCell網絡工具對細胞類型富集評分，并根據富集評分中位數將PDAC患者分成高分、低分2組，結果顯示，在64種細胞類型中只有Th1細胞和角質形成細胞在TCGA和GEO數據集中具有相同的預后價值(圖1)。Th1細胞富集評分中位數在TCGA數據集和GEO數據集中分別為0.04155(0.00000～0.22100)、0.2183(0.00000～0.64560)，角質形成細胞富集評分中位數在TCGA數據集和GEO數據集中分別為0.12850(0.01030～0.39400)、0.2700(0.00000～0.50350)，Th1細胞和角質形成細胞高分組患者的總生存率中位數均顯著低于低分組，差異均有統計學意義(P值均<0.05)。

二、與PDAC患者預后相關的DEGs篩選

TCGA數據庫的142例患者Th1細胞富集高分組與低分組比較，高分組有379個基因上調，294個基因下調(圖2A)；角質形成細胞富集高分組與低分組比較，高分組有481個基因上調，298個基因下調(圖2B)。 通過單變量生存分析并做韋恩圖(圖2C)，共篩選出216個與PDAC患者預后相關的DEGs。

三、216個DEGs的功能富集分析

對216個DEGs進行了GO分析和KEGG通路富集分析，結果顯示，以-log(矯正P值)>3設為閾值，21個BP、3個CC和5個MF被鑒定為差異有統計學意義(P值均<0.05)。21個BP中有9個(淋巴細胞分化、抗原受體介導的信號通路調控、淋巴細胞活化的調控、免疫系統過程的負調節、趨化因子介導的信號通路、B細胞活化、T細胞活化、對趨化因子應答、細胞對趨化因子的應答)與免疫過程密切相關(圖3A)。5條KEGG通路中4條通路(病毒蛋白與細胞因子和細胞因子受體的相互作用、細胞因子-細胞因子受體相互作用、趨化因子信號通路、原發性免疫缺陷)與免疫過程密切相關(圖3B)。

圖1 從GEO數據集篩選的Th1細胞(1A)和角質形成細胞(1B )富集評分以及從TCGA數據集篩選的Th1細胞(1C)和角質形成細胞(1D)富集評分的患者生存曲線

圖2 Th1細胞高低富集評分組(2A)和角質形成細胞高低富集評分組(2B)的DEGs火山圖以及預后DEGs的韋恩圖(2C)

四、216個DEGs的PPI網絡和模塊化分析

為了更好地理解216個DEGs之間的相互作用關系，使用STRING網絡工具建立了一個包含158個節點和431個邊的PPI網絡(圖4A)，并使用Cytoscape MCODE插件鑒定了一個MCODE得分>5分的模塊(圖4B)，結果顯示，CCR7、CD27、CD5、CXCL13、ZAP70、MS4A1和CCL19基因的連通度顯著高于該模塊中的其他基因。

五、216個DEGs在GEO數據集中的驗證

Kaplan-Meier單變量分析結果顯示，15個基因在GEO和TCGA數據集中具有相似的預后價值，其中6個基因(CNTNAP4、DLGAP1、GIMAP7、LGI1、SFRP1、TPH1)的高表達與生存期延長有關，9個基因(FAM83A、HMGA2、ITGA3、ITGB6、KRT7、KRT13、PTGES、S100A2、UBE2C)高表達與生存期縮短有關(圖5)。

圖3 參與免疫相關的9個生物過程(3A)以及參與5條KEGG通路(3B)的基因

圖4 216個與預后相關DEGs的PPI網絡(4A)及MCODE為11.333分的模塊圖(4B)。從藍色到紅色的梯度代表基因從下調到上調的變化過程，節點大小代表指定蛋白的連通度

圖5 TCGA(5A)和GEO(5B)數據集具有相似預后價值的15個基因的Kaplan-Meier生存曲線圖

七、共同差異表達基因與免疫相關因素的相關性

在TISIDB數據庫檢索上述15個共同基因，結果發現GIMAP7與胰腺癌的免疫相關因素淋巴細胞、趨化因子受體、趨化因子、免疫抑制因子、免疫刺激因子、主要組織相容性復合體(major histocompatibility complex, MHC)均呈顯著正相關(圖6)。

討 論

本研究通過TCGA數據庫中142例PDAC患者和GEO的168例PDAC患者資料篩選出與預后相關的Th1細胞及角質形成細胞，依據富集評分高低分組比較來確定與PDAC預后相關的腫瘤微環境基因。通過DEGs分析和進一步的生存分析明確與PDAC預后有關的216個 DEGs。對這些DEGs的功能分析結果顯示在生物學過程GO條目和KEGG通路中有近一半的GO條目和大多數KEGG通路與免疫反應密切相關。216個DEGs的PPI網絡顯示CCR7、CD27、CD5、CXCL13、ZAP70、MS4A1和CCL19具有更高的連通度，被鑒定為與PDAC密切相關的中樞基因。這7個基因參與了上述9個免疫相關的生物過程和4條免疫相關的KEGG通路，表明這些基因不僅與PDAC的預后密切相關，而且廣泛參與了免疫調節。

上述7個基因在腫瘤分化、轉移和預后中的作用均已有相關研究報道[5-11]。CCR7、CCL19、CD27和CXCL13是參與功能富集最多的4個基因，最有可能成為治療靶點。趨化因子CCL19及其受體CCR7可趨化樹突狀細胞進入腫瘤組織[12]，誘導T細胞殺傷腫瘤細胞[13]及免疫細胞浸潤腫瘤組織并釋放細胞因子，還能通過ERK和PI3K/AKT信號轉導上調Twist基因的表達，從而調控上皮-間充質轉化(epithelial-mesenchymal transition, EMT)過程，抑制腫瘤的增殖、侵襲和轉移。因此，CCR7/CCL19被廣泛推薦為癌癥治療的前瞻性和實用性靶點[6]。CD27基因是腫瘤壞死因子超家族的一員，其與CD28協同作用，促進T細胞的初始激活和擴增，還能維持外周效應細胞的擴張和存活，并增強其細胞溶解和細胞因子功能[9]。近年來有報道表明，抗人CD27單克隆抗體與PD-1阻滯劑聯合應用可提高腫瘤疫苗和小鼠免疫治療的效果[14-15]。CXCL13基因通過其受體CXCR5發揮生物學功能，CXCL13/CXCR5軸通過整合多個信號事件包括PI3K/AKT途徑的激活、基質金屬蛋白酶表達的增加等來促進癌細胞的增殖、存活和遷移。CXCL13/CXCR5軸還可通過下調T效應細胞直接幫助腫瘤細胞逃避宿主免疫監視，及誘導IL-10途徑間接幫助腫瘤細胞逃避T效應細胞免疫[16]。最近，針對CXCL13或CXCR5的敲除或抑制基因療法、抑制性抗體和載藥納米粒在其他癌癥中取得了一些進展[16-18]，然而針對上述基因的PDAC的治療還鮮有報道。

本研究從GEO數據庫下載的168例PDAC患者的基因表達譜數據對上述216個DEGs 進行交叉驗證，最終確定了15個與預后相關的腫瘤微環境相關基因，包括CNTNAP4、DLGAP1、FAM83A、GIMAP7、HMGA2、ITGA3、ITGB6、KRT7、KRT13、LGI1、PTGES、S100A2、SFRP1、TPH1和UBE2C。這些基因已被多項研究證明與腫瘤的發生發展有關，或被用作預測不同腫瘤預后的生物學標志物。

EMT是一個由多種腫瘤微環境因子觸發的生理過程，導致癌細胞從上皮到間充質表型發生形態學和遺傳學改變，它是胰腺癌細胞高轉移潛能和耐藥性的基礎[19]。先前研究表明，HMGA2、ITGA3、KRT13、S100A2和UBE2C可在許多癌癥中誘發EMT，而SFRP1可抑制EMT[20-23]。本研究結果顯示，SFRP1高表達與較長的生存期相關，而HMGA2、ITGA3、KRT13、S100A2和UBE2C基因高表達與較短的生存期相關，該結果與它們在EMT中的作用機制是一致的。

本研究運用TISIDB數據庫分析上述15個共同DEGs與免疫系統的關系，顯示GIMAP7與多種免疫相關因子包括淋巴細胞、趨化因子、趨化因子受體、免疫抑制因子、免疫刺激因子和MHC分子存在更高的相關性。進一步觀察發現，GIMAP7還與鑒定的中樞基因包括CCR7、CCL19、CXCL13和CD27也密切相關。GIMAP7蛋白是一種GTP酶，屬于免疫相關蛋白(group of immunity-associated proteins，GIMAP)家族，定位于內質網和高爾基體。GIMAP家族蛋白在調節T細胞的發育、選擇、存活、自身免疫和內環境穩定方面起著關鍵作用[24-25]。Schwefel等[25]報道，GIMAP7可刺激GIMAP2的GTP水解。Dion等[26]報道，在不同胸腺亞群中，雙陰性(double negative，DN)和雙陽性(double positive，DP)個體之間的大鼠GIMAP7表達顯著增加，并且在從CD4+CD8+DP胸腺細胞過渡到CD4+或CD8+單陽性胸腺細胞期間增加。但GIMAP7參與免疫調節的機制尚不清楚，需要進一步研究。

總之，本研究鑒定了一組與PDAC預后相關的腫瘤微環境基因及其7個中樞基因，并提供了CCR7、CCL19、CD27、CXCL13和GIMAP7 5個新的PDAC免疫治療潛在靶點。

利益沖突所有作者均聲明不存在利益沖突