We identified a total of 372 differentially expressed AF genes in cancer and normal samples (with a threshold of adj.P.Val<0.01 & |log (FC) |≥1) (Figures 1A,B). GO and KEGG functional enrichment analyses of the differentially expressed AF genes were then performed (Figures 1C–F). The enriched GO terms of DEGs were classified into three categories: molecular functions, cellular components, and biological processes. The results revealed that these genes were enriched for GO terms such as regulation of vasculature development, regulation of angiogenesis, ameboidal-type cell migration, and positive regulation of vasculature development, epithelial cell proliferation, and tissue migration. The KEGG pathway enrichment showed the enrichment of critical pathways involved in tumorigenesis and metastasis, including pathways in cancer, focal adhesion, the MAPK signaling pathway, the chemokine signaling pathway, the TGF-β signaling pathway, and renal cell carcinoma. The top 15 highly enriched KEGG pathways are presented in Figure 1F.

We performed a univariate Cox regression analysis of these differentially expressed AF genes and identified 58 AF-related genes that were associated with the prognosis of NSCLC. We performed survival analyses of the top five genes in terms of the p-value (Figures 2A–E). The low expression of these five genes was associated with a worse prognosis (Figure 2F).

A total of five AF genes significantly associated with prognosis (p < 0.001) in the univariate Cox regression were further selected for lasso regression (Supplementary Figure S2). We first used cross-validation to identify the minimal lambda, i.e. lambda min, and then selected the four most significant genes using lambda min to develop the prognostic risk model. The optimized model was: risk score = 0.104 * SLC2A1 + 0.138 * FSTL3 + 0.069 *C1QTNF6 + 0.046 * PTX3. We calculated the risk scores of each sample using this formula and classified all the samples into high- and low-risk groups according to the median for further analysis.

To validate the performance of our model, we plotted the Kaplan–Meier survival curves of the high- and low-risk groups (Figure 3A). A significant association was shown between the risk group and survival (p < 0.0001), suggesting that the model had a high prognostic value. Time-dependent ROC curves were further plotted, which showed AUC>0.6 in the 1-year, 3-years, and 5-years ROC curves. This indicated that the model had good prediction ability (Figure 3B). Based on the optimistic cutoff, the patients were divided into high AF risk score and low AF risk score groups (Figures 3C–F).

We used the GSE4573 and GSE68465 datasets to validate the model (Figures 4A–D). We combined the two datasets and removed the batch effect. We selected the prognostic genes in the datasets (C1QTNF6 was not identified) and calculated the risk score using the coefficients in the model for validation. The Kaplan–Meier plot showed that the samples in the high-risk group had a worse prognosis with a p-value < 0.05, which indicated that our model had high accuracy.

We analyzed the difference in AF risk scores of each group that was stratified by clinical characteristics. The risk score of LUSC was significantly higher than that of LUAD (Figure 5A). The risk score of the samples with EGFR mutations was significantly lower than that of samples without EGFR mutations (Figure 5B). The risk score also differed significantly across the different tumor stages and TNM stages, which was consistent with the process of carcinogenesis (Figures 5C–F). The patients with a smoking history also had significantly higher risk scores than those who never smoked (Figure 5G).

