基于单细胞分化轨迹的胶质母细胞瘤预后评估模型构建及验证

宁静源，范小晴，刘发强，郭玉梅

基于单细胞分化轨迹的胶质母细胞瘤预后评估模型构建及验证

宁静源1，2，3，范小晴3，刘发强1，2，郭玉梅1，2

1 石家庄市疾病预防控制中心微生物检验所，石家庄 050000；2 河北省疑难细菌研究重点实验室；3 河北医科大学免疫教研室

基于单细胞分化轨迹筛选胶质母细胞瘤（GBM）细胞分化相关基因（DRG），构建预后模型并加以验证。通过基因表达综合数据库中18例GBM样本的单细胞测序数据分析细胞分化轨迹，筛选DRG，并进行富集分析；结合癌症基因组图谱数据库中169例GBM患者的生存时间和生存状态，用单因素Cox回归分析和多因素Cox回归分析对DRG进行筛选并建立预后模型和动态列线图；对模型进行验证和内部特征分析。GBM细胞分化轨迹中存在7种不同状态的细胞亚群，在细胞发育初期，主要的细胞为间充质干细胞（MSC），中期部分MSC停止分化，部分MSC分化为神经元细胞（Neurons），晚期主要存在星形胶质细胞（Astrosyte）和部分Neurons细胞。对7组不同状态的细胞进行组间差异分析后共获得了518个DRG，京都基因与基因组百科全书（KEGG）富集分析结果显示，DRG主要富集在局部粘连、冠状病毒感染-COVID-19、p53信号传导途径、核糖体、细胞衰老、ECM与受体的相互作用等信号通路；基因本体富集结果显示，DRG主要富集在伤口愈合、轴突发育、细胞外基质、含胶原蛋白、焦点黏附、结构成分等；疾病本体富集分析结果显示，DRG与周围神经系统、自主神经系统、恶性胶质瘤、神经母细胞瘤、星形细胞瘤、多形性GBM密切相关。用4个DRG（SPP1、IFI44L、CD81、AEBP1）建立预后模型，经过验证该预后模型具有良好的准确性。根据风险分数的中位值将所有患者分为高风险组和低风险组，与低风险组相比，高风险组的总生存期短，高风险组和低风险组在免疫微环境中存在着大量差异。根据DRG的表达值和患者临床特征建立的动态列线图具有极强的预测性能和临床实用性。GBM细胞的复杂分化过程影响患者预后，DRG相关的预后模型和动态列线图可以准确预测GBM患者的预后。

胶质母细胞瘤；预后；列线图；单细胞分化轨迹

胶质母细胞瘤（GBM）是常见的侵袭性神经系统恶性肿瘤之一，在成人中具有极高的致死率［1］。因GBM复杂的肿瘤异质性，其治疗仍然非常困难［2-3］。常规放化疗并不能降低GBM的复发率［4-5］。免疫疗法使GBM患者的预后状态得到一定程度改善［6］，但并非所有患者都能从免疫治疗中获益［7-9］。目前对于GBM免疫疗法的研究大多处于Ⅰ、Ⅱ期实验，必须经过Ⅲ期实验才能确定该疗法的价值。GBM的异质性与肿瘤的发生、发展和复发密切相关［11］。准确揭示GBM的异质性，从而制定个性化的治疗方案已经成为精准医学关注的焦点［12］。混合RNA测序在细胞研究中取得了一定成果，但并不能揭示单个细胞的表达模式。随着单细胞技术的发展，这使得在单个细胞水平揭示更加复杂全面的表观基因组成为了可能。本研究用单细胞测序数据揭示GBM细胞分化轨迹，同时结合混合RNA测序数据进行分析，构建预后模型和动态列线图并加以验证，为GBM患者的诊断和治疗提供潜在靶点。

1 资料与方法

1.1临床资料从基因表达综合数据库（GEO，https：//www.ncbi.nlm.nih.gov/geo/）中的GSE159416中获得18例GBM患者的单细胞测序RNA数据，从癌症基因组图谱数据库中（TCGA，https：//portal.gdc.cancer. gov/）获得169例GBM患者的每百万映射读数的转录本片段数据。在R语言环境中，将FPKM值转化为每百万转录本。

