这是你要的5+单细胞挖掘思路吗？
导语
今天和大家分享的是今年7月份发表在Cancer Immunology, Immunotherapy 杂志（IF=5.442)上的一篇文章，“Integrated analysis of singlecell RNAseq and bulk RNAseq unravels tumour heterogeneity plus M2like tumourassociated macrophage infltration and aggressiveness in TNBC”。文章中作者用scRNA-seq和bulk RNA-seq数据鉴定肿瘤间异质性，应用随机森林和SVM找到与M2型TAM最相关的基因，并构建基于神经网络的深度学习框架以预测患者对免疫治疗的反应。
Integrated analysis of singlecell RNAseq and bulk RNAseq unravels tumour heterogeneity plus M2like tumourassociated macrophage infltration and aggressiveness in TNBC
单细胞RNA-seq和bulk RNA-seq联合分析揭示TNBC中肿瘤异质性以及M2型肿瘤相关的巨噬细胞侵袭特征
一、研究背景
三阴性乳腺癌（TNBC）具有广泛的肿瘤内异质性。尽管近年来乳腺癌治疗的有所进展，但很大一部分的TNBC患者进展到转移阶段且临床表现不佳，通过消除大多数肿瘤组织细胞可能无法有效抑制TNBC。单细胞和组织转录组分析的组合可以表征肿瘤异质性和免疫状况与临床结果之间的联系，并确定免疫状况与临床结果的关联。
二、结果解读
1 TNBC的肿瘤细胞异质性
作者首先进行单细胞RNA-seq分析，TAMs相关基因标记被用来预测生存和免疫治疗反应。UMAP方法用于非线性降维，结果显示，每个患者的肿瘤细胞高度聚集。
表1 TNBC患者(单细胞RNA-seq)信息
TNBC的肿瘤复发与几个侵袭性特征密切相关，如细胞干性、血管生成和EMT，因此作者在单细胞水平上分析了这三个特征。在TNBC细胞中，细胞干性、血管生成和EMT互相呈正相关。作者还分析了三者中分数均较高的细胞（散点图中的红色标记），表明其在肿瘤进展，转移和复发中潜在的协同作用。
图1 TNBC在单细胞水平上的异质性
2 TNBC的免疫细胞异质性
在UMAP降维下免疫细胞被分成三组，这三个集群与免疫细胞类型(T细胞，B细胞，和巨噬细胞)匹配良好。热图显示了识别这三个集群最关键的基因。
图2肿瘤免疫细胞亚群聚类
作者下一步分别对免疫细胞进行分析。TNBC中与T细胞激活和功能状态相关的特征基因集从单细胞转录组数据中提取，如T细胞调节、共刺激、天然、细胞毒性和衰竭。结果显示，BC09患者比其他患者具有更多的细胞毒性T细胞。相比之下，PT058和PT084患者的细胞毒性t细胞活性平均水平较低。具有高细胞毒活性的T细胞表达高水平的CST7、GZMA、GZMB、IFNG、NKG7、PRF1和TNFSF10。此外，TNBC组织中的T细胞即使在单个患者中也表现出混合的幼稚T状态。幼稚t细胞状态可能与t细胞在TNBC中的细胞毒性活性有关。
图3 TNBC单细胞中T细胞信号特征
作者也分析了肿瘤相关B细胞的功能状态，如抗凋亡增殖、相关基因表达特征。B免疫细胞在TNBC患者和同一TNBC患者中表现出异质性分布。总的来说，通过单个细胞的转录组数据， T细胞和B细胞在单个细胞水平上表现出明显差异，无论是在某个患者还是在不同的患者之间。
图4 TNBC单细胞中B细胞水平特征
总的来说，在TNBC中，M2样TAMs是巨噬细胞的主要类型，如下图。作者在巨噬细胞中发现了高M2活性。在多种癌症类型中，M2样TAMs通过诱导血管生成、EMT和干性来促进肿瘤发生。巨噬细胞中显著表达的基因也参与炎症和血管生成功能。M2型基因如TGFBI、MS4A6A和CD163在TNBC的大多数巨噬细胞中表达上调。此外，已知的促进血管生成、EMT和肿瘤进展的基因，如IL8和PLAUR，在TNBC巨噬细胞中过表达。
图5 单细胞水平上整合TNBC免疫细胞浸润分析
3 M2型TAMs与TNBC患者生存期相关
作者收集TCGA TNBC患者的Bulk-RNA-seq数据，利用CIBERSORT算法计算M2样TAMs的丰度。TCGA TNBC患者信息汇总见表2。
表2 TNBC患者(bulk RNA-seq)信息
Cox比例风险回归模型显示，TNBC患者中大量的TAMs与较差的无复发生存(RFS)相关。紧接着作者做了WGCNA分析，共生成7个非灰色模块。在这些模块中，与M2样TAM分数相关性最高的蓝色模块被认为是“TAM模块”。基因显著性代表了单个基因与TAM密度的关联，量化了模块特征值与基因表达谱之间的相关性。散点图表明该模块与M2型TAMs高度相关。
图6 Bulk-RNA seq显示TNBC患者高TAMs与不良预后相关
作者使用PCA方法计算每个患者的保护评分Anti-TAMpca。Anti-TAMpca保护评分根据各基因表达水平和PCA评分计算。Kaplan-Meier图显示，在TCGA TNBC队列中，保护评分高的患者比保护评分低的患者预后更好。作者基于TAM相关基因模型的保护评分应用于GEO外部验证队列，Cox回归分析显示TAM相关模型在GEO队列中的预测价值。通过ssGSEA来确定TAM相关模型的潜在机制，发现胆汁酸代谢、雌激素反应和P53通路与TAM相关保护评分呈高度负相关。
4 通过TAM相关基因进行分型
研究进一步根据共识聚类，将TNBC根据TAM相关基因特征的表达模式划分为不同的亚组。累积分布函数(CDF)图显示了每个k的共识矩阵的CDF，量化了共识矩阵在0到1范围内的分布情况。选择k = 2作为最优参数，将训练队列划分为不同的亚组。热图显示了两个集群之间差异表达的基因。Kaplan-Meier图显示，组1患者的预后优于组2患者。小提琴曲线图显示，与第1组相比，第2组TAM评分显著升高。在TBNC的第1组和第2组中进行GSEA分析，结果表明，上调的通路包括早期雌激素反应相关通路、晚期雌激素反应相关通路、M2巨噬细胞以及p53通路。下调的通路包括E2F靶点、G2M检查点和MYC靶点相关的通路。
图7 根据TAM相关基因特征进行分子亚型鉴定
5 随机森林算法和SVM算法揭示与M2型TAM相关的关键基因
作者通过随机森林算法和支持向量机算法均采用了特征重要性排序方法，提取与M2型TAM评分相关的前50个最重要基因并取交集。显示随机森林算法和SVM算法共有的19个基因，其中VIPR1、ABCD4、SLC40A1和SCN1A是不良预后因素，其余15个基因与更好的RFS相关。
6 TAM相关基因模型预测免疫治疗反应
作者建立了一个基于神经网络框架来预测TAM相关基因模型对免疫治疗的反应。神经网络框架如图8所示。TCGA的TNBC队列作为训练集，GEO队列作为测试集来评估准确性。在训练过程中，测试集的损失值随着迭代的增加而减小，说明训练过程的效率很高。当迭代数达到1000进行训练时，使用混淆矩阵对测试集中的所有样本进行正确识别， ROC曲线下面积(AUC)达到100%。
图8 基于神经网络的深度学习框架构建TAM相关基因模型
三、小结：
作者通过单细胞RNA-seq和Bulk RNA-seq数据，分析了肿瘤细胞和肿瘤免疫浸润细胞之间的异质性，在单细胞水平上发现了血管生成、干性和EMT之间高度相关，并鉴定出了同时高表达的细胞亚群。在TNBC患者中，M2型TAMs与不良预后相关。随后，作者通过随机森林算法和SVM算法发现了几个与M2型TAMs评分高度相关的基因，建立了一个神经网络来预测M2型TAMs相关预后模型，思路清晰。