Bioinformatics analysis of the clinical value and potential mechanisms of AHNAK2 in papillary thyroid carcinoma

Diagnostic value of AHNAK2

The mRNA expression level of AHNAK2 was significantly higher in PTC tissues than in normal tissues (P<0.0001) (Figure 1A). By performing receiver operating characteristic (ROC) curve analysis to discriminate PTC from normal thyroid tissues, we found that the area under the curve (AUC) of the AHNAK2 expression level was 0.88, suggesting that it could be a good diagnostic biomarker (Figure 1B).

Figure 1. Diagnostic and prognostic value of AHNAK2 in PTC. (A) Comparison of AHNAK2 mRNA expression levels between PTC and normal tissues. (B) Diagnostic efficacy of the ROC curve of AHNAK2. (C) Comparison of immunohistochemistry images of AHNAK2 between PTC and normal thyroid tissues with four different AHNAK2 antibodies (HPA000878, HPA004145, HPA002940 and HPA020111) based on the Human Protein Atlas. (D) K-M survival analysis was performed to determine differences in PFS between the H-AHNAK2 and L-AHNAK2 groups.

Additionally, immunohistochemistry performed with four different AHNAK2 antibodies showed deeper staining levels in PTC tissues than in normal thyroid tissues, suggesting higher protein expression in PTC (Figure 1C).

Prognostic value of AHNAK2

Kaplan-Meier (K-M) survival analysis showed that the high AHNAK2 (H-AHNAK2) group was associated with worse progression-free survival (PFS) than the low AHNAK2 (L-AHNAK2) group (P=0.0005) (Figure 1D). Multivariate Cox regression analysis suggested that AHNAK2 could become an independent predictor after adjusting for other parameters, including age, gender, American Joint Committee on Cancer (AJCC) stage, BRAF mutations, RAS mutations and tumor mutation burden (TMB) (Table 1) (PFS as the ending indicator). In summary, AHNAK2 could be an important tool for distinguishing PTC patient clinical outcomes.

Clinical value of AHNAK2

We analyzed the relationships between AHNAK2 and clinical parameters, including age, gender, stage, metastasis, N classification, T classification, pathologic type, and mutations in BRAF and RAS (Table 2, Supplementary Figure 1). AHNAK2 expression was significantly high in patients with an advanced stage, an advanced T classification, lymph node metastasis, increased BRAF mutations, decreased RAS mutations, and lower follicular PTC and higher tall-cell PTC proportions.

Interestingly, in the study of the relationship between AHNAK2 expression and gene mutations in PTC (Supplementary Figure 1A), we found that AHNAK2 expression was positively correlated with BRAF mutations (Supplementary Figure 1B) and negatively correlated with RAS (NRAS, HRAS and KRAS) mutations (Supplementary Figure 1C–1E).

Screening modules and genes related to AHNAK2

To construct a weighted coexpression network and screen for modules and genes related to AHNAK2, 2696 differentially expressed genes (DEGs) between PTC and normal tissues from The Cancer Genome Atlas (TCGA) were selected and subjected to weighted gene coexpression network analysis (WGCNA) (475 patients with complete clinical information were selected) (Figure 2A and 2B). After a series of adjustments for WGCNA parameters, the DEGs were divided into 11 modules by average linkage hierarchical clustering (Figure 2C–2F). The blue module, which contains 407 genes, exhibited the highest correlation with AHNAK2 expression (Pearson’s correlation coefficient =0.65, P<0.0001) (Figure 2G). Fifty-seven genes in the blue module were selected as hub genes (absolute module membership [MM] > 0.5 and absolute gene significance [GS] > 0.5) (Figure 2H).

