加权基因共表达网络分析结肠癌免疫相关生物标志物

王苗苗 张睿哲 徐磊 武寒 吴淑华

滨州医学院附属医院病理科，滨州 256600

结直肠癌（colorectal cancer，CRC）是世界范围内常见的恶性肿瘤之一，其发病率位居全球癌症发病率第3位，病死率位居第2位［1］。当前，微阵列技术已被广泛用于探索癌症基因表达的变化［2-3］。新一代测序和高通量技术的快速进展使许多生物标志物被发现，对于进一步探讨结肠癌的潜在机制具有重要意义［4-7］。

研究表明，免疫系统在预防结肠癌的发生、发展和转移扩散方面发挥着重要作用［8］。目前，免疫疗法已成为结肠癌治疗的重要方法之一［9］。但是，免疫微环境对结肠癌的影响仍处于探索阶段，且缺乏切实可靠的免疫相关生物标志物。因此，了解免疫反应在结肠癌中的作用和机制、寻找新的免疫相关生物标志物对改善结肠癌的诊断、预后和预测，以及开发新的治疗方法具有重要意义。

加权基因共表达网络分析（weighted gene coexpression network analysis，WGCNA）作为一种常用的基因模块分析技术，已被广泛用于复杂疾病的分子标记或药物靶点的识别和筛选［10］。本研究筛选了癌症基因组图谱（The Cancer Genome Atlas，TCGA）数据集［11］和基因表达综合数据库（Gene Expression Omnibus database，GEO）数据集［12］中最关键的免疫共表达模块，并从中得到免疫相关生物标志物，进一步分析其在组织中的表达情况和可能行使的免疫相关功能，为临床研究及进一步试验提供参考。

资料与方法

1.方法

1.1.数据下载与整理 从GEO数据库中下载GSE41657数据集（https：//www.ncbi.nlm.nih.gov/geo/query/acc.cgi？acc=GSE41657）。GSE41657数据集包含12份正常上皮细胞样本和25份结肠癌样本。所选TCGA（https：//portal.gdc.cancer.gov/）结肠癌数据集共有437份样本，包括398份肿瘤组织和39份相应的正常邻近正常黏膜组织，同时下载了385例肿瘤患者的临床数据。从IMMPORT数据库（https：//www.immport.org/）下载免疫相关基因列表（共1 929个免疫相关基因）。用R软件对以上数据进行归一化处理，并提取GSE41657和TCGA数据集免疫相关基因表达矩阵。本试验所用数据资料研究时间为2022年11月至12月。

1.2.WGCNA分析 分别对GSE41657和TCGA数据集中免疫相关基因进行WGCNA。使用“WGCNA”R包构建结肠癌样本和正常黏膜样本中所有免疫相关基因的共表达网络。首先对基因表达数据进行过滤以减少异常值，选取合适的软阈值β加权系数构建基因网络，利用基因间的相位关系构建层次聚类树。随后利用基因的加权相关系数按表达模式进行分类，将具有相似表达模式的基因归为同一个模块；根据基因的表达模式将其分为不同的模块进行下一步分析。利用该加权相关系数将相关矩阵转化为邻接矩阵，进而转化为拓扑重叠矩阵（topological overlap matrix，TOM）。

最后，利用WGCNA算法计算各模块基因和样本组成的模块特征基因（module eigengenes，ME）的皮尔逊相关系数和P值。采用Pearson相关系数测定不同模块与临床特征的关系，选取适当相关系数模块进行后续分析。

1.3.共表达模块富集分析 分别对GSE41657和TCGA数据集中最显著差异模块进行功能［基因本体（gene ontology，GO）］与通路［京都基因与基因组百科全书（Kyoto Encyclopedia of Genes and Genomes，KEGG）］富集分析，寻找结肠癌中免疫相关差异基因所涉及的功能和通路。KEGG是一个数据库资源，用于大规模的分子数据收集，以了解生物系统和功能［13］。GO是分析基因的生物学过程（biological process，BP）、分子功能（molecular function，MF）和细胞成分（cellular components，CC）的生物信息学工具［14］。P<0.05为差异有统计学意义。

1.4.目标基因筛选 利用R软件对两个数据库中结肠癌与正常样本中全部基因进行差异基因筛选，调整后P<0.05和|log2 FC|≥2.0认为差异有统计学意义。利用Venn图在线网站（http：//bioinformatics.psb.ugent.be/webtools/Venn/）分别将得到的差异基因与两数据库中最显著共表达模块基因做交集，得到共表达免疫相关差异基因。

