# paper.pdf ## Metadata - PDFFormatVersion=1.4 - IsLinearized=true - IsAcroFormPresent=false - IsXFAPresent=false - IsCollectionPresent=false - IsSignaturesPresent=false - CreationDate=D:20260417121440+05'30' - Creator=Adobe InDesign CS6 (Windows) - ModDate=D:20260417121541+05'30' - Producer=Adobe PDF Library 10.0.1 - Trapped=False - xmp:createdate=2026-04-17T12:14:40+05:30 - xmp:metadatadate=2026-04-17T12:15:41+05:30 - xmp:modifydate=2026-04-17T12:15:41+05:30 - xmp:creatortool=Adobe InDesign CS6 (Windows) - xmpmm:instanceid=uuid:caf0aef1-cf1d-4413-9070-ae90a1e1ecfd - xmpmm:originaldocumentid=xmp.did:B1A72632D5D1E311A39EEF241A36A0B8 - xmpmm:documentid=xmp.id:169EC1A0283AF1118599D61EFAB2D8DC - xmpmm:renditionclass=proof:pdf - xmpmm:derivedfrom=xmp.iid:4DDF3AE4273AF1118E62BF250972158Dxmp.did:12AA28B591A8EE11A7769276598EFA6Exmp.did:B1A72632D5D1E311A39EEF241A36A0B8default - xmpmm:history=convertedfrom application/x-indesign to application/pdfAdobe InDesign CS6 (Windows)/2026-04-17T12:14:40+05:30 - dc:format=application/pdf - pdf:producer=Adobe PDF Library 10.0.1 - pdf:trapped=False ## Contents ### Page 1 ReseaRch aRticlescSurvival: Single-Cell Survival Analysis of Clinical Cancer Cohort Data at Cellular ResolutionTao Ren1,2, Faming Zhao1,2, Canping Chen1,2, Le Zhou1,2, Ling-Yun Wu3,4, Gordon B. Mills2,5, Lisa M. Coussens2,6, and Zheng Xia1,2,7 ### Page 2 AACRJournals.org932 | CANCER DISCOVERY MAY 20261Biomedical Engineering Department, Oregon Health & Science University, Portland, Oregon. 2Knight Cancer Institute, Oregon Health & Science Uni- versity, Portland, Oregon. 3Academy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing, China. 4School of Mathematical Sciences, University of Chinese Academy of Sciences, Beijing, China. 5 Division of Oncological Sciences, Knight Cancer Institute, Oregon Health & Science University, Portland, Oregon. 6 Cell, Development and Cancer Biology Department, Oregon Health & Science University, Portland, Oregon. 7Center for Biomedical Data Science, Oregon Health & Science University, Portland, Oregon. T. Ren and F. Zhao contributed equally to this article.Corresponding Author: Zheng Xia, Knight Cancer Institute, Oregon Health & Science University, 3181 SW Sam Jackson Park Rd., Portland, OR 97239. E-mail: xiaz@ohsu.edu Cancer Discov 2026;16:931–52doi: 10.1158/2159-8290.CD-25-0965 This open access article is distributed under the Creative Commons Attribution- NonCommercial-NoDerivatives 4.0 International (CC BY-NC-ND 4.0) license.©2026 The Authors; Published by the American Association for Cancer ResearchintRoductionSurvival analysis has long been a cornerstone of cancer research for investigating time-to-event outcomes such as overall survival (OS) and progression-free survival (PFS). Bulk RNA sequencing (RNA-seq) studies have demonstrated the versatility of survival analysis by modeling gene expression as covariates, constructing multigene risk scores, and stratifying patients into molecular subgroups for prognostic evaluation and assessment of therapeutic efficacy (1–4). Meanwhile, single- cell RNA-seq (scRNA-seq) provides unprecedented resolution for characterizing the cellular heterogeneity of the tumor microenvironment (TME; refs. 5–8). Leveraging this capa- bility, recent advancements and the widespread adoption of scRNA-seq have led to a growing number of cohort-level cancer studies that profile single-cell transcriptomes from hundreds of patients alongside clinical survival data, offering new opportunities to elucidate how cellular heterogeneity influences disease progression and patient outcomes. Because individual cells can influence patient outcomes differently, identifying survival-associated subpopulations is crucial for improving risk assessment and developing more precise and targeted therapies beyond bulk-level analyses. In early studies, because of the limited sample size of scRNA-seq data in cancer studies, the predominant strategy for identifying risk-associated cell populations has involved transferring survival information from existing bulk RNA- seq datasets to single-cell data by integrating single-cell and bulk data, such as Scissor (9), scSurv (10), and DEGAS (11). Although these approaches have enabled the identification of cell populations associated with survival outcomes, they all rely on bulk transcriptomics and its associated survival outcomes as a bridge without employing the true survival information of the single-cell data. As cohort-level cancer studies increasingly generate scRNA- seq datasets from hundreds of patients with matched survival information, two simple strategies can be adopted to adapt traditional Cox regression frameworks to single-cell data: (i) aggregating cell-level gene expression into pseudobulk pro- files for each patient and (ii) annotating major cell types and using their proportions as patient-level covariates. Although both approaches allow the use of standard survival models, they suffer from major information loss. The pseudobulk strategy treats all cells equally and obscures the contributions of rare but prognostically important subpopulations. In con- trast, the cell type proportion strategy assumes homogeneity within annotated cell types and ignores intratype functional diversity. These methods offer a coarse approximation but fall short of realizing the full potential of single-cell data in capturing fine-grained risk structures. To our knowledge, directly performing survival prediction at single-cell resolution,Survival analysis is fundamental to cancer research. Advances in technology have enabled an increasing number of cohort-level cancer studies to incorporate single- cell sequencing alongside clinical survival data. However, no effective strategy currently exists for directly modeling survival outcomes from single-cell data. To address this gap, we present scSurvival, an attention-based multiple-instance Cox regression framework that models each tumor sample as an ensemble of cells to predict survival outcomes at both the patient and single-cell levels. To handle high dimensionality, sparsity, and batch effects, scSurvival integrates a variational autoencoder–based feature extraction module with generative modeling to enhance feature robustness and cross-batch generalizability. Comprehensive simulations demonstrate scSurvival’s superior performance and scalability. In melanoma and liver cancer single-cell RNA sequencing (scRNA-seq) cohorts, scSurvival accurately predicts patient outcomes and identifies the cell subpopulations most critical to survival. Overall, scSurvival enables robust prediction of patient survival while uncovering survival-associated cell subpopulations, advancing single-cell survival analysis in cancer research.SignifiCAnCe: Survival analysis is central to clinical oncology, yet no effective tools currently model survival outcomes directly from single-cell data. scSurvival bridges this gap by predicting patient out- comes and identifying key subpopulations from scRNA-seq with survival information, enabling scalable analyses and promoting broader adoption of cohort-level single-cell profiling in cancer research.aBstRact ### Page 3 RESEARCH ARTICLESingle-Cell Survival AnalysisMAY 2026 CANCER DISCOVERY | 933while preserving intrapatient heterogeneity and relying solely on single-cell cohort data, remains an unmet challenge, and no dedicated computational tools have yet been developed to address this need. Constructing risk prediction models such as the Cox pro- portional hazards model directly from single-cell cohort data also fundamentally differs from conventional bulk data–based Cox regression. In traditional bulk settings, each patient’s covariate values form a single feature vector that can be directly used for model construction. In contrast, single-cell expres- sion data are structured as a cell-by-gene matrix rather than a simple per-patient vector, with individual cells contributing differently to patient outcomes. Therefore, there is a clear need for a framework that can simultaneously predict patient sur- vival outcomes and identify the most informative cell subpop- ulations. Attention-based multiple-instance learning (12–14) offers a compelling solution for this task. In this framework, each tumor sample is treated as a collection of multiple instances (i.e., cells). The model learns to predict patient-level labels while simultaneously assigning weights to individual instances, thereby highlighting the cell subpopulations most critical to the survival outcomes. In this study, we conceptualize each tumor sample as an ensemble of cellular instances and thus develop a new tool, scSurvival, specifically designed to address heterogeneous cell subpopulations in survival analysis of single-cell cancer cohort data. To tackle the high dimensionality, sparsity, noise, and severe batch effects inherent in single-cell datasets, we design a feature extraction module based on variational autoen- coders (VAE; arXiv 1312.6114) and generative modeling (15, 16). Coupled with attention-based multiple-instance Cox regression (AMICR), scSurvival enables survival modeling of patient cohorts spanning multiple experimental batches. scSurvival not only integrates single-cell expression data with optional clinical covariates to construct accurate risk predic- tion models but also identifies outcome-associated cell sub- populations and characterizes their risk tendencies, thereby facilitating more refined downstream biological and clinical analyses.ResultsOverview of scSurvivalWe developed scSurvival, a single-cell survival outcome prediction tool based on attention-based multiple-instance learning, which accounts for the varying contributions of individual cells to survival outcomes. It is designed to construct survival prediction models from single-cell cancer cohort data while simultaneously identifying cell subpopulations that are strongly associated with patient risk. The core principle involves compressing each patient’s high-dimensional single- cell data into low-dimensional features, aggregating these features into patient-level representations through attention weights, and subsequently building a survival model at the patient level using Cox regression analysis. The attention mechanism preserves cellular heterogeneity within each patient, allowing subpopulations with higher attention scores to be more closely linked to survival probability. By aligning cell- level features with aggregated patient-level representations, the survival prediction model can directly infer the risk tendencies of specific cell subpopulations, thereby enabling the identification of risk-associated cells. Specifically, the main inputs to scSurvival include single-cell expression matrices measured for each patient (Fig. 1A and B) and patient survival information, i.e., survival event status and event time (Fig. 1C). Conceptually, scSurvival comprises two key modules: a cell feature extraction module (Fig. 1D) and an AMICR module (Fig. 1E). The feature extraction mod- ule is a VAE-based generative model that captures technical dropout using a zero-inflated Gaussian (ZIG) distribution (17, 18) and models the log-normalized expression matrix. LayerNorm [LN; (19)] and squeeze-and-excitation (SE; ref. 20) modules are integrated into this architecture to enable adap- tive standardization and reweighting of gene features. In addition, the feature extraction module supports batch label inputs to extract cross-batch single-cell features (Fig. 1D). The AMICR module aggregates cell-level features into patient-level representations using a multi-head attention mechanism and inputs the aggregated features into a risk scorer to perform hazard scoring and Cox regression analysis (Fig. 1E), in which the structures of the cell-level attention mechanism, gene- reweighting SE module, and risk scorer are shown in Fig. 1F. Furthermore, patient-level clinical covariates (21) can be incorporated into the risk scorer for joint analysis (Fig. 1G), enabling appropriate control of confounding influences and thereby improving the accuracy of survival prediction models and the identification of risk-associated cell subpopulations. The final loss function of the model comprises three com- ponents: the negative log-likelihood of Cox regression, the reconstruction likelihood of the ZIG-VAE, and an entropy regularization term (22) to control the sparsity of the cell-level attention weights. Training is performed in two stages: pretrain- ing and fine-tuning (Fig. 1H). In the pretraining stage, only the VAE module is trained to extract stable cell features. During the fine-tuning stage, all three loss components are jointly optimized to facilitate multiple-instance Cox regression, pro- ducing optimized cell representations and attention weights that are adaptively aligned with survival outcomes (“Methods”). The resulting outputs of scSurvival are the attention-adjusted hazard score for each cell (hazard_adj; see “Methods”, Fig. 1I), along with patient-level risk scores. This represents a funda- mental distinction from traditional survival analysis methods: scSurvival not only provides risk stratification at the patient level but also infers intrapatient cellular risk heterogeneity. Given the hazard scores of individual cells estimated by scSurvival, we can then perform downstream analyses of heterogeneous cell states within each cell type using other tools, such as identify- ing signature genes and altered pathway activities by comparing higher- and lower-hazard subpopulations (Fig. 1J). Further- more, scSurvival functions as a survival prediction model based on single-cell expression data, enabling risk assessment for new patients (Fig. 1K).Benchmarking scSurvival for Hazard Score Prediction at Single-Cell ResolutionTo systematically evaluate the effectiveness of scSurvival, we generated a series of simulated datasets and assessed its performance in controlled experiments. We first used Splatter ### Page 4 RESEARCH ARTICLEAACRJournals.orgRen et al.934 | CANCER DISCOVERY MAY 2026(23) to simulate a single-cell dataset comprising three groups (good.survival, bad.survival, and other) as the ground truth (Fig. 2A). Using these cells as risk-driving populations, we sampled single-cell expression data for 100 patients by varying the proportions of good.survival and bad.survival cells and simulated the corresponding survival times based on theseCellsGenes Attention layer Status Event Censored Feature extractor LN SE Genes SoftMax Scatter product Projection Multiply Sigmoid Tanh Sigmoid Patients SE blockbCellsA D F G H I J K E CBHazard scorer Multiply Add Linear Rare-gene weighting Multi-head attention Sum pooling Survival timex~ZIG (ρ,μ,σ) Generative model Patient ID Concat Hazard scorer Linear Linear Linear Tanh Linear Linear Linear Linear Dropout Relu and Predict patient-level hazard scores High riskP < 0.0001 0 0.00 0.25 0.50 Survival probability 0.75 1.00 500 1,000 Time (days) Risk group Low risk Pathway enrichment analysis Adjusted hazard score Pretrain: L VAEFinetune: loss = Lcox + λ1LVAE + λ2LEntgender, age, ... Differential expression profile Lower hazard Lower hazard Higher hazard Higher hazard GO KEGG GSEA Differential expression genes 10 2 3 4 5 Enrichment score Cells –15 –10 1bz2 3 4 5 –5 0 5 10 –10 –5 0 UMAP 1 UMAP 2 Pathway terms 5 10 10 5 0 –5 –10 AMICRμ σX cov X covFigure 1. Framework of scSurvival. A, Generation of single-cell cohort data: each patient undergoes scRNA-seq profiling, with corresponding clinical follow-up records. B, Single-cell transcriptomic expression matrix for each patient. C, Survival information for each patient, including event status and event time. D–H, Model architecture of scSurvival. D, A single-cell feature extraction module based on a ZIG-VAE, which takes log-normalized expression matrices as input and extracts batch-integrated single-cell features. LN and SE blocks are incorporated to adaptively normalize gene feature scales.e, A multi-head AMICR module that aggregates single-cell features to the patient level and computes relative hazard scores. f, The detailed architectures of the attention layer, the SE block in the VAE encoder, and the hazard scorer. g, The hazard scorer supports the optional inclusion of additional covariates (e.g., gender and age) if available. H, The overall loss function of scSurvival and the two-stage training strategy: the first stage involves pretraining only the VAE module to extract stable features and the second stage jointly optimizes the entire model to perform multiple-instance Cox regression and risk cell selection. i, Cell hazard scores estimated by scSurvival after adjustment using attention weights. J, Differential expression and pathway enrichment analyses of higher-hazard vs. lower-hazard cells identified by scSurvival. K, scSurvival can predict hazard scores for new patients in independent datasets to stratify patients into high- and low-risk groups. ### Page 5 RESEARCH ARTICLESingle-Cell Survival AnalysisMAY 2026 CANCER DISCOVERY | 935Bad.survival Bad.survival Status Status Group Young Old Event Censored Event Censored Patient ID 10 0 UMAP 2 Time (years)Time (years) –10 0.0 0 2 4 6 8 Bulk 1 Bulk 6 Bulk 11 Bulk 16 Bulk 21 Bulk 26 Bulk 31 Bulk 36 Bulk 41 Bulk 46 Bulk 51 Bulk 56 Bulk 61 Bulk 66 Bulk 71 Bulk 76 Bulk 81 Bulk 86 Bulk 91 Bulk 96 Patient ID Bulk 1 Bulk 6 Bulk 11 Bulk 16 Bulk 21 Bulk 26 Bulk 31 Bulk 36 Bulk 41 Bulk 46Bulk 51 Bulk 56 Bulk 61 Bulk 66 Bulk 71 Bulk 76 Bulk 81 Bulk 86 Bulk 91 Bulk 96 2.5 5.0 7.5 –10 –5 0 5 UMAP 1 10 10 0 UMAP 2UMAP 2 –10 –20 –10 0 10 UMAP 2 –10 –5 0 5 10 –10 –10 0 10 0 UMAP 1 UMAP 1 –10 –5 0 105 UMAP 1 10 20 Good.survival Batch 1 Batch 2 Good.survival Risk-cell identification (F1-score) Survival prediction (C-index, five-fold CV) Sim.0 0.968 0.331 0.317 0.262 0.941 0.909 0.977 0.951 0.953 0.904 0.971 0.865 0.814 0.7640.338 0.361 0.893 0.311 0.915 0.831 0.607 0.565 0.5090.687 0.858 0.845 0.602 0.682 0.945 0.806 0.848 0.595 0.843 0.507 0.799 0.604 0.857 0.7070.833 0.773 0.694 0.559 0.489 0.515 0.781 0.528 0.788 0.686 0.698 0.530 0.730 0.639 0.605 0.836 0.700 0.819 Sim.1 Sim.2 Sim.3 Sim.4 Sim.5 Sim.6 Overall Sim.0 Sim.1 Sim.2 Sim.3 Sim.4 Sim.5 Sim.6 Overall Good.survival Good.survival Bad.survival Bad.survival Bad.survival Other Other Other Other Other Good.survival 10 5 0 UMAP 2 –5 –10 10 5 0 UMAP 2 –5 –10 10 Status Survival time : 2 months Censor time : 10 months Survival time : 10 months Survival time : 40 months Event Censored 5 0 UMAP 2 –5 –10 –10 –5 0 UMAP 1 5 10 –10 –5 0 UMAP 1Benchmark Robustness Simulation Sim.1 Sim.3 Sim.5 Sim.6 Sim.4 Sim.2 Method Method scSurvival Cell type Coxph Pseudo-bulk Coxph scSurvival (VAE) scSurvival (PCA) scSurvival (NMF) scSurvival (HVG2000)Sim.1 Risk cell percentage Sim.3 Patient outcome distribution Sim.5 Patient group heterogeneity Sim.6 Risk subcluster heterogeneity Sim.4 Cell composition heterogeneity Sim.2 Batch imbalancehazard_adj A C F G J K H I D E B5 10 15 10 5 0 –5 –10 10 5 0 UMAP 2 –5 –10 –10–15 –5 0 UMAP 1hazard_adj(from Sim.6) 5 10 10 5 0 –5 –10 –10 –5 0 UMAP 1 5 10 10 5 0 UMAP 2 –5 –10 –10 –5 0 UMAP 1 5 10 90:10 40:60 10:90 10 5 0 UMAP 2 –5 –10 –10 –5 0 UMAP 1 5 10 10 5 0 UMAP 2Actual–5 –10 Other Higher LowerDetectedInattentive Good.survival Bad.survival 1,507 8 20 1,497 9 19 6,929 4 7 Confusion matrix Precision: 0.987 Recall: 0.990 F1-score: 0.989 Proportion 0.75 0.50 0.25ActualOther Higher LowerDetectedInattentive Good.survival Bad.survival 408 0 73 462 78 0 8,939 14 26 Confusion matrix Precision: 0.900 Recall: 0.966 F1-score: 0.931 Proportion 0.75 0.50 0.25 0.00 –10 –5 0 UMAP 1 5 10 10 5 0 UMAP 2 –5 –10 –10 –5 0 UMAP 1 5 10 10 5 0 UMAP 2 –5 –10 0.00 0 25 50 Time (months) 75 100 0.25 P < 0.0001 C-index = 0.942 ± 0.021 0.50 0.75 1.00 Survival probability –10 –5 0 UMAP 1 KM curve of risk groups (fivefold CV) Risk group High risk Low risk 5 10 Survival time : 90 months Survival time : 99 monthsFigure 2. Performance evaluation of scSurvival on simulated datasets. A, UMAP showing simulated scRNA-seq data with survival-driving cell subpop- ulations as the ground truth, including three groups: good.survival, bad.survival, and other background cells. B, UMAP plots of simulated scRNA-seq data from patients with varying survival times and differing proportions of cell subpopulations. Patients with longer survival times had higher proportions of good.survival cells and lower proportions of bad.survival cells. C, UMAP showing cell hazard scores estimated by scSurvival. D, Confusion matrix and performance evaluation for ground-truth cell subpopulation identification based on adjusted hazard scores predicted by scSurvival. e, Fivefold CV for patient-level survival prediction: merged relative risk percentiles across folds were used to evaluate risk stratification, with significance assessed by log-rank tests; the mean C-index and its SD across folds were also reported. f, Overview of six simulation scenarios designed to evaluate the performance and robustness of scSurvival under different conditions. g, UMAP or survival annotation illustrations for each simulation setting (Sim.1–Sim.6), highlighting the distinct designs and features used to introduce biological or technical variability. H, UMAP showing cell-level hazard scores (hazard_adj) estimated by scSurvival in Sim.6 (risk subcluster heterogeneity), in which true risk cells only form a small subcluster within a larger cluster. i, Confusion matrix showing the identification performance of scSurvival in Sim.6. J, Comparison of F1-score (mean of 10 repeats) for ground-truth risk cell identification across different simulation scenarios and feature extraction strategies, including scSurvival with the default VAE and alternative representations based on PCA, HVG2000, and NMF. Sim.0 refers to the baseline simulation setting shown in A and B, used as the basic scenario without added complexity. K, Comparison of patient-level survival prediction performance (mean C-index from five-fold CV) between scSurvival and two baseline models: Cell type Coxph and pseudobulk Coxph across all simulation settings. Sim.0 refers to the baseline setting illustrated in A and B. KM, Kaplan–Meier. ### Page 6 RESEARCH ARTICLEAACRJournals.orgRen et al.936 | CANCER DISCOVERY MAY 2026proportions, with random censoring applied to a subset of patients to mimic event truncation (Fig. 2B; “Methods”). We referred to this basic simulation as Sim.0. The cell attention scores and hazard scores inferred by scSurvival (Supplemen- tary Fig. S1A; Fig. 2C) accurately delineated the ground-truth risk-driving cell subpopulations. In one experiment, cell grouping based on attention and hazard scores (Supplemen- tary Fig. S1B; “Methods”) identified 99% of the ground-truth driver cells, with an overall accuracy of 98.7% and an F1-score of 0.989 (Fig. 2D). With regard to patient-level risk prediction, fivefold cross-validation (CV) showed that the hazard scores predicted by scSurvival significantly stratified patient survival outcomes (log-rank test, P value < 0.0001), achieving a mean concordance index (C-index) of 0.942 (SD, 0.021) across folds (Fig. 2E; Supplementary Fig. S1C). To assess whether scSurvival can accurately handle batch effects, we used Splatter to simulate ground-truth single-cell datasets from two distinct batches (Supplementary Fig. S1D; “Methods”). scSurvival successfully produced cell-level risk assessments consistent with the ground truth (Supplemen- tary Fig. S1E), demonstrating the capability of its embedded feature extraction module to correct batch effects (Fig. 1D). Uniform Manifold Approximation and Projection (UMAP) visualization based on the latent embeddings learned by scSurvival confirmed effective batch effect removal and accu- rate estimation of risk tendencies (Supplementary Fig. S1F and S1G). In this batch effect simulation example, scSurvival achieved an accuracy of 97.9%, a recall of 98.2%, and an F1-score of 0.980 for ground-truth risk-driving cell iden- tification (Supplementary Fig. S1H). Repeated evaluations across 20 runs further confirmed the model’s robustness in both batch-free and batch-affected settings, with high average F1-scores (0.964 vs. 0.923; Supplementary Fig. S1I). The performance of a survival prediction model is always affected by the sample size (cohort size) and the strength of association between covariates and survival time (risk effect size). To evaluate this, we simulated datasets with varying cohort sizes and risk effect sizes (“Methods”). The results showed that when the cohort size reached 10 patients and the effect size was moderate or stronger, the model consistently identified the risk-driving cell subpopulations, with F1-scores significantly higher than those observed under weaker effect sizes or smaller sample sizes (Supplementary Fig. S1J). In addition, when the cohort size was small, increasing the number of cells per patient also improved identification per- formance (Supplementary Fig. S1K). Building upon Sim.0, we further designed six more complex simulation scenarios (Sim.1–Sim.6, Fig. 2F), systematically introducing various biological and technical confounding factors at different levels of perturbation to evaluate the robustness and performance of scSurvival across diverse data conditions (Fig. 2F and G; “Methods”). The six scenarios included variation in the percentage of risk cells (Sim.1), interbatch imbalance in risk cell distribution (Sim.2), differ- ences in patient outcome distributions (varying effect size and censoring rate, Sim.3), cell composition heterogeneity (Sim.4), patient group heterogeneity in which subgroups respond differently to risk effects (Sim.5), and subcluster-level risk in which the risk-driving cells are embedded within larger cell clusters (Sim.6). Across most settings, scSurvival consistently identified the true risk-driving cell subpopulations, with aver- age F1-scores typically above 0.8 (Supplementary Fig. S2A–S2F). In certain challenging scenarios, such as when the risk cell proportion was as low as 1% (Sim.1), performance declined but rapidly recovered as signal strength increased, demon- strating strong trend stability and generalizability. Notably, in the most challenging setting (Sim.6), in which the true risk cells were nested within a large cluster and lacked clear bound- aries, scSurvival still accurately recovered the target population (Supplementary Fig. S2F; Fig. 2H), achieving a precision of 90%, recall of 96.6%, and F1-score of 0.931 in a representative example (Fig. 2I), highlighting its ability to perform fine- grained analysis without relying on predefined clustering. In patient-level survival prediction tasks, scSurvival also main- tained consistent and stable performance across all simulation conditions (Supplementary Fig. S2G–S2L), with C-index values remaining high across fivefold CV (mostly >0.8). Finally, we constructed benchmark comparisons on repre- sentative examples from each simulation setting (“Methods”) to systematically evaluate both the feature extraction module and the overall algorithmic performance of scSurvival. For feature representation, we compared our VAE-based joint learning framework with three alternative strategies: princi- pal component analysis (PCA), highly variable gene selection (HVG2000), and nonnegative matrix factorization (NMF). Joint learning of VAE consistently achieved the highest F1-scores across all simulation scenarios (Fig. 2J), highlighting its advan- tage in capturing complex expression structures and risk signals. At the algorithmic level, we further compared scSurvival with two commonly used baseline approaches: Cox regres- sion based on pseudo-bulk expression profiles (pseudobulk Coxph) and Cox regression based on cell type proportions (cell type Coxph). Because both baseline methods lack cell- level resolution, comparisons were restricted to patient-level survival prediction. Across most scenarios, scSurvival out- performed the baselines, achieving higher C-index values in fivefold CV (Fig. 2K). Together, these comprehensive simula- tion results demonstrate that scSurvival offers superior capa- bilities in both risk subpopulation identification and survival prediction, while fully preserving single-cell resolution.scSurvival is fast and ScalableCohort-level cancer studies often involve sequencing large numbers of patients, generating datasets with millions of cells. Performing survival analysis at this scale poses significant computational challenges for algorithm design. Benefiting from the supervised learning architecture of scSurvival, we can fully leverage the computational power of high-performance graphics processing units to process datasets of this magni- tude. To evaluate its computational performance, we simulated a series of datasets and measured runtime, memory usage, and GPU memory consumption. In these tests, early stopping was disabled, and the fine-tuning phase of each run was set to com- plete all 500 maximum epochs (“Methods”), suggesting that in practical applications, scSurvival could handle similarly sized datasets in even shorter times. We simulated two scenarios to evaluate the scalability of scSurvival, keeping the number of genes fixed at 2,000 in all cases, consistent with the top 2,000 most variable gene ### Page 7 RESEARCH ARTICLESingle-Cell Survival AnalysisMAY 2026 CANCER DISCOVERY | 937(MVG2000) feature selection strategy commonly used in real analyses. First, we fixed the cohort size at 100 patients and gradually increased the total number of cells from 5,000 to 1 million, distributing cells evenly across patients (Fig. 3A). The runtime and memory usage of scSurvival exhibited linear growth relative to the number of cells (x-axis in log scale), completing the analysis of 1 million cells in approximately 17.5 minutes, with around 30 GB of memory and 50 GB of GPU memory usage. We also compared cases with and without batch effect adjustment, treating each patient as a distinct batch. The results showed that batch correction did not intro- duce significant additional computational burden. Next, we fixed the total number of cells at 100,000 and simulated cohort sizes ranging from 10 to 500 patients (Fig. 3B). Runtime again exhibited linear scaling with respect to the number of patients (x-axis in log scale), whereas memory usage remained nearly constant. GPU memory usage increased slightly when batch correction was included but remained stable without it. In the largest simulated cohort of 500 patients, scSurvival completed the analysis in approximately 33 minutes, with a peak memory usage of about 3 GB and a peak GPU memory usage of about 6 GB. All tests were performed on an Arm64 CPU and H100 GPU with 96 GB VRAM. As different hardware configurations can lead to substantial performance variations, we also bench- marked scSurvival on the A100 platform. The computational cost remained within an acceptable range, approximately 2.5 hours for 1 million cells across 100 patients and about 1.85 hours for 100,000 cells across 500 patients (Supplementary Fig. S3A and S3B). These results demonstrate the scalability of scSurvival for processing ultra-large single-cell datasets.scSurvival identifies Risk-Associated Cell Subpopulations in a Melanoma CohortFor real-world evaluation, we first applied scSurvival to a melanoma cohort dataset (24) treated with immunotherapy, comprising 48 samples from 32 patients (Supplementary Fig. S4A; Supplementary Table S1) with a total of 16,291 im- mune cells annotated according to the original study (Fig. 4A; ref. 24). The key cell populations, together with their attention scores and hazard scores inferred by scSurvival, are shown in Fig. 4B and Supplementary Fig. S4B and S4C. Among them, B cells and plasma cells generally exhibited lower hazard scores (Fig. 4C), a finding potentially linked to the formation of tertiary lymphoid structures (TLS; refs. 26, 27). This associa- tion was further supported by spatial transcriptomics valida- tion, in which lower-hazard regions identified by the trained scSurvival model corresponded to TLS locations on two spatial slices (Supplementary Fig. S5A and S5B; ref. 28). By contrast, monocytes and macrophages were largely retained and showed higher overall hazard scores, but after further subgrouping of these cells (“Methods”), a bimodal pattern of hazard scores was observed, characterized by a clear polarization into higher- hazard and lower-hazard cell subsets (Fig. 4C). This observa- tion was of particular interest as monocytes/macrophages are known to polarize into two major functional states—M1-like (antitumor) and M2-like (protumor) phenotypes (29). How- ever, recent single-cell studies suggest that M1 and M2 phe- notypes define the possible extremes of polarized cells in vitrobut may not capture the full spectrum of tumor-associated macrophages (TAM) in vivo (25). We therefore sought to evaluate whether scSurvival could reliably distinguish pro- tumor and antitumor monocyte/macrophage populations by stratifying them into higher- and lower-hazard cell subsets. Reclustering analysis of the monocyte/macrophage com- partment revealed that higher- and lower-hazard cells tend to cluster separately (Fig. 4D), suggesting transcriptional and potentially lineage similarity within each risk group. Interest- ingly, we also observed that some subclusters, such as cluster 1, contained a mixture of both higher- and lower-hazard cells (Fig. 4D), indicating intracluster heterogeneity in risk states. Notably, the distribution of M1 and M2 gene signature scores showed strong concordance with our predicted risk states, fur- ther supporting the biological relevance of the subpopulations identified by scSurvival (Fig. 4D; Supplementary Fig. S5C). To further explore transcriptional differences between these cell populations, we performed differential expression analysis comparing higher-hazard and lower-hazard cells. As a result, 169 genes were significantly upregulated and 397 genes were significantly downregulated in higher-hazard cells relative to lower-hazard cells (Fig. 4E; Supplementary Table S2). Notably, several key M2-associated markers (30), including SPP1 andMSR1, were among the upregulated genes in higher-hazard cells. Conversely, several canonical M1-associated markers (e.g., IDO1 and CXCL9) and MHC class II genes (e.g., HLA-DOB) were significantly downregulated in these cells. Of particular interest, SPP1 and CXCL9 ranked among the top upregulated and downregulated differentially expressed genes (DEG), respectively (Fig. 4E; Supplementary Table S2). A recent study (25) has identified CXCL9 and SPP1 as key markers of macrophage polarization. Specifically, the CXCL9:SPP1 ratio was shown to reflect antitumor immune abundance and responsiveness to immunotherapy (25). Consistent with this, we observed that lower-hazard monocytes/macrophages pre- dominantly exhibited an SPP1− CXCL9+ phenotype, whereas higher-hazard cells were characterized by an SPP1+ CXCL9−profile (Fig. 4F; Supplementary Fig. S5D). Furthermore, Gene Ontology (GO) enrichment analysis revealed that DEGs upregulated in higher-hazard cells were enriched in pathways associated with neutrophil migration and chemotaxis and several lipid metabolism signals such as lipid localization and transport (Fig. 4G). In contrast, DEGs upregulated in lower- hazard cells were enriched in pathways related to immune activation, including lymphocyte differentiation and activation, type II interferon response, and MHC class II antigen pre- sentation (Fig. 4G). These enriched pathways are consistent with known features of TAMs, in which lipid metabolism and neutrophil-recruiting chemokines are commonly linked to pro- tumor phenotypes, whereas antigen presentation and interferon responses are hallmarks of antitumor macrophage activation (bioRxiv 2024.12.18.629070; refs. 30–32). Importantly, these biological distinctions between protumor and antitumor myeloid states were also recapitulated in an independent melanoma dataset treated with immunotherapy [GSE221553; ref. (33); Supplementary Fig. S6A–S6H], further supporting the robustness of scSurvival in integrating survival informa- tion to distinguish functionally distinct monocyte/macrophage populations. ### Page 8 RESEARCH ARTICLEAACRJournals.orgRen et al.938 | CANCER DISCOVERY MAY 2026Besides, based on morphologic assessments of immuno- therapy response, the original study classified the 48 samples from 32 patients into 17 responders and 31 nonresponders (24). As scSurvival enables the assessment of cellular risk het- erogeneity within each patient, it revealed intrapatient variability in risk states. Overall, the key cell subpopulations defined by scSurvival (excluding inattentive cells) and their associated hazard scores showed strong concordance with immunother- apy response: key cells in responder samples had significantly lower hazard scores than those in nonresponder samples (Fig. 4H; Supplementary Fig. S7A and S7B). The proportions of lower- and higher-hazard cells within each sample further demonstrated that responder samples contained significantly more lower-hazard cells, consistent with the findings of the original study (Fig. 4I; ref. 24). Figure 4J provides an intuitive visualization of this refined cellular risk heterogeneity. Notably, among the four patients who had both responder and nonre- sponder samples, three patients (P28, P5, and P4) exhibited cellular risk assessments consistent with their response status (Fig. 4J; Supplementary Fig. S7C). By preserving risk hetero- geneity at the cellular level, scSurvival also enables the recon- struction of sample-level summaries that integrate immune composition, survival information, and model attention. Sup- plementary Figure S7D summarizes these multidimensional outputs across all samples, providing an in-depth overview that highlights key immune features underlying patient-level risk. At the patient level, beyond assessing immunotherapy response, we also conducted leave-one-out CV to evalu- ate the predictive performance of the scSurvival-derived risk model based on single-cell expression data. The results demonstrated that the predicted patient hazard scores sig- nificantly stratified patient survival outcomes (P = 0.033) and accurately ranked risks (C-index = 0.812; Fig. 4K), indicating that scSurvival established a reliable survival prediction model at the patient level. Furthermore, univariate Cox regression confirmed that scSurvival-derived hazard scores were signifi- cantly associated with OS (P = 0.022), with a hazard ratio comparable with cytotoxicity and more significant than other conventional immune features, including the average expres- sion ratio of CXCL9 and SPP1 (Supplementary Fig. S8A). The achieved C-index of 0.812 also outperformed benchmark models based on cell type fractions, pseudobulk expres- sion, or the CXCL9/SPP1 expression ratio (Supplementary Fig. S8B).scSurvival identifies Prognostic T-Cell States andenables Survival Prediction in MelanomaT cells play a crucial role in tumor immunotherapy. To better highlight the role of T cells, we isolated T cells (10,685 cells) from the previous melanoma dataset (Fig. 5A) and analyzed them with scSurvival. To more rigorously assess the predictive power of scSurvival, we further evaluated a17.5 30 3,000 6,000 5,000 4,000 3,000 2,000 2,750 2,500 2,250 2,000 10 30 100 Number of patients 300 10 30 100 Number of patients 300 10 30 100 Number of patients 300 20 10AB30,000 50,000 GPU memory (MiB) 40,000 30,000 20,000 10,000 0 Memory (MiB) 20,000 10,000 0 With batch Without batch With batch Without batch 1e+04 1e+05 Number of cells 1e+06 1e+04 1e+05 Number of cells 1e+06 1e+04 1e+05 Number of cells 1e+06 Elapsed time (minutes)Elapsed time (minutes) Memory (MiB) GPU memory (MiB) 15.0 12.5 10.0 7.5Figure 3. Runtime and resource usage of scSurvival on the H100 high-performance GPU computing platform. A, Runtime and memory usage of scSur- vival for simulated cohorts of 100 patients across varying total numbers of cells, with or without simulated batch effects. B, Runtime and memory usage of scSurvival for simulated cohorts with a fixed total of 100,000 cells across varying numbers of patients, with or without simulated batch effects. ### Page 9 RESEARCH ARTICLESingle-Cell Survival AnalysisMAY 2026 CANCER DISCOVERY | 939Cell type Cell typeCytotoxicity Exhausted/HS CD8+T cells Exhausted CD8+T cells Lymphocytes exhausted/cell cycle Monocytes/macrophages Dendritic cells –15 –10 0 –10 –5 UMAP 1 UMAP 2 –10 2.5 0.0 –2.5 –5.0 –0.6 –0.3 0.0 0.3 10 5 0 –5 –10 1050–5–10 0 UMAP 2 Avg log 2 FC hazard_adj hazard_adj Expression level 0 5 –15 –10 –5 UMAP 1 0 5 –10 –5 hazard_adj hazard_adj Cell risk groupRisk enrichedResponderResponderNonresponder 10 1.000.75 0.50 Survival probability 0.25 P = 0.033 C-index = 0.812 0.00 0 250 500 Time (days) Number at risk Risk group High risk 16 9 6 3 1 0 16 15 9 1 1 1Low risk 750 1,000 1,250 0 250 500 Time (days) 750 1,000 1,250 5 0 –5 –10 P22 P12 P15P28P11 P2P5P16 P31 P3 P27 P14P1 P4P6 P7 P8 P10 P13 P17 P18 P19 P20 P21 P23 P24 P25 P26 P29 P30 P33 P35NonresponderHigher Higher- vs. lower-hazard cells Wilcoxon, P = 1.5e–10 Wilcoxon, P < 2.2e–16 High Low Enrichment analysis of DEGs (higher vs. lower) Lower Higher Inattentive Lowerhazard_adjPercentage 0 0% 25% 50% 75% 100% 20 10 0 −Log 10 (P value) 10 5 10 Lymphocytes Memory T cells Regulatory T cells B cells Plasma cellsA CG KJH I B E F DB–cells (1,455) Cytotoxicity (2,165) Dendritic cells (290) Exhausted CD8+ T cells (2,222) Exhausted/HS CD8+ T cells (1,656) Lymphocytes (2,165) Lymphocytes exhausted/cell cycle (1,129) Memory T cells (1,773) Monocytes/macrophages (1,391) Plasma cells (305) Regulatory T cells (1,740) Proportion of risk cells in each cell type Monocytes/macrophages Monocytes/macrophages M1-like score UMAP 1 UMAP 2 M2-like score High Low Cell risk group Lower Inattentive Higher 0 1 2 3 4 5 B cells Cytotoxicity Dendritic cells Exhausted CD8 +T cells Exhausted/HS CD8 +T cells Lymphocytes Lymphocytes exhausted/cell cycle Memory T cells Monocytes/macrophages Plasma cells Regulatory T cellsSPP1SPP1 CXCL9/SPP1 CXCL9 Wilcoxon, P < 2.2e–16 CXCL9/SPP1 6 4 Lower Higher 2 0 Ratio of CXCL9/SPP1 expression 6 4 2 0 Expression level 3 10P < 2.22e–16 5 0 –5 –10 100 Lower Higher Cell risk group Patient (high hazard –> low hazard) Lower Higher ** ns Risk group in samples (%) 75 50 25 0 Responder Nonresponder Lower Higher Lower Higher 2 1 0CTSD MSR1 MMP19 CTSB FN1 HLA–DOA CD40 CTSA CXCL10 CXCL9 CXCL11 HLA–DOB Up DownIDO1Difference (pct.1 – pct.2) Phagocytosis Neutrophil migration Granulocyte migration Neutrophil chemotaxis Apoptotic cell clearance Regulation of lipid localization Cholesterol transport Cholesterol storage Positive regulation of lipid transport Regulation of angiogenesis Upregulated in lower-hazard cells Upregulated in higher-hazard cells Dendritic cell migration B-cell activation Dendritic cell differentiation Response to type II interferon Regulation of immune effector process Lymphocyte differentiation Leukocyte cell–cell adhesion Regulation of T-cell activation MHC protein complex assembly Antigen processing and presentation KM curves by risk group (LOO-CV) Risk group High risk Low riskP adjust = 0.05Figure 4. scSurvival analysis of a melanoma scRNA-seq cohort. A, UMAP visualization showing diverse cell types as defined in the original article (24, 25). Cell numbers are indicated in parentheses. B, UMAP showing adjusted hazard scores estimated by scSurvival for each cell. C, Distribution of adjusted hazard scores and proportions of higher-hazard, lower-hazard, and inattentive cells within each cell type. Red and blue arrows indicate cell types enriched for higher-hazard and lower-hazard populations, respectively. Arrow direction is determined by the position of the maximum density estimated by Kernel Density Estimation (KDE) relative to zero (threshold set to 1; only directions with an absolute maximum hazard greater than 1 are shown). D, UMAP visualization showing the distribution of six monocyte/macrophage clusters, higher-hazard and lower-hazard subsets, and enrichment score for M1-like and M2-like signatures. e, Scatter plot of DEGs between higher-hazard and lower-hazard monocytes/macrophages. Red dots represent genes upregulated in higher-hazard cells, blue dots represent genes upregulated in lower-hazard cells, and gray dots denote nonsignificant genes. f, Violin plots showing the expression of SPP1, CXCL9, and their expression ratio (CXCL9/SPP1) in higher-hazard versus lower-hazard monocytes/macrophages. UMAP plot with cells colored by the CXCL9/SPP1 expression ratio. g, GO Biological Process enrichment analysis of DEGs between higher-hazard and lower-hazard monocytes/macrophages. H, Violin plot comparing adjusted hazard scores of cells between responders and nonresponders in the melanoma cohort. P value calculated using a one-sided Wilcoxon test. i, Boxplots comparing the proportions of higher-hazard versus lower-hazard groups within responders and nonresponders. P value calculated using a one-sided Wilcoxon test. **, P < 0.01; ns, not significant. J, Boxplots showing adjusted hazard scores of each patient, ranked from higher to lower patient-level hazard. K, Leave-one-out (LOO) CV at the patient level, evaluating the stratification ability of estimated patient-level hazard scores and reporting the overall C-index. P value calculated using the log-rank test. KM, Kaplan–Meier.T cell–based survival model across multiple independent scRNA-seq melanoma datasets containing T-cell expression profiles (Supplementary Table S1). In the T cell–only analysis, scSurvival identified a large number of key cells, with their attention and hazard scores shown in Fig. 5B and Supplementary Fig. S9A. Overall, cytotoxic ### Page 10 RESEARCH ARTICLEAACRJournals.orgRen et al.940 | CANCER DISCOVERY MAY 2026AD F I J K L M N G H E B CMelanoma T cells hazard_adj Cell risk group5 0 UMAP 2 –5 5 20 100 100 80 60 40 20 75 50 25 10 0 –10 NR (34) 0 20 40 R (57) NR (41) R (29) 10 0 –10 0 –5 6 7.5 10 5 0 –5 –10 5.0 2.5 0.0 –2.5 –5.0 3 0 –3 –10 –5 0 UMAP 1 5 10 –15 –10 –5 UMAP 1 0 5 –10 5 10 0 hazard_adj Avg log 2 FC Signature scoreOS (%) Survival probability PFS (%) –5 –10 P22 P12 –0.2 –2 –1 0 1 2 0.0 0.2 0.4 P15 P11 P2 P5 P16 P31 P3 P14 P28 P26 P29 P7 P35 P24 P27 P1 P4 P6 P8 P10 P13 P17 P18 P19 P20 P21 P23 P25 P30 P33 5 10 5 0 –5 –10 0 UMAP 2Cell type–5 –10 –8 –4 0 UMAP 1 Patient (high hazard –> low hazard)GZMH CCL4T cytotoxicity P.adjust Count 20 40 60 0.005 Wilcoxon, P = 0.00023 CancerCell.Gide.2019 CancerCell.Gide.2019 Log-rank Log-rankP = 5.16e–04 P = 1.03e–05 Low (46) High (45) Low (46) High (45) Time (months) 0 20 40 Time (months) Wilcoxon, P = 5e–06 0.010 0.015 Activation:effector function Cytotoxicity Stress Antigen processing and presentation TCR signaling Chemokine signaling pathway TNFD signaling via NF-NB UV response DN Naïve NESHSPA1A HSPA1B HSPH1 DNAJB1 TNFRSF25 SELL IL7R TCF7 CCR7 LEF1Up DownLAG3 CCR5 GZMA NKG7 HAVCR2Difference (pct.1 – pct.2)IFNG PRF1Higher vs. lower hazard (T cells) Higher vs. lower hazardK411 (NR) GSE134388 (R)hazard_adj hazard_adjK468 (NR)hazard_adj Multivariable Cox model Sex Age Patient hazardsPRJNA679099 PRJEB23709 Allen and Hugo Patient survival predictionhazard_adj Risk group High risk Low risk Nonresponders Responder 4 5 0 UMAP 2 UMAP 2 UMAP 2 10 5 4 1.00 0.75 0.50 0.25 P = 0.008 C-index = 0.757 0.00 0 Risk group 250 500 Time (days) 750 1,000 0 –4 –8 13 patients (76,112 cells) 5 0 –5 –10 5 0 –5 –10 5 0 –5 –10 5 10 0 –5 –10 0 –5 –10 –5P < 2.22e–16 GSE134388 (R) K411 (NR) K468 (NR)P < 2.22e–16 0 UMAP 1 5 –5 0 UMAP 1 Reference Reference F(N = 4)M(N = 9) 0.25(0.015–4.4) 0.344 0.3150.04*3.63(0.294–45.0)1.39(1.015–1.9)<60(N = 7)=60(N = 6) (N = 13)5 0.01 0.050.1 1 5 10 500.5 HR (95% CI) UMAP 2 hazard_adj UMAP 2 –5 –10 –8 –10 –5 0 hazard_adj hazard_adj 5 10 –4 0 UMAP 1 4 –8 –4 0 UMAP 1 4 Cytotoxicity (2,165) Higher Risk enriched Lower Higher Lower Inattentive Cytotoxicity Exhausted CD8+T cells Exhausted/HS CD8+ T cells Lymphocytes exhausted/cell cycle Memory T cells Regulatory T cells Exhausted CD8 +T cells (2,222) Exhausted/HS CD8 +T cells (1,656) Lymphocytes exhausted/cell cycle (1,129) Memory T cells (1,773) Regulatory T cells (1,740) –10 –5 0 5 10 Exhausted CD8T Number at risk High risk 3 2 1 10 0 250 500 Time (days) 750 1,000 9 7 1 7 0 6Low riskFigure 5. scSurvival analysis and prediction using T cells from a melanoma cohort. A, UMAP visualization generated using MVG2000, showing diverse T-cell subtypes. Cell numbers are indicated in parentheses. B, UMAP showing T-cell hazard scores estimated by scSurvival. C, Distribution of cell hazard scores across T-cell subtypes, with the overall risk trend indicated by red and blue arrows. D, Boxplot showing predicted hazard scores for each patient based on T cells, ranked by patient-level hazard scores. e, UMAP of all T cells colored by risk grouping. f, Scatter plot of DEGs between higher-hazard and lower-hazard T-cell groups. Red dots represent genes upregulated in higher-hazard cells, blue dots represent genes upregulated in lower-hazard cells, and gray dots denote nonsignificant genes. g, GSEA results comparing higher-hazard and lower-hazard T cells. H, A gene signature derived from DEGs between high-hazard versus low-hazard groups, evaluated in independent bulk RNA-seq melanoma datasets. Boxplots show signature scores in responder versus nonresponder groups across two independent datasets, with P values calculated using the Wilcoxon test. Kaplan– Meier survival curves compare patients with high versus low signature scores (median cutoff) in two independent datasets, with P values calculated using two-tailed log-rank tests. i, UMAP plots of predicted cell hazard scores in two independent nonresponder single-cell datasets using the trained scSurvival T-cell model. J, UMAP showing predicted cell hazard scores in an independent responder scRNA-seq sample using the trained scSurvival T-cell model. K, Violin plot comparing predicted hazard scores of cells from three new scRNA-seq samples. P value calculated using a one-sided Wilcoxon test. L, UMAP of predicted cell hazard scores of T cells from an independent melanoma cohort (n = 13 patients) using the trained scSurvival T-cell model. M, Risk stratification and C-index evaluation based on predicted patient-level hazard scores in this independent cohort. P value calcu- lated using the log-rank test. n, Multivariable Cox regression analysis incorporating sex, age, and scSurvival-derived patient-level hazard scores. Data are presented as hazard ratio (HR) ± 95% confidence interval (CI) and P value was calculated using a Wald test. NES, normalized enrichment score; NR, nonresponder; R, responder; TCR, T-cell receptor. ### Page 11 RESEARCH ARTICLESingle-Cell Survival AnalysisMAY 2026 CANCER DISCOVERY | 941T cells, exhausted CD8+ T cells, and regulatory T cells tended to be associated with longer survival, whereas memory T cells were linked to shorter survival (Fig. 5C). Notably, a subset of cytotoxic CD8+ T cells in this cohort expressed high levels of stress- and exhaustion-related markers (e.g.,TIGIT, HSPA1A, and HSPA1B; Supplementary Fig. S9B), which may explain why these cells were associated with unfavorable outcomes in the Cox regression model (Sup- plementary Fig. S8A). Intrapatient cellular hazard scores remained largely consistent with immunotherapy response status (Fig. 5D). Furthermore, lower-hazard T cells were enriched in effective immune checkpoint blockade (ICB) responses, whereas higher-hazard T cells predominated in nonresponders (Fig. 5E; Supplementary Fig. S9C and S9D). Differential expression analysis revealed that lower-hazard T cells displayed elevated expression of genes associated with T-cell longevity and stem-like memory potential (34), such as TCF7, CCR7, and IL7R (Fig. 5F; Supplementary Fig. S9E; Supplementary Table S3). In contrast, higher-hazard T cells upregulated multiple immune-inhibitory receptors, including HAVCR2 (TIM-3) and LAG3. Intriguingly, higher- hazard T cells also expressed a broad stress response program, with marked induction of heat shock genes such as HSPA1Aand HSPA1B (Fig. 5F; Supplementary Fig. S9E). A recent pan-cancer study has described a distinct intratumoral T-cell subset, termed TSTR cells, characterized by pervasive heat shock gene expression and linked to poor clinical outcomes in immunotherapy (35). In line with these findings, gene set enrichment analysis (GSEA) revealed strong enrichment of cellular stress pathways, T-cell exhaustion signatures, and cytotoxic dysfunction programs in higher-hazard T cells (Fig. 5G; Supplementary Fig. S9F), reinforcing the notion that stress-adapted, dysfunctional T-cell states underlie poor therapeutic responses (35). Furthermore, the DEG-derived signature score was significantly higher in ICB nonre- sponders than in responders across multiple independent ICB datasets, and it effectively stratified both OS and PFS in the PRJEB23709 bulk RNA-seq melanoma immunotherapy cohort (Fig. 5H; ref. 36). The trained scSurvival model can also be applied to predict patient survival outcomes in independent melanoma single- cell datasets. We first applied the model to three independent patients with melanoma treated with ICB, including two nonresponder samples (Fig. 5I; ref. 37) and one responder sample (Fig. 5J; ref. 38). The results showed that T cells in the responder sample had significantly lower hazard scores com- pared with those in the two nonresponder samples (Fig. 5K). We further validated the model on an independent melanoma cohort with survival outcomes (PRJNA679099; ref. 39), which includes 13 pretreatment samples comprising 76,112 cells (Supplementary Fig. S9G and S9H; Fig. 5L). The predicted patient hazard scores achieved a C-index of 0.757, and risk stratification (using optimal grouping) reached statistical significance (P = 0.008; Fig. 5M; Supplementary Fig. S9I). Univariate and multivariate Cox regression analyses based on the continuous patient hazard scores further confirmed that the predicted risks were significantly associated with survival outcomes (P = 0.04; Supplementary Fig. S9J; Fig. 5N). In addition, given that the training dataset included both pretreatment and posttreatment samples, we performed stratified analyses. The results showed that the features and signals identified by scSurvival were conserved across both pretreatment and posttreatment subsets (Supplementary Fig. S10A and S10B).scSurvival Reveals Survival-Associated Cell Subpopulations in a Liver Cancer CohortWe next applied scSurvival to a large liver cancer scRNA-seq atlas comprising 189 samples from 124 patients, including both tumor and TME cells with annotated cell types (40). Among these, survival information was available for 121 patients (Supplementary Fig. S11A; Supplementary Table S1). UMAP projection revealed diverse cell populations, including tumor/epithelial cells, multiple immune cell subtypes, and stromal components for a total of 1,092,172 cells (Fig. 6A). scSurvival computed cell-level hazard scores across all cells (Fig. 6B), which were further stratified into higher-hazard, lower-hazard, and inattentive groups via an attention-guided risk grouping strategy (Fig. 6C). Distinct cell types exhibited different distributions of hazard scores (Fig. 6D; Supplemen- tary Fig. S11B). Consistent with previous results from the melanoma immunotherapy cohort analysis, scSurvival iden- tified higher-hazard macrophages characterized by elevatedSPP1 and reduced CXCL9 expression (Supplementary Fig. S11C). This expression pattern is also consistent with recent findings in liver cancer, in which SPP1+ TAMs promote tumor progression and stemness (41, 42). In the tumor cell compartment, approximately one third of the epithelial cells were identified as higher hazard, with the remainder as lower hazard (Fig. 6D), highlighting sub- stantial heterogeneity among tumor cells that contribute to patient prognosis. To dissect transcriptional differences, we compared gene expression profiles of higher- and lower- hazard tumor cells. As a result, 1,672 genes were upregulated and 997 genes were downregulated in higher-hazard tumor cells (Fig. 6E; Supplementary Table S4). Notably, the higher- hazard cells exhibited upregulation of genes linked to hypoxia (HIF1A), stemness (PSCA and APLP1), and epithelial– mesenchymal transition (EMT; S100A4 and SPP1; ref. 43). GO and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses confirmed activation of EMT-related pathways, including wound healing, epithelium migration, and hypoxia-inducible factor-1 signaling pathway, in high- er-hazard tumor cells (Fig. 6F; Supplementary Fig. S11D). Furthermore, GSEA revealed significant enrichment of EMT, TNFα signaling via NF-κB, hypoxia, and TGFβ signaling in the higher-hazard group, indicating a more invasive phe- notype (Fig. 6G and H). In contrast, lower-hazard tumor cells showed elevated expression of liver-enriched metabolic genes such as ARG1 and ALDOB. Correspondingly, pathway enrichment analysis highlighted activation of metabolic programs in the lower-hazard group (Fig. 6F and G; Sup- plementary Fig. S11D), suggesting that these cells exhibit heightened metabolic activity and a hepatocyte-like transcrip- tional profile. This observation aligns with previous reports that a subset of tumor cells can display a well-differentiated phenotype with hepatocyte-like features (44), which are typically associated with more favorable clinical outcomes in liver cancer. ### Page 12 RESEARCH ARTICLEAACRJournals.orgRen et al.942 | CANCER DISCOVERY MAY 20265.0 2.5 0.0 –2.5 –5.0 –7.5 0 10 20 30Cell type A C D B G H K L J E I Fhazard_adj Cell risk groupEndothelial Tumor/endothelial 10 0 UMAP 2 –10 –10 0 UMAP 1 10 10 10 5.0 2.5 FOXC2 Higher vs. lower hazard (tumor cells) Wound healing Higher vs. lower hazard (tumor cells) TNFD signaling via NF-NB Branch 1 0.08 0.16 0.12 0.08 0.04 0 10 20 Pseudotime 30 0.09 0.07 0.05 0.03 Branch 1 0 10 20 30 Pseudotime 0 10 20 30 Pseudotime EMT EMT NES = 1.78P adj = 4.01e–04 NES = 1.89P adj = 8.2e–04 Fatty acid metabolism Cell risk group Lower Component 2 Component 1 Inattentive Higher 21 21 Adjusted hazard score TGFE signaling TGFE signaling GSE14520 (n = 445) TCGA LIHC (OS) GSE14520 (OS) TCGA LIHC (PFS) GSE14520 (PFS) KM curves by risk group (fivefold CV)P < 0.0001 C-index = 0.719 r0.098 Risk group High risk Low risk Higher vs. lower 100 OS (%) PFS (%) 80 60 40 20 100 OS (%) 90 80 70 60 100 PFS (%) 80 60 40 50 100 1.00 0.75 0.50 Survival probability 0.25 0.00 0 10Time (months) Time (months) Number at risk 0 10 20 30 Risk group 60High risk Low risk 61 59 53 7 49 31 7 20 30 80 60 40 20 0 0 20 40 Time (months) 60 0 20 40 Time (months) 60 25 50 75 Time (months) 100 125 0 25 50 75 Time (months) 100 125 Log-rank scLCSS Low (187) High (186) scLCSS Low (111) High (110) scLCSS Low (111) High (110) scLCSS Low (187) High (186)P = 0.01 Log-rankP = 2.94e–03 Log-rankP = 0.009 Log-rankP = 0.005 GSE36376 (n = 433) GSE36411 (n = 84) TCGA LIHC (n = 400) GSE76427 (n = 166) Nontumor Tumor Adjacent Tumor scLCSS (VIPER) P < 2.22e–16 P < 2.22e–16 3.6e–06 4.4e–10 0.0025 0.00022 0.007 9.8e–091.8e–05 1.4e–05 0.014 4 2 0 –2 4 6 5.0 2.5 0.0 –2.5 4 2 0 –2 6 4 2 0 –2 2 0 –2 Branch 2 Branch 1 Branch 2 Branch 1 Branch 2 Branch 1 Branch 2R = 0.26, P < 2.2e–16R = –0.79, P < 2.2e–16 R = 0.18, P < 2.2e–16R = 0.44, P < 2.2e–16R = 0.51, P < 2.2e–16R = –0.81, P < 2.2e–16 R = 0.33, P < 2.2e–16R = 0.56, P < 2.2e–16 AUCell scoreAUCell score Adjusted hazard score AUCell score 0.07 0.06 0.05 0.04 Branch 2 Pseudotime 0 10 20 30P.adjust Hypoxia TGFE signaling 0.002 Apical junction 0.004 IFND response 0.006 EMT IFNJ response Count –2 –1 0 1 2 IL2 stat5 signaling 25 Kras signaling up 50 Inflammatory response 75 Mitotic spindle 100 P53 pathway 125 Glycolysis Adipogenesis Fatty acid metabolism Peroxisome Oxidative phosphorylation Coagulation Xenobiotic metabolism Bile acid metabolism NES Epidermal growth factor receptor signaling pathway ERBB signaling pathway Cell–substrate adhesion Epithelium migration Response to hypoxia Endothelial cell migration Epithelial cell proliferation Positive regulation of angiogenesis Mesenchymal cell differentiation Proteinogenic amino acid metabolic process Carboxylic acid biosynthetic processUpregulated in higher-hazard cells Upregulated in lower-hazard cellsOrganic acid biosynthetic process Amino acid metabolic process Amino acid catabolic process Alpha–amino acid metabolic process Fatty acid metabolic process Small-molecule catabolic process Carboxylic acid catabolic process Organic acid catabolic process GOBP enrichment analysis of DEGs (higher vs. lower) −Log10 (P adjust) –60 –40 –20 0 10 20 PSCA LAMB3 APLP1 S100A4 PLAUR CCND2 EGFR ITGA2 ITGA3 SPP1 VEGFA HP LAMC2 HIF1APHGDH RBP4 ALDOB CYP2E1 CFHR1 PCK1 ARG1 GLUL SLC22A1 GNMT 0.0 Avg log 2 FC –2.5 –5.0 –0.4 –0.2 0.0 Difference (pct.1 – pct.2) Down NS Up 0.2 0.4 5 0 –5 –10 0 UMAP 2 –10 –10 0 UMAP 1 10 10 0 UMAP 2 –10 –10 0 Higher 100% 75% 50% Percentage 25% 0% Proportion of risk cells in each cell type Lower Inattentive UMAP 1 10 Mesenchymal CD4 + T B CD8 + T Treg NK gdT Mast Monocytes Neutrophils Dendritic Macrophages Endothelial Tumor/epithelial Mesenchymal CD4+ T B CD8 + T Normal CL Diff Undiff Adjacent I II III Adjacent I II III–IV Treg NK gdT Mast Monocytes Neutrophils Dendritic Macrophages B (37,597) Cell risk group Higher Lower Inattentive CD4 + T (162,752) CD8 + T (215,172) Dendritic (19,033) Endothelial (84,864) gdT (14,229) Macrophages (152,999) Mast (2,784) Mesenchymal (55,975) Monocytes (19,768) Neutrophils (34,307) NK (60,626) Treg (38,189) Tumor/epithelial (193,877)Figure 6. scSurvival analysis of a liver cancer single-cell cohort with survival data. A, UMAP visualization of a liver cancer scRNA-seq dataset colored by cell types with cell numbers indicated in parentheses. B, UMAP visualization of cell hazard scores estimated by scSurvival. C, UMAP visualization of cell risk groups, including higher-hazard, lower-hazard, and inattentive populations. D, Bar plots showing the proportions of higher-hazard, lower-hazard, and inattentive cells across cell types. e, Scatter plot of DEGs between high-hazard and low-hazard tumor cells. Orange dots indicate upregulated (up) genes, blue dots indicate downregulated (down) genes, and gray dots represent nonsignificant (NS) genes. (continued on following page) ### Page 13 RESEARCH ARTICLESingle-Cell Survival AnalysisMAY 2026 CANCER DISCOVERY | 943Pseudotime analysis revealed a trajectory from lower-hazard tumor cells branching into two higher-hazard tumor states (Fig. 6I). Notably, tumor cells from branch 2 displayed elevated EMT activity and showed a stronger correlation with the adjusted hazard score (Fig. 6I), suggesting that this subpopu- lation is linked to greater invasiveness and poorer prognosis. As a downstream analysis, we extended the scSurvival- identified higher-hazard and lower-hazard tumor cell states into bulk RNA-seq datasets to further assess their biological and clinical relevance. First, we combined the 100 most upregulated and 100 most downregulated DEGs into a uni- fied single-cell liver cancer unfavorable survival signature (scLCSS; Supplementary Table S4). The scLCSS signature scores calculated on bulk RNA-seq datasets robustly distin- guished tumor samples from nontumor or adjacent tissues, with scores increasing progressively with advancing tumor pathologic stage (Fig. 6J). More importantly, across four inde- pendent liver cancer cohorts with clinical follow-up, patients classified as scLCSS-high exhibited shorter OS and PFS than their scLCSS-low counterparts (Fig. 6K; Supplementary Fig. S11E). These findings highlight the ability of scSurvival to identify survival-associated cell subpopulations that are clinically relevant to prognosis. Finally, we evaluated the predictive performance of scSurvival at the patient level. We conducted five-fold CV using patients from the liver cancer cohort. The predicted patient-level hazard scores in the testing folds achieved a C-index of 0.719 ± 0.098. Kaplan–Meier survival analysis, based on stratifying patients into high- and low-risk groups by the median of their predicted hazard score percentiles, further revealed that the high-risk group had significantly worse prognosis than the low-risk group (P < 0.0001; C-index = 0.719 ± 0.098; Fig. 6L). Consis- tently, univariate Cox regression showed that the scSurvival- derived hazard scores were significantly associated with survival (P < 0.001), outperforming other conventional immune features (Supplementary Fig. S8A). The achieved fivefold CV C-index (0.719) also exceeded those of benchmark models using cell type fractions, pseudo-bulk expression, or the CXCL9/SPP1 ratio (Supplementary Fig. S8B). Together, these results demonstrate that the risk prediction model constructed by scSurvival using single-cell expression data can provide clinically meaningful sur- vival predictions in real-world patient cohorts.discussionIn this study, we develop scSurvival, a new framework that simultaneously constructs survival prediction models and identifies cell subpopulations associated with disease progression and patient outcomes. Unlike existing approaches that rely on bulk data as an intermediary for indirect associ- ation or that simply aggregate single-cell expression into pseudo-bulk profiles or cell type fractions, scSurvival adopts a multiple-instance learning strategy, modeling each patient’s tumor profile as an ensemble of multiple cellular instances to enable cellular resolution. By introducing attention mech- anisms and a multiple-instance Cox regression framework, scSurvival not only enables efficient and accurate hazard score prediction but also preserves cellular-level heterogeneity, allowing for the identification of key cell subpopulations closely associated with survival outcomes. Such survival- associated cell subpopulations identified by scSurvival could also be used to facilitate other forms of survival analysis such as pseudo-bulk analysis of those selected cell popula- tions to build multigene risk scores. Our results demonstrate that scSurvival achieves robust performance across multiple simulated datasets and real clinical cancer cohorts. Notably, it successfully identified high-risk macrophage/monocyte subpopulations in the melanoma immunotherapy dataset and specific tumor cell subpopulations in the liver cancer cohort, providing valuable insights into the tumor immune microenvironment (TIME) and prognostic mechanisms. Moreover, a survival prediction model constructed using T cells from the melanoma cohort exhibited robust pre- dictive accuracy when validated on independent datasets, further confirming the stability and broad applicability of scSurvival. Interestingly, our analysis revealed that a subset of CD8+ cytotoxic T cells was unexpectedly associated with poorer prognosis in this cohort. Although counterintuitive, this likely reflects heterogeneous CD8+ functional states, in which exhausted or stress-adapted subsets can override their canonical effector role (35). Similar context-dependent associations of CD8+ T cells with unfavorable outcomes have been reported in certain cancer settings, underscoring the importance of considering their functional state rather than abundance alone when evaluating prognostic impact (24, 35). In addition, we observed context-dependent asso- ciations of B cells across tumor types. In melanoma, B cells were enriched in the lower-hazard group and frequently localized within TLS, consistent with their reported role in enhancing responses to ICB. By contrast, in liver cancer, B cells were enriched in the higher-hazard group, which may reflect the tolerogenic hepatic microenvironment in which B cells often acquire immunoregulatory phenotypes that suppress antitumor immunity. These findings highlight the heterogeneous and context-dependent roles of B cells in cancer prognosis (45, 46).Figure 6. (Continued) f, GO Biological Process (GOBP) enrichment analysis of DEGs in high-hazard versus low-hazard tumor cells. Adjusted P values (P_adjust) calculated using the hypergeometric test and corrected for multiple testing using the FDR. g, GSEA of high-hazard tumor cells versus low-hazard tumor cells. H, GSEA enrichment curves of pathways altered in higher-hazard tumor cells versus lower-hazard tumor cells. i, Monocle 2 trajectory analysis of tumor cell differentiation, revealing two main divergent trajectories. Cells are colored by pseudotime (top) and by risk groups (bottom). Pearson correlation analysis between pseudotime scores and pathway activity or adjusted hazard score of cells. Smooth curves in each plot, generated using locally weighted smoothing method. J, Boxplots showing the distribution of scLCSS activity scores [inferred by the Virtual Inference of Protein Activity by Enriched Regulon analysis (VIPER) algorithm] across tumor and nontumor samples in five independent cohorts. K, Kaplan–Meier survival curves showing OS or PFS based on unfavorable survival signature scores in two liver cancer datasets. Patients were stratified into high- and low-risk groups according to the median risk scores. Statistical significance was determined using the two-tailed log-rank test. L, Five-fold CV at the patient level, evaluating the stratification ability of risk percentiles and reporting the mean C-index and SD. P value calculated using the log-rank test. LIHC, liver hepatocellular carcinoma; NES, normalized enrichment score; TCGA, The Cancer Genome Atlas; Treg, regulatory T cells. ### Page 14 RESEARCH ARTICLEAACRJournals.orgRen et al.944 | CANCER DISCOVERY MAY 2026Nevertheless, several limitations remain to be addressed in future work. First, although scSurvival incorporates mech- anisms for batch effect correction, its robustness against extreme batch heterogeneity arising from different sequenc- ing platforms or centers could be further improved. Future efforts may focus on enhancing batch effect modeling to facilitate more effective cross-batch data integration. Second, although the SE module currently serves to balance feature scales, its biological significance has not yet been fully explored. Future work could investigate the interpretability of SE weights or introduce alternative strategies for identifying key survival-related genes, thereby enhancing the biological transparency of the model. Third, although scSurvival sup- ports the integration of clinical covariates at the patient level, we did not systematically evaluate how covariate adjustment influences the identification of risk-associated cell popu- lations. In our initial evaluation, including covariates in scSurvival for the melanoma datasets produced similar results. Further studies are needed to dissect the interactions between covariate effects and intrinsic cellular risk features. Lastly, emerging spatial transcriptomics and proteomics technol- ogies in cancer research provide rich molecular information coupled with spatial coordinates, offering an additional layer of biological context. The scSurvival framework could be extended to these spatial modalities by incorporating posi- tional relationships among cells or spots, thereby enabling the identification of survival-associated, spatially co-localized pat- terns that extend beyond cell subpopulations. In summary, scSurvival provides a scalable and systematic framework for survival prediction from cohort-level cancer single-cell data. It not only enables accurate prediction of patient survival outcomes but also provides a powerful tool for identifying survival-associated cell subpopulations to facilitate downstream mechanistic investigations at cellular resolution. By preserving cellular-level heterogeneity and leveraging attention-based modeling, scSurvival enhances our ability to link molecular and cellular features to patient out- comes. We believe that scSurvival will accelerate the broad adoption of cohort-level single-cell profiling for survival analysis and deepen our understanding of how specific cell subpopu- lations in the TME influence cancer outcomes, positioning it as a valuable tool for translational cancer research.MethodsscSurvival FrameworkThe input to scSurvival originates from a patient cohort (Fig. 1A), with a total of m patients. Each patient undergoes scRNA-seq, result- ing in an individual gene expression matrix (Fig. 1B). The collection of all single-cell expression matrices was denoted as ( ) { } 1m i i X X == , in which the expression matrix for patient i is represented as ( ) ii n p X × ∈R , comprising in cells and p genes. Each patient is also associated with survival information ( ),Y T E= , in which ( )1, 2, , m mT t t t += … ∈Rrecords the event times and ( ) { }1, 2 , , 0,1 mmE e e e= … ∈ indicates the event statuses, with 1ie = denoting the occurrence of the event (e.g., death) and 0ie = representing censoring (Fig. 1C; refs. 47, 48). To account for potential batch effects arising from different data sources, scSurvival also supports automatic batch effect correction, in which the batch labels corresponding to individual cells are denoted as ( ) { } 1m i i B b == , in which ( ) ii n b ∈R . The scSurvival model primarily consists of two modules: a cell feature extraction module and a multiple-instance Cox regression module. The feature extraction module is a VAE-based generative model that learns unified cell embeddings across batches (Fig. 1D). In the multiple-instance Cox regression module, these cell embed- dings are aggregated at the patient level through a multihead atten- tion mechanism and subsequently passed into a scoring network to compute relative hazard scores and perform standard Cox regression (Fig. 1E). If patient-level covariates need to be considered, scSurvival also supports incorporating covariates simultaneously during Cox regression (Fig. 1G). In addition, we apply an entropy regularization (22) term on the attention weights to control the sparsity of attention and enable the identification of key cells. The overall loss function of the model is given by1 2,cox VAE Entλ λ= + +L L L Lin which coxL , VAEL , and EntL denote the loss functions for Cox regres- sion, VAE reconstruction, and entropy regularization, respectively, and 1λ and 2λ are hyperparameters balancing these components. Ultimately, scSurvival is capable of providing survival predictions at the patient level, and based on the attention weights and the hazard scoring network, it can also identify key cells associated with higher or lower risk (with optional adjustment for covariate effects; Fig. 1I–K).VAE Model ArchitectureDuring the feature extraction phase, we train the model based on the complete set of single-cell expression data from all patients. In addition to standard quality control, preprocessing of the single-cell data before input to scSurvival includes total normalization within each cell, log1p transformation, and selection of highly variable genes (by default, scSurvival uses MVG2000 as the input; refs. 49, 50). Sep- arating these preprocessing steps from the model facilitates reducing the model burden when making predictions on independent data- sets. We assume that the preprocessed expression data are generated according to the following process:( ) ( )0, ,c cp z z I= N( ) ( )2 | , | , , ,cg c c cg cg cg gp x z b ZIG xρ μ σ=( ), ; ,g cg c cf z bρ ρ ϕ=( ), ; .g cg c cf z bμ μ ϕ=Here, modeld cz ∈R represents the low-dimensional latent state of cell c, in which modeld denotes the embedding dimension. Its prior dis- tribution is assumed to follow the standard normal distribution,( )0,cz IN . The gene expression cgx is modeled by a ZIG distribu- tion (17, 18) with parameters cgρ , cgμ , and 2gσ , in which cgρ is the inflation probability, and cgμ and 2gσ are the mean and variance of the nonzero Gaussian component, respectively. The ZIG distribu- tion is capable of capturing technical dropouts in gene expression data and distinguishing them from biologically meaningful low expression. Its likelihood function is defined as( ) ( ) ( ) ( )2 2 | , , 1 | , ,cg cg cg g cg cg cg cg cg gZIG x x xρ μ σ ρ δ ρ μ σ= + − Nin which ( )xδ is the Dirac delta function (51), which equals 1 at zero and 0 elsewhere. The parameters cgρ and cgμ are computed by a param- eterized decoder ( ), ;c cf z bϕ• , which takes as input the latent state czand the corresponding one-hot encoded batch label cb . The decoder is implemented as a fully connected neural network with parame- ters denoted byϕ . The variance term 2gσ is a trainable, gene-specific parameter independent of the input. ### Page 15 RESEARCH ARTICLESingle-Cell Survival AnalysisMAY 2026 CANCER DISCOVERY | 945According to the principle of VAE (arXiv 1312.6114), the posterior distribution of the latent state cz is assumed to follow a Gaussian dis- tribution of the form( ) ( ) ( )( )2 | , | , , , ,c c c c c c c cq z x b z x b x bθ θμ σ= Nin which the parametersθμ and 2θσ lie in modeldR and are computed by an encoder that takes as input the expression vector cx of cell c and its encoded batch label cb . In practice, the latent state cz is sampled via the following differentiable reparameterization trick:( ) ( ) ( ), , , 0, .c c c c cz x b x b Iθ θ μ σ ε ε= + ∼⊙ NThe encoder consists of three modules: LN, SE, and a feature extraction module. LN (19) normalizes the input during training. The SE module, short for squeeze-and-excitation module (20), is a lightweight attention mechanism defined as( ) ( )( )SE sigmoid ·relu · ,e sx W W x x= ⊙in which d p sW × ∈R projects the high-dimensional input features into a lower-dimensional space for global feature summarization, with d being the reduced dimension and p the original feature dimen- sion (i.e., the number of genes). The matrix p d sW × ∈R maps the com- pressed features back to the original space to produce feature-wise reweighting coefficients. The SE module enables adaptive recalibra- tion of features during training. The feature extraction module within the encoder is implemented as a fully connected neural network, denoted as ( ), ;Eθ• • • with param- etersθ . It takes as input the concatenated recalibrated gene features and batch labels to compute the parameters of the latent posterior distribution:( ) ( )( )( ), SE LN , ; ,c c c cx b E x bθ μ μ θ=( ) ( )( )( ), SE LN , ; .c c c cx b E x bθ σ σ θ=VAE Reconstruction Loss and Rare Gene WeightingFollowing the principle ofβ -VAE (52), we define the reconstruc- tion loss using the negative evidence lower bound as( ) ( ) ( ) ( )( )1 log | , ,p VAE c cg KL c c c c g L x ZIG x D q z x b p zβ = = − + ‖in whichβ is a hyperparameter controlling the strength of the Kullback-Leibler (KL) divergence (53) regularization. During training, we further employ the FreeBits technique (54) to control the influ- ence of KL divergence and prevent feature collapse, by imposing a minimum KL divergence thresholdτ independently on each latent dimension. In addition, to stabilize the variance estimates, we intro- duce a prior regularization (55) on 2gσ as( ) ( )2 2Gamma 2,2 ,g gpσ σ=and incorporate it into the VAE loss via maximum likelihood. The resulting overall loss function becomes( ) ( ) ( ) ( )( ) ( )2 1 1 log | , log ,p p VAE c cg KL c c c c g g g L x ZIG x D q z x b p z pσ β λ σ = = = − + − ‖in whichσλ is a hyperparameter controlling the strength of the variance prior regularization. Moreover, to balance the expression levels across genes and pre- vent important lowly expressed rare genes from being dominated by highly expressed genes, we introduce gene-specific weights into the loss function. The gene weight is defined by1 1 , for 1,..., ,n g cg c x x g p n = = =( )1 2median , .., ,px x x x=ɶ,g g x xα ω −+ =  + ɶε εin which n is the total number of cells andε is a small constant to avoid division by zero, andα is a hyperparameter controlling the strength of weighting differences. Thus, the final VAE training loss becomes( ) ( ) ( )( ) ( ) ( )( )2 1 log log | , .p VAE c g cg g KL c c c c g L x ZIG x p D q z x b p zσ ω λ σ β = = − + + ‖The average loss over all cells is computed as( )1 1 ,n VAE VAE c c L x n = = Lin which 1m ii n n==  denotes the total number of cells across all patients. In addition to extracting cell features using ZIG-VAE, we also im- plement a standard VAE with mean-squared error as the reconstruc- tion loss, given by( ) ( ) ( ) ( )( )2 1 | , ,pG VAE c g cg cg KL c c c c g L x x D q z x b p zω μ β− = = − − + ‖and PCA (56) for feature extraction. These alternative methods can handle simpler scenarios and provide rapid performance validation in practical analyses.Multi-Head Attention Aggregation MechanismWe use the mean vector ( ),c cx bθμ output by the decoder as the sta- ble embedding feature for each cell and subsequently use it as the input to the downstream model. For simplicity, we denote it as cμwhen there is no ambiguity. We then employ a multi-head attention mechanism (57) to aggregate the cell embeddings into patient-level features. In each attention head, the cell embedding is first projected to a lower-dimensional space through a linear mapper:,c h c hh W bμ= +in which h modeld d hW × ∈R and hd hb ∈R are the weight matrix and bias vector, respectively, and model h h d d n= is the reduced feature dimension, with hn denoting the number of attention heads. Next, we compute the gated attention coefficients (12) as( ) ( )( )tan sigmoid ,T c h h c h ca w h U h V h= ⊙in which hd hw ∈R and , h hd dh hU V × ∈R are trainable parameters. Subsequently, the dimension-reduced features and corresponding attention weights are grouped by patient, denoted as ( ) ( ) ( ) ( ){ }1 2 1 , ..., i m i i i n i h h h =and ( ) ( ) ( ) ( ){ }1 2 1 , ..., i m i i i n i a a a =, respectively. For each patient, the attention coefficients are first normalized via the SoftMax function:( ) ( ) ( ) ( ) ( )1 exp , for 1,.., , for 1,.., . expi i ci c in i cc a a c n i m a= = = =ɶThe patient-level feature is then obtained by the weighted aggrega- tion of the cell features as( ) ( ) ( )1 ,for 1 ,.., .in i i i c c c H a h i m= = =ɶFinally, the complete patient representation is formed by concate- nating the outputs from all attention heads as( ) ( ) ( ) ( ) ( )1 2H Concat , ,..., , for 1,.., .h i i i i nH H H i m= =in which ( )i jH denotes the aggregated patient feature from the jth attention head. ### Page 16 RESEARCH ARTICLEAACRJournals.orgRen et al.946 | CANCER DISCOVERY MAY 2026Although the gated attention mechanism inherently tends to pro- duce sparse attention, in practice, to better identify the most crit- ical hazard cells, we further introduce entropy regularization (22) to explicitly control the sparsity of attention. First, the normalized information entropy for each patient is computed as( ) ( ) ( ) ( ) ( ), ,1 1 log1 , logh i i j i j n n c ci j ch i a a Ent n n= = = −  ɶ ɶin which ( ),i j caɶ denotes the normalized attention coefficient for cellc in patient i under the jth head. The entropy regularization loss is then defined as( )1 1 max ,0 ,m i Ent i Ent mγ = = − −   Lin whichγ is a threshold for the normalized entropy, such that no penalty is applied when the entropy decreases belowγ .Cox RegressionThe aggregated patient-level features are subsequently used for standard Cox regression (47, 48). We define the hazard model as( ) ( ) ( ) ( ) ( )( )0H exp H ,i i h t h t s=in which ( )0h t denotes the baseline hazard and ( ) ( )H i s is the relative hazard score. The scoring function s is implemented as a simple two- layer neural network with the following structure:( ) ( )( )2 0.5 1 1 2· dropout relu · ,o h W W h b b= + +( ) ( )( )clamp , 10,10 ,s h o h= −in which 0.5dropout denotes dropout with a rate of 0.5 and clampconstrains the output within the range [ ]10,10− to improve numer- ical stability. The loss function is based on the partial likelihood of the Cox pro- portional hazards model, defined as( ) ( ) ( ) ( ) ( )( ): 1 1 H log exp H ,i i i i cox i e i R te s s m = ∈ ′ ′ = − −      Lin which 1m e i i m e= =  denotes the total number of patients who experienced the event, ie and it represent the event status and event time for patient i, respectively, and ( ) { }:i i jR t i t t′′= ≥ denotes the risk set, comprising patients whose survival time is greater than or equal to it . When incorporating patient-level covariates, we combine the rel- ative hazard contributions from the covariates using a linear model. The hazard model is modified as( ) ( ) ( ) ( ) ( ) ( ) ( )0H exp H ,i i iT covh t h t s Xη= +in which ( )i covX denotes the covariate vector for patient i andη rep- resents the trainable weight vector associated with the covariates. The corresponding loss function becomes( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ): 1 1 H log exp H .i i i i i iT T cox cov cov i e i R te s X s X mη η′ ′′= ∈ = − + − +      LTotal Loss Function and Training StrategyscSurvival trains the model by jointly optimizing the three loss components described above, with the total loss function defined as1 2,cox VAE Entλ λ= + +L L L Lin which 1λ and 2λ are hyperparameters. In practice, to achieve sta- ble cell embeddings, we adopt a pretraining followed by fine-tuning strategy. First, we pretrain the VAE by minimizing VAEL alone using batch gradient descent. Subsequently, based on the pretrained VAE parameters, we fine-tune the model by jointly optimizing the overall loss L using the gradient surgery technique (58). The principle of gra- dient surgery is to apply different learning rates to different modules, thereby preserving the pretrained feature extractor’s stability. During fine-tuning in scSurvival, we set the learning rate for the VAE module to be 1/10 of that for the multiple-instance Cox regression module. This joint optimization strategy enables the model to learn cell representations that are adaptively aligned with survival outcomes, yielding more informative and task-specific embeddings than fixed features while simultaneously mitigating overcorrection of batch effects. Nevertheless, scSurvival also provides interfaces that allow users to directly input custom-fixed features if desired. When handling extremely large-scale datasets (e.g., when the total number of cells exceeds one million), computing the Cox loss requires aggregating features from all cells, which in turn forces the VAE module, especially the decoder, to reconstruct all cells simulta- neously. Because of the decoder’s complexity and the reconstruction loss, this would result in constructing and storing extremely large computational graphs, causing significant memory overhead. To address this, for ultra-large datasets, we adopt an alternating optimization strategy during fine-tuning to approximate joint train- ing. Specifically, the steps are as follows: 1. Perform one gradient update step on 1 VAEλ L . As only the VAE is updated in this step, we can employ mini-batch training as in the pretraining phase, which automatically releases the computational graph after each update and prevents memory accumulation. To further accelerate training, it is even possible to update using only a single batch. Empirically, this procedure is sufficient to maintain the stability of the pretrained model in large-scale settings. 2. Perform one full-batch gradient update on 2cox Entλ+L L . Given the relative simplicity of the multiple-instance Cox module, full-batch updates remain manageable in terms of memory consumption. 3. Within each training iteration, execute steps 1 and 2 sequentially to implement alternating updates.Cell Number Balancing SamplingWhen the number of measured single cells varies significantly across patients, directly training on all data may lead to uneven dis- tributions of attention coefficients among patients. For example, patients with a larger number of cells may have overall smaller atten- tion values, introducing biases in cell-level comparisons. To address this issue, we introduce an optional cell number–balancing sampling strategy in the multiple-instance Cox regression module. Specifi- cally, we first compute the average number of cells across all patients, denoted as n. During each fine-tuning step, for patients with more than n cells, we perform uniform sampling without replacement; for patients with fewer than n cells, we perform uniform sampling with replacement. As a result, each patient retains exactly n cells for the computation of the Cox loss and entropy regularization loss, fol- lowed by gradient updates. As the sampling procedure is repeated at every training step, eventually, almost all cells participate in the train- ing process without causing data resource waste.Hyperparameter SettingsThe thresholdγ in the entropy regularization is the most criti- cal hyperparameter in scSurvival as it directly controls the sparsity of the attention coefficients and thus influences the number of key cells retained. We recommend adjustingγ within the range of 0.5 to 1. Other hyperparameters are kept consistent across all our experi- ments on both simulated and real datasets and can be categorized into model hyperparameters and training hyperparameters. The model hyperparameters include the input feature dimension, which is set to 2000 by default; in all analyses, we use MVG2000 genes selected ### Page 17 RESEARCH ARTICLESingle-Cell Survival AnalysisMAY 2026 CANCER DISCOVERY | 947by Seurat (49) as input features; the latent feature dimension modeld , which is set to 128; the hyperparametersβ andσλ for the VAE recon- struction loss, with default values of 1 and 0.1, respectively; the KL divergence thresholdτ for FreeBits (54), which is set to 0.2; the gene weight scaling factorα , which is set to 0.2; the number of attention heads headn , which is set to 8; and the loss weight hyperparameters 1λand 2λ , which are set to 0.01 and 1, respectively. The training hyperparameters include a learning rate of 3 1 1 0− × ; 200 pretraining epochs; and a maximum of 500 fine-tuning epochs, with early stopping if no further improvement in the loss function or validation metrics is observed. We use the Adam optimizer during pretraining and switch to AdamW during fine-tuning, with a weight decay of 0.01. For small-scale datasets, the pretraining batch size is set to the maximum number of cells among all patients; if balanced sampling is enabled, the batch size is set to the sampling number (i.e., the average number of cells). For large-scale datasets, the batch size during VAE training can be appropriately increased depending on GPU memory availability, with a default value of 4,096.Cell Importance and Cell HazardAt the cell level, the attention coefficients reflect which cells play a critical role in risk estimation. Given the large number of cells, the values of the attention coefficients tend to be small after SoftMax normalization. Therefore, we define the importance of each cell based on the following form of attention coefficients:( ), ,sigmoid , for 1,..,ˆ ,c j c j heada a j n= =( ), 1 ˆ ˆ.head c c j j n a max a≤ ≤ =in which ,c ja denotes the unnormalized attention coefficient of cellc at the jth attention head. The sigmoid function rescales the raw attention values to the range (0–1), providing an interpretable mea- sure of each cell’s individual contribution to the model. Furthermore, to evaluate the relative hazard score at the cell level, we first align each cell’s feature with the patient-level representation( ),1 ,1 , , , , ˆ ˆ ˆConcat , ..., , .. , ,ˆ. head headc c c c j c j c n c nh h a h a h a=in which ,c jh represents the linearly projected feature of cell c from the jth attention head and , ˆc ja selectively preserves the parts of the feature that make significant contributions to the model. Subse- quently, the trained hazard scoring network is used to predict a rela- tive hazard score for each cell as( )ˆ.c cs s h=Although the scoring function s is trained at the patient level, because of the enforced sparsity of the attention coefficients, the aggregated patient feature is effectively close to the features of a few key cells. As a result, the scoring function can still provide reasonable hazard estimates at the cell level. This cell-level hazard score, modulated by the attention weights, effectively reflects the contribution-adjusted risk of each cell and is therefore referred to as the adjusted hazard (abbreviated as hazard_adj). Throughout this study, unless otherwise specified, all reported cell- level hazard scores refer to these attention-adjusted values.Model PredictionWhen there is no significant batch effect within the training dataset or between the training and testing datasets, we can ignore the batch variables during model training and directly make reason- able predictions on the test set. However, when substantial batch effects are present, prediction performance may deteriorate notice- ably. The feature extraction module in scSurvival does not require clinical information from patients and can operate independently from the multiple-instance module. Therefore, scSurvival supports incorporating the test set into the feature extraction module during training, allowing the model to learn cross-dataset repre- sentations and enabling transfer learning, thereby improving pre- diction performance.Evaluation MetricsAt the patient level, we typically use the C-index (59) to evaluate the model’s predictive performance, defined as( ) ( ) ( ) ( ), ,1 2C index ,i i i ii i i i I s s I s s′ ′∈′ ∈ ′ > + = =  C CC‐in which ( ){ }, , 1, i i ii i i i e t t ′= ≠ = <′ ′C is the set of all comparable patient pairs, is denotes the risk score for patient i, it is the survival time for patient i, and ( )·I is the indicator function. When the test set contains a sufficiently large number of samples, directly computing the C-index on the test set is appropriate. How- ever, when the test sample size is small, forming enough comparable pairs may be challenging. This situation commonly arises during CV on small datasets. Therefore, we adopt the following procedure to compute evaluation metrics: 1. In each validation fold, predict the hazard scores for the test samples while also computing the hazard scores for the training samples. 2. To eliminate prediction bias across folds, transform the test sam- ple hazard scores into quantiles relative to the training set risk score distribution. 3. Collect the quantile-transformed hazard scores across all valida- tion folds and compute a unified C-index as the validation metric. Note that all training steps, including the selection of most vari- able genes, are performed exclusively within the training set of each validation fold to avoid information leakage. During training, we also introduce a variant of the C-index, termed the conditional C-index (CC-index), to monitor early stop- ping and improve generalization performance. The CC-index is defined as( ) ( ) ( ) ( ) ( ), ,1 2CC index A ,i i i ii i i i I s s I s s′ ′′ ′∈ ′ ′∈ > + = = ′ ‐ C CCin which A A= −′ ×C C excludes all comparable pairs formed solely within set A. For evaluation, we merge the hazard scores from the validation and training sets and compute the CC-index by set- ting A as the training set. The CC-index evaluates not only the rank- ing accuracy within the validation set but also the relative ranking of validation samples within the combined dataset. As the CC-index incorporates comparable pairs between the validation and training sets, it remains applicable even for small validation sets. Note that the CC-index is only used for model selection during training and is not used for final performance reporting. Additionally, we stratify patients based on the hazard score into high-risk and low-risk groups and evaluate the model risk stratifica- tion capability using Kaplan–Meier survival curves (60) and the log- rank test (61). At the cell level, we classify cells into three groups based on their attention scores and hazard scores: inattentive, higher risk, and lower risk. Specifically, we first cluster cells into two groups based on their attention scores, labeling the low-attention group as inattentive. Within the high-attention group, we further divide cells into high- er-risk and lower-risk groups based on whether their hazard scores are positive or negative. As ground-truth labels are available in the simulated datasets, we evaluate the performance of cell classification only on simulated data. We use confusion matrices, precision, recall, and F1-scores (62) to assess the agreement between the predicted cell group labels and the ground-truth labels. ### Page 18 RESEARCH ARTICLEAACRJournals.orgRen et al.948 | CANCER DISCOVERY MAY 2026Basic Simulation Data SetupTo simulate a cohort of single-cell datasets with survival-driven cell subpopulations, we first used Splatter (23) to generate a synthetic single-cell expression dataset containing 10,000 cells and 5,000 genes. The cells were simulated as belonging to three groups, with assign- ment probabilities set to 0.7, 0.15, and 0.15, respectively. The differen- tial expression probabilities for the three groups were set to 0.2, 0.06, and 0.06, respectively. Among these, the two smaller groups were des- ignated as cell subpopulations associated with good prognosis and poor prognosis, respectively, whereas the largest group was treated as background cells. Next, we generated a cohort of single-cell expres- sion matrices along with corresponding survival times by sampling from this synthetic dataset. To ensure that the simulated survival outcomes were driven by the designed cell subpopulations, we mod- eled a dynamic progression: as survival time advances, the number of cells associated with poor prognosis decreases, whereas the number of cells associated with good prognosis increases; the number of back- ground cells remains approximately constant throughout. Suppose we simulate m samples with monotonically increasing survival times. For the i-th sample, the corresponding single-cell expression matrix is generated according to the following steps: 1. Set the sampling ratio 1i r m= − . Randomly sample a proportion rof cells from the good-prognosis group and a proportion 1 r− of cells from the poor-prognosis group; 2. Mix the sampled cells with the background cells and then ran- domly sample 1,000 cells from the combined pool to form the single-cell expression matrix for sample i; 3. Assign the event status to 1 and the survival time to i for sample i. As the probabilities of the two survival-driving cell groups are initially balanced, their total cell numbers are approximately equal. Through the above sampling procedure, the total number of mixed cells remains nearly constant for each sample, whereas the proportion of good-prognosis and poor-prognosis cells changes progressively, and the proportion of background cells remains stable. In addition, for each sample, we simulated censoring with a probability of 0.1. If censoring occurs, we set the event status to 0 and randomly advance the survival time (i.e., the censoring time). Based on this process, we generated a cohort of 100 samples with associated single-cell expression matrices and survival outcomes for evaluating the performance of scSurvival in identifying survival-driving cells and predicting survival outcomes. To further validate the model’s ability to integrate single-cell data across different batches, we also simulated patients originating from different batches. Specifically, using Splatter, we simulated two batches of data containing 6,000 and 4,000 cells, respectively, each with 5,000 genes. Each batch contained three cell clusters corre- sponding to good prognosis, poor prognosis, and irrelevant (back- ground) groups. We then applied the same sampling strategy as above to generate 50 samples from each batch, each with single-cell expres- sion matrices, survival times, and censoring statuses.Extended Simulation SetupTo comprehensively evaluate the robustness and generalizability of scSurvival under various biological and technical conditions, we designed six types of extended simulation scenarios (Sim.1–Sim.6) based on the basic settings. Each scenario introduced a specific per- turbation mechanism, and multiple parameter settings were set per scenario to simulate varying degrees of complexity and noise. 1. Risk cell percentage (Sim.1): We varied the proportions of good. survival and bad.survival cells in the single-cell reference dataset used for patient sampling, thereby generating datasets with dif- ferent sizes of ground-truth risk cell populations. This setting was used to assess the stability of scSurvival when identifying rare risk subpopulations. 2. Batch imbalance (Sim.2): To simulate risk signal imbalance across batches, we generated two batches of single-cell data and downs- ampled good.survival or bad.survival cells asymmetrically between batches. The downsampling ratio controlled the degree of imbal- ance, with extreme settings in which one risk cell type was present in only a single batch. This setting tested the model’s ability to handle risk feature shift across batches. 3. Patient outcome distribution (Sim.3): To simulate more real- istic, noise-prone outcome distributions, we modeled patient survival times using an accelerated failure time (AFT) frame- work with a Weibull baseline distribution (63). The proportion of risk-associated cells, denoted as s, was used as the primary covariate. Specifically, the survival time T for each patient was generated according to the following formula:( ) ( )1 exp · · · ln , IQRs T U sγ β α   = − −       in which ( )IQR s denotes the IQR of s across all patients, used to normalize the effect scale across different settings; ( )ln TRβ = , in which TR is the time ratio, defined as the median survival time ratio between high- and low-risk patients, indicating the extent to which a covariate influences expected survival time;α is the scale parameter of the survival distribution, fixed at 10;γ is the shape parameter, set to 1.5 to simulate higher risk in later stages of disease progression, such as cancer; and U is a random variable sampled from ( )Uniform 0.2,0.8to introduce stochasticity into survival times. We tested three levels of effect size, corresponding to TR values of 1.3, 2, and 3, to reflect varying degrees of biological signal strength. Additionally, we included a noiseless idealized setting ( infTR = ), in which U was fixed at 0.5 to yield deterministic survival times. Moreover, to simulate event censoring, we introduced random censoring under three censoring probabilities: 0.2, 0.5, and 0.7. When censoring occurred, the censoring time was uniformly sampled between 0 and the simulated survival time generated by the AFT–Weibull model. 4. Cell composition heterogeneity (Sim.4): To evaluate the robustness of scSurvival under complex cellular compositions, we simulated more complex celltype compositions for each patient. Specifically, the proportions of background cell types were randomly varied across patients to introduce heterogeneity and noise. The degree of cell composition heterogeneity was controlled by adjusting the number of distinct background cell types, allowing us to assess whether scSurvival could consistently identify true risk-associated cells under varying levels of complexity. 5. Patient group heterogeneity (Sim.5): Beyond heterogeneity arising from cell composition, we also considered interpatient variability in how individuals respond to the same biological signals. For instance, elderly patients may exhibit stronger clinical responses to the same level of risk signal compared with younger patients. To simulate this group-wise covariate effect, we extended the AFT–Weibull model to incorporate a binary grouping variable as( ) ( ) ( )1 exp · · · ln , IQRs T g U sγ β δ α   = − + −       in which { }0,1δ ∈ indicates whether the patient belongs to a high-response group (e.g., elderly) and g represents the magnitude of effect shift between the two groups. We varied both the proportion of patients assigned to the high-response group and the value of g to systematically evaluate the ability of scSurvival to model patient-level heterogeneity, and accordingly, this simulation scenario was analyzed using the covariate-aware mode of scSurvival. 6. Risk subcluster heterogeneity (Sim.6). To highlight the advantage of scSurvival in performing survival analysis at true single-cell res- olution, we introduced a challenging simulation scenario in which patient survival time originates only from a small subpopulation ### Page 19 RESEARCH ARTICLESingle-Cell Survival AnalysisMAY 2026 CANCER DISCOVERY | 949embedded within a larger cell cluster. We varied two key param- eters: (i) the proportion of ground-truth risk-driving cells within the larger cluster and (ii) the degree of separation between the risk subpopulation and the rest of the cluster, achieved by adjusting the log fold change of DEGs in Splatter (23). This setting was designed to assess whether scSurvival can still effectively model survival outcomes and accurately identify relevant subpopulations even in cases in which conventional clustering algorithms might fail because of the rarity or subtlety of the risk signal. Finally, based on the simulation approach described in Sim.3, we fixed the censoring rate at 0.3 and generated single-cell cohorts with varying numbers of patients and cells per patient. This allowed us to evaluate, under different effect sizes, how many patients are required for scSurvival to reliably identify risk-driving cell subpopulations. Besides, for benchmarking and visualization purposes, we selected one representative setting from each simulation category, reflecting moderately challenging but nonextreme scenarios (Fig. 2G):• Sim.1: risk cell proportion = 3%;• Sim.2: downsampling rate = 0.2;• Sim.3: censoring rate = 0.7 and TR = 2;• Sim.4: number of background clusters = 5;• Sim.5: group effect 3g = and older group proportion = 0.5;• Sim.6: subcluster proportion = 0.3 and separation = 0.35.Benchmark StrategiesTo systematically evaluate the performance of scSurvival, we designed two complementary benchmarking strategies. The first benchmarking strategy focused on identifying the ground-truth risk-driving cell subpopulations in the simulated datasets. To isolate the impact of feature extraction, we fixed the downstream multi- instance Cox regression module and replaced only the upstream feature encoder. Specifically, we compared the default (ZIG-VAE with alternative feature representations, including PCA, NMF (imple- mented via sklearn.decomposition.NMF; ref. 64), and top 2000 highly variable genes (HVG2000) identified by Scanpy (RRID: SCR_018139; ref. 50), which were used directly without further dimensionality reduction. Note that all methods other than ZIG-VAE did not include any explicit batch effect correction modules. The resulting cell-level features were passed into the same multi-instance Cox model, in which risk-relevant cells were identified by jointly analyzing attention weights and estimated cell-level hazard scores. Model performance was evaluated by comparing predicted risk cell groupings against the known ground-truth labels using macro-averaged precision, recall, and F1-score under a three-class setting (good.survival, bad.survival, and other). The second benchmarking strategy aimed to evaluate patient-level survival prediction. We constructed Cox regression models using either cell type proportions or pseudo-bulk expression as input features, serving as baseline methods for comparison with scSurvival. All models were implemented using sksurv.linear_model.CoxnetSurvivalAnalysis (65) and evaluated via nested CV, with a fivefold outer loop to assess model performance and an inner loop for tuning the regularization parameter (α). For the cell type proportion model, the input consisted of the relative abundance of major cell types per patient. To approxi- mate ridge regression and improve stability under feature collinearity, the l1_ratio was set to 0.1 in simulation datasets and 1e−3 in real-world data. Cell types in simulation were derived using the Leiden clustering algorithm (66). For the pseudo-bulk model, the input was the average gene expression profile per patient, and the l1_ratio was set to 1 (lasso regression) to facilitate feature selection. In real-data analysis, we addi- tionally constructed a univariate Cox model based on the CXCL9/SPP1 expression ratio, a widely used immune-related prognostic marker. All models were evaluated by calculating the C-index on the outer CV loop to assess risk prediction accuracy.GSEADifferential expression analysis between higher-hazard and lower-hazard cells was carried out with significance thresholds of|log 2 fold change| > 0.5 and adjusted P < 0.05. As previously described (67), we conducted GO, KEGG, and rank-based GSEA using the R package clusterProfiler (v4.7.1; RRID: SCR_013042). For bulk tran- scriptomic data, we employed the single-sample GSEA method with the R package GSVA. For scRNA-seq data, the AUCell algorithm was employed to mitigate the impact of dropouts and inherent technical variability. Gene sets used in this study, including Hallmark (v7.2), GO (v7.2), and KEGG (v7.2), were retrieved from the Molecular Signa- tures Database (v7.2; RRID: SCR_016863).Kaplan–Meier Survival Curve ComparisonAs previously described (9), we conducted ordinary survival anal- yses to assess the clinical relevance of the scSurvival-derived gene signatures. The resulting survival signature comprised the 100 most upregulated and the 100 most downregulated DEGs from higher- hazard versus lower-hazard cells. Signature scores were computed using the Virtual Inference of Protein Activity by Enriched Regulon analysis algorithm (68), which integrates both the upregulated and downregulated scSurvival-derived gene sets. Based on the median value of the signature scores, samples were dichotomized into high- score and low-score groups. Kaplan–Meier survival curves were gen- erated for these groups, and statistical significance was determined using the log-rank test.Public Dataset Collection and ProcessingIn this study, we curated four publicly available human scRNA-seq datasets and one spatial transcriptomic dataset for both discov- ery and validation (Supplementary Table S1). For the PRJNA679099 scRNA-seq dataset, raw FASTQ files were downloaded from the NCBI BioProject database. Raw gene expression matrices were generated using CellRanger (v8.0.1; RRID: SCR_017344) with the refdata-gex- GRCh38-2020-A reference genome, following the standard 10x Genomics pipeline with default parameters unless otherwise speci- fied. The expression matrices for the other datasets were downloaded from their respective public repositories. After importing raw counts, we removed low-quality cells (fewer than 250 detected genes, fewer than 500 total transcripts, or more than 20% mitochondrial reads). Seurat (v 5.1.0; RRID: SCR_016341) was then applied to perform normalization, identification of variably expressed genes, and PCA. A shared nearest-neighbor graph built from the first 10 principal components was used to delineate commu- nities. Finally, all cells were visualized in two dimensions with UMAP. Human bulk transcriptomic datasets were obtained from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih. gov/geo/) and UCSC Xena database (https://xena.ucsc.edu/), and the detailed information is summarized in Supplementary Table S1. For the PRJEB23709 melanoma bulk RNA-seq dataset, raw FASTQ files were downloaded from the NCBI BioProject database. As described in the official documentation (69), raw gene expression matrices were generated using the Salmon pipeline (v1.4.0; RRID: SCR_017036) following the recommended preprocessing workflow (https://github. com/COMBINE-lab/salmon). The expression matrices for the other bulk transcriptomic datasets were downloaded from their respective public repositories. For all bulk RNA-seq data, the transcript per mil- lion value was calculated.Main scRNA-Seq Dataset DescriptionsIn the GSE120575 melanoma immunotherapy cohort, we analyzed 48 tumor samples from 32 patients, encompassing a total of 16,291 CD45+ immune cells profiled using an optimized full-length Smart- seq2 protocol. Cell type annotations were derived from the original ### Page 20 RESEARCH ARTICLEAACRJournals.orgRen et al.950 | CANCER DISCOVERY MAY 2026study (24). Based on radiologic assessments of immunotherapy response, the samples were classified into 17 responders and 31 non- responders (24). In the present study, this dataset was used as both a training set and discovery cohort for identifying survival-associated immune cell subpopulations and constructing predictive transcriptional signatures. The processed expression matrix was downloaded from GEO (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE120575), and only protein-coding genes were retained for all analyses per- formed in this study. The detailed preprocessing steps are described in the Methods section. Besides, we collected one spatial transcrip- tomic dataset (28) and three bulk RNA-seq datasets, including PRJEB23709 (36), GSE78220 (70), and phs000452.v3.p1 (71), to further evaluate the robustness of the derived signature. To assess the risk landscape of CD8+ T cells in melanoma under ICB, we analyzed scRNA-seq datasets from three patients with annotated treatment outcomes. Two nonresponding metastatic melanoma samples (K411 and K468) were obtained from GSE159252 (37). Quality control followed the original publication’s parameters in Seurat by excluding cells with fewer than three counts, fewer than 400 detected features, or more than 10% mitochondrial content. CD8+ T-cell subsets were extracted for downstream analysis, yielding 3,737 cells from K411 and 4,878 cells from K468. For the respond- ing case, we accessed the processed expression matrix and cell type annotations of sample GSE134388 (38) from the TISCH24 portal (https://tisch.compbio.cn/home/) and similarly extracted the CD8+T-cell subset. Hazard-adjusted values were projected onto UMAP for visualization and compared across patients. In addition, the PRJNA679099 melanoma scRNA-seq dataset with available survival information was used as an independent validation cohort (39). The PRJNA679099 dataset comprises scRNA-seq data from 13 immunotherapy-naïve patients with cutaneous melanoma. In this study, researchers performed 10x Genomics scRNA-seq on sorted immune cells from resected tumors and peripheral blood, pro- filing tumor-infiltrating lymphocytes. In the present work, this data- set was used as an independent validation cohort for evaluating the performance of the prognostic model. Raw FASTQ files were down- loaded from the NCBI BioProject database (PRJNA679099), and upstream processing was performed to generate expression matrices (see “Methods” section for details). Furthermore, we included the GSE221553 melanoma scRNA- seq dataset (33) as an another independent validation cohort. This dataset comprises longitudinal tumor samples from 13 patients with metastatic melanoma treated with tumor-infiltrating lymphocyte adoptive cell therapy. The study used CD45-sorted scRNA-seq to interrogate the TME across treatment time points. In our analysis, we downloaded raw count matrix from GEO, performed quality con- trol, normalization, and cell type annotation, and then applied our scSurvival pipeline to validate whether the molecular signatures, hazard-associated cell populations, and pathway enrichments identi- fied in our training melanoma dataset are recapitulated in this inde- pendent scRNA-seq dataset. The PRJCA007744 dataset comprises scRNA-seq data from 189 samples collected from 124 patients with liver cancer (40). This com- prehensive dataset was generated using the 10x Genomics platform and includes more than 1 million cells, encompassing various liver cancer subtypes such as hepatocellular carcinoma, intrahepatic chol- angiocarcinoma, and combined hepatocellular cholangiocarcinoma. In the corresponding study, researchers performed scRNA-seq anal- ysis to investigate the heterogeneity of the TIME in liver cancer (40). In the present study, this dataset was used as both a training set and discovery cohort for identifying survival-associated tumor cell sub- populations and constructing predictive transcriptional signatures. The processed raw count matrix was obtained directly from the authors (40), and we removed single-cell data from three patients with- out available survival information. Only protein-coding genes were retained, and ribosomal genes (RPL and RPS) and mitochondrial genes (MT-) were excluded. The detailed preprocessing steps are described in the “Methods” section. Because of the lack of publicly available scRNA-seq datasets with survival annotations, we addition- ally collected eight bulk transcriptomic datasets as external valida- tion cohorts (Supplementary Table S1), including GSE36411 (72), GSE14520 (73), TCGA-LIHC (74), GSE76427 (73), GSE40873 (75), and GSE116174 (76).Data AvailabilityThe public datasets analyzed in this study are described in the “Methods” section, with corresponding accession numbers pro- vided therein. For convenience, all accession numbers and dataset details are also compiled in Supplementary Table S1. The open- source scSurvival program and its tutorials are freely available at GitHub (https://github.com/cliffren/scSurvival), Zenodo (https://doi. org/10.5281/zenodo.15399777), and Code Ocean (https://codeocean. com/capsule/7948185/tree/v1).Authors’ DisclosuresG.B. Mills reports personal fees and nonfinancial support from Amphista Therapeutics, Astex Pharmaceuticals, AstraZeneca, BlueDot, Ellipses Pharmaceuticals, ImmunoMET Therapeutics, Ionis Pharma- ceuticals, Leapfrog Bio, Bruker/NanoString Technologies, NeoPhore, Nerviano Medical Sciences, Nuvectis Pharma, Pangea Pharmaceuticals, PDX Pharmaceuticals, Qureator, RyboDyn, SignalChem Lifesciences, Tarveda Therapeutics, Turbine, and Zentalis Pharmaceuticals; own- ership of BlueDot stock, Catena Pharmaceuticals stock, ImmunoMet Therapeutics stock, Nuvectis Pharma stock, RyboDyn stock, Signal- Chem Lifesciences stock, Tarveda Therapeutics stock, and Turbine stock; other support from Adelson Medical Research Foundation, Breast Cancer Research Foundation, Komen Research Foundation, Ovarian Cancer Research Foundation, and Prospect Creek Founda- tion; other support from AstraZeneca (sponsored research), Zentalis Pharmaceuticals (sponsored research), NanoString (sponsored research), and Ionis Pharmaceuticals (sponsored research); other support from AstraZeneca (clinical trial support), Genentech (clinical trial sup- port), GSK (clinical trial support), and Eli Lilly and Company (clin- ical trial support) during the conduct of the study; grants, personal fees, and nonfinancial support from Amphista Therapeutics, Astex Pharmaceuticals, AstraZeneca, BlueDot, Ellipses Pharmaceuticals, ImmunoMET Therapeutics, Ionis Pharmaceuticals, Leapfrog Bio, Bruker/NanoString, NeoPhore, Nerviano Medical Sciences, Nuvectis Pharma, Pangea Pharmaceuticals, PDX Pharmaceuticals, Qureator, RyboDyn, SignalChem Lifesciences, Tarveda Therapeutics, Turbine, and Zentalis Pharmaceuticals; ownership of BlueDot stock, Catena Pharmaceuticals stock, ImmunoMet Therapeutics stock, Nuvectis Pharma stock, RyboDyn stock, SignalChem Lifesciences stock, Tarveda Therapeutics stock, and Turbine stock, other support from Adelson Medical Research Foundation, Breast Cancer Research Foundation, Komen Research Foundation, Ovarian Cancer Research Foundation, and Prospect Creek Foundation, other support from AstraZeneca (sponsored research), Zentalis Pharmaceuticals (sponsored research), NanoString (sponsored research), and Ionis Pharmaceuticals (spon- sored research), and other support from AstraZeneca (clinical trial support), Genentech (clinical trial support), GSK (clinical trial sup- port), and Eli Lilly and Company (clinical trial support) outside the submitted work. No disclosures were reported by the other authors.Authors’ ContributionsT. Ren: Data curation, software, formal analysis, validation, investigation, visualization, methodology, writing–original draft, writing–review and editing. F. Zhao: Data curation, formal analy- sis, validation, investigation, visualization, methodology, writing– original draft, writing–review and editing. C. Chen: Data curation, ### Page 21 RESEARCH ARTICLESingle-Cell Survival AnalysisMAY 2026 CANCER DISCOVERY | 951writing–review and editing. L. Zhou: Data curation. L.-Y. Wu:Methodology, writing–review and editing. G.B. Mills: Investiga- tion, methodology, writing–review and editing. L.M. Coussens:Investigation, methodology, writing–review and editing. Z. Xia:Conceptualization, resources, data curation, supervision, funding acquisition, investigation, methodology, writing–original draft, project administration, writing–review and editing.AcknowledgmentsThis work was supported by the following funding sources: the Department of Defense Prostate Cancer Data Science Award HT94252410551, NIH R01GM147365 and R01CA283171, and the Silver Family Innovation Foundation Award (to Z. Xia); NIH U01CA253472, U01CA281902, and U24CA264128 (to G.B. Mills); and the Susan G. Komen Foundation, the National Foundation for Cancer Research, and Hildegard Lamfrom Endowed Chair in Basic Research (to L.M. Coussens). The research reported in this publica- tion used computational infrastructure supported by the Office of Research Infrastructure Programs, Office of the Director, of the NIH under Award Number S10OD034224. We acknowledge the The Cancer Genome Atlas Research Network for generating, curating, and providing access to the liver hepatocellular carcinoma dataset used in this study. We also acknowledge the use of publicly available scRNA-seq datasets deposited in the GEO and the Sequence Read Archive. We thank the investigators who generated and shared these valuable datasets. The content is solely the responsibility of the authors and does not necessarily represent the official views of the funders.noteSupplementary data for this article are available at Cancer Discovery Online (http://cancerdiscovery.aacrjournals.org/). Received June 3, 2025; revised October 8, 2025; accepted December 24, 2025; posted first April 21, 2026.RefeRenCeS1. Zhou J-G, Zhao H-T, Jin S-H, Tian X, Ma H. Identification of a RNA-seq-based signature to improve prognostics for uterine sarcoma. Gynecol Oncol 2019;155:499–507.2. Gu Y, Lu L, Wu L, Chen H, Zhu W, He Y. Identification of prognostic genes in kidney renal clear cell carcinoma by RNA-seq data analysis. Mol Med Rep 2017;15:1661–7.3. Raman P, Zimmerman S, Rathi KS, de Torrenté L, Sarmady M, Wu C, et al. A comparison of survival analysis methods for cancer gene expression RNA-sequencing data. Cancer Genet 2019;235–236:1–12.4. Huang Z, Johnson TS, Han Z, Helm B, Cao S, Zhang C, et al. Deep learning-based cancer survival prognosis from RNA-seq data: approaches and evaluations. BMC Med Genomics 2020;13(Suppl 5):41.5. Tang F, Barbacioru C, Wang Y, Nordman E, Lee C, Xu N, et al. mRNA-Seq whole-transcriptome analysis of a single cell. Nat Methods 2009;6:377–82.6. Tang F, Barbacioru C, Nordman E, Li B, Xu N, Bashkirov VI, et al. RNA-Seq analysis to capture the transcriptome landscape of a single cell. Nat Protoc 2010;5:516–35.7. Potter SS. Single-cell RNA sequencing for the study of development, physiology and disease. Nat Rev Nephrol 2018;14:479–92.8. Kolodziejczyk AA, Kim JK, Svensson V, Marioni JC, Teichmann SA. The technology and biology of single-cell RNA sequencing. Mol Cell 2015;58:610–20.9. Sun D, Guan X, Moran AE, Wu L-Y, Qian DZ, Schedin P, et al. Iden- tifying phenotype-associated subpopulations by integrating bulk and single-cell sequencing data. Nat Biotechnol 2022;40:527–38.10. Mizukoshi C, Kojima Y, Hayashi S, Abe K, Kasugai D, Shimamura T. scSurv: a deep generative model for single-cell survival analysis. Bioin- formatics 2026;42:btaf671.11. Johnson TS, Yu CY, Huang Z, Xu S, Wang T, Dong C, et al. Diagnostic Evidence GAuge of Single cells (DEGAS): a flexible deep transfer learning framework for prioritizing cells in relation to disease. Genome Med 2022;14:11.12. Ilse M, Tomczak J, Welling M. Attention-based deep multiple instance learning. In: Jennifer D, Andreas K, editors. Proceedings of the 35th Inter- national Conference on Machine Learning. Volume 80. Proceedings of Machine Learning Research. Long Beach (CA): PMLR; 2018. p. 2127–36.13. Maron O, Lozano-Pérez T. A framework for multiple-instance learn- ing. Advances in Neural Information Processing Systems. 1997. Available from: https://proceedings.neurips.cc/paper_files/paper/ 1997/file/82965d4ed8150294d4330ace00821d77-Paper.pdf14. Carbonneau M-A, Cheplygina V, Granger E, Gagnon G. Multiple instance learning: a survey of problem characteristics and applica- tions. Pattern Recognition 2018;77:329–53.15. Lopez R, Regier J, Cole MB, Jordan MI, Yosef N. Deep generative modeling for single-cell transcriptomics. Nat Methods 2018;15: 1053–8.16. Xu C, Lopez R, Mehlman E, Regier J, Jordan MI, Yosef N. Probabilis- tic harmonization and annotation of single-cell transcriptomics data with deep generative models. Mol Syst Biol 2021;17:e9620.17. Pierson E, Yau C. ZIFA: dimensionality reduction for zero-inflated single-cell gene expression analysis. Genome Biol 2015;16:241.18. Ghosh SK, Mukhopadhyay P, Lu J-CJ. Bayesian analysis of zero-inflated regression models. J Stat Plann Inference 2006;136:1360–75.19. Xu J, Sun X, Zhang Z, Zhao G, Lin J. Understanding and improving layer normalization. In: Proceedings of the 33rd International Con- ference on Neural Information Processing Systems. Red Hook (NY): Curran Associates Inc.; 2019. p. 394.20. Hu J, Shen L, Albanie S, Sun G, Wu E. Squeeze-and-excitation net- works. IEEE Trans Pattern Anal Mach Intell 2020;42:2011–23.21. Tukey JW. Use of many covariates in clinical trials. Int Stat Rev 1991;59:123–37.22. Yves G, Yoshua B. Entropy regularization. In: Chapelle O, Scholkopf B, Zien A, editors. Semi-Supervised Learning. Cambridge (MA): The MIT Press; 2006. p. 151–68.23. Zappia L, Phipson B, Oshlack A. Splatter: simulation of single-cell RNA sequencing data. Genome Biol 2017;18:174.24. Sade-Feldman M, Yizhak K, Bjorgaard SL, Ray JP, De Boer CG, Jenkins RW, et al. Defining T cell states associated with response to check- point immunotherapy in melanoma. Cell 2018;175:998–1013.e20.25. Bill R, Wirapati P, Messemaker M, Roh W, Zitti B, Duval F, et al. CXCL9:SPP1 macrophage polarity identifies a network of cellular programs that control human cancers. Science 2023;381:515–24.26. Lauss M, Donia M, Svane IM, Jönsson G. B cells and tertiary lymphoid structures: friends or foes in cancer immunotherapy? Clin Cancer Res 2022;28:1751–8.27. Zhang Y, Yu B, Ming W, Zhou X, Wang J, Chen D. SpaTopic: a statisti- cal learning framework for exploring tumor spatial architecture from spatially resolved transcriptomic data. Sci Adv 2024;10:eadp4942.28. Thrane K, Eriksson H, Maaskola J, Hansson J, Lundeberg J. Spatially resolved transcriptomics enables dissection of genetic heterogeneity in stage III cutaneous malignant melanoma. Cancer Res 2018;78: 5970–9.29. Mantovani A, Sozzani S, Locati M, Allavena P, Sica A. Macrophage polarization: tumor-associated macrophages as a paradigm for polarized M2 mononuclear phagocytes. Trends Immunol 2002;23: 549–55.30. Cheng S, Li Z, Gao R, Xing B, Gao Y, Yang Y, et al. A pan-cancer single-cell transcriptional atlas of tumor infiltrating myeloid cells. Cell 2021;184:792–809.e23.31. Zhang T, Zhao F, Hu Y, Wei J, Cui F, Lin Y, et al. Environmental monobutyl phthalate exposure promotes liver cancer via repro- grammed cholesterol metabolism and activation of the IRE1α-XBP1s pathway. Oncogene 2024;43:2355–70.32. Yang Z, Huo Y, Zhou S, Guo J, Ma X, Li T, et al. Cancer cell-intrinsic XBP1 drives immunosuppressive reprogramming of intratumoral myeloid cells by promoting cholesterol production. Cell Metab 2022; 34:2018–35.e8. ### Page 22 RESEARCH ARTICLEAACRJournals.orgRen et al.952 | CANCER DISCOVERY MAY 202633. Barras D, Ghisoni E, Chiffelle J, Orcurto A, Dagher J, Fahr N, et al. Response to tumor-infiltrating lymphocyte adoptive therapy is asso- ciated with preexisting CD8+ T-myeloid cell networks in melanoma. Sci Immunol 2024;9:eadg7995.34. Chen JH, Nieman LT, Spurrell M, Jorgji V, Elmelech L, Richieri P, et al. Human lung cancer harbors spatially organized stem-immunity hubs associated with response to immunotherapy. Nat Immunol 2024;25:644–58.35. Chu Y, Dai E, Li Y, Han G, Pei G, Ingram DR, et al. Pan-cancer T cell atlas links a cellular stress response state to immunotherapy resis- tance. Nat Med 2023;29:1550–62.36. Gide TN, Quek C, Menzies AM, Tasker AT, Shang P, Holst J, et al. Distinct immune cell populations define response to anti-PD-1 monotherapy and anti-PD-1/anti-CTLA-4 combined therapy. Cancer Cell 2019;35:238–55.e6.37. Pauken KE, Shahid O, Lagattuta KA, Mahuron KM, Luber JM, Lowe MM, et al. Single-cell analyses identify circulating anti-tumor CD8 T cells and markers for their enrichment. J Exp Med 2021; 218:e20200920.38. Li N, Kang Y, Wang L, Huff S, Tang R, Hui H, et al. ALKBH5 regulates anti-PD-1 therapy response by modulating lactate and suppressive immune cell accumulation in tumor microenvironment. Proc Natl Acad Sci U S A 2020;117:20159–70.39. Lu BY, Lucca LE, Lewis W, Wang J, Nogueira CV, Heer S, et al. Circulat- ing tumor-reactive KIR+CD8+ T cells suppress anti-tumor immunity in patients with melanoma. Nat Immunol 2025;26:82–91.40. Xue R, Zhang Q, Cao Q, Kong R, Xiang X, Liu H, et al. Liver tumour immune microenvironment subtypes and neutrophil heterogeneity. Nature 2022;612:141–7.41. Wang Y, Wang Q, Tao S, Li H, Zhang X, Xia Y, et al. Identification of SPP1+ macrophages in promoting cancer stemness via vitronectin and CCL15 signals crosstalk in liver cancer. Cancer Lett 2024;604:217199.42. Fan G, Xie T, Li L, Tang L, Han X, Shi Y. Single-cell and spatial analyses revealed the co-location of cancer stem cells and SPP1+ macrophage in hypoxic region that determines the poor prognosis in hepatocellular carcinoma. NPJ Precis Oncol 2024;8:75.43. Liberzon A, Birger C, Thorvaldsdóttir H, Ghandi M, Mesirov JP, Tamayo P. The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell Syst 2015;1:417–25.44. Barcena-Varela M, Lujambio A. The endless sources of hepatocellular carcinoma heterogeneity. Cancers (Basel) 2021;13:2621.45. Yang Y, Chen X, Pan J, Ning H, Zhang Y, Bo Y, et al. Pan-cancer single-cell dissection reveals phenotypically distinct B cell subtypes. Cell 2024;187:4790–811.e22.46. Villéga F, Fernandes A, Jézéquel J, Uyttersprot F, Benac N, Zenagui S, et al. Ketamine alleviates NMDA receptor hypofunction through syn- aptic trapping. Neuron 2024;112:3311–28.e9.47. Breslow NE. Analysis of survival data under the proportional hazards model. Int Stat Rev 1975;43:45–57.48. Cox DR. Regression models and life-tables. J R Stat Soc Ser B (Methodol) 1972;34:187–202.49. Hao Y, Hao S, Andersen-Nissen E, Mauck WM, Zheng S, Butler A, et al. Integrated analysis of multimodal single-cell data. Cell 2021;184: 3573–87.e29.50. Wolf FA, Angerer P, Theis FJ. SCANPY: large-scale single-cell gene expression data analysis. Genome Biol 2018;19:15.51. Balakrishnan V. All about the Dirac delta function (?) Resonance 2003;8:48–58.52. Higgins I, Matthey L, Pal A, Burgess C, Glorot X, Botvinick M, et al. beta-VAE: learning basic visual concepts with a constrained varia- tional framework. 2017 [cited 2024 Jan 12]. Available from: https:// openreview.net/forum?id=Sy2fzU9gl.53. Csiszár I. I-divergence geometry of probability distributions and min- imization problems. Ann Probab 1975;3:146–58.54. Kingma DP, Salimans T, Jozefowicz R, Chen X, Sutskever I, Welling M. Improved variational inference with inverse autoregressive flow. In: Lee D, Sugiyama M, Luxburg U, Guyon I, Garnett R, editors. Proceedings of the 30th International Conference on Neural Information Processing Systems. Barcelona (Spain): Curran Associates Inc.; 2016. p. 4743–51.55. Brown PJ, Griffin JE. Inference with normal-gamma prior distribu- tions in regression problems. Bayesian Anal 2010;5:171–88.56. Wold S, Esbensen K, Geladi P. Principal component analysis. Chemo- metrics Intell Lab Syst 1987;2:37–52.57. Vaswani A, Shazeer N, Parmar N, Uszkoreit J, Jones L, Gomez AN, et al. Attention is all you need. In: Guyon I, Von Luxburg U, Bengio S, Wallach H, Fergus R, Vishwanathan S, et al, editors. Proceedings of the 31st International Conference on Neural Information Processing Systems. Long Beach (CA): Curran Associates Inc.; 2017. p. 6000–10.58. Yu T, Kumar S, Gupta A, Levine S, Hausman K, Finn C. Gradient surgery for multi-task learning. Adv Neural Inf Process Syst 2020;33: 5824–36.59. Harrell FE Jr, Lee KL, Mark DB. Multivariable prognostic models: issues in developing models, evaluating assumptions and adequacy, and measuring and reducing errors. Stat Med 1996;15:361–87.60. Kaplan EL, Meier P. Nonparametric estimation from incomplete observations. J Am Stat Assoc 1958;53:457–81.61. Mantel N. Evaluation of survival data and two new rank order sta- tistics arising in its consideration. Cancer Chemother Rep 1966;50: 163–70.62. Hossin M, Sulaiman MN. A review on evaluation metrics for data classification evaluations. Int J Data Min Knowl Manag Process 2015;5:1.63. Saikia R, Barman MP. A review on accelerated failure time models. Int J Stat Syst 2017;12:311–22.64. Lee DD, Seung HS. Learning the parts of objects by non-negative matrix factorization. Nature 1999;401:788–91.65. Pölsterl S. Scikit-survival: a library for time-to-event analysis built on top of scikit-learn. J Machine Learn Res 2020;21:1–6.66. Traag VA, Waltman L, van Eck NJ. From Louvain to Leiden: guaran- teeing well-connected communities. Sci Rep 2019;9:5233.67. Zhao F, Zhang T, Wei J, Chen L, Liu Z, Jin Y, et al. Integrated single-cell transcriptomic analyses identify a novel lineage plasticity-related cancer cell type involved in prostate cancer progression. EBioMedicine 2024;109:105398.68. Lefebvre C, Rajbhandari P, Alvarez MJ, Bandaru P, Lim WK, Sato M, et al. A human B-cell interactome identifies MYB and FOXM1 as master regulators of proliferation in germinal centers. Mol Syst Biol 2010;6:377.69. Liu F, Zhang T, Yang Y, Wang K, Wei J, Shi JH, et al. Integrated anal- ysis of single-cell and bulk transcriptomics reveals cellular subtypes and molecular features associated with osteosarcoma prognosis. BMC Cancer 2025;25:280.70. Hugo W, Zaretsky JM, Sun L, Song C, Moreno BH, Hu-Lieskovan S, et al. Genomic and transcriptomic features of response to anti-PD-1 therapy in metastatic melanoma. Cell 2016;165:35–44.71. Łuksza M, Riaz N, Makarov V, Balachandran VP, Hellmann MD, Solovyov A, et al. A neoantigen fitness model predicts tumour response to checkpoint blockade immunotherapy. Nature 2017;551: 517–20.72. Wang YP, Yu GR, Lee MJ, Lee SY, Chu IS, Leem SH, et al. Lipocalin-2 negatively modulates the epithelial-to-mesenchymal transition in hepatocellular carcinoma through the epidermal growth factor (TGF-beta1)/Lcn2/Twist1 pathway. Hepatology 2013;58:1349–61.73. Roessler S, Jia HL, Budhu A, Forgues M, Ye QH, Lee JS, et al. A unique metastasis gene signature enables prediction of tumor relapse in early-stage hepatocellular carcinoma patients. Cancer Res 2010;70:10202–12.74. Cancer Genome Atlas Research Network. Electronic address wbe, cancer genome atlas research N. Comprehensive and integrative genomic characterization of hepatocellular carcinoma. Cell 2017;169: 1327–41.e23.75. Kudo A, Mogushi K, Takayama T, Matsumura S, Ban D, Irie T, et al. Mitochondrial metabolism in the noncancerous liver determine the occurrence of hepatocellular carcinoma: a prospective study. J Gastro- enterol 2014;49:502–10.76. Ning J, Ye Y, Shen H, Zhang R, Li H, Song T, et al. Macrophage-coated tumor cluster aggravates hepatoma invasion and immunotherapy resistance via generating local immune deprivation. Cell Rep Med 2024;5:101505.