基于单细胞数据研究泛癌相关坏死性凋亡的生物标志物

Biomarkers of Pan-Cancer-Associated Necroptosis Studied Based on Single-Cell Data

摘要: 目的：基于单细胞测序数据筛选泛癌相关坏死性凋亡的生物标志物，阐明生物标志物对肿瘤的预后和免疫治疗的价值。方法：采用高通量基因表达数据库中下载的23种癌症的单细胞转录组测序数据，基于基因名片数据库中下载的坏死性凋亡基因集对单细胞数据打分，筛选差异表达基因。通过机器学习模型筛选差异基因得到泛癌相关坏死性凋亡的基因集。采用TCGA数据库下载的33种癌症的转录组数据和临床数据，利用基因集变异分析对基因集打分。基于GO、KEGG和GSEA分析基因集的功能和通路、使用COX分析、Kaplan-Meier分析和深度生存模型分析生物标志物的预后价值，并利用深度学习生存分析模型验证基因集预后效果。使用cBioPortal数据库分析免疫浸润相关性。结果：通过单细胞数据筛选出包含34个基因的生物标志物。基因集在大多数肿瘤中表达升高，在坏死性凋亡通路上高表达，并且在葡萄膜黑色素瘤、胰腺腺癌、间皮瘤等多种癌症中与较差的预后相关。基因集的评分与多种免疫细胞浸润相关，包括巨噬细胞、T细胞。结论：基因集在多种肿瘤中高表达，并于预后不良相关，且与免疫细胞密切相关。基因集可能作为一种生物标志物和免疫治疗效果的预测因子，成为有潜力的新的治疗靶点。

Abstract: Objective: To screen pan-cancer-associated necroptosis biomarkers based on single-cell sequencing data, and to elucidate the value of biomarkers for tumor prognosis and immunotherapy. Methods: Single-cell transcriptome sequencing data of 23 cancers downloaded from a high-throughput gene expression database were used to screen differentially expressed genes by scoring the single-cell data based on the necroptosis gene set downloaded from the gene business card database. The gene set of pan-cancer-associated necroptosis was obtained by screening differential genes through machine learning model. Transcriptomic data and clinical data of 33 cancers downloaded from TCGA database were used to score the gene set using gene set variant analysis. The prognostic value of the biomarkers was analyzed based on GO, KEGG, and GSEA analysis of the gene set’s functions and pathways, using COX analysis, Kaplan-Meier analysis, and deep survival modeling, and the prognostic effect of the gene set was validated using a deep learning survival analysis model. Immune infiltration correlations were analyzed using the cBioPortal database. Results: Biomarkers containing 34 genes were screened by single-cell data. Gene sets had elevated expression in most tumors, were highly expressed on the necroptosis pathway, and were associated with poorer prognosis in several cancers, including uveal melanoma, pancreatic adenocarcinoma, and mesothelioma. Gene set scores were associated with infiltration of a variety of immune cells, including macrophages and T cells. Conclusion: Gene sets are highly expressed in a variety of tumors and are associated with poor prognosis and are strongly associated with immune cells. Gene sets may be potential new therapeutic targets as a biomarker and predictor of immunotherapy efficacy.

1. 引言