基于TCGA数据库临床信息，利用R软件筛选共表达免疫相关差异基因中的预后相关基因，通过TIMER数据库（https：//cistrome.shinyapps.io/timer/）研究其在结肠癌中与各免疫细胞、免疫标记物、免疫检查点的关系，并选取在结肠癌与免疫相关性最密切的基因作为目的基因用于进一步分析。

1.5.目标基因分析 根据TCGA数据库基因表达量对其进行差异表达可视化，通过GSEA数据库［15］（https：//www.gsea-msigdb.org/gsea/index.jsp）下载c2.cp.kegg.v7.1.symbols.gmt数据集，对其进行单基因KEGG富集分析，推测其可能在结肠癌中发挥免疫相关功能的机制。

2.统计学处理

使用R软件和R包ggplot2、corrplot、survminer等对数据进行可视化处理，P<0.05为差异有统计学意义。

结果

1.WGCNA分析

基因表达矩阵来自TCGA数据集和GSE41657数据集，并进行了标准化处理。用R软件删除免疫相关基因列表中的重复基因，提取免疫相关基因的表达矩阵，得到890个和671个免疫相关基因的表达矩阵，用于共表达网络的构建。本研究将所有样本纳入WGCNA。为了更符合无标度特性，选择TCGA数据集4和GSE41657数据集8作为构建共表达网络的β值，剪切高度设为0.25，对模块进行合并（图1A、D）。根据TOM检测，在两个数据集中都发现了3个基因模块（图1B、E）。从模块和性状的热图来看，TCGA数据库中的蓝色（307个基因）模块和GSE41657数据库中的蓝绿色模块（248个基因）与结肠癌的相关性最强。为了验证结果的可靠性，计算了两个模块中MM和GS基因的相关系数，即蓝绿色模块（COR=0.97，P=4.2E-153）和蓝色模块（COR=0.86，P=4.1E-91）（图1C、F）。因此，选择GEO的蓝绿色模块和TCGA的蓝色模块进行后续研究。

图1 GSE41657数据集和TCGA数据集的WGCNA。A：GSE41657数据集基因聚类树状图；B：GSE41657数据集临床特征与模块特征相关性热图；C：GSE41657数据集蓝绿色模块基因模块隶属度与基因显著性相关的散点图；D：TCGA数据集基因聚类树状图；E：TCGA数据集临床特征与模块特征相关性热图；F：TCGA数据集蓝色模块基因模块隶属度与基因显著性相关的散点图

2.共表达模块的富集分析

通过R软件对两个模块进行GO与KEGG 富集分析，并将KEGG富集结果的前30项及GO富集结果的前10项进行可视化，结果发现：在GO富集中，两个模块均富集在免疫细胞趋化、游走、活化、增殖、迁徙等重要免疫相关功能中；在KEGG富集中，两个模块均富集在免疫细胞受体信号通路、趋化因子信号通路、抗原处理和呈递、免疫细胞分化、MAPK信号通路、HIF-1信号通路、PI3K-Akt信号通路、NF-kappa B信号通路等重要的免疫相关通路。见图2。

图2 GSE41657数据集和TCGA数据集富集分析。A：GSE41657数据集GO富集分析；B：GSE41657数据集KEGG富集分析；C：TCGA数据集GO富集分析；D：TCGA数据集KEGG富集分析

3.目标基因筛选与免疫相关性研究

通过比较结肠癌样品和正常样品，筛选出所有的差异基因。在GSE41657数据集中，共筛选出4 712个差异基因，其中2 502个上调基因和2 210个下调基因。在TCGA数据集中，共筛选出3 572个差异基因，其中1 343个上调基因和2 229个下调基因。

通过对TCGA数据集和GSE41657数据集的差异基因和最显著模块中的基因进行交集，得到68个共表达的免疫相关差异基因（图3A）。基于TCGA数据库结肠临床数据进行预后筛选（Kaplan-Meier法），得到10个预后相关基因（NMB、SCG2、IL1A、ULBP2、INHBB、COLEC12、F2RL1、ANGPTL1、NR3C2、TNFRSF17）（图3B），通过TIMER数据库的筛选，最终选取结肠癌中与免疫最为密切的2个基因——COLEC12、ANGPTL1作为目标基因。COLEC12和ANGPTL1在结肠癌中与6种免疫细胞（B细胞、CD8+T细胞、CD4+T细胞、巨噬细胞、中性粒细胞和树突状细胞）均呈正相关关系，差异均有统计学意义（均P<0.05），其中COLEC12与CD4+T细胞、巨噬细胞、中性粒细胞和树突状细胞关系更为密切，ANGPTL1与巨噬细胞、树突状细胞的关系更为密切（图4A）。进一步对这两个基因在结肠癌中与免疫标志物的相关性分析发现，其结果与之对应的免疫细胞相关性基本相符（表1）。