Figure 2. Screening for modules and genes related to AHNAK2 in PTC. (A) Volcano plot of DEGs between PTC and adjacent tissues in TCGA data. (B) Heatmap of DEGs between PTC and adjacent tissues in TCGA data. (C) Calculation of the scale-free fit index of various soft-thresholding powers (β). (D) Analysis of the mean connectivity of various soft-thresholding powers (β). (E) Clustering dendrogram of 475 PTC patients. (F) A total of 2696 DEGs were clustered based on the dissimilarity measure (1-TOM) and were divided into 11 modules. (G) A correlation heatmap between module eigengenes and clinical parameters (AHNAK2 was used as the main research object) of PTC. (H) Scatter plot of blue module eigengenes.

Analysis of the potential mechanism of AHNAK2

Cell adhesion- and cell junction-related pathways were the most frequently noted pathways in the functional enrichment analysis of the 407 blue module genes. “Cell-cell adhesion via plasma membrane adhesion molecules”, “cell-cell junction”, and “cell adhesion mediator activity” were the most significant Gene Ontology (GO) terms for cellular components, biological processes and molecular functions, respectively (Figure 3A). “Cytokine−cytokine receptor interaction” was the most significant pathway in the Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis (Figure 3B). In addition, we constructed a protein–protein interaction (PPI) network with the 57 hub genes and found that ICAM1 and FN1 play key roles (Figure 3C).

Figure 3. Potential mechanisms of AHNAK2 in PTC. (A) GO analysis of 407 blue module eigengenes; the blue, red and green areas represent cellular components, biological processes and molecular functions, respectively. (B) KEGG analysis of blue module eigengenes: the four most significantly enriched pathways. (C) PPIs of hub genes. (D) Enriched pathways in the H-AHNAK2 group based on GSEA. (E) Enriched pathways in the L-AHNAK2 group based on GSEA.

Additionally, gene set enrichment analysis (GSEA) was used to identify the mechanism and functional differences between the L-AHNAK2 group and H-AHNAK2 group. The H-AHNAK2 group was enriched in the following terms: 1. Cell adhesion- and cell junction-related pathways: cell adhesion molecules, focal adhesion, adherens junction, and tight junction; 2. Immune-related pathways: cytokine-cytokine receptor interaction, JAK-STAT signaling pathway, natural killer cell-mediated cytotoxicity, chemokine signaling pathway, and leukocyte transendothelial migration; and 3. Cancer-related pathways: p53 signaling pathway, small cell lung cancer, and pathways in cancer (Figure 3D). Conversely, the L-AHNAK2 group was enriched in metabolism-related pathways: glycine serine and threonine metabolism, fatty acid metabolism, beta alanine metabolism, peroxisome, valine leucine and isoleucine degradation, oxidative phosphorylation, and tryptophan metabolism (Figure 3E).

Further analysis of the relationship between AHNAK2 and immune infiltration

The results of the ESTIMATE analysis suggested that the H-AHNAK2 group had a higher immune score and stromal score than the L-AHNAK2 group, while the tumor purity score was lower (Figure 4A–4D). These results indicate that the proportions of immune cells and stromal cells are higher, and that the proportion of tumor cells is lower in the H-AHNAK2 group than in the L-AHNAK2 group. The TIMER analysis showed that AHNAK2 has significant positive correlations with B cells, CD4⁺ T cells, macrophages, neutrophils and dendritic cells (DCs) in PTC (Figure 4E). The results of the TISIDB analysis suggested that AHNAK2 has positive correlations with 28 tumor-infiltrating lymphocyte (TIL) types (Figure 4F) and human leukocyte antigens (HLAs) (Figure 4G) across human cancers, which is particularly significant in THCA (papillary thyroid carcinoma in TCGA).

Figure 4. AHNAK2 is closely related to immunity in PTC. Comparison of (A) stromal scores, (B) immune scores, (C) ESTIMATE scores, and (D) tumor purity between the L-AHNAK2 and H-AHNAK2 groups based on ESTIMATE. (E) TIMER analysis of purity-corrected partial Spearman’s correlation between the expression of AHNAK2 and six immune cells in THCA. (F) Correlation analysis between the expression of AHNAK2 and 28 types of TILs across human cancers by TISIDB. (G) Correlation analysis between the expression of AHNAK2 and MHCs across human cancers by TISIDB. *P<0.05, ** P<0.01, *** P<0.001, **** P<0.0001 here and in the following figures.

