结肠癌核心基因和独立预后因子筛选的生物信息学分析

刘 迁,祁国萍,于华裔,戴宇阳,陆文斌,金建华

（江苏大学附属武进医院 徐州医科大学武进临床学院肿瘤科，江苏 常州 213017）

结肠癌是世界上致死率排名第三位的疾病［1］。目前存在一些临床常用的结肠癌标志物［2］，但用于诊疗结肠癌的标志物仍然较少，故探寻潜在的结肠癌标志物仍是结肠癌研究的热点和重点。近年来，高通量测序和基因芯片技术被广泛地应用于生命科学领域［3-4］。生物信息学是分析现存大量生物数据的重要工具，通过分析数据中存在的潜在重要核心基因或预后因子［5-6］，挖掘潜在肿瘤标志物或治疗靶点［7］。

当前采用生物信息学分析的方法系统挖掘肿瘤致病基因及其调控机制将对肿瘤研究起到巨大的推动作用。基于此，有研究［5］通过分析基因表达综合（Gene Expression Omnibus，GEO）数据库，联合多种生物信息学分析方法，对GEO基因芯片数据库分析筛选结肠癌枢纽基因，对结肠癌的发病机制研究提供理论基础。本研究通过对癌症基因组图谱（The Cancer Genome Atlas，TCGA）数据库和GEO数据库独立分析及采用2种不同的数据差异分析方法进行联合分析，获得13个共有基因，进一步行蛋白-蛋白互作 （protein-protein interaction，PPI）网络构建，并识别出具有相互联系的7个核心基因。本研究对这7个核心基因进行基因本体（Gene Ontology，GO）和京都基因和基因组百科全书（Kyoto Encylopaedia of Genes and Genomes，KEGG）富集分析，还对这7个核心基因进行泛癌（32个不同肿瘤样本）分析。通过对核心基因分别进行结肠癌总生存的预后风险因素分析（单因素和双因素COX分析），筛选出关键基因——AXIN2基因的表达对结肠癌患者总生存的预后具有显著意义，为结肠癌的研究提供依据。

1 资料与方法

1.1 原始数据来源从TCGA数据库（https：//portal.gdc.cancer.gov/）中下载结肠癌转录组数据和临床数据，并进行标准化处理。从GEO数据库（http：//www.ncbi.nlm.nih.gov/geo/）中下载高通量测序数据集GSE37182基因芯片数据。

1.2 差异分析筛选采用R/Bioconductor软件包中的EdgeR用于处理来源于TCGA的数据。采用R软件中的“limma”［9］、“edgeR”［10］和“ggplot2”［11］程序包进行数据处理分析和作图观察，倍数（fold change，FC）以∣logFC∣大于1及矫正后P＜0.05作为筛选条件进行筛选。GEO2R（http：//www.ncbi.nlm.nih.gov/geo/geo2r/）可用于分析GEO数据集中不同数据组的差异表达水平。采用GEO2R在线分析工具分别分析并筛选出GSE37182的差异表达基因 （differentially expressed genes，DEGs），以∣logFC∣＞1，调整后以P＜0.05作为识别DEGs的阈值。采用R软件中的“ggplot2”和“pheatmap”程序包分析DEGs并进行可视化，分别生成火山图和热图。

1.3 加 权 基因 共 表 达 网 络 分 析［12］（weighted correlation network analysis，WGCNA）采 用R Studio软件中“WGCNA”函数包，对数据进行预处理，采用“goodsamplegenes（）”函数检查数据是否有异常基因，选择合适的阈值剔除离群值。加载临床特征数据并绘制基因-特征热度聚类图，完成共表达网络的构建和模块的选择。设定的模块基因数为50个。表达相似的基因进行聚类形成树状图，树状图中的每片叶子代表1个基因，每个分支代表1组表达相近或相似的基因组，每组不超过50个基因。

1.4 维恩分析WGCNA方法分析获得TCGA数据中的magenta+salmon颜色模块对应的基因组合（TCGA_salmon+magenta）、GEO数据中的blue+pink颜色模块对应的基因组合（GEO_blue+pink）、差异分析获得的TCGA数据中的差异基因（TCGA_diff）和GEO数据中获得的差异基因（GEO_diff），将这4组基因取交集，获得13个共有基因。

