凡志祥 李坤 張彥彥 左世凱 黃冬梅

子宮頸癌是最常見的婦科惡性腫瘤，其發(fā)病率和死亡率在女性惡性腫瘤中位居第四。子宮頸鱗狀細(xì)胞癌（Cervical squamous cell carcinoma，CSCC）占子宮頸癌的80%～85%，是最常見的子宮頸癌亞型[1]。近年來，腫瘤免疫檢查點治療成為子宮頸癌治療的研究熱點。然而，尚未發(fā)現(xiàn)用于子宮頸癌免疫治療的特異性生物標(biāo)記物[2]。因此，探索免疫相關(guān)生物標(biāo)記物成為研究子宮頸癌免疫治療的重要方向。檢查點免疫療法已徹底改變了黑色素瘤的治療方式，并在肺癌、腎癌和頭頸癌中顯示出顯著的療效[3～6]。目前，子宮頸癌檢查點免疫治療的試驗尚處于早期階段，而將其作為一種治療選擇主要是因為子宮頸癌是一種病毒驅(qū)動的癌癥，免疫系統(tǒng)應(yīng)將其視為異物[7]?，F(xiàn)已證實CD8+T 細(xì)胞有助于腫瘤適應(yīng)性免疫，且在大多數(shù)實體瘤中，高度浸潤的CD8+T 細(xì)胞對腫瘤治療有益[8，9]。且與正常子宮頸組織相比，子宮頸癌中的CD8+T 細(xì)胞更為豐富[10，11]。因此，鑒定與CD8+T 細(xì)胞浸潤相關(guān)的新生物標(biāo)志物將有助于探索免疫浸潤機(jī)制和豐富免疫治療。

加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析（Weighted gene coexpression network analysis，WGCNA）是一種分析多樣本基因表達(dá)譜的方法。它可以通過對具有相似表達(dá)模式的基因進(jìn)行聚類，形成模塊，并分析模塊與特定特征（如患者的臨床信息）之間的關(guān)系，該算法已廣泛用于尋找生物標(biāo)記物[12]。通過估計轉(zhuǎn)錄本的相對子集進(jìn)行細(xì)胞類型識別（CIBERSORT）是另一種用于分析基因表達(dá)數(shù)據(jù)的生物信息學(xué)工具，該工具使用反卷積算法量化免疫細(xì)胞的組成，已成功用于估計各種癌癥中的免疫細(xì)胞浸潤水平，如乳腺癌和肝癌[13，14]。

為探索腫瘤微環(huán)境的特征并確定CSCC 的潛在生物標(biāo)志物，本研究使用CSCC 的基因表達(dá)數(shù)據(jù)進(jìn)行WGCNA，并利用CIBERSORT 算法計算樣品的T 細(xì)胞組成，然后鑒定出與CD8+T 細(xì)胞浸潤水平相關(guān)的關(guān)鍵模塊和中樞基因。

1 材料與方法

1.1 材料 通過TCGA 數(shù)據(jù)庫（https：//portal.gdc.cancer.gov）獲取304 個原發(fā)CSCC 樣本的基因表達(dá)譜，并選擇中位絕對偏差前10%的5 605 個基因進(jìn)行額外分析。從UCSC 數(shù)據(jù)庫（https：//xenabrowser.net/）中檢索285 例CSCC 患者的臨床數(shù)據(jù)，用于后續(xù)分析。

1.2 方法

1.2.1 腫瘤浸潤性免疫細(xì)胞（Tumor infiltrating immune cells，TIIC）的評估 使用R 包“CIBERSORT”估算TCGA-CSCC 數(shù)據(jù)集中每個樣本的22 種TIIC分?jǐn)?shù)[15]。并選擇每個樣本中T 細(xì)胞的7 個亞型分?jǐn)?shù)作為WGCNA 的性狀數(shù)據(jù)。

