Abstract

   This study tackles the persistent prognostic and management challenges
   of clear cell renal cell carcinoma (ccRCC), despite advancements in
   multimodal therapies. Focusing on anoikis, a critical form of
   programmed cell death in tumor progression and metastasis, we
   investigated its resistance in cancer evolution. Using single-cell RNA
   sequencing from seven ccRCC patients, we assessed the impact of
   anoikis-related genes (ARGs) and identified differentially expressed
   genes (DEGs) in Anoikis-related epithelial subclusters (ARESs).
   Additionally, six ccRCC RNA microarray datasets from the GEO database
   were analyzed for robust DEGs. A novel risk prognostic model was
   developed through LASSO and multivariate Cox regression, validated
   using BEST, ULCAN, and RT-PCR. The study included functional
   enrichment, immune infiltration analysis in the tumor microenvironment
   (TME), and drug sensitivity assessments, leading to a predictive
   nomogram integrating clinical parameters. Results highlighted dynamic
   ARG expression patterns and enhanced intercellular interactions in
   ARESs, with significant KEGG pathway enrichment in MYC + Epithelial
   subclusters indicating enhanced anoikis resistance. Additionally, all
   ARESs were identified in the spatial context, and their locational
   relationships were explored. Three key prognostic genes—TIMP1, PECAM1,
   and CDKN1A—were identified, with the high-risk group showing greater
   immune infiltration and anoikis resistance, linked to poorer prognosis.
   This study offers a novel ccRCC risk signature, providing innovative
   approaches for patient management, prognosis, and personalized
   treatment.

   Subject terms: Cancer, Computational biology and bioinformatics,
   Biomarkers, Medical research, Oncology, Urology

Introduction

   Renal cell carcinoma (RCC), recognized as the third most common
   urological malignancy, poses significant challenges in clinical and
   scientific domains for urologists. In 2020, over 400,000 new RCC cases
   were documented globally. ccRCC, the predominant subtype of RCC,
   constitutes about 70–80% of all RCC cases^[28]1,[29]2. Characterized by
   high metastasis and mortality rates, ccRCC is increasingly diagnosed in
   younger populations^[30]3. In 2022, China reported 77,410 new RCC cases
   and 46,345 deaths^[31]4. Although partial and radical nephrectomies are
   the primary ccRCC treatments, surgical intervention often doesn’t
   prevent recurrence in early-stage cases^[32]5,[33]6. The prognosis for
   ccRCC patients remains relatively poor. Early-stage ccRCC often
   presents without distinct clinical symptoms, resulting in approximately
   30% of patients receiving a diagnosis at advanced stages, characterized
   by distant metastasis and consequently missing opportunities for
   surgical intervention^[34]7. This underscores the critical need for
   innovative diagnostic and prognostic models, alongside the development
   of new biomarkers and molecular targets for ccRCC.

   Anoikis, a unique form of programmed cell death, is triggered when
   cells detach from adjacent cells or the extracellular matrix they
   typically adhere to. This mechanism is pivotal in the regulation of
   cellular survival, functioning by selectively eliminating cells that
   have abnormally detached from the adjacent structure^[35]8. Anoikis
   acts as a significant protective factor in both cancerous and
   non-cancerous diseases, particularly in inhibiting the dissemination of
   cancer cells that could lead to distant metastases^[36]9. Epithelial
   cells, which preserve normal tissue architecture through intercellular
   adhesion and engagement with the extracellular matrix, are particularly
   prone to anoikis^[37]10. Furthermore, epithelial cells constitute a
   critical component of tumor tissue. Tumor-Associated Epithelial Cells
   (TECs), in contrast to their normal counterparts, demonstrate
   resistance to anoikis, a critical precondition for tumor
   metastasis^[38]11. Anchorage-independent growth and the
   Epithelial-Mesenchymal Transition (EMT) are two pivotal hallmarks of
   anoikis resistance, significantly impacting tumor progression and the
   metastatic potential of cancer cells^[39]12. the activation of various
   signaling pathways plays a crucial role in imparting anoikis resistance
   to cancer cells, thereby enabling distant metastasis^[40]13. Extensive
   research has indicated that Anoikis-Related Genes (ARGs) are
   intricately linked with tumor progression and metastasis. The
   interaction between TIMP1 and CD63 activates the PI3K-AKT pathway,
   inducing resistance to anoikis in melanoma^[41]14, while CPT1A mediates
   fatty acid oxidation, promoting resistance to anoikis and inducing
   metastasis in colorectal cancer cells^[42]15.

   The objective of this research is to investigate the regulatory impact
   of ARGs on ccRCC TECs at a single-cell level, and to establish an
   innovative risk prognostic model grounded on ARGs for ccRCC
   Anoikis-related epithelial subclusters (ARESs). Utilizing multiple
   datasets, the study validates the predictive efficacy of this model.
   Additionally, it develops a nomogram related to clinical features to
   augment its applicability in clinical settings. Additionally, this
   research examines the disparities in pathway and functional enrichment
   among patients from distinct risk groups, as well as the variance in
   immune cell infiltration within the TME grounded on gene signatures.
   The goal is to potentially pioneer novel approaches for diagnosing,
   prognostically evaluating, and tailoring treatments for ccRCC patients.

Materials and methods