1.5 PPI和Hub基因筛选为了探讨维恩分析结果中PPI，本研究采用STRING数据库（https：//string-db.org/）进行分析比对，获得已知的蛋白之间的相互作用图，构建PPI网络。为了进一步挖掘模块中的关键基因节点，将PPI网络导入Cytoscape3.6.1软件可视化，通过cytoHubba插件中的最大团中心性 （maximal clique centrality，MCC）算法筛选Hub基因。

1.6 GO和KEGG富集分析通过在线数据库DAVID（http：//dabid.ncifcrf.gov/）对核心基因进行GO富集分析［13］和KEGG富集分析［14］，分析核心基因主要参与的生物学功能及参与的重要信号通路等。

1.7 泛癌基因的免疫细胞亚型分析首先在UCSC Xena网站下载TCGA数据库中的免疫亚型数 据 包 （https：//xenabrowser.net/datapages/？cohort=TCGA%20Pan-Cancer%20（PANCAN）&remove Hub=https% 3A% 2F% 2Fxena.treehouse.gi.ucsc.edu%3A 443），通过“perl”脚本对数据进行处理，采用R软件中“limma”、“ggplot2”［15］和“reshape2”程序包对数据进行分析，分析核心基因的表达与免疫细胞亚型的相关性［16］。

1.8 核心基因的泛癌热图和相关性分析采用R软件中的“pheatmap”程序包对处理后的数据进行分析获得热图。采用“corrplot”程序包分析核心基因之间相关性，通过每个基因在不同癌种中的相对表达量进行相关性分析，此程序里会删除正常样品，只分析不同肿瘤组织中的基因表达。

1.9 核心基因与肿瘤微环境的相关性分析采用R软件中的“estimate”和“corrplot”程序包进行数据分析统计，并分析基因表达与基质细胞评分之间的相关性、基因表达与免疫细胞评分之间的相关性和本次检测的基因表达与综合（基质细胞和免疫细胞的总和）评分之间的相关性［17］。

1.10 核心基因的风险回归模型分析下载TCGA数据库中临床数据样本，采用“perl”脚本进行数据处理，再通过R软件中的“survival”程序包进行基因的单因素或多因素COX分析，采用R软件中的“survminer”程序包绘制森林图［18］。

1.11 实时荧光定量PCR（real-time fluorescence quantitative PCR，RT-qPCR）、Western blotting和免疫组织化学法检测本研究对江苏大学附属武进医院收集的结肠癌组织和癌旁组织的临床标本分别进行组织RNA提取、反转录及RT-qPCR检测［19］；行组织蛋白提取，制备蛋白样品进行Western blotting法检测；免疫组织化学法检测所需的空白涂片由江苏大学附属武进医院病理科提供，上述实验具体实验方法参考文献［20-21］。

1.12 统计学分析采用SPSS 22.0统计软件进行统计学分析。结肠癌组织中NKD1和AXIN2表达水平多组间比较采用单因素方差分析，组间两两比较采用t检验。采用单因素和多因素Cox回归分析方法评估AXIN2基因的表达与结肠癌患者总生存（overall survival，OS）的关系。以P＜0.05为差异有统计学意义。

2 结 果

2.1 TCGA和GEO数据库基因分别从TCGA和GEO数据库中下载的结肠癌数据进行基因表达差异分析，以∣logFC∣＞1和矫正后的P＜0.05作为筛选条件进行筛选，获得TCGA数据中的差异表达基因和GSE37182芯片中的基因表达差异图，其中红点代表结肠癌组织中高表达基因，绿点代表低表达基因，而黑点代表无明显表达差异的基因。见图1。

图1 TCGA和GEO数据库中的结肠癌数据差异表达分析Fig.1 Differential expression analysis on colon cancer data in TCGA and GEO databases