程序性细胞死亡(Programmed cell death, PCD)途径可清除功能上可有可无的、受感染的细胞或潜在的肿瘤细胞，这凸显了它们在体内平衡、宿主抵御病原体、癌症和一系列其他病症中的重要作用。目前几种程序性细胞死亡途径已经被研究，主要包括细胞凋亡、细胞焦亡、坏死性凋亡等[1]，而坏死性凋亡被认为在癌症生物学的调控中的主要细胞死亡过程，包括肿瘤发生、癌症转移、癌症免疫和癌症亚型[2]。坏死性凋亡是当细胞凋亡受阻时，通过细胞外信号或细胞内信号被激活的细胞自我破坏的过程[3]。坏死性凋亡通常由细胞外刺激触发，包括由RIG-I样受体、TLRs和死亡受体(Death Receptors, DR)刺激启动的信号通路，配体(如TNFα)与细胞膜中的DR结合后，这些受体被激活后与衔接蛋白TRADD和TRAF2结合，从而导致RIP激酶的下游活化[4]。RIPK1激活RIPK3的磷酸化导致MLKL的募集和随后的磷酸化，并转移到质膜形成孔复合物。细胞膜破裂后细胞内容物溢出进入器官，导致炎症表型的出现和损伤相关分子模式(damage associated molecular patterns, DAMPs)的释放，如IL-1α、IL-β和IL-33等，从而引发免疫应答，最终导致细胞死亡。与传统的细胞凋亡不同，坏死性凋亡具有一定的促炎作用，会促进肿瘤的发展[5]。肿瘤细胞坏死性凋亡可调节肿瘤微环境(tumor microenvironment, TME)，促进肿瘤的发生和转移[6]。然而，关于坏死性凋亡与肿瘤之间相关性的进一步研究表明，这些特征并不是绝对的。RIPK1和RIPK3在胰腺导管腺癌(pancreatic ductal adenocarcinoma, PDA)中均高度表达。研究表明，体内RIPK3 缺乏或 RIPK1抑制可延缓小鼠胰腺导管腺癌的进展。这种表型与抗肿瘤免疫反应增强有关，表现为RIPK3缺乏的PDA中淋巴细胞浸润增加，免疫抑制的髓细胞浸润减少。由于发生坏死性凋亡的肿瘤细胞可能会释放促进瘤周免疫抑制的可溶性因子，因此在RIPK3缺失的PDA中发现了两种细胞因子CXCL1和SAP130下调。当RIPK3介导的肿瘤坏死性凋亡发生时，释放的可溶性因子与炎症细胞上的受体(如SAP130及其同源受体Mincle)结合，诱导免疫抑制，促进PDA的进展[7]。Ruan等[8]的研究表明MLKL低表达与原发性宫颈鳞状细胞癌患者的预后不良显著相关。由于MLKL是坏死性凋亡信号通路的关键介质，MLKL的低表达可能暗示患者坏死性凋亡信号通路的降低[9]。在宫颈鳞状细胞癌患者中，低MLKL表达与预后不良相关的一个潜在机制可能是这些患者的坏死性凋亡信号通路减少的结果[10]。肿瘤类型的不同，坏死对肿瘤产生的影响也不同。深入研究坏死性凋亡与不同癌症之间的相互作用有助于更好地理解癌症的本质，并为治疗提供新的策略。

单细胞转录组测序(Single-cell RNA sequencing, scRNA-seq)是解决上述问题的方法之一[11]。scRNA-seq使用从人体获得的样本进行分析，可以代表体内实验的结果。它是一种高分辨率的工具，克服了传统体测序的局限性，使我们能够在单个细胞水平上测量肿瘤的异质性，深入了解细胞的功能、特性和相互作用等方面的信息[12]。通过单细胞研究，可以对细胞在正常和疾病状态下的差异进行详细分析，揭示出更精确的细胞类型和组织发育过程的信息，对TME的研究很有价值。使用scRNA-seq可以用来评估坏死性凋亡和泛癌之间的关系。

肿瘤坏死性凋亡的具体机制尚不清楚，在不同的肿瘤中观察到相反的效果。可以基于泛癌来分析坏死性凋亡在癌症中的作用[13]。本文使用来自GEO数据库中泛癌的scRNA-seq数据。基于Genecards数据库中的636个坏死相关基因，筛选泛癌中的坏死型凋亡的标志基因。利用机器学习模型对坏死性凋亡基因进行筛选得到包含34个基因的坏死性凋亡基因集(NCPS.Sig)。针对NCPS.Sig进行了泛癌分析，研究NCPS.Sig的表达与TCGA数据库中33种肿瘤类型的预后之间的相关性，并通过通路富集分析探讨其可能的分子机制。此外，还研究了NCPS.Sig与肿瘤突变负担(tumor mutation burden, TMB)、微卫星不稳定性(microsatellite instability, MSI)与肿瘤的免疫学以及预后相关性。本研究为临床癌症的治疗提出了新的思路。

2. 材料与方法

2.1. 泛癌数据采集和预处理

