The 24 angiogenesis-related genes investigated in this study were obtained from the Gene Set Enrichment Analysis (GSEA) package on the WikiPathways website (https://www.gsea-msigdb.org/gsea/index.jsp). The TCGA datasets (https://portal.gdc.cancer.gov) were downloaded to obtain the CNVs, SNVs, and expression level changes in these genes in 32 different types of cancer. The data were analyzed with the Perl language and R Studio, and the results were visualized using Toolbox for Biologists (TBtools) (19). Additionally, RNA-Seq transcriptome data and associated clinical data on the patients with KIRC were obtained from TCGA, which contains 539 tumor samples and 72 normal samples. “Corrgram” in RStudio (MA, USA) was used to plot the co-expression relationships among the angiogenesis-related genes, with Pearson's correlation coefficient applied in the statistical analysis of the results.

Due to the large variations in gene expression profiles observed in the previously obtained datasets, we constructed an angiogenesis score model based on mRNA expression to show the differential expression levels between the samples. In brief, a single-sample gene set enrichment analysis (ssGSEA) was first used to evaluate the enrichment scores of genes in the angiogenesis pathway. The “Gplots” package in RStudio was used to perform the differential analysis, and the “pheatmap” package was used to draw the heat map of the cluster analysis results. Through comparisons with the mRNA expression levels of the genes in normal tissues, the mRNA expression statuses in the tumor tissues were classified into three categories: angiogenesis inactive (cluster 1 or C1), angiogenesis active (cluster 2 or C2), and normal angiogenesis (cluster 3 or C3). To further illustrate the relationships between the gene expression levels of these three clusters, we used the violin plot drawn with the “ggpubr” package to depict the enrichment score levels of the clusters. Finally, “pheatmap” in RStudio was used to draw a heat map of the relationship between two of the clusters and the clinicopathological characteristics of the patients with KIRC. The differences with a p < 0.05 were considered statistically significant.

Of the 266 drugs listed in the 2019 edition of the Genomics of Drug Sensitivity in Cancer (GDSC) database (https://www.cancerrxgene.org/), which is the largest of its type on pharmacogenomics, we selected 12 classical and novel targeted drugs used for the KIRC treatment for study. On the basis of the cell expression profiles indicated on the GDSC database, we used the “pRRopheticl” package to construct a ridge regression model to estimate the half-maximal inhibitory concentration (IC₅₀) of the drugs in the three clusters (20). We also used 10-fold cross-validation based on the GDSC training set to evaluate the prediction accuracy. Except for “combat” and “allSoldTumours” tissue types, all the parameters were set to default values, and the duplicate gene expression levels were summarized as the mean value. Finally, the “ggplot2” and “cowplot” packages were used to draw the boxplot. The results with a p < 0.05 were considered statistically significant.

To explore the potential regulatory mechanism of the angiogenesis pathway in KIRC, we examined the expression levels of various oncogenes and tumor suppressor genes in the three clusters in the form of a heat map, using the “string,” “pheatmap,” “gplots,” and “gird” packages in RStudio. One-way ANOVA was used to compare the expression levels of various oncogenes and tumor suppressor genes in KIRC in the different clusters. The differences with a p < 0.05 were considered statistically significant. Similarly, we used the same approach to demonstrate the differences among the three angiogenesis clusters in terms of the expression of sirtuins (SIRTs) and histone deacetylases (HDACs), which aside from being involved in histone modification also play an important role in regulating the angiogenesis pathway.

A ssGSEA was used to quantify the 29 immune-associated gene sets obtained from TCGA and there are a total of 707 genes, representing different immune cell types, functions, and pathways (21–23). The “ggplot2” and “dplyr” packages in R Studio were then used to draw a heat map of the correlation between the angiogenesis-related genes and immune cell infiltration, with Spearman's correlation coefficient being applied for the statistical analysis. A ssGSEA can be applied to gene signals expressed by the immune cell populations in a single sample. Therefore, on the basis of the ssGSEA results, we used the “ggstatsplot,” “data.table,” “dplyr,” “tidyr,” and “ggplot2” packages in RStudio to analyze and plot the correlation between the angiogenesis score and immune substances. In the plotted figure, the area of each sphere represents the degree of correlation and the color represents the p-value. Finally, we used the “ggscatterstats” package to generate a scatter plot for representing the correlations between the four classical immune cell populations (viz., T follicular helper (Tfh) cells, mast cells, neutrophils, and type II interferon (IFN) response cells) and the angiogenesis score. A p < 0.05 was considered to be statistically significant.

We used the “corrplot” package to show the co-expression relationship between any two of the angiogenesis pathway genes. The “glmnet” package in RStudio was used to perform a least absolute shrinkage and selection operator (LASSO)–Cox regression analysis to further identify the most useful prognostic genes and to establish a risk model. Next, we computed the risk score (RS) of each sample on the basis of the gene expression and coefficient values, using the following formula:

where coefi represents the coefficient and x_i represents the expression value of each selected gene. We used the “survminer” package to obtain the best cut-off value for dividing the samples into a high-risk group and a low-risk group as well as the “survival” package in R to compute the survival curves of these two groups. A receiver operating characteristic (ROC) curve was generated using the “survival-ROC” package in R to obtain the area under the ROC curve. To highlight the superiority and accuracy of our model, we carefully compared our prognostic model with another three prognostic signatures (an autophagy-related long non-coding RNA signature constructed by Yu et al., an m6A-related lncRNA signature constructed by Yu et al., and a seven-MDEG signature constructed by Hu et al.) (24–26). All the genes used for the construction of prognostic signatures of Yu and HU were obtained. The “timeROC” R package was further applied to compute the area under the curve (AUC) values of each model. Finally, based on this model, we analyzed the correlation of the RS with the clinicopathological characteristics of patients with KIRC and used a heat map to depict it. There were no N-stage data because of the large amount of Nx data for these patients in the dataset of TCGA. A p < 0.05 was statistically significant.

The univariate and multivariate Cox regression analyses were used to determine the correlations of the patient age, tumor stage, tumor grade, tumor size (T), and tumor metastasis (M) with the RS in the model. Finally, the “rms” package in RStudio was used to draw a nomogram from the Cox analysis results and clinical characteristics for evaluating the survival probability of patients with KIRC.

The renal cancer cell line 786-O was purchased from the Institute of Cell Research, Chinese Academy of Sciences (China). The cells were routinely cultured in Roswell Park Memorial Institute (RPMI) 1640 medium, containing 10% fetal bovine serum and penicillin-streptomycin, at 37°C under a 5% CO₂ atmosphere. Once the culture had reached the logarithmic growth phase, 2 × 10⁵ cells were seeded into each well of 6-well plates. On the following day, 100 μl of a mixture containing Lipofectamine 3000 (Invitrogen, Carlsbad, CA, USA) and plasmid fragments carrying the tissue inhibitor of metalloproteinase 3 (TIMP3) gene diluted in serum-free medium was added to each well. After 6 h of incubation at 37°C, the medium was changed to serum-containing medium and the plates were further incubated for 24 h. Finally, the cells were digested with trypsin and collected for cell proliferation analysis using the CCK-8 assay (Dojindo, Japan). In brief, 1 × 10³ cells were seeded into each well of 96-well culture plates (four replicate wells per group). Subsequently, 10 μl of the CCK-8 reagent was added to each well according to the instructions from the manufacturer, and the plates were incubated in a 37°C incubator for 2 h. Finally, the absorbance of each well was measured at 450 nm using a microplate reader (EL340; BioTek Instruments, Hopkinton, MA, USA).

First, we draw a corresponding flowchart (Figure 1) to show the process of this whole research more clearly. In total, 24 angiogenesis-related genes were obtained from the GSEA website and referenced to TCGA datasets to obtain their CNVs and SNVs in 32 different types of cancer (Supplementary Figures 1A–C, Supplementary Tables 4–3). Although there were CNVs and SNVs in the angiogenesis pathway genes of most of the cancer types, there were almost no CNV gains or losses in these genes in thyroid cancer (THCA), thymoma, pancreatic adenocarcinoma, KIRC, and prostate adenocarcinoma. There was a high frequency of SNVs in uterine corpus endometrial carcinoma (UCEC), skin cutaneous melanoma, and colon adenocarcinoma. By contrast, the frequency of SNVs was low in kidney chromophobe (KICH), pheochromocytoma and paraganglioma, THCA, and uveal melanoma. Next, we used the expansion package in the R-language program to draw an image showing co-expression between the genes, whereupon we observed that there was a highly positive correlation among the genes FLT1, PDGFB, VEGFA, KDR, and TEK, whereas there was a negative correlation between SRC and FLT1. There was also a negative correlation between SRC and KDR and between SRC and TEK (Figure 2).

To analyze the expression changes in the angiogenesis-related genes in different cancer types, we used log2(fold change) to represent the ratio of the gene expression levels in the cancer tissues to those in the corresponding normal tissues. We found that most of the gene expression levels in the cancer tissues were different from those in the normal tissues (Figure 3A, Supplementary Table 4). Next, we constructed a survival landscape of these genes according to the correlation between the patient survival rate and gene expression in TCGA (Figure 3B, Supplementary Table 5). The gene is regarded as a protective gene when the hazard ratio (HR) is <1 and as a risk gene when the value is >1. Since the angiogenesis pathway itself is a cancer-related pathway, the angiogenesis-related genes were present as risk genes in most cancer types. Usually, the protective genes are downregulated and the risk genes are upregulated in tumors. Interestingly, we found five genes that were all risk factors in UCEC but were downregulated, which seemed contradictory. The following could be possible reasons: (1) tumors are heterogeneous, and the same genes are expressed differently in different tumors. (2) These risk genes are not the drivers of tumors and cannot directly lead to their occurrence and onset, and their expression may be regulated by other genes. (3) What we analyzed above was the mRNA expression of related genes; however, these genes ultimately play roles through their encoded proteins. Sometimes, the post-translational modifications, epigenetics, negative feedback, and other factors may lead to the mRNA and protein expression levels being inconsistent. We also observed from the pan-cancer analysis that most angiogenesis-related genes were protective genes in the patients with KIRC, an interesting phenomenon that was not evident in many other cancers, such as UCEC and KICH. Considering that an abundant blood supply, vigorous lipid metabolism, and insensitivity to radiotherapy and chemotherapy are all characteristics of KIRC, we focused on the relationship between the angiogenesis pathway genes and this disease in the following study. Using the “survminer” package, we divided the genes into the high-expression and low-expression groups based on the best cut-off values. Subsequently, we plotted Kaplan–Meier curves for identifying the statistically significant angiogenesis pathway genes in the patients with KIRC (Figure 3C). The results were consistent with the constructed survival landscape (Figure 3B).

To further explore the expression of angiogenesis-related genes in KIRC, we plotted their related heat maps (Figure 4A). We observed that the expression levels of most of these genes were significantly different between the tumor tissues and normal tissues. Next, a univariate Cox regression analysis of these genes in KIRC was carried out (Supplementary Figure 2, Supplementary Tables 9, 11), with the forest plot showing the HRs with 95% CIs and p-values. It was found that a high expression level of TEK, MAKP1, PDGFB, PIK3CA, CREBBP, ANGPY1, FLT1, PTK2, NOS3, TIMP3, KDR, and ARNT correlated with better survival rates in the patients with KIRC. By contrast, the high expression of FGF2, SRC, and PDGFRA correlated with the worse survival rates. Next, we clustered the angiogenesis-related genes obtained from the dataset of TCGA and divided the patients into three clusters based on the final angiogenesis score and gene expression level. C1 comprised the angiogenesis-inactive tumor tissues, C2 the angiogenesis-active tumor tissues, and C3 the tumor tissues with normal angiogenesis (Figure 4B, Supplementary Tables 6, 7). The violin plot clearly showed that the enrichment scores for the three clusters were in the order of C2 > C3 > C1 (Figure 4C). Next, we plotted the survival curves for the three clusters to determine whether the clustering was reasonable. The patients in C2 had significantly higher overall survival rates than the patients in C1 and C3 (Figure 4D), indicating that the angiogenesis score represented a protective factor. Finally, we explored the correlation between two of the clusters and the clinicopathological features (Figure 4E) and found that a higher angiogenesis score was negatively correlated with the tumor size (T), grade, and stage and tumor metastasis (M).

Considering that molecularly targeted therapy is currently a common way to treat KIRC, we obtained related data from the GDSC database to assess the response of these three angiogenesis clusters to 12 types of drugs. These drugs primarily include both commonly used targeted therapeutics, especially for kidney cancer, and classical drugs in tumor research, such as metformin. For the treatment of recurrent or metastatic RCC, the targeted drug therapy is preferred. This includes the use of tyrosine kinase inhibitors (e.g., sunitinib, sorafenib, axitinib, and pazopanib) and mTOR inhibitors (e.g., Texiromu) (27). Metformin, an activator of AMP-activated protein kinase (AMPK), is primarily used for the management of type 2 diabetes. The studies have shown that metformin exerts antitumor effects against different types of cancer (e.g., breast, colon, liver, prostate, and kidney cancer) (28–30) primarily by reducing glycemia to cutoff the phosphoinositide 3-kinase (PI3K)/mitogen-activated protein kinase (MAPK) pathway or by activating the AMPK pathway (31). Therefore, it is necessary to explore the correlation between these targeted drugs and the angiogenesis pathway. Hence, we conducted a drug sensitivity analysis and obtained estimated IC₅₀ values of the drugs for each sample. A lower IC₅₀ value is indicative of better drug sensitivity. The ridge regression model revealed the different drug sensitivities among the angiogenesis clusters to be as follows: pazopanib: C3 > C2 > C1; sorafenib: C1 > C2; sunitinib: C3 > C2; nilotinib: C2 > C1; vorinostat: C3 > C2 > C1; axitinib: C2 > C1; gefitinib: C1 > C2; temsirolimus: C3 > C1; lapatinib: C2 > C1; metformin: C1 > C2; bosutinib: C3 > C2; and tipifarnib: C1 > C2 (Figures 5A–L).

To further explore the differential expression patterns of the oncogenes and tumor suppressor genes in the three angiogenesis clusters, their related heat maps were plotted. We found that the expression levels of the oncogenes CCND1, BRAF, AKT1, MYC, KRAS, MTOR, PIK3A, and VEGFA were significantly higher in C2 than in C1. Additionally, the expression levels of the tumor suppressor genes von Hipel-Lindau (VHL), TP53, and PTEN were significantly lower in C1 than in C2 (Figure 6A). Among these tumor suppressor genes, the mutation of VHL can result in the overexpression of hypoxia-inducible factor-1 alpha (HIF-1α) protein, which is considered a hallmark of KIRC (32, 33). Continuously activated HIF-1α is shown to be related to cell proliferation, angiogenesis, and epithelial–mesenchymal transition, leading to the progression of KIRC and its metastasis to other organs (32, 34). The above results indicate that the poor prognosis of C1 may be correlated with the inhibition of tumor suppressor gene expression, which may play a more important role than the activation of oncogenes in C2. Additionally, the expression level of the oncogene HRAS was significantly higher in C1 than in C2.

In recent years, many studies have shown that SIRTs are involved in a variety of biological processes related to tumorigenesis, such as the changes in both the tumor-related metabolic pathways and tumor microenvironment and uncontrolled cell proliferation. In different cancer types and under different experimental conditions, SIRTs are thought to play a complex role as oncogenes or tumor suppressors (35). In our study, the level of SIRT1 expression was significantly higher in the angiogenesis-active group than in the angiogenesis-inactive group. By contrast, the levels of SIRT4, SIRT6, and SIRT7 expression were significantly lower in the angiogenesis-active group (Figure 6B). One study showed that Patrinia scabiosaefolia induced the death of 786-O cells via metabolic disruptions mediated by SIRT1 and mTOR signaling (36). Thus, SIRT1 inhibitors may be more effective for the patients who belong to the angiogenesis-active group. However, another study showed that SIRT1, SIRT3, and SIRT6 function as the tumor suppressors in RCC (37). In summary, these results indicate that SIRTs have a strong correlation with the angiogenesis signaling pathway and may act synergistically to promote or inhibit the multiple processes in the progression of KIRC.

Histone deacetylases, which catalyze the removal of acetyl groups from the histones and non-histone lysine residues, play an important role in gene transcription regulation (38) and are closely related to tumorigenesis and tumor metastasis. HDAC inhibition has recently become a clinically validated strategy for the treatment of cancer (39). We found that the expression levels of HDAC1, HDAC2, HDAC3, HDAC4, HDAC5, and HDAC9 were significantly higher in the angiogenesis-active group than in the angiogenesis-inactive group. By contrast, the expression levels of HDAC8, HDAC10, and HDAC11 were significantly lower in the angiogenesis-active group (Figure 6C). These results may offer new directions for the future precision treatment of tumors. For example, because HDAC2 was almost exclusively expressed in the angiogenesis-active group, the use of HDAC2 inhibitors may be more beneficial for patients with an active angiogenesis pathway.

The tumor microenvironment is composed of immune cells, tumor cells, stromal cells, and secreted cytokines and chemokines (40) and regulates the occurrence and progression of cancer (41). As an important component of the tumor microenvironment, the immune cells are closely related to the clinical outcome of cancer and are effective targets for anticancer treatment (42). In addition, angiogenesis plays a key role in regulating the tumor immune microenvironment (43). To investigate the correlation between the angiogenesis pathway and immunity in the patients with KIRC, we performed a correlation analysis between the angiogenesis pathway and immune cell infiltration (Figure 6D) and found that many genes related to the angiogenesis signaling pathway were associated with the immune cell infiltration, especially PDGFRA, MMP9, MAPK14, FGF2, and FGFR2. Among these, PDGFRA, MMP9, MAPK14, and FGF2 were positively correlated—whereas FGFR2 was negatively correlated—with immune cell infiltration. Next, we analyzed the correlation of the angiogenesis score with immune cell infiltration (Figures 6E–I), whereupon it was found to be positively correlated with the infiltration of type II IFN response cells, mast cells, and neutrophils, and negatively correlated with that of Tfh cells.

To explore the relationship between the genes in the angiogenesis pathway, co-expression analysis of the 24 angiogenesis-related genes was carried out, and the results are shown in Figure 7A. To determine whether the angiogenesis-related genes could be used to establish a model for predicting the clinical outcomes of patients with KIRC, the LASSO–Cox regression analysis of the 24 genes was performed, and 15 were finally selected to build the risk score model (Figures 7B,C, Supplementary Tables 10, 12). The patients with KIRC were divided into a high-risk group and a low-risk group according to their RS values. The overall survival rates of patients in the high-risk group were significantly lower than those of patients in the low-risk group (Figure 7D, Supplementary Table 12). To explore the prognostic prediction efficiency of the new survival model in these same patients, a ROC curve analysis was performed. The areas under the ROC curves of the survival model for predicting 3-, 5-, 7-, and 10-year survival rates were 0.734, 0.752, 0.763, and 0.787, respectively, indicating that the risk model had a high predictive value (Figures 7E–H). Survival probability predicted by our angiogenesis prognostic signature was superior to the m6A-related lncRNA signature constructed by Yu et al. (25) and the seven-MDEG signature constructed by Hu et al. (26). For the patients with KIRC whose survival time was 1–4-years, our prognostic signature also showed a higher predictive accuracy compared with the autophagy-related long non-coding RNA signature constructed by Yu et al. (24). Additionally, for the patients whose survival time was 4–8-years, our prognostic signature showed similar predictive accuracy compared with the autophagy-related long non-coding RNA signature (24) (Supplementary Figure 3). A further statistical test was performed on the differences between the risk subgroups and a heat map was drawn to visualize the correlation between the RS and clinical data (Figure 7I). We found that the risk model was correlated to tumor metastasis (M) and the tumor size (T), stage, grade, and fustat, with patients in the high-risk group tending to have an advanced histological grade and to be in an advanced clinical stage.

First, we performed a univariate Cox regression analysis on the RS and other clinicopathological features of the patients with KIRC (Figure 8A, Supplementary Table 13). The forest plot showed that the patient age, tumor grade, stage, and size (T), tumor metastasis (M), and RS correlated with the overall survival of the patients. The multivariate Cox regression analysis (Figure 8B, Supplementary Table 14) revealed that the patient age, tumor grade and stage, and RS were independent risk factors correlated with the overall survival. On the nomogram based on the risk model (Figure 8C), the second to ninth lines represent the patient age, tumor grade, tumor stage, RS, total points, 5-, 7-, and 10-year survival, respectively. The total score in the sixth row is the sum of the scores for each item from the second to fifth lines. The 5-, 7-, and 10-year survival rates were predicated based on the total score. For example, if the total score is 100, then the 5-year survival rate is approximately 0.3.

The viability of TIMP3-overexpressing 786-O renal cancer cells was determined using the CCK-8 assay, which showed that the proliferation of these cells was significantly inhibited (Figure 9). These results were consistent with our previous bioinformatics analysis, which indicated that TIMP3 was a protective gene.