细胞焦亡+单细胞RNA测序联合分析，搭配实验验证，轻松发到5分+。

今天给同学们分享一篇单细胞+细胞焦亡+实验的生信文章“Exploring the role of pyroptosis in shaping the tumor microenvironment of colorectal cancer by bulk and single-cell RNA sequencing”，这篇文章于2023年5月18日发表在Cancer Cell Int期刊上，影响因子为5.8。

新兴研究显示，细胞焦亡在肿瘤的发展和治疗中起着不可忽视的作用。然而，细胞焦亡在结直肠癌（CRC）中的机制仍然不清楚。因此，本研究调查了细胞焦亡在CRC中的作用。

1. CRC和正常样本之间的差异

在这项研究中，在减少GEO和TCGA队列之间的批次效应后，保留了52个与细胞焦亡相关的基因，以供后续分析。通过比较这些基因在CRC和正常样本之间的表达，作者在GEO队列中确定了35个FDR < 0.05的基因。这些基因被称为DEGs，并以热图的形式显示。如图1A所示，CRC样本中有13个上调基因和22个下调基因。GO富集分析显示DEGs在一些重要的生物过程中发挥作用，如炎症反应调节、细胞焦亡、细胞因子产生、对细菌源性分子的反应、对脂多糖和白细胞介素产生的反应。除此之外，这些DEGs还具有复杂的分子功能，如内肽酶活性调节和细胞因子过程（图1B）。KEGG富集分析的结果显示，差异表达基因(DEGs)参与了许多细胞和生物功能，例如多种类型的感染、NOD样受体信号通路、结直肠癌和炎症性肠病(图1C)。根据富集分析的上述结果，作者推测DEGs与CRC的预后和免疫有关，表明炎症性细胞死亡调节基因的表达失衡可能通过上述途径和机制在CRC的发生和进展中起着关键作用。

图1 构建与细胞焦亡相关的风险模型

2.&nbsp;细胞焦亡相关风险模型的构建

根据纳入和排除标准，分别从GEO和TCGA队列中纳入了556和515个CRC样本，用于后续分析。基于上述差异表达基因（DEGs），进行了单变量Cox回归分析，筛选出与CRC样本预后相关的8个DEGs。NLRP1在CRC样本的生存中起到了风险因素的作用（HR = 1.477（1.058–2.062）；p = 0.022）。然而，CHMP6、GSDMD、GZMB、CYCS、CASP3、CASP1、CASP6的HR < 1（p < 0.05），代表了CRC预后的保护因素（图1D）。为了全面评估每个CRC样本，通过LASSO Cox回归分析，从这8个预后DEGs中建立了一个与细胞焦亡相关的风险模型（图1E，F）。根据该模型计算得出的每个基因的相应系数。如图1G所示，399个CRC样本中有36个（9.02%）发生了与细胞焦亡相关的基因突变。突变频率最高的基因是NLRP1，而CYCS在CRC样本中没有发生突变。就突变类型而言，错义突变是最常见的突变。最后，主成分分析显示高风险和低风险评分组之间存在明显的分离。这些结果表明，与炎症性细胞死亡相关的风险模型能够区分具有不同PRS的CRC样本。而且，在TCGA队列中进一步证实了这个结论。

3.&nbsp;不同PRS的CRC样本之间的生存分析

在培训组中，根据截断值（-7.052），将所有CRC样本分为高风险和低风险评分组。生存分析显示，高风险评分组的CRC样本的总生存期和无进展生存期较差（图2A，B）。随后，绘制了时间相关的ROC曲线，以测试PRS在预测总生存期和无进展生存期方面的准确性。如图2C所示，1年、3年和5年的AUC分别为0.644（p = 0.003）、0.676（p = 1.5e-08）和0.673（p = 2.4e-09）。这意味着PRS在预测总生存期方面具有价值。同样，图2D说明PRS可以作为可靠的指标来预测无进展生存期。预期地，PRS在验证组中仍然能够准确预测CRC样本的总生存期和无进展生存期。图2E和F显示了190个低风险CRC样本和180个高风险CRC样本的遗传突变概况。APC、TP53、TTN、KRAS以及MUC16是突变频率最高的五个基因。APC中的无义突变是最常见的突变类型，而其他基因中最常见的突变类型是错义突变。