从高通量基因表达数据库(Gene Expression Omnibus data base, GEO)数据库(https://www.ncbi.nlm.nih.gov/geo/)中下载了23种癌症的单细胞转录组测序(Single-cell RNA sequencing, scRNA-seq)数据集，包括头颈部腺样囊性癌(Adenoid Cystcarcinoma of the Head and Neck, ACC)胶质母细胞瘤(Glioblastoma, GBM)、非小细胞肺癌(Non-small cell lung cancer, NCSLS)、结肠癌(Colorectal cancer, COAD)、肾透明细胞癌(Kidney renal clear cell carcinoma, KIRC)、乳腺癌(Breast invasive carcinoma, BRCA)、肺腺癌(Lung adenocarcinoma, LUAD)、卵巢癌(Ovarian cancer, OV)、甲状腺癌(Thyroid carcinoma, THCA)、胃癌(Stomach adenocarcinoma, STAD)、食管鳞状细胞癌(Esophageal carcinoma, ESCA)、胰腺癌(Pancreatic adenocarcinoma, PAAD)、膀胱癌(Bladder Urothelial Carcinoma, BLCA)、肾脏嫌色细胞癌(Kidney Chromophobe, KICH)、口腔鳞状细胞癌(Oral squamous cell carcinoma, OSCC)、肾乳头状细胞癌(Kidney renal papillary cell carcinoma, KIRP)、肉瘤(Sarcomav, SARC)、子宫内膜样癌(Uterine Corpus Endometrial Carcinoma, UCEC)、睾丸癌(Testicular Germ Cell Tumors, TGCT)、胆管癌(Cholangiocarcinoma, CHOL)、肝细胞癌(Liver hepatocellular carcinoma, LIHC)、直肠癌(Rectum adenocarcinoma, READ)、弥漫性大B细胞淋巴瘤(Lymphoid Neoplasm Diffuse Large B-cell Lymphoma, DLBC)。33种癌症类型的The Cancer Genome Atlas (TCGA)转录组数据均从UCSC XENA网站(https://xenabrowser.net/datapages/)下载。仅选择有总生存期数据的患者进行进一步的数据分析。从Genotype-Tissue Expression(GTEx)数据库下载健康样本的转录组数据。利用R软件的Seurat包进行单细胞分析，排除基因表达数小于50的细胞以及线粒体基因表达占比大于30%的细胞。对质控后的细胞进行归一化和标准化。然后进行主成分分析(principal components analysis, PCA)，利用均匀流形近似和投影(Uniform Manifold Approximation and Projection, UMAP)对细胞聚类进行可视化。根据已知得细胞类型标记基因对每个细胞簇进行细胞注释。

2.2. 获取坏死性凋亡相关基因

在GeneCards数据库(https://www.genecards.org/)中，鉴定出636个与坏死性凋亡相关的基因。将相关系数阈值设置为大于或等于1.0，得到102个与坏死性凋亡相关的基因。

2.3. 筛选癌症相关坏死性凋亡基因

使用AddModuleScore函数基于Genecards下载的坏死性凋亡基因相关性大于1的102个基因进行打分，根据得分的平均数，将细胞分为高分组和低分组。利用FindMarkers函数筛选log2FoldChange > 0.5的差异表达基因。将筛选出的标志基因与已知的坏死性凋亡基因取交集，得到每个癌症的坏死性凋亡相关基因。

2.4. 基于机器学习筛选生物标志物

将每个scRNA-seq数据坏死性凋亡高分组和低分组细胞进行随机抽样，每个肿瘤数据高分组和低分组各随机抽取125个样本，合并数据共计2875个样本。利用R的sva包的ComBat函数将样本进行去批次，随后拆分数据为训练集(n = 4600, 80%)和测试集(n = 1150, 20%)。利用泛癌的scRNA-seq数据相关坏死性凋亡的基因特征构建预测模型，分别应用三种机器学习算法进行模型构建。三种机器学习算法是支持向量机递归特征消除(Support Vector Machine-Recursive Feature Elimination，SVM-RFE，R包e1071)、基于Boruta特征选择算法的随机森林(Random Forest，RF，R包caret)、最小绝对收敛和选择算子(Least Absolute Shrinkage and Selection Operator，LASSO，R包glmnet)。

Boruta特征选择算法与RF结合进行标志物筛选。Boruta通过比较原始特征和其随机变量之间的表现来确定对目标变量有显著影响的特征。利用Boruta筛选出的特征构建RF模型，RF可以评估每个特征在构建决策树时的重要性筛选特征。LASSO通过对目标函数添加L1正则化项，倾向于产生稀疏权重向量，从而使得部分特征的权重变为零，进行特征选择。SVM-RFE是基于SVM的特征选择方法，通过RFE的方式来选择最重要的特征子集。SVM通过找到一个最优的超平面来将不同类别的数据分隔开，同时最大化间隔，以提高泛化能力。使用SVM对所有特征进行建模评估，得到每个特征的权重，根据特征的重要性，逐步剔除权重较小的特征，重新构建模型，直到达到设定的特征数量或停止条件，得到最佳的特征子集。

基于5折交叉验证，通过ROC曲线下的面积(Area Under Curve, AUC)指标评估机器学习模型性能。AUC是ROC曲线下的面积，ROC曲线横轴表示假阳率(False Positive Rate, FPR)，纵轴表示真阳率(True Positive Rate, TPR)。ROC曲线展示了在不同分类阈值下，模型的灵敏度和特异度之间的权衡关系。AUC值越接近1，说明模型性能越好，具有更好的分类能力；AUC值为0.5时，则表示模型的分类效果等同于随机猜测。AUC作为评价二分类模型性能的重要指标，在实际应用中被广泛使用。将三种机器学习模型筛选结果取交集，得到泛癌相关坏死性凋亡的生物标志物。

2.5. NCPS评分的构建和验证

为进行NCPS.Sig基因集通路富集分析，基于TCGA泛癌普通转录组数据(bulk RNA-seq, bulk-seq)的归一化基因表达数据，利用R的GSVA包进行打分。基因集变异分析(Gene Set Variation Analysis, GSVA)是一种用于基因表达数据分析的方法，可用于评估基因集在不同样本中的得分。通过将基因表达数据转化为基因集水平的表达变化，可以更好地捕捉基因集的生物学含义和功能，帮助理解基因调控网络、信号通路等方面的生物学过程。通过基因表达数据计算基因集的得分，对得到的基因集得分进行统计分析，比较不同样本之间基因集的得分差异。坏死性凋亡基因特征NCPS.Sig的GSVA得分根据中位数，将泛癌bulk-seq数据集中的样本分为高坏死性凋亡肿瘤样本和低坏死性凋亡肿瘤样本。评估NCPS.Sig的GSVA得分在正常样本和泛癌肿瘤样本之间的表达差异。

2.6. 通路和功能富集分析

基因本体论(Gene Ontology, GO)是一种用于描述基因和蛋白质功能的标准化注释体系，将生物学过程(biological processes, BP)、细胞组分(cellular components, CC)和分子功能(molecular functions, MF)分为不同的术语和层级结构。通过R的cluster Profiler包对NCPS.Sig基因集进行功能富集分析，可视化p值小于0.05的相关术语。随后使用基因集富集分析(Gene Set Enrichment Analysis, GSEA)分析MsigDB下载的坏死性凋亡相关通路上基因集的表达情况。

2.7. 免疫细胞浸润的评估

TMB反映了肿瘤细胞基因组中的突变数量，与肿瘤的遗传多样性和进化速度相关[14]。高TMB的肿瘤通常伴随着更多的新抗原的产生，会被免疫系统识别并引发免疫反应，高TMB可能与较好的临床预后相关。MSI是肿瘤组织中DNA错配修复功能缺陷的结果，MSI高的肿瘤往往预示着更好的治疗效果[15]。TMB和MSI数据来自GDC (Genomic Data Commons, https://portal.gdc.cancer.gov/)数据库，利用TCGAbiolinks包下载泛癌的TMB数据，进行TMB与NCPS.Sig打分的相关性分析。MSI数据来自cBioPortal数据库(https://www.cbioportal.org/)预测与NCPS.Sig得分的相关性。对于免疫细胞浸润的评估，根据TCGA肿瘤队列的转录组数据，使用CIBERSORT算法来评估NCPS.Sig评分与TCGA队列中不同癌症类型的22个免疫细胞浸润的相关性。

2.8. 生存分析

利用NCPS.Sig评分和来自TCGA数据库的泛癌症临床数据，通过单变量Cox分析研究了NCPS.Sig与肿瘤预后之间的关系。计算了相对p值和危险比，并绘制了森林图。随后，采用R包survminer针对基因集评分生成了各种类型肿瘤的Kaplan-Meier生存曲线。

2.9. 深度生存分析模型

利用深度神经网络Cox模型DeepSurv分析生存数据[14]，验证NCPS.Sig在预后的预测效果。DeepSurv是一种深度前馈神经网络，具有多个隐藏层并加入了整流线性单位(ReLU)、批处理归一化、dropout层、权重衰减正则化等技术。隐藏层由全连接的节点层和dropout交叠组成，以增强模型的泛化性能和避免过拟合。加入ReLU非线性激活函数用于学习复杂的模式和特征，提高模型非线性拟合能力，它的形式为：

$f\left(x\right)=\{\begin{array}{c}x,x>0\\ 0,x\le 0\end{array}$ (1)

整合输出后得到风险率${\widehat{h}}_{\theta }\left(x_i\right)$ ，用于估计Cox模型中的交叉熵损失函数：

$l\left(\theta \right):=-\frac{1}{N_{E=1}}{\displaystyle \sum _{i:E_i=1}\left({\widehat{h}}_{\theta }\left(x_i\right)-\log {\displaystyle \sum _{j\in R\left(T_i\right)}{e}^{{\widehat{h}}_{\theta }\left(x_j\right)}}\right)}+\lambda \cdot {\Vert \theta \Vert }_2^2$ (2)

其中，$N_{E=1}$ 是患者样本数，$E$ 为生存状态，$T$ 为生存时间，$x$ 是特征变量，$\lambda$ 是L2正则化的参数，通过参数权重$\theta$ 预测患者基因表达量对风险率的影响。

将K-M生存分析显著的TCGA数据库的肿瘤类型包括UVM、LGG、PRAD、LIHC、LUAD、KIRC、HNSC、MESO的基因表达数据和临床数据输入DeepSurv模型，共4360个样本。将样本分为训练集(80%，3488个样本)和测试集(20%，872个样本)进行预测。参数设置3层隐藏层，每层由全连接层、Batch normalization、激活层(ReLU)、dropout层堆叠组成，学习率为1e−2，L2正则化参数为1e−4。采用适应性矩估计(Adaptive Moment Estimation, Adam)算法，根据每个参数的梯度情况来动态调整学习率，它维护了每个参数两个额外的变量：一阶矩估计和二阶矩估计，用来估计梯度的均值和方差。同时，采用Cox比例风险回归模型和随机森林生存分析模型与DeepSurv进行对比，Cox比例风险模型基于风险比例假设，即各个预测因子对风险的影响在时间上是恒定的。Cox风险回归模型如下：

$h\left(t|x_i\right)=h_0\left(t\right)\exp \left(x_i\beta \right)$ (3)

$x=\left(x_{i1},x_{i2},\cdots ,x_{in}\right)$ 代表预测因子，表示可能与生存时间相关的协变量。$h_0\left(t\right)$ 表示当所有危险因素均为0时的基线风险率。$\beta$ 表示Cox模型的偏回归系数，需根据实际样本数据进行估计。在Cox比例风险回归中，通常通过部分似然估计来预测模型的参数$\beta$ ，从而推断不同因素对生存时间的影响。随机森林生存分析模型是一种结合了随机森林和生存分析的方法，用于处理生存数据和预测生存时间。每棵决策树都会针对生存数据进行训练，考虑个体的生存时间和相关因素。利用AUC值比较模型性能，并验证NCPS.Sig基因集在预后的作用。

3. 结果

3.1. 肿瘤环境中细胞群的坏死性凋亡基因状态

假设泛癌中与坏死性凋亡的转录组变化是癌症的生存预后和免疫治疗靶点。为了证明，下载了23种肿瘤的scRNA-seq数据，以及坏死性凋亡相关基因集NCPS。将scRNA-seq数据质控后，进行手动细胞注释，并通过UMAP可视化，如图1(a)~(d)所示。利用NCPS基因集评估单细胞肿瘤数据，将细胞分为高表达坏死性凋亡基因组High和低表达坏死性凋亡组Low。随后基于AddModuleScore打分分析各细胞簇得分显示，T细胞得分高于单核细胞、B细胞、巨噬细胞、上皮细胞等细胞簇。因此t细胞表现出最高的坏死性凋亡水平。

3.2. 坏死性凋亡相关生物标志物的鉴定

从Genecards数据库下载的坏死性凋亡相关的基因集，共636个基因，筛选出相关性大于1的102个基因命名为NECR.Sig，对细胞进行打分。基于打分将细胞分为高低组，进行差异基因分析，将各肿瘤数据筛选出的差异表达基因与NECR.Sig取交集得到NECR.UP.Sig基因集，共58个基因。将NECR.UP.Sig基因集输入机器学习模型，如图2所示，LASSO回归模型从NECR.UP.Sig鉴别出46个基因，经过5折交叉验证，AUC值为0.74。通过Boruta特征选择后的随机森林模型识别出44个基因，AUC值为0.78。SVM-RFE识别出55个相关基因。将三个模型鉴别的基因通过韦恩图取交集，得到包含34个基因的泛癌相关坏死性凋亡的基因集，命名为NCPS.Sig (图2)。

Figure 1. Expression of necroptosis genes in scRNA-seq pan-cancer data. (a) Cellular annotation visualization of NSCLC, downscaled plots of high and low groupings after scoring of necroptosis genes, and mountain range plots of different cell cluster scores. (b) Cellular annotation visualization of COAD, downscaled plots of high and low groupings after scoring of necroptosis genes, and maps of different cell cluster scores mountains. (c) Cellular annotation visualization of ESCA, downscaled plots of high and low groupings after scoring of necroptosis genes, and maps of different cell cluster scores mountains. (d) Cellular annotation visualization of THCA, downscaled plots of high and low grouping after necroptosis gene scoring, and different cell cluster scoring mountain range plots

图1. scRNA-seq泛癌数据中坏死性凋亡基因的表达情况。(a) NSCLC的细胞注释可视化，坏死性凋亡基因打分后高低分组的降维图，以及不同细胞簇得分山峦图。(b) COAD的细胞注释可视化，坏死性凋亡基因打分后高低分组的降维图，以及不同细胞簇得分山峦图。(c) ESCA的细胞注释可视化，坏死性凋亡基因打分后高低分组的降维图，以及不同细胞簇得分山峦图。(d) THCA的细胞注释可视化，坏死性凋亡基因打分后高低分组的降维图，以及不同细胞簇得分山峦图

Figure 2. Machine learning model to identify pan-cancer and necroptosis-related genes. (a)~(c) LASSO regression model screening features with model AUC value of 0.74; (d)~(f) Random forest model based on Boruta feature selection with model AUC value of 0.78; (g) SVM-RFE model screening 56 features; (h) Venn diagram of LASSO, RF and SVM-RFE model screening features taken as intersection

图2. 机器学习模型识别泛癌与坏死性凋亡相关基因。(a)~(c) LASSO回归模型筛选特征，模型AUC值为0.74；(d)~(f) 基于Boruta特征选择的随机森林模型，模型AUC值为0.78；(g) SVM-RFE模型筛选出56个特征；(h) 将LASSO、RF和SVM-RFE模型筛选特征取交集的韦恩图

3.3. 在泛癌症水平上的坏死性凋亡景观的描述

为了探讨泛癌中NCPS.Sig的生物学特性。利用GSVA评估NCPS.Sig，计算TCGA数据库中33种癌症的每个患者的转录组特征得分。随后对GTEx数据库中的正常样本进行打分，分析正常样本和肿瘤样本的基因表达差异。由图3所示，肿瘤样本中GSVA高分组比例较高，正常样本中低分组占比更高。图3(a)显示，在33种癌症数据中，共有21个肿瘤数据正常样本和肿瘤样本基因表达差异显著，并且其中19个NCPS.Sig评分肿瘤样本都高于正常样本，包括HNSC、PAAD、LIHC等肿瘤类型，但LUSC和KICH的正常样本评分高于肿瘤样本。

对NCPS.SIG基因集进行功能通路富集分析，GO富集分析显示基因集在多种生物学过程、细胞组分和分子功能中富集，包括模式识别受体、调节凋亡信号通路、程序性细胞凋亡。由图4所示，NCPS.Sig中多个基因都参与程序性细胞凋亡和调节凋亡信号通路的生物学过程。

3.4. 肿瘤微环境中NCPS.Sig评分的免疫特性

高NCPS.Sig GSVA评分与肿瘤免疫浸润的相关性，这种转录组特征也增加了预测抗PD-L1/PD-1免疫检查点抑制剂治疗反应的可能性。TMB和MSI是免疫治疗反应的一个强有力的泛癌预测因子[16]。我

Figure 3. Differential expression of necroptosis genes in normal and cancer samples. (a) Differential expression violin plots of normal and cancer samples in 33 cancer types; (b) Heatmap of differential GSVA scores in normal and tumor samples

图3. 坏死性凋亡基因在正常样本和癌症样本的表达差异。(a) 正常样本和癌症样本在33种癌症类型种的差异表达小提琴图；(b) 正常样本和肿瘤样本的GSVA评分差异热图

Figure 4. NCPS.Sig gene set enrichment analysis. (a-b) GO enrichment analysis; (c) GSEA analysis of gene sets in the apoptotic pathway

图4. NCPS.Sig基因集富集分析。(a-b) GO富集分析；(c) GSEA分析基因集在凋亡通路的表现

们评估了NCPS.Sig与TMB之间的相关性。11种肿瘤与TMB显著相关，并且THYM、STAD、PAAD等7种癌症呈正相关。STAD与MSI也呈正相关。随后，使用CIBERSORT来评估GSVA评分与TCGA队列中不同癌症类型的22个免疫细胞浸润的相关性。图5显示，较高的GSVA评分与免疫细胞明显改变的肿瘤浸润有关。特别是，几乎所有癌症类型的GSVA评分与M2巨噬细胞肿瘤浸润呈正相关，与T细胞滤泡辅助性细胞肿瘤浸润呈负相关。因此，我们认为坏死性凋亡相关的转录组特征，可作为免疫细胞反应失衡的泛癌预后因子。

Figure 5. Tumor immune properties correlated with necroptosis score. (a-b) Radar plots of the correlation of NCPS.Sig score with TMB and MSI; (c) Correlation of NCPS.Sig score with 22 immune cells in pan-cancer. Red color represents positive correlation, blue color represents negative correlation, and larger points represent higher correlation. p-value < 0.05 is *, p-value < 0.01 is **, and p-value < 0.001 is ***

图5. 与坏死性凋亡评分相关的肿瘤免疫特性。(a-b) NCPS.Sig评分与TMB和MSI的相关性雷达图；(c) NCPS.Sig得分与泛癌的22种免疫细胞相关性。红色代表正相关，蓝色代表负相关，点越大代表相关性越高。p值 < 0.05为*, p值 < 0.01为**，p值 < 0.001为***

3.5. 生存分析

利用NCPS.Sig基因在TCGA泛癌数据中使用Cox比例风险回归模型(图6(a))。结果显示，共有10种癌症的GSVA评分与预后显著相关。UVM、THYM等9种癌症显示高评分的样本与不良预后相关。还可以观察到，NCPS.Sig评分对包括HNSC、KIRC、LGG、LIHC、LUAD、MESO、PAAD、UVM在内的各种其他癌症类型的总生存率具有很强的预后能力(图6(b)~(i))，表明我们选择的基因集评分在各种癌症中是可靠的预后指标。

Figure 6. Prognostic value of the NCPS.Sig score. (a) Forest plot of pan-cancer COX analysis of NCPS.Sig. (b)~(i) Survival analysis plot of NCPS.SIG scores in pan-cancer

图6. 坏死性凋亡评分的预后价值。(a) NCPS.Sig的泛癌COX分析森林图。(b)~(i) NCPS.Sig打分在泛癌的生存分析曲线图

3.6. 与坏死性凋亡相关的泛癌预后模型的构建和验证

为了验证所获得的NCPS.Sig基因特征对预测泛癌患者预后风险的能力，在8种癌症的TCGA队列上构建了一个DNN分类模型DeepSurv。在该模型中，以NCPS.Sig基因集作为输入，以生存状态作为分类标签。测试集模型(TCGA队列)的AUC为0.894。此外，在随机生存森林和COX比例风险回归模型中验证基因特征的可靠性，AUC值分别为0.894和0.649 (图7)。结果表明基因集可作为预测患者预后风险的有效生物标志物。

Figure 7. Necroptosis-associated pan-cancer prognostic modeling. ROC plots of the NCPS.Sig gene set in the DeepSurv model, the RF survival model, and the COX risk regression model

图7. 坏死性凋亡相关的泛癌预后模型。NCPS.Sig基因集在DeepSurv模型、RF生存模型以及COX风险回归模型的ROC曲线图

4. 讨论

坏死性凋亡在肿瘤发生和进展中的作用是多方面的。虽然许多癌症的坏死性凋亡特征及对TME的影响已被广泛探索，但泛癌中坏死性凋亡相关预后的单细胞机制仍不确定[17]。因此，本文整合单细胞和批量转录组数据，全面分析泛癌的坏死性凋亡特征。对单细胞打分可看出，T细胞高表达坏死性凋亡相关基因。由于坏死性凋亡本质上是高度免疫原性的，因此T细胞会富集坏死性凋亡相关基因[18]。然而RPA1作为坏死性凋亡关键基因，可通过RPA1耗竭激活ZBP1-RIPK3信号转导，从而导致T细胞坏死性凋亡，这种由RPA1缺陷诱导的坏死性凋亡性T细胞死亡允许释放DAMP，进而募集白细胞并加剧炎症反应[19]。随后，利用RF、SVM-RFE和LASSO回归算法研究了坏死性凋亡基因特征，筛选出包含34个基因的生物标志物。基因集中TNF肿瘤坏死因子是先天免疫和促炎反应的关键调节因子。TNF促使RIP1和RIP3形成一个促坏死样凋亡复合物，进而招募并激活下游效应蛋白MLKL，使MLKL转位到质膜上，诱导细胞死亡[20]。而MYC抑制TNF-α诱导的坏死性凋亡，阻止RIPK1-RIPK3复合物的形成[21] [22]。而基因集中CTSB被内质网应激，促进细胞坏死性凋亡，从而通过PKCα-JNK-cJun通路诱导TNFα分泌[23]。基于由34个基因组成的坏死性凋亡特征，构建了一个针对泛癌患者的评分系统。利用GSVA对TCGA和GTEx样本进行打分，分为高低坏死性凋亡得分组。并对评分进行验证，富集分析结果显示，高NCPS.Sig得分组不仅富集到regulation of apoptotic、programmed Necrotic等凋亡相关通路，且富集到多条免疫相关通路包括regulation of Inflammatory response等通路[24]。NCPS.Sig评分与TMB、MSI相关性分析表明，基因集可作为肿瘤免疫的治疗靶点[25]。利用COX分析基因集在癌症预后的作用。高分组的患者具有更广泛的坏死性凋亡特征和较低的生存率。基因集评分在预测患者生存方面也表现良好，表明基因集在各种癌症中是可靠的预后指标。K-M生存曲线看出，HNSC、LIHC与LGG等癌症中基因集高评分与预后不良相关，但KIRC的高评分与预后改善相关。Urbano等人研究表明KIRC的细胞表现出高水平的RIPK1和RIPK3，它们更容易受到TNFR1引发的坏死性凋亡的影响，坏死性凋亡途径参与肿瘤坏死，导致肿瘤自发性坏死性凋亡[20] [26]，这与我们的研究结果一致。建立了一个基于深度学习的生存模型，证实了该基因特征作为预测患者预后的生物标志物的有效性。坏死性凋亡基因集可作为一种生物标志物和免疫治疗效果的预测因子和预后指标，成为有潜力的新的治疗靶点。然而，由于目前能收集到的scRNA-seq数据有限，后续的研究应该结合多组学数据，进一步挖掘坏死性凋亡在癌症的分子机制。同时，后续的研究应该进一步结合相关临床数据和实验对上述研究结果进行验证。

NOTES

^*第一作者。

^#通讯作者。