基于GEO对多发性骨髓瘤关键基因生物信息学分析及免疫浸润模式与验证

侯 丽，张 丽，唐 婧，牛 瑶，张 媛，朱有森（.新疆医科大学第一附属医院医学检验中心，乌鲁木齐 830054；.喀什地区第二人民医院检验科，新疆喀什 844000）

多发性骨髓瘤(multiple myeloma，MM) 是一种浆细胞恶性肿瘤，其特征表现为骨髓中克隆性浆细胞的恶性增殖，临床上主要表现为高血钙症、肾损害、贫血及骨破坏[1]。据报道，多发性骨髓瘤约占血液系统恶性肿瘤的10%，男性发病率高于女性[2]。2020年全球多发性骨髓瘤的年龄标准化发病率为1.78 / 10万人，死亡率为1.14/10万人[3]。由于多发性骨髓瘤存在较高的异质性，容易出现耐药和复发频繁，多发性骨髓瘤仍然无法完全治愈，五年生存率约为50%[4]。目前临床上的治疗方案，主要是通过早期化疗等方法控制病情，但治疗效果并不理想。一些研究发现免疫细胞参与了骨破坏和骨形成，并且指出骨细胞和免疫细胞之间存在双向调节的作用[5]，但各类免疫细胞，如巨噬细胞、粒细胞和非特异性免疫途径对其的影响尚不清楚。本研究运用生物信息学、机器学习算法等方法对公共数据库所基因表达综合数据库(gene expression omnibus，GEO)中多发性骨髓瘤相关的数据集进行挖掘，分析多发性骨髓瘤患者和健康对照组的基因表达微阵列数据，通过差异表达基因分析、功能富集和生存分析筛选关键基因，探索多发性骨髓瘤相关基因通路和功能的变化，揭示与其发生发展的分子机制提供理论依据。

1 材料与方法

1.1 资料来源 以“multiple myeloma”为关键词，从GEO数据库中检索诊断为多发性骨髓瘤的患者和健康对照样本的芯片，见表1。选取基因芯片GSE118985，GSE133346，GSE146649作为训练集，芯片数据生成于同一分析GPL570[HG-U133_Plus_2]Affymetrix Hunan Genome U133 Plus 2.0 Array，GSE7116以作为验证集。

表1 MM患者和健康对照基因芯片基本信息

1.2 方法与统计学分析

1.2.1 差异表达基因筛选：根据各数据集相应的平台探针信息注释基因，采用R 软件（4.2.2）的“limma”,“sva”软件包对三个训练集的测序原始数据进行批次校正均一化处理，并分析多发性骨髓瘤与健康对照组差异表达基因进行初筛，标准为P<0.05，且|Log2FC|≥1，在将各数据集筛选出的差异表达基因用稳健排序整合(robust rank aggregation，RAA)法排序，筛选出差异基因，使用R软件中的“pheatmap”软件包构建差异表达基因（differentially expressed genes，DEGs）的热图。

1.2.2 DEGs功能富集分析：基于筛选出的共同差异表达基因，通过R软件“stringi”“ggplot2”“clusterProfiler”“enrichplot”等软件包对获取的差异基因进行分析，基因本体论（gene ontology，GO）数据库对共同差异表达基因进行生物学功能注释；京都基因与基因组百科全书（kyoto encyclopedia of genes and genomes，KEGG）通路数据库进行共同差异表达基因信号通路的富集。

1.2.3 蛋白质互作网络构建：利用STRING 11.5（https://cn.string-db.org/）构建共同差异表达基因编码蛋白的互作网络，筛选置信度得分≥0.4的互作蛋白，去除网络中的游离蛋白，并用Cytoscape软件中CytoHubba插件依据项目进行打分筛选出其中的关键基因（hub gene）。

1.2.4 应用受试者工作特性曲线评价：利用GSE7116中的数据构建疾病-对照的模型验证集，评价hub基因与多发性骨髓瘤的关联程度。利用R软件绘制应用受试者工作特性（receiver operator characteristic,ROC）曲线并通过计算ROC曲线下面积对筛选出的核心基因进行评价。

2 结果

2.1 差异表达基因筛选结果 见图1。通过筛选共获得101个差异表达基因，其中25个下调基因，76个上调基因，并将排名前20的上调及下调差异基因绘制差异表达基因的热图。

图1 多发性骨髓瘤患者与健康对照差异表达基因热图

