骨肉瘤（osteosarcoma,OS）是最常见的原发性恶性骨肿瘤，由间充质细胞株发展而来，多发于四肢长骨干骺端，其特征为增殖的肿瘤细胞直接形成骨或类骨组织[1]。此外，OS能迅速破坏周围组织进行转移，具有较高的局部侵袭和转移倾向。近年来，随着手术切除结合全身化疗放疗治疗方案的应用，OS患者的远期生存率已经提高到65%～70%，但OS患者复发的概率依然很高，仍是患者临床死亡的主要原因之一[2]。OS基因组的复杂性和不稳定性及对OS发生的确切机制缺乏了解使OS治疗进展受到阻碍[3]。因此，亟需进一步研究OS的分子机制以探索OS的新治疗靶点。

生物信息学主要是处理遗传信息的学科，涉及计算工具和方法的研究、开发或应用，以获取、存储、可视化和解释医学或生物数据。近年来，高通量技术已视为医学肿瘤学具有巨大应用前景的工具，使许多基因表达数据得以整合，现已被广泛应用于探索OS的诊断和预后生物标志物。结合生物信息学的微阵列技术可以更全面分析OS发生和发展过程中的mRNA表达变化。有学者运用微阵列技术将早期生长反应1(early growth response1,EGR1）、CXC趋化因子配体10(C-X-C motif chemokine ligand 10,CXCL10）、骨髓细胞瘤癌基因（myelocytomatosis oncogene,MYC）和CXC趋化因子受体4(CXC chemokine receptor 4,CXCR4）确定为OS预后的潜在生物标志物，并进一步探索OS发展的分子调节机制[4]。通过微阵列技术了解OS的发病机制对于早期诊断和筛选治疗靶点具有重要意义。

本研究从基因表达综合（Gene Expression Omnibus,GEO）数据库中下载2个OS原始数据集GSE16088和GSE42572，通过矫正和分析后筛选出687个差异表达基因（differentially expressed genes,DEGs），其中包括523个上调基因和164个下调基因。将加权基因共表达网络分析方法（weighted gene co-expression network analysis,WGCNA）分析的最相关模块基因与DEGs取交集。对交集基因进行疾病本体论（Disease Ontology,DO）、基因本体论（Gene Ontology,GO）和京都基因与基因组百科全书分析（Kyoto Encyclopedia of Genes and Genomes,KEGG）富集分析，以及通过蛋白互作（protein-protein interaction,PPI）网络分析筛选出Hubbe基因并绘制受试者工作特征（receiver operating characteristic,ROC）曲线，使用另一组OS数据集GSE19276对Hubbe基因进行验证，筛选出OS关键基因。对关键基因进行浸润性免疫细胞相关性分析。通过实时荧光定量聚合酶链式反应（real-time quantitative polymerase chain reaction,qRT-PCR）检测关键基因在组织样本中的相对表达水平，从而进一步探索OS的潜在候选生物标志物，为寻找OS的药物治疗靶点和深入了解其发病机制奠定基础。

1、资料和方法

1.1 数据下载和数据处理

使用关键词osteosarcoma在GEO数据库中检索公开日期为2009-01-01/2022-01-01的数据，下载3个与本研究相符的OS原始数据集，即GSE16088、GSE42572和GSE19276。选择标准：（1）包含正常骨组织样本和OS组织样本；（2）单个芯片样本数量>10。GSE16088[5]包括9个正常骨组织样本和14个OS组织样本。GSE42572[6]包含5个正常骨组织样本[男性3例，女性2例；年龄11～50岁，（23.6±15.93）岁]和7个OS组织样本[男性3例，女性4例；年龄8～15岁，（12±2.58）岁]。GSE19276[7]包括5个正常骨组织样本[男性3例，女性2例；年龄49～79岁，（66±12.88）岁]和44个OS组织样本[男性23例，女性21例；年龄7～76岁，（23.64±17.15）岁]。将GSE16088和GSE42572作为训练集，GSE19276则为验证集。使用Perl软件注释每个数据集系列矩阵文件的探针，将微阵列探针的名称转换为基因名。对数据集进行log2转换和标准化。使用Perl软件合并GSE16088和GSE42572数据集，将每个样本中的基因纳入研究，获得1个整合的基因表达矩阵（14个正常骨组织样本和21个OS组织样本）。

1.2 DEGs筛选

使用R软件limma软件包在整合的基因表达矩阵中筛选OS组织与正常骨组织的DEGs。根据错误发现率（false discovery rate,FDR）调整P值（校正后P值），以校正后P<0.05和|log2FC|>1作为截取标准，使用pheatmap包绘制整合基因集中DEGs的可视化热图。

1.3 加权共表达基因网络分析

使用R软件中的WGCNA包和limma包对合并后的数据集进行WGCNA分析。WGCNA是一种非常有用的系统生物学方法，有助于构建自由尺度的共表达基因网络和检测基因模块。首先，定义一个相关性幂（软阈值参数）以确保标准的无标度网络，随之建立一个加权邻接矩阵显示基因之间的强相关性，去除弱相关性的基因。将邻接关系转化为拓扑重叠矩阵（topological overlap matrix,TOM）来度量基因的网络连通性，TOM将相邻基因的网络基因比率相加，并计算相应的相异性。基于TOM不相似性度量的平均连锁层次聚类法，对表现出与基因模块相似表达谱的基因进行聚集，构建聚类树状图。最后，选择r值绝对值最大的模块作为关键模块。此外，基于TOM及其聚类树状图的不同，将与OS最相关的模块基因可视化为热图。使用R软件中的VennDiagram包筛选出DEGs和WGCNA关键模块的交集基因用作后续分析。

1.4 交集基因的DO富集分析

DO是一个社区驱动的开源本体，旨在通过疾病概念链接不同的数据集，将高通量数据中的分子发现转化为临床相关性，进而从疾病的角度对基因进行注释。使用R中的DOSE和clusterProfiler软件包对交集基因进行DO富集分析，以了解交集基因与疾病的相关性。

1.5 交集基因GO功能与KEGG途径富集分析

GO富集分析是注释基因与基因产物和识别高通量基因组或转录组数据的特征生物学属性的常用方法，可根据定义的特征提供有关基因组产品基因功能的全面信息。GO富集分析包括生物学过程（biological process）、细胞成分（cellular component）和分子功能（molecular function）分析。KEGG是一个用于系统分析基因功能和基因组信息的数据库，目的是将基因组中的基因与通路中的基因产物加以联系。本研究通过org.hg.eg.db R包获得每个DEGs的Entrez-ID后，对DEGs进行GO功能富集和KEGG途径富集分析，以探索DEGs的功能。将P<0.05作为截取标准并认为差异具有统计学意义。

1.6 PPI网络构建和分析

PPI网络是了解细胞功能、疾病机制和药物设计的重要工具。通过STRING在线工具，选择交集基因的PPI得分（中位数置信度）>0.9为标准构建PPI网络，通过Cytoscape软件（https://cytoscape.org/，3.9.0版本）将PPI网络中的节点蛋白可视化，并隐藏断开的节点。最后将该网络中每个节点Degree值排名前10的蛋白认为是该网络中的关键蛋白。

1.7 ROC曲线

为了研究基因表达对OS疾病状态的影响，使用GraphPad Prism软件对10个OS的Hubbe基因绘制ROC曲线，并采用另一个数据集GSE19276进行验证，根据曲线下面积（area under the curve,AUC)>0.8且P<0.05的标准，筛选出具有较高的临床诊断价值的OS相关基因。其可作为OS的潜在生物标志物。

1.8 细胞免疫浸润分析

CIBERSORT算法是一种从基因表达谱中定量复杂组织细胞分数的通用反卷积算法，利用基因表达数据估计混合细胞中成员细胞类型的浸润情况，可以计算22种类型的免疫细胞的分布。为了评估免疫浸润的情况，将基因表达矩阵数据上传到CIBERSORT，获得免疫细胞浸润矩阵。然后根据P<0.05的标准筛选免疫细胞基质，用R软件包鉴定正常骨组织和OS组织样本中22种肿瘤浸润免疫细胞（tumor infiltrating immune cells,TIICs）在OS微环境中的含量，最后使用vioplot包绘制小提琴图可视化OS与正常组织样本免疫细胞浸润的差异。

1.9 关键基因与浸润免疫细胞的相关性分析

利用R软件中的limma、reshape2、ggpubr和ggExtra包等相关分析探讨关键基因与浸润免疫细胞水平的关联，使用ggplot2软件包绘制相关性棒棒糖图将其结果可视化。

1.10 qRT-PCR检测

为了证实生物信息学分析的结果，收集2020年9月1日至2022年6月30日期间于广西中医药大学第一附属医院住院手术的4例OS患者的OS组织及其癌旁组织用于qRT-PCR验证，其中男性3例，女性1例；年龄11～19岁，（15.75±3.5）岁；≥15岁2例，<15岁2例。本研究方案已通过广西中医药大学第一附属医院伦理委员会批准，所有患者均签署知情同意书。用Trizol提取组织样本中总RNA，通过HiFiScript cDNA第一条链合成试剂盒（江苏省泰州市康为世纪公司）将总RNA的RNA样本反转录成cDNA，使用CFX ConnectTM荧光定量PCR仪（Bio-Rad，美国）进行qRT-PCR。基因引物序列见表1。采用2-△△CT法计算各样本中目的基因相对mRNA表达。

1.11统计学分析

采用R.4.2.3和GraphPad Prism软件进行数据分析。计量资料采用均数±标准差表示。组间比较采用单因素方差分析。以P<0.05为差异具有统计学意义。

表1 引物序列

2、结果

2.1 DEGs的鉴定

根据矫正后P<0.05和—log2FC—>1的标准，从OS组织样本和正常骨组织样本中筛选出687个DEGs，其中包括523个上调的DEGs和164个下调的DEGs（图1）。

2.2 加权共表达基因网络分析

使用R中的WGCNA软件包分析得出10个模块基因，选择r值绝对值最大的黄色模块（MEyellow，|r|=0.65）作为与OS最相关的基因模块（图2A～2B）。该模块中共有2 338个基因。DEGs与WGCNA的交集基因共有545个（图2C）。

2.3 DO富集分析结果

DO富集分析结果表明，交集基因（DEGs与WGCNA最相关模块基因的交集）富集于泌尿系统癌症、肾癌、生殖细胞癌、肌肉骨骼系统癌症和胚胎癌等多种癌症疾病中（均P<0.05，图3）。

2.4 GO和KEGG富集分析结果

GO富集分析结果显示，交集基因（DEGs与WGCNA最相关模块基因的交集）在生物学过程中主要参与骨化的形成、细胞外基质组织、细胞外结构组织、外部封装结构组织和生物矿物组织发育；在细胞成分中，交集基因主要富集在含胶原蛋白的细胞外基质、内质网腔、核包膜、组蛋白脱乙酰酶复合物和胶原蛋白三聚体中；在分子功能中，交集基因在细胞外基质结构成分、整合素结合、赋予拉伸强度的细胞外基质结构成分、细胞外基质结合和血小板衍生生长因子结合中显著富集（图4A）。KEGG富集分析显示，交集基因主要富集在磷脂酰肌醇3激酶（phosphatidylinositide 3-kinase,PI3K）/蛋白激酶B(protein kinase B,PKB,又名Akt）信号通路、过氧化物酶体增殖物激活受体（peroxisome proliferators-activated receptor,PPAR）信号通路、黏着斑、蛋白质消化和吸收以及流体剪切应力和动脉粥样硬化等信号通路（图4B）。

2.5 PPI网络分析结果

用Cytoscape软件对上述分析所得到的的文本文件构建一个PPI网络，并使用Cytoscape软件中的Cytohubba插件对PPI网络中的重要节点基因进一步分析，获得排名前10位的Hubbe基因：磷脂酰肌醇4,5-二磷酸3-激酶催化亚基α(phosphatidylinositol-4,5-bisphosphate 3-kinase catalytic subunit alpha,PIK3CA）、脯氨酰4-羟化酶亚单位α1(prolyl 4-hydroxylase subunit alpha 1,P4HA1）、整合素αⅤ（integrin alphaⅤ,ITGAⅤ）、组蛋白去乙酰化酶2(histone deacetylase 2,HDAC2）、连环蛋白β1(catenin beta 1,CTNNB1）、Ⅲ型胶原蛋白α1(collagen typeⅢalpha 1,CO-L3A1）、Ⅰ型胶原蛋白α2(collagen typeⅠalpha 2,COL1A2）、转导蛋白β样1X相关蛋白1(transducin beta-like 1X-related protein 1,TBL1XR1）、小核核糖核蛋白多肽G(small nuclear ribonucleoprotein polypeptide G,SNRPG）和Ras相关核蛋白（Ras-related nuclear protein,RAN；图5）。

2.6 Hubbe基因的ROC曲线分析

对Hubbe基因（PIK3CA、P4HA1、ITGAV、H D A C 2、C T N N B 1、C O L 3 A 1、C O L 1 A 2、TBL1XR1、SNRPG和RAN）绘制ROC曲线以评估这10个基因的表达水平对OS的诊断效能。根据AUC≥0.8且P<0.05的标准筛选出7个关键基因（TBL1XR1、PIK3CA、COL3A1、COL1A2、CTNNB1、P4HA1和ITGAV，图6），并通过数据集GSE19276对这7个基因进行再次验证（图7）。结果显示，P4HA1和ITGAV的ROC曲线在训练集和验证集中均满足AUC≥0.8且P<0.05的标准，表明其可能是OS的潜在诊断标志物或治疗靶点。

2.7 免疫细胞浸润分析

免疫细胞浸润差异的小提琴图显示，与正常骨组织比较，OS组织中幼稚B细胞、CD8+T细胞和休眠期CD4+记忆性T细胞表达相对较低（均P<0.05）；而OS组织中γδT细胞、M0巨噬细胞和M2巨噬细胞的表达高于正常骨组织（均P<0.05，图8）。

2.8 关键基因与浸润性免疫细胞的相关性分析

ITGAV与M0巨噬细胞（r=0.74,P<0.01）、M2巨噬细胞（r=0.69,P<0.01）和γδT细胞（r=0.46,P=0.005）呈正相关，与CD8+T细胞（r=-0.75,P<0.01）、幼稚B细胞（r=-0.57,P<0.01）和休眠期CD4+记忆性T细胞（r=-0.41,P=0.015）呈负相关（图9A）。

P4HA1与M0巨噬细胞（r=0.62,P<0.01）和M2巨噬细胞（r=0.52,P=0.001）呈正相关，与CD8+T细胞（r=-0.61,P<0.01）和幼稚B细胞（r=-0.59,P<0.01）呈负相关（图9B）。

2.9 qRT-PCR验证

对P4HA1和ITGAV进行qRT-PCR实验验证显示，P4HA1和ITGAV mRNA表达量与微阵列分析结果一致，即P4HA1和ITGAV在OS组织中均呈高表达（均P<0.01，图10）。

图1 GSE16088和GSE42572数据集中DEGs的鉴定

Fig.1 Identification of DEGs in GSE16088 and GSE42572 datasets

图2 GSE16088和GSE42572数据集中与OS最相关基因模块的WGCNA分析

图3 交集基因（OS组织和正常骨组织的DEGs与WGCNA分析中与OS最相关模块基因的交集）的DO富集分析结果

3、讨论

OS是一种起源于原始转化细胞的侵袭性恶性骨肿瘤。据相关研究表明，手术治疗后OS患者复发率高达80%，严重影响患者的预后[8]。了解OS发生和发展的分子机制对于OS的诊断、治疗和预后至关重要。随着诸多学者对GEO数据库的深入挖掘，已有很多关于筛选OS生物标志物的研究。

本研究运用生物信息学技术和WGCNA对OS基因数据集分析发现，OS的DEGs和WGCNA的交集基因与肌肉骨骼系统癌症、肾癌和胚胎癌等肿瘤的发生和发展相关。GO分析表明，这些交集基因的生物学功能主要富集于细胞外基质组织。细胞外基质是细胞外微环境的主要组成部分，可以直接与细胞相互作用，通过整合素或其他细胞表面受体调节细胞生长、迁移、增殖、分化和代谢等功能[9]。有研究显示，肿瘤基质的细胞外基质分子失调会改变肿瘤细胞与白细胞和成纤维细胞的相互作用、促进血管生成和激活肿瘤细胞增殖与侵袭，从而促进肿瘤进展[10]。本研究结果与此基本一致，表明细胞外基质可能是OS发生和发展的关键驱动因素。

此外，OS的DEGs和WGCNA的交集基因主要富集在PI3K-Akt和PPAR信号通路。研究表明，调节OS细胞中PI3K-Akt信号通路的表达水平可抑制OS细胞的增殖、侵袭和迁移，从而抑制肿瘤的生长[11]。有学者报道三重基序蛋白46(tripartite motif 46,TRIM46）介导PPARα的表达水平可调节OS细胞的生长和凋亡，证实PPAR信号通路在抑制OS细胞增殖和诱导肿瘤细胞凋亡中发挥重要作用[12]。本研究通过KEGG分析发现，交集基因主要富集于PI3K-Akt信号通路，表明该信号通路可能与OS的发生和发展密切相关。

图4 交集基因（OS组织和正常骨组织的DEGs与WGCNA分析中与OS最相关模块基因的交集）的GO和KEGG富集分析

图5 交集基因（OS组织和正常骨组织的DEGs与WGCNA分析中与OS最相关模块基因的交集）的PPI网络分析结果和Hubbe基因的鉴定

注蓝色节点表示交集基因的普通基因，其他颜色节点表示Hubbe基因，其中颜色越深表示degree排名越靠前，灰色线条表示蛋白质之间的相互作用

图6 GSE16088和GSE42572数据集中PPI网络前10个Hubbe基因的表达水平对OS的诊断效能的ROC曲线分析

图7 GSE19276数据集中PPI网络Hubbe基因的表达水平对OS的诊断效能的ROC曲线验证

免疫浸润分析发现，OS组织中M0巨噬细胞和M2巨噬细胞表达较高，而CD8+T细胞和休眠期CD4+记忆性T细胞表达较低，说明巨噬细胞是OS组织中主要的浸润性免疫细胞，参与OS的增殖与进展。肿瘤相关巨噬细胞（tumor-associated macrophage,TAM）是肿瘤微环境中的一个主要影响成分，是参与炎性反应和组织稳态的免疫细胞[13]。研究表明，与原发性OS比较，在OS肺转移灶中发现更高密度的M2型巨噬细胞；M2型巨噬细胞增强肿瘤的侵袭力[14-15]。

使用OS数据集GSE19276对PPI网络排名前10位的Hubbe基因进行验证，筛选出P4HA1和ITGAV作为具有诊断价值的关键基因。为了进一步验证本研究结果，分别对4例OS患者的OS组织及其癌旁组织进行qRT-PCR验证，发现OS组织中P4HA1和ITGAV均呈高表达（均P<0.01），且与微阵列分析结果一致。P4HA1和ITGAV可能是OS的潜在生物标志物和治疗靶点。

图8 GSE16088和GSE42572数据集的正常骨组织与OS组织的免疫浸润分析

图9 ITGAV和P4HA1基因与浸润性免疫细胞的相关性分析

P4HA1是胶原蛋白合成过程中的关键酶，由2个相同的α亚单位和2个β亚单位组成，通过增加4-羟脯氨酸来增强胶原蛋白的修饰，调节细胞外基质代谢过程，参与肿瘤纤维化，增加肿瘤组织硬度，调节肿瘤免疫，促进肿瘤转移[16]。P4HA1上调可导致细胞外基质中基质金属蛋白酶（matrix metalloproteinases,MMPs）过表达，进而促进肿瘤细胞增殖、侵袭和迁移[17]。OS发生和发展的主要原因之一是OS细胞外基质的降解，而MMPs过表达可促进OS细胞外基质成分降解，从而导致OS细胞通过血液循环向远处发生转移[18-19]。同时，研究发现，通过阻断PI3K-Akt信号通路激活抑制OS细胞外基质的降解可能是OS的潜在治疗策略[20]。此外，也有研究表明，P4HA1过表达与OS的不良预后密切相关[21]。OS组织中P4HA1表达上调促进OS细胞的迁移和侵袭[22]。因此，推测P4HA1可能为OS的潜在生物标志物和治疗靶点并通过调控细胞外基质的功能对OS的发生和发展产生影响，然而其具体作用机制有待进一步深入研究。

图1 0 qRT-PCR实验检测OS组织和癌旁组织中P4HA1和ITGAV的mRNA表达水平

ITGAⅤ是整合素α链家族的成员，其编码产物是整合素的一种成分，可与5个整合素β亚基结合形成αvβ1、αvβ3、αvβ5、αvβ6和αvβ8整合素受体；其表达水平在诱导肿瘤新血管生成和促进肿瘤迁移与侵袭中发挥重要作用[23]。血管生成是指从现有血管系统形成新血管的过程。其在OS生长和转移过程中至关重要。血管内皮生长因子（vascular endothelial growth factor,VEGF）是血管生成的重要调节因子，可调节内皮细胞基因的激活形式，诱导内皮细胞基因表达组织的相关因子，并通过诱导间质胶原酶和蛋白水解酶刺激血管生成，促进OS细胞侵袭[24]。研究表明，ITGAV可通过介导VEGF诱导内皮细胞迁移和成管腔样结构发育，进而生成新血管，促进肿瘤细胞的转移[25]。此外，ITGAV的表达水平在抑制或促进肿瘤细胞的迁移和侵袭中发挥关键作用。相关研究证实，ITGAV在OS组织中表达量的增加促进OS细胞的增殖、侵袭与迁移[26]。本研究中，GEO数据库分析和实验验证显示，ITGAV在OS组织中的mRNA表达量高于癌旁组织，与国内外研究报道基本一致[26-28]。综上所述，ITGAV可能通过调控肿瘤内皮细胞血管生成参与OS的发生和发展。

总之，本研究通过生物信息学分析和qRT-PCR实验验证筛选出对OS具有诊断和治疗价值的关键基因P4HA1和ITGAV，同时也从生物信息学的角度验证细胞外基质、PI3K-Akt和PPAR信号通路等在OS发生和发展中的重要意义，有利于更加深入理解OS发展的分子调节机制，为后续OS的诊断和治疗提供新的研究方向。P4HA1和ITGAV可能是OS的新的生物标志物和治疗靶点。然而本研究仍存在一定的不足，单一的OS组织样本缺乏广泛的代表性，需要更多研究来阐明P4HA1和ITGAV在OS诊断和治疗中的作用。

参考文献:

[15]林跃玮,刘文渤,胡志锋,等.基于SEER数据库使用支持向量机和列线图对骨肉瘤诊断时肺转移风险预测模型的建立与验证[J].实用肿瘤杂志, 2023, 38(3):251-257.

[19]赖续文,罗锦霞,崔文志,等.经典型骨肉瘤中Hes1蛋白的表达及其临床意义[J].实用肿瘤杂志, 2022,37(4):315-319.

[25]吕凤,郭晓静,付丽.整合素表达变化在乳腺癌发生及发展中的作用研究进展[J].中华病理学杂志, 2015,44(4):294-296.

[27]范兆阳,尹彩星,鲜文峰,等.利用蛋白组学技术探讨骨肉瘤肺转移患儿血清中蛋白表达情况[J].中国卫生检验杂志, 2023, 33(3):261-265.

[28]刘伟,赵宝辉,刘晓峰. MiR-136通过下调靶基因ITGAV表达抑制骨肉瘤细胞增殖及迁移的分子机制[J].临床与病理杂志, 2022, 42(5):1013-1021.

基金资助:国家自然科学基金项目(81860793); 广西自然科学基金项目(2020JJA140375);

文章来源:李威材,秦刚,苏国威,等.基于生物信息学分析骨肉瘤的关键基因和qRT-PCR实验验证[J].实用肿瘤杂志,2024,39(04):344-355.