多发性骨髓瘤(Multiple Myeloma，MM) 是一种以骨髓中单克隆浆细胞增殖为特征的血液系统恶性肿瘤。MM细胞产生的大量单克隆免疫球蛋白使内质网聚集了大量未折叠或错误折叠的蛋白。近年来，随着蛋白酶体抑制剂等新型药物的应用，MM的生存率有了显着提高^[1]。硼替佐米、来那度胺和地塞米松是目前新诊断MM的一线治疗标准^[2]。蛋白酶体抑制剂、免疫调节物质和经典化疗药物是复发性MM主要的治疗措施^[3]。但即使在这些治疗下，几乎所有MM患者最终都会复发，而且不同患者间的预后存在高度异质性。因此，对MM患者进行风险分层，寻找可靠的预后生物标志物，可以更好地提高预后的准确性，指导临床治疗具有重要意义。当前针对MM的预后评估系统较多，包括国际分期系统(International Staging System，ISS)、R-ISS分期、梅奥mSMART分期、Durie-Salmon (DS)分期等。近年来也有学者提出将炎症指数评分及营养状况评分纳入到MM患者的预后评估中^[4-5]。但上述评分系统尚不能反映MM患者疾病状态的全貌。自噬作为一种细胞自我清除的机制，可使细胞质内错误折叠的蛋白质和受损细胞器自动降解，对MM细胞的生存至关重要^[6]。此外，营养缺乏、缺氧、DNA损伤或活性氧等不利于细胞生存的条件均会触发自噬，克服细胞能量缺乏，清除多余物质，达到维持细胞存活的目的^[7]。研究发现，通过敲除Beclin-1基因抑制自噬或使用自噬抑制剂3-甲基腺嘌呤和氯喹可诱导MM细胞凋亡，阻止自噬小体的形成，但同时也发现过度自噬可以导致MM细胞死亡^[8-9]。因此，在MM的发生发展中，自噬似乎是一把双刃剑。最近的研究显示，通过识别自噬相关基因(Autophagy-Related Genes，ARGs)可以有效识别不同风险级别的类风湿性关节炎和结肠癌等^[10-11]。目前，国内尚无基于ARGs评分与MM预后关系的相关研究，本研究拟采用生物信息学分析方法，通过MM患者转录组与临床数据，系统分析MM ARGs的差异表达并预测MM预后的生物标志物，建立MM预后模型。

