Abstract
Objective. This study aimed to analyze the cuproptosis-related long non-coding RNA (lncRNA) in patients with bladder urothelial carcinoma (BLCA), construct a prognostic model, and screen its potential drugs. Methods. The transcriptome expression and clinical and mutation burden data related to BLCA were downloaded from The Cancer Genome Atlas database. The prognostic lncRNAs were screened using univariate Cox and Lasso regression analyses, and then included in the multifactor risk ratio model. The risk score of each sample was calculated based on the prognostic model formula, and the patients were divided into high- and low-risk groups for survival difference analysis. Clinically relevant receiver operating characteristic (ROC) curve, C-index principal component analysis, and clinical data statistics were used to evaluate the predictive power of the model. The risk-differential lncRNAs were functionally enriched. We calculated the tumor mutation burden of risk lncRNAs, and survival and the Tumor Immune Dysfunction and Exclusion analyses for high- and low-risk groups. Finally, immunocorrelation analysis and potential drug screening were performed. Results. Eleven lncRNAs with independent prognostic significance were screened out to construct the prognostic model. Survival analysis showed a significant difference in survival between the high- and low-risk groups. The areas under the ROC curve at 1, 3, and 5 years were 0.711, 0.679, and 0.713, respectively. The discrimination between the lncRNA high- and low-risk groups in the constructed model was the most obvious. The risk-differential lncRNAs were closely related to immunity. The treatment drugs with high sensitivity were screened based on the IC50 value. Conclusion. The 11 cuproptosis-related lncRNAs may serve as molecular biomarkers and therapeutic targets for BLCA.
1. Introduction
Bladder cancer (BC) is a malignant tumor occurring on the bladder mucosa, among which transitional cell carcinoma is the most common one and is called bladder urothelial carcinoma (BLCA). BC is the most common malignant tumor of the urinary system in China and ranks 10th in the incidence of all types of cancer in the world [1, 2]. Despite medical developments, improved diagnostics, and increased health awareness in recent years, the incidence has not been stable throughout the world over time, nor will it be in the near future. Countermeasures should be taken to deal with the adverse effects of aging [3].
The main treatments of BC at present include surgery, radiotherapy, chemotherapy, immune support, and so on. Invasive BC accounts for about 70% of all BCs, and the muscle layer is generally involved in the resection of transurethral bladder tumor, but the possibility of recurrence and progression is higher. Muscularis invasive BC accounts for about 30% of all BCs, and radical resection and bladder pelvic lymph node cleaning treatment are the gold standard. However, still, a high recurrence rate or the possibility of progression exists [4, 5]. Currently, authoritative guidelines recommend immune checkpoint inhibitors as the second-line treatment in patients with BC who have developed metastases and lost the chance of surgical resection, and as the first-line treatment in patients with programmed death-ligand 1 (PD-L1) who are ineligible for platinum-based chemotherapy [6].
A variety of prespecified and precisely controlled programmed cell death occur, such as apoptosis, necroptosis, pyroptosis, and ferroptosis, during the development of multicellular organisms. The mechanism of copper toxicity has been demonstrated to be different from all other known mechanisms regulating cell death, and this previously uncharacterized cell death mechanism is termed cuproptosis [7].
Therefore, this study used The Cancer Genome Atlas (TCGA) database to predict the long non-coding RNA (lncRNA) associated with copper mortality in BC, so as to construct the prognostic model of related lncRNA and screen potential drugs in patients with BLCA.
2. Materials and Methods
2.1. Data Downloading and Sample Sorting
The BLCA-related transcriptome expression data, clinical data, and the mutation load were downloaded from the National Cancer Institute GDC Data Portal website (https://portal.gdc.cancer.gov), which included the control and the patient’s clinical data, such as age, sex, survival time and survival state, tumor classification, T stage, N stage, and M installment. The Perl software was used to collate the transcriptome data and transformation ID, and separate lncRNA and microRNA.
2.2. Expression of Cuproptosis-Related lncRNA and Co-Expression Analysis
The limma package of the R-Studio software was used to collate transcriptome data and perform cuproptosis-related gene (NFE2L2, NLRP3, ATP7B, ATP7A, SLC31A1, FDX1, LIAS, LIPT1, LIPT2, DLD, DLAT, PDHA1, PDHB, MTF1, NFE2L2, NLRP3, ATP7B, ATP7A, SLC31A1, FDX1, LIAS, LIPT1, LIPT2, DLD, DLAT, PDHA1, PDHB, MTF1, GLS, CDKN2A, DBT, and GCSH DLST) analysis. The extraction of the lncRNA expression quantity associated with cuproptosis via the limma package quantity of cuproptosis-related gene expression was analyzed via lncRNA isolation (set conditions: corFilter = 0.4 and ).
2.3. Construction of the Prognostic Model
The limma package of the R-Studio software was used to merge cuproptosis-related lncRNA expression data with the clinical data (survival time and survival state) via the survival pack, the caret bag, glmnet package, survminer merged data and time receiver operating characteristic (ROC) package, and univariate Cox regression analysis (filter criteria: ) to screen the lncRNA related to the prognosis in patients with BLCA. Further screening analysis was performed by Lasso regression to reduce overfitting of the data and to screen the key cuproptosis-related lncRNAs. Cross-validation was used in Lasso regression to select parameters, and the Lasso regression coefficient spectrum was drawn. Finally, multivariate Cox regression analysis was performed to establish a valuable lncRNA model for the prognosis of BLCA and presented in the form of a nomogram. The lncRNA model was constructed based on the multivariate Cox regression analysis, and the risk score equation was: . The data were randomly divided into two groups: Train group and Test group.
2.4. Evaluation and Clinical Value Analysis of the Prognostic Model
The R-Studio software was used to perform clinical statistical analysis between the Train and Test groups. Limma, reshape2, tidyverse, and ggplot2 packages were used for correlation analysis of cuproptosis-related gene expression data, cuproptosis-related lncRNA expression data, and risk gene data. The survival, survMiner, and timeROC packages were used to analyze the OS and progression-free survival (PFS) of risk gene data. Independent prognostic and ROC analyses were performed on risk gene and clinical data. The C-index analysis was performed using dplyr, survival, rms, and pec packages. Risk genetic data and clinical data were analyzed by using regplot, survival and rms packages, and then histogram was drawn to predict the survival of BLCA patients by nomogram. Survival was analyzed using the R-Studio software package and survminer package data for risk genes associated with clinical data, and the validation model built was applicable to patients with different clinical groups. The principal component analysis (PCA) was performed using the limma and Scatterplot3d packages to verify whether the lncRNA involved in the model could distinguish the patients in the high- and low-risk groups.
2.5. Screening of Risk Genes with Significant Differences and Enrichment and Immune-Related Functional Analyses
The R-Studio software was used to conduct risk difference analysis on risk gene data and gene expression data, and to screen the differential risk genes with significant differences. The Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed. Then, risk gene data and gene expression data were analyzed for immune-related functions, and statistically significant differences in immune-related functions were observed between the high- and low-risk groups.
2.6. Comparison between High and Low Risk of BC Risk Gene Mutation Frequency and Mutation Load Survival Analysis
The downloaded mutation burden data were processed using the Perl software, and the mutation burden of BLCA was calculated. The Perl software was used to process the tumor mutation data and risk gene data to obtain the mutation gene data of high- and low-risk groups. The MafTools package of the R-Studio software was used to analyze the mutation frequency of risk genes between high- and low-risk groups. Reuse the survival of cancer mutations and the survminer package load data and risk genes in the tumor mutation load of survival analysis.
2.7. Analysis of Immune Evasion and Immunotherapy in BLCA
The transcriptome data were uploaded to the Tumor Immune Dysfunction and Exclusion (TIDE) database (http://tide.dfci.harvard.edu/) to obtain the TIDE rating of the transcriptome data. The R-Studio software limma package and ggpubr package were used to analyze risk genes data. The high- and low-TIDE scores were analyzed to find the difference between the risk groups.
2.8. BC Screening of Potential Drugs
Gene expression data and risk genes were analyzed using the limma, ggpubr, pRRophetic bags, and ggplot2 packages, and screening indicated a significant difference between high- and low-risk groups of drugs (filter condition: ).
3. Results
3.1. Extraction of Cuproptosis-Related lncRNA Expression and Co-Expression Analysis
Using the Perl transcriptome data processing software, we obtained 16,876 lncRNAs from the data set of RNA-seq of 412 samples of BLCA and 19 corresponding tissue adjacent to carcinoma samples after isolating lncRNA expression data. The R-Studio software was used to analyze the transcriptome data of cuproptosis-related genes (NFE2L2, NLRP3, ATP7B, ATP7A, SLC31A1, FDX1, LIAS, LIPT1, LIPT2, DLD, DLAT, PDHA1, PDHB, MTF1, GLS, CDKN2A, NFE2L2, NLRP3, ATP7B, ATP7A, SLC31A1, FDX1, LIAS, LIPT1, LIPT2, DLD, DLAT, PDHA1, PDHB, MTF1, GLS, CDKN2A, DBT, and GCSH DLST). The expression levels of cuproptosis-related genes in each sample were obtained through analysis, and then combined with the obtained 16,876 lncrnas for co-expression analysis, a total of 762 copper-death co-expression lncrnas were screened (corFilter = 0.4, ; Table 1). The R-Studio software was used to draw the mulberry map of the co-expression data, and the correlation between cuproptosis-related genes and lncRNA could be intuitively observed (Figure 1).

3.2. Construction of Prognostic Model
The R-Studio software was used to merge the lncRNA expression files related to cuproptosis with clinical data (survival time and survival status). Furthermore, 67 lncRNAs related to the prognosis of patients with BLCA were initially screened using univariate Cox regression analysis (; Figure 2(a)), and then, Lasso regression analysis was used to reduce the overfitting of data. Seventeen cuproptosis-related lncRNAs that were more valuable for the prognosis of patients with BLCA were screened (Figures 2(b) and 2(c)). Finally, 11 lncRNAs with prognostic values for copper mortality were screened out via multivariate Cox regression analysis (Table 2), and the prognostic model was constructed based on the risk scores of the screened-out 11 lncRNAs. The samples were randomly divided into two groups (Train and Test groups). The samples in the Train and Test groups were divided into high- and low-risk groups. The risk score was obtained using the following formula:

(a)

(b)

(c)
3.3. Evaluation and Clinical Value Analyses of the Prognostic Model
Combined with the clinical data of patients with ccRCC, R-Studio software was used to conduct clinical statistical analysis of the training group and the Test group. No significant difference was observed in each clinical trait between the Train and Test groups, indicating no deviation in clinical traits between the random grouping of the samples (; Table 3). Correlation analysis was performed on the expression data of cuproptosis-related genes, cuproptosis-related lncRNA expression data, and risk gene data to know the correlation between lncRNA involved in the model construction and cuproptosis-related genes (Figure 3). The OS analysis of the risk gene data indicated significant differences in the survival rate between the high- and low-risk groups in the risk gene data, in both the Train and Test groups, and the survival of patients in the high-risk group in the risk gene data was shorter in both the Train and Test groups compared with those in the low-risk group. The PFS analysis of risk gene data showed a significant difference in the progression-free survival between the high- and low-risk groups, and the high-risk group had significantly shorter progression-free survival than the low-risk group (Figure 4). The R-Studio software was used to Train set and Test set the risk of genetic data is analyzed, through the risk curve can be intuitive grouping situation of both high- and low-risk groups (median of risk score), through the survival state diagram can be found that patients with increased risk, the cases of death also along with the increase, and by heat maps can be observed, CASC20, NR2F2-AS1, AC124283.3, AP000593.3 are high-risk lncRNA, secretory carrier membrane protein 1 (SCAMP1)-AS1, AC108066.2, AC110611.1, AC022165.1, MIR181A2HG, AC010132.4, and AC078778.1 were low-risk lncRNAs (Figure 5). Univariate and multivariate Cox regression analyses, via independent prognostic analysis of risk gene data and clinical relevant data (sex, age, grade, and stage), indicated that the constructed model could be used as an independent prognostic factor in patients with BLCA independent of other clinical traits (risk score ; Figure 6). Combining the ROC curve with clinical relevant data revealed that the constructed model predicted the survival time of patients Area Under Curve (AUC) = 0.711) better than age, sex, grade, and stage. The ROC curve showed that the constructed model predicted the survival time in patients with BLCA with high sensitivity and accuracy (AUC after 1, 3, and 5 years was 0.711, 0.679, and 0.713, respectively). Furthermore, the C-index curve showed that the model constructed had a high accuracy in predicting the survival of patients (Figures 7(a), 7(b), and 7(c)). The survival nomogram of patients with BLCA was drawn by combining the risk gene data with clinical relevant data, and the risk of patients could be scored by the nomogram to predict the survival rate after 1, 3, and 5 years. The calibration chart showed that the predicted probability was basically consistent with the actual probability (Figures 7(d) and 7(e)). The analysis of risk gene data and clinical relevant data (stages) revealed that the constructed model was suitable for survival prediction in both Stages I and II and Stages III and IV patients (; Figure 8). PCA showed that the discrimination level between the lncRNA high- and low-risk groups in the constructed model was the most obvious, indicating that the patients in the high- and low-risk groups could be distinguished by the lncRNA involved in the model construction (Figure 9).


(a)

(b)

(c)

(d)

(a)

(b)

(c)

(d)

(e)

(f)

(a)

(b)

(a)

(b)

(c)

(d)

(e)


(a)

(b)

(c)

(d)
These results indicated that the prediction model based on 11 lncRNAs associated with copper mortality was significantly better than the clinical factors, such as age, sex, grade, and tumor stage, in predicting the prognosis of patients, and the risk score was significantly correlated with the progression of BLCA.
3.4. Screening and Enrichment Analysis of Risk Genes with Significant Differences
The R-Studio software was used for the differential analysis of gene expression data and risk gene data, and a total of 1042 genes with significant differences between the high- and low-risk groups were obtained (Table 4). The GO and KEGG enrichment analyses were further performed on the differential risk genes to analyze the enrichment pathways of differential risk genes in tumor tissues. The GO enrichment analysis showed that in the biological process of GO, differential risk genes were enriched in humoral immune response, extracellular matrix organization, extracellular structure organization, defense response to bacterium, phagocytosis, humoral immune response mediated by circulating immunoglobulin, complement activation classical pathway, collagen-containing extracellular matrix, immunoglobulin complex, circulating immunoglobulin complex, antigen binding, glycosaminoglycan binding, extracellular matrix structural constituent, immunoglobulin receptor binding, sulfur compound binding, cytokine receptor binding, cytokine activity, chemokine activity, and so on. The differential risk genes were enriched on GO:0003823, GO:0005539, GO:0005201, GO:0006959, GO:0045229, GO:0030198, GO:0043062, GO:0042742, GO:0006909, GO:0062023, and GO:0045229. The enrichment on GO:0019814 and GO:0009897 was significant (Figure 10(a)). The KEGG enrichment analysis showed that differential risk genes were enriched in cytokine–cytokine receptor interaction, viral protein interaction with cytokine and cytokine receptor, interleukin (IL)-17 signaling pathway, osteoclast differentiation, extracellular matrix–receptor interaction, protein digestion and absorption, phosphoinositide 3-kinase–Akt signaling pathway, chemokine signaling pathway, cell adhesion molecules, necrosis factor–kappa B signaling pathway, focal adhesion, proteoglycans in cancer, tumor necrosis factor signaling pathway, AGE–RAGE signaling pathway in diabetic complications, Janus kinase–signal transducer and activator of transcription signaling pathway, and transcriptional misregulation in cancer. The significant enrichment was observed on HSA05150, HSA04060, HSA04061, HSA05323, HSA05146, HSA04145, and HSA04640 (Figure 10(b)).

(a)

(b)
The immune-related function analysis of risk gene data and tumor gene expression data showed APC_co_inhibition, T_cell_co-inhibition, checkpoint, T_cell_co-stimulation, cytolytic_activity, inflammation promotion, Human Leukocyte Antigen (HLA), APC_co_stimulation, Chemokine receptors (CCR), Major Histocompatibility Complex (MHC)_class_I, and parainflammation. Significant differences were observed in Type_I_IFN_Reponse between the high- and low-risk groups (). Meanwhile, the heat map showed that the aforementioned immune-related functions were more active in the high-risk group (Figure 11).

3.5. Survival Analysis to Compare the Mutation Frequency and Mutation Burden of Risk Genes between High- and Low-Risk Groups of BLCA
The downloaded mutation burden data were processed using the Perl software, and the mutation burden of BLCA was calculated (Table 5). Risk genes for high- and low-risk groups of data of mutated genes, using the R-Studio software analysis between high- and low-risk group risk gene mutation frequency, and draw a waterfall figure, can be found through the Perl software for tumor mutation data and data processing. The mutation frequency of risk genes in the high-risk group was mostly higher than that in the low-risk group (Figure 12(a)). Survival analysis of tumor mutation load data and risk genes showed that there was a significant difference in the survival time of patients with high and low mutation load, and the survival time of the group with high BLCA mutation load was longer than that of the group with low mutation load (; Figure 12(b)), high mutation load groups of high- and low-risk patients survival time has significant differences. The survival time of high-risk patients was shorter than that of low-risk patients (). A significant difference was observed in the survival time between high- and low-risk patients in the low-mutation burden group, with high-risk patients having a shorter survival time than low-risk patients (; Figure 12(c)).

(a)

(b)

(c)

(d)
3.6. Analysis of Immune Evasion and Immunotherapy in BLCA
The analysis of immune escape and immunotherapy of risk gene data indicated significant differences in TIDE scores between the high- and low-risk groups, and the TIDE scores in the high-risk group were significantly higher than that in the low-risk group, indicating that the immune escape potential in patients with BLCA in the high-risk group was greater, and hence, the effect of immunotherapy was worse (; Figure 13).

3.7. Screening of Potential Drugs for BLCA
The gene expression data and risk data for potential drug screening were analyzed using the R-Studio software (screening conditions: ). From the IC50 values of drugs in the high-risk and low-risk groups, it could be found that bryostatin 1, doxorubicin, epothilone B, cyclopamine, dasatinib, FMK (Fluoromethyl Ketone), Genentech Cpd 10, obatoclax mesylate, paclitaxel, rapamycin, ruxolitinib, parthenolide, pazopanib, saracatinib, STF-62247, and other drugs indicated significant differences between the high- and low-risk groups, were more sensitive to the high-risk group, and could be found according to the correlation, IC50 value concentration of the drug, and negatively correlated with the risk score, the higher risk score, IC50 value concentration decreases (Figure 14).

4. Discussion
BLCA is one of the most common cancers in the urinary system. Although the diagnosis and treatment of BLCA have definite curative effects in the clinic, the disease may recur even after surgery. Therefore, screening of new BLCA biomarkers in the high-risk group is urgently needed for early and individualized treatment of patients with BLCA using potential drugs to improve the survival rate of patients. Cuproptosis is copper-dependent programmed cell death. The progress in human genome sequencing data indicated that lncRNAs play an indispensable role in the diagnosis, treatment, and prognosis of tumors [8–10]. Recent studies have shown that lncRNAs are involved in the occurrence and development of BLCA [11–13]. BLCA has been found to be associated with ferroptosis [14]. However, only a few studies on cuproptosis-related genes have been conducted. Therefore, it was of significant help to establish a prediction model of cuproptosis-related lncRNA based on the TCGA database to predict the prognosis of patients with BLCA and find new biomarkers for individualized treatment of high-risk patients.
In this study, 11 lncRNAs with independent prognostic significance were screened using univariate Cox regression analysis, Lasso regression analysis, and multivariate Cox regression analysis to construct a prognostic model. Previous studies have shown that CASC20 is involved in the carcinogenesis of gastric cancer by competitively binding with Mir-143-5p and regulating the MEMO1 expression, thereby inducing epithelial–mesenchymal transition [15]. The SCAMP1 is a protein-coding gene. Among its related pathways is an innate immune system that plays an important role in the occurrence and development of ovarian cancer and glioma [16, 17]. SCAMP1-AS1 could inhibit the proliferation and migration of esophageal cancer TE-13 cells by targeting miR-483-5p, and SCAMP1-AS1 was downregulated in esophageal cancer cells [18]. Studies have shown that AC078778.1 is associated with the occurrence and development of BLCA [19]. NR2F2 antisense RNA 1 (NR2F2-AS1) expression is associated with renal clear cell carcinoma, gastric cancer, thyroid cancer, and other diseases. In thyroid cancer, NR2F2-AS1 promotes the proliferation and migration of thyroid cancer cells and inhibits cell death by regulating the miRNA-338-3p/CCND1 axis [20–22]. Studies have shown that downregulation of lncRNA MIR181A2HG (MIR181A2 host gene) by high glucose impairs vascular endothelial cell proliferation and migration through the dysregulation of the miRNAs/Akt2 axis [23]. At present, the remaining six lncRNAs (AC108066.2, AC110611.1, AC022165.1, AC010132.4, AC124283.3, and AP000593.3) associated with cuproptosis have not been reported in tumors; hence, further studies are needed to explore their significance. The BLCA prognostic model constructed based on 11 cuproptosis-related lncRNAs showed that the OS and PFS in the low-risk group were significantly better than those in the high-risk group in the survival analysis of the BLCA samples in the TCGA database. The survival status and risk score diagram revealed that the number of patients who died and the risk score continuously increased with the increase in the risk value. The analysis of the ROC curve and calibration chart indicated that the model had high sensitivity and accuracy in predicting the prognosis of BLCA, which could provide a potential direction for clinical research.
Most of the pathways were found to be correlated with immunity via enrichment analysis and immune-related function analysis of differentially expressed risk genes, and significant differences were observed in the immune-related functions. As parainflammation characterizing 50% of major cancer types has been shown to be correlated with p53 mutations and defective p53 pathways, parainflammation may serve as a driver of p53 mutagenesis and help in cancer prevention when treated with nonsteroidal anti-inflammatory drugs [24, 25]. APC_co_stimulation, MHC_class_I, APC_co_inhibition, HLA, and Type_I_IFN_Reponse were significantly different among the BLCA subtypes. The study of the intratumoral immune microenvironment may provide a new perspective for BLCA treatment [26, 27]. Some studies found that T-cell co-stimulation in combination with targeting Focal adhesion kinase (FAK) drives enhanced anti-tumor immunity [28]. At present, no conclusive evidence suggests that cuproptosis is directly correlated with the occurrence of BLCA. However, this study will help explore the mechanism of cuproptosis-related lncRNAs.
Tumor mutational burden (TMB) was defined as the total number of detected somatic gene coding errors, base substitutions, and gene insertions or deletions per megabase [29]. PD-1 or its ligand (PD-L1) has achieved remarkable clinical efficacy in treating a variety of tumors. TMB, the newest marker for evaluating the efficacy of PD-1 antibodies, has been demonstrated in treating colorectal cancer with a deficiency in mismatch repair [30, 31]. The TMB level in tumor tissue can predict the efficacy of targeted therapy in patients with advanced non-small-cell lung cancer with driver mutations [32]. This study indicated that TMB could be considered an independent prognostic factor in BLCA.
Collectively, our results suggest that the proposed model, particularly the identified lncRNAs, could be promising biomarkers for estimating the prognosis of BC patients and could help clinicians to stratify patients into high- and low-risk groups, which could lead to more personalized treatment strategies. Furthermore, these lncRNAs could be investigated as therapeutic targets, as more in-depth understanding on their functional roles in BC and their involvement in cuproptosis could lead to the development of novel therapeutic approaches targeting these molecules or their downstream effectors. Thus, further research could explore the relationship between cuproptosis, lncRNAs, and BC in more detail, including studies on their underlying molecular mechanisms. Additionally, the proposed prognostic model could be validated using additional independent datasets, which would strengthen its potential use in clinical practice.
In conclusion, this study successfully constructed a prognostic model based on 11 promising lncRNAs with independent prognostic significance for patients with BC. The model effectively discriminated between high- and low-risk groups and demonstrated a strong association with immunity. Moreover, we identified potential treatment drugs with high sensitivity based on IC50 values, which provide important references for the potential development of personalized treatment strategies for BLCA patients.
Data Availability
Transcriptome expression data, clinical data, and mutation load data related to BLCA were downloaded from TCGA database (https://portal.gdc.cancer.gov). The transcriptome data uploaded to TIDE database (http://tide.dfci.harvard.edu/).
Conflicts of Interest
The author(s) declare(s) that they have no conflicts of interest.
Acknowledgments
We would like to thank TCGA database and TIDE database for providing high-quality data on BLCA. This work was supported by the National Natural Science Foundation of China (81800668) and Natural Science Foundation of Shanxi Province (201801D221276).
Supplementary Materials
Supplementary Materials. SupplementaryTable S1: The specific data description of Figure 1(a).
Supplementary Materials. SupplementaryTable S2: All data obtained after GO enrichment analysis of differential risk genes.
Supplementary Materials. SupplementaryTable S3: All data obtained after KEGG enrichment analysis of differential risk genes.