1.2单细胞数据预处理和细胞分化轨迹分析在R语言环境，应用“Seurat R”包（版本 3.2.1）来预处理单细胞测序数据。首先通过PercentageFeatureSet函数计算线粒体基因的百分比，并通过相关性分析阐明测序深度与线粒体基因序列和（或）总细胞内序列之间的关系。然后使用带默认参数的NormalizeData函数对单细胞表达式矩阵进行归一化，FindVariableFeatures函数和默认的“vst”方法将所有基因进行排序，选择排序后的前1 500个基因作为高变基因。使用FindClusters函数识别不同的细胞簇，分辨率为0.5。使用RunUMAP函数可视化低维度细胞景观。预处理结束以后应用“monocle”包的标准流程进行细胞分化轨迹分析。

1.3富集分析在R语言环境中，使用“clusterProfiler”包中的enrichKEGG和enrichGO函数进行富集分析，用“ggplot2”包、“ggpubr”包内置绘图函数对富集结果进行可视化。

1.4预后模型构建通过单因素Cox回归分析，获得与患者总生存期显著相关的分化相关基因（DRG），使用“glmnet”包进行最小绝对收缩和选择算子（LASSO）分析，得到最佳DRG。通过多因素Cox回归分析建立基于最佳DRG的预后模型，并基于以下公式计算所有患者的风险分数：

Coefi代表通过多因素Cox回归分析得到的每个DRG的系数值，xi代表选定的每个DRG的表达值。

1.5动态列线图生成在R语言环境中，预先载入该过程所需的依赖包（“shiny”包、“plotly”包、“compare”包、“stargazer”包），根据构建预后模型产生的多因素Cox回归分析结果，使用“DynNom”包中的DynNom函数生成动态列线图，用Dnbuilder函数导出本地动态列线图脚本。最后将导出的脚本文件上传至shinyapp网站（https：//www.shinyapps.io/），为动态列线图建立网页连接。

1.6统计学方法所有数据均在R4.1.2相应的软件包和Strawberry Perl软件中进行分析。<0.05为差异有统计学意义。

2 结果

2.1单细胞分化轨迹分析结果对GSE159416中18例GBM患者的单细胞RNA测序数据进行质量控制和标准化，共纳入4 073个细胞。测序深度与线粒体基因序列间不存在相关性，与细胞内总序列存在正相关（=0.92）。共分析了15 034个基因，其中13 534个基因在细胞间的变异度较低，1 500个基因被确定为高变异基因。4 073个细胞聚集为16个簇，其中0、1、2、3、4、8、10、13簇均为间充质干细胞（MSC），5、11、12、9为星形胶质细胞（Astrocyte），6、7、14为神经元细胞（Neurons），15为巨噬细胞（Macrophage）细胞。每种细胞的标志基因见图1A。拟时序分析结果显示，GBM细胞分化过程中存在7种不同状态。状态1和7主要由Astrocyte、Neurons组成，状态2和3主要由Neurons组成，状态4、5、6主要由MSC组成，见图1B。

注：A：每种细胞类型的前10个标志基因表达的热图；B：拟时序分析。

2.2DRG及其富集分析结果对7组不同状态的细胞进行组间差异分析（|LogFC|>1，FDR<0.05）后共获得了518个DRG。京都基因与基因组百科全书（KEGG）富集分析结果显示，DRG主要富集在局部粘连、冠状病毒感染-COVID-19、p53信号传导途径、核糖体、细胞衰老、ECM与受体的相互作用等信号通路。基因本体富集结果显示，DRG主要富集在伤口愈合、轴突发育、细胞外基质、含胶原蛋白、焦点黏附、结构成分等。疾病本体富集分析结果显示，DRG与周围神经系统、自主神经系统、恶性胶质瘤、神经母细胞瘤、星形细胞瘤、多形性GBM密切相关。

2.3预后模型构建及验证结果对DRG进行单因素Cox回归分析得到65个预后相关DRG（均<0.05），依次通过LASSO回归分析和多因素Cox回归分析筛选得到4个DRG，分别为分泌性磷蛋白1（SPP1）、扰素诱导样蛋白 44（IFI44L）、CD81抗原（CD81）、脂肪细胞增强子结合蛋白 1（AEBP1）。这4个DRG在不同细胞类型中的表达见图2A，SPP1在Macrophage和Astrocyte中高表达，IFI44L在Neurons中高表达，CD81在Macrophage中高表达，AEBP1在MSC中高表达。