图2 训练组中高风险和低风险评分组的预后差异

4.&nbsp;PRS与临床特征之间的关联

在这项研究中，作者调查了不同临床亚组中PRS的分布情况。在GEO队列中，作者发现PRS在年龄和性别上的分布没有显著差异（图3A，B）。然而，结直肠癌样本的PRS与肿瘤进展有关。晚期肿瘤样本的PRS明显高于早期肿瘤样本（图3C）。作者还发现，AJCC-TNM分期较高的结直肠癌样本具有较高的PRS（图3D-F）。TCGA队列中也得到了相同的结果。这意味着具有高PRS的结直肠癌样本往往处于晚期阶段，这也在一定程度上解释了为什么高风险评分组的结直肠癌样本预后较差。CIN（+）亚组中PRS也显著增加（p = 0.049；图3G）。与熟练的MMR（pMMR）亚组相比，缺陷的MMR（dMMR）亚组的结直肠癌样本PRS较低，这表明具有低PRS的结直肠癌可能对免疫治疗敏感（p = 0.0013；图3H）。KRAS突变的结直肠癌样本PRS较高（p = 0.0089；图3I）。然而，作者没有观察到CIMP、BRAF/TP53突变亚组之间PRS的显著差异。除此之外，作者基于GEO队列中的RNA表达数据预测了CRC样本的CMS。160个样本被分为CMS2，其特点是微卫星稳定性（MSS）。89个样本被分为CMS3，其特点是分化和脂肪酸代谢。160个样本被分为CMS4，其特点是上皮间质转化（EMT）。图3J显示了4个CMS亚组中PRS的分布情况。在GEO队列中，单变量Cox回归分析的结果显示年龄、病理分期、AJCC-TNM分期和PRS是CRC的风险因素（图4A）。随后的多变量Cox回归分析表明，只有年龄、AJCC-TNM分期和PRS可以作为独立的预后指标（图4B）。此外，如附加文件1：图S6所示，PRS也被证实是TCGA队列中的独立预后指标。基于这些独立的预后指标，作者开发了一个图表来评估结直肠癌样本的1年、3年和5年生存率（图4C）。时间依赖的校准曲线显示，该图表在预测结直肠癌预后方面具有很高的准确性（图4D-F）。此外，1年、3年和5年的AUC表明，该图表比单个预后指标具有更好的预测能力（图4G-I）。

图3 GEO队列中PRS与临床特征的关联

图4 CRC的生存预测的等高线图

5.&nbsp;CRC的生存预测的等高线图

在GSEA之后，作者发现DNA修复、E2F靶点、G2M检查点、α干扰素反应以及γ干扰素反应在低风险评分组中富集（图5A）。凝血、上皮间质转化、缺氧、肌肉生成以及对紫外线反应下调在高风险评分组中富集（图5B）。这些结果表明高风险评分组与低风险评分组之间可能存在免疫特征上的差异。如图5C所示，作者发现几乎所有与免疫相关的基因在低风险评分组的CRC样本中高表达。虽然高风险评分组和低风险评分组之间的免疫评分、估计评分和肿瘤纯度没有显著差异，但作者发现高PRS的CRC样本具有更高的基质评分（图5D-G）。基质评分代表肿瘤样本中基质细胞的丰度。因此，高基质评分有利于肿瘤细胞的生长和转移，也会影响抗肿瘤免疫效应。同时，作者还探索了PRS与基质相关活动之间的相关性，包括血管生成和泛纤维母细胞TGF-β反应特征（Pan-F-TBRS）。如图5H和I所示，PRS与血管生成呈正相关（r = 0.17，p = 4.29e−05），与Pan-F-TBRS呈正相关（r = 0.35，p = 1.8e−17）。在免疫细胞浸润评估中，作者发现高风险评分组与低风险评分组之间的免疫细胞浸润存在显著差异（图5J）。在低风险评分组的CRC样本中，活化的CD8+ T细胞、活化的CD4+ T细胞、活化的B细胞和记忆B细胞浸润更多。与此同时，高风险评分组中天然杀伤细胞、天然杀伤T细胞、浆细胞样树突状细胞、未成熟树突状细胞、巨噬细胞和肥大细胞的浸润更高。这些结果表明，具有不同PRS的CRC样本可能对免疫治疗的敏感性不同。

