基于高通量数据的骨肉瘤相关核心基因的挖掘和鉴定分析

Mining and Identification of Osteosarcoma Related Hub Genes Based on Microarray Data

摘要: 背景：骨肉瘤是源于间充质细胞的恶性肿瘤，起病隐匿，好发于青少年，早期症状仅表现为关节肿胀和疼痛，且进展迅速，出现远处转移时其生存率大幅降低。为了筛选骨肉瘤发生发展的潜在基因，本研究从GEO数据库中获得GSE99671和GSE197158进行生物信息学分析。方法：首先，利用R包进行去批次处理，再利用R包“limma”识别差异表达基因(DEGs)，通过GO (Gene Ontology)和KEGG (Kyoto Encyclopedia of Genes and Genomes)分析对差异表达基因(DEGs)进行功能注释。利用STRING工具构建蛋白–蛋白相互作用(PPI)网络，挖掘出最重要的模块和核心基因。结果：共鉴定出275个差异表达基因。差异表达基因的定位主要集中在细胞外基质和细胞外骨架中，主要定位于细胞外基质胶原纤维中和分泌小泡内，功能变化主要集中在细胞识别，细胞膜、骨架等结构形成和多种酶或转运体如糖胺聚糖结合和无机阳离子跨膜转运体活性等方面。CTSG、MPO、ELANE、DEFA4、CAMP、LTF、PRTN3、MMP8、LCN2、BPI基因被鉴定为核心基因。结论：综上所述，本研究中发现的差异表达基因和核心基因有可能成为潜在的诊断和治疗靶点。

Abstract: Background: Osteosarcoma is a malignant tumor derived from mesenchymal cells. Its onset is oc-cult. It often occurs in adolescents. Its early symptoms are only joint swelling and pain, and its pro-gress is rapid. Its survival rate is greatly reduced when distant metastasis occurs. In order to screen the potential genes for the occurrence and development of osteosarcoma, GSE99671 and GSE197158 were obtained from geo database for bioinformatics analysis. Methods: Firstly, the R package was used for de batch processing, and then the R package “limma” was used to identify dif-ferentially expressed genes (DEGs). The function of the differentially expressed genes (DEGs) was annotated by go (gene ontology) and KEGG (Kyoto Encyclopedia of genes and genes) analysis. Using string tool to construct protein protein interaction (PPI) network, the most important modules and core genes were mined. Results: 275 differentially expressed genes were identified. The location of differentially expressed genes is mainly concentrated in the extracellular matrix and extracellular skeleton, mainly in the extracellular matrix collagen fibers and secretory vesicles. The functional changes are mainly concentrated in cell recognition, the formation of cell membrane, skeleton and other structures, the binding of a variety of enzymes or transporters such as carbohydrates and the activity of inorganic cationic transmembrane transporters. CTSG, MPO, ELANE, DEFA4, CAMP, LTF, PRTN3, MMP8, LCN2, BPI were identified as hub genes. Conclusion: In conclusion, the differentially expressed genes and core genes found in this study may become potential diagnostic and therapeu-tic targets.

1. 引言

骨肉瘤好发于长骨，偶见发于软组织，是临床上青少年最常见的恶性肿瘤之一，其特征为肿瘤组织内部大量的类骨质和未成熟骨 [1] [2]。尽管目前手术治疗、化疗、靶向治疗等多种手段被应用于骨肉瘤的治疗当中 [3] [4]，但由于患者临床表现不典型，发现分期较晚以及肿瘤恶性程度高易出现侵袭转移等原因，虽然治疗措施在近几年有很多变化，但是生存率并没有明显改变 [5]。因此深入研究骨肉瘤的分子机制尤为重要。

现代生物信息学是一门新兴交叉学科，主要利用计算机技术协助下一代测序数据分析 [6]，随着数据量的指数爆炸性增长，我们的研究手段也日新月异。数据来源主要有基因组学和蛋白质组学，常见的分析方法有GO分析、GSEA分析、WCGNA分析等，未来生物信息学则向全基因组关联研究(Genome-wide association study, GWAS) [7]，以及单细胞测序分析 [8] 等方向发展，未来前景广阔，而我们的研究则使用生物信息学的研究方法寻找差异表达基因，指导骨肉瘤的后续诊疗。