用4个DRG建立预后模型并计算所有患者的风险分数，风险分数=（0.176×SPP1表达值）+（-0.242×IFI44L表达值）+（0.740×CD81表达值）+（0.472×AEBP1表达值）。根据风险分数的中位值将所有患者分为高风险组和低风险组。SPP1、CD81、AEBP1在高风险组中高表达，IFI44L在低风险组中高表达，见图2B。主成分分析（PCA）结果显示，参与风险分数计算的4个DRG区分患者的效果优于所有基因和所有DRG，见图2C。生存分析结果显示，高风险组的总生存期（OS）明显缩短（<0.05），见图2D。单因素Cox回归分析和多因素Cox回归分析证明，风险分数是影响患者预后的独立因素（均<0.05），见图2E、2F。在所有患者队列中，风险分数预测预后的接受者操作特征曲线（ROC）的曲线下面积为0.766，高于其他临床特征的曲线下面积，见图2G。另外，风险分数在1、3、5年的曲线下面积分别为0.766、0.711、0.824，见图3H。将所有患者随机分为训练集和验证集进行内部验证。结果与预期一致，训练集（图2I）和验证集（图2J）中高风险组OS均短于低风险组。

注：A：SPP1、IFI44L、CD81、AEBP1在单细胞测序数据不同细胞类型中的表达情况；B：SPP1、IFI44L、CD81、AEBP1在高低风险患者中的表达情况；C：PCA分析；D：生存分析；E：临床特征和风险分数的单因素Cox回归分析；F：临床特征和风险分数的多因素Cox回归分析；G：临床特征和风险分数的ROC曲线；H：风险分数在1、3、5年的ROC曲线；I：训练集生存曲线；J：验证集的生存曲线。

2.4预后模型特征基因集富集分析（GSEA）显示，高风险组细胞黏附分子、糖胺聚糖生物合成硫酸角质素、白细胞跨内皮迁移、补体和凝血级联、焦点黏附等通路被激活（图3A）。肿瘤免疫微环境分析结果表明，高风险组肿瘤纯度较低，其中基质细胞和免疫细胞含量均高于低风险组（图3B）。免疫细胞浸润结果显示，高风险组的DC细胞、B细胞、巨噬细胞、Treg细胞浸润较多（图3C）。免疫功能分析结果显示，高风险组抗原提呈能力较高，T细胞激活和Ⅱ型干扰素反应较为敏感（图3D）。免疫检查点结果显示高风险组中包括CD274在内的多个免疫检查点表达较高（图3E）。

2.5预测模型的临床验证和动态列线图的建立使用基因表达谱交互分析数据库（GEPIA，http：//gepia.cancer-pku.cn/index.html）中TCGA合并GTEx的转录组数据对参与模型构建的4个DRG（SPP1、IFI44L、CD81、AEBP1）进行了外部数据集RNA水平验证。结果显示，与健康组相比，肿瘤组中4个DRG表达上调（图4A）。通过人类蛋白质图谱（HPA，http：//www.proteinatlas.org/）数据库收集GBM患者和健康者临床免疫组织化学染色结果，与健康者相比，GBM患者4个DRG蛋白表达高（图4B）。用4个DRG和临床特征发布了交互式动态列线图（https：//medicalcalculator.shinyapps.io/GBMNomapp/）。该列线图界面如图所示（图4C），用户仅需输入患者基因表达值和相应的临床特征，点击Predict便可以快速对患者预后状态进行评估。矫正曲线证明该列线图具有良好的准确性（图4D）。

3 讨论

GBM的细胞分化过程是决定肿瘤恶性程度以及促进肿瘤侵袭、转移的最主要因素，肿瘤细胞的分化程度越高的GBM，恶性程度较轻，且不容易发生转移［13］；而分化程度越低的肿瘤细胞恶性程度越高，也越易发生侵袭转移，患者往往预后差［4］。充分了解GBM中的异质性对于研究新的治疗方法必不可少［3］。本研究在单细胞水平探讨了GBM细胞分化过程中的异质性。根据单个细胞转录组数据，确定了GBM细胞分化过程中存在的7种不同状态。在细胞发育初期，主要的细胞为MSC（状态5），中期部分MSC停止分化（状态6），部分细胞群分化为Neurons细胞（状态2、3和4），晚期主要存在Astrosyte细胞和部分Neurons细胞（状态1和7）。以上结果表明，GBM中MSC存在着不同的亚型，部分亚型具有转分化为Astrosyte细胞的能力

