加權(quán)基因共表達網(wǎng)絡(luò)分析結(jié)腸癌免疫相關(guān)生物標志物

王苗苗 張睿哲 徐磊 武寒 吳淑華

濱州醫(yī)學院附屬醫(yī)院病理科，濱州 256600

結(jié)直腸癌（colorectal cancer，CRC）是世界范圍內(nèi)常見的惡性腫瘤之一，其發(fā)病率位居全球癌癥發(fā)病率第3位，病死率位居第2位［1］。當前，微陣列技術(shù)已被廣泛用于探索癌癥基因表達的變化［2-3］。新一代測序和高通量技術(shù)的快速進展使許多生物標志物被發(fā)現(xiàn)，對于進一步探討結(jié)腸癌的潛在機制具有重要意義［4-7］。

研究表明，免疫系統(tǒng)在預(yù)防結(jié)腸癌的發(fā)生、發(fā)展和轉(zhuǎn)移擴散方面發(fā)揮著重要作用［8］。目前，免疫療法已成為結(jié)腸癌治療的重要方法之一［9］。但是，免疫微環(huán)境對結(jié)腸癌的影響仍處于探索階段，且缺乏切實可靠的免疫相關(guān)生物標志物。因此，了解免疫反應(yīng)在結(jié)腸癌中的作用和機制、尋找新的免疫相關(guān)生物標志物對改善結(jié)腸癌的診斷、預(yù)后和預(yù)測，以及開發(fā)新的治療方法具有重要意義。

加權(quán)基因共表達網(wǎng)絡(luò)分析（weighted gene coexpression network analysis，WGCNA）作為一種常用的基因模塊分析技術(shù)，已被廣泛用于復(fù)雜疾病的分子標記或藥物靶點的識別和篩選［10］。本研究篩選了癌癥基因組圖譜（The Cancer Genome Atlas，TCGA）數(shù)據(jù)集［11］和基因表達綜合數(shù)據(jù)庫（Gene Expression Omnibus database，GEO）數(shù)據(jù)集［12］中最關(guān)鍵的免疫共表達模塊，并從中得到免疫相關(guān)生物標志物，進一步分析其在組織中的表達情況和可能行使的免疫相關(guān)功能，為臨床研究及進一步試驗提供參考。

資料與方法

1.方法

1.1.數(shù)據(jù)下載與整理 從GEO數(shù)據(jù)庫中下載GSE41657數(shù)據(jù)集（https：//www.ncbi.nlm.nih.gov/geo/query/acc.cgi？acc=GSE41657）。GSE41657數(shù)據(jù)集包含12份正常上皮細胞樣本和25份結(jié)腸癌樣本。所選TCGA（https：//portal.gdc.cancer.gov/）結(jié)腸癌數(shù)據(jù)集共有437份樣本，包括398份腫瘤組織和39份相應(yīng)的正常鄰近正常黏膜組織，同時下載了385例腫瘤患者的臨床數(shù)據(jù)。從IMMPORT數(shù)據(jù)庫（https：//www.immport.org/）下載免疫相關(guān)基因列表（共1 929個免疫相關(guān)基因）。用R軟件對以上數(shù)據(jù)進行歸一化處理，并提取GSE41657和TCGA數(shù)據(jù)集免疫相關(guān)基因表達矩陣。本試驗所用數(shù)據(jù)資料研究時間為2022年11月至12月。

1.2.WGCNA分析 分別對GSE41657和TCGA數(shù)據(jù)集中免疫相關(guān)基因進行WGCNA。使用“WGCNA”R包構(gòu)建結(jié)腸癌樣本和正常黏膜樣本中所有免疫相關(guān)基因的共表達網(wǎng)絡(luò)。首先對基因表達數(shù)據(jù)進行過濾以減少異常值，選取合適的軟閾值β加權(quán)系數(shù)構(gòu)建基因網(wǎng)絡(luò)，利用基因間的相位關(guān)系構(gòu)建層次聚類樹。隨后利用基因的加權(quán)相關(guān)系數(shù)按表達模式進行分類，將具有相似表達模式的基因歸為同一個模塊；根據(jù)基因的表達模式將其分為不同的模塊進行下一步分析。利用該加權(quán)相關(guān)系數(shù)將相關(guān)矩陣轉(zhuǎn)化為鄰接矩陣，進而轉(zhuǎn)化為拓撲重疊矩陣（topological overlap matrix，TOM）。