ELANE基因属于丝氨酸蛋白酶家族，编码具有细胞毒性的丝氨酸蛋白酶即中性粒细胞弹性蛋白酶(NE)。有研究发现，肺癌中血清中NE升高不仅与疾病状态呈正相关，而且与疾病进展呈正相关 [9] [10]。CTSG基因属于肽酶S1家族成员之一，编码组织蛋白酶G，参与先天免疫、炎症、感染、凝血和细胞凋亡等的调控 [11]。有研究表明，CTSG可促进口腔鳞癌发生发展 [12]，还对急性髓性白血病的免疫逃避有一定意义 [13]。然而目前这两种基因与肺癌之间的关系尚不清楚。

所以本研究拟以R语言为研究基础，通过不同功能的R包的对来自数据库的高通量测序数据做处理，去除不同样品的批次效应，寻找差异表达基因并行GO和KEGG富集分析，找到差异表达基因的定位及功能，并最终通过PPI网络找到差异表达基因中的核心基因，旨在指导骨肉瘤后续的临床诊疗发展。

2. 材料和方法

2.1. 数据集

我们从gene Expression Omnibus (GEO)数据库下载了两个数据集，GSE136043和GSE146460 (https://www.ncbi.nlm.nih.gov/geo/)。GSE99671包括18个骨肉瘤样本和18个正常骨组织样本，GSE197158包括11个骨肉瘤样本和3个正常骨组织样本。

2.2. 差异表达基因(DEGs)标识

我们应用了基于R语言开发的R包“BiocManager”，用于消除不同数据集间的批次效应，合并为一个数据集。之后我们应用R包“limma”识别合并后数据集的骨肉瘤和正常组别之间的差异表达基因。临界标准为P值 < 0.05，log(FC) > 1或log(FC) < −1。

2.3. 差异表达基因注释

Metascape (http://metascape.org/gp/index.html)网站和基于R语言开发的R包“org.Hs.eg.db”，“clusterProfiler”，“ggplot2”，“enrichplot”是对基因进行生物过程(BP)、细胞成分(CC)、分子功能(MF)和KEGG (Kyoto Encyclopedia of Genes and Genomes)分析等的强有力注释工具。我们通过上述R包和Metascape网站注释了共同差异转录基因的功能。

2.4. 蛋白质相互作用网络的构建

利用网站STRING (http://string-db.org)，我们可以检索相互作用基因，并形成相应的蛋白–蛋白相互作用(PPI)网络。该网站可将经limma包处理的差异表达基因列表转换为相应的蛋白质分子列表，用于构建PPI网络。我们将STRING网站输出的PPI网络输入到Cytoscape (3.8.1版)中，行可视化处理后进行下一步操作。

2.5. 核心基因的鉴定与表达

分子复合物检测工具(MCODE)可在Cytoscape中被开放获取，主要被用于识别PPI网络中联系更为紧密的部分，从而可以用于识别核心基因模块。筛选标准是node score cut-off = 0.2，maxdepth = 100，MCODE score > 2，cutoff = 2，kscore = 2，。此外，我们还应用cytoHubb，另一个可在Cytoscape中可开放获取的插件，可用于从MCODE的核心模块中进一步筛选，并用MCC算法将排名前十的基因排序，寻找最核心的基因。最后对得到的排名前十的基因在去批次化处理后的数据集中作热图，可视化观察核心基因在肿瘤组织和非肿瘤组织间的表达差异。

3. 结果

3.1. 差异转录基因(DEGs)

火山图显示了GSE99671和GSE197158中的差异转录基因(图1)。共有276个差异转录基因。

Figure 1. The volcano map of DEGs after batch removal of GSE197158 and GSE99671

图1. GSE197158和GSE99671去批次后的差异基因火山图

3.2. 差异转录基因注释

应用metascape进行的富集分析显示于图2，应用R语言进行的生物过程(BP)、细胞成分(CC)、分子功能(MF)富集分析显示于图3，而KEGG富集分析则显示于图4。与生物过程(BP)相关的差异表达基因变异主要集中在细胞外基质的合成、细胞骨架的合成、白细胞迁移、对细菌的反应、对真菌的杀伤作用、Killing细胞对其他组分的响应中(图2)。与细胞成分(CC)相关的差异转录基因主要集中在细胞外基质胶原、细胞质囊泡、特定克里囊泡、分泌囊泡等多种囊泡中(图2(b))。与分子功能(MF)相关的差异转录基因的变异主要集中在细胞外基质合成、信号识别激活调控、糖原聚糖结合、四吡咯结合、肝素结合等(图2(c))。与KEGG富集分析相关的差异转录基因变异主要集中在细胞因子–细胞因子受体相互作用、PI3K-Akt信号通路、癌症中的转录失调、PPAR信号通路、中性粒细胞胞外受体陷阱形成、心肌病等。

3.3. 蛋白质相互作用网络与核心基因

蛋白-蛋白相互作用(PPI)如图5所示。MCODE分析的核心模块如图5(a)所示。cytoHubb显示了在MCC算法中筛选出的前10个基因(图5(b))，包括CTSG、MPO、ELANE、DEFA4、CAMP、LTF、PRTN3、MMP8、LCN2、BPI。热图显示了将GSE99671和GSE197158数据集去批次化处理之后的合并数据集的核心基因表达情况(图6)。在骨肉瘤中，除DEFA4外，所有核心基因表达均下调。

Figure 2. Enrichment analysis of DEGs using metascape website. (a) Using metascape website, all DEGs were enriched into the corresponding items according to different enrichment items, and the p value was expressed by color depth. The darker the color, the smaller the p value, and the higher the reliability; (b) For a network composed of multiple enrichment items, different cluster identities of enrichment items give different colors to nodes. Nodes with the same cluster identity are closer and more connected; (c) For a network composed of multiple enrichment items, the color depth is defined by the p value. The more genes, the darker the color, the smaller the p value, and the higher the reliability

图2. 利用Metascape网站对DEGs进行富集分析。(a) 利用Metascape网站，将所有差异基因按不同的富集项分别富集到相应项目中，用颜色深浅表示p值高低，颜色越深，p值越小，可信度越高；(b) 由多个富集项组成的网络，以富集项的群集标识不同赋予节点不同的颜色，其中群集标识相同的节点相距更近，联系也越多；(c) 由多个富集项组成的网络，颜色的深浅由p值定义，包含的基因数目越多，颜色越深，p值越小，可信度越高

Figure 3. The DEGs were annotated by GO enrichment analysis based on r-packet. The number of genes was expressed by the size of the point, and the color was expressed by the q value

图3. 利用基于R包的GO分析对差异表达基因做出注释，用点的大小表示基因数量，颜色表示q值

Figure 4. The DEGs were annotated by KEGG enrichment analysis based on r-packet. The number of genes was expressed by the size of the point, and the color was expressed by the q value

图4. 利用基于R包的KEGG分析对差异表达基因做出注释，用点的大小表示基因数量，颜色表示q值

Figure 5. Protein protein interaction (PPI) network and core genes. (a) The core module of MCODE analysis in protein-protein interaction (PPI) networks; (b) The top 10 genes screened by MCC in cytoHubb

图5. 蛋白–蛋白相互作用(PPI)网络和核心基因。(a) 蛋白–蛋白相互作用(PPI)网络中的MCODE分析的核心模块；(b) 应用MCC在cytoHubb中筛选出的前10个基因

Figure 6. Heatmap of expression analysis of core genes

图6. 核心基因的表达分析热图

4. 讨论

骨肉瘤是严重威胁青少年健康的常见恶性肿瘤之一，发病率约占所有恶性肿瘤的0.2% [1]。其临床进展迅速，未出现转移时联合治疗生存率较高可达60%~70% [14]，但出现远处转移时其生存率大幅降低。且目前骨肉瘤的危险因素尚不明确，唯一明确的环境因素是电离辐射 [15]，因此从基因分子水平挖掘与骨肉瘤发生发展相关的基因对骨肉瘤治疗及改善预后具有重要意义。本研究利用生物信息学技术分析骨肉瘤和正常组织，筛选核心基因。本研究主要结果是CTSG和ELANE基因在骨肉瘤中低表达。

4.1. ELANE基因

ELANE基因属于丝氨酸蛋白酶家族，编码嗜中性粒细胞弹性蛋白酶(NE)，可以水解多种蛋白质底物，如细胞表面蛋白(G-CSF-Receptor、VCAM、c-kit和CXCR4)，另外其具有抗菌活性，参与细胞因子、趋化因子和生长因子的加工，包括肿瘤坏死因子α (TNF-α)和基质细胞衍生因子-1α (SDF-1α) [16] [17]。ELANE基因突变可导致周期性中性粒细胞减少症(CyN)和严重的先天性中性粒细胞减少症(SCN) [18]。Lerman等人提出在前列腺癌中，中性粒细胞可通过NE表达升高浸润组织引起癌症进展 [10]。有趣的是，另有研究提出NE能水解释放与组蛋白H1亚型相互作用的CD95死亡域(DD)，从而选择性地诱导多种癌细胞死亡，表达出抗癌活性，且可触发由CD8+T细胞介导的异位效应来防止远端转移 [19]。提出ELANE可作为选择性抗癌治疗的药物靶点 [20]。在本研究中我们发现ELANE基因在骨肉瘤中低表达与cui等人关于ELANE基因的抗癌作用相一致。因此我们推测ELANE基因表达下调可能使得该死亡域难以暴露，在癌细胞产生后难以诱导其水解，从而引起癌细胞继续进展，而癌细胞摄取ELANE减少，难以激活CD8+T细胞活化，原位进展和远处转移均失去一定的抑制。

4.2. CTSG基因

CTSG属于肽酶S1家族成员之一，编码组织蛋白酶G，具有特征性Ile16 N末端，在合成后以活性形式储存。其参与分解细胞外基质和血浆蛋白，与白细胞的分离及其迁移有关，还可作为单核细胞的趋化因子，促进炎症细胞的浸润，放大炎症反应 [11] [21] [22] [23]，有研究指出将CTSG作为一种关键的免疫基因，在膀胱上皮癌细胞组织中表达水平较低，从这可能帮助癌症细胞建立免疫逃避，引起膀胱上皮癌进展 [24]。另有报道在t (8; 21)急性髓性白血病(AML)中，白血病细胞可能通过AML1-ETO靶向抑制CTSG阻碍外源蛋白降解，从而逃脱免疫系统的监视 [13]。有系列研究将CSTG确定为口腔鳞状细胞癌(OSCC)中潜在的被确定为潜在的独立免疫生物标志物，也可作为治疗靶点 [12] [25]。还有研究利用蛋白质组学的方法解释了CTSG可用于区分恶性胸腔积液和结核性胸腔积液的潜在诊断性标志物 [26]。因此我们推测，CTSG基因表达水平的降低，可能影响免疫系统对癌症细胞的识别和清除，帮助癌症组织建立免疫逃避机制，从而引起骨肉瘤的发生和进展，并可作为潜在的诊断和治疗标志物，在未来作为骨肉瘤临床诊疗上做出帮助。

4.3. 研究的局限性

我们的研究局限性在于仅利用可开放获取的生物信息学数据库GEO数据库，对两个不同的骨肉瘤数据集行去批次化处理后再寻找相关核心基因。虽可信度比单一数据库更高，更尤其是对目前关于ELANE的免疫治疗对癌症的抑制有验证意义。但是并没有使用生物学技术进行验证，也没有研究相关通路。我们后续可以利用生物技术，在骨肉瘤动物模型上处理，对CTSG和ELANE基因行过表达或敲除操作，观察不同组小鼠的肿瘤是否有进展或控制，并推断基因在基因通路中的意义，从而指导后续药物开发临床治疗等一系列问题。

5. 结论

综上所述，ELANE和CTSG基因可能在临床骨肉瘤诊疗上有一定帮助，我们的研究可以为骨肉瘤诊治提供新思路。

NOTES

^*共同第一作者。

^#通讯作者。