1.2.2 加權(quán)基因共表達(dá)網(wǎng)絡(luò)的構(gòu)建 為了解5 605個基因之間的相互關(guān)系，并確定基因模塊和影響腫瘤免疫細(xì)胞浸潤的中樞基因，使用R Studio 中的“WGCNA”軟件包構(gòu)建加權(quán)基因共表達(dá)網(wǎng)絡(luò)[16]。首先，基于成對基因之間的Pearson 相關(guān)系數(shù)，將單個基因的表達(dá)水平轉(zhuǎn)換為相似矩陣。使用表達(dá)式矩陣構(gòu)造一個分層聚類樹來檢測異常值。然后根據(jù)網(wǎng)絡(luò)拓?fù)浞治龊瘮?shù)計算出無標(biāo)度拓?fù)鋽M合指數(shù)作為軟閾值功率的函數(shù)。接著使用拓?fù)渲丿B不相似性度量（TOM）計算具有相似表達(dá)譜基因之間的模內(nèi)連接性。最后采用一種動態(tài)混合切割方法將這些基因分成不同的模塊，每個模塊至少有30 個基因。

1.2.3 構(gòu)建模塊-特征關(guān)系 應(yīng)用Pearson 檢驗法計算模塊特征基因與T 細(xì)胞浸潤水平之間的相關(guān)性，選擇P值顯著并且相關(guān)系數(shù)高的T 細(xì)胞亞型和模塊，定義其為中心模塊。為確定中心模塊中基因的功能，使用R 包“clusterProfiler”進(jìn)行富集分析[17]。

1.2.4 hub 基因的識別 根據(jù)中心模塊中每個基因的模塊連通性和臨床特征關(guān)系選擇候選hub 基因。模塊連通性定義為基因間Pearson 相關(guān)性（Module Membership）的絕對值。臨床性狀關(guān)系定義為每個基因與性狀之間Pearson 相關(guān)性（Gene Significance）的絕對值。將候選hub 基因的Module Membership設(shè)置為＞0.8，Gene Significance 設(shè)置為＞0.4。同時選擇中心模塊中的所有基因，應(yīng)用檢索交互基因的搜索工具（Search Tool for the Retrieval of Interacting Genes，STRING；https：//STRING db.org/）構(gòu)建蛋白質(zhì)互作（Protein protein interaction，PPI）網(wǎng)絡(luò)并尋找中心節(jié)點[18]。節(jié)點連通性＞22 的基因被認(rèn)為是中心節(jié)點。使 用Cytoscape（https：/Cyto-scape.org/）展示網(wǎng)絡(luò)并使用R 包“VennDiagram”進(jìn)行Venn 分析以比較PPI 網(wǎng)絡(luò)中的中心節(jié)點和候選hub 基因[19]。

1.2.5 hub 基因的驗證 使用TCGA 數(shù)據(jù)庫驗證hub基因的免疫相關(guān)性。首先，根據(jù)CIBERSORT 的計算結(jié)果提取出每個樣本的CD8+T 細(xì)胞比例。使用R 包“ggpubr”計算CD8+T 細(xì)胞浸潤水平與hub基因表達(dá)之間的Spearman 相關(guān)性，并使用R 包“ggstatsplot”對結(jié)果進(jìn)行比較。同時檢索腫瘤-免疫系統(tǒng)相互作用數(shù)據(jù)庫（Tumor Immune System Interactions Database，TISIDB；http：//cis.hku.hk/TISIDB）以確定hub 基因和TIIC 之間的Spearman相關(guān)性[20]。最后應(yīng)用R 包“heatmap”將這些結(jié)果以熱圖形式呈現(xiàn)出來。

1.2.6 免疫因子鑒定 hub 基因和不同免疫因子之間的Spearman 相關(guān)性來自TISIDB 數(shù)據(jù)庫，其中包括免疫抑制因子、免疫刺激因子、趨化因子和受體。然后使用R 包“heatmap”構(gòu)建熱圖，并使用STRING和Cytoscape 選擇與相關(guān)免疫因子的平均相關(guān)系數(shù)大于0.6 的hub 基因構(gòu)建網(wǎng)絡(luò)[19]。

