单细胞转录组测序分析结肠癌微环境中免疫细胞亚群与癌症进程的关系

周彦彤，李春晓，王劲松，孙芳洲，许东奎，张雪燕，袁芃，钱海利*

1国家癌症中心/国家肿瘤临床医学研究中心/中国医学科学院北京协和医学院肿瘤医院分子肿瘤学国家重点实验室，北京 100021；2国家癌症中心/国家肿瘤临床医学研究中心/中国医学科学院北京协和医学院肿瘤医院肿瘤结直肠外科，北京 100021；3国家癌症中心/国家肿瘤临床医学研究中心/中国医学科学院北京协和医学院肿瘤医院肿瘤特需医疗部，北京 100021

根据2020年世界卫生组织(WHO)发布的《2020年世界癌症进展报告》，结直肠癌发病率在世界范围内居第3位，病死率居第2位，是严重威胁人类健康的疾病[1]。为了攻克癌症，新的治疗技术层出不穷，其中免疫治疗展现出不俗的治疗效果，于2013年被Science杂志评为该年度最重要的科学突破。为了使免疫治疗这一手段让更多的患者受益，针对肿瘤微环境中的免疫细胞与肿瘤细胞的相互调控及其影响肿瘤细胞行为的研究成为热点[2]。有研究显示，肿瘤微环境中的免疫细胞具有多种协助肿瘤发生发展的生物学行为，包括参与肿瘤局部血管生成、协助肿瘤细胞发生转移及免疫逃逸等[3-6]。但目前在探究免疫细胞与肿瘤细胞相关调控行为的研究中，对免疫细胞表型的划分较为笼统，虽然基于单细胞测序及分析技术的提升，已经能在结直肠癌免疫微环境中识别更细分的免疫亚群，但这些免疫亚群的分布及其在大样本临床队列中的潜在价值尚不清楚。CIBERSORTx算法是一种机器学习方法，可用于估计细胞比例及推断细胞类型特异性基因表达谱[7-8]。本研究以单细胞转录组测序(scRNAseq)结果为参考，利用CIBERSORTx算法构建特征表达谱，测算癌症基因组图谱(TCGA)结肠癌队列中肿瘤的免疫细胞浸润情况，并评价其对肿瘤患者预后的预测价值。

1 资料与方法

1.1资料数据来源 通过基因表达数据库(GEO)获得GSE146771[9]中43 817个单细胞转录组，分析CD45+细胞样本的每千碱基记录本(transcripts per kilobase million，TPM)数据，调整样本特征(nFeatures)数值，进一步筛选获得更为匀质的13 193个细胞。通过10×测序平台，获得这些细胞的转录组数据。结肠癌单细胞转录组测序数据包含来自9例结肠癌患者的外周血、术中切除的正常肠黏膜组织、结肠腺瘤及腺癌组织的单细胞测序基因表达谱；从加州大学圣克鲁兹分校癌症功能基因组数据库(UCSC XENA)获取TCGA结直肠癌队列，该队列包含临床结直肠癌样本513例，收集其临床资料如年龄、性别、TNM分期、随访时间及生存状态等，并将每千个碱基的转录每百万映射读取片段(fragments per kilobase per million，FPKM)值的转录组测序(RNAseq)数据转换为TPM数据，以便与scRNA-seq组数据保持一致。