2.2 差异表达基因功能富集及通路分析 见表2。对筛选出的101个共同差异表达基因进行GO和KEGG通路富集分析，结果显示，差异基因BP主要富集在细胞因子的正向调节、B 细胞活化、T细胞、免疫效应的调节及蛋白质磷酸化；CC主要富集在内质网腔，细胞外基质；MF主要富集在信号受体激活，细胞因子激活及生长因子激活等；共富集到KEGG通路12个，主要涉及细胞因子受体相互作用、细胞外基质受体相互作用、粘着斑、Hedgehog信号通路、致心律失常性右心室心肌病(ARVC)等通路。

表2 共同差异表达基因GO和KEGG通路富集分析结果

2.3 蛋白质互作网络的构建及hub基因的筛选 将101个差异表达基因输入STRING进行互作网络分析，去除游离的蛋白后，共得到由89个节点（靶点蛋白）和61条边（蛋白质互作）构成的蛋白质互作网络，见图2。再进一步应用Cytohubba计算打分，筛选出蛋白质互作网络中的关键基因，皆为上调基因，分别为多配体蛋白聚糖1(Syndecan,SDC1)、干扰素调节因子4（interferon regulatory factor，IRF4）、分化抗原簇38(Cluster of Differentiation 38，CD38)、肿瘤坏死因子受体超家族(tumor necrosis factor receptor superfamily，TNFRSF17)和细胞周期素(cyclin D1，CCND1)，并将hub基因绘制差异表达基因的热图，见图3。

图2 蛋白质互作网络关键基因簇分析

2.4 Hub基因的诊断意义评价 用筛选出的核心基因对验证集GSE7116绘制ROC曲线，见图4。SDC1，IRF4，CD38，TNFRSF17和CCND1的曲线下面积(area under curve，AUC)分别为0.714，0762，0.714，0.590，0.686，核心基因对验证模型联合诊断的AUC为0.933，说明筛选出的核心基因对MM具有较好的诊断意义，

图4 Hub基因对MM诊断的ROC曲线

2.5 免疫细胞浸润分析 根据具有22种免疫细胞的基因表达矩阵，用R软件计算出三个训练集中MM组和健康对照组的免疫细胞浸润情况，结果显示，M2型巨噬细胞和单核细胞相互作用关系最明显，见图5A。与正常人群相比，多发性骨髓瘤患者有12类免疫细胞差异具有统计学意义，其中幼稚性B淋巴细胞、记忆性B淋巴细胞、幼稚性CD4，γδT细胞、单核细胞及中性粒细胞免疫浸润差异程度较大，差异具有统计学意义（均P＜0.01），见图5B。

图5 22种免疫细胞浸润结果分析

3 讨论

多发性骨髓瘤(MM)是一种从骨髓中克隆性浆细胞的异常增殖为特征的血液学恶性肿瘤，主要临床表现为血细胞减少，血液或尿液中的单克隆蛋白增多，骨溶解病变[6-7]。多发性骨髓瘤是一种无法根治的浆细胞克隆性恶性肿瘤，尽管使用蛋白酶体抑制剂和免疫调节剂改善了多发性骨髓瘤的治疗，但几乎所有最初治疗存活下来的患者最终都会复发，需要进一步治疗[8]。因此，筛选与多发性骨髓瘤相关的早期诊断标志物和调控途径非常重要。

