单细胞分析如何过渡到预后模型？
导语
今天和大家分享的是2020年10月份发表在Aging-US杂志上的一篇文章(IF=4.831)“Glioblastoma cell differentiation trajectory predicts the immunotherapy response and overall survival of patients”。文章中作者根据GBM样品的整合单细胞RNA测序数据进行分析，通过轨迹分析找到了与分化相关的基因，继而通过TCGA数据库以这些基因为基础分别构建了分子分型以及预后预测模型。
Glioblastoma cell differentiation trajectory predicts the immunotherapy response and overall survival of patients
胶质母细胞瘤分化轨迹可预测患者的免疫治疗反应和总体存活率
一、研究背景
胶质母细胞瘤（GBM）是最常见的原发性恶性脑肿瘤，占所有肿瘤的14.6％和中枢神经系统（CNS）的恶性肿瘤的48.3％。由于GBM具有高度侵略性，GBM患者在数天或数月内会迅速发生神经功能缺损，癫痫发作和颅内高压症状。尽管在包括手术切除，放疗和化学疗法在内的各种治疗方法的开发方面取得了长足的进步，但在过去的几十年中，在延长生存期和改善预后方面进展甚微。癌症干细胞（CSC）的分化过程中在肿瘤微环境影响癌细胞多种因素，导致异源细胞分化状态和细胞命运。单细胞转录组学分析已成为一种强大的方法，可通过以显微镜分辨率同时研究整个肿瘤样品的基因组的全面性质，来提供表征细胞状态及其转变的机会。将这种综合的组成肿瘤的细胞排列成轨迹，有助于根据分化状态了解肿瘤细胞亚群，并揭示伴随细胞命运规范的遗传级联和相关致瘤途径。然而，尚不清楚GBM细胞是否处于不同的分化状态，以及基于细胞分化轨迹的GBM患者新分类是否与肿瘤生物学行为相关，并在预测患者存活率和免疫治疗反应中发挥作用。
二、研究思路
三、结果解读
1、使用scRNA-seq数据确定13个细胞簇
根据质量控制标准和GBM scRNA-seq数据的标准化，排除了194个低质量细胞，来自GBM核心的2149个细胞纳入分析。检测到的基因数量与测序深度显着相关。总共包括19,752个对应基因，方差分析显示1,500个高度变异基因。采用主成分分析(PCA)确定可用维度并筛选相关基因。但是PCA结果并未证明GBM中细胞间存在明显的分离。选择了20个P值<0.05主要成分（PCs）以进行后续分析。
图1. 基于单细胞RNA-seq数据的13个细胞簇
应用tSNE算法，并将人类GBM中的细胞成功分类为13个独立簇。差异表达分析共鉴定自13个聚类的8,025个标记基因。根据标记基因的表达模式，这些簇被标记为singleR和CellMarker。包含518个单元的簇0注释为GBM CSC；包含878个细胞的簇1、2、6和10标注为GBM癌细胞或GBM细胞；包含196个细胞的簇3标注为星形胶质细胞；包含44个细胞的簇11被标注为少突胶质细胞。包含319个细胞的簇4、5和9被标注为肿瘤相关的巨噬细胞；包含77个细胞的簇8被注释为典型的M1巨噬细胞；包含81个细胞的簇7被注释为典型的M2巨噬细胞；将包含36个单元的簇12标记为T细胞。
2、将GBM细胞可分为两个不同分化亚群，其GDRG与免疫和代谢途径相关
通过轨迹分析将GBMs中的所有细胞投射到一个根和两个分支上。结果表明，GBM CSCs主要位于根中，而GBM细胞主要位于分支中。有趣的是，分支I中的GBM细胞（这里定义为I型GBM细胞）全部来自簇2、6和10（434个细胞），而分支II中的GBM癌细胞（这里定义为II型GBM细胞）全部来自簇1（444个单元）。I型和II型GBM亚群中细胞的分化程度差异很大。将I型和II型GBM细胞的分支依赖性标记基因为GDRGs。进行差异分析，并将265个标记基因鉴定为I型GDRG，将193个标记基因鉴定为II型GDRG。因此，GBM中最终鉴定出总共498个GDRG。
图2. 不同分化模式的细胞亚群的GO注释、轨迹分析和GSEA
进行基因集富集分析（GSEA），以鉴定具有不同分化模式的GBM细胞的相关分子机制和途径。结果显示，I型GDRGs与免疫过程的调节呈负相关，例如抗原加工和免疫细胞分化，II型GDRGs与代谢相关途径（如碳代谢）呈显着正相关（氨基酸的生物合成，糖酵解和糖异生）。GDRGs的这些发现表明，处于不同分化状态的GBM细胞表现出独特的肿瘤生物学特征，这可能为GBM的分子标记（包括内在特性和相关通路的调控）提供新的证据。
然后确定是否可以使用大量RNA-seq数据识别观察到的GBM细胞亚群。如在所示图3A - 4C，相关性分析表明，无论是TCGA数据库还是CGGA数据库，I型GDRGs的scRNA-seq数据与bulk RNA-seq数据高度相关II型GDRGs。这些发现表明，I型和II型GBM细胞也可以通过大量RNA-seq数据通过GDRG的表达来识别，因为这些高度相关的基因可能表明一个共同的细胞起源。
图3. 两种GDRG亚型的相关性分析和体细胞突变分析
3、两个GBM细胞亚群的GDRG在功能上相关且大部分发生突变
为确定来自不同GBM细胞亚群的基因谱是否在功能上相关，作者利用metagens来代表相应基因谱的总体表达模式。I型和II型元基因(分别由I型和II型gdrg组成)的表达量是由组成基因表达的加权平均数得出的。结果显示I型和II型metagene表达之间之存在很强的相关性，并在scRNA-seq和RNA-seq数据II型基因表达谱表达。表明I型和II型GDRGs功能相关,功能上相关的具有不同区分模式的不同GBM子集。
作者还分析了TCGA队列中GDRG的体细胞突变状态。大多数GDRG（90.8％，452/498）存在突变。表皮生长因子受体（EGFR）表现出最高的突变频率（53％），其次是CDK4（16％）和TSPAN31（16％）。I型GDRGs中有246个基因（92.8％）具有突变，II型GDRGs中有262个基因（89.4％）具有突变。两组之间的突变频率无统计学意义。但是，排名前9位的GDRG中有8个是II型GDRG，而只有1个是I型GDRG。这些结果表明，GDRGs的突变状态具有高度的异质性，提示GDRGs在GBMs的发生发展过程中发挥着关键作用。
4、基于GDRG的GBM患者分类与不同的OS结局和临床病理特征相关
为建立基于GDRGs的表达模式的GBM分类，对来自TCGA数据库的151名GBM患者进行了基于机器学习的无监督共识聚类。根据累积分布函数（CDF）曲线和共识热图下面积的相对变化，确定最佳簇数为两个（k值= 2），并且在CDF下面积未观察到明显的增加。曲线（图5A – 5C）。因此，所有GBMs分为两组：分子簇1 (MC1)的80例(53.0%)和分子簇2 (MC2)的71例(47.0%)。Kaplan-Meier生存分析表明，MC1 GBM的患者的OS显着低于MC2 GBM的患者。
图4.GBM患者基于GDRG的分类的识别和验证
之后在CGGA队列中验证了GDRG的分类。如在所示图5G - 5I，也被确定的簇的最佳数目是两个（k值= 2），并且将患者还归类为MC1（265例，75.7％）和MC2（85例，占24.3％）。Kaplan-Meier生存分析还表明，与MC2 GBM相比，MC1 GBM的患者的OS较差。
此外比较了TCGA队列中患者的两个MC之间GDRGs的表达模式和临床病理特征。I型GDRG基因表达谱的表达水平显著升高，II型GDRG基因表达谱的表达水平MC1 GBM的比MC2 GBM患者显著降低。在CGGA验证队列中也观察到了相同的发现。因此假设MC1患者主要对应于I型GBM子集的功能特性，而MC2患者主要对应于II型GBM子集的功能特性。
MC1型GBMs患者中I型GDRG metagene的表达水平显著高于MC2型GBMs患者， II型GDRG metagene的表达水平显著低于MC2型GBMs患者。在CGGA验证队列中也发现了同样的结果。因此，因此假设MC1患者主要符合I型GBM亚群的功能特性，而MC2患者主要符合II型GBM亚群的功能特性。
表1.TCGA训练队列和CGGA验证队列的MC1和MC2亚组GBM患者的基线资料
表1显示了MC1和MC2 GBMs患者的人口学特征和临床病理特征。与MC2患者相比，TCGA队列和CGGA队列中的MC1患者均明显年轻(P=0.007)。然而，患者两种MCs在其他变量上无显著差异(P均> 0.05)。综上所述，上述结果表明，基于GDRG的GBM患者分类在不同人群中是稳健和可靠的，并且可以根据该分类明确区分不同的生存结局。
5、基于GDRG的GBM患者分类与免疫检查点的表达模式与免疫治疗反应的不同可能性相关
确定了6个主要免疫检查点的表达，分别是PDCD1（PD1），CD274（PDL1），PDCD1LG2（PDL2），CTLA4，CD80和CD86，并比较了GBM患者的两个GBM细胞亚群和两个MC 。根据scRNA-seq数据，PD1，PDL1和PDL2在I型GBM细胞中高表达，而CTLA4，CD80和CD86在II型GBM细胞中高表达。就大量RNA序列数据而言，在TCGA数据库和CGGA数据库中，MC1 GBM患者中PD1，PDL1和PDL2高表达，而MC2 GBM患者中CTLA4，CD80和CD86高表达。
图5.GBM患者的免疫治疗反应预测
之后使用肿瘤免疫功能障碍和排斥（TIDE）算法来预测免疫治疗反应的可能性。根据TCGA训练队列的结果，MC2 GBM患者（43.7％，31/71）比MC1患者（20.0％，16/80，P = 0.003）更有可能对免疫疗法产生反应。同样，在CGGA验证队列中，MC2 GBM患者比MC1患者更有可能对免疫疗法作出反应。
然后，进行SubMap分析以预测在GBM患者的两个MC中对PD1和CTLA4抑制剂产生临床反应的可能性。TCGA和CGGA队列中的MC1 GBM患者对抗PD1治疗更敏感，而MC2 GBM患者对抗CTLA4治疗的更敏感。
6、FN1，APOE，RPL7A和GSTM2是人类GBM中4种最显著的的预测生存GDRG
进行单变量Cox分析并在TCGA训练集中确定了45个与预后相关的GDRG。然后进行最小绝对收缩和选择算子（LASSO），然后进行多变量Cox分析，并鉴定出4个显着预测生存的GDRG：纤连蛋白1，载脂蛋白E，核糖体蛋白L7a和谷胱甘肽S-转移酶mu 2。FN1在GBM细胞和T细胞中显着上调（一般）；在GBM细胞，CSC和巨噬细胞中，APOE明显上调；RPL7A在GBM细胞，CSC和巨噬细胞中显着上调；而GSTM2在GBM细胞和星形胶质细胞中显着上调。使用包括163个GBM（TCGA）样品和207个正常（GTEx）样品的基因表达谱交互式分析（GEPIA）数据库验证了这4个GDRG的表达。发现与正常组织相比，GBM中所有4种可预测生存的GDRGs均上调。
图6.预测模型中四个GDRG的表达和生存分析
Kaplan-Meier生存分析表明，GBM患者中FN1的高表达和APOE，GSTM2和RPL7A的低表达与GBM患者的OS差有关。
7、以GDRG的预后风险评分模型的生成和验证
基于上述4个GDRG，使用以下公式开发了预后风险评分模型：风险评分= Exp FN1 ×1.66 + Exp APOE ×（-0.93）+ Exp RPL7A ×（-1.30）+ Exp GSTM2 ×（-0.90）。计算TCGA训练集中所有患者的风险评分，并使用风险评分的中位数作为临界值，将患者分为高风险（高评分）组或低风险（低评分）组值。Kaplan-Meier生存分析表明，高危组患者的OS显着低于低危组。用于OS预测的GDRG签名的C指数为0.781。随时间变化的接收器工作特性（ROC）分析还表明，GDRG签名在预测0.5年，1年，2年和3年OS率方面表现出卓越的性能，曲线下面积（AUC）值为0.767， 0.712、0.752和0.776。
图7. 基于GDRG的风险评分模型的生存分析和风险评分
然后使用CGGA队列以类似方式验证预测公式，所有GBM患者均分为高危组或低危组。与TCGA训练集的结果一致，生存分析也显示高风险组患者的OS显着低于低风险组患者。GDRG验证的C指数为0.715，基于时间的ROC分析还提出了在CGGA验证集中预测OS的有良好价值。这些结果表明，基于GDRG的预后风险评分模型可以作为GBM患者不同人群的可靠预后指标。
8、开发验证的GDRG和临床病理参数可用于预后诺模图
研究GDRG签名的预后意义是否独立于其他临床病理变量来预测GBM患者的生存，进行了单变量和多变量Cox回归分析，结果表明GDRG验证在两个患者中均独立于OS。TCGA和CGGA队列（表2）。
表2.TCGA GBM培训队列和CGGA GBM验证队列中临床病理变量和GDRG验证的单变量和多变量Cox比例风险分析
最后成功开发了预后列线图，为个体OS预测提供了临床上可应用的定量方法。年龄，药物治疗，放疗，IDH突变状态，MGMT启动子甲基化状态和GDRG特征包含在最终的OS预测模型中，预后列线图的C指数为0.896。基于时间的ROC分析显示，对于0.5年，1、2年和3年OS率具有良好的预测能力，AUC值分别为0.734、0.771、0.864和0.919。校正图显示了预测的0.5年，1年和3年OS率与TCGA队列中的实际观察值之间的极好的一致性。然后在CGGA队列中验证了预后模型，其C指数为0.729。0.5年，1年，2年和3年OS率随时间变化的AUC与预后诺模图分别为0.725、0.696、0.694和0.701。以上发现表明，预后诺模图方法可用于OS预测，在不同人群的GBM患者中具有较高的可靠性。
图8. 预后列线图可预测GBM患者的0.5年，1年和3年OS
四、小结
文章中作者首先使用单细胞RNA测序（scRNA-seq）数据通过轨迹分析确定了处于不同分化状态的两个GBM细胞亚群，探索了GBM细胞分化相关基因（GDRGs）的生物学功能，随后使用4个可预测OS的关键可预测GDRG（FN1，APOE，RPL7A和GSTM2）和其他临床病理变量为GBM患者构建了分子分型以及预后预测模型进行验证。最后确定了不同的肿瘤内GBM细胞分化状态，根据GBM细胞的分化特征对患者进行分类可以预测肿瘤的免疫治疗反应和患者生存率，对免疫治疗反应具有重要作用。