1.2.7 hub 基因的預(yù)后分析 通過R 包“survminer”找到hub 基因表達(dá)量的最優(yōu)截斷值，并將患者分為高表達(dá)組和低表達(dá)組。然后通過R 包“survival”進(jìn)行Kaplan-Meier 分析。同時采用log-rank 檢驗對hub 基因進(jìn)行單因素生存分析，當(dāng)P＜0.05 時差異具有統(tǒng)計學(xué)意義。并將有統(tǒng)計學(xué)意義的hub 基因進(jìn)行多因素Cox 生存分析，鑒別出對CSCC 預(yù)后有獨立影響因素的基因。

1.2.8 基因集富集分析（GSEA） GSEA 是一種基于基因集的富集分析方法[21]。根據(jù)基因表達(dá)的中值，將樣本分為兩組，并進(jìn)行“c2.cp.kegg.v7.0.symbols”基因集富集分析，以P＜0.05 和q＜0.05 作為統(tǒng)計顯著性指標(biāo)。使用R 包“ggplot2”和“clusterProfiler”對富集途徑進(jìn)行可視化。在圓圖中，本研究只顯示富集過程的部分核心基因。

2 結(jié)果

2.1 CSCC 的加權(quán)基因共表達(dá)網(wǎng)絡(luò) 使用R 包“WGCNA”構(gòu)建5 605 個基因的共表達(dá)網(wǎng)絡(luò)，并對304 個CSCC 樣本進(jìn)行聚類。為構(gòu)建無標(biāo)度網(wǎng)絡(luò)，選擇β=4（無標(biāo)度R2=0.85）作為軟閾值功率（見圖1A、1B），并使用動態(tài)混合切割技術(shù)構(gòu)建了層次聚類樹，其中樹上的每片葉子代表一個基因，具有相似表達(dá)數(shù)據(jù)的基因靠一起，形成樹的一個分支，代表一個基因模塊，共生成10 個模塊（見圖1C）。

圖1 選擇合適的β值構(gòu)造層次聚類樹

2.2 識別中心模塊并進(jìn)行富集分析 模塊特征基因與T 細(xì)胞浸潤的相關(guān)性見圖2A。在這10 個模塊中黑色模塊與CD8+T 高度相關(guān)（r=0.4，P=2.3e-13）。黃綠色模塊與未活化記憶性CD4+T 淋巴細(xì)胞也有相對較高的相關(guān)性（r=0.33，P=4.6e-9）。其他模塊與T 細(xì)胞之間的相關(guān)性小于0.3。將與CD8+T 細(xì)胞高度相關(guān)的黑色模塊作為中心模塊，并使用R 包“GOplot”分析該模塊中所含基因，以進(jìn)行通路和過程富集分析。前20 個最高富集項中絕大多數(shù)與免疫相關(guān)（見圖2B），其中前3 個最高富集項為免疫系統(tǒng)的過程（Immune system process）、對刺激的調(diào)節(jié)反應(yīng)（Regulation of response to stimulus）和免疫應(yīng)答（Immune response）。

圖2 中心模塊和功能

2.3 hub 基因的識別和驗證 黑色模塊中的高度連接基因被視為與CD8+T 細(xì)胞浸潤水平相關(guān)的潛在因素。基于PPI 網(wǎng)絡(luò)，應(yīng)用intera-ction score＞0.7和degree＞22 的截止值（節(jié)點）從黑色模塊中確定34 個基因作為中心節(jié)點，同時使用Cytoscape將其可視化（見圖3A）。根據(jù)閾值標(biāo)準(zhǔn)（Module Membership＞0.8 和Gene Significance＞0.4），從黑色模塊中篩選68 個基因作為候選hub 基因（見圖3B）。最終選擇在以上兩種方法中均出現(xiàn)的10 個基因為hub 基因（CCL5、CCR5、CD2、CD3D、CD3E、CD8A、CXCR3、IFNG、PDCD1、PRF1），見圖3C。為研究這些hub 基因與CD8+T 細(xì)胞之間的關(guān)系，使用R 軟件分析TCGA-CSCC 數(shù)據(jù)集中hub 基因的表達(dá)數(shù)據(jù)與CD8+T 細(xì)胞的相關(guān)性。結(jié)果顯示，10 個基因的mRNA 表達(dá)量與CD8+T 細(xì)胞的浸潤水平呈顯著正相關(guān)（所有基因的相關(guān)系數(shù)最小為0.47），見圖4A。作為一個示例，圖4B 顯示了CD3E 表達(dá)量與CD8+T細(xì)胞浸潤水平的關(guān)系。查詢TISIDB 數(shù)據(jù)庫，以獲得腫瘤浸潤淋巴細(xì)胞的豐度與hub 基因表達(dá)之間的Spearman 相關(guān)值。結(jié)果顯示hub 基因與腫瘤浸潤淋巴細(xì)胞呈正相關(guān)?；罨疌D8+T 細(xì)胞（Act CD8）和效應(yīng)器記憶CD8+T 細(xì)胞（Tem CD8）的相關(guān)值最高（見圖4C）。此外，以TISIDB 數(shù)據(jù)庫為基礎(chǔ)，篩選出33 個免疫因子，這些免疫因子與hub 基因的平均相關(guān)性大于0.6?；趆ub 基因與33 個免疫因子，又構(gòu)建了免疫浸潤互作網(wǎng)絡(luò)，以此來探索hub基因的免疫浸潤機(jī)制（見圖4D）。而相關(guān)性熱圖詳細(xì)展示了hub 基因與不同免疫因子之間的關(guān)聯(lián)性（見圖5）。