2.2 WGCNA分析采用WGCNA分析对TCGA和GEO数据库数据进行基因表达差异分析。分析TCGA数据获得的树状图及对应的颜色模块见图2 A，其中洋红和鲑鱼色颜色模块对应肿瘤中高表达的基因组合，而蓝色和黄色颜色模块对应正常结肠组织中高表达的基因组合（图2B）；GEO数据库分析获得的树状图和对应的颜色模块见图2 C，其中蓝色和粉色颜色模块对应肿瘤组织中高表达的基因组合，而绿松石色、红色、黄色及棕色颜色模块对应正常结肠组织中高表达的基因组合（图2D）。

图2 WGCNA法分析TCGA和GEO数据库中的结肠癌数据Fig.2 Colon cancer data in TCGA and GEO databases analyzed by WGCNA method

2.3 维恩分析和核心基因识别对WGCNA分析得到的两组高表达模块基因进行维恩分析，共获得13个共有基因（图3A）。通过在线STRING数据库对这13个基因进行PPI网络构建，检测这些蛋白之间的相互结合关系，发现ARID3A、SALL 4、NKD1、KND2、AXIN 2、CA 9和TACSTD2 7个基因之间存在相互关系，认为这7个基因为结肠癌潜在的核心基因（图3 B）。

图3 维恩分析和PPI网络Fig.3 Venn analysis and PPI network

2.4 核心基因的GO和KEGG富集分析对上述7个核心基因进行基因功能GO富集分析显示：在生物过程（biological process，BP）方面，NKD2/AXIN2/NKD1对经典/非经典Wnt信号通路进行正/负调控［1］；细胞成分（cellular component，CC）方面，基底外侧质膜（NKD2/CA 9/TACSTD2），质膜外侧（NKD2/TACSTD2）；分子功能（molecular function，MF）方面，NKD2/AXIN2主要结合泛素蛋白连接酶或结合泛素样蛋白连接酶。此外，对这7个核心基因进行KEGG富集分析显示：NKD2/AXIN2/NKD1主要集中在Hippo信号通路和Wnt信号通路。见图4。

图4 核心基因的GO和KEGG富集分析Fig.4 GO and KEGG enrichmental analysis on core genes

2.5 泛癌基因分析对上述7个核心基因进行泛癌（32种不同肿瘤组织）分析，结果发现TASCTD2基因在泛癌中的表达量高于其他6个核心基因（图5）。热图分析显示：7个核心基因在结肠癌中的表达相对于其他癌种均明显升高，其中NKD1和AXIN2在17个不同肿瘤组织中的表达量相似（图6）。在线UALCAN数据库分析显示：AXIN2基因在结肠癌组织中的表达量明显高于正常结肠组织中的表达，而AXIN 2在肿瘤组织中的表达量明显低于正常组织（图7）。NKD1和AXIN2蛋白在32种不同肿瘤组织中的表达相关性最高（相关性系数为0.69），此外，NKD2还分别与NKD1、AXIN2和ARID3A具有较高的相关性，相关系数分别为0.30、0.38和0.32（图8）。为分析核心基因的表达与免疫细胞亚型之间的关系，将肿瘤样品中的免疫细胞分为6个亚型：C1为伤口愈合型，C2为IFN-γ主导型，C3为炎性型，C4为淋巴细胞缺失型，C5为免疫静默型，C6为TGF-beta主导型等六个亚型。其中NKD2基因在C1、C2和C6亚型中高表达，而在C3、C4和C5亚型中的表达相对较低；NKD1基因在C5亚型中高表达，而在其他亚型中呈低表达；AXIN2基因在C2和C4中呈低表达，而其他亚型中呈高表达（图9）。核心基因的表达与肿瘤微环境的相关性分析结果显示：NKD2和NKD1基因的表达在泛癌中与综合（样品中基质细胞和免疫细胞总和）评分大多呈正相关关系。NKD2、AXIN 2、ARID3A和NKD1等基因的表达与结肠腺癌组织的综合评分呈负相关关系，即NKD2等基因的表达越高，结肠癌组织中肿瘤细胞数越多（图10）。

图5 泛癌分析法分析7个核心基因的相对表达量Fig.5 Relative expression amounts of 7 core genes detected by pan-cancer analysis method

图6 热图分析核心基因在不同癌组织中的表达量Fig.6 Heat map analysis on expression amount of core genes in different kinds of cancer tissues