图5 高风险评分组与低风险评分组之间的生物特征差异

6.&nbsp;评估和预测CRC的治疗策略

根据高风险评分组和低风险评分组之间的IC50值比较，高PRS的CRC样本的5-氟尿嘧啶IC50值较高（图6A）。这个结果表明基于5-氟尿嘧啶的化疗方案在低风险评分组中效果更好。与高PRS的CRC样本相比，低PRS的样本对抗PD1免疫疗法更敏感（Bonferroni校正p < 0.01）（图6B）。IMvigor 210是关于atezolizumab在局部晚期或转移性尿路上皮膀胱癌患者中的II期试验，atezolizumab是一种改良的抗PDL1抗体。在IMvigor 210队列中，作者发现“炎症”免疫表型患者的PRS显著低于“沙漠”和“排斥”免疫表型患者（图6C）。此外，对atezolizumab有反应的患者（n = 68）的PRS显著低于无反应的患者（n = 230）（图6D）。此外，本研究通过CTRP和PRISM预测了对具有不同PRS的CRC样本可能有效的某些化合物/药物。根据化合物/药物与CRC样本的PRS之间的相关性，筛选出与高PRS或低PRS的CRC治疗潜在药物的前6个具有正或负相关系数的化合物/药物。在CTRP中，KU 0060648、BRD-K50799972、C6-酸酐、VU0155056、BRD-K13999467、BRD-K52037352是高PRS的CRC的前6个潜在药物。而硝苯地尔和异氟氧氮等药物可能用于低PRS的CRC治疗（图6D）。在PRISM中，一些药物如瑞诺酮和异氟氧氮被预测为治疗高PRS的CRC的潜在药物，而替莫普利、乙炔雌酚等可能有助于低PRS的CRC治疗（图6E）。

图6 CRC治疗策略的预测

&nbsp;

7.&nbsp;单细胞水平下与焦亡相关风险模型的探索

作者从GSE132257数据集中收集了5个CRC样本的scRNA-seq数据。在数据过滤和标准化后，选择了方差最大的3000个基因进行后续的细胞分类。这3000个基因的表达水平通过PCA降维到PC1-20。然后，对PC1-20进行t分布随机邻域嵌入（t-SNE）以将所有细胞分类为21个亚型。然后，作者对每个亚型进行了细胞注释。T细胞、B细胞、自然杀伤细胞（NK细胞）、上皮细胞、树突状细胞（DC）、巨噬细胞、单核细胞、成纤维细胞、共同髓系祖细胞（CMP）以及内皮细胞是这些亚型中的主要细胞类型（图7A）。不同细胞亚群中凋亡相关基因的表达和百分比显示在图7B中。此外，作者还应用AUCell评分算法将所有细胞分为低AUC组和高AUC组。作者在所有细胞的AUC值中发现了两个峰值，当AUC值阈值设定为0.039时，4614个细胞显示出相对较高的AUC值。为了进一步研究低和高AUC组中免疫细胞亚群的比例，作者提供了更详细的免疫细胞亚群注释，显示在图7C中。低和高AUC组之间免疫细胞亚群的比例差异显示在图7D中。与高AUC组相比，低AUC组中B细胞、T细胞和DC的比例更高，这可能至少部分解释了低PRS样本对抗PD1免疫疗法更敏感的原因。

图7 通过单细胞RNA测序数据和RT-qPCR探索与细胞焦亡相关的基因

8.&nbsp;正常细胞系和结直肠癌细胞系中与细胞焦亡病相关基因表达水平的比较

为了进一步验证与细胞焦亡相关的风险模型的临床实用性，作者进行了RT-qPCR实验，以量化这些风险模型中8个基因的表达水平，包括正常细胞系（NCM-460）和结直肠癌细胞系（HT-29，HCT 116）。图7E-L显示了每个基因的相对表达水平，用均值和标准偏差（SD）表示。误差线的顶部表示均值+SD，底部表示均值-SD。浅蓝色代表NCM-460，珊瑚色代表HT-29，黄色代表HCT 116。ns表示无显著性差异；*表示p < 0.05；**表示p < 0.01；***表示p < 0.001；****表示p < 0.0001。与NCM-460相比，CASP1和NLRP1在HT-29和HCT 116中下调，而CASP3、CASP6、CHMP6、CYCS、GSDMD和GZMB在HT-29和HCT 116中上调。