圖3 hub 基因的鑒定

2.4 預(yù)后生物標(biāo)志物的鑒定 使用單因素Cox 回歸分析探索10 個hub 基因與CESC 患者預(yù)后之間的關(guān)系，發(fā)現(xiàn)除PRF1 與IFNG 外，其余8 個hub 基因均與患者預(yù)后有顯著的相關(guān)性（見圖6A）。而多因素Cox 回歸分析結(jié)果表明，僅CCR5（HR=1.61，P=0.03）與CD3E（HR=0.5，P=0.04）與預(yù)后顯著相關(guān)（見圖6B）。Kaplan-Meier 生存分析結(jié)果進(jìn)一步表明，與低表達(dá)組相比，10 個hub 基因高表達(dá)組患者有更好的總生存期（OS），見圖7，這與單因素Cox回歸分析結(jié)果基本相符。通過觀察發(fā)現(xiàn)多因素分析結(jié)果顯示CCR5 為CSCC 的危險因素，這與單因素分析和Kaplan-Meier 生存分析的結(jié)果完全相反。因此剔除CCR5 后選擇CD3E 作為進(jìn)一步分析的預(yù)后生物標(biāo)志物。

圖4 hub 基因驗證及蛋白質(zhì)互作網(wǎng)絡(luò)圖構(gòu)建

圖5 免疫因子和hub 基因之間相關(guān)性熱圖

圖6 hub 基因的Cox 回歸分析森林圖

圖7 10 個hub 基的Kaplan-Meier 曲線

2.5 CD3E 的基因集富集分析 根據(jù)CD3E 表達(dá)中值，將TCGA 中的CSCC 樣本分為高表達(dá)組和低表達(dá)組進(jìn)行基因集富集分析。結(jié)果表明，CD3E高表達(dá)組樣本富集到的免疫相關(guān)通路共21 個（FDR＜0.05，q＜0.05）。其中4 個富集最顯著的途徑分別為細(xì)胞黏附分子、趨化因子信號通路、細(xì)胞因子-細(xì)胞因子受體相互作用和T 細(xì)胞受體信號通路（見圖8A、8B）。而低表達(dá)組沒有明顯的富集途徑。

圖8 GSEA 富集分析

3 討論

CSCC 是子宮頸癌的主要組織學(xué)亞型，早期CSCC 的預(yù)后較好，而對于已發(fā)生遠(yuǎn)處轉(zhuǎn)移的晚期子宮頸癌或復(fù)發(fā)病例，傳統(tǒng)治療手段的效果并不理想[22]。免疫檢查點抑制劑在CSCC 中的成功應(yīng)用增加了研究人員對探索免疫治療中潛在特異性免疫相關(guān)因子的興趣[2]。CD8+T 細(xì)胞在免疫治療中起著核心作用。在本研究中，我們鑒定了10 個hub基因，其表達(dá)量與CD8+T 細(xì)胞浸潤水平相關(guān)，提示這些基因可能是影響CESC 免疫浸潤的潛在因素。在這些hub 基因中，CD3E 被確定為CSCC 潛在的預(yù)后生物標(biāo)志物和治療靶點。

