今年9月底发表在Aging-US（IF=4.831）上的一篇纯生信分析文章（“Glioblastoma cell differentiation trajectory predicts theimmunotherapy response and overall survival of patients”），该文章基于GEO的单细胞测序数据发现了具有不同分化特征的胶质母细胞瘤（GBM）细胞，继而进行差异表达分析找到分化相关的基因（GDRG），最后根据这些基因分别构建了分子分型与预后预测模型。
研究背景介绍
GBM是最常见同时也是死亡率最高的颅内原发肿瘤，5年生存率仅5-6.8%，综合治疗手段，诸如手术、放疗、化疗、免疫治疗等都未能有效延长患者生存时间。肿瘤微环境中的多种因素会影响肿瘤细胞与肿瘤干细胞（CSC）的分化，从而导致细胞分化状态的异质性，而单细胞转录组学分析为我们提供了一个强大的方法，用于研究整个肿瘤样品的基因组本质，继而提供了定义和明确细胞状态及其分化状态的机会。
图1
结果
1.基于单细胞测序数据确定细胞亚群
从GSE84456获取了来源于4个GBM标本的3589个细胞测序结果，仅有2343个来源于肿瘤核心部位的细胞被纳入进一步研究。Seurat R包被应用于单细胞数据分析。经过质控及标准化数据，194个低质量细胞被排除（图2A-B）。总共纳入19752个基因，方差分析显示了1500个高度变异的基因（图2C）。主成分分析PCA未见GBM细胞显著分离趋势，我们最终选了20个主成分行进一步分析（图2D-E）。之后，tSNE分析将GBM细胞分为13个亚群，差异表达分析从所有13个亚群中鉴定出总共8025个标记基因（图2F-G）。根据标记基因的表达量，我们使用singleR和CellMarker对各个亚群进行注释（图3A）。
图2
2.肿瘤分化轨迹分析并确定GBM分化相关基因（GDRG）
使用Monocle 2算法进行时序及轨迹分析，发现肿瘤干细胞主要位于树根部（root），而2个树枝（branch）分布有不同的GBM细胞，branch I包含了434个GBM细胞，而branch II包含了444个细胞（图3B）。差异表达分析得到265个I型GDRG和193个II型GDRG。基因集富集分析（GSEA）发现I型GDRG与免疫反应调控相关通路显著负相关，而II型GDRG与代谢相关通路显著正相关（图3C-D）。
图3
3.TCGA中GDRG的分析
为探索能否在TCGA转录组测序结果中找到具有不同分化特征的GBM亚群，相关性分析表明我们可以通过传统转录组测序数据中GDRG表达来鉴定I型和II型GBM细胞（图4A-C），因为这些高度相关的基因可能表明共同的细胞起源。为了确定来自不同GBM细胞亚群的基因是否在功能上相关，我们利用metagene来代表相应基因的整体表达模式，即I型和II型的metagene分别由I型和II型GDRG组成。图4D-F显示I型和II型的metagene在单细胞及传统转录组数据中高度相关，表明I型和II型GDRG在功能上是高度相互关联的。此外，90.8%的GDRG出现了体细胞突变。
图4
5.免疫检查点分子（ICM）及免疫治疗反应性分析
PD1、PDL1、PDL2在I型GBM细胞亚群（单细胞测序数据）和MC1患者（传统转录组测序）中相对高表达，而CTLA4、CD80、CD86在II型GBM细胞亚群和MC2患者中相对高表达（图6A-C）。根据TIDE算法，MC2患者对免疫治疗的反应性优于MC1。此外，SubMap分析发现MC1对PD1抑制剂反应性更好，而MC2对CTLA4抑制剂反应性更好（图6D-E）。
图6
6.构建基于GDRG的预后风险评分模型
套路化的单因素COX回归，再接着LASSO、多因素COX回归，继而构建一个四基因风险评分模型。分析评分计算公式如下：
根据风险评分的中位数，将所有患者被分为2组，即高危组和低危组，K-M生存曲线提示高危组患者OS显著差于低危组。ROC曲线及C指数均提示较好的预测价值。同样的结论在CGGA验证集中得到验证。
图7
7.构建预后相关的列线图Nomogram
套路化的单因素+多因素COX回归，继而构建Nomogram，C指数、ROC曲线、校正曲线均提示较好的预后预测价值。同样的结论在CGGA验证集中得到验证。
图8
总结
该文章首先挖掘GEO数据库里的单细胞数据，通过分化轨迹分析找到了与分化相关的基因，继而通过TCGA数据库以这些基因为基础分别构建了分子分型以及预后预测模型。该文章思路新颖，可谓生信分析中的一股清流，令编辑和审稿人耳目一新（仅用1个月即接收），大家不妨按照相似的套路模仿一下，想必不难发表高分文章。