图7 结直肠癌组织和癌旁组织中AXIN2基因的表达量Fig.7 Expression amounts of AXIN 2 gene in cancer and adjacent cancer tissues

图8 7个核心基因的相关性分析Fig.8 Correlation analysis among 7 core genes

图9 核心基因表达与免疫细胞亚型的相关性分析Fig.9 Correlation analysis between expression of core genes and immune cell subtypes

图10 核心基因的表达与肿瘤微环境的相关性分析Fig.10 Correlation analysis between expression of core genes and microenvironment of tumor

2.6 NKD1与AXIN2在癌组织和癌旁组织中的表达水平为了验证生物信息学分析结果，采用临床样本分别检测NKD1和AXIN 2在结肠癌组织和癌旁组织中的表达情况。RT-qPCR检测结果显示：AXIN2（图11）和NKD1（图12）在结肠癌组织中明显高表达，并且二者在结肠肿瘤组织中的表达水平具有显著正相关关系（图13）。通过免疫组织化学法检测进一步确认AXIN2、NKD1和NKD2在结肠癌组织中高表达，且表达水平具有正相关关系（图14）。对癌和癌旁组织进行Western blotting法检测发现：AXIN2和NKD1在结肠癌组织中明显高表达，且表达水平具有显著的正相关关系（图15）。

图11 在结肠癌和癌旁组织中AXIN2基因表达水平Fig.11 Expression levels of AXIN2 gene in colon cancer and adjacent tissues

图12 结肠癌和癌旁组织中NKD1基因表达水平Fig.12 Expression levels of NKD1 gene in colon cancer and adjacent tissues

图13 结肠癌组织中NKD1和AXIN2表达水平相关关系Fig.13 Correlation between expression levels of NKD1 and AXIN2 in colon cancer tissue

图14 结肠癌组织中AXIN2（A）、NKD2（B）和NKD1（C）表达情况（免疫组织化学,×20）Fig.14 Expressions of AXIN 2（A）,NKD2（B）,and NKD1（C）in colon cancer tissue(Immunohistochemistry,×20)

图15 Western blotting法检测结肠癌和癌旁组织中NKD1、NKD2和AXIN2表达电泳图Fig.15 Electrophoregram of expressions of NKD 1,NKD2，and AXIN2 in colon cancer and adjacent cancer tissues detected by Western blotting method

2.7 AXIN2基因的单因素和多因素预后风险分析对7个核心基因分别进行单因素和多因素预后风险分析结果显示：单因素预后风险分析和多因素预后风险分析均显示：只有AXIN2表达差异有统计学意义（P=0.026），并且AXIN2表达量（0.82）＜1，故AXIN2基因的表达可以作为结肠癌患者OS的低风险预后因子。此外，年龄（P=0.001）和原发肿瘤（T）（P=0.045）具有显著意义，并且二者表达量均大于1，故可以作为结肠癌患者OS的高风险独立预后因子。见表1。

表1 AXIN 2基因表达与结肠癌患者OS的单因素和多因素风险分析Tab.1 Univariate and multivariate risk analysis on expression of AXIN 2 gene and OS of colon cancer patients

3 讨 论

随着近两年生物信息学的发展，有较多潜在重要的结肠癌基因已被检测发现，但应用于临床的较少，其中存在的问题主要有：利用单一数据库存在系统误差、单一的差异分析方法导致假阳性等误差存在［25］以及对生物信息学分析结果进行实验验证。为了寻找结肠癌特异性的肿瘤标志物或潜在的重要的药物靶点，同时为了尽可能减少系统误差或假阳性等结果，本研究对2个不同的独立数据库（TGCA和GEO数据库）中结肠癌数据进行分析，避免单一数据库的系统误差出现。本研究对结肠癌基因表达数据采用2种不同的数据差异分析方法进行分析：基因表达差异分析方法和WGCNA法，通过2种不同的差异分析方法分析会较大程度上减少假阳性等误差的出现。本研究随后对4组差异分析数据取交集进行维恩分析，获得了13个共有基因。本研究通过在线STRING数据库对这13个基因进行蛋白网络构建并识别了具有相互关系的7个潜 在 的 核 心 基 因（ARID3A［26］、SALL 4［27］、NKD1［21，28］、NKD2、AXIN2［29］、CA 9［30］和TACSTD2［31］），这些结果与上述文献报道一致。

