多发性骨髓瘤(multiple myeloma,MM)是骨髓浆细胞恶性克隆性增殖的血液系统肿瘤。其治疗首选蛋白酶体抑制剂(proteasome inhibitors,PIs)、免疫调节剂和地塞米松
[1]。近几十年来,随着治疗药物的开发,患者预后显著改善
[2-3],但MM的遗传异质性导致患者的化疗敏感性存在差异,复发难治仍是MM治疗的难点
[4]。硼替佐米(bortezomib,BTZ)作为新诊断多发性骨髓瘤(newly diagnosed multiple myeloma,NDMM)患者的首选PIs类药物,与其他药物联合广泛应用于MM治疗
[5]。然而,在BTZ治疗的过程中,很多患者会出现获得性耐药,甚至部分初次治疗的患者也表现出耐药性。MM患者在疾病不同阶段均可复发,经过多轮治疗后产生的耐药性和治疗相关并发症可能导致复发难治MM患者死亡。其中,初始诱导治疗无效MM患者的预后通常比初始缓解者更差
[6]。这部分初始诱导治疗失败患者即为原发难治性多发性骨髓瘤(primary refractory multiple myeloma,PRMM)患者。PRMM患者通常在未检测到传统高危因素(如细胞遗传学异常和β2微球蛋白异常等)的情况下,也表现出较差的生存结局
[6]。由于现有治疗药物(如BTZ、环磷酰胺、来那度胺和地塞米松等)的联合应用对大部分初诊患者具有较好的疗效,因此在临床上发现的PRMM患者相对较少,当前的临床研究也缺乏针对这部分患者生存结局的系统评估
[6]。
研究
[7]揭示代谢紊乱可能在MM的进展和耐药中发挥关键作用。代谢组学技术和生物信息学分析技术的联合应用可为挖掘这些代谢变化提供有力的技术支持。此外,结合公共数据库的测序数据可分析MM进展过程中的分子事件,并为MM筛选新的预后生物标志物
[8]。因此,本研究聚焦于对含BTZ化疗方案诱导治疗无反应或反应较差的PRMM患者,基于不同疗效MM患者代谢组学的分析结果,结合数据库转录组测序数据,利用生物信息学分析及机器学习算法,构建预测NDMM患者疗效的机器学习模型,旨在为PRMM患者的治疗提供新的思路,同时为耐药机制的研究提供理论依据。
1 对象与方法
1.1 伦理声明
本研究已获得中南大学湘雅三医院伦理委员会批准(审批号:22298)。患者均签署知情同意书。
1.2 对象
收集2022年8月至2023年7月在中南大学湘雅三医院血液内科住院的MM患者。在中南大学湘雅三医院检验科收取患者检测后的剩余血清,离心采集后即刻置于冰上,储存于-80 ℃冰箱。在电子病历系统中收集患者临床信息,判断患者疗效。
纳入标准:接受多药联合化疗的MM患者,初始诱导方案为VCD(BTZ+环磷酰胺+地塞米松)或VRD(BTZ+来那度胺+地塞米松)。化疗2~3个疗程或4~6个疗程。
排除标准:化疗方案不含BTZ及更换化疗方案后仍含BTZ患者;多次复发,化疗方案复杂的患者;行自体造血干细胞移植的患者。
定义:NDMM,初次诊断为MM,暂未经治疗。治疗有效MM(treatment response multiple myeloma,TRMM),初次化疗起使用含BTZ诱导方案,疗效达非常好的部分缓解(very good partial response,VGPR)。PRMM,初次化疗起使用含BTZ诱导方案,对初始治疗无反应或反应差,疗效未达VGPR或早期复发。
1.3 代谢组学检测
液相色谱采用Thermo Vanquish(Thermo Fisher Scientific,USA)超高效液相系统。质谱采用Thermo Orbitrap Exploris 120质谱检测器(Thermo Fisher Scientific,USA)。测得患者血清代谢物测序的原始数据后,首先通过Proteowizard软件包(v3.0.8789)中MSConvert工具将原始质谱下机文件转换为mzXML文件格式。采用R软件XCMS包进行峰检测、峰过滤和峰对齐处理,得到物质定量列表。采用基于质控样本的方法实现数据校正,消除系统误差。然后保留质控样本中变异系数小于30%的物质进行后续分析。评估代谢物对样本分类的影响力和解释力,根据统计检验计算P值,采用正交偏最小二乘法判别分析(orthogonal partial least squares discriminant analysis,OPLS-DA)降维计算变量投影重要度(variable importance in projection,VIP),差异倍数分析计算组间差异倍数。以P<0.05且VIP>1筛选差异代谢物。采用“MetaboAnalyst”软件包进行功能富集分析。
1.4 转录组数据处理
在基因表达综合(Gene Expression Omnibus,GEO)数据库中下载数据集:GSE68871(n=118)、GSE55145(n=67)、GSE116324(n=44)、GSE159426(n= 51)、GSE9782(n=239)和GSE136337(n=426)。使用R软件处理芯片数据集和高通量测序数据集。GSE68871和GSE55145为同一平台芯片测序数据,使用“ComBat”包处理批次效应,合并两数据集用于差异表达基因(differential expression genes,DEGs)筛选。采用主成分分析(principal component analysis,PCA)观察合并数据集处理批次效应前后的样本分布情况。GSE68871、GSE55145、GSE116324和GSE159426中的患者均采用含BTZ化疗方案,初诊时采集。GSE9782中的患者为复发MM患者,参与了BTZ 2期和3期临床试验,试验前采集。GSE116324、GSE159426和GSE9782用于外部验证。GSE136337用于预后分析。
1.5 DEGs筛选
基于GSE68871和GSE55145合并后的转录组数据,使用“limma”包筛选PRMM组和TRMM组DEGs,筛选条件为P≤0.05。在京都基因和基因组数据库(Kyoto Encyclopedia of Genes and Genomes,KEGG)及分子标签数据库(Molecular Signatures Database,MSigDB)中获取氨基酸代谢相关基因(amino acid metabolism genes,AAMGs)。AAMGs与DEGs取交集,筛选不同疗效患者间差异表达AAMGs,用于构建模型。
1.6 预后分析
使用“limma”包从GSE136337中提取差异表达AAMGs的表达量。“survival”和“survminer”包分析生存数据,筛选与预后相关的DEGs(P<0.05)。得到Kaplan-Meier(K-M)生存曲线分析结果,绘制基因的K-M生存曲线和预后网络图。
1.7 模型构建及验证
采用“glmnet”包进行最小绝对收缩和选择算子(least absolute shrinkage and selection operator,LASSO)回归分析,初步筛选预后相关的差异表达AAMGs。基于机器学习算法构建疗效预测模型,采用“randomForest”包构建随机森林(random forest,RF)模型,“kernlab”包构建支持向量机(support vector machine,SVM)模型,“xgboost”包构建极致梯度提升(extreme gradient boosting,XGB)模型,“glm()”函数构建广义线性模型(generalized linear model,GLM)。采用“pROC”包绘制受试者操作特征(receiver operating characteristic,ROC)曲线,评价模型的预测性能。采用“rms”包构建所筛模型基因的列线图,“rmda”包构建模型的临床决策曲线。
1.8 富集分析与免疫分析
使用“clusterProfiler”包进行KEGG和GO富集分析,其中基因本体(Gene Ontology,GO)数据库注释DEGs,KEGG分析DEGs涉及的功能和通路。使用基因集变异分析(Gene Set Variation Analysis,GSVA)计算每个样本中特定基因集合的富集程度。为进一步挖掘模型基因在MM中的潜在作用机制,利用GSVA的富集分析鉴定基因高、低表达组中差异富集的通路。GSVA分析能够识别因基因表达水平变化而差异调节的途径,进而确定受干扰的生物学过程。使用“GSVA”包进行GSVA富集分析。使用“CIBERSORT”包进行免疫浸润分析,计算不同队列的免疫功能评分。
1.9 统计学处理
采用SPSS 27.0统计学软件进行数据分析。计数资料以例表示,计量资料以均数±标准差表示。采用R 4.3.1软件进行生物信息学分析。DEGs筛选采用Wilcoxon检验。预后分析采用Cox比例风险回归模型(proportional-hazards model)和K-M生存曲线分析。基因高表达和低表达患者的总生存期(overall survival,OS)比较采用K-M生存曲线分析和双侧log-rank检验。P<0.05为差异具有统计学意义。
2 结 果
2.1 临床资料
共收集61例MM患者,其中22例NDMM患者、23例TRMM患者和16例PRMM患者,患者均为60岁左右的中老年人,免疫球蛋白(immunoglobulin,Ig)A型和IgG型MM占比较高,临床数据见
表1。TRMM和PRMM组化疗方案及疗效见
表2。
2.2 3组MM患者血清代谢物分析
代谢组学结果显示:NDMM组与TRMM组比较,有46个差异代谢物上调,34个差异代谢物下调(
图1);NDMM组与PRMM组比较,有37个差异代谢物上调,37个差异代谢物下调(
图2);TRMM组与PRMM组比较,有21个差异代谢物上调,24个差异代谢物下调(
图3)。
整体来看,3组患者的代谢物差异较显著,共筛选到3组总体差异代谢物70个,主要富集在氨基酸代谢通路,其中在NDMM组上调的代谢物,在PRMM和TRMM组表现为下调,且两治疗组代谢物水平较一致(附
图1,
https://doi.org/10.57760/sciencedb. xbyxb.00103)。
2.3 不同疗效NDMM患者血清代谢物分析
因药物治疗对患者体内代谢物产生了较大影响,故查询了22例NDMM患者的后续疗效,追踪到16例患者,其中12例疗效为完全缓解(complete response,CR)或VGPR,归为初诊_治疗有效组(newly diagnosed_treatment response,ND_TR);4例疗效未达VGPR,归为初诊_原发难治组(newly diagnosed_primary refractory,ND_PR)。2组间筛选出23个差异代谢物,在ND_PR组中,柠檬酸等17个代谢物下调,色氨酸等6个代谢物上调(
图4A)。差异代谢物主要富集在氨基酸代谢和癌症中心碳代谢等通路(
图4B)。上述结果说明,不同疗效MM患者体内的氨基酸代谢可能存在显著差异。然而,由于本研究收集的NDMM患者数量较少,代谢组学结果不适合构建模型。因此在GEO数据库下载了转录组测序数据用于模型构建。
2.4 差异表达AAMGs筛选及功能注释
从数据库下载了NDMM患者转录组测序数据。GSE55145和GSE68871的患者分为TRMM组(108例)和PRMM组(77例)。去除批次效应后,两数据集样本均匀分布。最后筛选了2 749个DEGs,数据库获取了1 745个AAMGs。DEGs与AAMGs取交集后,得到了不同疗效患者间差异表达的232个AAMGs,其中PRMM组有86个基因上调,146个基因下调。232个基因参与的生物过程(biological process,BP)包括氨基酸代谢和嘌呤核苷酸代谢等过程;细胞组分(cellular component,CC)层面主要定位在线粒体;分子功能(molecular function,MF)层面富集在糖基转移酶活性和磷酸酯水解酶活性等功能层面。此外,232个基因主要富集在嘌呤代谢、精氨酸和脯氨酸代谢及氧化磷酸化等通路(附
图2,
https://doi.org/10.57760/sciencedb. xbyxb.00103)。
2.5 预后分析
在这232个基因中,117个基因与预后显著相关(
P<0.05),其中76个为预后高风险基因(
HR>1,
P<0.05),基因的高表达与预后较差相关;41个为预后低风险基因(
HR<1,
P<0.05),基因的低表达与预后较好相关。基于K-M生存分析的结果,绘制了前35个基因的预后网络图,包括27个预后高风险基因和8个预后低风险基因(
图5A)。根据基因表达水平将患者分为高表达组和低表达组,图
5B~
5D展示其中3个预后相关基因的生存曲线。
2.6 模型构建及验证
基于117个AAMGs,采用LASSO回归筛选了48个预测疗效的重要特征基因(图
6A、
6B)。后续基于48个基因构建疗效预测模型。模型残差分析显示:SVM、RF和XGB 3种模型表现出较好的性能,GLM模型则相对较差(图
6C、
6D);且SVM、RF和XGB模型ROC的曲线下面积(area under the curve,AUC)均大于0.8,GLM模型的AUC值相对较小,为0.714(
图6E)。在外部验证中,XGB模型表现出更优秀的预测性能,3个外部验证队列的AUC值分别为0.815、0.708和0.709(图
7D~
7F)。因此最终选择XGB算法构建的预测模型。
用合并后的建模数据集(GSE68871和GSE55145)构建XGB模型列线图(包含基因
SMPD3、
PDE5A、
PDHB、
ITPKB、
YOD1)(
图7A)。列线图使用基因表达量计算预测NDMM患者为PRMM患者的风险,显示模型具有较好的预测性和准确性(图
7B、
7C)。外部队列的ROC曲线验证了模型的外推性(图
7D~
7F)。其中复发队列(GSE9782)的AUC为0.709(
图7F),表明模型还可有效预测复发患者对含BTZ化疗方案的疗效。
2.7 模型基因的表达和预后
2.8 模型基因富集分析和免疫分析
在
ITPKB高表达组中,心肌收缩通路和Hedgehog信号通路等活性较高;而在低表达组中,磷脂酰肌醇信号系统通路和谷胱甘肽代谢通路等活性较高。在
PDE5A高表达组中,牛磺酸和次牛磺酸代谢及RNA聚合酶通路活性较高;而在低表达组中,内吞作用通路及体内多种物质如脂肪酸和氨基糖等代谢通路的活性较高(附
图5,
https://doi.org/10.57760/sciencedb. xbyxb.00103)。在
PDHB高表达组中,NOTCH信号通路和MAPK信号通路活性较高;在低表达组中,氧化物酶体和柠檬酸循环等代谢通路活性较高。在
SMPD3高表达组中,自噬调控、泛素介导的蛋白质水解通路和p53信号通路活性较高;在低表达组中,亚油酸和亚麻酸代谢等通路活性较高。而在
YOD1高表达组中,药物代谢细胞色素p450、酪氨酸和α-亚麻酸代谢通路活性较高;在低表达组中,RIG-I样受体信号通路、Toll样受体信号通路和p53信号通路活性较高(附
图5,
https://doi.org/10.57760/sciencedb.xbyxb. 00103)。
免疫功能评分结果(附
图6,
https://doi.org/10. 57760/sciencedb.xbyxb.00103)显示:
ITPKB和
PDE5A高表达患者及
PDHB和
YOD1低表达患者的免疫功能评分较高。这表明这些患者的肿瘤微环境中可能存在更强的抗肿瘤免疫反应。此外,
ITPKB、
PDE5A高表达患者和
PDHB、
YOD1低表达患者的预后较好。这在一定程度上表明较强的抗肿瘤免疫反应可能给这部分患者带来了更好的预后。
3 讨 论
近年来越来越多的研究开始关注氨基酸代谢在癌症代谢中的作用
[9]。有研究
[10-12]指出氨基酸代谢参与了MM的发生、发展和耐药形成过程,氨基酸代谢研究可能为MM治疗和管理提供新的思路。本研究结果显示色氨酸在ND_PR组中上调,且多个色氨酸代谢产物水平存在差异,支持上述观点。色氨酸及其代谢产物吲哚在免疫调节中具有重要价值,主要代谢途径的抑制剂已用于癌症治疗
[13]。在MM患者中,色氨酸代谢异常与来那度胺治疗的不良结局相关
[14]。骨髓微环境中功能失调的浆细胞样树突状细胞与MM细胞的相互作用可上调色氨酸分解代谢途径的酶,导致免疫抑制和MM细胞增殖
[15]。此外,柠檬酸在ND_PR组中下调,且所筛差异代谢物中含有谷氨酸代谢产物,差异代谢物富集的通路包括谷氨酸代谢通路。柠檬酸是三羧酸循环的中间体,代谢酶的改变会导致循环代谢物失调,可能与肿瘤发生有关
[16]。由于肿瘤的代谢重排,许多肿瘤细胞对谷氨酰胺的依赖增加,以满足三羧酸循环,获取增殖生长所需的能量
[17]。一旦剥夺谷氨酰胺,肿瘤细胞会停止生长甚至死亡。在MM中,恶性浆细胞进入骨髓微环境会导致成骨细胞分化严重失衡,MM细胞会消耗大量谷氨酰胺
[18]。色氨酸代谢和谷氨酰胺代谢可能是MM治疗的重要靶点。
研究认为色氨酸代谢是癌症、神经退行性变性等疾病的治疗靶点
[14, 19];并且色氨酸及其代谢产物吲哚在免疫调节中具有重要价值,其主要代谢酶吲哚胺 2,3-双加氧酶1(indoleamine 2,3-dioxygenase 1,IDO1)抑制剂已用于癌症治疗
[13]。色氨酸的分解代谢途径在MM的治疗中也具有重要价值。在MM患者中,色氨酸代谢异常与来那度胺治疗的不良结局相关
[14]。骨髓微环境中功能失调的浆细胞样树突状细胞(plasmacytoid dendritic cells,pDCs)与MM细胞的相互作用可诱导犬尿氨酸途径(色氨酸分解代谢途径)的犬尿氨酸3-单加氧酶(kynurenine-3-monooxygenase,KMO)上调,导致免疫抑制和MM细胞的增殖
[15]。此外,在pDC-T细胞-NK细胞-MM细胞共培养模型中,KMO的抑制剂Ro61-8048处理可以激活pDCs,并触发MM特异性细胞毒性T淋巴细胞(cytotoxic T lymphocytes,CTL)和自然杀伤(natural killer,NK)细胞对肿瘤细胞的杀伤活性,单独使用抑制剂和与抗程序性死亡受体1(programmed death receptor 1,PD-1)抗体联合使用都可以增强MM细胞毒性并恢复抗MM免疫
[15]。另有研究
[20]发现谷氨酰胺代谢上调与浆细胞骨髓瘤(plasma cell myeloma,PCM)的进展和PIs耐药相关,PCM细胞系的生长和生存依赖于所摄取的谷氨酰胺。敏感PCM细胞系和耐药PCM细胞系的细胞膜上均表达谷氨酰胺转运体ASCT2,并且都对ASCT2抑制剂V9302敏感,V9302可协同增强PIs的细胞毒性。目前,研究
[21]认为谷氨酰胺代谢抑制剂的利用是非常有前景的癌症治疗策略。肿瘤细胞及肿瘤微环境的氨基酸代谢是MM等多种癌症发生、发展和耐药的关键一环,对氨基酸的利用及其代谢的调控可能是重要的癌症治疗新策略。
本研究构建的疗效预测模型具有较好的预测性和外推性。目前MM疗效预测的研究主要是基于转录组测序数据,预测患者疗效或细化与疗效相关的分子特征。有研究
[22-23]基于NDMM患者的转录组测序数据构建了预测PAD/VCD化疗方案疗效的机器学习模型。另有研究
[24]使用单细胞测序鉴定DEGs,结合临床变量和DEGs开发预测BTZ治疗反应的联合模型。这些研究可在一定程度上助力MM治疗,为耐药机制研究提供思路。
磷脂酰肌醇3激酶/蛋白激酶B(phosphoinositide 3-kinase/glucocorticoid regulated kinases,PI3K-AKT)信号通路是一种细胞内信号转导途径。肌醇-1,4,5-三磷酸激酶B(inositol-trisphosphate 3-kinase B,ITPKB)在这一通路中发挥重要作用。ITPKB可将肌醇-1,4,5-三磷酸转换为可溶性拮抗剂肌醇-1,3,4,5-四磷酸,进一步抑制通路的信号转导
[25]。ITPKB诱导产生的 肌醇-1,3,4,5-四磷酸可抑制B细胞、T细胞和中性粒细胞的钙离子摄取,调控细胞的存活和功能
[26]。调节能量代谢的丙酮酸脱氢酶β(pyruvate dehydrogenase beta,PDHB)的基因位于线粒体,限制谷氨酰胺饮食与抑制丙酮酸代谢相结合可有效治疗肝细胞癌,敲除丙酮酸脱氢酶α、PDHB和丙酮酸羧化酶均可诱导三羧酸循环的代谢重编程,破坏线粒体功能,抑制谷氨酰胺耗竭适应下的肿瘤细胞增殖
[27]。此外,研究
[28]揭示PDHB的表达与多个免疫细胞的浸润和免疫微环境的状态密切相关,PDHB可能是预测癌症患者免疫反应的潜在分子标志物。本研究同样发现:在MM患者中,PDHB低表达患者的免疫评分更高,且预后更佳。这表明PDHB的低表达可能与更活跃的免疫微环境相关,具有更强的抗肿瘤免疫反应。YOD1是一种去泛素化酶,其表达与维持蛋白质稳态和介导蛋白质降解相关。在急性髓细胞性白血病(acute myeloid leukemia,AML)中,YOD1的表达显著降低,可诱导肿瘤抑制因子p53的降解,进而抑制AML细胞的凋亡
[29]。本研究发现YOD1在PRMM队列中低表达,其低表达与不良预后相关,并且YOD1低表达队列中p53信号通路活性高。因此,在MM中可能存在与AML类似的机制,即YOD1低表达所导致的预后不良可能与肿瘤抑制因子p53的降解有关。
本研究存在以下局限性:1)由于数据限制,未分析临床特征与疗效的相关性。2)所筛模型基因和通路未做分子实验验证。3)未使用临床标本及数据对模型性能做验证。4)由于代谢组学和转录组学的检测样本非同一批患者的纵向检测样本,故未做联合分析。
总之,本研究充分利用代谢组学、转录组学和生物信息学分析构建的疗效预测模型可助力临床医师对PRMM患者的管理,相关分析也为探讨耐药机制,挖掘新的治疗靶点提供了理论支持。氨基酸代谢是MM发生、发展和耐药的关键一环,对氨基酸的利用及其代谢的调控可能是重要的癌症治疗新策略。
国家自然科学基金(81870166)