今天给同学们分享一篇单细胞分析+肿瘤预后模型构建+验证实验思路的生信文章“Identification of a novel immune-related gene signature for prognosis and the tumor microenvironment in patients with uveal melanoma combining single-cell and bulk sequencing data”，这篇文章于2023年1月30日发表在Front Immunol 期刊上，影响因子为7.3。
葡萄膜黑色素瘤（UVM）是成人中侵袭性最强的眼内恶性肿瘤，预后较差。越来越多的证据表明，免疫相关基因与肿瘤的发生和预后有关。本研究旨在构建葡萄膜黑色素瘤的免疫相关预后特征，并明确分子和免疫分类。
1.葡萄膜黑色素瘤聚类的构建和验证
在分析中，作者从 TCGA 数据库中获得了 80 个 UVM 样本。然后，采用ssGSEA方法对UVM样本的RNA-seq数据进行量化，以评估29种免疫细胞的浸润水平。计算并得到每个 UVM 样本的 ssGSEA 分数。生成的热图说明了肿瘤微环境（TME）中免疫细胞景观之间不同的关联模式（图 1A）。根据ssGSEA的得分，通过无监督分层聚类算法，作者发现了两个具有不同免疫浸润模式的聚类，包括Immunity_H聚类（n = 10）和Immunity_L聚类（n = 70）（图1B，C）。为了验证上述分组方法的可行性和实用性，根据每个 UVM 样本的表达谱，使用ESTIMATE 算法计算了两个群的肿瘤纯度、基质得分、免疫得分和ESTIMATE得分。与免疫_L群相比，免疫_H群的基质得分、免疫得分和ESTIMATE得分更高，但肿瘤纯度则相反（图1D）。小提琴图还显示，免疫_H 群组的基质得分、免疫得分和ESTIMATE得分均高于免疫_L 群组（p < 0.001，图 1E）。此外，作者还发现，与免疫_L群相比，免疫_H群中大多数HLAs的表达量都变高了（均 p < 0.05，图 1F）。在对Immunity_H和Immunity_L簇间基因的KEGG分析中，富集了多种生物过程，如自然杀伤细胞介导的细胞毒性、T细胞受体信号通路、趋化因子信号通路、B细胞受体信号通路、NOD样受体信号通路、Toll样受体信号通路等（图1G）。更具体地说，在所测试的 22 种免疫细胞类型（包括 APC 细胞、aDC 细胞、T 细胞和 Th1 细胞）中，免疫_H 群组的免疫得分明显高于免疫_L 群组（图 1H）。
图1 葡萄膜黑色素瘤聚类的构建与验证
2.免疫力高、低群组间差异表达基因及免疫相关基因的鉴定
根据|logFC| > 1和adj p < 0.05的阈值，作者探索了TCGA数据库中Immunity_H和Immunity_L集群之间的DEGs。随后，作者得到了3122个DEGs，包括2588个上调基因和534个下调基因（图2A）。热图显示了免疫_H和免疫_L群中DEGs的表达水平（图2B）。此外，还从 ImmPort 数据库中下载了 2498 个免疫相关基因。热图显示了免疫相关基因在 Immunity_H 和 Immunity_L 群组中的表达水平（图 2C）。此外，作者还根据 DEGs 和 ImmPort 数据库中的免疫相关基因进行了双向 Venn 分析，在两个基因集中都发现了 316 个基因。
图2 差异表达的免疫相关基因分析
3.葡萄膜黑色素瘤免疫相关基因预后特征的鉴定
在将临床病理信息整合到基因表达图谱中后，作者获得了80个紫外线瘤样本。为了研究能预测 UVM 预后的免疫相关基因，作者对 316 个交叉基因的作用进行了单变量 Cox 回归分析。根据 p < 0.001 的标准，28 个基因与 OS 相关（图 3A）。此外，为了确定 28 个免疫相关基因与转录因子之间的相关性，作者进行了相关性分析。作者绘制了冲积图来表示 28 个免疫相关基因与转录因子之间的关系（图 3B）。此外，作者还通过 STRING 数据库建立了这些基因的相互作用网络，并通过 Cytoscape 显示出来（图 3C）。为了缩小预后特征的范围，对28个免疫相关基因采用了LASSO回归算法，并通过1000轮交叉验证确定了参数的最佳值。LASSO分析结果表明，当预后模型包括MMP9、S100A13和SEMA3B等三个免疫相关基因时，预后模型的性能最佳（图3D、E）。结合 LASSO 算法和 SVM-RFE 算法所剔除的基因，发现这两种算法同时选择了 3 个基因（图 3F、G），这 3 个基因被认为是潜在的分类和预后特征。如图 3H 所示，MMP9 与 B 细胞记忆相关（p, < 0.05），SEMA3B 与巨噬细胞 M1 和浆细胞相关（分别为 p < 0.01 和 p < 0.05）。利用 PCA 图验证了从 DEGs 中筛选出的三个预后基因在不同免疫集群中的分布情况。这表明训练队列中的预后基因在免疫_H 群组和免疫_L 群组之间呈反方向倾斜（图 3I）。根据这三个基因的表达水平，构建了一个预后特征来计算每个样本的风险得分。风险得分=（MMP9 的表达量*0.602435）+（S100A13 的表达量*1.333449）+（SEMA3B 的表达量*-1.092650）。根据中位风险评分，TCGA 数据库中的 80 例 UVM 患者被分为高风险组和低风险组。该 PCA 图显示，作者从 DEGs 中筛选出的预后基因的分布有两个方向，即低风险组和高危组（图 3J）。此外，在检测的 22 种免疫细胞类型（包括 APC 细胞、B 细胞、MHC 细胞、T 细胞、Th1 细胞和 Th2 细胞）中，高风险组的免疫评分明显高于低风险组（图 3K）。这些研究结果发现了三个与免疫相关的基因，它们是紫外线损伤患者高度敏感和特异的预后指标。
图3 葡萄膜黑色素瘤免疫相关基因预后特征的鉴定
4.葡萄膜黑色素瘤免疫相关基因预后特征的验证
柱状图显示，高危组患者的死亡比例高于低危组（图 4A）。Kaplan-Meier 生存曲线和对数秩检验表明，高风险组患者的 OS 明显低于低风险组（图 4E）。预测 1 年和 3 年 OS 率的 ROC 曲线下面积（AUC）值均高于 0.8（图 4I）。图 4M 和 5Q 显示了 UVM 样本的风险评分和生存状态分布。作者还观察到，低风险组的免疫抑制细胞因子也出现了下调（图 4U）。
图4
图5
为了验证免疫相关基因预后特征在预测 UVM 患者 OS 方面的有效性和稳健性。作者将 GSE44295GSE84976 数据库中的 57 个 UVM 样本和 28 个 UVM 样本作为验证集。根据上述免疫相关预后模型，作者按照中位风险评分将验证集中的患者分为高风险组和低风险组。如图 4B 和 C 所示，高风险组患者的死亡率高于低风险组。Kaplan-Meier 生存曲线显示，高危组比低危组的生存结果更差（p = 0.009 和 p < 0.001，图 4F、G）。此外，ROC 曲线显示，在 GSE44295 中，风险评分预测 1 年、3 年和 5 年生存率的 AUC 分别为 0.833、0.777 和 0.631（图 4J）。在 GSE84976 队列中，1 年生存率的 AUC 为 0.844，3 年生存率的 AUC 为 0.850（图 4K）。此外，预后特征的风险评分和生存状态也显示在散点图中（图 4N、O、R、S）。作者还观察到，低风险组的免疫抑制细胞因子也出现了下调（图 4V，W）。作者的结果表明，风险评分是预测 UVM 患者 OS 的良好模型。
5.抗肿瘤治疗临床反应、免疫状态和免疫微环境分析
作者使用 R 软件包 "pRRophetic "探讨了风险评分与对几种免疫治疗药物和一些化疗药物的临床反应之间的联系。以抗肿瘤药物的半数最大抑制浓度（IC50）来衡量，高风险组和低风险组对前 15 种化学疗法或靶向疗法的敏感性存在显著差异（图 5A）。大多数免疫检查点在高危组中更为活跃。此外，作者还发现免疫治疗的几个免疫检查点基因，如 CD70、CD48、CTLA4 和 TNFSF4 基因表达在高危组中的增加，显示它们在各组中的影响不同（图 5B）。这表明，根据风险评分，作者可以为 UVM 患者选择最佳的检查点抑制剂。根据 CIBERSOFT 算法，作者发现风险评分与 CD4 T 细胞、B 细胞、NK 细胞、树突状细胞和巨噬细胞的免疫细胞浸润水平呈正相关（图 5C）。作者还评估了风险评分与免疫疗法的免疫通路之间的关系，如致癌通路、与靶向疗法相关的基因特征和辐射反应基因特征。大约所有与抗癌免疫疗法相关的标记物在风险评分和富集评分之间都有正相关性（图 5D）。
6.TMB、预后因素和免疫相关基因预后特征提名图的评估
为了评估免疫相关基因预后模型是否可以作为一个独立的预测因子，研究人员进行了单变量和多变量 Cox 回归分析。上述结果表明，新的预后模型可以作为一个独立的预后因素来影响紫外线吸收患者的生存结果（图 6A、B）。为了预测 UVM 患者 1 年、3 年和 5 年的 OS 率，作者在 TCGA 数据库的基础上建立了一个提名图。作者选择年龄、性别、分期、TMN 状态和风险评分作为变量（图 6C）。
图6
为了确定免疫相关基因预后特征与肿瘤突变负荷（TMB）之间的潜在相关性，作者分析了高危组和低危组的TMB水平。结果显示，高危组和低危组的 TMB 水平无明显差异（图 6D）。Kaplan-Meier生存分析也表明，高危组患者的OS概率与低危组无明显差异（图6E）。接下来，作者分析了免疫相关基因的预后特征与 TMB 的相关性。为了确定TMB状态的作用，作者对风险组进行了生存分析，分别为低风险组/低TMB、低风险组/高TMB、高风险组/低TMB和高风险组/低TMB。如图 6F 所示，这四组之间存在显著差异（P < 0.001）。总之，这些结果表明风险评分与体细胞突变模式相关。
7.单细胞测序数据分析
在应用质量控制标准后，作者使用 UMAP 方法降低了这些数据集的维度，以说明单细胞测序图谱的分布情况（图 7A）。通过无监督分类，细胞被有效地分为 9 个簇（图 7B、C）。根据 CellMarker 标记的表达模式，作者将这些细胞簇人工标注为以下 9 种细胞类型：1) 天真 CD4+ T 细胞（表达 IL7R 和 CCR7）；2) CD14+ 单核细胞（表达 CD14 和 LYZ）；3) 记忆 CD4+ 细胞（表达 IL7R 和 S100A4）；4) B 细胞（表达 MS4A1）；5) CD8+ T 细胞（表达 CD8A）；6) 单核细胞（表达 FCGR3A、MS4A7 和 FCGR3A）；7) NK 细胞（表达 GNLY 和 NKG7）；8) DC 细胞（表达 FCER1A 和 CST3）；9) 血小板细胞（表达 PPBP）。然后，作者开始研究UVM中壁细胞与免疫细胞之间的交流方式。作者使用 CellcRhat 软件包，利用现有数据库预测蛋白质与蛋白质之间的相互作用，从而了解细胞内的相互作用。净图显示，在所有细胞类型中，CD4+ T 细胞和 CD8+ T 细胞的相互作用最多（图 7D）。作者利用单细胞测序数据研究了建模基因在不同细胞类型中的表达情况。MMP9 主要在 NK 细胞中表达，SEMA3B 主要在 Naive CD4+ T 细胞和血小板细胞中表达，S100A13 主要在 Naive CD4+ T 细胞、记忆 CD4+ T 细胞和 CD8+ T 细胞中表达（图 7E, F ）。最后，作者研究了特定信号通路与细胞集群之间的相关性。作者发现 ROS 通路主要富集在 NK 细胞和血小板细胞中（图 7G）。
图7
8.体外和体内实验验证
为了进一步证实MMP9、S100A13和SEMA3B的表达模式，作者分析了MMP9、S100A13和SEMA3B在UVM细胞系（C918）和正常成人视网膜色素上皮细胞系（ARPE-19；附图1A）中的mRNA和蛋白表达水平。结果显示，MMP9、S100A13 和 SEMA3B 的 mRNA 表达水平在 C918 细胞系中升高。作者还研究了 UVM 细胞系与正常细胞系之间 MMP9、S100A13 和 SEMA3B 蛋白表达水平的差异。为了证实某些基因在 C918 细胞增殖和迁移中的作用，作者进一步对预后标志物进行了实验研究。由于 S100A13 基因在早期建立的模型中具有很高的系数值和抗逆性，因此作者在额外的实验中研究了 S100A13 的致癌作用。随后，作者进行了功能缺失实验，在 C918 细胞中沉默 S100A13，以研究 S100A13 在紫外线转移过程中的生物学作用。qRT-PCR 和 WB 验证了 siRNA 在 C918 细胞中的敲除效果，由于 sii-S100A13-3 的基因沉默效果最好，因此作者选择了 sii-S100A13-3 进行接下来的测试（图 8A，B）。CCK-8 检测显示，S100A13 基因敲除后，细胞活力在 72 小时内显著降低（图 8C）。EdU 染色试验也被用于评估细胞增殖。与 siRNA 阴性对照组相比，S100A13 表达减少的细胞的 EdU 阳性细胞数量明显减少（图 8D）。此外，经孔侵袭和划痕实验结果表明，敲除 S100A13 后，UVM 细胞的迁移和侵袭能力大大降低（图 8E、F）。与CCK-8结果一致，各组的活/死染色检测结果显示，敲除S100A13后，C918细胞的死亡数量急剧增加（图8G）。这些结果表明，S100A13 基因可能在 UVM 细胞存活过程中发挥了重要作用。根据 GSEA 结果，S100A13 的表达与 ROS 通路特征显著正相关（图 8H）。氧化应激在紫外线损伤中起着重要作用。荧光染色显示，S100A13 基因敲除后，C918 细胞中的 ROS 生成明显增加（图 8I）。同时，荧光染色发现，与对照组相比，S100A13敲除后的C918细胞中DHE意外上调（图8J）。qRT-PCT显示，与对照组相比，S100A13敲除的C918细胞中Nrf2和NQO1表达上调（图8K），而Western印迹检测到的HO-1表达下调（图8L）。这些结果表明，S100A13 可能通过缓解氧化应激在紫外线损伤进展过程中发挥作用。
图8
为了确定在湘雅医院眼科中心接受治疗的患者的 UVM 组织和邻近非肿瘤组织中具有代表性的 S100A13 蛋白水平，对其进行了免疫组化（IHC）标记。与之前的研究结果一致，UVM 组织的 S100A13 蛋白水平远高于邻近的非肿瘤组织（图 8M，N）。
总结
作者的研究通过系列生物信息学策略确定了免疫相关基因预后特征，它不仅在三个相对较大的独立患者队列中得到了验证，而且在预测患者生存结果方面具有显著的临床价值和可靠的表现。作者还构建了主要免疫相关基因-TFs 网络，以揭示免疫基因调控的可能机制。此外，作者还建立了多变量 Cox 预测模型，以确定预测紫外线乳腺癌预后的独立因素。新的免疫相关基因预后特征为紫癜患者的预后评估提供了依据。