9.&nbsp;共识聚类

在这项研究中，作者在高风险评分组和低风险评分组之间确定了6719个FDR < 0.01的差异表达基因（DEGs）。基于这些DEGs，GEO队列中的结直肠癌样本被分成了4个亚型，其中亚型1有202个样本，亚型2有102个样本，亚型3有153个样本，亚型4有99个样本（图8A）。图8B显示了当亚型的数量为2-9时的累积分布函数（CDF）曲线。根据图8C，作者可以发现当k = 4时，CDF面积变化的拐点最显著。t-SNE的结果显示，四个亚型之间有明显的差异，表明不同亚型中的样本可能具有不同的预后和其他生物学特征（图8D）。生存分析显示，亚型2具有明显的生存优势，而亚型4的预后较差（图8E）。此外，这些亚型之间PRS的比较表明，亚型4的PRS显著高于其他亚型（图8F）。如图所示。8G，99个样本中有83个样本（84%）属于cluster4，而这些样本中的高风险评分组占了很大一部分。相比之下，cluster1（67个样本；33%）和cluster2（39个样本；38%）中只有少部分来自高风险评分组。这些结果解释了为什么cluster4的样本有最差的预后。此外，作者还探索了这四个cluster与CMS亚组之间的关系。如图8H所示，作者推测cluster1对应CMS2，cluster2对应CMS1，cluster3对应CMS3，cluster4对应CMS4。通过alluvial图（图8I），作者可视化了个体样本中PRS、CMS和clusters的变化。

图8 共识聚类和CMS亚群

为了更深入地了解不同聚类中样本的特征，作者进行了GSVA评估生物功能和免疫相关活动的富集情况。GSVA的结果显示，聚类2和聚类4与免疫相关活动显著相关。然而，聚类4中也富集了基质激活，如上皮间质转化和血管生成。聚类3中观察到富集了生物代谢，如异物代谢、胆酸代谢以及脂肪酸代谢。聚类1中激活了一些经典的肿瘤通路（图9A）。随后的分析揭示了聚类2中CD8 T效应细胞和免疫检查点的显著激活；聚类4中基质活动，包括血管生成、全面免疫炎症反应评分和上皮间质转化1/2/3，得到了显著增强（图9B）。这些结果进一步确认了聚类与CMS亚组之间的对应关系。

图9 4个亚型之间生物特征的探索

10.&nbsp;高风险和低风险评分组之间的关键基因

从差异表达基因（DEGs）和免疫基因中确定了338个重叠基因，用于构建PPI网络。基因的颜色代表基因的PPI得分（表示相互作用程度），越接近红色得分越高，越接近紫色得分越低。高风险评分组中的上调基因标记为红色，高风险评分组中的下调基因标记为蓝色。根据度数，STAT1、FGF2、HRAS、IGF1、CXCL10、CD8A、EGF、IL17A、CCL5以及IFNG被认为是前10个中心基因。此外，图9C显示了这些中心基因在高风险评分组、低风险评分组和临床亚组中的表达水平。经进一步分析，作者发现CXCL10和FGF2的表达水平在CRC样本的生存中起着重要作用。与高表达CXCL10的样本相比，低表达CXCL10的样本预后较差（图9D）。相反，低表达FGF2的CRC样本比高表达FGF2的样本预后更好（图9E）。作者还发现，CXCL10和FGF2的表达水平与CRC中的B细胞、CD8+ T细胞、CD4+ T细胞和树突状细胞等免疫细胞的浸润呈正相关。

总结

通过细胞焦亡相关风险模型计算的PRS与CRC样本的病理分期、AJCC-TNM分期和预后相关。此外，具有不同PRS的CRC样本在免疫和基质特征上也存在差异，从而对化疗和免疫治疗的敏感性也不同。这些发现表明PRS在CRC中起着关键作用，为肿瘤预后和治疗预测提供了一种新颖而高效的视角。