最后，利用WGCNA算法計算各模塊基因和樣本組成的模塊特征基因（module eigengenes，ME）的皮爾遜相關(guān)系數(shù)和P值。采用Pearson相關(guān)系數(shù)測定不同模塊與臨床特征的關(guān)系，選取適當相關(guān)系數(shù)模塊進行后續(xù)分析。

1.3.共表達模塊富集分析 分別對GSE41657和TCGA數(shù)據(jù)集中最顯著差異模塊進行功能［基因本體（gene ontology，GO）］與通路［京都基因與基因組百科全書（Kyoto Encyclopedia of Genes and Genomes，KEGG）］富集分析，尋找結(jié)腸癌中免疫相關(guān)差異基因所涉及的功能和通路。KEGG是一個數(shù)據(jù)庫資源，用于大規(guī)模的分子數(shù)據(jù)收集，以了解生物系統(tǒng)和功能［13］。GO是分析基因的生物學過程（biological process，BP）、分子功能（molecular function，MF）和細胞成分（cellular components，CC）的生物信息學工具［14］。P<0.05為差異有統(tǒng)計學意義。

1.4.目標基因篩選 利用R軟件對兩個數(shù)據(jù)庫中結(jié)腸癌與正常樣本中全部基因進行差異基因篩選，調(diào)整后P<0.05和|log2 FC|≥2.0認為差異有統(tǒng)計學意義。利用Venn圖在線網(wǎng)站（http：//bioinformatics.psb.ugent.be/webtools/Venn/）分別將得到的差異基因與兩數(shù)據(jù)庫中最顯著共表達模塊基因做交集，得到共表達免疫相關(guān)差異基因。

基于TCGA數(shù)據(jù)庫臨床信息，利用R軟件篩選共表達免疫相關(guān)差異基因中的預(yù)后相關(guān)基因，通過TIMER數(shù)據(jù)庫（https：//cistrome.shinyapps.io/timer/）研究其在結(jié)腸癌中與各免疫細胞、免疫標記物、免疫檢查點的關(guān)系，并選取在結(jié)腸癌與免疫相關(guān)性最密切的基因作為目的基因用于進一步分析。

1.5.目標基因分析 根據(jù)TCGA數(shù)據(jù)庫基因表達量對其進行差異表達可視化，通過GSEA數(shù)據(jù)庫［15］（https：//www.gsea-msigdb.org/gsea/index.jsp）下載c2.cp.kegg.v7.1.symbols.gmt數(shù)據(jù)集，對其進行單基因KEGG富集分析，推測其可能在結(jié)腸癌中發(fā)揮免疫相關(guān)功能的機制。

2.統(tǒng)計學處理

使用R軟件和R包ggplot2、corrplot、survminer等對數(shù)據(jù)進行可視化處理，P<0.05為差異有統(tǒng)計學意義。

結(jié)果

1.WGCNA分析