Analysis of AHNAK2 methylation in PTC

Our data show that the methylation level of AHNAK2 in PTC tissues is significantly lower than that in normal tissues (Figure 5A). A significant negative correlation was identified between the methylation level of AHNAK2 and the expression level of AHNAK2 (Spearman r=-0.678, P<0.0001) (Figure 5B). Among the DNA methylation sites of AHNAK2, cg06799735, cg01513078, cg11138227 and cg23385208 had the most significant negative correlations with the AHNAK2 expression level (Table 3).

Figure 5. Analysis of AHNAK2 methylation in PTC. (A) AHNAK2 methylation levels were compared between PTC and normal tissues. (B) Correlation analysis between AHNAK2 methylation and AHNAK2 expression in PTC. (C) Correlation analysis between AHNAK2 methylation and tumor purity in PTC based on LinkedOmics.AHNAK2 methylation levels among PTC patients with different (D) stages, (E) T classifications, (F) N classifications, and (G) histological types. (H) AHNAK2 methylation was highly correlated with genes identified by the Pearson test in the THCA cohort based on LinkedOmics. (I) Heatmaps of the top 50 genes positively and negatively correlated with AHNAK2 in TCHA based on LinkedOmics. (J) Significantly enriched GO biological process analysis and (K) KEGG pathways of AHNAK2 methylation in the THCA cohort.

Moreover, we found significant connections between AHNAK2 methylation and clinical parameters such as tumor purity, stage, T classification, N classification and histological type (Figure 5C–5G).

To understand the biological significance of AHNAK2 methylation in THCA, a functional module of LinkedOmics was used to examine AHNAK2 coexpression pattern in the THCA cohort. Based on RNAseq, we screened 9,944 genes related to AHNAK2 methylation (false discovery rate (FDR) <0.01) (Figure 5H). The top 50 significant genes that were positively and negatively correlated with AHNAK2 methylation are presented as two heatmaps in Figure 5I. The GO (biological process) analysis results derived by GSEA were significant; the results indicated that AHNAK2 methylation coexpressed genes participate primarily in metabolism-related pathways, while genes related to activities such as the adaptive immune response, granulocyte activation, leukocyte migration, regulation of leukocyte activation, and positive regulation of the defense response were inhibited. KEGG pathway analysis also showed that genes related to the complement and coagulation cascades, cytokine-cytokine receptor interaction, cell adhesion molecules, and Epstein–Barr virus infection, among other pathways, were inhibited. These results illustrate the extensive effects of AHNAK2 methylation on transcriptome and immune-related pathways.

Gene expression omnibus (GEO) verification of the diagnostic value of AHNAK2

By analyzing mRNA expression in 8 GEO validation cohorts (GSE3467, GSE3678, GSE5364, GSE27155, GSE33630, GSE58545, GSE53157, and GSE60542), we verified the efficacy of AHNAK2 in distinguishing PTC from normal tissues (Figure 6A–6P).

Figure 6. GEO verification of the diagnostic value of AHNAK2 in PTC. (A–H) mRNA expression levels of AHNAK2 in 8 GEO verification cohorts (GSE3467, GSE3678, GSE5364, GSE27155, GSE33630, GSE58545, GSE53157, and GSE60542). (I–P) ROC curves were generated to verify the diagnostic performance of AHNAK2 expression in the verification cohorts.

Generalization value of AHNAK2 in pan-cancer

To investigate whether AHNAK2 has broad value, we performed a series of studies on AHNAK2 across all cancers. Gene Expression Profiling Interactive Analysis (GEPIA) showed that AHNAK2 expression status varies in different cancers (Figure 7A). TISIDB showed that the high expression of AHNAK2 in pan-cancer tended to be accompanied by an advanced tumor stage (Figure 7B) and short overall survival (OS) (Figure 7C). K-M survival analysis showed that the high AHNAK2 group had significant associations with short OS in bladder urothelial carcinoma (BLCA), glioblastoma multiforme (GBM), lung adenocarcinoma (LUAD), mesothelioma (MESO), pancreatic adenocarcinoma (PAAD), skin cutaneous melanoma (SKCM) and uveal melanoma (UVM) (Figure 7D–7J). AHNAK2 is an extremely significant gene in the survival statistics of a large number of pan-cancer samples (N = 9497, HR = 1.4, P<0.0001) (Figure 7K).

