Transcriptome and clinical data from 33 cancers were downloaded from TCGA database (https://tcgadata.nci.nih.gov/tcga/) using UCSC Xena (https://xena.ucsc.edu/), a genome browser for visualizing gene and variant information. The 33 tumor datasets from the TCGA database were obtained by clicking on the “Launch Xena” and “DATA SETS” options on the UCEC website. Gene expression and clinical data from any tumor can be downloaded in the “gene expression RNAseq” and “phenotype” sections of the tumor datasets, respectively. We further downloaded the RNA sequencing data of 31 different tissues from the GTEx biobank (https://commonfund.nih.gov/GTEx) to obtain more normal tissues. We clicked on “Downloads”—“Open Access Data” from the GTEx database homepage. GTEx Analysis V8 (dbGaP Accession phs000424.v8.p2) version was chosen. The annotation file and gene expression data were downloaded in the “Annotations” and “RNA-Seq Data” sections, respectively. The data from TCGA and GTEx was merged. The RNA sequencing data was transformed into transcripts per million (TPM) reads and normalized by log2 transformation for the following analysis.

An analysis of ANXA2 expression in tumor and normal tissues was performed using R statistical software (version 4.0.3). Box plots were created using the R package “ggplot2”. The HPA database (https://www.proteinatlas.org/) was used to investigate the gene expression of ANXA2 in 27 types of normal tissues, the protein expression of ANXA2 in 20 cancers and the protein location of ANXA2 in cells. We downloaded immunohistochemical images of colon adenocarcinoma (COAD), glioblastoma multiforme (GBM), liver hepatocellular carcinoma (LIHC), prostate adenocarcinoma and stomach adenocarcinoma (STAD) with the corresponding normal tissues from HPA to evaluate the differential expression of ANXA2. The GEPIA database (http://gepia.cancer-pku.cn/) was used to verify the differential expression of ANXA2 in the aforementioned tumors and corresponding normal tissues.

The receiver operating characteristic (ROC) curve was calculated to evaluate the diagnostic ability of ANXA2 by the R software (version 4.0.3). The curves were drawn using the R packages “pROC” and “ggplot2”. The area under the curve (AUC) > .8 was considered of adequate diagnostic value (Mandrekar, 2015). The survival analysis was conducted to analyze the prognosis of differential expression of ANXA2 in 33 types of cancer by using the packages “survival” and “survminer”. The overall survival (OS), disease-specific survival (DSS) and progression-free interval (PFI) were used to evaluate the prognostic value of ANXA2. The Kaplan-Meier method (Kaplan-Meier curve plus univariate Cox regression analysis) was used for survival analysis. Also, a pan-cancer survival analysis of ANXA2 gene was performed using the GEPIA database to explore the prognostic value of ANXA2 in all 33 tumors from the TCGA database.

The correlation between ANXA2 expression and tumor staging was tested with the R software (version 4.0.3). Box diagrams were drawn using the R package “ggplot2”.

Tumor types with statistical significance in both survival analysis and tumor staging correlation analysis were chosen to draw the nomograms using R statistical software (version 4.0.3), so as to highlight the value of ANXA2 expression in cancer prognosis and reduce the interference of other clinical factors. R package “rms” and “survival” were used to construct nomograms and calibration curves which were used to verify the accuracy of the nomograms.

The RNA sequencing expression profile was used to investigate the infiltration of 24 types of immune cells in tumors with the R software (version 4.0.3). The correlation between ANXA2 expression and immune cell infiltration in 33 types of cancer was analyzed by the R package “GSVA”. The correlation between ANXA2 expression and StromalScore, ImmuneScore and ESTIMATEScore was analyzed by the R package “estimate”, |r| > .3 was considered statistically relevant. In addition, we conducted a co-expression analysis to explore the correlation between expression of ANXA2 and 47 immune checkpoints and immunomodulators (including immune-activating genes, immune-suppressing genes, chemokine ligands, chemokine receptors and MHC genes) by using the R package “ggplot2”. The analyses were performed by using the R software (version 4.0.3).

The correlation analysis was performed to investigate the relationship between ANXA2 expression and DNA methylation, TMB, MSI and MMR by the R software (version 4.0.3). The package “ggplot2” was used for visualizing data. The DNA methylation analysis was based on Illumina methylation 450 data and cg02395965 probe. The cBioPortal website (http://www.cbioportal.org/) was applied to analyze the genetic characteristics of ANXA2 in 30 types of cancer. The genetic variation characteristics of ANXA2 included mutation, structural variants, amplification, deep deletion and multiple alterations in pan-cancer.

The Gene Set Enrichment Analysis (GSEA) was performed to predict the ANXA2-related signaling pathways by using the R software (version 4.0.3). The package “clusterProfiler” was used to obtain the enrichment maps. The adjusted p-value (<.05), normalized enrichment score (|NES| > 1), and False Discovery Rate (FDR, q value < .25) were used as criteria to determine a statistically significant phenotype.

The differential expression of ANXA2 was estimated using the Wilcoxon rank sum testand Wilcoxon signed rank test. A log-rank test was used to perform the survival analysis of ANXA2. The Kruskal–Wallis test was used to conduct the correlation analysis between ANXA2 expression and tumor TNM staging. The Spearman’s rank correlation coefficient was used to estimate the correlation between ANXA2 expression and infiltration of 24 types of immune cells, expression of 47 immune checkpoints, StromalScore, ImmuneScore, ESTIMATEScore, immunomodulators, DNA methylation, TMB, MSI and MMR. A p-value < .05 was considered statistically significant (*p < .05, **p < .01, ***p < .001).

In this study, we aimed at analyzing the expression of ANXA2 in various tumors. The sample sizes for each comparison were added to Supplementary Table S1. Based on the TCGA database, the expression of ANXA2 in cervical squamous cell carcinoma (CESC), CHOL, COAD, esophageal carcinoma (ESCA), GBM, head and neck squamous cell carcinoma (HNSC), kidney renal clear cell carcinoma (KIRC), kidney renal papillary cell carcinoma (KIRP), LIHC, STAD, thyroid carcinoma (THCA) and uterine corpus endometrial carcinoma (UCEC) was higher compared with corresponding normal tissues. The expression of ANXA2 in breast invasive carcinoma (BRCA), KICH, lung adenocarcinoma (LUAD) and PRAD was lower compared with corresponding normal tissues (Figure 1A). Due to the shortage of normal samples in the TCGA database, we combined the data of TCGA and GTEx databases for expression analysis. Compared with normal tissues, ANXA2 was upregulated in bladder urothelial carcinoma (BLCA), BRCA, CESC, CHOL, COAD, lymphoid neoplasm diffuse large B cell lymphoma (DLBC), ESCA, GBM, HNSC, KIRC, KIRP, brain lower grade glioma (LGG), LIHC, lung squamous cell carcinoma (LUSC), ovarian serous cystadenocarcinoma (OV), pancreatic adenocarcinoma (PAAD), rectum adenocarcinoma (READ), STAD, testicular germ cell tumor (TGCT), THCA, thymoma (THYM), UCEC and uterine carcinosarcoma (UCS). ANXA2 was downregulated in adrenocortical carcinoma (ACC), KICH, acute myeloid leukemia (LAML), PRAD and skin cutaneous melanoma (SKCM) (Figure 1B). In paired samples from TCGA database, ANXA2 was upregulated in CHOL, ESCA, HNSC, KIRC, KIRP, LIHC, STAD, THCA compared with normal tissues. ANXA2 was downregulated in KICH, LUAD, PRAD compared with normal tissues (Figure 1C). Based on the TCGA database alone, the expression of ANXA2 was lower in BRCA compared with corresponding normal tissues, the average expression levels of ANXA2 were 9.675 and 9.887 in tumor and normal tissues, respectively. When TCGA and GTEx databases were combined together, the expression of ANXA2 was higher in BRCA compared with corresponding normal tissues, the average expression levels of ANXA2 were 9.675 and 9.545 in tumor and normal tissues, respectively.

We further investigated the gene and protein expression levels of ANXA2 through the HPA database. Among 27 normal tissues, protein expression level of ANXA2 was highest in the esophagus (Supplementary Figure S1A). Among 20 different types of cancer, the protein level of ANXA2 was highest in renal cell carcinoma (Supplementary Figure S1B). We speculated that ANXA2 is more likely to exert significant physiological or pathological functions in normal tissues and cancers with high levels of ANXA2 protein expression. The ANXA2 protein was mainly located in the plasma membrane and cytoplasmic matrix (Supplementary Figure S1C). To evaluate the expression of ANXA2 in more detail, we analyzed the immunohistochemistry results obtained through the HPA database. In addition, we compared the results with the ANXA2 gene expression from the GEPIA database. The findings were consistent between the two databases (Figures 2A–E). Normal brain, liver and stomach tissues showed a weak or no ANXA2 staining, whereas GBM, LIHC and STAD tissues showed a moderate staining. Normal colon tissues showed a moderate staining, while COAD tissues showed a strong staining. Normal prostate tissues showed a strong staining, while PRAD tissues showed a weak or no staining.

The ROC curve was used to evaluate the diagnostic efficacy of ANXA2 in 33 types of cancer. ANXA2 showed a high diagnostic efficiency (AUC > .8) in the following 12 types of cancer: CHOL (AUC = .997), GBM (AUC = .986), KICH (AUC = .864), KIRP (AUC = .874), LIHC (AUC = .893), OV (AUC = .917), PAAD (AUC = .972), PRAD (AUC = .828), READ (AUC = .931), STAD (AUC = .900), TGCT (AUC = .880) and THYM (AUC = .883) (Figures 3A–L).

A survival analysis was conducted to examine the association between ANXA2 expression and prognosis of 33 types of cancer. Based on the Cox proportional hazards model, the ANXA2 expression level were related to OS of BLCA (p = .013), CESC (p = .014), HNSC (p = .007), LGG (p < .001), LIHC (p = .021), LUAD (p < .001), mesothelioma (MESO) (p < .001), OV (p = .008), PAAD (p < .001) and uveal melanoma (UVM) (p = .004) (Figure 1D). The Kaplan-Meier survival analysis showed that the high expression of ANXA2 was associated with a short OS in the aforementioned tumors (Figure 4A). The expression level of ANXA2 correlated with DSS of BLCA (p = .016), CESC (p = .031), HNSC (p = .006), KIRC (p = .028), LGG (p < .001), LUAD (p = .003), MESO (p < .001), OV (p = .011), PAAD (p < .001) and UVM (p = .003) (Figure 1E). The Kaplan-Meier survival analysis showed that the high expression of ANXA2 was associated with a poor prognosis in the ten types of cancer described above (Figure 4B).

According to the forest plot, ANXA2 expression was related to PFI in ACC (p = .043), BLCA (p = .018), GBM (p = .044), HNSC (p = .014), KIRC (p = .008), KIRP (p = .027), LGG (p < .001), LIHC (p = .022), LUAD (p = .019), MESO (p = .009), PAAD (p < .001), UCEC (p < .001) and UVM (p = .002) (Figure 1F). The Kaplan-Meier survival analysis showed that the high expression of ANXA2 was associated with a poor PFI in ACC, BLCA, GBM, HNSC, KIRC, KIRP, LGG, LIHC, LUAD, MESO, PAAD and UVM. The low expression of ANXA2 was associated with a poor PFI in UCEC (Figure 4C). The pan-cancer survival analysis of ANXA2 gene based on the GEPIA database showed that high expression of ANXA2 was associated with lower overall survival in all 33 types of cancer (p = 0) (Supplementary Figure S2).

The gene expression of ANXA2 was related to the clinical T stage of ACC, BRCA, HNSC, KIRC and THCA (Figures 5A–E), the clinical M stage of BLCA, CESC, CHOL, LUAD and THCA (Figures 5F–J) and the clinical N stage of ACC, CESC, KIRP and LUSC (Figures 5K–N) in 33 types of cancer.

BLCA, HNSC and LUAD were selected for the construction of nomograms. The clinical factors included in the nomogram for BLCA were: ANXA2 expression, age, clinical T, M, N stage and gender (Supplementary Figure S3A). The clinical factors included in the nomogram for HNSC were ANXA2 expression, age, clinical T, N stage, gender and radiation therapy (Supplementary Figure S3C). The clinical factors included in the nomogram for LUAD were ANXA2 expression, age, clinical T, N stage, gender and race (Supplementary Figure S3E). Among the three nomograms, ANXA2 showed the best prognostic prediction ability, and the calibration curves confirmed the accuracy of the nomograms (Supplementary Figures S3B, D, F). The nomograms further showed the good prognostic value of ANXA2 in different tumor types.

The Spearman’s correlation coefficient was used to indicate the association between the gene expression of ANXA2 and immune infiltration in 33 types of cancer. The association between ANXA2 expression and infiltration of six important immune cells (B cells, dendritic cells, macrophages, neutrophils, NK cells and T cells) was analyzed firstly. ANXA2 expression was positively correlated with the immune cells mentioned before in pheochromocytoma and paraganglioma (PCPG) (Figure 6A) and PRAD (Figure 6B) (all r > .3, all p < .05). We also found that the ANXA2 expression positively correlated with the StromalScore of 12 types of cancer (Supplementary Figure S4A) (all r > .3, all p < .05), the ImmuneScore of 11 types of cancer (Supplementary Figure S4B) (all r > .3, all p < .05) and the ESTIMATEScore of 11 types of cancer (Figure 6C) (all r > .3, all p < .05). We further explored the association between the ANXA2 gene expression and infiltration of 24 types of immune cells in pan-cancer. ANXA2 gene expression was associated with varying levels of immune infiltration in 32 types of cancer (except CHOL) (Figure 6D) according to the findings. We explored the correlation between ANXA2 expression and 47 immune checkpoints in pan-cancer. ANXA2 expression was variably related to the immune checkpoints in all 33 types of cancer, with the most relevant correlation being with CD276 (31/33) (Figure 6E). CHOL was associated with two immune checkpoints, CD276 and TNFSF18.

We conducted a gene co-expression analysis in 33 types of cancer exploring the correlation between ANXA2 expression and immunomodulators such as immune-activating genes, immune-suppressing genes, chemokine ligands, chemokine receptors and MHC genes (Figures 7A–E). ANXA2 was variably associated with immunomodulators studied in all the tumors (p < .05) and positively associated with most immunomodulators in 30 types of cancer (except CESC, HNSC and LUSC).

We investigated the correlation between the ANXA2 expression and DNA methylation in pan-cancer (Table 1). The expression of ANXA2 positively correlated with DNA methylation in BLCA (r = .162, p = .001), BRCA (r = .075, p = .035), CESC (r = .135, p = .019), ESCA (r = .257, p = .001), HNSC (r = .146, p = .001), LUSC (r = .114, p = .029), SKCM (r = .102, p = .027), TGCT (r = .567, p < .001), UCEC (r = .210, p < .001) and UCS (r = .363, p = .006) and negatively correlated with DNA methylation in sarcoma (SARC) (r = −.198, p = .001), STAD (r = −.244, p < .001) and THYM (r = −.355, p < .001) (Figure 8A).

We further investigated the correlation between ANXA2 expression and TMB, MSI and MMR in 33 types of cancer. The ANXA2 expression positively correlated with TMB in eight types of cancer (all r > 0): COAD (p < .001), LGG (p < .001), PAAD (p < .001), SARC (p < .001), SKCM (p = .002), STAD (p < .001), THYM (p = .0017) and UCEC (p < .001). The ANXA2 expression negatively correlated with TMB in five types of cancer (all r < 0): CESC (p = .022), LAML (p = .0087), LUAD (p = .0014), LUSC (p = .0013), and PRAD (p < .001) (Figure 8B). The expression of ANXA2 positively correlated with MSI in seven types of cancer (all r > 0): COAD (p < .001), KIRC (p = .032), READ (p = .018), SARC (p < .001), STAD (p < .001), TGCT (p = .0095) and UCEC (p < .001). The ANXA2 expression negatively correlated with MSI in three types of cancer (all r < 0): GBM (p = .027), LUAD (p = .0012) and PRAD (p = .022) (Figure 8C). The ANXA2 expression was variably associated with five MMR genes (MLH1, MSH2, MSH6, PMS2 and EPCAM) in 32 types of cancer and was closely correlated with MMR in BLCA, KIRP and LIHC (Figure 8D). ANXA2 showed the highest frequency of variation in UCEC (2.9%). Gene mutations were the most common type of genetic variations. The highest mutation frequency of ANXA2 was reported in SKCM (2.48%) (Figure 8E).

The GSEA was conducted to further explore the ANXA2-related pathways in tumorigenesis and tumor development. ANXA2 was associated with immune-related pathways in CESC, CHOL, ESCA, HNSC, KICH, LAML, LUAD, LUSC, MESO, PAAD, PCPG, PRAD, READ, SKCM and UVM. In CESC, ANXA2-related immune pathway were adaptive immune response based on somatic recombination of immune receptors bulit from immunoglobulin superfamily domains and humoral immune response (Figure 9A). In CHOL, ANXA2-related immune pathways were regulation of lymphocyte activation, immune response regulating signaling pathway and regulation of immune effector process (Figure 9B). In ESCA, ANXA2-related immune pathway was complement activation (Figure 9C). In HNSC, ANXA2-related immune pathway were adaptive immune response based on somatic recombination of immune receptors bulit from immunoglobulin superfamily domains, antigen binding, antigen receptor mediated signaling pathway, B cell activation and B cell mediated immunity (Figure 9D). In KICH, LAML and PRAD, ANXA2-related immune pathway were regulation of lymphocyte activation and immune response regulating signaling pathway (Figures 9E,F,L). In LUAD, the ANXA2-related immune pathway were humoral immune response mediated by circulating immunoglobulin and immunoglobulin complex (Figure 9G). In LUSC, ANXA2-related immune pathway were adaptive immune response based on somatic recombination of immune receptors bulit from immunoglobulin superfamily domains, B cell activation and B cell mediated immunity (Figure 9H). In MESO, the ANXA2-related immune pathway was the humoral immune response (Figure 9I). In PAAD, ANXA2-related immune pathway was immune response regulating signaling pathway (Figure 9J). In PCPG, ANXA2-related immune pathway was adaptive immune response based on somatic recombination of immune receptors bulit from immunoglobulin superfamily domains (Figure 9K). In READ, ANXA2 related immune pathway were humoral immune response, adaptive immune response based on somatic recombination of immune receptors bulit from immunoglobulin superfamily domains and lymphocyte mediated immunity (Figure 9M). In SKCM, ANXA2-related immune pathways were immune response regulating signaling pathway, humoral immune response, adaptive immune response based on somatic recombination of immune receptors bulit from immunoglobulin superfamily domains, humoral immune response and lymphocyte mediated immunity (Figure 9N). In UVM, the ANXA2-related immune pathway was leukocyte migration (Figure 9O). The above-mentioned immune-related pathways were not observed in other tumors (Supplementary Figure S5).