Potential role of glucosamine-phosphate N-acetyltransferase 1 in the development of lung adenocarcinoma

The expression of GNPNAT1 in LUAD

We obtained data regarding the expression of GNPNAT1 in LUAD and normal lung tissues from the Cancer Genome Atlas (TCGA) and the Gene Expression Omnibus (GEO) database. There were five cohort meta-analyses of the differential expression of GNPNAT1 in the Lung Cancer Explorer (LCE) database, showing that the mRNA GNPNAT1 expression level was significantly higher in LUAD tumor tissues than in normal tissues (Figure 1A). We obtained similar results for the mRNA GNPNAT1 expression features using the Gene Expression Profiling Interactive Analysis (GEPIA) database (Figure 1B) and Oncomine database (Figure 1C).

We used the UALCAN database to identify the mRNA GNPNAT1 expression in TCGA LUAD samples, grouped by characteristics such as age, sex, stage, race, and smoking status. (Figure 2A). In all subgroups, the mRNA GNPNAT1 expression level was higher in LUAD tumor tissues than in normal tissues. Therefore, GNPNAT1 might be a potential biomarker for the diagnosis of LUAD.

To further assess GNPNAT1 expression levels in LUAD, we also detected GNPNAT1 protein levels in LUAD tissues and normal tissues by immunohistochemical staining in the Human Protein Atlas (HPA) database. As shown in Figure 2B, GNPNAT1 protein expression in cell cytoplasmic/membranous staining intensity was moderate or low in LUAD tissues, but it was not detected in normal tissues.

GNPNAT1 expression associated with survival in LUAD

This study aimed to determine the GNPNAT1 expression associated with survival in patients with LUAD. Kaplan-Meier survival curves were used to identify survival outcomes in the multiple LUAD cohorts. The median value of GNPNAT1 expression was the cutoff value, and each LUAD cohort was divided into a high and a low GNPNAT1 expression group. Compared to the low GNPNAT1 expression group, the high GNPNAT1 expression group had significantly poorer overall survival (OS) (log-rank test, p < 0.05) in the LUAD TCGA cohort (Figure 3A) and in the GSE72094 cohort (Figure 3B). From the LCE database, we got a 17 cohorts meta-analysis of GNPNAT1 expression associated with survival in LUAD; Figure 3C shows that the high GNPNAT1 expression group had poorer OS (test for all cohorts, p < 0.01, HR = 1.27) compared to the low GNPNAT1 expression group.

The co-expression networks of GNPNAT1 in LUAD

We used LinkedOmics to obtain the GNPNAT1 co-expression networks in the LUAD TCGA cohorts. There were 4039 positively co-expressed genes and 6654 negatively co-expressed genes with GNPNAT1 (FDR < 0.01); in Figure 4A, the dark red dots represent positively correlated genes and the dark green dots represent negatively correlated genes. We drew heatmaps with the 50 most significant positively and negatively co-expressed genes, respectively (Figure 4B). All these significantly correlated genes are shown in Supplementary Table 1.

The five most significant genes positively associated with GNPNAT1 expression were cyclin dependent kinase inhibitor 3 (CDKN3), cyclin B1 (CCNB1), DLG associated protein 5 (DLGAP5), NIMA related kinase 2 (NEK2), and mitotic arrest deficient 2 like 1 (MAD2L1). The five most significant negatively associated genes were chromobox 7(CBX7), Golgi-associated, gamma adaptin ear-containing, and ARF-binding protein 2 (GGA2), calpain 3 (CAPN3), solute carrier family 27 member 1 (SLC27A1), and enhancer of zeste 1 polycomb repressive complex 2 subunit (EZH1). We used GEPIA to identify the 10 genes associated with OS in LUAD. Kaplan-Meier survival curves are shown in Figure 4C. All the five positive high expressions showed significant risk genes in LUAD (p < 0.05); conversely, the five negative high expressions showed protective genes in LUAD (p < 0.05; except CAPN3, p = 0.072).

