Biomarkers associated with cell-in-cell structure in kidney renal clear cell carcinoma based on transcriptome sequencing

Data acquisition

The current research encompasses datasets from four aspects. Firstly, RNA-Seq data for The Cancer Genome Atlas (TCGA)-KIRC were downloaded using the TCGA Genomic Data Commons (GDC) application programming interface (API). The FPKM values were converted to TPM and then log2-transformed, with a total of 513 primary tumor samples and 72 adjacent normal control samples retained. Secondly, expression profiles for The Renal Cell Cancer-European Union (RECA-EU)/Renal cell carcinoma were obtained from the International Cancer Genome Consortium (ICGC) database, encompassing 91 primary tumor samples and 45 adjacent normal control samples. Furthermore, the KIRC single-cell dataset GSE224630 was obtained from the Gene Expression Omnibus (GEO) database. The GSE224630 dataset includes tumor samples from 6 patients with untreated clear cell renal cell carcinoma. Finally, 101 CIC-correlated genes were sourced from previous literature (Song et al., 2022; Ren et al., 2024).

Identification and enrichment analysis of DEGs

The limma package (Ritchie et al., 2015) was utilized in the TCGA-KIRC and ICGC datasets to screen for DEGs between KIRC patients and normal controls, and to identify common upregulated genes across both datasets (Song et al., 2023a). The criteria for selecting DEGs were a p < 0.05 and |log2FC| > 1. Additionally, the clusterProfiler package in R (Yu et al., 2012) was employed to enrich the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) of these DEGs. Adjusted p < 0.05 indicated significantly enriched pathways (Song et al., 2023b).

Machine learning for key gene selection, diagnostic model development, and validation

The overlapping genes obtained from the intersection of DEGs and CIC genes were further screened using machine learning techniques. The rfe function from the R package caret (Kuhn, 2008) was employed, utilizing the svmlinear method. During the selection of the number of features through recursive feature elimination (RFE), the cross-validation (CV) accuracy of the support vector machine (SVM) model was used to screen for key disease genes. Feature selection was also conducted using Least Absolute Shrinkage and Selection Operator (LASSO) regression from the glmnet package (Friedman et al., 2021), with key parameters set to nfolds = 10 and family = ‘binomial’. Finally, the feature genes selected by both the SVM-RFE and LASSO methods were intersected to obtain the biomarkers for this study. Subsequently, the R package e1071 (Meyer et al., 2019) was used to construct a diagnostic model using the SVM method. The accuracy of the biomarkers selected by machine learning was tested by plotting receiver operating characteristic (ROC) curves for KIRC samples and normal controls in the TCGA-KIRC training set. A larger the area under the curve (AUC) indicated a higher accuracy of considering the genes as hub genes. The validity of these genes was further verified using the ICGC validation set using the same method.

Correlation analysis of biomarkers with immune infiltration

The CIBERSORT package (Newman et al., 2015) was employed to quantify different immune cells in KIRC samples and control samples. Spearman correlation coefficients were used to perform correlation analyses between biomarkers and immune cells.

Processing and analysis of scRNA-seq data from KIRC

The Read10X function from the Seurat package (Hao et al., 2024) was utilized to read the scRNA-seq data for each sample in GSE224630, retaining cells with a gene count between 200 and 5000 and a mitochondrial gene proportion of less than 10%. Subsequently, the SCTransform function (Hao et al., 2024) was applied for normalization. After Principal Component Analysis (PCA) dimensionality reduction, the harmony package (Korsunsky et al., 2019) was used to remove batch effects among different samples. The RunTSNE function (Hao et al., 2024) was then employed for t-Distributed Stochastic Neighbor Embedding (TSNE) dimensionality reduction. Finally, the FindNeighbors and FindClusters functions (Hao et al., 2024) were used for clustering, with parameters set to dims = 1:25 and resolution = 0.1. Cell subpopulations were annotated based on marker genes provided by the CellMarker2.0 database (Hu et al., 2022). The expression patterns of the biomarkers obtained in this study across different cell types were subsequently investigated.