基因表達矩陣來自TCGA數(shù)據(jù)集和GSE41657數(shù)據(jù)集，并進行了標準化處理。用R軟件刪除免疫相關(guān)基因列表中的重復(fù)基因，提取免疫相關(guān)基因的表達矩陣，得到890個和671個免疫相關(guān)基因的表達矩陣，用于共表達網(wǎng)絡(luò)的構(gòu)建。本研究將所有樣本納入WGCNA。為了更符合無標度特性，選擇TCGA數(shù)據(jù)集4和GSE41657數(shù)據(jù)集8作為構(gòu)建共表達網(wǎng)絡(luò)的β值，剪切高度設(shè)為0.25，對模塊進行合并（圖1A、D）。根據(jù)TOM檢測，在兩個數(shù)據(jù)集中都發(fā)現(xiàn)了3個基因模塊（圖1B、E）。從模塊和性狀的熱圖來看，TCGA數(shù)據(jù)庫中的藍色（307個基因）模塊和GSE41657數(shù)據(jù)庫中的藍綠色模塊（248個基因）與結(jié)腸癌的相關(guān)性最強。為了驗證結(jié)果的可靠性，計算了兩個模塊中MM和GS基因的相關(guān)系數(shù)，即藍綠色模塊（COR=0.97，P=4.2E-153）和藍色模塊（COR=0.86，P=4.1E-91）（圖1C、F）。因此，選擇GEO的藍綠色模塊和TCGA的藍色模塊進行后續(xù)研究。

圖1 GSE41657數(shù)據(jù)集和TCGA數(shù)據(jù)集的WGCNA。A：GSE41657數(shù)據(jù)集基因聚類樹狀圖；B：GSE41657數(shù)據(jù)集臨床特征與模塊特征相關(guān)性熱圖；C：GSE41657數(shù)據(jù)集藍綠色模塊基因模塊隸屬度與基因顯著性相關(guān)的散點圖；D：TCGA數(shù)據(jù)集基因聚類樹狀圖；E：TCGA數(shù)據(jù)集臨床特征與模塊特征相關(guān)性熱圖；F：TCGA數(shù)據(jù)集藍色模塊基因模塊隸屬度與基因顯著性相關(guān)的散點圖

2.共表達模塊的富集分析

通過R軟件對兩個模塊進行GO與KEGG 富集分析，并將KEGG富集結(jié)果的前30項及GO富集結(jié)果的前10項進行可視化，結(jié)果發(fā)現(xiàn)：在GO富集中，兩個模塊均富集在免疫細胞趨化、游走、活化、增殖、遷徙等重要免疫相關(guān)功能中；在KEGG富集中，兩個模塊均富集在免疫細胞受體信號通路、趨化因子信號通路、抗原處理和呈遞、免疫細胞分化、MAPK信號通路、HIF-1信號通路、PI3K-Akt信號通路、NF-kappa B信號通路等重要的免疫相關(guān)通路。見圖2。

圖2 GSE41657數(shù)據(jù)集和TCGA數(shù)據(jù)集富集分析。A：GSE41657數(shù)據(jù)集GO富集分析；B：GSE41657數(shù)據(jù)集KEGG富集分析；C：TCGA數(shù)據(jù)集GO富集分析；D：TCGA數(shù)據(jù)集KEGG富集分析

3.目標基因篩選與免疫相關(guān)性研究

通過比較結(jié)腸癌樣品和正常樣品，篩選出所有的差異基因。在GSE41657數(shù)據(jù)集中，共篩選出4 712個差異基因，其中2 502個上調(diào)基因和2 210個下調(diào)基因。在TCGA數(shù)據(jù)集中，共篩選出3 572個差異基因，其中1 343個上調(diào)基因和2 229個下調(diào)基因。

通過對TCGA數(shù)據(jù)集和GSE41657數(shù)據(jù)集的差異基因和最顯著模塊中的基因進行交集，得到68個共表達的免疫相關(guān)差異基因（圖3A）?；赥CGA數(shù)據(jù)庫結(jié)腸臨床數(shù)據(jù)進行預(yù)后篩選（Kaplan-Meier法），得到10個預(yù)后相關(guān)基因（NMB、SCG2、IL1A、ULBP2、INHBB、COLEC12、F2RL1、ANGPTL1、NR3C2、TNFRSF17）（圖3B），通過TIMER數(shù)據(jù)庫的篩選，最終選取結(jié)腸癌中與免疫最為密切的2個基因——COLEC12、ANGPTL1作為目標基因。COLEC12和ANGPTL1在結(jié)腸癌中與6種免疫細胞（B細胞、CD8+T細胞、CD4+T細胞、巨噬細胞、中性粒細胞和樹突狀細胞）均呈正相關(guān)關(guān)系，差異均有統(tǒng)計學意義（均P<0.05），其中COLEC12與CD4+T細胞、巨噬細胞、中性粒細胞和樹突狀細胞關(guān)系更為密切，ANGPTL1與巨噬細胞、樹突狀細胞的關(guān)系更為密切（圖4A）。進一步對這兩個基因在結(jié)腸癌中與免疫標志物的相關(guān)性分析發(fā)現(xiàn)，其結(jié)果與之對應(yīng)的免疫細胞相關(guān)性基本相符（表1）。

