卵巢癌发生发展关键基因的生物信息学筛选

余展鹏， 彭 亮， 姜 巧

卵巢癌是妇科最为常见的恶性肿瘤之一，其发病率逐年上升，病死率位居女性恶性肿瘤的第8位[1]。卵巢癌起病隐匿，大部分患者确诊时已为晚期，预后不佳，5年生存率仅为30%～40%[2,3]。为此，寻找与其发生、发展相关联的基因，对其诊断、治疗和预后评估具有重要意义。本研究拟使用生物信息学方法，整合分析高通量基因表达(Gene Expression Omnibus，GEO)数据库和癌症基因图谱(The Cancer Genome Atlas，TCGA)数据库及KMplot(Kaplan-Meier Plotter)等多个数据库中与卵巢癌的相关数据集，以筛选卵巢癌发生发展的关键基因。现报告如下。

1 资料与方法

1.1数据来源与差异基因(differentially expressed ge-nes,DEGs)筛选 (1)在GEO数据库(https://www.ncbi.nlm.nih.gov/geo/)中，按以下设置检索：(“ovarian neoplasms”[MeSH terms] or ovarian cancer[all fields]) and “homo sapiens”[porgn] and (“gse”[filter] and “expression profiling by array”[filter])。根据检索结果下载得到GSE54388、GSE105437、GSE40595数据集，3份数据均由GPL570平台完成、上传。它们分别包含6例、5例及8例正常卵巢上皮样本(正常组)和16例、10例及31例卵巢癌样本(卵巢癌组)。使用R 3.5.2软件进行数据预处理：以log FC>1.5或log FC<-1.5且校正P<0.05为标准，筛选正常组与肿瘤组的DEGs，并选出其中的共有部分，作为本研究的候选基因。(2)从Xena数据库[4](https://xena.ucsc.edu/)下载TCGA与GTEx(Genotype-Tissue Expression)的合并数据集(卵巢癌组=426，正常组=88)，分析两组候选基因的表达情况，进一步确定正常卵巢组织与卵巢癌组织的DEGs。

1.2DEGs与患者预后关系分析 使用KMplot数据库(https://kmplot.com/analysis/)卵巢癌数据集分析DEGs对患者预后的影响，以其自动计算各DEGs的cut-off值将患者分为高表达组和低表达组，并通过数据库缺省设置分别绘制各组Kaplan-Meier生存曲线，分析DEGs与患者生存预后的关系，筛选卵巢癌关键基因。

1.3卵巢癌关键基因的表达数据分析 从TCGA数据库(https://portal.gdc.cancer.gov/)中检索“Ovarian cancer”，下载380例卵巢癌组织的RNA-Seq数据集及与之对应的临床信息。剔除石蜡包埋标本，整理数据并进行标准化处理，分析卵巢癌关键基因表达水平与患者临床特征的关系。

1.4卵巢癌关键基因调控机制的初步分析

1.4.1 基因集富集分析 选取TCGA下载的卵巢癌组织RNA-Seq数据，运行GSEA 3.0(Gene Set Enrichment Analysis 3.0)软件，以两个卵巢癌关键基因[ARX(Aristaless Related Homeobox)，FAM150B(Family with Sequence Similarity 150 Member B)]的中位表达值为界，分为高表达组(190例)和低表达组(189例)，设置Hallmark Gene Sets为参照基因集，进行KEGG(Kyoto Encyclopedia of Genes and Genomes)信号通路聚类分析。

1.4.2 卵巢癌关键基因与免疫相关基因的关联分析 使用在线工具TIMER(https://cistrome.shinyapps.io/timer/)，选取“Ovarian Serous Cystadenocarcinoma”数据集，设置“tumor purity”为校正变量，分别使用GSEA富集得到的免疫系统相关信号通路的关键分子与ARX和FAM150B基因的表达水平进行Spearman相关分析。

1.4.3 卵巢癌关键基因的相互作用网络分析 使用STRING数据库(https://string-db.org/)分别检索“ARX”和“FAM150B”及与之相关的蛋白(即1.4.2步骤获得的信号通路的关键分子)，构建ARX和FAM150B的蛋白相互作用网络。

2 结果

2.1卵巢癌DEGs的筛选结果 从GEO数据库下载的GSE54388、GSE105437、GSE40595数据集中分别筛选到450、925、1097个DEGs。根据其上调、下调情况进行韦恩图分析，得到10个下调的候选基因及4个上调的候选基因。见表1。使用TCGA和GTEx数据集中的卵巢癌组织及正常卵巢组织的RNA-Seq数据对DEGs进行验证，结果显示有13个基因为卵巢癌的候选关键基因。见图1。

2.2卵巢癌DEGs与患者生存预后关联性分析结果 使用KMplot数据库分析以上13个DEGs与卵巢癌患者总生存期的关系，按检验水准α=0.01，结果显示，在下调的9个DEGs中，ADH1B低表达组的生存预后显著优于高表达组(P<0.01)；而ARX、FAM150B低表达组预后显著差于高表达组(P<0.01)。在上调的4个DEGs中，COL4A1高表达组预后显著优于低表达组(P<0.01)。见图2。

表1 卵巢癌DEGs基因

图1 卵巢癌DEGs验证结果图

图2 卵巢癌DEGs与患者生存预后的生存曲线图

2.3ARX和FAM150B表达水平与卵巢癌患者临床病理特征的关联性分析结果 结合2.2的结果，发现ARX和FAM150B表达量的下调，可能会导致患者的预后不良，鉴此，本研究选取ARX和FAM150B作进一步研究。结果显示，ARX的表达量与卵巢癌患者肿瘤直径具有关联性(P<0.05)，肿瘤直径越大，ARX表达量越低。FAM150B的表达量与患者的年龄及肿瘤直径具有关联性(P<0.05)，年龄或肿瘤直径越大，FAM150B表达量越低。见图3。

图3 ARX和FAM150B表达水平与卵巢癌患者临床病理特征的关联性分析结果图

2.4ARX及FAM150B调控机制的初步分析结果 使用TCGA卵巢癌RNA-Seq数据进行GSEA PATHWAY富集分析，结果显示，ARX低表达组较高表达组富集的通路有24条；FAM150B低表达组较高表达组富集的通路有14条。进一步分析发现，ARX与FAM150B低表达组富集到的通路中分别有12条和8条与免疫系统相关。整理以上免疫系统相关信号通路的关键分子，并使用TIMER数据库分析其与ARX和FAM150B表达量的相关关系。Spearman相关分析结果显示，ARX与PIK3R2(Phosphoinositide-3-Kinase Regulatory Subunit 2)的表达量呈正相关(P<0.05)，与IL-6及PIK3CD的表达量呈负相关(P<0.05)。见表2。FAM150B与PIK3R1的表达量呈正相关，与VAV1、PIK3CA、RAC2、MAPK1的表达量呈负相关(P<0.05)。见表3。使用STRING数据库分别构建ARX与IL-6、PIK3R2、PIK3CD，以及FAM150B与VAV1、PIK3CA、RAC2、PIK3R1、MAPK1的互作网络，预测结果显示，ARX可能可以通过EGFR对IL-6、PIK3R2及PIK3CD进行调控，FAM150B可能可以通过PLXNA2对VAV1、PIK3CA、RAC2、PIK3R1、MAPK1进行调控。见图4。

表2 ARX与免疫相关基因的Spearman相关分析结果

表3 FAM150B与免疫相关基因的Spearman相关分析结果

图4 ARX和FAM150B与免疫相关基因的互作网络图

3 讨论

3.1本研究通过联合GEO、TCGA及GTEx数据库分析发现ARX与FAM150B基因在卵巢癌中显著下调，且它们的表达量与肿瘤体积及患者预后密切相关，表达量越低，肿瘤体积越大，患者预后越差(P<0.01)。ARX基因定位于Xp22.1-Xp21.3之间，其包含5个外显子，编码含562个氨基酸的ARX蛋白，该蛋白为转录调控因子。既往研究[5]表明其在神经系统的发育中起重要作用，近来也有研究显示ARX基因与胰腺神经内分泌瘤、前列腺癌、卵巢癌等疾病的发生发展密切相关。FAM150B基因定位于2p25.3，其包含11个外显子，编码含152个氨基酸的FAM150B蛋白，该蛋白为ALK(anaplastic lymphoma receptor tyrosine kinase)和LTK(leukocyte receptor tyrosine kinase)的配体。既往研究[6]表明其在细胞的增殖、分化中起信号转导作用，但目前国内外与之相关的研究较少，仅发现其与成神经细胞瘤和鼻咽癌的发生、发展相关。

3.2本研究使用TCGA卵巢癌数据集中的数据，将获取数据分为高表达组与低表达组进行GSEA PATHWAY富集，结果发现两个基因的低表达组均大比例地富集到与免疫系统相关的信号通路。因此，我们推测这两个基因的下调可能会影响机体的免疫功能。鉴此，本研究进一步使用TIMER数据库分析ARX基因及FAM150B基因和与之相关免疫信号通路关键分子的表达情况。结果显示，在ARX基因下调的状态下，IL-6和PIK3CD基因的表达量上调，PIK3R2基因的表达量下调；而在FAM150B基因下调的状态下，VAV1、RAC2、MAPK1及PIK3CA基因的表达量上调，PIK3R1基因的表达量下调。其中PIK3R1和PIK3R2基因编码的蛋白为PI3K的调节亚基，PIK3CA和PIK3CD基因编码的蛋白为PI3K的催化亚基[7]。调节亚基的下调或催化亚基的上调都会导致PI3K的活性增强，从而减弱抗肿瘤免疫效应[8]。Vav1基因编码的蛋白是一种鸟嘌呤核苷酸转换因子，能促进Rho家族的Rac2表达量增加[9]，而Rac2又可以进一步上调MAPK1的表达[10]。MAPK1基因编码蛋白ERK2作为MAPK信号通路上的重要分子，其上调可使MAPK信号通路激活，继而上调PD-L1的表达量，使肿瘤细胞产生免疫耐受[11,12]。IL-6基因编码的IL-6作为炎症因子，在肿瘤微环境中不但能抑制树突细胞的迁移能力，减少其与T细胞的有效接触[13]，还能下调Th1细胞MHC Ⅱ类分子的表达,并减少IFN-γ和IL-2的释放，从而降低细胞毒性T细胞的活性[14],使肿瘤细胞能够逃避抗肿瘤免疫反应。

3.3为了探讨ARX和FAM150B调控上述基因的机制，本研究使用STRING数据库绘制了它们的互作网络，结果显示ARX可能通过TUBA1A及EGFR对IL-6、PIK3R2、PIK3CD进行调控；FAM150B可能通过PLXNA2对VAV1、PIK3CA、RAC2、PIK3R1、MAPK1进行调控。EGFR作为肿瘤发生、发展过程中的至关重要的基因，目前已被广泛研究[15]，而对于TUBA1A和PLXNA2的研究则相对较少，但近来亦有研究[16,17]表明TUBA1A和PLXNA2也是肿瘤发生发展中的关键调控基因。TUBA1A基因编码的α-微管蛋白是微管蛋白家族中重要的一员，是参与细胞分裂和运动所必需的蛋白。Zhu等[18]研究发现，在宫颈癌中TUBA1A可作为miR-15a-5p和miR-15b-5p的靶点，参与肿瘤细胞增殖及侵袭过程的调控。另外，PLXNA2既往一直只被认为在机体神经系统的发育过程中发挥关键调控作用，但目前也有越来越多的研究强调了其在肿瘤细胞侵袭和迁移过程中的调控功能。Tian等[19]分析基因芯片发现，与原发或局限性的前列腺癌相比，在转移性的前列腺癌组织中，PLXNA2基因显著上调。该研究还通过体外实验证实了PLXNA2参与了前列腺癌细胞迁移和侵袭的调控。因此，通过对上述结果的分析，我们推测在卵巢癌中下调的ARX和FAM150B可能可以通过调控TUBA1A、EGFR、PLXNA2，进而对IL-6、PIK3R2、PIK3CD、VAV1、PIK3CA、RAC2、PIK3R1及MAPK1进行调控，从而使机体抗肿瘤免疫反应减弱，最终导致肿瘤细胞出现免疫逃逸或者产生免疫耐受，促进卵巢癌的发生发展。

综上所述，本研究通过生物信息学分析发现ARX和FAM150B在卵巢癌组织中呈低表达状态，且它们的表达量和患者的预后密切相关，经过进一步分析，我们还发现ARX和FAM150B能对一系列基因进行调控，从而减弱机体的抗肿瘤免疫反应。在此研究中，所有用于分析的基因芯片原始数据均来自于同一平台，研究方法一致，数据客观，因此认为本研究能在一定程度上为ARX和FAM150B在卵巢癌中的进一步研究提供理论基础。