梁一豪，赖颖君，袁燕文，袁 炜，张锡波，张拔山，卢志锋
南方医科大学第十附属医院（东莞市人民医院）1检验科，2消化内科，3病理科，广东 东莞523059
2022年国家癌症中心发布2016年我国癌症数据，统计结果显示我国胃癌发病率为28.68/10万，死亡率为20.87/10 万，均位居恶性肿瘤发生和死亡人数的第3位［1］。截至2020年我国胃癌新发和死亡病例仍持续上升，分别占全球的43.94%和48.62%，均处于全球较高水平，其中晚期患者占新发病例的比例较高，亟需加强防控和改善流行状况［2-4］。因此，早期筛查和靶向治疗对于胃癌的早发现和诊断治疗至关重要。目前主流的胃癌筛查方法主要包括胃镜检查、X线钡餐造影、血清幽门螺旋杆菌抗体以及胃肿瘤相关标志物检测，如癌胚抗原（CEA）、细胞角蛋白19片段（CYFRA21-1）、胃肠道肿瘤相关抗原（CA72-4）、血清胃蛋白酶原和生长激素释放肽等［5,6］。确诊方法是通过内窥镜或手术获取的病理组织活检，临床常用的胃癌标志物可作为辅助诊断和判断疗效预后的依据，但其特异性和灵敏度仍有一定的局限性［7］。分子标志物是肿瘤发生发展过程中细胞内部或间质中异常表达的基因或蛋白质组分，如EGFR、p53、HER2和FGFR2等基因在胃癌中普遍存在过度表达和激活，靶向上述基因的药物治疗也逐步进入临床试验［8-10］。但由于胃癌具有细胞分化程度低、增殖和侵袭速度快、肿瘤微环境复杂等性质，潜在的致病机制如免疫浸润和炎症反应等过程尚未阐述清楚，因此在依赖靶基因进行早期诊断和药物治疗方面，胃癌相关分子标志物仍有广阔的探索前景［11］。
随着基因芯片技术和生物信息学分析方法的成熟，全基因组筛查和分析验证是发现新型肿瘤分子标志物的有效方法之一。通过整合分析GEO和TCGA数据，目前已发现GPX8、BIRC5、NPM1等靶基因在胃癌中的表达显著上调且与胃癌患者的生存预后显著相关，表明基因芯片技术和生物信息学分析方法是探究胃癌等疾病发生机制有效方法之一［12-14］。但由于胃癌的发病机制复杂，涉及多个生物学过程，如细胞凋亡、细胞自噬和炎症细胞浸润等，因此与胃癌相关的核心基因仍需进一步筛选和验证［15］。因此，本研究拟从GEO数据库下载胃癌相关的数据集，寻找差异表达基因（DEGs），并对DEGs进行生物信息学分析筛选得到核心基因。利用癌症基因图谱数据库中的胃腺癌（TCGA-STAD）基因组测序数据进行分析验证。最后利用qRT-PCR检测核心基因的表达水平，为胃癌的早期诊断、靶向治疗及预后分析挖掘潜在的有效靶点。
1.1 芯片数据的获取
根据样本来源为人胃癌与正常胃黏膜组织、每组样本例数≥10，以及芯片数据引用率高等条件，从NCBI的GEO 数据库（https://www.ncbi.nlm.nih.gov/geo/）中筛选下载3组胃癌相关的基因芯片数据集（GSE118916、GSE54129、GSE79973），保证芯片数据的代表性和可靠性。其中数据集GSE118916基于GPL15207平台，而数据集GSE54129、GSE79973 基 于GPL570 平 台。GSE118916 包含GC 组织和癌旁组织各15 例；
GSE79973 包含GC 组织和癌旁组织各10 例；
GSE54129包含111例GC组织和21例正常胃组织。
1.2 检验芯片数据的质量
利用R包FactoMineR和factoextra进行主成分分析（PCA），对各组基因芯片中所有的样本表达数据进行降维分析，在二维坐标平面内得出样本降维后的分布和离散情况，验证胃癌和对照组样本的分组合理性。根据基因的表达数据，对样本进行层次聚类，并且挑选方差值排行前1000个基因绘制分组聚类热图，观察基因在GC与对照组样本间的表达差异。
1.3 筛选GC相关的DEGs
利用GEO数据库的在线分析工具GEO2R处理和分析每组表达芯片数据，筛选阈值为差异倍数|log2FC|>2.0，校正后P值<0.05，筛选在GC中表达差异的基因，其中log2FC>2.0 的基因被标记为上调基因，log2FC<-2.0 的基因被标记为下调基因。利用R 包tinyarray可视化差异基因的表达情况，并通过绘制韦恩图筛选3组数据集中共有DEGs的作为共表达差异基因集。
1.4 GO功能和KEGG富集分析
为进一步探究共表达差异基因集的功能，利用R包ClusterProfiler对其进行基因本体（GO）功能注释和京都基因与基因组百科全书（KEGG）信号通路的富集分析，其中GO包括3部分：生物过程（BP）、分子功能（MF）和细胞组分（CC），GO和KEGG分析均以P<0.05为标准筛选结果。通过气泡图对显著富集的GO和KEGG条目进行可视化分析，同时结合共表达差异基因的log2FC值绘制圈图和弦图，进一步分析DEGs在GC发生发展中参与的关键生物过程和信号通路。
1.5 构建蛋白互作网络及筛选核心基因
为探索DEGs潜在的互作关系，利用String（http://string-db.org/）数据库构建DEGs蛋白互作网络（PPI），以互作评分>0.4作为阈值条件，将关联基因的互作数据导入到Cytoscape 软件进行可视化分析，通过Cytohubba的最大互相关方法（MCC）分析每个节点基因的网络特征并赋值，依据分值高低筛选PPI中重要的核心基因。
1.6 核心基因的验证分析
TCGA是一个开源性的肿瘤数据库，收录众多肿瘤和癌旁样本的基因组测序数据和临床随访信息，旨在为科研工作者提供一个广阔的探索平台。从TCGA中下载407例STAD和32例癌旁组织的RNA-seq数据和临床随访信息，基于Mann-Whitney U检验和配对样本t检验分析对比核心基因在STAD和癌旁组织中的表达水平；
通过绘制受试者工作曲线（ROC）评估核心基因的诊断效能。基于Kaplan-Meier plotter（http://kmplot.com/analysis/）数据库绘制生存曲线评估核心基因的预后价值。人类蛋白质图谱（THPA）数据库是涵盖肿瘤与正常组织蛋白组学信息的整合型数据库，收录了胃癌与正常胃粘膜组织的免疫组化芯片数据，用于验证核心基因的在组织上的表达水平。上述研究方法见流程图（图1）。
图1 生信分析流程图Fig.1 Flowchart of bioinformatics analysis.
1.7 细胞培养和qRT-PCR验证基因表达
人正常胃粘膜上皮细胞系GES-1、人胃癌细胞系MKN45、MKN74、AGS、HGC-27和KE-39在含10%的胎牛血清的RPMI 1640培养基中培养，置于37 ℃恒温和含5%CO2细胞培养箱中传代增值至细胞密度约为80%。利用Trizol在4 ℃低温下提取各总RNA，按照逆转录试剂盒（PrimeScripeTMRT Master Mix,Takara）说明配置体系，合成cDNA。根据荧光定量PCR试剂（TB Green®Premix Ex Taq™,Takara），选用GAPDH作为内参基因，进行基因的定量扩增和检测，其余差异核心基因引物见表1。以GES-1 中各差异核心基因的mRNA 水平为对照，计算其在不同胃癌细胞系中的相对表达含量。
表1 核心基因引物序列Tab.1 Primer sequence of the hub genes for RT-PCR
1.8 统计学方法
利用SPSS 22.0软件对实验数据进行统计分析，计量数据以均数±标准差表示。利用Shapiro-Wilk 和Levene法检验相关数据是否符合正态分布和满足方差齐性。当数据符合正态分布且方差齐时，多组比较采用单因素方差分析；
方差不齐时采用非参数检验。本研究中qRT-PCR数据符合正态分布，采用t检验比较各组胃癌细胞株与胃黏膜上皮细胞的差异，以P<0.05为差异有统计学意义，所有实验均独立重复3次。
2.1 数据质量控制
由3组胃癌芯片的原始数据，提取出表达矩阵、样本分组信息和芯片对应的平台注释ID（表2）。对每组基因芯片中所有样本的表达数据进行PCA，降维结果显示，在二维坐标轴内GSE54129和GSE118916数据集中的胃癌组和正常组织对照组样本没有产生交集，而GSE79973虽有部分区域重叠，但整体可区分（图2）。分组聚类热图显示，3组基因芯片表达数据在胃癌与正常对照组织中存在显著差异（图3）。
表2 胃癌表达芯片数据集信息Tab.2 GC-related microarray information in GEO dataset
图2 胃癌基因芯片数据的主成分分析Fig.2 Principal component analysis of GC microarray data.
图3 胃癌基因芯片数据的分组聚类热图Fig.3 Grouping clustering heatmap of GC microarray data.
2.2 胃癌相关的DEGs及共表达差异基因
通过GEO2R分析GSE118916、GSE54129、GSE79973 数据集得到表达上调和下调的差异基因。绘制火山图显示GSE54129共筛选出456个DEGs，包括上调基因202个，下调基因254个；
GSE79973共筛选出379个DEGs，包括上调基因81个，下调基因298个；
GSE118916 共筛选出231 个DEGs，包括上调基因56个，下调基因175个（图4）。绘制韦恩图获取77个共同的DEGs（图5）。
图4 差异表达基因的火山图Fig.4 Volcano map of the DEGs.Blue color represents the downregulated genes(log2FC<-2.0,P<0.05)and red the upregulated genes(log2FC>2.0,P<0.05)in GC tissues.
图5 共表达差异基因的韦恩图Fig.5 Venn diagram of co-expressed differential genes.
2.3 共有DEGs涉及的生物学过程与信号通路
GO功能富集分析结果显示，77个共有DEGs主要位于胶原纤维丰富的ECM，在细胞的基底膜和顶端分泌部位均有分布。在调节生物过程方面主要参与机体消化和激素代谢过程，且在应对外源性刺激反应和维持组织稳态中发挥重要作用。在分子功能方面主要具有作用于CH-OH基团的氧化还原酶活性，同时具备ECM受体和配体的活性。KEGG分析结果显示，77个共有DEGs主要在视黄醇代谢通路、胃酸分泌途径、ECM受体互作途径、蛋白质的消化和吸收等信号通路富集（图6，表3）。
表3 GC相关DEGs的GO和KEGG富集分析Tab.3 GO and KEGG enrichment analysis of GC-related DEGs
图6 GO和KEGG分析显示DEGs显著富集的条目Fig.6 GO and KEGG enrichment analysis of the DEGs.Y axis:GO terms and KEGG pathways;X axis:GeneRatio (The ratio of the number of target genes enriched to the target pathway genes to the total number of target genes,including the gene set and total genes).
结合77个DEGs的差异表达倍数，对其中显著富集的条目进行圈图和弦图的可视化，结果显示机体消化途径、作用于供体的CH-OH基团的氧化还原酶活性、细胞激素的代谢过程、对异生刺激的反应、组织稳态、视黄醇代谢通路以及胃酸分泌等生理过程和信号通路在GC中的调节作用显著下降，涉及的大多数基因均表达下调；
而在ECM结构成分、ECM-受体的相互作用通路、蛋白质代谢等肿瘤相关部位和信号通路在GC中的调节作用显著上升，其所涉及的大多数基因均表达上调（图7）。
图7 圈图和弦图可视化DEGs的表达水平和GO和KEGG的富集结果Fig.7 Circle and chord diagram visualization of DEGs expression and enrichment results.A:Enrichment results between DEGs and entries significantly enriched for GO and KEGG. B: Relationship between DEGs with differential expression fold DEGs and significantly enriched entries for GO and KEGG.P value:Height of the inner circle fan diagram;LogFC value:Height of the gene in the exosphere.Z-score=(Up-regulated genes-downregulated genes)/total genes.
2.4 PPI的构建及核心基因的获取
利用String 数据库对77 个共有DEGs 构建PPI（图8），将关联基因的互作数据导入Cytoscape软件进行可视化（图9），分析相互作用关系，发现多个网络关键节点。利用Cytohubba插件的MCC算法对PPI中的各个基因节点进行赋值排名，共获得9个核心基因，分别是：SPARC、TIMP1、THBS2、COL6A3、TFF1、GKN1、THY1、TFF2、PGC（图10）。
图8 Sting数据库中共有差异基因的蛋白互作网络图Fig.8 PPI of common differential genes in String database.
图9 Cytoscape可视化PPI网络图Fig.9 Cytoscape visualization of the PPI network.
图10 Cytohubba分析PPI中的核心基因节点Fig.10 Cytohubba analysis of the core gene nodes in PPI network.
2.5 核心基因的预后验证分析
TCGA-STAD数据验证核心基因的表达水平结果显示，SPARC、TIMP1、THBS2、COL6A3和THY1在GC中的表达显著上调，而TFF1、TFF2、GKN1和PGC在GC 中表达下调；
在配对GC 样本中SPARC、TIMP1、THBS2、COL6A3、THY1和GKN1可得到相同的结果，而TFF1、TFF2和PGC差异无统计学意义（图11）。生存曲线结果显示，SPARC、TIMP1、THBS2、COL6A3、THY1高表达的GC患者生存时间下降（P<0.05），TFF1、TFF2低表达与胃癌患者总体生存时间有关（P<0.05），其他基因的异常表达对GC 患者生存时间无明显差异（图12）。AUC结果显示，TIMP1、THY1、SPARC、THBS2对胃癌有较高的诊断价值（P<0.05，图13）。免疫组化结果显示，SPARC、TIMP1、THBS2、COL6A3在胃癌患者肿瘤组织中的表达上调（图14）。分析SPARC、TIMP1、THBS2、COL6A3在各个病理分型样本中的表达水平，发现SPARC和THBS2在I期和II期GC样本中差异具有统计学意义；
TIMP1在I期和II期、III期、IV期样本的比较中的差异具有统计学意义；
COL6A3则在I期与II期、III期中具有显著差异（P<0.05，图15）。
图11 核心基因在TCGA数据库中的配对和非配对GC样本的表达水平Fig.11 Expression levels of the hub genes in paired and unpaired GC samples in the TCGA database analyzed using paired sample T test(A)and Mann-Whitney U test(B).*P<0.05;**P<0.01;***P<0.001.
图12 GC核心基因的预后分析Fig.12 Prognostic analysis of the hub genes in GC.
图13 核心基因的ROC曲线Fig.13 ROC curves of the hub genes in GC.
图14 THPA数据库验证核心基因在GC组织中的表达水平Fig.14 Verification of expressions of the hub genes in GC tissues by immunohistochemical staining from THPAdatabase.
图15 TCGA-STAD病理分型数据验证核心基因在各期中的表达水平Fig.15 Verification of hub genes in different pathological stages from TCGA-STAD.*P<0.05;**P<0.01;***P<0.001.
2.6 qRT-PCR验证
qRT-PCR结果显示，相对于正常人胃粘膜上皮细胞，SPARC在AGS、HGC-27和KE-39细胞株中表达上调（P<0.05）；
THBS2在HGC-27和KE-39细胞株中表达上调（P<0.05）；
TIMP1在5种候选胃癌细胞株中均表达上调（P<0.05）；
COL6A3在MKN45细胞株中表达显著上调（n=3，P<0.05，图16）。
图16 核心基因SPARC、THBS2、TIMP1和COL6A3在人胃粘膜上皮细胞和胃癌细胞系中的差异表达Fig.16 Differential expressions of the hub genes SPARC,THBS2,TIMP1 and COL6A3 in normal human gastric epithelial cells and gastric cancer cell lines.*P<0.05 vs GES-1 group.
早期诊断和手术切除是防治胃癌的有效方式，然而对于早期无症状患者，胃内窥镜属于侵入性检查，人群检出率较低；
对于晚期患者，受身体素质等限制，难以获得合适的手术条件，且术后转移和复发也是不可忽略的危险因素［16］。分子标志物是可以通过分子生物学技术检测的肿瘤相关基因、蛋白质、代谢产物等生物大分子，可用于胃癌的早期诊断、预后评估和靶向治疗［17］。因此，利用生物信息学方法挖掘胃癌相关基因芯片的数据，筛选疾病进展相关的核心基因，是实现早期诊断和靶向治疗的关键。
本研究通过分析来自GEO数据库的3个胃癌相关的基因芯片数据，发现77个共表达差异基因，可能通过不同的生物学途径或信号通路影响胃癌的发生和进展。对其进行GO和KEGG富集分析，结果显示差异基因主要表达在ECM 和基底膜部位，具有作用于带有CH-OH基团的氧化还原酶和ECM受体配体的活性，参与细胞激素代谢、视黄醇代谢、胃酸分泌和ECM受体互作等胃恶性肿瘤相关的信号通路。研究表明ECM是正常细胞维持组织稳态的重要结构，在胃癌组织中异常表达的胶原和沉积蛋白容易破坏ECM结构促进肿瘤微环境形成，在调控癌症的侵袭和转移过程中发挥关键作用［18,19］。胶原蛋白是细胞外基质的主要成分，研究发现与正常组织相比，胃癌间质中胶原蛋白的沉积显著增加，而且胶原蛋白的含量、成熟度、形态和结构等在胃癌中改变明显，可作为胃癌潜在的辨别指标［20］。核心基因SPARC、TIMP1、THBS2、COL6A3 和THY1 等在肿瘤细胞外基质中均有分布和表达，此外本研究还发现在维持组织稳态的基因TFF1、TFF2 和GKN1 表达显著下调，表明胃癌破坏正常胃粘膜的微环境。激素代谢异常是引起相关癌症发生的危险因素，如胃泌素代谢异常是幽门螺旋杆菌感染介导胃癌发生的原因之一［21］。视黄醇是人体正常胃粘膜组织中高表达的活性分子，发挥胃粘膜的免疫调节功能，在视黄醇代谢的富集条目中，靶基因表达显著下调，表明筛选的共表达差异基因可能通过视黄醇代谢途径调节胃粘膜的免疫功能［22］。因此，本研究的差异基因富集分析结果与胃癌密切相关，在胃癌的组织稳态调节通路方面仍具有探索空间。
进一步通过STRING数据库对共表达差异基因进行蛋白质网络互作分析，通过cytohubba计算前4个得分最高的关键基因，分别是SPARC、TIMP1、THBS2和COL6A3，提示其在GC 调控中发挥关键作用。早在2010年SPARC被发现与胃癌进展及患者的不良预后息息相关［23］，其作为分泌蛋白在胃癌病理组织中表达显著升高，经多因素分析表明SPARC是评估肿瘤侵袭深度和淋巴结转移的独立危险因素［24］。此外SPARC联合GLT8D2和VCAN建立关于肿瘤成纤维细胞的风险预测模型，提示高风险组患者预后不良且耐受免疫治疗［25］，佐证本研究发现SPARC在ECM相关的条目中显著富集，在胃癌中表达显著升高。TIMP1具有抑制基质金属蛋白酶的活性，可以破坏基底膜和ECM成分，是促进GC 细胞侵袭和转移的关键［26］。研究发现外周血TIMP1水平与肿瘤病理分期相关，且高表达的TIMP1水平提示不良预后，是与胃癌进展相关的理论依据［27］。THBS2属于内皮间质细胞转换的相关基因，研究发现其与胃癌的微卫星不稳定性亚型相关，可用于预测总体生存率和区分不同模式的胃癌免疫微环境［28］。本研究也发现TIMP1和THBS2在细胞外基质与细胞间的相互作用中显著富集，与胃癌的总体生存率显著相关。此外有研究发现ECM相关蛋白THBS2联合胶原蛋白家族COL6A3，二者的高表达与肿瘤的病理T分期相关，同时与患者的不良预后及肿瘤的免疫浸润程度相关［29］。此外胶原蛋白家族的其他成员也通过基因芯片数据鉴定为胃癌相关的核心基因，表明本研究挑选的COL6A3靶标有一定研究依据［30］。随着基因芯片和测序技术的发展，对胃癌靶基因的筛查逐步细化到肿瘤细胞黏附、侵袭和免疫浸润等具体进程，因此本研究仍需关注SPARC、TIMP1、THBS2、COL6A3等核心基因在具体通路中的作用［31,32］。
综上所述，本研究应用生物信息学分析方法筛选出可能与胃癌病程调控相关的4 个核心基因SPARC、TIMP1、THBS2和COL6A3。经表达水平和诊断预后的验证分析，表明其可能在GC的进展中发挥关键作用，为胃癌的早期诊断、靶向治疗及预后分析挖掘潜在有效的靶点提供一定数据参考。但本研究仍存在不完善之处：核心基因对胃癌的作用仍需要进一步的体内外实验验证；
另一方面，核心基因的具体调控和作用机制是后续研究的关注重点。