目前利用動物腫瘤模型、體外細(xì)胞系和臨床樣本技術(shù)研究子宮頸癌的分子基礎(chǔ)已取得重大進(jìn)展。然而，CSCC 具有十分復(fù)雜的微環(huán)境，需要利用更大的數(shù)據(jù)集來進(jìn)一步分析。我們使用基因表達(dá)矩陣構(gòu)建共表達(dá)網(wǎng)絡(luò)，計算T 細(xì)胞的浸潤水平，并確定相關(guān)性，以篩選與CD8+T 細(xì)胞最相關(guān)的基因。對所選模塊的基因集富集分析表明，它是一個與免疫高度相關(guān)的模塊。WGCNA 與PPI 網(wǎng)絡(luò)之間關(guān)聯(lián)最緊密的基因被認(rèn)為是hub 基因（CCL5、CCR5、CD2、CD3D、CD3E、CD8A、CXCR3、IFNG、PDCD1、PRF1）。在TISIDB 數(shù)據(jù)庫中查詢這10 個基因與免疫細(xì)胞之間的關(guān)系，發(fā)現(xiàn)這些基因的表達(dá)與免疫細(xì)胞，尤其是CD8+T 細(xì)胞顯著正相關(guān)。以TISIDB 和STRING 數(shù)據(jù)庫為基礎(chǔ)，又構(gòu)建hub 基因和相關(guān)免疫因子之間的CD8+T 細(xì)胞浸潤網(wǎng)絡(luò)，以探索CSCC的免疫機(jī)制。分析結(jié)果證實hub 基因與CD8+T 細(xì)胞浸潤水平密切相關(guān)，并在免疫微環(huán)境中發(fā)揮重要作用。我們還對10 個hub 基因進(jìn)行了單因素、多因素Cox 回歸分析和Kaplan-Meier 分析。其中CCL5、CCR5、CD2、CD3D、CD3E、CD8A、CXCR3、PDCD1 與CSCC 患者預(yù)后有顯著的相關(guān)性，且hub基因高表達(dá)時患者的預(yù)后較好。而多因素Cox 回歸分析僅證實CD3E 低表達(dá)是CSCC 不良預(yù)后的獨立影響因素。結(jié)合這三種分析結(jié)果，我們選定CD3E為CSCC 患者潛在的預(yù)后指標(biāo)和治療靶點。

CD3 分子可以與T 細(xì)胞受體（TCR）結(jié)合，形成CD3/TCR 復(fù)合物并介導(dǎo)TCR 信號傳導(dǎo)和T 細(xì)胞分化[23，24]。由CD3E 基因編碼的CD3ε是CD3 的一個亞單位，與嚴(yán)重免疫缺陷有關(guān)，經(jīng)常被用作CD3 抗體的蛋白質(zhì)靶點[25，26]。研究發(fā)現(xiàn)，CD3E mRNA 表達(dá)水平較低的頭頸部鱗狀細(xì)胞癌患者復(fù)發(fā)風(fēng)險較高[27]。在肝癌、乳腺癌中也觀察到類似的結(jié)果，即CD3E 的高表達(dá)傾向于有更好的預(yù)后[28，29]。GSEA 結(jié)果還表明，CD3E 的高表達(dá)顯著富集于基質(zhì)相關(guān)信號通路，如T 細(xì)胞受體信號通路和趨化因子信號通路。這些結(jié)果表明，CD3E 可能還參與了腫瘤微環(huán)境從免疫型向代謝型的轉(zhuǎn)變。因此CD3E 有望成為CSCC 的預(yù)后指標(biāo)和新的免疫治療靶點。

總之，本研究嘗試使用WGCNA 和CIBERSORT算法來識別CSCC 潛在的CD8+T 細(xì)胞相關(guān)生物標(biāo)志物，并發(fā)現(xiàn)10 個與CSCC 預(yù)后相關(guān)的hub 基因。通過生物信息學(xué)驗證，CD3E 被確定為CSCC 免疫治療的潛在生物標(biāo)記物和免疫治療靶點。