Identification and validation of shared gene signature of kidney renal clear cell carcinoma and COVID-19

Data sources and preprocessing

We obtained the expression profiles of COVID-19 patients from the GEO database. Two cohorts, namely GSE196822 (GPL20301) and GSE211979 (GPL16791) consisting of a total of 47 samples (33 patient samples and 14 normal samples), were chosen. If multiple probes matched a single gene, the gene’s average expression was calculated. The study encompassed a total of 22,235 genes. We employed the “limma” R package to identify differentially expressed genes (DEGs) between normal and COVID-19 patient samples. These DEGs were named COVID-19-DEGs and filtered based on logFC ≥1 and adjusted p-value < 0.05.

We downloaded the expression profiles of KIRC patients from the TCGA database KIRC cohort (downloaded on 2023.01.03), which comprised 542 tumor tissue samples and 72 adjacent noncancerous tissue samples. If a gene had multiple probes, we calculated the average expression of that gene. A total of 60,550 genes were acquired. We utilized the “limma” R package to obtain DEGs between tumor and normal groups. These DEGs were named KIRC-DEGs and were filtered based on logFC ≥ 1.5 and FDR < 0.05. We also downloaded the corresponding mutation data and clinical information for KIRC patients from the TCGA database. We excluded samples with a survival time of 0 or missing information. Ultimately, we selected a final cohort of 533 patients for analysis.

To validate the prognostic risk model, we acquired the ICGC-RECA-EU cohort from the ICGC database and selected 91 KIRC patients with reliable prognostic information for external validation.

Identification and functional annotation of shared genes

To identify the most relevant genes associated with COVID-19, we employed weighted correlation network analysis (WGCNA) to classify the COVID-19-DEGs into distinct modules. We considered the module genes significantly associated with COVID-19 as candidate genes. The shared genes were obtained by intersecting these candidate genes with the KIRC-DEGs. To perform function enrichment analysis on the shared genes, we utilized the “clusterProfiler” R package, focusing on Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO).

Machine learning-based variable screening and model construction

The shared genes in the TCGA-KIRC training set were analyzed by univariate Cox regression to select the significant shared genes (p < 0.05), followed by Lasso-Cox regression analysis, which is a novel algorithm for screening parameters in high-dimensional data (Goeman, 2010). Through Lasso-Cox regression analysis, we further reduced the number of shared genes. Finally, using multivariate Cox regression analysis, we narrowed down the genes and constructed a risk model consisting of hub genes, which represent the signature genes of the model.

The risk score is calculated using the following formula:

$$RiskScore=\sum _{i=0}^n\left({\beta }_i\ast x_i\right)$$where n is the number of genes, β is the coefficient and x is the expression level of the gene.

Validation of hub genes expression

Through the previous analysis, we identified differential expression of hub genes in both COVID-19 and KIRC. To further validate our findings, we examined the expression of these hub genes in RCC cells.

Assessment of prognostic features and validated predictive power

The TCGA-KIRC training set was divided into low- and high-risk groups based on the median value of the hub genetic risk score. The survival rate between these groups was assessed using the Kaplan-Meier method and log-rank test with the “survminer” R package. Receiver operating characteristic (ROC) curves were generated using the “survivalROC” R package, and the area under the curve (AUC) was calculated to evaluate the predictive power of prognostic characteristics and other clinical parameters for patients’ overall survival (OS) at 1, 3, and 5 years. Furthermore, the accuracy and robustness of the risk model were confirmed in the TCGA-KIRC testing set, the complete TCGA-KIRC cohort, and the external validation set (ICGC-RECA-EU cohort).

Principal component analysis and development and validation of prediction models for nomograms

The discriminative capability of the risk model between low- and high-risk patients was assessed by generating Principal Component Analysis (PCA) scatterplots using the “scatterplot3d” R package. Clinical parameters, including gender, age, histological grade, pathological stage, and risk score, were incorporated to develop predictive nomograms for estimating the 1-year, 3-year, and 5-year survival rates of KIRC patients in the training set, testing set, entire set, and external validation set using the “rms” R package. Calibration plots were generated to evaluate the predictive performance of prognostic nomograms for patients’ OS at 1-year, 3-year, and 5-year.

Immune landscapes for risk model

The "ESTIMATE" algorithm was employed to evaluate the tumor microenvironment and determine the immune score, stromal score, and estimated score for each tumor sample using gene expression values from the entire TCGA-KIRC cohort (Yoshihara et al., 2013) Subsequently, we examined the profiles of 22 tumor-infiltrating immune cells using the “CIBERSORT” R package and determined their relative proportions in the entire TCGA-KIRC cohort (p < 0.05) (Newman et al., 2015). Additionally, we analyzed the associations between the risk model and immune function, as well as immune checkpoints.

Relationship between the risk model and tumor mutation burden

Tumor Mutational Burden (TMB) is a potential biomarker used to predict the response to immunotherapy, and patients with a high TMB often experience greater benefits from immunotherapy (Goodman et al., 2017; Yin et al., 2020). Mutation data were obtained from the TCGA-KIRC cohort. The TMB score was calculated using the formula: (total mutations/total covered bases) × 10⁶. We compared the TMB scores of the low- and high-risk groups and evaluated the survival differences between risk scores and TMB scores.