本研究对7个核心基因进行泛癌（32种不同肿瘤样本）分析，首先对7个核心基因在32个不同肿瘤组织中的相对表达量进行分析，发现TACSTD5基因在泛癌中的表达量高于其他6个核心基因。进行热图分析时，将正常组织样本量少于5个样本的肿瘤样本进行删除，使统计的数据更有意义，因此热图结果显示了17个不同的肿瘤组织。热图结果显示：7个核心基因在结直肠癌中的表达量均明显升高，此外，TACSTD2还在ESCA、THCA、BLCA和UCEC中高表达，CA 9基因还在BLCA、UCEC、HNSC、LUSC及GBM中高表达。整体而言，NKD1和AXIN 2基因在所有的肿瘤中的表达量均存在较高的相似度。本研究通过在线UALCAN数据库分析AXIN2基因的泛癌表达情况发现：AXIN2在结直肠癌组织中表达量明显高于癌旁的正常结肠组织，而在其他肿瘤样本中情况相反，即AXIN2在正常组织中的表达量明显高于肿瘤组织，该结果与本研究通过R语言分析结果一致。本研究对7个核心基因进行相关性分析显示：NKD1和AXIN2基因在不同肿瘤组织中存在较高的相关性（相关系数为0.69），此外，NKD2还与NKD1、AXIN 2和ARID3A存在较高的相关性。核心基因与免疫细胞亚型之间的相关性分析结果显示：TACSTD2基因在免疫细胞C1、C2、C3和C6亚型中表达量明显高于其在C4和C5亚型中的表达量。NKD1基因在免疫细胞C5亚型中的表达量明显高于其在C1、C2、C3、C4和C6亚型中的表达量。虽然本研究探讨了核心基因与免疫细胞亚型之间的相关性，但是由于篇幅有限，未对所有基因进行一一分析说明，仅代表性地以NKD1和TACSTD2为例进行说明。7个核心基因与肿瘤微环境的相关性分析结果显示：TACSTD2基因的表达与泛癌多数样品中综合（基质细胞和免疫细胞总和）评分呈显著正相关关系，即TACSTD2基因的表达越高，GBM、KICH、PCPG和UVM等肿瘤样品中的基质细胞和免疫细胞的总数越多，则肿瘤细胞数就会越少。AXIN 2基因的表达与ACC、COAD、GBM、MESO、READ、SARC和TGCT等肿瘤样品综合评分呈负相关关系，即AXIN2表达越高，这些肿瘤样品中的肿瘤细胞总数越多。为了验证生物信息学筛选结果，本研究通过实验进一步验证了AXIN2、NKD1和NKD2在结肠癌组织中高表达，且其表达具有显著正相关关系。本研究未对所有基因进行逐一验证说明，这也是本研究不足之处。此外，本研究对这7个核心基因进行结肠癌生存风险概率分析，分别对其表达与结肠癌患者的临床参数，如年龄、性别、肿瘤分期、原发肿瘤、区域淋巴结和远处转移等因素进行单因素及多因素分析，结果发现仅AXIN2基因的表达与结肠癌患者OS存在显著相关性；多因素分析结果显示：AXIN2的表达量为0.82，P=0.026，因此AXIN2为结肠癌患者OS的低风险预后因子。研究［32］显示：AXN2和β-catenin在结直肠癌组织中较高表达，且β-catenin的表达可以作为结直肠癌患者OS的独立预后因子，但是尚未提及AXIN 2与预后的相关性［32］。研究［33］显示AXIN2可能是结肠癌的一个风险标志物。但目前尚无文献报道AXIN 2可作为结肠癌OS的低风险独立预后因子，本研究是对该领域的重要补充。

本研究通过STRING数据库筛选了具有相互联系的7个核心基因，但这不代表另外6个基因无研究价值。本研究通过生物信息学数据分析发现这6个基因不能作为独立因素起作用。本研究通过生物信息学方法筛选潜在重要的核心基因及独立预后因子，这些结果需要更多数据进行进一步确认，为后续深入研究提供一定的数据支持。