生物信息學鑒定甲狀腺乳頭狀癌腫瘤微環境中免疫相關基因

馬強，張霞*，吳敬醫，丁彩云，谷冰雨，徐四娟

（1.皖南醫學院弋磯山醫院超聲醫學科，安徽 蕪湖 241001；2.重癥醫學科；3.婦產科）

甲狀腺癌作為最常見的內分泌惡性腫瘤，患病率逐年增加，其中甲狀腺乳頭狀癌（PTC）是最常見的亞型［1］。PTC的進展和整個發生發展過程涉及到基因改變，引起細胞增殖、細胞生長、血管生成和分化損失［2-3］。其中腫瘤微環境在PTC的侵襲和轉移過程中發揮關鍵作用［1］。為探究PTC腫瘤微環境相關基因的作用機制，我們研究了癌癥基因組圖譜（TCGA）項目成果，進行了全面和多平臺的方法來分析496個PTC的多維分子數據（不包括差和未分化的癌細胞，匯編所有腫瘤基因的改變），并將其與臨床相關參數聯系起來［4-5］。本文通過生物信息學方法對TCGA中PTC進行數據挖掘及分析，找出與PTC腫瘤微環境免疫相關的關鍵基因，為臨床上從基因上預測甲狀腺癌的臨床發展狀態，以及為治療各期甲狀腺癌提供新的靶點。

1 材料與方法

1.1 數據庫PTC患者的基因表達譜及臨床資料（如性別、年齡、分期、等級、生存率和最終結局）可從TCGA數據庫（https://portal.gdc.cancer.gov/）獲得。PTC患者的免疫評分和基質評分通過使用R軟件（R Foundation for Statistical Computing，Vienna，Austria）中的ESTIMATE包來計算。根據免疫、基質評分的中位值，將所有患者分別分為高、低免疫細胞評分組以及高、低基質細胞評分組。

1.2 鑒定差異基因分別對基質及免疫細胞高、低分組，使用limma R軟件包根據以下臨界值識別出基質細胞和免疫細胞中差異表達基因：錯誤發現率（FDR）＜0.05和|log2倍數變化（FC）|＞2。再分別將基質細胞打分中上調、下調基因與免疫細胞打分中上調、下調基因取交集，得出腫瘤微環境相關的差異基因（DEGs）。

1.3 DEGs的GO和KEGG富集分析為了解DEGs在THCA中作用機制以及潛在生物學功能，在R軟件中使用clusterProfiler軟件包［6］進行GO分析［7］和KEGG富集分析［8］，包括生物學過程（BP）、分子功能（MF）和細胞成分（CC），P＜0.05為差異有統計學意義，最后將結果用ggplot2軟件包繪制氣泡圖和圈圖。

1.4 PPI網絡的建立與模塊分析為了了解不同基因之間的相互作用，使用STRING（https://stringdb.org/）在線數據庫［9］分析DEGs的蛋白質-蛋白質相互作用（PPI）網絡，并考慮了臨界標準（最低要求互動得分＞0.90）。使用3.8.0版（https://cytoscape.org）的cytoscape重建該網絡，并使用分子復雜檢測（MCODE）插件以識別PPI網絡中的所有模塊。

1.5 核心基因的鑒定與臨床相關性分析使用cytoHubba（一個Cytoscape插件）中的度數算法［10］選擇了排名前3位的核心基因，并與單因素Cox回歸分析結果取交集得出差異的核心基因，并從中選擇部分基因行單基因分析，通過R軟件分析該基因在正常樣本與甲狀腺癌腫瘤樣本間的表達差異，以及與甲狀腺癌的預后以及臨床分期的相關性。

1.6 GSEA富集分析借助GSEA 4.0.3軟件對核心基因進行GSEA富集分析。

1.7 免疫細胞浸潤的差異分析及相關性分析先對下載的基因表達數據進行歸一化和過濾。然后利用CIBERSORT算法推斷甲狀腺癌和正常甲狀腺中22種侵襲性免疫細胞的豐度，包括B細胞、樹突狀細胞、嗜酸性粒細胞、巨噬細胞、肥大細胞、單核細胞、中性粒細胞、自然殺傷（NK）細胞、漿細胞、T細胞及其不同亞型。通過使用蒙特卡羅采樣，CIBERSORT得到了每個樣本的P值，這為結果提供了一個置信度度量。以核心基因分組分別得出各免疫細胞與核心基因的相關性，繪制小提琴圖以及散點圖。

2 結果