Construction of hub gene-based protein-protein interaction network and gene set enrichment analysis

We constructed a protein-protein interaction (PPI) network based on hub genes and their interacting genes using the GeneMANIA portal (https://genemania.org). This network was used to predict correlations among co-localization, co-expression, shared protein structural domains, and common pathways (Warde-Farley et al., 2010). Gene Set Enrichment Analysis (GSEA) is a visual method used to determine statistically significant and consistent differences in gene sets among different functional phenotypes (Subramanian et al., 2005). By comparing the enriched biological processes in the low and high-risk groups of the entire TCGA-KIRC cohort, we determined the potential biological functions of hub genes. The reference gene collection used in the GSEA was downloaded from the Molecular Signature Database (MSigDB), including c2.cp.kegg.v6.0.symbols.gmt and c5.go.v7.4.symbols.gmt. The critical criterion was set at p < 0.05.

Drug screening based on hub genes

To facilitate therapeutic development, we explored potential drug molecules significantly associated with hub genes using the Drug Signature Database (DSigDB) and cellMiner Database (Yoo et al., 2015; Shankavaram et al., 2009).

Statistical analysis

Statistical analyses were conducted using R software version 4.2.2 (R Core Team, 2022). Survival comparisons were performed using the Kaplan-Meier method and log-rank test. A significance level of p < 0.05 was used unless otherwise stated.

Identification and functional annotation of shared genes

Initially, we conducted a differential expression analysis on the COVID-19 dataset to identify genes with differential expression (COVID-19-DEGs), based on logFC ≥ 1 and an adjusted p-value < 0.05 (Fig. 1A). The COVID-19-DEGs were then partitioned into eight modules using WGCNA. The modules MEbrown, MEblue, MEpurple, MEblack, MEpink, and MEgrey, which demonstrated higher module significance (MS) indicating a stronger correlation with the disease, were selected for further analysis. These modules collectively comprised 556 genes (Fig. 1B). Subsequently, we performed an intersection between these genes and differentially expressed genes in kidney renal clear cell carcinoma (KIRC-DEGs) (n = 3,106) to identify shared genes (Figs. 1C and 1D). Heatmaps of the COVID-19-DEGs and KIRC-DEGs are presented in Figs. S1 and S2, respectively, while Fig. S3 illustrates a segment of the WGCNA analysis process.

Next, we conducted an enrichment analysis of the shared genes using KEGG and GO enrichment analysis. The GO enrichment analysis revealed predominant enrichment of immune-related functions in these genes, including phagocytosis, complement activation, immunoglobulin complexes, antigen binding, and immunoglobulin receptor binding, among others (Fig. 1E). The KEGG analysis indicated that the enriched pathway associated with these genes was the cell cycle (Fig. 1F). These findings suggest a close association between the shared genes and immune function, implying that immune-related functions are likely to be impacted during the onset of the disease.

Construction of the risk model

A random assignment was performed on the entire TCGA-KIRC cohort with a 1:1 ratio, followed by univariate Cox and Lasso-Cox regression analyses conducted on 156 genes within the training set (n = 265). Initially, univariate Cox regression analysis was employed to explore the association between the expression levels of selected genes in the training set and OS time. Utilizing a significance threshold of p < 0.05 in Cox analysis, we identified 83 genes as potential risk factors associated with OS (Table S1). Subsequently, this gene set underwent refinement through the calculation of regression coefficients using the LASSO regression algorithm (Figs. 2A and 2B). This process identified the nine most valuable predictive genes (Table S2), from which a risk scoring system was established using the formula described earlier. Further refinement of these nine genes was conducted through multivariate Cox regression analysis. Ultimately, we identified four hub genes based on a Cox p value of < 0.01, namely GTSE1, CEACAM4, HECW2, and KCNMA1, as features for constructing the risk model (Fig. 2C). Specifically, GTSE1 and CEACAM4 were identified as risk factors (HR > 1), whereas HECW2 and KCNMA1 were characterized as protective factors. The risk score for each KIRC patient could be calculated by utilizing the expression values of these four genes along with the coefficients obtained from the multivariate Cox regression, employing the following formula:

Risk Score = (0.786708751 * GTSE1 exp) + (0.716034712 * CEACAM4 exp) +

(−0.489951732 * HECW2 exp) + (−0.414148249 * KCNMA1 exp).

Validation of the risk model

The Risk Score was individually computed for each sample in the TCGA training set (n = 265) using the gene expression data and the risk score formula. Subsequently, the samples were divided into low and high-risk groups based on the median risk scores to evaluate the predictive performance of the risk model. Kaplan-Meier curves demonstrated a significantly improved prognosis for low-risk patients, with higher risk scores corresponding to shorter survival time and increased mortality rates (Fig. 2D). The effectiveness of the risk model in predicting patient survival at 1, 3, and 5 years was evaluated using ROC curves, revealing AUC values of 0.791, 0.755, and 0.782 for the respective time points (Figs. 2E and 2F). Survival analysis of the TCGA-KIRC testing set demonstrated a markedly worse prognosis for high-risk patients (p < 0.001), with corresponding AUC values of 0.698, 0.667, and 0.669 for 1, 3, and 5 years, respectively (Figs. 3A–3C). Comparable findings were observed across the entire TCGA-KIRC cohort and the external validation set (Figs. 3D–3I). Additionally, the distributions of risk scores, survival status, and gene expression of the model genes for each dataset were presented, further substantiating the robustness and accuracy of the risk model. The baseline characteristics for each dataset are provided in Table 1.

Development and validation of nomogram and PCA analysis

Nomograms were employed to validate the clinical applicability and utility of the risk models. Specifically, nomograms were developed to predict OS at 1, 3, and 5 years by integrating various clinical parameters, including age, gender, pathological stage, histological grade, and risk score (Fig. 4A). Calibration plots were utilized to assess the predictive accuracy of the nomogram prediction model, demonstrating similar trends to the ideal model in predicting patients’ OS at 1, 3, and 5 years, affirming a good fit between the predicted values of the risk model and the actual outcomes (Fig. 4B). Similarly, consistent trends were observed in the TCGA-KIRC testing set, the entire TCGA-KIRC cohort, and the external validation set (Figs. 4C–4H). The PCA scatter plot results across various cohorts indicated the risk model’s ability to distinctly differentiate between low- and high-risk patients (Figs. 4I–4L). Collectively, these indications suggest the accuracy and clinical applicability of the risk model.

Immune landscape between the low-risk and high-risk groups

Initially, we evaluated the tumor microenvironment within the entire TCGA-KIRC cohort using the ESTIMATE method. The analysis revealed substantially higher estimated and immune scores in the high-risk group compared to the low-risk group (p < 0.001) (Fig. 5A). Subsequently, we employed the CIBERSORT algorithm to determine the distribution of 22 immune cell types within the TCGA-KIRC cohort. The composition of these immune cell types was visualized as a histogram, depicting the relative proportions of distinct immune cells in each sample (Fig. S4). We further utilized the “ggplot2” R package to calculate correlations between risk scores and immune cells. The analysis demonstrated significant positive associations between risk scores and specific immune cell populations, including CD8⁺T cells, activated CD4⁺ memory T cells, and Tregs (Fig. 5B). Moreover, we examined the relationship between the expression levels of prevalent immune checkpoint molecules and the risk score. Our findings revealed a robust correlation between the risk score and these immune checkpoints (p < 0.05). Notably, CTLA4, LAG3, TIGIT, PDCD1, SIRPA, TNFRSF9, and LILRB1 exhibited significantly higher expression levels in high-risk KIRC patients compared to low-risk patients (p < 0.001) (Fig. 5C). Additionally, our analysis revealed a correlation between high-risk status and heightened immune-related functions, including APC_co_stimulation, T_cell_co-stimulation, and inflammation promotion (Fig. 5D). Furthermore, we investigated the association between TMB and the risk score. The findings indicated a higher TMB in high-risk patients compared to low-risk patients (p = 0.023). Additionally, we noted that individuals exhibiting both high TMB and high risk displayed the poorest prognosis, whereas those with low TMB and low risk exhibited the most favorable prognosis (Figs. 5E–5G).

Biological functions of hub genes

To investigate the biological functions of hub genes, we performed a PPI analysis of hub genes and their interacting genes using GeneMANIA, and the network illustrates that the four genes are enriched in DNA damage response, potassium channel complex, cell cycle checkpoint, microtubule polymerization, transmembrane transporter complex, regulation of ion transmembrane transport, detection of chemical stimulus, etc., (Fig. 6A). Then, we further performed GSEA pathway enrichment analysis for each gene in the low and high-risk groups of the entire TCGA-KIRC cohort. In the high-risk group, the top five GO terms were considered to be significantly enriched representatives. The results showed that hub genes were involved in immunoglobulin complex, humoral immune response mediated by circulating immunoglobulin, complement activation, phagocytosis recognition, and humoral immune response (Fig. 6B). The top five KEGG pathways in the high-risk group showed that hub genes were involved in complement and coagulation cascades, cytokine receptor interaction, retinol metabolism, oxidative phosphorylation, and drug metabolism cytochrome p450 (Fig. 6C).

Expression validation of hub genes

Analysis of public databases revealed significant overexpression of hub genes in both COVID-19 and KIRC (Figs. 6D and 6E). Additionally, qPCR results demonstrated variable degrees of overexpression of the four hub genes in RCC cell lines (Figs. 6F–6I). These findings further substantiate the accuracy and reliability of the risk model.

Drug screening

To identify potential drugs for the treatment of COVID-19 and KIRC co-morbidity, we conducted drug target enrichment analysis with four hub genes (GTSE1, CEACAM4, HECW2, and KCNMA1) as genetic targets associated with the construction of risk models via cellMiner and DsigDB databases (Tables S3, S4). The results revealed that Etoposide, Fulvestrant and Topotecan were significantly positively correlated with the expression level of the four hub genes (Table 2). These drugs serve as important references for the development of treatment strategies for comorbidity.