Figure 7. Generalization value of AHNAK2 in pan-cancer. (A) Comparison of AHNAK2 mRNA expression between pan-cancer cancerous and paracancerous tissues based on GEPIA. (B) Associations between AHNAK2 expression and stage across human cancers based on TISIDB. (C) Associations between AHNAK2 expression and overall survival across human cancers based on TISIDB. K-M survival analysis of (D) BLCA, (E) GBM, (F) LUAD, (G) MESO, (H) PAAD, (I) SKCM, (J) UVM, and (K) pan-cancer in the low AHNAK2 and high AHNAK2 groups based on TISIDB and GEPIA. NS: no significance.

Materials

The TCGA THCA dataset, which contains 58 normal thyroid samples (N) and 512 PTC samples (T), was selected as the discovery cohort. The normalized level-three RNA-seq data (FPKM), simple nucleotide variation data (VarScan) and DNA methylation data were downloaded from the TCGA GDC (https://portal.gdc.cancer.gov/). The clinical information on TCGA THCA was downloaded from the University of California at Santa Cruz (UCSC) Xena browser (https://xena.ucsc.edu/).

Eight gene expression microarray datasets, each containing both PTC and normal thyroid tissue samples, were downloaded from the National Center for Biotechnology Information GEO database (http://www.ncbi.nlm.nih.gov/geo): GSE3467 (N=9, T=9), GSE3678 (N=7, T=7), GSE5364 (N=16, T=35), GSE27155 (N=4, T=51), GSE33630 (N=45, T=49), GSE58545 (N=18, T=27), GSE53157 (N=3, T=15), and GSE60542 (N=30, T=33).

Comparison of AHNAK2 expression between THCA and normal thyroid tissues

The mRNA expression levels of AHNAK2 in PTC and normal thyroid tissues were analyzed using TCGA data. We defined the top quarter as the H-AHNAK2 group (based on mRNA expression rank in the THCA dataset) and the remainder as the L-AHNAK2 group.

ROC curve analysis was used to evaluate the diagnostic efficacy of AHNAK2. The AUC followed the criteria 0.50–0.60 = fail, 0.60–0.70 = poor, 0.70–0.80 = fair, 0.80–0.90 = good and 0.90–1 = excellent.

Additionally, a direct comparison of AHNAK2 protein expression between PTC and normal thyroid tissues was performed by analyzing immunohistochemistry images from The Human Protein Atlas (https://www.proteinatlas.org) [28].

WGCNA reveals key modules and hub genes related to AHNAK2

First, DEGs between PTC and normal tissues in the TCGA dataset were identified with the “limma” package in R (adjusted p-value <0.05 and |log₂FC| >1).

WGCNA aims to identify coexpressed gene modules and to explore associations between gene networks and phenotypes of interest, as well as core genes in the network [29]. We divided the DEGs into different gene modules using the “WGCNA” package in R (minModuleSize = 30).

We defined the module with the highest absolute module significance as the key module. We defined the hub genes as those that satisfied the following criteria:

1. In the key module;

2. An absolute value of MM >0.5, which represents Pearson’s correlation between the gene and the module; and

3. An absolute value of GS >0.5, which represents Pearson’s correlation between the gene and the clinical parameter.

Analysis of the potential mechanism of AHNAK2 in PTC

GO and KEGG enrichment analyses

The genes in the key module were subjected to GO and KEGG enrichment analyses with DAVID 6.7 (https://david-d.ncifcrf.gov/).

PPI network analysis

The hub genes were input into the STRING database (version 11.0, https://string-db.org/) to analyze the interaction network of the hub gene-encoded proteins. The minimum required interaction score was set to 0.4 (medium confidence), and the results were visualized using Cytoscape (v3.7.2).

GSEA

GSEA software (Version 2.0.1) (http://www.broad.mit.edu/gsea) was used to explore the potential biological function of AHNAK2 in PTC. A normalized enrichment score (NES) ≥ 1.0, a p-value ≤ 0.05, and an FDR q-val ≤ 0.25 were used as criteria to sort the GSEA results (c2.cp.kegg.v6.2.symbols.gmt was used as the reference gene set) into the L-AHNAK2 and H-AHNAK2 groups after 1000 permutations for each analysis.

Further analysis of the relationship between AHNAK2 and immune infiltration in PTC

ESTIMATE

ESTIMATE [30], a method that uses gene expression signatures to infer the fractions of stromal and immune cells in tumor samples, was used to evaluate the levels of immune cell infiltration (immune score), the stromal content (stromal score), the stromal-immune comprehensive score (ESTIMATE score) and tumor purity for each THCA sample.

TIMER

The TIMER online database [31], a web server for comprehensive analysis of tumor-infiltrating immune cells, was used to analyze and visualize associations between AHNAK2 expression and the abundance of 6 tumor-infiltrating immune cell subtypes (B cells, CD4⁺ T cells, CD8⁺ T cells, macrophages, neutrophils, and DCs).

TISIDB

We used TISIDB, a comprehensive repository portal for tumor-immune system interactions, to determine the Spearman correlations between AHNAK2 expression and 28 TIL types across human cancers [32].

Analysis of AHNAK2 methylation in PTC

Based on DNA methylation data downloaded from TCGA, we compared the methylation levels of AHNAK2 between PTC and normal tissues. We analyzed the correlation between AHNAK2 methylation and AHNAK2 expression in PTC. We also performed a correlation analysis between each methylation site of AHNAK2 and AHNAK2 expression.

LinkedOmics is an online platform for analyzing multiomics data across 32 cancer types (http://www.linkedomics.org/) [33]. We downloaded and visualized the correlation results between AHNAK2 methylation and clinical parameters from LinkedOmics. Pearson correlations were used to analyze the AHNAK2 methylation-induced coexpression of genes, and the results are displayed with a volcano map and a heatmap. The functional modules of LinkedOmics are based on AHNAK2 methylation data, and GSEA was performed with data from GO (biological process) and KEGG pathway databases. We set the minimum number of genes (Size) =10 and simulations=1000 and selected a weighted set cover to reduce redundancy.

GEO verification of the diagnostic value of AHNAK2

mRNA expression and ROC curves were used to verify the diagnostic efficacy of AHNAK2 in 8 GEO datasets (GSE3467, GSE3678, GSE5364, GSE27155, GSE33630, GSE58545, GSE53157, and GSE60542).

Generalization of conclusions in pan-cancer

GEPIA (http://gepia.cancer-pku.cn/) is used as a tool to generalize conclusions in pan-cancer [34]. A comparison of AHNAK2 gene transcripts between PTC and normal thyroid tissues in pan-cancer was performed using TCGA and Genotype-Tissue Expression (GTEx) data (p-value <0.05 and |log₂FC| >1). K-M curves in pan-cancer were generated using the cohort of TCGA data. We used TISIDB to analyze relationships between AHNAK2 expression and tumor stage in pan-cancer, as well as associations between AHNAK2 expression and OS.

Statistical analysis

K-M survival analysis was performed with the "survival" package (PFS was used as the endpoint). The chi-square test was used to test differences in clinical parameters between the L-AHNAK2 and H-AHNAK2 groups. The Spearman method was used to test for correlations. The Mann-Whitney test was used to compare two groups. The log-rank method was used to calculate significant P values for survival. R (v3.6.0) and SPSS version 25.0 software programs were used for statistical processing. Data were visualized with GraphPad Prism V.8.0 software and R.