本研究发现，调控GBM细胞分化的差异基因共518个。富集分析结果显示，DRG与周围神经系统、自主神经系统、恶性胶质瘤、神经母细胞瘤、星形细胞瘤、多形性GBM密切相关，这进一步证明DRG对于患者疾病进展的重要性。另外大量的DRG富集在新型冠状病毒感染通路，该结果提示GBM的发生可能与新型冠状病毒肺炎之间存在着一定的串扰基因。GBM患者作为极为脆弱的群体，对于新型冠状病毒感染的抵抗能力似乎更差，而目前对于新型冠状病毒感染阳性的GBM患者并没有明确的治疗策略，分析GBM进展和COVID-19感染之间的标志性串扰基因可能是突破口。

结合GBM患者生存时间和生存状态，最终确定了4个DRG（SPP1、IFI44L、CD81、AEBP1），建立了预后预测模型并证明了模型的稳定性和有效性。KIJEWSKA等［14］的研究发现，SPP1主要参与免疫反应、组织重塑；在许多癌症中过表达，并通过促进肿瘤细胞迁移、侵袭和癌症干细胞自我更新来调节肿瘤进展。胶质瘤中SPP1过表达可以增加胶质瘤细胞的干性，使肿瘤细胞具有更低的分化程度，从而增加了肿瘤细胞的恶性程度。SPP1过表达的胶质瘤细胞更易发生侵袭、转移。临床研究验证了SPP1过表达的肿瘤细胞具有更高的恶性程度，SPP1过表达肿瘤患者具有更差的预后［15］。WEI等［16］研究发现，SPP1介导GBM相关巨噬细胞浸润。这也与本研究结果相似，SPP1的表达与GBM的免疫细胞浸润显著相关［17］。另外研究也表明，CD81高表达为GBM细胞提供了放射治疗抗性［18］。AEBP1也是GBM免疫疗法的潜在治疗靶点［19］。但是IFI44L在GBM中的作用目前并无相关研究。

本研究根据风险分数的中位值将患者分为高风险组和低风险组，进行模型特征分析。免疫细胞浸润结果显示，高风险组的DC细胞、B细胞、巨噬细胞、Treg细胞浸润较多。免疫功能分析结果显示，高风险组抗原提呈能力较高，T细胞激活和Ⅱ型干扰素反应较为敏感。免疫检查点结果显示，高风险组包括CD274在内的多个免疫检查点表达较高，提示高风险组可能更易从PD-L1治疗中获益。这与之前的研究［20-21］结果相似。再次说明两种不同分化程度的GBM具有不同的免疫浸润状态，这也可能是影响患者生存的因素之一。分化程度更低的GBM成免疫抑制状态，Treg细胞浸润较多，但同时CD274在内的多个免疫检查点表达较高。这也反映该类型的肿瘤虽然恶性程度高预后较差，但其免疫治疗效果可能更好。重要的是，用患者的临床特征和参与模型构建的4个DRG发布了动态列线图网页。据搜索，该研究是GBM患者中第一个通过多因素Cox回归分析构建的DRG相关的列线图，并且证明了该列线图具有良好的准确性。总之动态列线图使研究结果实现了更好的临床转换，为临床使用提供了更便利的诊断程序。

注：A：GSEA分析；B：肿瘤纯度分析；C：免疫细胞浸润分析；D：免疫细胞功能分析；E：免疫检查点分析。

注：A：GEPIA数据库中健康样本和肿瘤样本RNA表达对比；B：HPA数据库中健康样本和肿瘤样本蛋白表达对比；C：动态列线图网页界面；D：校正曲线。

总之，本研究在单细胞水平鉴定了GBM细胞分化过程中的特征基因，这些基因的表达变化决定了GBM的异质性程度，与肿瘤的侵袭转移有关，强调了GBM细胞分化在患者预后中的意义，为GBM治疗提供了方向。但本研究也存在不足，需要进一步评估该列线图的稳健性。另外，TCGA数据库GBM患者的临床特征变量并不完善，因此该列线图还具有一定的优化空间。