表1 COLEC12、ANGPTL1與免疫標志物相關(guān)性分析

圖3 Venn圖和預(yù)后相關(guān)免疫基因。A：最顯著模塊基因與差異基因Venn圖；B：預(yù)后相關(guān)免疫基因生存曲線

圖4 ANGPTL1、COLEC12免疫相關(guān)性分析。A：ANGPTL1、COLEC12與免疫細胞的相關(guān)性；B：ANGPTL1、COLEC12與免疫檢查點的相關(guān)性

此外，我們分析了COLEC12、ANGPTL1在結(jié)腸癌中與免疫檢查點（CD274/PDCD1/CTLA4）的相關(guān)性，以探討該基因作為免疫治療相關(guān)基因的可能性，結(jié)果發(fā)現(xiàn)，COLEC12和ANGPTL1與免疫檢查點均呈正相關(guān)關(guān)系，差異均有統(tǒng)計學意義（均P<0.05）。見圖4B。

4.目標基因分析

根據(jù)TCGA數(shù)據(jù)庫目標基因的表達情況，我們發(fā)現(xiàn)COLEC12、ANGPTL1在癌組織中的表達均低于正常組織，這一結(jié)果在GSE41657數(shù)據(jù)集中得到了驗證（圖5A、B）。通過單基因富集分析發(fā)現(xiàn)，COLEC12與ANGPTL1均富集到了諸多免疫相關(guān)和癌癥相關(guān)的通路，如JAK STAT信號通路、Toll樣受體信號通路、趨化因子信號通路、白細胞跨內(nèi)皮遷移、黑色素瘤、自然殺傷細胞介導(dǎo)的細胞毒性、MAPK信號通路等，將關(guān)鍵通路進行可視化處理（圖5C、D）。

圖5 ANGPTL1、COLEC12表達及其機制分析。A、B：正常組織與癌組織中的基因表達；C、D：單基因KEGG富集分析

討論

結(jié)腸癌是最常見的惡性腫瘤之一。目前，免疫療法作為一種新的治療方法已被用于部分結(jié)腸癌的治療［16-17］。但是，目前免疫療法對結(jié)腸癌的作用和影響尚不完全明確，且缺乏可靠的相關(guān)免疫生物標志物。因此，尋找新的可靠的免疫靶點和生物標志物顯得尤為重要。

在本研究中，結(jié)腸癌組織和正常腸組織之間的差異基因來自兩個mRNA微陣列數(shù)據(jù)集，GSE41657數(shù)據(jù)集和TCGA結(jié)腸癌數(shù)據(jù)集構(gòu)建WGCNA用于結(jié)腸癌分析，此外，TCGA數(shù)據(jù)集還被用于結(jié)腸癌免疫相關(guān)的分析［18-19］。這充分說明兩個數(shù)據(jù)庫中數(shù)據(jù)的可靠性，但GSE41657數(shù)據(jù)集中的數(shù)據(jù)尚未被用于結(jié)腸癌的免疫相關(guān)分析，由于其樣本量較小，采用TCGA數(shù)據(jù)集與其聯(lián)合分析，進一步增加其準確性和可靠性，本研究采用與之前學者不同的研究方法，得到了具有價值的新結(jié)果。

