单细胞+氧化应激+转录组预后模型

今天给同学们分享一篇单细胞+氧化应激+预后模型的生信文章“Oxidative Stress Response Biomarkers of Ovarian Cancer Based on Single-Cell and Bulk RNA Sequencing”，这篇文章于2023年1月27日发表在Oxid Med Cell Longev期刊上，影响因子为7.31。

卵巢癌（OV）的发生和发展在很大程度上受到氧化应激（OS）副产物水平升高和缺乏抗氧化应激修复系统的影响。因此，有必要探索与卵巢癌氧化应激相关的标记物，以帮助预测卵巢癌患者的预后和免疫治疗反应。

一、细胞亚群的鉴定

1. 肿瘤单细胞图谱

从 GEO 数据库检索到的单细胞测序数据集 GSE146026 由六名 OV 患者的肿瘤组织组成。在对原始数据进行整合和过滤后，剩余的细胞和基因被用于亚组鉴定和注释。图 1（a）显示了 OV 患者 TME 中细胞亚群的分布和差异。根据单细胞测序数据集的聚类结果和利用SingleR获得的注释信息，利用UMAP降维显示单细胞的表达模式。结果显示，细胞可分为 0 至 16 的 17 个亚群，并被注释为六种主要细胞类型（图 1(b) 和 1(c)）。根据每个细胞亚群的前两个基因绘制了小提琴图。小提琴图显示了每个细胞亚群特定标志物的表达和分布（图 1(d)），从而表明肿瘤细胞具有异质性。根据每个细胞亚群中特异性强的标记基因与操作系统反应相关基因的交叉点，共筛选出 28 个基因。气泡图显示了这些基因在各细胞亚群中的表达情况（图 1(e)）。结果显示，与操作系统反应相关的基因对第 12 组的特异性最强。第 12 组细胞被注释为祖细胞。

图1 从单细胞 RNA 测序数据库中识别细胞亚群和标记基因的表达

2. OS 亚群的鉴定

将与 OS和细胞亚群相关的特定基因交集，共得到 56 个基因。这 56 个基因被定义为 ROS 标记。根据这 56 个 ROS 标记的表达确定的活跃细胞亚群被用来研究单细胞水平上 OS基因的表达模式。使用最佳阈值确定细胞活性，结果显示 OS 亚群中有 1751 个活性细胞（图 2(a)）。图 2(b) 和 2(c) 显示了活跃细胞的 UMAP 图和累积分布直方图。结果显示，祖细胞和树突状细胞是活跃细胞。图 2(d) 显示了活跃细胞亚群中前十个标记基因的表达情况。

图2 识别活跃细胞亚群

3. 确定 OS 反应细胞亚群的功能

利用 GSEA，根据特异性强的标记基因的表达确定活跃 OS 细胞群的功能。对标记基因组进行了 GO 和 KEGG 通路富集分析，并选取前 10 个显著富集的通路绘制气泡图。结果显示，标记基因组明显富集了与细胞外组织相关的功能以及与粘连斑块和癌症中蛋白多糖相关的通路。使用基于 log2FC 值的 HALLMARK 通路进行了 GSEA。共有八条通路被显著富集。在干扰素反应因子和肿瘤中的免疫调节之间观察到了密切的相关性。

4. OV 样本的批量 RNA 序列分析及基因表达模式

从 TCGA 数据库中获得的 OV 数据不包含正常样本，因此作者进行了批量校正以消除批量效应。将基因型-组织表达中的正常样本和TCGA-OV样本进行整合，得到352个OV样本和88个正常样本，用于差异表达分析。共获得 2928 个 DEGs。OV患者的临床信息见表1。将 DEGs 与活跃细胞亚组的标记基因进行交叉，共得到 151 个差异表达活跃的标记基因。使用 HALLMARK 通路对 DEGs 进行了 GSEA，结果显示 DEG 显著富集了 HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION 和 HALLMARK_INTERFERON_GAMMA_RESPONSE。

二、预后风险模型的构建与验证

1. 识别 OV 的预后基因特征

采用单变量 Cox 分析来识别活跃细胞群中差异表达的标记基因（P < 0.05）。最终共鉴定出 12 个与 OV 预后相关的基因。随机抽样选取 7/10 个 TCGA-OV 样本（n = 352）作为训练集（n = 246）。对训练集进行 LASSO 回归分析以去除冗余基因，并设定种子 = 1110。结果如图 3(a)-3(c) 所示。通过 Cox 回归分析确定了这九个基因，并以每个基因的中位表达值为临界值，将患者分为高危和低危两组。绘制了TCGA-OV数据集中患者的KM生存曲线，并选择了能显著预测患者生存期的基因作为代表。在这些基因中，观察到八个基因之间的 KM 曲线存在明显差异。与低风险组相比，高风险组患者的预后更差（图 3(d)-3(f)）。

图3 TCGA-OV 数据集的 Cox 和 LASSO 回归分析

2. TCGA-OV 数据集的 Cox 和 LASSO 回归分析

为了确定利用九个基因特征建立的模型的稳健性，以 RS 中位值为阈值将样本分为高风险组和低风险组。构建了 TCGA 训练队列（图 4(a)）、整个 TCGA 队列（图 4(c)）、GSE17260 队列（图 4(e)）和 GSE26712 队列（图 4(g)）中不同组别患者的 KM 生存曲线。结果显示，与所有队列中的低风险组患者相比，高风险组患者的预后明显较差。ROC 曲线用于确定模型预测患者预后的效果（图 4(b)-4(h)）。在TCGA训练队列中，1年、3年和5年生存率的AUC分别为0.647、0.711和0.756（图4(b)）。整个 TCGA 队列 1 年、3 年和 5 年生存率的 ROC 分别为 0.624、0.674 和 0.723（图 4(d)）。在 GSE17260 队列中，1 年的 ROC 为 0.668，3 年的 ROC 为 0.662，5 年的 ROC 为 0.681（图 4(f)）。GSE26712 组群 1 年、3 年和 5 年的 AUC 分别为 0.626、0.722 和 0.679（图 4(h)）。这些结果表明，作者的预后模型在不同队列中的表现良好。