1.2建立scRNA-seq特征矩阵 使用CIBERSORTx在线分析平台(https://cibersortx.stanford.edu/)，在不进行物理分离的情况下，分析细胞类型特异性基因表达谱。按照指令上传单细胞表达矩阵，在默认参数保持不变的条件下，用CIBERSORTx计算获得38个免疫细胞类型的签名矩阵。

1.3CIBERSORTx计算细胞组分 依据CIBERSORTx的提示，上传TCGA及GEO，获得结直肠癌表达谱数据集，选取之前以scRNA-seq为参考数据得到的签名矩阵，选择S模式进行批量校正。将排列数设置为500，其他参数保持默认值。运行CIBERSORTx后，获得了每个样本中38个肿瘤浸润免疫细胞亚群的相对比例，将P＜0.05的样本纳入进一步研究。

1.4免疫细胞亚群与转移的相关性分析 为了探究免疫细胞亚群与肿瘤转移的关系，将513例TCGA结直肠癌队列分为转移组(216例)与未转移组(297例)，比较两组免疫细胞亚群分数的差异，并分析免疫细胞亚群与淋巴结转移数的相关性。

1.5基于免疫细胞亚群的生存分析 为探究不同免疫细胞亚群对结肠癌预后的影响，使用Kaplan-Meier法绘制患者的生存曲线，用log-rankχ2检验分析患者总生存期与免疫细胞亚群细胞分数的关系。

1.6免疫细胞亚群功能分析 采用功能富集分析数据库(the database for annotation, visualization and integrated discovery，DAVID)以及单细胞转录组分析和注释工具(single-cell transcriptome analysis of pathway activation signatures，scTPA)进行功能聚类。将CIBERSORTx分析获得的CD4+T亚群的评分与其对应的表达谱的基因列表进行相关性分析，获得的相关性列表剔除在先前评估细胞分数中使用的基因列表，在剔除后的列表中取正负相关性最高的前500个基因进行功能聚类。对非免疫反应通路进行比较，观察不同CD4+T亚群对结肠癌细胞功能调控的倾向。在聚类分析结果中，除由于免疫细胞亚群差异导致的免疫功能相关通路的差异，还着重关注免疫细胞亚群差异对结肠癌肿瘤细胞生物功能的影响，关注其对肿瘤细胞周期、细胞生长、细胞增殖、能量代谢、氧化磷酸化、缺氧应激等方面的调控差异。

1.7基于免疫细胞亚群的预后模型的构建 采用LASSO算法筛选免疫细胞亚群，并使用Cox回归分析构建风险预后模型。

1.8统计学处理 采用SPSS 23.0软件进行统计学分析。正态分布的计量资料组间比较采用独立样本t检验。不符合正态分布的计量资料比较采用Wilcoxon秩和检验。P＜0.05为差异有统计学意义。

2 结 果

2.138个免疫细胞亚群的签名矩阵 取GSE146771检测基因数目nFeatures_RNA值在800～1200的匀质细胞的测序数据作为后续分析的参考集，基于CIBERSORTx算法创建包括38个免疫细胞集群共3801个基因的特征矩阵。如图1所示。

2.2不同TNM分期的免疫细胞亚群的比例差异

与TCGA结肠癌未转移患者比较，转移患者中hT04_CD4-TCF7分数上调(P＜0.01)，hT11_CD4-CTLA4、hT12_CD8-LEF1、hB01_PlasmaB-IgG、hM10_Macro-IL1B分数下调(P＜0.05，图2A)，且上述hT04_CD4-TCF7、hT12_CD8-LEF1、hB01_PlasmaB-IgG、hM10_Macro-IL1B免疫细胞亚群与淋巴结转移数存在明显相关关系(P＜0.05，图2B)。

2.3免疫细胞亚群与总生存期的相关性 如图3所示，hT04_CD4-TCF7、hT03_CD4-GNLY、hT09_CD4-CXCL13、hT16_CD8-CD6、hT17_CD8-CD160、hT14_CD8-CX3CR1、hI02_NK-GZMK、hB03_FollicularB-IgD、hB01_PlasmaB-IgG、hB02_GALTB-IgA、hM01_Mast-TPSAB1、hM12_TAMC1QC、hM04_cDC1-BATF3、hM03_cDC2-CD1C、hM05_Mono-CD14与患者的总生存期明显相关(P＜0.05)。

2.410个CD4+T细胞亚群对结肠癌细胞调控功能的异质性 CIBERSORTx分析共获得10个CD4+T亚群的评分。其中CD4-GNLY、CD4-ANXA1、CD4-CXCR6对结肠癌肿瘤的调控作用优先聚集在血管生成相关通路，CD4-IL23R优先聚集在能量代谢相关通路，CD4-GZMK优先聚集在缺氧应激相关通路，CD4-ANXA1、CD4-TCF7、CD4-CTLA4优先聚集在细胞增殖相关通路，CD4-ANXA1、CD4-TCF7、CD4-CTLA4、CD4-CXCL13优先聚集在细胞黏附及迁移相关通路(图4A)。此外，CCR7、CXCR5两个亚群富集结果显示，这两个亚群在促癌生物学过程中的作用较其他亚群不显著。

比较转移组与非转移组间存在明显差异的免疫细胞亚群，CD4+T细胞亚群中T04_CD4-TCF7与转移呈正相关，T11_CD4-CTLA4与转移呈负相关。使用scTPA线上工具比较两个亚群细胞内通路激活的差异，发现两者在细胞黏附及免疫激活或抑制相关通路中的激活状态存在明显差异(P＜0.05，图4B)。生存相关性分析结果显示，T04_CD4-TCF7与患者生存呈正相关关系，T09_CD4-CXCL13与患者生存呈负相关关系(P＜0.05)。与CD4-TCF7相比，CD4-CXCL13在抗肿瘤及TCR相关通路中的激活上明显上调(P＜0.05，图4C)。生存分析结果显示，T17_CD8-CD160对TCR相关通路、Fc受体及免疫刺激细胞因子相关的通路激活明显高于T16_CD8-CD6细胞亚群(P＜0.05，图4D)。

2.5基于免疫细胞亚群反映患者免疫状态的预后风险模型的构建 通过LASSO算法从38个免疫细胞亚群中选择4类免疫细胞亚群CD160-CD8+T细胞、CX3CR1-CD8+T细胞、BATF3-cDC1细胞、IgG-PlasmaB细胞(图5A)。其中CD160-CD8+T细胞(P=0.0064)、IgG-PlasmaB细胞与风险评分呈负相关，CX3CR1-CD8+T细胞(P=0.0098)、BATF3-cDC1细胞(P=0.011)与风险评分呈正相关(P＜0.05)。采用Cox分析建立预后风险评分模型，根据整个TCGA结肠癌队列获得的风险评分，患者可分为高风险组与低风险组(图5B)。Kaplan-Meier曲线表明，在训练集、验证集及整组中，高风险组患者的死亡风险较高(图5C–E)。在免疫细胞亚群功能差异的分析中，CD160-CD8+T细胞比CD8+T细胞表现出更高的Fc受体(Fc receptors，FcR)激活水平及T细胞抗原受体(T cell receptor，TCR)激活水平(P＜0.05)。高风险组患者的T细胞功能失调评分高于低风险组，差异有统计学意义(P＜0.05，图5F)。

3 讨 论

免疫治疗在结肠癌中的应用较少，尚处于较初级的探索阶段。鉴别结肠癌的肿瘤微环境特征，特别是免疫细胞亚群及功能状态的特点，以及可能为临床精准实施免疫治疗提供参考的可靠预测标志物，对于提高结肠癌患者的免疫治疗效果，提出可行的辅助治疗方案具有重要意义。近年来，单细胞测序技术的发展使研究者能在更高分辨水平研究肿瘤微环境中的各种细胞成分，其中包括免疫细胞亚群在基因表达及功能状态上的异质性。一系列已经发表的单细胞测序分析结果显示，不同癌种的免疫细胞亚群谱存在差异[10-12]，癌种间免疫治疗效果的差异可能与其免疫细胞谱系的差异相关，目前相关研究主要通过不同的表面标记分子，结合免疫组化空间定位分析免疫细胞亚群的功能及预后价值[13-14]。免疫细胞调控肿瘤细胞功能及肿瘤发生发展的研究更多着眼于传统的免疫细胞分群，而探讨免疫微环境中免疫细胞异质性的单细胞测序分析多侧重于免疫细胞本身的基因表达特性及功能特性。为进一步揭示结肠癌肿瘤微环境中各种免疫细胞亚群对结肠癌细胞功能调控的差异以及对肿瘤进展的影响，本研究利用CIBERSORTx算法，基于单细胞测序分析定义的高分辨结肠癌肿瘤免疫微环境免疫细胞谱对TCGA的结肠癌患者队列进行研究，以期探究结肠癌微环境中特征免疫细胞谱对结肠癌肿瘤细胞的调控作用及其与肿瘤进程的关系，结果显示，hT04_CD4-TCF7、hT11_CD4-CTLA4、hT12_CD8-LEF1、hB01_PlasmaB-IgG、hM10_Macro-IL1B细胞亚群与结肠癌淋巴结转移相关，hT04_CD4-TCF7、hT03_CD4-GNLY、hT09_CD4-CXCL13、hT16_CD8-CD6、hT17_CD8-CD160、hT14_CD8-CX3CR1、hI02_NK-GZMK、hB03_FollicularBIgD、hB01_PlasmaB-IgG、hB02_GALTB-IgA、hM01_Mast-TPSAB1、hM12_TAM-C1QC、hM04_cDC1-BATF3、hM03_cDC2-CD1C、hM05_Mono-CD14细胞亚群与患者的总生存期相关。

本研究还探索了不同CD4+T细胞亚群对结肠癌肿瘤细胞调控功能的差异，发现对于肿瘤血管生成、肿瘤细胞生长增殖、能量代谢、氧化磷酸化、缺氧应激等通路的调节，不同CD4+T细胞亚群表现出不同的优先调控方向，从而影响肿瘤的发展进程，导致患者在转移和预后结局方面的不同。

在转移与预后分析中，本研究发现，CD4+和CD8+T细胞的不同亚群对转移及预后具有相反的相关关系，通过比较这些突出亚群的细胞内通路激活状态，发现这些细胞亚群间通路激活状态的差异很好地解释了其对肿瘤发生发展方向的影响。先前的分析只是单独地看待这些免疫细胞亚群的功能差异，而经过构建模型，可以看到这种显著影响患者生存的肿瘤微环境差异主要来自于哪些免疫细胞亚群，从而可以从整体观察免疫细胞亚群对患者疾病进程的影响。肿瘤部位基于免疫细胞亚群的免疫细胞构成的异质性及免疫状态，对于患者的生存具有较高的影响，该模型在训练集、验证集及整组中显示出较好的基于免疫状态的死亡风险预测能力，相较单纯基于差异基因簇，基于免疫状态预测患者的死亡风险具有更高的临床参考价值，展现了患者检测节点的免疫状态。在这一基于免疫细胞亚群构建的预后模型中，笔者发现，低风险组与高风险组中有4个免疫细胞亚群的构成存在差异，这可能表明某些免疫细胞亚群之间存在免疫应答过程中的相互作用，降低肿瘤微环境中BATF3-cDC1树突细胞的浸润可降低CX3CR1-CD8+T细胞的浸润，或提升IgG-PlasmaB细胞的浸润，进而提升CD160-CD8+T细胞的浸润，改善CD8+T细胞群的功能活性，从而改善患者生存，或辅助免疫检测点治疗。

综上所述，结肠癌肿瘤微环境中的免疫细胞构成及免疫微环境状态能够通过对结肠癌细胞功能的差异调控影响结肠癌的进展，进而影响患者的预后。本研究基于单细胞测序数据，使用CIBERSORTx算法在TCGA队列中研究免疫细胞谱对结肠癌进展及预后的影响，为进一步明确各免疫细胞亚群对肿瘤微环境的影响提供了方向和线索。但本研究尚有一定的局限，如研究结果基于公共数据和算法分析，未能阐述各免疫细胞亚群之间的相互关系网络，且需要进一步的细胞生物学实验对结果进行验证。