We also visualized the association of the risk score with tumor mutational burden (TMB), homologous recombination deficiency (HRD), neoantigen burden, chromosomal instability (CIN), and stemness index (mRNAsi) (Figures 6A–D). TMB is a marker for genomic instability measured by sequencing the whole tumor genome and has been shown to correlate with immunotherapy (Gibney et al., 2016). Therefore, TMB is emerging as a predictor of immunotherapeutic responses. For all indexes, the highest correlation was obtained for TMB (Figure 6A). This further illustrates that the interaction of AFs can affect immunotherapy. The discovery of homologous recombination deficiency (HRD) in lung cancer is of great importance for patients who will benefit from poly ADP-ribose polymerase inhibitor (PARPi) (Weston et al., 2010). However, we did not find a correlation between HRD and AFs (Figure 2B). Neoantigens are another important index for predicting the clinical response to immunotherapy. The current studies of neoantigen sources mainly focus on single nucleotide variants (SNVs), such as small insertions and deletions (indels), somatic copy number variations (SCNVs), and large scale transition (LSTms). Similarly, we found no significant differences in these parameters (Figures 6C–H). The stemness index (mRNAsi) is used to measure the tumor development and evaluate the reliability of stem cell indexes as shown in Figures 6A,I significant positive correlation was found between AFs and mRNAsi. These results confirmed that AFs were related to biological processes, cancer metastasis, and the immune microenvironment.

The immune infiltration status was highly associated with the prognosis of NSCLC. We used the CIBERSORT algorithm to calculate and compare the proportion of immune infiltration in the high- and low-risk groups based on TCGA data (Figure 7A). The proportions of naive B cells, memory activated CD4 T cells, gammadelta T cells, and resting dendritic cells were significantly increased in the low-risk group, while the proportions of memory B cells, and macrophages. M0, macrophages. M2, and activated mast cells was significantly higher in the high-risk group, which indicated that the immune infiltration status was different in the high- and low-risk groups.

We also found that the stroma score (p = 7.8e-16), immune score (p = 0.012), and tumor purity (p = 1.7e-08) were significantly higher in the high-risk group than in the low-risk group (Figures 7B–D).

We further investigated the difference in the mutation profiles between the high- and low-risk groups based on TCGA data. The mutation rate of the high-risk group was slightly higher than that of the low-risk group (92.81 vs. 90.91%). The mutation rate of TP53 was the highest in both the high- and low-risk groups. Additionally, missense mutations were the most dominant among all mutation types. Single nucleotide polymorphisms (SNPs) occurred more frequently in the high-risk group than in the low-risk group (Figures 8A,B).

We also investigated the difference in CNV between the high- and low-risk groups (Figures 8C–E). The copy numbers of amplification and deletion were distributed differently in the same position. Significant differences in distribution could be observed in the figures (Figures 8C,D). We analyzed the Z-score of the high- and low-risk groups (Figure 8E) by t-test. The results showed a significant difference between them (p < 2.22e-16).

Immunotherapy offers a new approach to cancer treatment. For a long period of time, immunotherapy approaches targeting PD1, PDL1, and ctla-4 have all been successfully applied in cancer, with largely promising outcomes. Tumor immune dysfunction and exclusion (TIDE) is a gene expression biomarker developed for predicting the clinical response to immune checkpoint blockade. We used the TIDE score to assess the performance of the risk score to predict the response to immunotherapy and visualized it in R software. A significant difference in the TIDE score was demonstrated between the high- and low-risk groups (p = 0.0027) (Figure 9A), while its prediction accuracy was lower than that of the risk score (Figure 9B).

To assess the effect of the risk score on prognosis, we performed univariate and multivariate Cox regression analyses of the above clinical characteristics and validated the risk model using validation datasets (Figures 10A–D). The risk score showed a significant effect on the prognosis in both univariate and multivariate regression analyses.

Finally, we developed nomograms using the risk score and clinical characteristics and validated them with calibration plots (Figure 11A). The risk score showed the highest accuracy of prediction (Figure 11B). The 1-year, 2-years, and 3-years calibration plots demonstrated the highest accuracy of our nomograms (Figures 11C–E).