［1］ WIRSCHING H G， GALANIS E， WELLER M. Glioblastoma［J］. Handb Clin Neurol， 2016，134：381-397.

［2］ DAVIS M E. Glioblastoma： Overview of disease and treatment［J］. Clin J Oncol Nurs， 2016，20（5 Suppl）：2-8.

［3］ FRIEDMANN-MORVINSKI D. Glioblastoma heterogeneity and cancer cell plasticity［J］. Crit Rev Oncog， 2014，19（5）：327-336.

［4］ LAH T T， NOVAK M， BREZNIK B. Brain malignancies： Glioblastoma and brain metastases［J］. Semin Cancer Biol， 2020，60：262-273.

［5］ CAMPOS B， OLSEN L R， URUP T， et al. A Comprehensive profile of recurrent glioblastoma［J］. Oncogene， 2016，35（45）：5819-5825.

［6］ HUTÓCZKI G， VIRGA J， BIRKÓ Z， et al. Novel Concepts of glioblastoma therapy concerning its heterogeneity［J］. Int J Mol Sci， 2021，22（18）：10005.

［7］ LE RHUN E， PREUSSER M， ROTH P， et al. Molecular targeted therapy of glioblastoma［J］. Cancer Treat Rev， 2019，80：101896.

［8］ DAUBON T， HEMADOU A， ROMERO GARMENDIA I， et al. Glioblastoma Immune landscape and the potential of new immunotherapies［J］. Front Immunol， 2020，11：585616.

［9］ RIVERA M， BANDER E D， CISSE B. Perspectives on Microglia-Based Immune Therapies Against Glioblastoma［J］. World Neurosurg， 2021，154：228-231.

［10］ MEDIKONDA R， DUNN G， RAHMAN M， et al. A review of glioblastoma immunotherapy［J］. J Neurooncol， 2021，151（1）：41-53.

［11］ DECORDOVA S， SHASTRI A， TSOLAKI A G， et al. Molecular heterogeneity and immunosuppressive microenvironment in glioblastoma［J］. Front Immunol， 2020，11：1402.

［12］ BROWN D V， STYLLI S S， KAYE A H， et al. Multilayered Heterogeneity of Glioblastoma Stem Cells： Biological and Clinical Significance［J］. Adv Exp Med Biol， 2019，1139：1-21.

［13］ SENGA S S， GROSE R P. Hallmarks of cancer-the new testament［J］. Open Biol， 2021，11（1）：200358.

［14］ KIJEWSKA M， KOCYK M， KLOSS M， et al. The embryonic type of SPP1 Transcriptional regulation is reactivated in glioblastoma［J］. Oncotarget， 2017， 8（10）：16340-16355.

［15］ CHEN P， ZHAO D， LI J， et al. Symbiotic macrophage-glioma cell interactions reveal synthetic lethality in PTEN-null glioma［J］. Cancer Cell， 2019，35（6）：868-884，e6.

［16］ WEI J， MARISETTY A， SCHRAND B， et al. Osteopontin mediates glioblastoma-associated macrophage infiltration and is a potential therapeutic target［J］. J Clin Invest， 2019，129（1）：137-149.

［17］ SZULZEWSKY F， PELZ A， FENG X. Glioma-associated microglia/macrophages display an expression profile different from M1 and M2 polarization and highly express Gpnmb and Spp1［J］. PLoS One， 2015，10（2）：0116644.

［18］ ZHENG W， CHEN Q， LIU H， et al. CD81 enhances radioresistance of glioblastoma by promoting nuclear translocation of Rad51［J］. Cancers， 2021，13（9）：1998.

［19］ LIU M， YU Y， ZHANG Z， et al. AEBP1 as a potential immune-related prognostic biomarker in glioblastoma： A bioinformatic analyses［J］. Ann Transl Med， 2021，9（22）：1657.

［20］ CHEN Z， HAMBARDZUMYAN D. Immune microenvironment in glioblastoma subtypes［J］. Front Immunol， 2018，9：1004.

［21］ DE FELICE F， PRANNO N， MARAMPON F， et al. Immune check-point in glioblastoma multiforme［J］. Crit Rev Oncol Hematol， 2019，138：60-69.

（2022-06-17）

郭玉梅（E-mail： 325069675@qq.com）

10.3969/j.issn.1002-266X.2022.28.015

R739.4

A

1002-266X（2022）28-0063-06