We used gene set enrichment analysis (GSEA) to annotate GNPNAT1 co-expressed genes with a significant GO term. It showed that GNPNAT1 co-expressed genes are mainly involved in chromosome segregation, ncRNA processing, organelle fission, ribonucleoprotein complex biogenesis, rRNA and tRNA metabolic process, etc., and in inhibited ion transmembrane transport, transporter activity, metal ion transport, neuron projection guidance, export across plasma membrane, and so forth (Figure 4D and Supplementary Table 2). KEGG pathway analysis showed enrichment in the ribosome, proteasome, RNA transport, ribosome biogenesis in eukaryotes, spliceosome, and so forth (Figure 4E and Supplementary Table 3).

GNPNAT1 networks of kinase, miRNA, and transcription factor targets in LUAD

To further explore the targets of GNPNAT1 in LUAD, we analyzed the kinase, miRNA, and transcription factor target networks of GNPNAT1 co-expressed genes. The top five most significant kinase-target networks were cyclin-dependent kinase 1 (CDK1), polo-like kinase 1 (PLK1), Aurora kinase B (AURKB), cyclin-dependent kinase 2 (CDK2), and Aurora kinase A (AURKA) (Table 1 and Supplementary Table 4). Except for CDK2, the expression of these genes in LUAD tissue was higher compared with that in normal tissue. Furthermore, except for CDK2, for all the kinase genes, higher expression was associated with poor OS of LUAD (Supplementary Figure 1). CDK1 is a diagnostic biomarker and a prognostic biomarker in LUAD [10, 11].

We obtained miRNA-target networks of GNPNAT1 co-expressed genes by GSEA, but there was no statistical significance. (Supplementary Table 5). The transcription factor-target networks related primarily to the E2F transcription factor family, including E2F-Q6, E2F1-Q6, E2F-Q4, E2F1DP1_01, and E2F1DP2_01 (Supplementary Table 6). High expression levels of E2F family genes are associated with an unfavorable prognosis in LUAD [12, 13].

Genomic alterations of GNPNAT1 in LUAD

We used the cBioPortal to determine the types and frequency of GNPNAT1 alterations in LUAD based on sequencing data from LUAD patients in the TCGA database. There were 39 out of 522 (7.47%) LUAD patients with GNPNAT1 alterations (Figure 5A). The alterations were mRNA upregulation in 23 cases (4.41%), amplification in 11 cases (2.11%), mutation in 1 case (0.19%), deep deletion in 1 case (0.19%), and multiple alterations in 3 cases (0.57%). Although GNPNAT1 alterations in LUAD were not frequent, mRNA upregulation and amplification were the most common types of GNPNAT1 copy number variation (CNV) in LUAD. We estimated the frequency distribution of GNPNAT1 CNV patients in different stage groups; as shown in Figure 5B, stage I and stage II patients had a high occurrence of GNPNAT1 CNV alteration in LUAD. We divided the patients with LUAD into a GNPNAT1 CNV alteration group and a no GNPNAT1 CNV alteration group; GNPNAT1 CNV alteration group patients had poorer OS than patients in the no GNPNAT1 CNV alteration group (p < 0.05) (Figure 5C). The median survival time was 32.82 months and 49.31 months, respectively.

Gene co-occurrence of GNPNAT1 alterations in LUAD

Gene co-occurrence reflected common genetic risk factors constituting functional relationships. We subsequently identified the co-occurrence genes with GNPNAT1 amplification in LUAD. There were 120 significant co-occurrences with GNPNAT1 amplification genes, as shown in Figure 5D (Supplementary Table 7). The top 10 alterations were solute carrier family 25 member 21 (SLC25A21), DDHD domain containing 1(DDHD1), endoplasmic reticulum oxidoreductase 1 alpha (ERO1A), fermitin family member 2 (FERMT2), G-protein-coupled receptor 137C (GPR137C), proteasome 26S subunit, ATPase 6 (PSMC6), serine/threonine/tyrosine interacting protein (STYX), thioredoxin domain containing 16 (TXNDC16), long intergenic non-protein coding RNA 517 (LINC00517), bromodomain adjacent to zinc finger domain 1A (BAZ1A), shown in Figure 5E. Enriched KEGG pathway analysis indicated that these co-occurrent genes were mainly enriched in Hippo signaling pathway-multiple species (Figure 5F). This signal pathway is primarily associated with the proliferation and apoptosis of tumor cells. GO term analysis of these genes showed enrichment in the G1/S transition of the mitotic cell cycle (Figure 5G and Supplementary Table 8), which also indicated that these genes participate in the growth of tumors.

