馬強(qiáng)，張霞*，吳敬醫(yī)，丁彩云，谷冰雨，徐四娟

（1.皖南醫(yī)學(xué)院弋磯山醫(yī)院超聲醫(yī)學(xué)科，安徽 蕪湖 241001；2.重癥醫(yī)學(xué)科；3.婦產(chǎn)科）

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

1 材料與方法

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

1.2 鑒定差異基因分別對基質(zhì)及免疫細(xì)胞高、低分組，使用limma R軟件包根據(jù)以下臨界值識(shí)別出基質(zhì)細(xì)胞和免疫細(xì)胞中差異表達(dá)基因：錯(cuò)誤發(fā)現(xiàn)率（FDR）＜0.05和|log2倍數(shù)變化（FC）|＞2。再分別將基質(zhì)細(xì)胞打分中上調(diào)、下調(diào)基因與免疫細(xì)胞打分中上調(diào)、下調(diào)基因取交集，得出腫瘤微環(huán)境相關(guān)的差異基因（DEGs）。

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

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

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

1.6 GSEA富集分析借助GSEA 4.0.3軟件對核心基因進(jìn)行GSEA富集分析。

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

2 結(jié)果

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

圖1 免疫、基質(zhì)總評分與臨床參數(shù)的相關(guān)性分析散點(diǎn)圖

2.2 鑒定DEGs為了更好地了解基因表達(dá)與免疫、基質(zhì)評分之間的相關(guān)性，我們分析了510個(gè)PTC樣本的基因表達(dá)譜。對于免疫評分，高分組與低分組相比，上調(diào)了690個(gè)基因，下調(diào)了60個(gè)基因（|log2倍數(shù)變化（FC）|＞2，F(xiàn)DR＜0.05）。對于基質(zhì)評分，高分組與低分組相比，495個(gè)基因上調(diào)，8個(gè)基因下調(diào)（|log2倍數(shù)變化（FC）|＞2，F(xiàn)DR＜0.05）。分析免疫和基質(zhì)評分的高分組中共同的上調(diào)和下調(diào)基因，發(fā)現(xiàn)了434個(gè)上調(diào)和7個(gè)下調(diào)的基因。獲得的這441個(gè)基因被視為DEGs，并應(yīng)用于進(jìn)一步分析。

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

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

2.4 DEGs的PPI網(wǎng)絡(luò)的構(gòu)建和核心基因的鑒定為了進(jìn)一步探討DEGs及其調(diào)控機(jī)制，我們利用STRING在線數(shù)據(jù)庫分析了DEGs，并構(gòu)建了PPI網(wǎng)絡(luò)。在PPI網(wǎng)絡(luò)中，我們確定了139個(gè)基因（圖3A）和7個(gè)功能模塊（使用MCODE功能）。為更好地鑒定PTC微環(huán)境的核心基因，我們使用了Cytoscape軟件中的cytohubba度數(shù)算法進(jìn)一步分析了DEGs的PPI網(wǎng)絡(luò)。根據(jù)它們的程度分?jǐn)?shù)篩選了16個(gè)核心基因。排在前3位的核心基因如下：SAA1、CXCL10、CXCL11（圖3B）。

圖3 DEGs的PPI網(wǎng)絡(luò)構(gòu)建（A）和核心基因的鑒定（B）

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

圖4 SAA1臨床相關(guān)性的箱線圖與生存曲線分析

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

圖5 GSEA富集分析圖

2.7 免疫細(xì)胞浸潤的相關(guān)性分析通過上述分析結(jié)果可得出腫瘤微環(huán)境的核心基因SAA1作用于免疫機(jī)制影響甲狀腺癌的臨床狀態(tài)。

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

圖6 免疫細(xì)胞相關(guān)性分析圖以及SAA1與免疫細(xì)胞的小提琴圖

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

3 討論

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

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

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

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