Single-cell and spatial analysis of aging-related cancer-supporting cells across cancers reveals their molecular features and predicts outcomes in neuroblastoma using multi-omics and machine learning

Analyzing the bulk-seq data, we obtained five transcriptome datasets of NB: GSE49710↗ and GSE85047↗ sourced from GEO database, TARGET-NB from TARGET database, and E-MTAB-8248 and E-MTAB-179 from ArrayExpress database. We utilized the log2 (x+1) algorithm to normalize transcriptomic data, and conducted the combat function of the “sva” R package to solve batch effects (14). After excluding patients with incomplete follow-up data, we incorporated 1617 patients in total. We used the GSE49710↗ cohort as the training cohort to develop both prognostic and diagnostic signatures, contrastingly, the remaining four cohorts served as validation cohorts. Comprehensive patient information in bulk-seq cohorts was listed in Supplementary Table S1. Moreover, scRNA-seq datasets of NB (GSE137804↗, GSE192906↗ and GSE140819↗) were acquired from the GEO database, with detailed clinic information of NB patients in GSE137804↗ included in Supplementary Table S2. Another NB scRNA-seq cohort was download from https://www.neuroblastomacellatlas.org/↗. There is no requirement for ethical approval or patient consent, as the bioinformatics data in our analysis is publicly available. The general workflow of our study was illustrated in Supplementary Figure S1A. For pan-cancer single cell data, GSE176078↗, GSE203612↗, GSE138709↗, GSE142784↗, GSE166555↗, GSE181919↗, GSE149614↗, GSE131907↗, GSE184880↗, GSE215120↗ and GSE139829↗ datasets were sourced from the GEO database for analysis. The scRNA-seq data of PRAD from Chen et al. were downloaded from http://www.pradcellatlas.com/↗. Genes related to senescence were sourced from Molecular Signatures Database (MSigDB) (https://www.gsea-msigdb.org/gsea/msigdb/cards/FRIDMAN_SENESCENCE_UP.html↗). This search identified 77 senescence related genes, detailed in Supplementary Table S3.