Moreover, we constructed a GNPNAT1 co-occurrence gene protein-protein interaction (PPI) network using lung-specific data collected from the DifferentialNet database (Figure 5H). Proteasome 20S subunit alpha (PSMA) 3, PSMA4, PSMA5, tripartite motif containing 9 (TRIM9), and GTP cyclohydrolase 1(GCH1) were the top 5 hub genes. PSMA is a proteasome subunit alpha type associated with the occurrence and development of multiple cancers [14, 15].

Finally, we constructed the GNPNAT1 co-occurrence genes TF-miRNA coregulatory interactions using the RegNetwork repository (Figure 5I). Forkhead box A1(FOXA1), Kelch-like family member 28 (KLHL28), cofilin 2 (CFL2), SIX homeobox 4 (SIX4), and glia maturation factor beta (GMFB) were the top five TFs. Many studies have shown that FOXA1 participates in the development of lung cancer, prostate cancer, and several other types of cancers [16–19].

Association between GNPNAT1 expression and immune infiltration level in LUAD

We used the TISIDB database to assess whether GNPNAT1 expression was significantly correlated with immune cell infiltration level in LUAD, as shown in Table 2 (Figure 6). There were 20/28 significantly correlated immune cells, including B cells, CD4 T cells, CD8 T cells, CD56dim cells, eosinophils, IMM B cells, macrophages, mast cells, MDSC cells, neutrophils, NK cells, NKT cells, pDC cells, Tem CD8 cells, Tfh cells, Tgd cells, Th1 cells, Th2 cells, Th17 cells, and Treg cells. We also observed associations between GNPNAT1 expression and immune cell infiltration levels across different cancer types (Supplementary Figure 2A).

Immunomodulators included immunoinhibitors, immunostimulators, and MHC molecules, which regulate immune system functions. We found that GNPNAT1 expression was significantly negatively correlated with immunomodulators (Table 3). The five most highly correlated immunoinhibitors (Figure 7A) were galectin 9 (LGALS9), adenosine A2a receptor (ADORA2A), transforming growth factor beta 1 (TGFB1), colony-stimulating factor 1 receptor (CSF1R), B and T lymphocyte associated (BTLA). The five most highly correlated immunostimulators (Figure 7B) were transmembrane protein 173 (TMEM173), TNF receptor superfamily member 14 (TNFRSF14), CD40 ligand (CD40LG), TNF receptor superfamily member 13B (TNFRSF13B), and PVR cell adhesion molecule (PVR). The five most highly correlated MHC molecules (Figure 7C) were major histocompatibility complex, class II, DP beta 1 (HLA-DPB1), major histocompatibility complex, class II, DO alpha (HLA-DOA), major histocompatibility complex, class II, DM alpha (HLA-DMA), major histocompatibility complex, class II, DR beta 1 (HLA-DRB1), and major histocompatibility complex, class II, DP alpha 1 (HLA-DPA1). We also observed associations between GNPNAT1 expression and immunomodulators across different cancer types (Supplementary Figure 2B–2D).

Chemokines expression levels are key factors in controlling immune cell infiltration. We identified the correlation between GNPNAT1 expression and chemokines (Table 4). The five most highly correlated chemokines (Figure 8A) were C-C motif chemokine ligand (CCL)-14, C-X3-C motif chemokine ligand 1 (CX3CL1), CXC motif chemokine ligand 8 (CXCL8), CCL17, CCL19. The top five chemokine receptors (Figure 8B) were C-X3-C motif chemokine receptor 1 (CX3CR1), CC motif chemokine receptor (CCR)-6, CCR7, CCR4, CXC motif chemokine receptor 5 (CXCR5). Moreover, we described the correlation of GNPNAT1 expression with chemokines or receptors across different cancer types (Supplementary Figure 3).