为确定可预测MM患者预后的ARG，本研究以：(1)数据来源于初次诊断MM患者；(2)至少10个生物学重复；(3)数据为原始fastq数据采用HTSeq2比对定量后的FPKM文件，三点为标准，以“Multiple Myeloma/MM”为关键词，在Gene Expression Omnibus(GEO, https://www.ncbi.nlm.nih.gov/)、The Cancer Genome Atlas(TCGA, https://portal.gdc.cancer.gov/)、National Cancer Institute center for cancer Research(NCICCR, https://portal.gdc.cancer.gov/repository)、Human Protein Atlas(HPA, https://www.proteinatlas.org/)和Multiple Myeloma Research Foundation(MMRF, https://themmrf.org/)数据库进行检索。

应用R语言“limma”软件包，以| log₂FC |>1(FC表示两样品间表达量的比值)，错误发现率(False Discovery Rate，FDR)＜0.05为条件筛选多发性骨髓瘤患者骨髓样本与正常骨髓样本中的差异ARGs(Differential Expression Autophagy-Related Genes，DEARGs)，绘制差异表达火山图。

使用R软件“clusterprofiler”程序包对DEARGs进行GO和KEGG功能富集分析，以adj.P＜0.05为标准对结果进行筛选。

应用R语言“survival”软件包中的单因素COX回归算法以P＜0.05为筛选条件，分析与MM患者预后相关ARGs，并绘制森林图。随后对与预后生存相关的ARGs构建多因素COX模型，计算每个MM样本的风险值，并依据风险值的中位值将所有样本分为高风险组和低风险组，绘制两组生存曲线。

应用单因素和多因素COX分析方法评价患者年龄、性别、分期、家族史、风险评分对MM患者预后的影响，绘制森林图。应用多指标受试者操作特征(Receiver Operating Characteristic，ROC) 曲线评价这些因素独立预测MM患者预后的准确性。

从MMRF数据库(https://themmrf.org/) 下载764例MM患者的骨髓组织转录组和临床数据，结合4例来自人类蛋白质图谱(Human Protein Atlas，HPA，https://www.proteinatlas.org/)的正常骨髓样本。基于人类自噬数据库(Human Autophagy Database，HADb) (http://www.autophagy.lu/)中的232个ARGs，构建MM ARGs表达谱，用于后续分析。

应用R语言“limma”软件包从764例MM骨髓样本及4例正常骨髓样本中提取差异表达ARGs，同时提取临床资料完整的744例MM患者的年龄、性别、分期、家族史等临床信息(见表 1)。结果显示：有104个基因的表达在MM样本中具有显著差异(|log₂FC|>1，FDR＜0.05)，其中上调基因46个，下调基因58个(见图 1)。

表1(Table 1)

表 1 纳入研究患者基本资料 Table 1 Basic data of patients included in the study 临床特征 例(%) 年龄 ≤65 461(61.962) ＞65 283(38.038) 性别 男 438(58.871) 女 306(41.129) ISS分期 Ⅰ期 262(35.215) Ⅱ期 270(36.290) Ⅲ期 212(28.495)

表 1 纳入研究患者基本资料 Table 1 Basic data of patients included in the study

图 1(Figure 1)

图 1 MM患者骨髓样本与正常对照骨髓样本中ARGs差异基因火山图及热图 Figure 1 Volcano map and heat map of ARGs differential genes in bone marrow samples of MM patients and normal controls 注：红点代表上调的差异ARGs，绿点代表下调的差异ARGs，黑点表示在MM骨髓样本中未差异表达的ARGs.

对104个DEARGs进行功能富集分析，GO富集分析显示，差异表达ARGs主要集中在巨自噬、自噬调节、细胞对外部刺激的反应等本体学注释(见图 2)。KEGG通路显示，DEARGs主要集中在自噬、细胞凋亡、NOD样受体信号通路、PI3K-Akt信号通路等(见图 3)。

图 2(Figure 2)

图 2 GO功能富集分析结果 Figure 2 Results of GO enrichment analysis 注：图中ID号码为GO通路的编号，红色为富集通路中上调基因，蓝色为下调基因。

图 3(Figure 3)

图 3 KEGG功能富集分析结果 Figure 3 Results of KEGG enrichment analysis 注：横坐标代表富集在各通路的基因数目，纵坐标代表通路名称，P＜0.05。

对744个临床信息完整的样本进行生存分析。使用上述104个DEARGs单因素COX分析共发现33个ARGs与MM患者整体生存明显相关，其中13个ARGs是风险基因(HR＞1)。进一步使用多因素COX回归筛选出13个预后相关ARGs(NKX2-3、NCKAP1、BIRC5、PEX3、HGS、RUBCN、PARP1、ARSA、DNAJB9、HSP90AB1、EEF2、FKBP1B、CD46)，详见表 2。根据每个基因的mRNA表达水平和风险系数计算每个患者的风险值，风险值计算公式为：风险值=(1.911) ×NKX2-3表达值+(-0.212)×NCKAP1表达值+(0.189)×BIRC5表达值+(-0.253)×PEX3表达值+(0.353)×HGS表达值+(-0.493)×RUBCN表达值+(0.406)×PARP1表达值+(-0.192)×ARSA表达值+ (-0.199)×DNAJB9表达值+(0.318)×HSP90AB1表达值+(-0.29)×EEF2表达值+(-0.246)×FKBP1B表达值+ (-0.162)×CD46，将MM患者分为高风险组和低风险组，Kaplan-Meier生存曲线分析显示，高、低风险组患者的生存率比较，差异有统计学意义(P＜0.001) (见图 4)。

表2(Table 2)

表 2 单因素和多因素COX分析结果 Table 2 Results of univariate and multivariate COX analyses 基因 单因素分析结果 多因素分析结果 HR(95%CI) P value HR(95%CI) P value 系数 NKX2-3 3.476(1.113，10.856) 0.032 6.762(1.697，26.946) 0.006 1.911 FKBP1B 0.590(0.427，0.817) 0.001 0.782(0.558，1.097) 0.155 -0.246 NCKAP1 0.748(0.651，0.858) ＜0.001 0.809(0.701，0.934) 0.003 -0.212 BIRC5 1.582(1.390，1.802) ＜0.001 1.208(1.038，1.406) 0.015 0.189 PEX3 0.808(0.660，0.989) 0.038 0.777(0.592，1.019) 0.068 -0.253 HGS 1.436(1.055，1.956) 0.021 1.423(0.981，2.066) 0.063 0.353 RUBCN 0.695(0.550，0.879) 0.002 0.611(0.470，0.793) ＜0.001 -0.493 PARP1 1.435(1.127，1.826) 0.003 1.500(1.135，1.983) 0.004 0.406 ARSA 0.732(0.602，0.891) 0.002 0.825(0.666，1.024) 0.080 -0.192 CD46 0.795(0.649，0.975) 0.028 0.850(0.683，1.058) 0.146 -0.163 DNAJB9 0.623(0.496，0.783) ＜0.001 0.819(0.623，1.077) 0.153 -0.200 HSP90AB1 1.612(1.178，2.208) 0.003 1.375(0.977，1.935) 0.068 0.318 EEF2 0.700(0.578，0.849) ＜0.001 0.746(0.613，0.907) 0.003 -0.293

表 2 单因素和多因素COX分析结果 Table 2 Results of univariate and multivariate COX analyses

图 4(Figure 4)

图 4 高风险组和低分险组MM患者Kaplan-Meier生存曲线 Figure 4 Kaplan-Meier survival curve of MM patients in high risk group and low risk group

单因素分析显示，MM患者的年龄、性别、分期和风险值与OS有关(P＜0.05)。此外，根据临床病理特征年龄、性别、分期、家族史和风险值进行多因素COX回归分析表明，年龄、分期、风险值可作为评估MM患者预后的独立因素(P＜0.05，见图 5)，ROC曲线下面积分别为0.561、0.685、0.719(见图 6)。

图 5(Figure 5)

图 5 多因素COX回归分析影响745例患者总体生存期的因素 Figure 5 Analysis of factors affecting overall survival of 745 patients by multivariate COX regression

图 6(Figure 6)

图 6 5年生存多指标ROC曲线 Figure 6 Five-year survival multi-indicator ROC curve 注：横坐标为假阳性率，纵坐标为真阳性率(灵敏度).

DNA测序数据显示，不同多发性骨髓瘤患者或同一患者随发病时间均存在巨大的遗传异质性^[12]。因此，当前对MM的生物学异质性及其对治疗的影响仍有很大的研究空间。尽管近年来出现的蛋白酶体抑制剂、免疫调节药物和单克隆抗体等许多新的治疗选择提高了MM患者的中位生存率，但仍有25%的新诊断患者生存率低于3年，且化疗耐药和复发仍然是MM临床治疗中的主要问题^[13]。Deng等在血液系统恶性肿瘤中研究发现，自噬在化疗耐药方面起着至关重要的作用^[14]。本研究使用MMRF数据库中764例MM患者的骨髓组织表达谱和744例MM患者临床信息，结合HADb数据库中的232个ARGs，最终得到NKX2-3、NCKAP1、BIRC5、PEX3、HGS、RUBCN、PARP1、ARSA、DNAJB9、HSP90AB1、EEF2、FKBP1B和CD46等13个预后相关基因，通过多因素COX分析和多指标ROC分析验证了年龄、分期和风险值可作为独立预后标志。

丝氨酸/苏氨酸蛋白激酶AMPK是细胞中的能量传感器，在分解代谢的激活和合成代谢的失活中起关键作用。研究表明，在各种生理和病理条件下，AMPK被上游激酶磷酸化并与AMP或ADP结合后活化，之后通过磷酸化mTORC1，ULK1和PIK3C3/VPS34复合物中的自噬相关蛋白来直接促进自噬，或通过调节转录因子如FOXO3，TFEB和BRD4等下游的ARGs的表达来间接促进自噬^[15]。本研究中富集分析也显示大量差异表达ARGs富集到丝氨酸/苏氨酸蛋白激酶活化信号通路。NKX家族成员在淋巴器官发育中起着至关重要的作用^[16]。尽管目前研究主要集中于NKX2-3与炎症性肠病的关系，但已有研究表明B细胞中过表达NKX2-3的转基因小鼠可通过NF-κB和PI3K-AKT途径促进B细胞向脾和结外组织迁移，导致边缘区淋巴瘤的进展^[17]。本研究发现NKX2-3与MM患者预后不良高度相关，此外KEGG富集分析也发现了差异表达ARGs在PI3K-AKT信号通路的大量富集.研究表明PI3K-Akt信号通路和自噬之间存在密切关联，PI3K-Akt信号通路的激活可能会在多种癌症中诱导自噬^[18-19]。Zhong等^[20]研究发现，NCKAP1通过增强肝细胞癌中的Rb1/p53激活改善患者预后并抑制肝癌细胞生长，本研究中也发现NCKAP1对于MM患者的预后是一种保护性因素，但具体机制还有待进一步研究。MM治疗中产生的耐药是影响MM患者预后的重要因素，研究表明，BIRC5表达的上调与骨髓基质细胞诱导的MM耐药高度相关^[21]。缺氧诱导因子1(Hypoxia inducible factor 1，HIF-1)是研究最多的ARGs之一，缺氧是癌症发生和发展的主要特征。据报道，MM中的缺氧会增强肿瘤的发生，上调原始B细胞阶段基因的表达，参与蛋白酶体抑制剂耐药的发生^[22-23]，刺激恶性浆细胞迁移^[24]，促进溶骨性骨破坏及血管生成^[25]。本研究中发现的HIF-1信号通路也再次表明HIF-1信号通路可能是MM治疗的潜在靶点。除自噬外，ARGs可能还具有多种功能，本研究进行的GO和KEGG通路分析也发现了除自噬外的某些其他注释，包括细胞对刺激的反应，慢性髓系白血病，微生物感染和神经退行性疾病等。

虽然目前已有几种包含临床特征和分子生物学的MM危险分层预后模型如Durie-Salmon分期系统和ISS分期系统等，但尚无模型能够有效预测MM患者的预后。以往研究表明，年龄是MM预后的独立危险因素^[26]，但在我们的研究中，年龄相较于由ARGs得出的风险值预后意义较小。而ISS分期的ROC值则更接近于风险值。此外，患者的肾功能状态，是否存在意义未名的单克隆免疫球蛋白血症病史，是否存在髓外疾病及影像学指标等也均与MM患者预后高度相关^[27]。近几年来，研究人员尝试利用ARGs评估MM患者的预后情况。Jung G等^[28]检测了89例MM患者的骨髓活检标本，发现骨髓中自噬标志物Beclin 1和微管相关蛋白轻链3的高表达与更好的预后相关，但并未建立相应的预后评分系统。Zhu等^[29]则利用GEO数据库中559例MM患者数据建立了包含16个ARGs的MM患者预后评分系统，5年生存率的ROC曲线下面积与本研究类似，同时找到了与本研究中相同的FKBP1B、NCKAP1和DNAJB9基因。与上述研究相比，本研究中利用的MMRF数据库样本量更大，对上述研究进行了有效补充。

综上，虽然当前对于MM的治疗及预后仍存在极大的挑战，但新的治疗方法同样为MM患者带来了极大的希望。因此，我们看待MM预后不良的问题也需要综合多方面因素综合考虑。有理由相信，随着越来越多的治疗方案问世，高危MM中的某些亚型在日后将有可能成为低危MM。而基于RNA的定量转录组分析也可以发现更多与突变无关的，涉及多基因，多途径的致病机制。DNA和RNA的综合分析也将进一步为我们提供疾病的生物学进展及预后情况。

NKX2-3、NCKAP1、BIRC5、PEX3、HGS、RUBCN、PARP1、ARSA、DNAJB9、HSP90AB1、EEF2、FKBP1B和CD46等可能是与MM患者预后相关对ARGs，据此建立的预后模型对于临床具有一定指导意义，但需进一步验证。