The expression profile combined with patients’ clinical and annotation information in LUAD and LUSC datasets were downloaded from UCSC (https://xenabrowser.net/). Next, we averaged the expression level of genes with the same name and removed the genes with expression levels lower than 30%. We merged the two expression profiles after processing and converted the data type from FPKM to TPM. The samples from patients aged >18 years were extracted, and batch effects were removed (Figures 12A,B). We then searched the NCBI database (https://www.ncbi.nlm.nih.gov/gene/?term=angiogenic) using “angiogenic factor (AF)” as the keyword and extracted AF expression data of 1,054 samples from the downloaded expression profile.

We used the R package “limma” to identify AF-related differentially expressed genes (DEGs) (threshold: adj.P.Val<0.01 & |log (FC) |≥ 1) in 372 cancer and normal samples. Next, gene ontology (GO) enrichment analysis (p-value cutoff < 0.05) and KEGG pathway enrichment analysis (p-value cutoff < 0.05) of differentially expressed genes were performed using the R package “clusterProfiler”.

Other data of cancer samples were further extracted, and a univariate Cox regression analysis of DEGs associated with overall survival was performed using the R packages “survival” and “survminer” with a threshold of p < 0.05. DEGs associated with prognosis were identified after screening.

Lasso regression was performed using the R package “glmnet” for downscaling prognostic genes. We first screened lambda by cross-validation, and then selected the model with lambda. min. Next, the expression matrix of the selected genes for the model was extracted, and the risk score of each sample was calculated using the following formula: Riskscore i = ∑ i = 1 n exp ji ∗ β j .

It represented the expression level of gene j in sample i, and represented the coefficient of gene j in the lasso regression model. All the samples were stratified into high- and low-risk groups according to the median-risk score.

Kaplan–Meier survival curves were plotted according to high- or low-risk groups. The ROC curves were drawn based on the predicted risk score of each sample.

The samples with AF risk scores were grouped according to clinical characteristics. We used the R package “ggplot2” to show the distribution of AF risk scores in each group and the R package “ggpubr” to illustrate the significant difference between groups.

We calculated tumor mutational burden, homologous recombination deficiency (HRD) (from technical support), tumor neoantigen burden (according to the literature The Immune Landscape of Cancer), chromosome instability (CIN) (according to the literature The Immune Landscape of Cancer), and stemness index (according to the literature Machine Learning Identifies Stemness Features Associated with Oncogenic Dedifferentiation) based on AF risk scores and performed association analyses.

The proportion of 22 immune cells in the samples can be derived using the CIBERSORT algorithm based on the expression of certain genes. We assessed the difference in immune infiltration between the high- and low-risk groups by t-test with a significance threshold of p < 0.05.

We analyzed the differences in the immune score, stromal score, and tumor purity of AF in high- and low-risk groups using the R package “ESTIMATE” and assessed the differences in immune infiltration in the high- and low-risk groups by t-test with a threshold of p < 0.05.

MAF files of NSCLC were downloaded from the GDC database, and we extracted the mutation information of AF from the somatic mutation profile. The mutation profile of AF in high- and low-risk groups was demonstrated with the help of the “oncoplot” function, using the R package “maftools”.

The copy number variation (CNV) data of LUSC and LUAD were downloaded from UCSC. The copy numbers of the high- and low-risk groups were extracted to generate files and plotted on the gadget using the “CNV distribution chart - bar graph section”.

The expression profiles of immune genes were extracted from the processed data of endometrial cancer samples, and the immune gene sets were obtained from the ImmPort database (https://www.immport.org/) and InnateDB database (https://www.innatedb.ca/). The expression profiles of the immune gene sets were subsequently normalized. The predicted TIDE scores of the samples were calculated using the TIDE online database. The distribution of TIDE scores in the high- and low-risk groups was illustrated with box plots using ggpubr, and the significance was tested by t-test.

To validate whether the risk score was an independent prognostic factor, univariate Cox regression analyses of the candidate prognostic factors using TCGA sample data were first performed, including risk score, age, tumor stage, and TNM stage. A multivariate Cox regression analysis was subsequently performed to assess the effect size of the risk score. We used the function cph in the R package “rms” to plot the nomograms and calibration plots for visualization.

All statistical analyses were performed using R software version 4.0.3. p < 0.05 was set as the significance criterion.