Data collection and research design

   Single-cell RNA sequencing (scRNA-seq) data were collected from tumor
   samples of seven patients with clear cell renal cell carcinoma (ccRCC)
   to explore the regulatory roles of ARGs in TECs. Additionally, spatial
   RNA sequencing (stRNA-seq) data from two ccRCC patients were downloaded
   to identify ARESs within the spatial context. We conducted thorough and
   extensive data mining to pinpoint robust DEGs between normal and tumor
   samples at the tissue level. The datasets included in this study
   adhered to the following criteria: each dataset was required to contain
   a minimum of 10 ccRCC tumor samples and corresponding adjacent normal
   samples. The entire scRNA-seq and stRNA-seq datasets, along with the
   microarray datasets, were obtained from the Gene Expression Omnibus
   (GEO) database, accessible at [43]www.ncbi.nlm.nih.gov under the
   accession numbers [44]GSE159115, [45]GSE210041, [46]GSE53757,
   [47]GSE36895, [48]GSE15641, [49]GSE66272, [50]GSE68417, and
   [51]GSE40435. The bulk RNA sequencing data and clinical details for
   TCGA-KIRC were acquired from the TCGA database, accessible at
   [52]https://cancergenome.nih.gov/. A supplementary validation dataset,
   E-MTAB-1980 was obtained from the ArrayExpress database, which can be
   accessed at [53]https://www.ebi.ac.uk/arrayexpress/. All of the ARGs
   were retrieved from Gene Card database
   ([54]https://www.genecards.org/), the protein coding genes and
   relevance score > 0.3 as filtering criteria. All data analyzed or
   generated in this study are freely accessible from prior publications
   or public databases.

ccRCC scRNA-seq data processing

   The 'Seurat' R package (version 4.4.0) was employed to preprocess the
   seven ccRCC scRNA-seq sample data. Seurat objects were generated for
   each sample, derived from their respective scRNA-seq gene expression
   matrix. Rigorous quality control was implemented, eliminating cells
   with gene expression counts below 200 or above 6000, raw counts under
   1000, or mitochondrial gene expression surpassing 20%. Data
   normalization was subsequently carried out using the NormalizeData
   function, employing the LogNormalize method with its standard settings.
   The top 2000 variable genes were calculated using the
   FindVariableFeatures function. Employing these genes, the ScaleData and
   RunPCA functions were executed on the Seurat objects, extracting the
   top 30 principal components (PCs) for subsequent analysis. To remove
   batch effects across samples, the ‘Harmony’ R package was employed,
   followed by Uniform Manifold Approximation and Projection (UMAP) to
   visualize distinct cell clusters in each scRNA-seq dataset. Ultimately,
   major cell types within the ccRCC TME were annotated and visualized,
   leveraging ccRCC cell annotation data derived from previous
   studies^[55]16.

Anoikis score and identify anoikis related TECs subcluster

   The expression score for ARGs was calculated using normalized data from
   each cell cluster, employing the 'Ucell' R package. To investigate the
   regulatory impact of ARGs on TECs, the non-negative matrix
   factorization (NMF) algorithm was utilized, specifically version 0.26
   of the NMF R package. This approach entailed conducting dimensionality
   reduction analysis on ARGs' expression in TECs, thereby categorizing
   distinct cell subtypes through the scRNA expression matrix. This
   analysis conformed to the methodologies outlined in previous relevant
   studies, guaranteeing consistency and comparability^[56]17,[57]18. DEGs
   within each NMF subcluster were calculated via the FindAllMarkers
   function, which was employed with its standard parameters. Subclusters
   exhibiting an average log2 fold change (log2FC) in ARGs exceeding 1,
   and manifested expression in over 70% of cells, were classified as
   ARESs.

Pseudotime analysis, cell communication, functional enrichment analysis and
metabolic analysis of ARESs

   The Monocle R package (version 2.30.0) was employed to analyze
   scRNA-seq data of ARESs, with the goal of elucidating the relationship
   between ARGs and cell pseudotime trajectories^[58]19. Dimensionality
   reduction via the DDRTree method enabled the visualization of ARGs'
   dynamic expression changes in ARECs' pseudotime trajectories in ccRCC,
   employing the plot_pseudotime_heatmap and plot_cell_trajectory
   functions for this purpose. Employing the Cellchat R package (version
   2.1.0), tailored for scRNA-seq data analysis across various cell
   clusters, predictions of ARESs' cellular interactions were conducted
   using the CellChatDB.human database^[59]20. Subsequently, the circle
   plot was utilized to visualize the intensity and extent of
   intercellular communication networks among various ARESs. The
   'netAnalysis_signalingRole_heatmap' function depicted the signaling
   pathway and their respective signaling input–output patterns. The
   scMetabolism function was applied to estimate cellular metabolism
   analysis, grounded on the related metabolic pathway enrichment results
   of ARESs’ DEGs. Finally, KEGG pathway enrichment analysis was performed
   on the DEGs between ARESs, employing the ClusterProfiler R package
   (version 4.10.0). Gene sets with an adjusted p-value less than 0.05
   were categorized as having significant enrichment, and the top three
   pathways with the highest enrichment were showcased.

Identifying ARESs in the spatial context of ccRCC patients

   To identify within the spatial context and address the challenge of
   inter-patient heterogeneity, we utilized stRNA-seq data from two
   patients with ccRCC. The cell types identified at the single-cell level
   in our study served as references for spatial context. The