We evaluated these 25 immune genes correlated with GNPNAT1, including 23 negatively correlated genes and 2 positively correlated genes, and we found that higher expression of 18/23 negatively correlated genes was associated with favorable prognosis in patients with LUAD, while higher expression of 2/2 positively correlated genes was associated with unfavorable prognosis. (Table 5 and Figure 9).

LCE database analysis

LCE is a lung cancer-specific database including expression data and clinical data from over 6700 patients in 56 studies [44]. We could easily obtain an overview analysis of tumor versus non-malignant tissue (normal tissue) differential gene expression and expression–survival association by meta-analyses. In addition, we obtained an individual data set-based survival analysis, comparative analysis, and correlation analysis.

GEPIA database analysis

The GEPIA database (http://gepia.cancer-pku.cn/) is an interactive website that contains 9736 tumor samples and 8587 normal samples from TCGA and GTEx datasets [45]. We used GEPIA to generate OS curves, based on gene expression with the log-rank test and the Mantel-Cox test in LUAD. We obtained tumor versus non-malignant tissue (normal) gene differential expression map. The threshold values were p-value of 0.05 and fold change of 1.0.

Oncomine database analysis

The Oncomine database (https://www.oncomine.org/) is a cancer microarray database and web-based data-mining platform. The gene GNPNAT1 expression level in LUAD was examined in the Oncomine 4.5 database. We used a p-value of 0.05, fold change of 1.2, and gene ranking of all as the threshold values.

UALCAN database analysis

The UALCAN (http://ualcan.path.uab.edu) database is a comprehensive, user-friendly, and interactive web resource for analyzing cancer omics data [46]. We used the UALCAN database to obtain the gene GNPNAT1 expression analysis across LUAD and normal samples in various tumor sub-groups based on TCGA data.

HPA database analysis

The HPA (http://www.proteinatlas.org) database maps human proteins in cells, tissues, and organs using the integration of various omics technologies [47]. GNPNAT1 protein expression in LUAD tissues and normal lung tissues was mapped by immunohistochemistry.

LinkedOmics database analysis

The LinkedOmics database (http://www.linkedomics.org/login.php) is a publicly available portal that includes multi-omics data from all 32 TCGA cancer types [48]. We used LinkedOmics to gain the GNPNAT1 co-expression assessed by Pearson’s correlation coefficient statistical analysis, presented in volcano plots and heat maps. We used the GSEA function modules in the LinkedOmics database to obtain Gene Ontology biological process (GO_BP), KEGG pathways, kinase-target enrichment, miRNA-target enrichment, and transcription factor-target enrichment analysis. FDR < 0.05 and 1000 simulations were the standard in this enrichment analysis.

c-BioPortal database analysis

The cBioPortal (http://cbioportal.org) includes multidimensional cancer genomics [49]. We used the c-BioPortal tool to analyze GNPNAT1 mutation, gene co-occurrence, and CNV in LUAD (TCGA, Firehose Legacy). We also obtained the GNPNAT1 alterations overview per sample via this portal.

NetworkAnalyst database analysis

NetworkAnalyst 3.0 (https://www.networkanalyst.ca/) [50] is a tool that allows to create cell-type or tissue-specific PPI networks, gene regulatory networks, enrichment networks, and gene co-expression networks. All the GNPNAT1 co-occurrence gene networks were built using this tool.

TISIDB database analysis

TISIDB is a web portal for tumor and immune system interaction, which integrates data regarding 30 cancer types from TCGA (http://cis.hku.hk/TISIDB/) [51]. We used TISIDB tools to analyze the correlation of GNPNAT1 with 28 tumor-infiltrating lymphocytes (TILs), 45 immunostimulators, 24 immunoinhibitors, 21 MHC molecules, 41 chemokines, and 18 receptors.

Statistical analysis

We used Student’s t-tests to identify the different GNPNAT1 expression levels. Kaplan-Meier curves and the log-rank test were used to compare the OS of various gene expression subgroups. The correlation between GNPNAT1 expression, immune cell infiltration, and immune genes was evaluated by Spearman’s method. In these analyses, p < 0.05 was considered statistically significant.