图4 模型在不同队列中的表现

三、OV 风险评分模型与临床特征之间的相关性

为了验证 RS 是否可作为独立的预后因素，作者对患者的临床特征（如年龄、FIGO 分期和分级）进行了单变量和多变量 Cox 回归分析。结果显示，无论采用哪种类型的 Cox 回归分析进行统计检验，RS 都是患者的独立预后风险因素（图 5(a)和 5(b)）。多变量 Cox 回归分析用于构建提名图，结果表明 RS 可显著预测临床结果（图 5（c））。接下来，作者检验了根据 RS 建立的预测模型能否预测不同特征患者的预后。KM生存曲线显示，与不同临床特征的低风险组患者相比，高风险组患者的预后较差（图5(d)-5(g)，年龄＜60岁，P＜0.0001；年龄≥60岁，P＜0.0001；III/IV期，P＜0.0001；3/4级，P＜0.0001）。

图5 RS是临床特征的独立预后因素

四、基因突变的差异

基因突变在癌症的发生和发展中起着重要作用。因此，加强对 OV 基因突变的了解将有助于开发靶向药物和新的肿瘤疗法。因此，作者对高危组和低危组患者中突变频率最高的前 20 个基因进行了分析，并绘制了瀑布图。结果显示，高危组和低危组患者的基因突变频率存在差异（图 6(a)和 6(b)）。此外，还计算了 TMB，并根据上四分位数将患者分为高 TMB 组和低 TMB 组（表 S12）。绘制了高 TMB 组和低 TMB 组患者的 KM 生存曲线。结果显示，高 TMB 组患者的预后良好（图 6(c)）。散点图显示 RS 与 TMB 呈显著负相关（图 6(d)，R = -0.13，P < 0.05）。

图6 模型组之间的单核苷酸变异（SNV）

五、分析模型的免疫细胞浸润特征

TME中的免疫细胞浸润会影响癌症的发生和发展。因此，作者使用CIBERSORT算法探讨了高危组和低危组患者TME的差异。作者使用 CIBERSORT 估算了渗入 OV 的 22 个免疫细胞的比例。分析结果显示，高危组和低危组患者肿瘤中浸润的免疫细胞比例存在明显差异。结果显示，在高风险组/低风险组患者中，M1 巨噬细胞、活化记忆 CD4 T 细胞、T 滤泡辅助细胞和γ δ T 细胞等浸润免疫细胞的比例存在明显差异（图 7(a)）。图 7（b）和 7（c）显示了高风险组/低风险组患者中 22 种免疫细胞的分布情况。免疫检查点是免疫细胞表达的分子，可以调节免疫激活的程度。它们在自身免疫性疾病的发生中起着重要作用。因此，作者分析了五类免疫调节剂、细胞因子[ 26]和九个模型基因表达之间的相关性。结果显示，CXCL10 与免疫检查点之间存在很强的相关性（图 7(d)）。

图7 分析模型的免疫细胞浸润特征

六、不同风险组患者的富集途径差异

根据 OV 患者的基因表达模式计算了 HALLMARK 和 KEGG 通路富集分析得分。高危组和低危组富集通路的差异将有助于破译可能影响患者预后的潜在致癌机制。研究人员使用 "limma "R软件包（调整后的P值小于0.01）计算了高风险组和低风险组之间基于富集得分的差异，并绘制了热图。结果显示，高风险组与低风险组之间有 46 条通路的富集得分存在显著差异（图 8(a)）。在富集的通路中，WNT/β-catenin 信号通路在介导肿瘤发生和氧化应激之间的串联过程中发挥了重要作用。WNT/β-catenin 信号通路的异常激活可刺激 ROS 生成和慢性炎症激活。此外，根据 HALLMARK 和 KEGG 通路富集分析所富集的通路，计算 RS 与富集得分之间的相关性。绘制了热图，结果如图 8（b）和图 8（c）所示。

图8 高风险组和低风险组之间的通路富集差异

七、预测药物敏感性和对免疫疗法的反应

根据TCGA-OV样本的表达数据，预测了GDSC数据库中138种药物的IC 50 值。小提琴图显示了高风险组和低风险组之间 IC 50 值分布的差异。图 9(a)-9(h)显示了有显著差异的 8 种化疗药物的结果。计算了模型基因的表达与化疗药物 IC 50 值之间的相关性。结果显示，CXCL10 与化疗药物的 IC 50 值有很强的相关性（图 9(i)）。此外，为了确定该模型能否预测免疫治疗反应，作者根据 RS 的中位值将 IMvigor210 数据集分为高风险组和低风险组。利用累积分布图比较了高风险组和低风险组免疫应答样本（完全应答（CR）+部分应答（PR））和非免疫应答样本（疾病稳定（SD）+疾病进展（PD））的分布差异（图 9(j)）。与高危组患者相比，低危组患者对免疫疗法的反应更好。此外，SD/PD 组患者的 RS 明显高于 CR 组患者（图 9(k)），这表明低风险组患者对免疫疗法的反应更好。

图9 RS 模型可预测 OV 患者的治疗结果

总结

总之，作者利用单细胞和大量 RNA 序列数据创建了一个风险评分模型。用于构建模型的基因是基于OV中与OS反应和通路相关的基因以及与OS相关的基因的交叉。九基因特征预后模型可用于预测OV患者的预后和免疫反应。