Cell culture and transfection

Two cell lines, human embryonic kidney 293T (CRL-3216) and human renal clear cell adenocarcinoma cell 786-O (CRL-1932), were all obtained from the American Type Culture Collection (Manassas, MD, USA). These cells were cultured in DMEM (11965092, Gibco, Waltham, MA, USA) or RPMI 1640 medium (11875093, Gibco, Waltham, MA, USA) with the supplementation of 10% fetal bovine serum (FBS) (S9020, Solarbio Lifesciences, Beijing, China) and 1% penicillin-streptomycin (15140148, Gibco, Waltham, MA, USA). All cells were tested via short tandem repeat profiling and incubated in the incubator at 37 °C with 5% CO₂.

For the liposome transfection, the small interfering RNA against CTSS (si-CTSS) and the control small interfering RNA (si-NC) were all purchased from GenePharma (Shanghai, China) and transfected into 786-O cells with the use of lipofectamine 2000 transfection reagent (11668027, Invitrogen, Carlsbad, CA, USA) as per the manuals. The sequences applied for the transfection were 5′-TCACATATAAGTCAAACCCTA-3′.

CCK8 assay

During the logarithmic phase, 786-O cells were plated in a 96-well plate with a density of 1  × 10⁴ cells per well and incubated at 37 °C in an atmosphere containing 5% CO₂ for 0, 24, or 48 h. Subsequently, 10 µL of CCK-8 reagent was added to the medium, and the samples were incubated at 37 °C for 2 h. To generate the CCK-8 curve, the absorbance was measured at 450 nm, which served as the ordinate, while time was represented on the abscissa. The results were averaged from three independent experimental repeats.

Cell migration assay

The transfected 786-O cells (5  × 10⁵ cells/well) were cultivated in a six-well plate with serum-free media. When they achieved complete confluence, a 200-µL sterile pipette tip was used to create an artificial scratch on the monolayer. Forty-eight h later, the cells were photographed by an inverted optical microscope (DP27, Olympus, Tokyo, Japan), and the wound closure (%) was quantified correspondingly to assess the migration of KIRC cells (Zhang et al., 2024). Wound closure (%) = (Initial scratch width–scratch width at measurement time point)/Initial scratch width ×100%.

Cell invasion assay

For the invasion assay, 786-O cells (1  × 10⁵/100 µL) were suspended in 200 µL serum-free medium and plated in the upper Transwell chamber (3422, Corning, Inc., Corning, NY, USA) coated with matrix gel (C0372, Beyotime, China), while the lower chamber was filled with 700 µL culture media containing 10% bovine calf serum. After 48 h, the invaded cells were fixed by 4% paraformaldehyde (P0099, Beyotime, Shanghai, China) and stained with 0.1% crystal violet (C0121, Beyotime, Shanghai, China) for 30 min. Then, three random fields were observed under an inverted optical microscope (DP27, Olympus, Japan), and the number of invaded cells was quantified (Wang et al., 2023).

QRT-PCR experiment

Following the instructions, total RNA was isolated from 293T and 786-O cells using the TriZol total RNA extraction kit (15596026, Invitrogen, Carlsbad, CA, USA). Subsequently, the concentration of the isolated RNA was determined. Then, complementary DNA was synthesized by reverse transcription with a relevant assay kit (D7178S, Beyotime, Shanghai, China). After that, SYBR Green qPCR Mix (D7260, Beyotime, Shanghai, China) was used for the PCR assay according to the protocols. Finally, the relative level was calculated by the 2^−ΔΔCT method with GAPDH as the reference gene. The qRT-PCR primers used in this study were designed according to National Center for Biotechnology Information (NCBI) sequences using Primer Premier 6 software. The sequences of the primers used were presented in Table S1.

Statistical analysis