表1 COLEC12、ANGPTL1与免疫标志物相关性分析

图3 Venn图和预后相关免疫基因。A：最显著模块基因与差异基因Venn图；B：预后相关免疫基因生存曲线

图4 ANGPTL1、COLEC12免疫相关性分析。A：ANGPTL1、COLEC12与免疫细胞的相关性；B：ANGPTL1、COLEC12与免疫检查点的相关性

此外，我们分析了COLEC12、ANGPTL1在结肠癌中与免疫检查点（CD274/PDCD1/CTLA4）的相关性，以探讨该基因作为免疫治疗相关基因的可能性，结果发现，COLEC12和ANGPTL1与免疫检查点均呈正相关关系，差异均有统计学意义（均P<0.05）。见图4B。

4.目标基因分析

根据TCGA数据库目标基因的表达情况，我们发现COLEC12、ANGPTL1在癌组织中的表达均低于正常组织，这一结果在GSE41657数据集中得到了验证（图5A、B）。通过单基因富集分析发现，COLEC12与ANGPTL1均富集到了诸多免疫相关和癌症相关的通路，如JAK STAT信号通路、Toll样受体信号通路、趋化因子信号通路、白细胞跨内皮迁移、黑色素瘤、自然杀伤细胞介导的细胞毒性、MAPK信号通路等，将关键通路进行可视化处理（图5C、D）。

图5 ANGPTL1、COLEC12表达及其机制分析。A、B：正常组织与癌组织中的基因表达；C、D：单基因KEGG富集分析

讨论

结肠癌是最常见的恶性肿瘤之一。目前，免疫疗法作为一种新的治疗方法已被用于部分结肠癌的治疗［16-17］。但是，目前免疫疗法对结肠癌的作用和影响尚不完全明确，且缺乏可靠的相关免疫生物标志物。因此，寻找新的可靠的免疫靶点和生物标志物显得尤为重要。

在本研究中，结肠癌组织和正常肠组织之间的差异基因来自两个mRNA微阵列数据集，GSE41657数据集和TCGA结肠癌数据集构建WGCNA用于结肠癌分析，此外，TCGA数据集还被用于结肠癌免疫相关的分析［18-19］。这充分说明两个数据库中数据的可靠性，但GSE41657数据集中的数据尚未被用于结肠癌的免疫相关分析，由于其样本量较小，采用TCGA数据集与其联合分析，进一步增加其准确性和可靠性，本研究采用与之前学者不同的研究方法，得到了具有价值的新结果。

通过WGCNA对两个数据集免疫相关基因表达矩阵进行共表达分析，分别筛选出两个最显著功能模块。我们对其进行富集分析，结果显示，这两个模块中的基因富集到的功能与通路相近，主要与免疫细胞趋化、游走、活化、增殖、迁徙等重要免疫相关功能相关，与免疫细胞受体信号通路、趋化因子信号通路、抗原处理和呈递、免疫细胞分化、MAPK信号通路、HIF-1信号通路、PI3K-Akt信号通路、NF-kappa B信号通路等重要免疫信号通路相关。近期研究显示，MAPK靶向疗法可能与免疫细胞产生协同作用，为开发新的肿瘤联合疗法提供了依据［20］。缺氧诱导因子1（hypoxia inducible factor-1，HIF-1）在缺氧肿瘤微环境中起着不可或缺的作用［21］。Shay等［22］发现，通过抑制HIF-1的表达可以阻止免疫功能正常小鼠中已建立的结肠炎相关结肠癌的进展。PI3K-Akt-mTOR抑制剂可以靶向癌细胞生物学过程，减弱免疫细胞效应功能，调节肿瘤微环境［23］。NF-κB被证明是可以将慢性炎症与癌症风险增加联系起来的关键分子［21］。由此可见，这些功能和通路均与癌症免疫密切相关。