通過WGCNA對兩個數(shù)據(jù)集免疫相關(guān)基因表達矩陣進行共表達分析，分別篩選出兩個最顯著功能模塊。我們對其進行富集分析，結(jié)果顯示，這兩個模塊中的基因富集到的功能與通路相近，主要與免疫細胞趨化、游走、活化、增殖、遷徙等重要免疫相關(guān)功能相關(guān)，與免疫細胞受體信號通路、趨化因子信號通路、抗原處理和呈遞、免疫細胞分化、MAPK信號通路、HIF-1信號通路、PI3K-Akt信號通路、NF-kappa B信號通路等重要免疫信號通路相關(guān)。近期研究顯示，MAPK靶向療法可能與免疫細胞產(chǎn)生協(xié)同作用，為開發(fā)新的腫瘤聯(lián)合療法提供了依據(jù)［20］。缺氧誘導(dǎo)因子1（hypoxia inducible factor-1，HIF-1）在缺氧腫瘤微環(huán)境中起著不可或缺的作用［21］。Shay等［22］發(fā)現(xiàn)，通過抑制HIF-1的表達可以阻止免疫功能正常小鼠中已建立的結(jié)腸炎相關(guān)結(jié)腸癌的進展。PI3K-Akt-mTOR抑制劑可以靶向癌細胞生物學過程，減弱免疫細胞效應(yīng)功能，調(diào)節(jié)腫瘤微環(huán)境［23］。NF-κB被證明是可以將慢性炎癥與癌癥風險增加聯(lián)系起來的關(guān)鍵分子［21］。由此可見，這些功能和通路均與癌癥免疫密切相關(guān)。

為了得到關(guān)鍵模塊中的預(yù)后相關(guān)差異基因，將關(guān)鍵模塊與差異基因交集，得到68個共表達的免疫相關(guān)差異基因，并基于TCGA臨床生存數(shù)據(jù)對其進行了篩選，結(jié)果得到10個預(yù)后相關(guān)基因：NMB、SCG2、IL1A、ULBP2、INHBB、COLEC12、F2RL1、ANGPTL1、NR3C2、TNFRSF17。通過TIMER數(shù)據(jù)庫，選取結(jié)腸癌中與免疫最為密切的基因——COLEC12、ANGPTL1。本研究發(fā)現(xiàn)，這兩個基因與各種免疫細胞均呈正相關(guān)關(guān)系，其中COLEC12與CD4+T細胞、巨噬細胞、中性粒細胞和樹突狀細胞關(guān)系更為密切，ANGPTL1與巨噬細胞、樹突狀細胞關(guān)系更為密切。有研究發(fā)現(xiàn)，腫瘤浸潤的CD20+B淋巴細胞在結(jié)腸癌中具有良好的預(yù)后價值［24］，而CD8+T細胞作為腫瘤殺傷細胞，能抑制腫瘤的發(fā)展［25］。也有研究表明，CD4+T細胞在癌癥免疫治療中起著重要作用［26］，而中性粒細胞與淋巴細胞的比例是結(jié)腸癌患者的明確預(yù)測指標，腫瘤相關(guān)中性粒細胞能夠促進腫瘤的發(fā)展［27］。腫瘤相關(guān)巨噬細胞是指被募集到腫瘤微環(huán)境中的巨噬細胞，其在腫瘤發(fā)生發(fā)展中發(fā)揮著重要作用［28-29］。Zhang等［30］發(fā)現(xiàn)，M2巨噬細胞通過調(diào)控COLEC12基因的表達來促進腎透明細胞癌的進展，從而影響腎透明細胞癌患者預(yù)后。本研究發(fā)現(xiàn)，預(yù)后相關(guān)基因COLEC12、ANGPTL1與這些免疫細胞密切相關(guān)，同時在這兩個基因與免疫標志物相關(guān)性分析中驗證了這一點。有研究顯示，COLEC12可促進胃癌細胞增殖、遷移和侵襲，抑制胃癌細胞凋亡，其可能在抑制胃癌腫瘤免疫應(yīng)答中發(fā)揮作用［31］。由此推測，其可能通過調(diào)控腫瘤免疫，在結(jié)腸癌發(fā)生發(fā)展過程中發(fā)揮作用。