本研究检索GEO数据库中多发性骨髓瘤相关的数据集，通过生物信息学分析筛选多发性骨髓瘤患者与健康人群差异表达基因，筛选出的差异基因主要富集在免疫反应过程中，包括细胞因子的正向调节、B 细胞活化、T细胞、免疫效应的调节及蛋白质磷酸化等，这与血液系统肿瘤患者存在严重的免疫功能紊乱、细胞周期调控机制被破坏[9]相符合。SDC1是一种硫酸肝素蛋白聚糖，属于多配体蛋白聚糖家族，是维持细胞形态所必需的细胞表面黏附分子，通过硫酸肝素链与黏附分子、基质组分、生长因子、酶、酶抑制剂等一系列配体结合，与周围微环境相互作用。SDC1主要在上皮细胞中表达，因此，SDC1表达异常会促进细胞增殖、转移、侵袭和血管生成而导致癌症的发展[6,10]，有研究表明，SDC1在多发性骨髓瘤中已被普遍用作浆细胞标志物[11]。IRF4位于6p25.3位点，是淋巴细胞发育并分化成分泌IgG的浆细胞的重要因子，在细胞增殖、凋亡、致癌易感性和T细胞免疫反应中起主要作用[12]。有研究表明[13]，骨髓瘤细胞的生存依赖IRF4，IRF4基因缺失会导致浆细胞凋亡，抑制IRF4的表达或干扰其转录可作为治疗各种亚型多发性骨髓瘤药物的新靶点。外周血中90%的浆细胞为CD38阳性，自然杀伤细胞和单核细胞中约有60%表达CD38[14]，与正常淋巴细胞和骨髓细胞相比，CD38在多发性骨髓瘤患者浆细胞中呈现出高特异性和高表达的状态，该位点可能成为多发性骨髓瘤的潜在治疗靶点，因此，一些CD38的单克隆抗体药物已经被开发出来如达雷妥尤单抗，并且已经应用到了多发性骨髓瘤的临床治疗。TNFRSF17是B细胞表面分子，也叫B细胞成熟抗原，属于肿瘤坏死因子受体家族，主要表达于成熟的B淋巴细胞，对于B细胞的成熟和自身免疫反应发挥重要作用[15]。TNFRSF17在幼稚性B细胞、造血干细胞或正常的非血液组织中检测不到，它与B细胞活化因子及其受体和跨膜激活剂、钙调节剂和亲环蛋白配体相互作用调节B细胞增殖、成熟、存活以及向浆细胞分化。CCND1是一种编码细胞周期调节蛋白的基因，是正常细胞和癌细胞由G1期向S期发展所必需的蛋白，CCND1在细胞周期进展中的作用已得到广泛认可，但在癌症发生、增殖、迁移、侵袭、转移中的作用尚未完全阐明，它在多发性骨髓瘤中作为癌基因被激活，并调节细胞周期和促进增殖。已有研究表明，CCND1的高表达与乳腺癌、结直肠癌、食管癌预后也存在一定的关联[16]。破骨细胞为介导多发性骨髓瘤的主要效应细胞，它起源于单核-巨噬细胞系，是一类多核且高度分化的对骨质具有再吸收作用的细胞，骨髓内破骨细胞受到来自恶性浆细胞或骨髓微环境中其他细胞的刺激出现增殖且功能活跃，导致溶骨亢进。本研究显示，单核细胞是多发性骨髓瘤患者中最重要的免疫浸润细胞，M2型巨噬细胞和单核细胞相互作用关系最明显，这与该疾病的主要效应细胞的起源存在一致性，其次幼稚性B淋巴细胞、记忆性B淋巴细胞、CD4+幼稚性T淋巴细胞、γδ T细胞及中性粒细胞免疫浸润均存在明显的差异，本研究训练集基于同一分析平台三个数据集，经过数据校正后，筛选出的5个核心基因并没有集中富集在某几条通路上，SDC1主要参与肌细胞分化，细胞外基质受体相互作用；IRF4主要参与免疫效应过程的调节，细胞因子的正向调节，T细胞活化的调节；CD38主要参与B细胞活化；TNFRSF17主要参与细胞因子受体相互作用；CCND1主要参与hedgehog信号通路、粘着斑通路，均富集在免疫反应过程，这可能由于在分子水平上，多发性骨髓瘤并不是一种单一的疾病，而是一种特定的基因表达和染色体易位反复发生的多亚型表现的异质性疾病[17-18]。筛选出的5个核心基因均有研究报道在多发性骨髓瘤中的特异性表达，并有可能成为多发性骨髓瘤的潜在治疗靶点，在外部数据集验证了核心基因在多发性骨髓瘤与正常人群中的差异性，ROC曲线显示在外周血单个核细胞中核心基因对多发性骨髓瘤患者同样也具有较好的诊断意义。

综上所述，本研究采用生物信息学方法分析现有的多发性骨髓瘤基因芯片数据，通过分析预测了多发性骨髓瘤可能的发病机制并筛选出可能参与多发性骨髓瘤发生发展的核心基因，筛选出5个关键基因，可能成为多发性骨髓瘤辅助诊断的新的标志物或新药物治疗的靶点。本研究完全基于GEO数据库中的数据，筛选出的关键基因与多发性骨髓瘤的相关性及相关发病机制仍需在临床样本中进行验证。