为了得到关键模块中的预后相关差异基因，将关键模块与差异基因交集，得到68个共表达的免疫相关差异基因，并基于TCGA临床生存数据对其进行了筛选，结果得到10个预后相关基因：NMB、SCG2、IL1A、ULBP2、INHBB、COLEC12、F2RL1、ANGPTL1、NR3C2、TNFRSF17。通过TIMER数据库，选取结肠癌中与免疫最为密切的基因——COLEC12、ANGPTL1。本研究发现，这两个基因与各种免疫细胞均呈正相关关系，其中COLEC12与CD4+T细胞、巨噬细胞、中性粒细胞和树突状细胞关系更为密切，ANGPTL1与巨噬细胞、树突状细胞关系更为密切。有研究发现，肿瘤浸润的CD20+B淋巴细胞在结肠癌中具有良好的预后价值［24］，而CD8+T细胞作为肿瘤杀伤细胞，能抑制肿瘤的发展［25］。也有研究表明，CD4+T细胞在癌症免疫治疗中起着重要作用［26］，而中性粒细胞与淋巴细胞的比例是结肠癌患者的明确预测指标，肿瘤相关中性粒细胞能够促进肿瘤的发展［27］。肿瘤相关巨噬细胞是指被募集到肿瘤微环境中的巨噬细胞，其在肿瘤发生发展中发挥着重要作用［28-29］。Zhang等［30］发现，M2巨噬细胞通过调控COLEC12基因的表达来促进肾透明细胞癌的进展，从而影响肾透明细胞癌患者预后。本研究发现，预后相关基因COLEC12、ANGPTL1与这些免疫细胞密切相关，同时在这两个基因与免疫标志物相关性分析中验证了这一点。有研究显示，COLEC12可促进胃癌细胞增殖、迁移和侵袭，抑制胃癌细胞凋亡，其可能在抑制胃癌肿瘤免疫应答中发挥作用［31］。由此推测，其可能通过调控肿瘤免疫，在结肠癌发生发展过程中发挥作用。

此外，本研究还分析了其与癌症中免疫检查点（CD274/PDCD1/CTLA4）的相关性，探讨该基因作为免疫治疗相关基因的可能性，结果发现，这两个基因与免疫检查点密切相关并均呈正相关关系，提示其有成为免疫治疗相关基因的潜力。

目标基因COLEC12、ANGPTL1在结肠癌中的表达低于正常组织，其高表达提示不良预后。为了推测其具体功能和机制，进行单基因KEGG富集分析发现，COLEC12、ANGPTL1与多种免疫相关和癌症相关的通路——JAK STAT信号通路、Toll样受体信号通路、趋化因子信号通路、白细胞跨内皮迁移、黑色素瘤、天然杀伤细胞介导的细胞毒性、MAPK信号通路等——有关。TET2介导干扰素γ-JAK-STAT信号通路以控制趋化因子和细胞程序性死亡-配体1表达、淋巴细胞浸润和癌症免疫［32］。Toll样受体家族是肿瘤免疫研究的热点之一［33-34］。外泌体ANGPTL1通过调节Kupffer细胞分泌和阻止MMP9诱导血管渗漏来降低结直肠癌的肝转移［35］，COLEC12也被用于预测结肠腺癌患者的预后［36］。由此推测，COLEC12、ANGPTL1可能通过这些免疫相关通路调控结肠癌患者免疫微环境，并最终影响结肠癌患者的进展及预后。

综上所述，本研究基于TCGA和GEO数据库联合分析，运用WGCNA等方法，最终筛选出2个与结肠癌免疫最相关的基因——COLEC12、ANGPTL1，其可能在肿瘤免疫微环境中起着重要作用。本研究为结肠癌的发生发展机制以及肿瘤免疫环境的研究提供了新的研究思路和参考，筛选出的预后相关基因可能成为诊断结肠癌的生物标志物。但是，本文只是从生物信息水平进行了初步分析，这些基因在结肠癌中具体的功能机制尚不清楚，仍需扩大样本量分析并从多层面开展分子机制、细胞实验、动物实验等更深入的探究，进一步验证。同时，这也将成为课题组下一步研究的方向。

利益冲突所有作者均声明不存在利益冲突

作者贡献声明王苗苗：酝酿和设计试验，实施研究，采集数据，分析/解释数据，起草文章，对文章的知识性内容作批评性审阅，统计分析；张睿哲：酝酿和设计试验，实施研究，分析/解释数据，起草文章，对文章的知识性内容作批评性审阅，统计分析；徐磊、武寒：酝酿和设计试验，实施研究，采集数据，对文章的知识性内容作批评性审阅，统计分析；吴淑华：酝酿和设计试验，分析/解释数据，对文章的知识性内容作批评性审阅，统计分析获取研究经费，行政、技术或材料支持，指导，支持性贡献