All statistical analysis was conducted using R (version 3.6.0). The t-test was performed on continuous variables between two groups. Data normality was assessed using the Shapiro–Wilk test. If normality assumptions were violated, the Wilcoxon rank-sum test was performed as a supplementary analysis. To investigate the correlation between gene expression and immune cell fractions, Spearman’s rank correlation test was employed. A p < 0.05 signified a significant level. Sangerbox (http://sangerbox.com/) provided assistance with this study (Shen et al., 2022).

Differential gene selection and enrichment analysis

Comprehensive differential gene expression analysis was conducted on tumor and control samples within the TCGA-KIRC and ICGC datasets (Figs. 1A–1B). Subsequently, upregulated genes common to both the TCGA and ICGC datasets were identified, totaling 1,256 (Fig. 1C). Functional enrichment analysis, including KEGG pathway analysis and GO analysis, was performed on these commonly upregulated genes. The top 10 enriched KEGG pathways included Human T-cell leukemia virus 1 infection, PI3K-Akt signaling pathway, Cytokine-cytokine receptor interaction, Phagosome, Epstein-Barr virus infection, and cell adhesion molecules (CAMs) among others (Fig. 1D). The results of the GO functional enrichment analysis covered biological process (BP), cellular component (CC), and molecular function (MF) aspects. The top five GO functional enrichments indicated that these genes were significantly involved in critical processes related to immune system regulation, such as regulation of lymphocyte activation, MHC protein complex, regulation of T cell activation, peptide antigen binding, T cell activation, regulation of leukocyte activation, and leukocyte cell–cell adhesion (Fig. 1E).

Machine learning for biomarkers screening

Twelve overlapping genes were obtained by intersecting the DEGs with CIC genes. When selecting the number of features using the RFE method, the trend of the CV accuracy of the SVM model with the number of features can be observed. It can be seen that when the number of selected features is six, the model’s CV accuracy reaches its maximum (Fig. 2A). To further screen the genes, LASSO regression analysis was employed. Figure 2B displayed the changes in regression coefficients for different gene features in LASSO regression as the penalty parameter (λ) varies, showing the deviance at the optimal λ value determined by 10-fold CV. By intersecting the genes obtained from both the RFE-SVM and LASSO methods, five hub genes were ultimately identified as biomarkers in this study, namely CDKN2A, VIM, TGFB1, CTSS, and CDC20 (Fig. 2C). Figures 2D–2E demonstrated the differential expression of these hub genes between tumor samples and normal controls, showing that all biomarkers are expressed higher in the tumor group than the control group in both the TCGA training set and the ICGC validation set (p < 0.0001).

Construction and verification of the diagnostic model

The five biomarkers were integrated to establish a diagnostic model, and the predictive ability of this model in the TCGA training set was assessed using ROC curves. Each curve represents a hub gene, and the AUC value was used to explore the gene’s ability to distinguish between diseases. The results showed that the AUC > 0.8 (Fig. 3A), indicating that these genes had high predictive power. Figure 3B displayed the ROC curve for the overall predictive ability of the constructed SVM model, with an AUC of 0.962, suggesting that the model exhibits very high predictive performance in distinguishing between KIRC and normal samples. The confusion matrix of the classification model demonstrated the classification results for the tumor group and the control group (Fig. 3C), further validating the robustness of the model in distinguishing between the two groups of samples. Subsequently, the same method was used for validation in the ICGC validation set. Both the AUC values for individual genes and the overall AUC value for the SVM model were above 0.8, indicating that the diagnostic model in this study has high predictive performance (Figs. 3D–3F).

Correlation between biomarkers and immune infiltration

To investigate the relationship between immune cells and key biomarkers, Figs. 4A–4B presented the distribution of different types of immune cells in normal and KIRC tissues. Overall, immune cells such as macrophages (M0, M1, M2) and T-cells (including CD8⁺ T-cells, regulatory T-cells), and others were significantly more abundant in tumor tissues compared to normal tissues (p < 0.05), while B-cells (B cells naive) were more prevalent in normal tissues (p < 0.05). These differences showed slight variations between the TCGA training set and the ICGC validation set, for example, there was a significant difference in M2 macrophages between the tumor group and the control group in the TCGA training set, but not in the ICGC validation set. However, the overall results were consistent. Finally, the correlation between hub genes and immune cell scores in tumor samples was calculated and presented using correlation coefficients (ranging from −1 to 1) and significance levels (indicated by asterisks). The results indicated that these genes have significant correlations with multiple types of immune cells, and these correlations were particularly concentrated in T-cells and macrophages (Figs. 4C–4D). Among them, T cells CD4 memory activated and T cells CD4 memory resting showed significant differences across 5 biomarkers in the TCGA training set (p < 0.05).

Relationship between hub genes and clinicopathological characteristics of KIRC

Next, we observed at the expression levels of the five hub genes in relation to clinical stage and grade. As shown in Fig. 5A, the expression levels of CDKN2A were significantly different in different stages and grades (p < 0.01). The expression of VIM and CTSS did not change significantly across different stages or grades (Figs. 5B, 5D). The trend in the expression level of the TGFB1 gene across different stages and grades was not significant, although there was a slight increase in Stage III and IV, it did not reach statistical significance (Fig. 5C). The expression level of CDC20 gene produced significant changes with the progression of clinical stage and grading, and its expression level was significantly up-regulated in Stage IV and G4 (Fig. 5E, p < 0.001). These results imply that these hub genes may be closely associated with tumor progression in KIRC.

Distribution and expression of hub genes in the kidney

Six major cell clusters were identified from public scRNA-seq datasets. These clusters included endothelial cells, Natural killer T (NKT) cells, fibroblasts, B cells, tumor cells, and proliferative tumor cells (Fig. 6A). Further gene expression characteristics were used to demonstrate the marker genes of various cell types. For instance, genes highly expressed in B cells included CD79A and MS4A1, endothelial cells specifically expressed PLVAP, VWF, and EMCN, fibroblasts specifically expressed COL3A1 and COL1A2, NKT cells specifically expressed GZMA, NKG7, and GNLY, while proliferative tumor cells and tumor cells expressed specific tumor marker genes such as MKI67 and EPCAM, respectively (Fig. 6B). Comparison of the proportions of various cell types in the samples revealed that tumor cells accounted for the highest proportion, reaching 42%, followed by endothelial cells and B cells, each accounting for 19%, fibroblasts accounting for 17%, and NKT cells and proliferative tumor cells accounting for the smallest proportions (Fig. 6C). Finally, the expressions of key hub genes in different cell types was further analyzed. These results demonstrated that TGFB1, CDKN2A, and CDC20 were significantly expressed in proliferative tumor cells, CTSS was most prominently expressed in B cells, while VIM was expressed in multiple cell types (Fig. 6D).

Cellular validation based on in vitro experiment

The possibilities of CDKN2A, VIM, TGFB1, CTSS and CDC20 as biomarkers of KIRC were further verified through in vitro experiments. The mRNA expression levels of these five genes were calculated, and it was found that the expressions of all these five genes in 786-O KIRC cells were significantly higher than those in 293T cells (Fig. 7A, p < 0.05). Since we observed a significant up-regulation of CTSS expression in KIRC cell lines relative to the other key genes screened, and its less studied in CTSS in KIRC. For this reason, the impact of CTSS silencing on the biological function of KIRC cells was investigated (Fig. 7B). CCK-8 results showed that silencing the expression of CTSS significantly reduced the proliferative capacity of 786-O cells (Fig. 7C, p < 0.001). Furthermore, the relevant results demonstrated that the silencing of CTSS led to a reduction in the migration and invasion capabilities of KIRC cells (Figs. 7D–7E, p < 0.001).