2.1 免疫評分和基質評分與臨床參數顯著相關我們從TCGA數據庫下載了496例PTC患者的基因表達譜和臨床信息。通過R中的ESTIMATE包對這些患者的免疫和基質評分進行匹配。共獲得了510例具有完整免疫和基質評分的樣本；基質評分從-1 709.75至1 617.75，免疫評分從-1 298.62至3 244.00。使用中位免疫、基質評分作為臨界值，將所有PTC樣本分為高、低免疫評分組和高、低基質評分組。使用R中的生存軟件包對患者免疫或基質評分與總體生存率之間的關系進行分析，結果顯示，高、低免疫和基質評分組患者預后差異無統計學意義。此外，樣本免疫和基質總評分與臨床參數之間的關系分析發現，兩者的總評分在不同的臨床參數（臨床分級、TNM分期）中存在顯著的相關性（圖1）。

圖1 免疫、基質總評分與臨床參數的相關性分析散點圖

2.2 鑒定DEGs為了更好地了解基因表達與免疫、基質評分之間的相關性，我們分析了510個PTC樣本的基因表達譜。對于免疫評分，高分組與低分組相比，上調了690個基因，下調了60個基因（|log2倍數變化（FC）|＞2，FDR＜0.05）。對于基質評分，高分組與低分組相比，495個基因上調，8個基因下調（|log2倍數變化（FC）|＞2，FDR＜0.05）。分析免疫和基質評分的高分組中共同的上調和下調基因，發現了434個上調和7個下調的基因。獲得的這441個基因被視為DEGs，并應用于進一步分析。

2.3 DEGs的GO功能和KEGG途徑富集分析為了解這441個DEGs與PTC相關性的潛在機制，使用R軟件對DEG進行了GO功能注釋和KEGG途徑富集分析。DEGs的GO功能分析分為生物過程（BP）、分子功能（MF）和細胞成分（CC），結果分析表明，主要的BP富集通路是體液免疫反應。在CC方面，免疫球蛋白復合物通路富集顯著。對于MF，最顯著的為抗原結合通路（圖2A）。KEGG途徑富集分析的結果表明，DEGs主要富集于包括細胞因子-細胞因子受體相互作用、造血細胞譜系、病毒蛋白與細胞因子及其受體的相互作用、趨化因子信號通路（圖2B）。

圖2 DEGs的GO功能和KEGG途徑富集通路圖

2.4 DEGs的PPI網絡的構建和核心基因的鑒定為了進一步探討DEGs及其調控機制，我們利用STRING在線數據庫分析了DEGs，并構建了PPI網絡。在PPI網絡中，我們確定了139個基因（圖3A）和7個功能模塊（使用MCODE功能）。為更好地鑒定PTC微環境的核心基因，我們使用了Cytoscape軟件中的cytohubba度數算法進一步分析了DEGs的PPI網絡。根據它們的程度分數篩選了16個核心基因。排在前3位的核心基因如下：SAA1、CXCL10、CXCL11（圖3B）。

圖3 DEGs的PPI網絡構建（A）和核心基因的鑒定（B）

2.5SAA1單基因分析在下載的樣本數據中對篩選出的3個核心基因分別進行表達差異分析，發現唯有SAA1存在表達差異，P=0.002 7（圖4A）。對SAA1進行臨床相關性的分析，包括性別、年齡、TNM分期以及分級等臨床信息，發現SAA1在PTC不同的臨床狀態如腫瘤的大小、淋巴結轉移以及腫瘤的分級存在較顯著的差異，當SAA1高表達時，PTC患者更易發生淋巴結轉移以及處于更高的分期。然而生存分析顯示SAA1與PTC的預后無明顯相關性（圖4），此結果可能與PTC患者通常預后良好，生存率較高有關。

圖4 SAA1臨床相關性的箱線圖與生存曲線分析

2.6 GSEA富集分析按SAA1表達水平的中位數將樣本分為SAA1高表達組和SAA1低表達組，分別導入GSEA 4.0.3軟件進行GSEA富集分析，結果表明，SAA1高表達樣本中有22個富集通路（P＜0.05），選擇其中與免疫相關密切的抗原加工提呈通路、細胞黏附分子通路、趨化因子信號通路、細胞因子相互作用通路等9個通路繪制GSEA富集分析圖（圖5）。結果表明，免疫相關通路與SAA1高表達的腫瘤樣本呈正相關。

圖5 GSEA富集分析圖

2.7 免疫細胞浸潤的相關性分析通過上述分析結果可得出腫瘤微環境的核心基因SAA1作用于免疫機制影響甲狀腺癌的臨床狀態。

為了進一步研究通過哪些免疫細胞發揮效能，我們對TCGA數據進行歸一化和過濾。然后利用CIBERSORT算法推斷甲狀腺癌和正常甲狀腺中22種侵襲性免疫細胞的豐度（圖6A）以及免疫細胞之間的相關性（圖6B）。再根據SAA1的表達計算各個免疫細胞在SAA1高表達組和低表達組中的含量，結果發現：CD4記憶T細胞、調節性T細胞、γ-δ T細胞、靜息的樹突狀細胞、活化的樹突狀細胞、靜息的肥大細胞差異存在統計學意義（圖6C）。