此外，本研究還分析了其與癌癥中免疫檢查點（CD274/PDCD1/CTLA4）的相關(guān)性，探討該基因作為免疫治療相關(guān)基因的可能性，結(jié)果發(fā)現(xiàn)，這兩個基因與免疫檢查點密切相關(guān)并均呈正相關(guān)關(guān)系，提示其有成為免疫治療相關(guān)基因的潛力。

目標基因COLEC12、ANGPTL1在結(jié)腸癌中的表達低于正常組織，其高表達提示不良預(yù)后。為了推測其具體功能和機制，進行單基因KEGG富集分析發(fā)現(xiàn)，COLEC12、ANGPTL1與多種免疫相關(guān)和癌癥相關(guān)的通路——JAK STAT信號通路、Toll樣受體信號通路、趨化因子信號通路、白細胞跨內(nèi)皮遷移、黑色素瘤、天然殺傷細胞介導(dǎo)的細胞毒性、MAPK信號通路等——有關(guān)。TET2介導(dǎo)干擾素γ-JAK-STAT信號通路以控制趨化因子和細胞程序性死亡-配體1表達、淋巴細胞浸潤和癌癥免疫［32］。Toll樣受體家族是腫瘤免疫研究的熱點之一［33-34］。外泌體ANGPTL1通過調(diào)節(jié)Kupffer細胞分泌和阻止MMP9誘導(dǎo)血管滲漏來降低結(jié)直腸癌的肝轉(zhuǎn)移［35］，COLEC12也被用于預(yù)測結(jié)腸腺癌患者的預(yù)后［36］。由此推測，COLEC12、ANGPTL1可能通過這些免疫相關(guān)通路調(diào)控結(jié)腸癌患者免疫微環(huán)境，并最終影響結(jié)腸癌患者的進展及預(yù)后。

綜上所述，本研究基于TCGA和GEO數(shù)據(jù)庫聯(lián)合分析，運用WGCNA等方法，最終篩選出2個與結(jié)腸癌免疫最相關(guān)的基因——COLEC12、ANGPTL1，其可能在腫瘤免疫微環(huán)境中起著重要作用。本研究為結(jié)腸癌的發(fā)生發(fā)展機制以及腫瘤免疫環(huán)境的研究提供了新的研究思路和參考，篩選出的預(yù)后相關(guān)基因可能成為診斷結(jié)腸癌的生物標志物。但是，本文只是從生物信息水平進行了初步分析，這些基因在結(jié)腸癌中具體的功能機制尚不清楚，仍需擴大樣本量分析并從多層面開展分子機制、細胞實驗、動物實驗等更深入的探究，進一步驗證。同時，這也將成為課題組下一步研究的方向。

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

作者貢獻聲明王苗苗：醞釀和設(shè)計試驗，實施研究，采集數(shù)據(jù)，分析/解釋數(shù)據(jù)，起草文章，對文章的知識性內(nèi)容作批評性審閱，統(tǒng)計分析；張睿哲：醞釀和設(shè)計試驗，實施研究，分析/解釋數(shù)據(jù)，起草文章，對文章的知識性內(nèi)容作批評性審閱，統(tǒng)計分析；徐磊、武寒：醞釀和設(shè)計試驗，實施研究，采集數(shù)據(jù)，對文章的知識性內(nèi)容作批評性審閱，統(tǒng)計分析；吳淑華：醞釀和設(shè)計試驗，分析/解釋數(shù)據(jù)，對文章的知識性內(nèi)容作批評性審閱，統(tǒng)計分析獲取研究經(jīng)費，行政、技術(shù)或材料支持，指導(dǎo)，支持性貢獻