For pan-cancer scRNA-seq data, raw gene expression matrices were loaded into a Seurat object and read into R software via Seurat R package (15). We removed low-quality cells according to criteria of >40,000 UMI per cell, <500 genes per cell, >5,000 genes per cell and >20% mitochondrial genes. For scRNA-seq data retrieved from inDrop platform, we filtered out cells with UMI counts >40,000, gene counts <200, gene counts >5000, and mitochondrial gene count >20%. Doublets were excluded by DoubletFinder R package (16). The “Harmony” R package was used to remove technical batch effects while keeping biological variance after batch consolidation (17). The local inverse Simpson’s Index (LISI) was utilized to evaluate corrections of batch effects (17). For analysis of NB, scRNA-seq data of GSE137804↗ was loaded as Seurat object via “Seurat” R package (15). We performed quality control to delete low-quality cells with feature counts above 7500 or below 300, or mitochondrial counts above 20%. The function FindVariableFeatures was used to identify the top 2000 genes with most cell variants (18), at which PCA was conducted. The FindClusters function was employed to identify different clusters at a resolution of 0.5. We performed the RunUMAP function to reduce dimensions and manually annotated major cell types based on canonical markers (5). With the makers of CAF subtypes previously reported (19, 20), various clusters of CAFs were identified. Obtaining markers of CAF subtypes was performed with the criteria of log2FC >0.25 and p value <0.05 by the FindAllMarkers function, illustrated in Supplementary Table S4. Seven single cell scoring algorithms (AUCell in “AUCell” R package, Ucell in “Ucell” R package, ssGSEA and GSVA in “GSVA” R package, singscore in “singscore” R package, AddModuleScore and PercentageFeatureSet in “Seurat” R package) were applied to perform enrichment scoring and pinpoint senescence related cell type (15, 21–23), visualized by the “irGSEA” R package (https://chuiqin.github.io/irGSEA/index.html↗). Functional enrichment analysis via GO and KEGG database in scRNA data was conducted by “ClusterGVis” R package (https://github.com/junjunlab/ClusterGVis↗). To explore the specific subgroup preferences of CAF subpopulations, we calculated the odds ratios (OR) via the computational method of Zhang et al. (24). A numeric vector, predicting cellular status from least (1.0) to most (0.0) differentiated was generated from the RNA matrix by “CytoTRACE” R package (25). The “slingshot” and “Monocle3” R package was utilized to infer cell lineages and pseudotime states (26, 27).

To explore biological landscapes of CAFs in pan-cancer spatial transcriptomics (ST) resolution, we sourced ST slide data of various cancer types form 10x database (https://www.10xgenomics.com/cn/datasets↗), as well as GSE176078↗, GSE179572↗, GSE203612↗ and GSE181300↗ from GEO database. To accurately acquire the cell composition at each spot on the 10x ST slide, deconvolution analysis was applied (28), which is based on ST and scRNA-seq data, with particular consideration given to the corresponding tumor type. We first obtained scRNA-seq data of various samples with the same tumor type in Tumor Immune Single-cell Hub 2 (TISCH2) (29), and then constructed a comprehensive scRNA-seq reference library. Ensuring the reliability of analysis, we applied strict quality control measures to the scRNA-seq data according to numbers of expressed genes, counts of UMI, and percentage of mitochondrial RNA. In terms of screening parameters, we refer to the relevant studies of sRNA-seq data sources to ensure the scientific and accurate screening criteria. We then constructed a signature score matrix via computing the average expressions of the top 25 specifically expressed genes of various cell types in each site’s scRNA reference. Subsequently, by get_enrichment_matrix and enrichment_analysis in the “Cottrazm” R package (30), we successfully generated the enrichment score matrix, which can provide a powerful support for subsequent cell composition analysis. The enrichment score for each cell type is visualized using the SpatialFeaturePlot function via “Seurat” R package. If the score of Malignant cells in the microregion is 1, the Malignant group is defined. If it is 0, the Normal group is defined; otherwise, the Mixed group is defined. Wilcoxon Rank Sum Test was utilized to assess the statistical difference in gene expression between the three predicted areas. Meanwhile, to obtain the spatial coordinates of CAFs, we performed combined analysis of the scRNA-seq data and ST data via CellTrek R package (31) with its default parameters. We used the run_kdist function to calculate the spatial k-distance between different CAF subtypes and cell subpopulations.

Utilizing marker genes of senes CAF with prognostic value, we established senes CAF related signature (SCRS) for diagnosing INSS stage 4 NB and predicting prognosis in NB. We unified 10 prognostic machine learning (ML) methods, involving random survival forest (RSF), elastic network (Enet), Lasso, Ridge, stepwise Cox, CoxBoost, partial least squares regression for Cox (plsRcox), supervised principal components (SuperPC), generalized boosted regression modeling (GBM), and survival support vector machine (survival-SVM) for prognosis prediction. Subsequently, 12 diagnostic ML procedures were applied, such as random forest (RF), Lasso, Ridge, elastic net (Enet), stepwise Glm, GlmBoost, LDA, partial least squares regression for Glm (plsRglm), GBM, XGB, SVM and Naives Bayes for the diagnosis of stage 4 NB among all NB patients. A total of 101 prognostic ML methods and 113 predictive ML algorithms were trained within the training cohort, via the leave-one-out cross-validation (LOOCV) framework to establish the prognostic and diagnostic models. Models with fewer than five genes were removed. GSE49710↗ cohort served as the training set, while GSE85047↗, TARGET-NB, E-MTAB 8248 and E-MTAB 179 cohorts served as testing sets. The concordance index (C-index) and the average area under the curve (AUC) was calculated in each ML combination across the five sets (13). AUC, Area under precision-recall curve (PRAUC), accuracy, sensitivity, specificity, precision, cross-entropy, Brier scores, balanced accuracy and F1 Score, calculated by “mlr3” R package, were used to select the best diagnostic model (32). We employed precision-recall curve (PRC) to evaluate the performance of classification models in handling imbalanced datasets. Logarithmic loss, recall and decision calibration were utilized to select the best prognostic model, among the top five prognostic models with the highest C-index, via “mlr3proba” R package (33). Risk scores were retrieved with a linear combination method for each prognostic ML combination by including gene expression data from multiple feature selection patterns. Likewise, the most powerful diagnostic model was utilized to calculate the possibility of stage 4 NB.

Comprehensive validation methods were carried out to validate the superior precision, stability, and replicability of SCRS. With the prognostic SCRS, patients were divided into high or low-risk groups by the median risk score from the training set. Validations was performed based on Kaplan-Meier (KM) survival methods and a log-rank test, with the “survival” and “survminer” R packages. For the diagnostic model, a confusion matrix was utilized to test accuracy of SCRS by “cvms” R package. Receiver Operating Characteristic (ROC) curves, calibration curves and decision curve analysis (DCA) were applied to assess the precision, differentiation, and clinic utility of both diagnostic and prognostic SCRS. Subsequently, we compared the predictive power of the prognostic signature with conventional clinic features in time-dependent ROC curves. Moreover, univariate and multivariate Cox analysis were used to confirm the independent predictive advantage of SCRS.

With altogether 16 model genes in diagnostic and prognostic SCRS, we performed unsupervised clustering analysis in three NB cohorts (GSE49710↗, E-MTAB 8248, and TARGET-NB). We used the “ConsensusClusterPlus” R package and the k-means algorithm, which involved 1,000 repetitions and sampled 80% of the total data in each instance, to explore the molecular landscapes of senescent myofibroblast makers (34). Principal components analysis (PCA) was performed to depict the heterogeneity between clusters. The effectiveness of clustering analysis was appraised by comparing differences in clinicopathological characteristics and gene expression levels between clusters, utilizing the “ComplexHeatmap” R package. Moreover, survival analysis was employed to illustrate survival outcome disparities between clusters.

Differentially expressed genes (DEGs) were identified from two clusters divided by consensus molecular clustering and from two risk groups categorized by prognostic SCRS. DEGs were discovered with the “limma” R package, by a False-discovery rate (FDR) threshold of <0.05 and an absolute log2fold change (FC) of >1. The function enrichment of DEGs was performed in Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) terms with the “clusterprofiler” R package (35). Gene set variation analysis (GSVA) was conducted by KEGG terms with “GSVA” R package (23), based on the “h.all.v7.4.symbols.gmt” gene set from MSigDB. Gene set enrichment analysis (GSEA) was applied to explore the molecular pathways associated with different clusters and risk groups (36), with a threshold of p < 0.05 and Normalized Enrichment Score (NES) > 1.

We utilized several immune calculation algorithms by the “IOBR” R package and the single sample gene set enrichment analysis (ssGSEA) to assess the immune infiltration levels between groups and clusters by “wilcox” test (37–45). Immune cell marker genes for ssGSEA were sourced from literature (46). Next, the correlation analysis was applied in the Spearman method to reveal the relation between SCRS risk scores, model gene expressions, and immune cell infiltrations. Further, utilizing immune function marker genes, we performed ssGSEA to compare the immune function level between risk groups with “wilcox” test (47). Subsequently, ssGSEA was employed to assess the seven steps of cancer immunity cycle, with marker genes in Tracking Tumor Immunophenotype (TIP) (http://biocc.hrbmu.edu.cn/TIP/↗) (48). Lastly, we appraised the gene expression levels of immune checkpoint genes in two risk groups. We conducted immune subtype analysis which could identify immune response subtypes and predict immunotherapy reactions (49). Five immunological patterns were revealed in GSE49710↗ cohort, covering wound healing (C1), IFN-gamma dominant (C2), inflammatory (C3), lymphocyte depleted (C4) and TGF-β dominant (C6). We then compared the proportions of immunological subtypes between groups and clusters.

Obtaining the somatic mutational data in cBioPortal database (https://www.cbioportal.org/↗), we appraised the mutational types and frequencies of SCRS model genes with “maftools” R package (50). Moreover, we analyzed the tumor mutation burden (TMB) via calculating the aggregate count of somatic mutations per megabase (MB) within the exonic region of the human genome. Gene mutations were stratified into two types namely synchronous or nonsynchronous mutation. The latter comprised Frame_Shift_Del, Frame_Shift_Ins, In_Frame_Del, In_Frame_Ins, Missense, Nonsense, Nonstop, Splice_Site, and Translation_Start_Site aberrations. Key mutation regions were identified from the copy number variation (CNV) data in cBioPortal database by GISTIC 2.0 (51). Besides, the gene frequency of somatic CNV was visualized in “bubble plot”, while the chromosomal location of gene mutation was visualized in “circle plot” via “RCircos” R package (52).

To assess effects of SCRS to forecast immunotherapy response, we compared immune dysfunction and exclusion (TIDE, http://tide.dfci.harvard.edu/↗) scores between risk groups. Subsequently, we utilized submap algorithm to evaluate immunotherapy responsiveness of patients based on an immunotherapy cohort (53, 54). Meanwhile, we assessed the power of SCRS to predict responses to immunotherapy in immunotherapy datasets (IMvigor210, GSE78220↗, GSE135222↗, and GSE91061↗). Lastly, we acquired the sensitivity to chemotherapy drugs of human cancer cell lines in Cancer Therapeutics Response Portal (CTRP, https://portals.broadinstitute.org/ctrp↗) database and Profiling Relative Inhibition Simultaneously in Mixtures (PRISM, https://depmap.org/portal/prism/↗) database. Cell line which is more responsive to a chemotherapy agent could get a lower AUC, which help explore potential therapeutics for high-risk patients (55).

Six scRNA scoring algorithms were utilized to assess the specific enriched cells of SCRS model genes. We utilized the ‘Scissor’ R package to pinpoint the particular cell populations responsible for SCRS status variations (56), which capitalizes on both aggregate data and phenotypic data, enabling the autonomous selection of cell subpopulations from single-cell datasets, attributing to distinct phenotypes. Pseudotime trajectory analysis was conducted by “Monocle” and “Monocle3” R package, depicting maps of development trajectories devoid of preexisting knowledge about differentiation commencement (27, 57). We used “InferCNV” R package to determine CNVs of neuroendocrine (NE) cells, Schwann cells, endothelial cells, and fibroblasts, referred to T cells, B cells, monocytes and macrophages (58). “CellChat” and “NicheNet” R package was utilized to explore intercellular communication networks in high-SCRS and low-SCRS cells, respectively (59, 60). Subsequently, we utilized “pySCENIC” (version 0.11.2) with Python (version 3.7) to explore enriched transcription factor (TF) and regulon activities of each cell type, which build TF regulatory network and identify steady cell state (61).

We performed pan-cancer analysis by “TCGAplot” R package to reveal the commonalities and disparities in genomic and cellular modifications of SCRS hub gene across various cancer types, focusing on gene expression, tumor mutation burden (TMB), microsatellite instability (MSI), immunological microenvironment, and prognosis value (62). We compared gene expression in tumor and normal samples by “wilcox” test and paired samples t-test. Spearman correlation analysis was performed to reveal associations between expression of hub gene and immune infiltrations. Immune cell ratio data was obtained in The Immune Landscape of Cancer (https://api.gdc.cancer.gov/data/b3df502e-3594-46ef-9f94-d041a20a0b9a↗), and immune score was calculated by ESTIMATE method.

To verify the different expressions of JAK1 between stage 4 and other NB samples, as well as its prognostic value, we performed immunohistochemistry (IHC) staining in tissue samples from 24 tissue specimens of stage 4 NB and 16 NB tissue specimens of other stages. The study was approved by the ethics committee of Children’s Hospital of Chongqing Medical University. NB tissues were paraffin-embedded and sectioned into 4 mm slices. After dewaxing, hydration, and antigen retrieval, the samples were incubated overnight at 4°C with primary antibody: Anti-JAK1 (Proteintech, Cat No: 66466-1-Ig). Subsequent steps included incubation with a Goat anti-Rabbit IgG secondary antibody (ZENBIO, China), DAB staining (ZENBIO, China), and blocking. Staining was observed under a microscope. Each sample was evaluated for staining intensity (0: none, 1: mild, 2: moderate, 3: strong) and the percentage of positive cells (0: 0%, 1: 1–25%, 2: 26–50%, 3: 51–75%, 4: 76–100%). The final IHC score was the sum of intensity and percentage scores.

To construct a pan-cancer single-cell and spatial transcriptomics landscape, we downloaded pan-cancer scRNA-seq data from GEO database of 13 prevalent cancer types, along with ST data from 10x website and GEO database. The cancer types of scRNA-seq data that we retrieved from database involved: breast cancer (BRCA), colorectal cancer (CRC), liver hepatocellular carcinoma (LIHC), lung adenocarcinoma (LUAD), cholangiocarcinoma (CHOL), ovarian cancer (OVCA), prostate adenocarcinoma (PRAD), head and neck squamous cell carcinoma (HNSC), neuroblastoma (NB), skin cutaneous melanoma (SKCM), stomach adenocarcinoma (STAD), uveal Melanoma (UVM) and uterine corpus endometrial carcinoma (UCEC) (Supplementary Figure S1B). To diminish the batch effects among different scRNA-seq datasets, we independently analyzed each dataset and annotated cells with canonical markers of major cell types. We extracted CAFs from each scRNA-seq dataset and merged them into a Seurat object to establish an integrated scRNA-seq dataset of pan-cancer CAFs. After quality control (Supplementary Figure S1C) and batch effect correction by Harmony (Supplementary Figure S1D), a total of 35,048 cells in the pan-cancer scRNA-seq data were remained for the following analysis. We next explored the heterogeneity of CAFs in pan-cancer landscape and identify seven CAF subpopulations, as well as pericytes and smooth muscle cells (SMCs) (Figure 1A). To evaluate Harmony for batch effect removement, we calculated local inverse Simpson’s Index (LISI) of four batch correction methods, indicating well batch correction after Harmony, as well as the superior ability of batch correction in Harmony algorithm (Supplementary Figure S1E).

In pan-cancer CAF atlas (Figure 1A), we annotated the CAF subpopulation which highly expressed chemokines (CCL19, CCL21 and CXCL2, Figure 1B) as inflammatory CAF (iCAF), likely to the previously reported iCAFs (63). A CAF subpopulation which highly expressed extracellular matrix (ECM) remodeling genes (MMP11, CTHRC1 and COL1A2, Figure 1B) was annotated as matrix CAF (mCAF), resemble the previously reported mCAFs in cancer data (64). To reveal the molecular landscape of senescence in CAF, we utilized six scRNA scoring algorithms and senescence geneset to conduct enrichment scoring, investigating senescence related CAF subset with the highest enrichment scores of senescence-related genes (Figures 1D, E). Finally, we discovered senescence related CAF (senes CAF) subset which ranked first in senescence enrichments scores among all CAF subpopulations (Figure 1E; Supplementary Figure S1F), as well as highly expressed senescence related genes retrieved from MSigDB database (OPTN, RAB31 and IFI16, Figure 1F). Moreover, we found a CAF subset highly expressed MHC-II-associated antigen presentation genes (HLA-DRA, HLA-DRB1 and CD74, Figure 1F), which we annotated it as antigen presenting CAF (ap CAF) (65). Then, we observed a CAF subpopulation which was associated with glycolytic process and ATP generation from ADP in GO BP terms (Figure 1C), annotating it as metabolic CAF (meta CAF) which was similar to previous research (66). Subsequently, a CAF subpopulation which highly expressed cycle-related genes (CENPF, NUSAP1, and PTTG1, Figure 1F) was annotated as proliferative CAF (prolif CAF), consistent with previous pan-cancer research (67). Besides, we annotated progenitor like CAF (pro CAF) subtype with marker genes of IGF1, OGN and GSN (Figure 1F) (68). We utilized top marker genes of each CAF subtype to perform GO (Figure 1C) and KEGG (Supplementary Figure S1G) functional enrichment analysis, which further verified our CAF annotation results. In pan-cancer CAF landscape, we found different CAF subtypes displayed significantly different preferences of cancer types (Figures 1G, H). We found that iCAFs and mCAFs were most abundant in CAF subpopulations, while iCAFs were preferred in LUAD, whereas mCAFs were preferred in OVCA (Figures 1G, H). And we observed that senes CAFs, were preferred in UCEC and HNSC (Figure 1H). The diverse cellular characteristics of CAF subtypes stem from their varied origins, including both tissue-resident fibroblasts and pericytes (69). We utilized CytoTRACE analysis to discover that prolif CAFs owned highest CytoTRACE scores, inferring a possibility to occur in the earlier state of prolif CAFs (Figure 1I). With prolif CAFs set as the start point, we used Monocle3 (Figure 1J) and slingshot (Figure 1K) respectively to conclude differentiation lineages and pseudotime scores among CAF subtypes, revealing the complexity and heterogeneity of CAF differentiation. Meanwhile, we explored the differentially essential motifs within the subtypes of CAFs through SCENIC analysis (Figure 1L), with TCF12 motif associated with ECM remodeling ranked highly in mCAFs (70), as well as interferon regulatory factors (IRF) family associated with inflammaging enriched in senes CAFs (71).

CellTrek is a computational tool designed to map individual cells directly to the corresponding spatial positions in tissue sections by integrating scRNA-seq and ST data (31). Unlike traditional ST deconvolution techniques, this method transfers ST coordinates to single cells, enabling resolution at the single-cell level. We utilized CellTrek on high-quality scRNA-seq and ST datasets across various cancer types to reconstruct spatial single-cell atlases, involving CRC and CRC liver metastasis (72), as well as BRCA (73), BRCA brain metastasis (74), OVCA (73) and LIHC (75) (Figures 2A–H). Remarkably, even in the absence of corresponding scRNA-seq data from the same patient, the ST datasets were largely represented by scRNA-seq data through co-embedding analysis (Supplementary Figure S1H). From left to right, the first figure depicts the spatial distribution of various cell types by RCTD deconvolution analysis in tumor tissue section (Figures 2A–H). The second figure displays the spatial characteristics of malignant areas, mixed areas and normal areas divided by Cottrazm analysis, indicating the spatial distribution of tumor border and tumor immune barrier in tumor tissue section (Figures 2A–H). The third figure reveals the spatial positions of CAF subtypes in tissue sections by CellTrek mapping analysis, which demonstrates an abundant infiltration of various CAF subpopulations and spatial existence of senes CAFs in tumor microenvironment (Figures 2A–H). For the fourth figure, we calculated the spatial k-distance between all major cell types and all CAF subtypes in every tumor tissue section, suggesting that CAF subpopulations exhibit the minimum spatial k-distance among themselves (Figures 2A–H). For major cell types, we observed that senes CAF exhibit close spatial k-distance to endothelial cells, especially in CRC liver metastasis (Figure 2B) and OVCA (Figure 2E), while senes CAF also show close spatial k-distance to B cells in CRC (Figure 2A) and LIHC (Figure 2G).

In the GSE137804↗ scRNA cohort, we finally retrieved 172,564 cells after proper quality control (Supplementary Figure S2A). We performed “harmony” algorithm for batch effect removement, with a fine correction before and after integration (Supplementary Figure S2B). Based on cell markers in literature, we manually annotated 10 major cellular subtypes, including NE cells, T cells, B cells, NK cells, endothelial cells, macrophages, monocytes, Schwann cells, fibroblasts, and plasmacytoid DC cells (Figure 3A), with subpopulation of macrophages, fibroblasts, B cells and T cells illustrated respectively (Figure 3A). Figure 3B depicted top 5 marker genes of every major cell type. We visualized top marker genes of each major cell type in heatmap (Figure 3C), and utilized GO and KEGG terms to functionally annotated each major cell type, which further verified our single cell annotation (Figure 3C). We utilized “Dimplot” to visualize expression levels of several marker genes (CD79A for B cells, PHOX2B for NE cells, COL1A1 for fibroblasts, CD7 for T cells, SPP1 for macrophages and LYZ for monocytes, Figure 3D). Figure 3E demonstrated the different cell proportions in NB patients of INSS stage 1, 3, 4 and 4S, respectively, with NE cells predominating. Figure 3F depicted UMAP view of each cell subsets (top) and cell density (bottom) showing cell distribution across four stages. Downsampling was applied for four tissue groups. High relative cell density is shown as bright magma.

We subset 3105 fibroblasts and created a new Seurat object for the following scRNA analysis. Supplementary Figure S2C showed a well quality control of fibroblast subtypes and Supplementary Figure S2D displayed a fine correction of batch effects with Harmony. Based on CAF markers retrieved from literature review (20, 68), we manually annotated 10 CAF subtypes in NB single cell data, namely myofibroblasts (ACTA2, TAGLN, MYH11), proliferative CAF (MKI67), matrix CAF (POSTN, COL3A1), inflammatory CAF (CFD, C3), tumoral CAF (TMEM158), heat-shock protein CAF (HSP90AA1), reticular CAF (CCL19, CCL21), antigen-presenting CAF (CD74, HLA-DRA, HLA-DRB1), vascular CAF (NOTCH3), and interferon CAF (IL32) (Figure 4A; Supplementary Figure S2E). Figure 4B depicted top 5 marker genes of every major cell type. Figure 4C illustrated different compositions of CAF subtypes in different INSS stages. Previous literature has reported about senescent myofibroblasts that localize near tumor ducts and accumulate with tumor progression (76), so we utilized seven scRNA scoring algorithms to conduct enrichment scoring and discover senescence related CAF subset, which highlighted myofibroblasts closely associated with senescence (Figure 4D). Ranking results of single cell scoring confirmed that myofibroblasts were enriched in senescence and could possibly be identified as the senes CAFs (Figure 4E; Supplementary Figure S2F). We have included additional three independent scRNA-seq datasets of NB in our analysis to discover the consistent findings (Supplementary Figures S2G, H). We visualized top marker genes of each CAF subtype in heatmap (Figure 4F), and utilized GO and KEGG terms to functionally annotated each cell type, which further verified our annotation (Figure 4F). Interestingly, marker genes of iCAFs were enriched in aging, shedding light on its potential associations with senescence. With the FindAllMarkers function, we identified 1088 significant senescence related CAF markers with a criterion of absolute log2 (fold change) > 0.25 and p-value (Padj) < 0.05 (Supplementary Table S4). To sum up, we succeeded in identifying senes CAF subtype, namely myofibroblasts, which is significantly related to senescence. In CAF landscape of NB, different CAF subtypes exhibited different preferences of INSS stages (Figure 4G). Subsequently, we utilized CytoTRACE and pseudotime analyses to infer cellular differentiation states and dynamics of lineage specification. Heat-shock protein CAF owned highest CytoTRACE scores, with a likelihood to occur in the earlier state (Figure 4H). Setting heat-shock protein CAF as the beginning, we then utilized Monocle3 (Figure 4I) and slingshot (Figure 4J) respectively to infer several differentiation lineages among CAF subtypes, which showed a relatively senescent state of senes CAF. Moreover, we analyzed the differentially essential motifs within the subpopulations of CAFs with SCENIC analysis (Supplementary Figure S2I), with motif SRF, highly associated with cellular senescence, expressed in senes CAFs (77).

We performed batch removements across five NB cohorts, with PCA plots visualizing the well behaviors of batch correction (Supplementary Figure S3A). Marker genes of senes CAF were used for feature selection and model development. Ultimately, 13 marker genes with diagnosis significance were involved in the development of diagnostic ML model, while 8 marker genes with prognosis significance were used for establishment of prognostic ML model (Supplementary Figure S3B). We performed 101 prognosis ML combinations and 113 predictive ML combinations based on LOOCV framework, to select the best ML combination (Supplementary Figure S3C). The model with the highest AUC was established with RF in feature selection and model development, scoring the top average AUC (0.862) in five cohorts (Figure 5A). AUC of every ML combination was computed in every test cohort (Supplementary Table S5). We utilized Area under precision-recall curve (PRAUC), accuracy, sensitivity, specificity, precision, cross-entropy, Brier scores, balanced accuracy, F1 Score and precision-recall curve (PRC) to reveal that RF model is the best diagnostic model with superior performance (Supplementary Figures S3D, E). In GSE49710↗ cohort (Figure 5B) and E-MTAB 8248 cohort (Figure 5C), confusion matrix exhibited a fine precision of SCRS with RF model. ROC curves in five datasets displayed a fine discrimination of SCRS (Figure 5D). We calculated AUC of SCRS, clinic variable and logistic regression (LR) model including clinic variables and SCRS, showing that SCRS and LR model outperformed (Figure 5E). Calibration curves in five datasets displayed well alignments between SCRS predicted possibility and observed possibility (Figure 5F), indicates that the model’s probability estimates are reliable and well-calibrated, as it ensures that the risk estimates provided by the model can be trusted to reflect the true likelihood of patient outcomes. DCA curves revealed the great clinical benefit of SCRS and LR model, outperforming other clinic factors (Figure 5G). We visualized the diagnostic LR model with nomogram to help clinical decision-making (Figure 5H). We obtained feature importance of 11 model genes finally selected by RF, which revealed that hub gene EPN2 was the most powerful in RF algorithm (Figure 5I).

The ML combination with the highest C-index was constructed by SuperPC in feature selection and model development, obtaining the highest average C-index (0.763) in five cohorts (Figure 6A). C-index of every ML combination was computed in every test cohort (Supplementary Table S6). Logarithmic loss, recall and decision calibration were computed to reveal that SuperPC model had the best performance in calibration and precision (Supplementary Figure S3F). In GSE85047↗ dataset, the low-risk patients had a longer overall survival (OS) and progression-free survival (PFS) than the high-risk patients (Figures 6B, C). ROC curves in 1-, 3- and 5-year OS displayed fine specificity in SCRS (Figure 6D). AUC values of 3-year OS indicated that SCRS and cox regression model including SCRS and clinic factors were more discriminative to predict prognosis than other clinic factors (Figure 6E). Time dependent ROC curves revealed that SCRS and cox regression model behaved better than common clinic factors in discriminative ability (Figure 6F). Calibration curves (Figure 6G) and DCA curves (Figure 6H) demonstrated that SCRS is powerful in accuracy and clinical benefit, which implies that using the SuperPC model to guide clinical decision-making would result in more effective identification of patients. Multivariate Cox regression analysis indicated that SCRS risk score was an independent prognosis variable in GSE85047↗ cohort (P < 0.001) (Figure 6I). We utilized univariate cox regression to depict the prognosis value of each model gene (Figure 6J). We visualized the feature importance of 8 model genes chosen by SuperPC, with CKS2 being the most influential (Figure 6K). These metrics of model evaluation consistently proved that SCRS exhibited precision and robustness in model performance.

Further explorations were undertaken to elucidate the underlying mechanism of SCRS model genes. PCA showed significant division between two risk groups, split by prognostic SCRS (Figure 7A). Utilizing DEGs sourced from differential analysis between two risk groups, we employed function enrichment analysis through GO and KEGG databases, which revealed that DEGs were enriched in outer kinetochore in GO terms, and were enriched in DNA replication in KEGG terms (Figures 7B, C). GSEA discovered that Ribosome and Motor proteins were elevated in high-risk patients, and cell adhesion molecules were decreased in high-risk patients (Figures 7D, E). Then, GSVA with”h.all.v7.4.symbols.gmt” gene set in MSigDB website demonstrated that high-risk patients were elevated in myc_targets_v1, and decreased in hedgehog_signaling (Figure 7F). The diverse expressions of SCRS model genes and the variation of clinic factors between two risk groups were obvious (Figure 7G). The hub gene JAK1 expressed significantly lower in high-risk patients, exhibiting significant protective prognosis value. Spearman correlation analysis indicated tight relations (correlation p value < 0.001) across SCRS model genes (Figure 7H). Utilizing CNV data to plot, we displayed that INPP4B had the highest somatic CNV frequency in diagnosis model genes, and VGLL4 had the highest somatic CNV frequency in prognosis model genes (Figures 7I, J). Moreover, we portrayed the loci of mutations within SCRS model genes on the chromosome (Figure 7K).

To assess the discrimination of prognostic SCRS in immune infiltration, we analyzed the immune cell abundance in two groups using eight immune algorithms. We used “ComplexHeatmap” R package to revealed the significantly lower immune cell infiltrations in high-risk patients (Figure 8A). Meanwhile, the spearman correlation method showed associations between immune cell abundances and SCRS risk scores, as well as relationships between immune cell abundances and model gene expressions (Figure 8B). Additionally, low-risk patients had higher expression profiles of immune checkpoint genes, prompting a sensitivity to immunotherapy (Figure 8C). Moreover, ssGSEA based on immunological function signatures revealed that the low-risk group was significantly more infiltrated in immune cells (Figure 8D). Besides, ssGSEA results of six key steps in cancer immunity cycle was significantly higher in low-risk patients (Figure 8E). The spearman correlation method showed inverse associations between immune function levels and SCRS risk scores (Figure 8F). We utilized top 10 maker genes of myofibroblasts and conducted ssGSEA to obtain immune abundance of senescence related myofibroblasts. Then we conducted spearman correlation analysis to reveal tight relations (correlation p value < 0.001) across immune cells and CAFs in TME, as well as their cox p value (Figure 8G). Given the well-established diverse metabolic inclinations and dependencies (78), we acquired various metabolic pathways in KEGG terms to explore associations between risk scores and cancer metabolic pathways, which revealed underlying metabolic functions of SCRS (Figure 8H).

We demonstrated the landscapes of somatic mutation in high-risk patients (Figure 9A) and low-risk patients (Figure 9B) with mutational information in TARGET cohort. High-risk patients owned more tumor mutation burden (TMB) than low-risk patients with no significance (Figure 9C). Meanwhile, TMB was positively related to SCRS risk scores according to spearman correlation analysis with no significance (Figure 9D). Then, we classified patients of TARGET cohort into high TMB and low TMB subgroup based on median TMB. Survival analysis revealed that high-risk patients with low TMB showed the shortest OS and EFS, and low-risk patients with high TMB showed the longest EFS, with significance (Figures 9E, F).

Utilizing TIDE and submap method, the responsiveness of immunotherapy was assessed in high-risk and low-risk patients. In GSE49710↗ dataset, high-risk group showed higher TIDE scores, higher TIDE dysfunction and exclusion scores, which is prone to exhibit immune escape during immunotherapy (Figure 9G). In IMvigor210 dataset, low-risk patients responsive to immunotherapy showed the longest OS, while high-risk patients not responsive to immunotherapy owned the shortest OS, with significance (Figure 9H). Meanwhile, in E-MTAB 8248 dataset, low-risk patients were more likely to respond to immunotherapy (Figure 9I). Subsequently, submap method revealed that low-risk group were more responsive to CTLA4 inhibitors (p = 0.007) and not responsive to PD1 inhibitors (Figure 9J). Subsequently, we compared therapy effects in four immunotherapy datasets between risk groups, which revealed that patients responsive to immunotherapy had lower risk scores, with significance (Figure 9K). Ultimately, to discover novel chemotherapy agents for high-risk patients identified by SCRS, we forecasted the medicine reaction utilizing drug sensitivity data sourced in CTRP and PRISM. By cross correlating the two pharmacogenomics datasets, we triumphantly obtained five possible medicines or compounds (BI-2536, GSK461364, SB-743921, ispinesib and talazoparib), which exhibited therapeutic effectiveness in high-risk patients (Figure 9L).

To deeply explore the expression profiles of SCRS model genes, GSE49710↗, E-MTAB 8248 and TARGET-NB were involved to perform consensus molecular clustering. We used SCRS model genes to perform unsupervised clustering analysis in every dataset, which revealed k = 2 with outstanding discrimination (Figure 10A). Meanwhile, PCA displayed indispensable disparities between two clusters (Figure 10B). Then, survival analysis showed that the cluster 2 had shorter OS (Figure 10C) and EFS (Figure 10D), with significance. Next, the expression landscapes of SCRS model genes and the clinic variables between two clusters were significantly different (Figure 10E). Using DEGs obtained between two clusters, we employed function enrichment analysis to found that DEGs were enriched in negative T cell selection in GO, and were enriched in ABC transporters in KEGG (Figures 10F, G). GSEA demonstrated that biosynthesis of cofactors was elevated in cluster 2, and NK-kappa B signaling pathway were decreased in cluster 2 (Figures 10H, I). GSVA, with “h.all.v7.4.symbols.gmt” gene set in MSigDB, demonstrated that cluster 2 was elevated in myc_targets_v2, and decreased in hedgehog_signaling (Figure 10J). Subsequently, we utilized eight immune infiltration methods to appraise the immunological infiltration variations in two clusters, as well as depicting Cox P value of every cell type (Figure 10K).

To validate the superior prognostic prediction ability of SCRS, we gathered gene coefficients of 39 released public NB prognosis signatures. Then, we compared C-index of every prognostic signature with SCRS in five NB datasets. Ultimately, we revealed that SCRS outperformed most of previous signatures in five datasets in prediction performances (Figure 11A), which qualified SCRS as a meaningful NB prognostic signature. Moreover, we demonstrated the relationships between risk groups, clusters and clinic factors via “sankey plot” (Figure 11B). After defining immunological subtypes of patients in GSE49710↗, we found significant differences of respective proportions of immunological subtypes between risk groups and clusters, indicating more wound healing (C1) subtype in high-risk patients and cluster 2 (Figure 11C). Comparisons of clinical factors between two risk groups revealed that low-risk patients showed longer prognoses and better clinic status in GSE49710↗ cohort (Figure 11D).

Verifying the risk-stratify ability of SCRS in single cell landscape, we calculated the SCRS enrichment scores via scRNA scoring algorithms, which revealed that higher-SCRS cells were mostly abundant in NE cells (Figures 12A, B). To elucidate the cellular origins supporting high-SCRS clinical manifestations, we employed the “scissors” R package to make correlations between bulk-seq data and scRNA data, which independently singles out cells having remarkable alignment with the desired phenotype. We labeled high-SCRS and low-SCRS states of patients as our foremost phenotypes, thus enabling the thorough detection of 2120 high-SCRS cells (scissors+ cells) and 1875 low-SCRS cells (scissors- cells) (Figures 12C, D). The SCRS values of scissors+ cells were significantly higher than scissors- cells and other cell types (Figures 12D, E), which indicated successful discovery of scissors+ cells representing diverse SCRS status. Subsequently, we performed pseudotime trajectory analysis via Monocle 2 algorithm to explore the temporal sequences of cellular differentiation in NE cells, fibroblasts, Schwann cells and endothelial cells, with part of NE cells being poorly differentiated (Figures 12F, G). We utilized SCRS scores calculated by singscore and Ucell algorithm to divided cells into high-SCRS and low-SCRS according to median SCRS score. Comparisons of pseudotime scores revealed that high-SCRS NE cells had an earlier differentiation trajectory than low-SCRS NE cells, which indicated that immature NE cells scored higher SCRS points (Figures 12H–J). Additionally, setting NE cells as differentiation starting point in Monocle 3 algorithm, we verified that immature NE cells scored higher SCRS points, which could serve as malignant cells (Figure 12K).

Validating the power of SCRS in single cell level, we used inferCNV analysis to explore the clone structures of four cell types mentioned above, which indicated that NE cells with chromosome 17q gain were probably malignant cells, while high-SCRS cells owned more CNVs (Figures 13A, B). In cell chat analysis, we utilized “circle plot” to display the communication frequencies and communication strengths between every cell type, as well as the integrative cell chatting networks in high-SCRS cells and low-SCRS cells (Figure 13C). We visualized the cell chat landscapes in high-SCRS cells, and demonstrated over-expressed ligand-receptor pairs and communication profiles between B cells, myeloid cells, Schwann cells and fibroblasts in high-SCRS cells (Figures 13D, E). Distinct cell types would generate dissimilar contributive cues affecting the total, inbound and outbound signals among high-SCRS and low-SCRS cells, with macrophages, monocytes, Schwann cells, fibroblast, and endothelial cells noting exceptional significance (Figures 13F, G; Supplementary Figure S3G). The cellular communication between CAF subtypes holds significant biological importance, we then explore the communication network of CAF subpopulations in NB and pan-cancer landscape (Supplementary Figure S3H). We calculated the correlations between the hub gene of senes CAFs and cell-cell communication strength, finding a significant association between JAK1 and communication strength of CAF subtypes in NB, which offered a research value of JAK1 in modulating cellular communication (Supplementary Figure S3I). Meanwhile, we assessed the association among regulon (TFs and target genes) activities and cell types in high-SCRS and low-SCRS cells via SCENIC, revealing that regulons of SMARCA4 and PBX3 were more active in high SCRS cells (Figure 13H), and regulons of VEZF1 and PDLIM5 were more active in low SCRS cells (Figure 13I).

Interestingly, we revealed that hub gene JAK1, with favorable prognosis value in NB, was abundant in oncological research value, which has been proved as oncogene in breast cancer and non-small cell lung cancer (79). Therefore, we conducted pan-cancer analysis to explore the heterogeneity of JAK1expression in tumor and normal samples in 33 cancer types (Figure 14A). Besides, the links between JAK1 and TMB, MSI, immune cells, and immunological scores underscored importance of JAK1 in TME, immune cell invasions, and immune checkpoints (Figures 14B–F). The protective prognostic value of JAK1 was seen in KIRC, similar to NB. The risky prognostic value of JAK1was seen in LUSC (Figure 14G). Furthermore, we conducted pan-cancer spatial transcriptomics analysis to comprehensively explore the expressions of JAK1 in malignant cells and malignant spots. In tumor types of BRCA, CRC, HNSC, CESC, GIST, KIRC, LUSC and OVCA (Figures 14H, I; Supplementary Figures S4A–F), hub gene JAK1 solidly exhibited positive correlations with the abundance of malignant cells across eight cancer types, while expressing higher in malignant areas than in normal areas in multiple spatial transcriptomics slides, indicating the underlying oncogenic role of JAK1. Interestingly, hub gene JAK1 abnormally exhibited negative correlations with the abundance of malignant cells in LIHC spatial transcriptomics slide, while expressing lower in malignant areas than in normal areas (Figure 14J). This result suggests a unique functional role for JAK1 in this particular cancer type, which highlights the complexity of cancer biology and the importance of context in gene function.

Our bioinformatics analysis has revealed that three hub genes, which were both diagnostic model genes as well as prognostic model genes in SCRS, were expressed significantly lower in stage 4 NB (Figure 15A) and were protective prognosis genes (Figure 15B) in GSE49710↗ cohort. To validate the elevated protein levels of JAK1 in other stages NB tissues compared to stage 4 NB tissues, as well as its protective prognosis value, we conducted IHC staining to corroborate our bioinformatics findings (Figure 15C). Our scoring results showed that protein levels of JAK1 was significantly higher in other stages NB tissues than stage 4 NB tissues (Figure 15D). Subsequently, we categorized 40 NB patients into high and low expression groups based on their median IHC scores of JAK1. Survival analysis indicated that patients with high JAK1 expression had better OS than those with low JAK1 levels, without significance (p = 0.15, Figure 15E).