圖6 免疫細胞相關性分析圖以及SAA1與免疫細胞的小提琴圖

2.8 作用免疫細胞將SAA1基因高、低表達分組后，免疫細胞含量存在差異；通過spearman方法檢驗SAA1與各個免疫細胞的相關性，再將兩者得出的免疫細胞取交集可得出，5個SAA1作用的免疫細胞。分別為CD4記憶T細胞、調節性T細胞、γ-δ T細胞、靜息的樹突狀細胞、活化的樹突狀細胞。

3 討論

甲狀腺癌是最常見的內分泌腺癌，PTC是最常見的甲狀腺癌。由于目前對各種形式的甲狀腺癌遺傳背景的了解還遠遠不夠，因此大多數研究都集中在遺傳基礎上［11］。目前已有PTC生物信息學分析，結果顯示在GEO平臺中數據挖掘分析發現LPAR5、TFPI和ENTPD1與PTC患者的發展和預后相關［12］，這些基因指標常與免疫浸潤相關。炎癥目前被認為是惡性腫瘤的重要組成部分［13］。如今腫瘤微環境已吸引了腫瘤研究界的關注，并在腫瘤發生和發展、治療和預后中已被確認為一個關鍵因素［14-16］。越來越多的證據表明，腫瘤微環境的各種成分，如免疫細胞、可溶性因子和改變的細胞外基質，都對癌癥的進展有積極作用［17］。而與腫瘤微環境相關的基因與甲狀腺癌之間的聯系尚未得到充分闡明。

本文采用生物信息學方法研究TCGA數據庫中PTC基因表達譜及臨床信息與腫瘤微環境之間的聯系。首先使用ESTIMATE算法計算出樣品中免疫、基質細胞的評分以及對臨床狀態如TNM分期等進行相關性分析，發現腫瘤微環境影響了PTC的發病，并作用于PTC的臨床分期。本研究使用limma R軟件包識別基質細胞高、低分組與免疫細胞高、低分組中DEGs，并對免疫、基質細胞高評分組和低評分組取交集得出441個DGEs。GO和KEGG分析得出主要的細胞因子受體相互作用等富集通路。進一步分析了DEGs的PPI網絡，得出了SAA1等16個基因為重要連接站點。并對連接數前三的基因SAA1、CXCL10、CXCL11進行單基因分析發現，唯有SAA1在TCGA正常樣品與腫瘤樣品中的表達差異存在統計學意義，在配對的腫瘤樣品和癌旁組織中，SAA1在腫瘤樣品中表達量顯著增高（P=0.002 7）。為進一步研究SAA1作用機制，我們對其行GSEA富集分析，發現結果里有22條通路上調，其中9條通路與免疫密切相關。進一步分析發現，SAA1可能通過影響免疫通路以及CD4記憶T細胞、調節性T細胞、γ-δ T細胞、靜息的樹突狀細胞、活化的樹突狀細胞，最終影響甲狀腺癌患者的臨床狀態。

目前研究表明SAA1在胰腺癌的發生發展中具有重要的調節作用［18］。高表達的SAA1具有作為晚期腎透明細胞癌患者的診斷和預后生物標志物的潛力，其可能作為晚期腎透明細胞癌患者的新型治療靶點［19］。另外有研究表明SAA1可介導NOX4/ROS通路并激活p38MAPK/NF-κB通路，促進炎癥介質的釋放［20］。SAA1與機體的炎癥反應存在一定關聯，這與本研究的結果相符。本研究發現SAA1與PTCTNM分期及臨床分級具有相關性，而與年齡、性別以及預后相關性無統計學意義。PTC腫瘤通常生長緩慢，且結合手術、甲狀腺激素和放射性碘治療效果較好，故PTC的預后一般都較好［21］，且在TCGA數據庫中的PTC數據是排除惡性度較高的癌癥樣本，所以可能會導致最終結果顯示SAA1對于PTC的預后并不存在預測價值。但研究表明一旦PTC惡化，更具侵襲性和致命性，故若想預測PTC的預后，仍需對PTC發展進行更深入的機制研究，也可聯合對比增強造影超聲來預測甲狀腺癌預后［22］。另外關于SAA1如何通過免疫反應通路影響PTC的發生發展仍需要更進一步的實驗室研究。

本文的局限性在于只納入了TCGA數據庫中PTC數據，沒有與其他數據庫進行聯合比較。本研究結果表明腫瘤免疫微環鏡對于PTC的發生發展具有重要意義，并發現了SAA1作為PTC潛在的免疫浸潤的生物標志物，在未來將有助于PTC的診斷、預后和靶向治療的進一步臨床應用。