Characteristics of hypoxic tumor microenvironment in non-small cell lung cancer, involving molecular patterns and prognostic signature.

BACKGROUND: The mechanisms of hypoxia or immune microenvironment in cancer have been studied respectively, but the role of hypoxia immune microenvironment in non-small cell lung cancer (NSCLC) still needs further exploration.
METHODS: By applying the K-means algorithm, 1,121 patients with NSCLC were divided into three categories. We evaluated the constructed signature in order to link it with the prognosis, which was constructed by univariate and least absolute shrinkage operator (LASSO) Cox regression analysis.
RESULTS: A total of three clusters were obtained by clustering five Gene Expression Omnibus (GEO) data sets. Gene Set Variation Analysis (GSVA) and immune infiltration analysis were performed to explore the biological behavior. Cluster one presented an activated state of oncogenic pathways, and compared with the other two clusters, the median risk score was the highest, which was the reason for its poor survival. Cluster three showed that the immune pathway was active and the median risk score was the lowest, so the survival was the best. However, cluster two presented a state in which both immune and matrix pathways were activate. This was manifested as mutual antagonism, and its risk score was in the middle. Its survival was in the middle.
CONCLUSIONS: This work revealed the role of hypoxia related genes (HRGs) modification in tumor microenvironment, which was conducive to our comprehensive analysis of the prognosis of NSCLC, and provided direction and guidance for clinical immunotherapy. 2021 Translational Lung Cancer Research. All rights reserved.

Introduction

Lung cancer is one of the most malignant tumors that threatens the health, and its incidence has increased mainly owing to the rise in smokers, a known risk factor (1,2). Hypoxia related genes (HRGs) play an important role in maintaining energy metabolism, angiogenesis and metastasis of tumor cells, whose essence is mainly due to the vigorous metabolic capacity of tumor cells, resulting in an increase in oxygen content requirements beyond the normal range. Hypoxic tumor microenvironment (TME) contributes to the promotion of tumor cell invasion, thus contributing to the study of cancer mechanisms, progression, and recurrence (3). Some researches have suggested that certain HRGs and their mediators, hypoxic-inducible factors, may serve as prognostic and therapeutic targets for certain cancers, such as colorectal cancer, breast cancer, and hepatocellular carcinoma (4). Genomic analysis has been the main method used internationally to identify new biological targets for lung cancer. Although this method does not reveal a new mechanism, some studies have revealed the significance of tumor-related structures and the signaling pathways of HRGs in the tumor microenvironment. According to the report, there are differences in cell composition within TME, which are specifically reflected in cytotoxic T cells, dendritic cells, helper T cells, mesenchymal stem cells, tumor-related macrophages and related inflammatory pathways (5,6). Therefore, comprehensive identification of the characteristics of TME cell infiltration mediated by multiple HRGs will be good for our realization of the immune regulation of TME. In this research, we integrated genomic information from 1,121 non-small cell lung cancer (NSCLC) samples to comprehensively assess the modification patterns of HRGs and associate these modifications with TME immune infiltration. Three different modification patterns of HRGs were revealed to indicate that HRGs modification patterns play an indispensable role in shaping the individual tumor microenvironment. To this end, we specifically constructed a scoring system to build a scoring system to quantify the modification patterns of HRGs. We presented the following article in according with the MDAR checklist (available at http://dx.doi.org/10.21037/tlcr-20-1314).

Methods

Data sources and data processing for NSCLC

Complete clinical note and gene expression data were downloaded from the Gene Expression Omnibus (GEO) and The Cancer Genome Atlas (TCGA) database. Patients with missing clinical information were deleted. There were six data cohorts involved in this paper, including GSE30219, GSE42127, GSE41271, GSE50081, GSE37745, and TCGA cohort. For the GEO datasets, the normalized matrix files were downloaded, while the RNA sequencing data (FPKM value) came from the TCGA dataset. The “ComBat” algorithm in the sva package was used to correct for batch effects.

Unsupervised clustering for seventy-six HRGs

Seventy-six HRGs were originated from GSEA database (https://www.gsea-msigdb.org/gsea/index.jsp). Unsupervised cluster analysis was applicable for identifying the modification patterns of different HRGs according to the expression of HRGs, and the patients were classified for next research. Consensus algorithm was applied to determine the number of cluster and stability. The steps described above were conducted by the ConsensuclusterPlus package, and they were repeated 1,000 times to ensure the stability of the cluster.

Gene Set Variation Analysis (GSVA)

To further explore the differences of HRGs in biological processes, GSVA enrichment analysis was carried out using GSVA package (7). In a nonparametric and unsupervised approach, GSVA was commonly used to estimate changes in path and bioprocess activity in a sample of an expression dataset. The gene set of “c2.cp.kegg.v6.2.-symbols” originated from MSigDB database, aiming at performing GSVA analysis. P value <0.05 after correction [false discovery rate (FDR)] is generally considered to be statistically significant, which was also the basis for gene selection for GSVA analysis.

Evaluation of immune cell infiltration

The ssGSEA (single-sample gene set enrichment analysis) was used to quantify the relative abundance of 28 immune cell infiltrates. The set of genes used to label each TME infiltrating immune cell type that was derived from the Charoentong research, which contained a variety of human immune cell subtypes, consisting of macrophages, activated dendritic cells, activated CD8 T cells, regulatory T cells, naturally killing T cells, etc. The enrichment scores calculated by ssGSEA analysis was used to reflect the relative abundance of each TME infiltrated cell in each sample (8).

Analysis of the differentially expressed genes (DEG) among HRGs clusters

To verify the reliability and stability of the clustering, we screened out the DEG among the clusters, and the limma package was used to determine DEG. An empirical Bayesian approach was used to estimate the change of gene expression, and the involved method was the moderated t-test. The criterion for screening differential genes was P<0.001.

Construction and validation of HRGs signature

In an integrated GEO training cohort, a univariate Cox regression analysis was performed in the first step to assess the relationship between HRGs expression levels and patient survival. The least absolute shrinkage and selection operator (LASSO) method was used to screen out prognostic genes with P value <0.05. Subsequently, HRGs-based risk scores were constructed to predict the prognosis of NSCLC. Riskscore = Exp (Gene1) × β1 + Exp (Gene2) × β2 + … + Exp (Genen) × βn. Exp represented the expression level of the gene, and β represented the regression coefficient of each gene calculated by the LASSO Cox regression. Then the patients were then divided into low and high risk groups according to the optimal cut-off value of risk score (9). The Kaplan-Meier survival curve was used to analyze the survival differences between low and high risk groups. The time-dependent receiver operating characteristic (ROC) curve analysis was used to evaluate the predictive accuracy of the HRGs-based prognostic signature. What’s more, multivariate regression analysis was applied to determine whether the signature was an independent prognostic factor of NSCLC.

Correlation between HRGs signature and other related biological processes

A set of gene sets were constructed by Mariathasan which stored genes related to some biological process, such as (I) antigen processing machinery; (II) CD8 T-effector signature; (III) immune-checkpoint; (IV) pan-fibroblast TGFb response signature (Pan-F-TBRS); (V) angiogenesis signature; (VI) epithelial-mesenchymal transition (EMT) markers consisting of EMT1, EMT2 and EMT3; (VII) Nucleotide excision repair; (VIII) mismatch repair; (IX) DNA damage repair; (X) WNT target; (XI) antigen processing and presentation; (XII) DNA replication (10).

Statistical analysis

One-way ANOVA and Kruskal-Wallis test, non-parametric and parametric methods, which devoted to the comparison of two or more groups. Spearman and distance correlation analysis were used to calculate the correlation coefficient. The optimal cutoff value for each dataset was evaluated based on the overall patient survival and risk score value. To identify genes in the analysis of differential genes, the Benjamini-Hochberg method converted the P value to FDR. Kaplan-Meier method generated survival curves for the subgroups in each data set, and log rank (Mantel-Cox) test was applied to determine the statistical significance of the differences. All heat maps are generated by the pheatmap function. P values were two-side, and P value of less than 0.05 were regarded statistically significant. The study was in accordance with the Helsinki Declaration (as revised in 2013).

Results

Cluster analysis of five GEO data sets based on HRGs

To further explore the new molecular classification of HRGs expression patterns, 1,121 NSCLC patients which came from a combination of five data sets from the GEO database were subjected to unsupervised consensus clustering (GSE30219, GSE37745, GSE41271, GSE42127 and GSE50081). According to the relative change of the area under the cumulative distribution function (CDF) curve and the uniform heat map, the optimal number of cluster was determined to be three (k value =3), which suggested three molecular patterns, and there was no significant change in the area under the CDF curve (, Figure S1A,S1B). There were 353 cases in the first cluster (31.5%), 409 cases in the second cluster (36.5%), and 358 cases in the third cluster (32%). Kaplan-Meier survival analysis showed that patients in the third cluster had the best overall survival rate (), while patients in the first cluster had the worst overall survival rate. Principal component analysis displayed that there were significant differences among the three different clusters (). What’s more, the heat map showed that the expression of differential genes between different clusters was significantly different, which proved the reliability and stability of the three clusters () (11).

Figure 1

Cluster analysis of HRGs. (A) Unsupervised cluster analysis of 76 HRGs from cBioPortal database, consensus matrices for k =3. (B) The survival analysis of the three HRG modification patterns involved 1,121 patients, of which cluster one included 353 cases (31.5%), cluster two included 409 cases (36.5%), and cluster three included 358 cases (32%). (C) The principal component analysis of the transcriptome profile of the three HRG modification patterns aimed to show the reliability of clustering. (D) In the heat map, cluster, status, stage, gender and age were regarded as patient annotations. Yellow represented high expression, blue represented low expression. HRG, hypoxia related gene; CDF, cumulative distribution function.

Analysis of immune infiltration among different clusters and GSVA enrichment

To explore the biological behavior among these different HRGs molecular patterns, immune infiltration and GSVA enrichment analysis were performed. The box plot displayed that the three molecular patterns were significantly different in twenty-eight kinds of immune cells (), while the heat map similarly showed the expression differences among molecular patterns () (12). Compare cluster one and cluster two, cluster one was obviously enriched in carcinogenic activation pathways such as P53 pathway, PI3K AKT MTOR signaling and WNT beta catenin signaling, while cluster two showed multiple enrichment in the immune pathway such as IL2 STAT5 signaling, IL6 JAK STAT3 signaling and inflammatory response. Compare cluster one and cluster three, cluster one was significantly enriched in hypoxia and carcinogenic activation pathways. For instance, hypoxia, MTORC1 signaling and epithelial mesenchymal transition () (13,14). Kruskal-Wallis test revealed that there were significant differences in HRG score between HRG clusters. In the stromal score, there were significant differences among the clusters. Meanwhile, cluster two showed higher stromal score than the other two clusters, indicating that cluster two was correlated with the matrix activation-related characteristics, which also helped explain that although cluster two was more active in immune infiltration, it did not survive as well as cluster three. Among the immune score, there was no significant difference between cluster one and cluster three. Although cluster two showed the highest immune score, the performance of immune-related characteristics is not obvious due to the characteristics of matrix activation. Furthermore, in the analysis of tumor purity, the differences between clusters were obvious, and cluster three showed higher tumor purity, which further verified that cluster three had the best survival () (15).

Figure 2

Immune infiltration analysis and GSVA Enrichment Analysis. (A) The abundance of each TME infiltrating cell among the clusters, involving twenty-eight kinds of immune cells. The box line of the box plot represented the median value, the black dots represented the outliers, and the asterisk represented P value (***, P<0.001). (B) Expression differences of three different modification patterns among twenty-eight immune cells. Yellow represented high expression, blue represented low expression. (C) Activation and inhibition of different hypoxia modification patterns in biological pathways. Yellow represented high expression, blue represented low expression. (D) Differences of three hypoxia modification patterns in immune microenvironment. Kruskal-Wallis test was used to compare the statistical differences among the three clusters. GSVA, Gene Set Variation Analysis; TME, tumor microenvironment.

Analysis of cluster differences and GO enrichment

By evaluating the stromal microenvironment among the clusters, matrix expression was highest in cluster one and worst in cluster three, which also indicated that cluster one had the worst survival (). Although cluster two played a certain role in the immune pathway, it reduced the influence of the immune pathway due to the activation of the matrix pathway, so the survival status was worse than cluster three. Moreover, there were significant differences among costimulus molecules among different clusters (). By comparing the three clusters in pairs, we used the genes with a P value of less than 0.001 as the difference genes between the two, and screened out 1,694 common difference genes through the Venn diagram (). Through enrichment analysis of differential genes, it was found that these differential genes were mainly enriched in hypoxia pathway, extracellular matrix pathway and immune pathway, such as response to oxygen levels, response to hypoxia, regulation of adaptive immune response, activation of innate immune response, negative regulation of immune system process, extracellular matrix, extracellular matrix organization and so on ().

Figure 3

Selection of differential genes and analysis of GO enrichments. (A) The differences in matrix activation pathways of the three hypoxia modification patterns. EMT, epithelial-mesenchymal transition. (B) The differences in costimulatory molecules of the three hypoxia modification patterns. (C) The condition for screening 1,964 common differential genes was P<0.001. (D) The GO enrichment was designed to functionally annotate hypoxia-related patterns. (*, P<0.05; ***, P<0.001). GO, gene ontology; EMT, epithelial-mesenchymal transition.

Construction of HRG prognostic signature

To further determine the value of HRGs in predicting the prognosis of patients with NSCLC, we constructed a prognosis signature of HRGs by a two-step approach. The first step was univariate Cox regression, and the second step was LASSO regression analysis (16). A total of 18 HRGs were selected to construct the signature. The HRG-based prognostic risk score model was established with the following formula: ExpADM × 0.053 + ExpBIK × 0.009 + ExpDDIT3 × (−0.003) + ExpENO1 × 0.17 + ExpEPAS1 × (−0.09) + ExpFGF3 × 0.051 + ExpGAPDH × 0.018 + ExpMIF × 0.044 + ExpNFKB1 × (−0.192) + ExpPFKP × 0.113 + ExpPGK1 × (−0.095) + ExpPLAUR × 0.0007 + ExpSPP1 × 0.078 + ExpSTC1 × 0.115 + ExpTEK × (−0.045) + ExpTFRC × (−0.02) + ExpTGFA × 0.035 + ExpXRCC6 × 0.0009 (). In the TCGA data set, there are significant differences in the expression of 18 genes in normal tissues and tumor tissues (Figure S1C). We have shown the correlation between genes and clinical factors, and we only list meaningful genes. Patients were divided into high and low risk groups by the optimal cutoff value (17). The survival analysis of signature genes showed that the HRGs of model construction had prognostic significance (). The area under the red curve in ROC curve was 0.69 (). Multivariate Cox regression signature analysis included patient age, sex, stage, and risk score, confirming that risk score was an independent prognostic factor for NSCLC patients (). To further validate the stability of the signature, we applied the HRGs signature constructed on the basis of the merged GEO dataset to the separate GEO database in order to verify its prognostic value (GSE30219, GSE37745, GSE41271, GSE42127, GSE50081) (). Prognostic analysis of HRGs showed that the genes in the signature could significantly predict the significance of patient survival (). In addition, we found significant differences in risk scores among the clusters through the box plot. Cluster one indicated the highest median score, while cluster three indicated the lowest median score (). The Sankey diagram fully demonstrated the association between prognostic signature, clusters, and clinical characteristics (). The network diagram of HRGs in signature was intended to visually show the combined appearance of the correlation between signature genes and their impact on the overall survival of NSCLC patients (). To better demonstrate the features of the HRGs signature, we also explored the link between the signature and the immune microenvironment (). In the three HRGs clusters, significant differences among signature genes were observed, which further verified the significance of clustering (). The signature presented significant differences in the matrix microenvironment (). The vast majority of costimulatory molecules showed significant differences between high and low risk groups (). Most immune cells are significantly different between the high and low risk groups (). Furthermore, the differences of somatic mutation distribution of lung adenocarcinoma (LUAD) and lung squamous carcinoma (LUSC) in TCGA were analyzed to further studying the characteristics of signature genes. The gene NFKB1 in the model had the highest mutation rate in LUSC, reaching 2%, while EPAS1 had the highest mutation rate in LUAD, reaching 2% (, Figure S1D).

Table 1

Gene coef. The coefficients of 18 genes constitute the signature

Gene	Coef.
ADM	0.0528654052451945
BIK	0.00851596297956771
DDIT3	−0.00328450940798743
ENO1	0.170228628163932
EPAS1	−0.0899785170147726
FGF3	0.0511361133771376
GAPDH	0.0176369963829078
MIF	0.0439351480327378
NFKB1	−0.191944037561544
PFKP	0.11322226986168
PGK1	−0.0952325026833525
PLAUR	0.000655402266938651
SPP1	0.0776793659864285
STC1	0.114889187197857
TEK	−0.04500445358283
TFRC	−0.0203908440853215
TGFA	0.0347500036062051
XRCC6	0.0009599431305934

Figure 4

Construction of prognostic signature. (A) The Kaplan-Meier curve aimed to analyze the survival significance of high and low risk groups. (B) The receiver operating characteristic curve aimed to reflect the predictive efficiency of signature. (C) Forest plot aimed to confirm that the signature was an independent prognostic factor for NSCLC. (D,E,F,G,H) Kaplan-Meier curve was used to analyze the meaning of signature in each data set (GSE30219, GSE37745, GSE41271, GSE42127, GSE50081).

Figure 5

Interrelationships among signature factors. (A) Forest plot was used to display the prognostic significance of factors in the signature. (B) Analysis of the differences of signature factors in three hypoxia modification patterns. (C) Alluvial diagram showed the connection and changes of status, stage, cluster and riskScore. (D) Interaction among signature factors. The size of circle represented the influence of each signature on the prognosis. The green dots in the circle represented prognostic risk factors, while the purple dots in the circle represented prognostic favorable factors. The lines connecting the signature factors showed their interaction, and the thickness indicated the mutual strength. Negative correlation was blue, while positive correlation was red. The signature factors were divided into three groups, marked with red, yellow and blue. (***, P<0.001).

Figure 6

Significance of signature factors. (A) The relationship among riskScore and signature factors were analyzed by spearman. (B) The Kruskal-Wallis test was used to compare the differences among three different hypoxia modification patterns and riskScore. (C) Analysis of the signature in the matrix microenvironment, and the asterisk represented P value (*, P<0.05; **, P<0.01; ***, P<0.001). (D) Difference analysis of costimulatory molecules in high and low risk groups. (E) Analysis of the difference between high and low risk groups in immune infiltration. (F,G) Waterfall plot of tumor somatic mutations involving signature factors in lung adenocarcinoma (LUAD) and lung squamous cell carcinoma (LUSC). (F) LUAD (G) LUSC.

Verification of the prognostic signature in the TCGA cohort

To further explore the reliability and stability of the signature construction, we applied the risk score to the TCGA cohort, so that the feasibility of the hypoxia signature was confirmed (). It was noting that multivariate regression analysis suggested that the prognostic signature was indeed an independent prognostic factor for NSCLS (). The risk heat map showed that as the risk score increased, survival was worse and the expression of signature factors were different in the high and low risk groups (). We build a nomogram that could predict 1-, 3- and 5-year OS by signature and other clinical factors. The 1-, 3- and 5-year OS probability calibration curve showed that the OS predicted by nomogram was in good agreement with the actual OS of NSCLC patients. The ROC curve in the nomogram showed that the 1-, 3- and 5-year forecast values were 0.646, 0.652, and 0.615 (). The correlation between the 18 hypoxia-related genes in the constructed model and clinical factors is shown in box plots, and we only show meaningful genes (Figure S2).

Figure 7

Validation of the signature in the TCGA cohort. (A) The Kaplan-Meier curve was designed to verify the prognostic significance of high- and low-risk groups in the TCGA cohort. (B) The forest plot was used to verify that the signature was indeed an independent prognostic factor in the TCGA cohort. (C) The risk plot was used to show that the survival status became worse as the risk value increased. (D) The nomogram contained age, gender, stage, T, N, signature. The x-axis of the calibration chart was the predicted recurrence probability result, and the y-axis was the actual recurrence probability. ROC analysis detected the accuracy of prediction and inspection. TCGA, The Cancer Genome Atlas; ROC, receiver operating characteristic.

Efficacy of the signature in immunotherapy

The data on the smoking history of patients came from the cBioPortal database. It divided patients into four categories, including current reformed smoker for < or =15 years, current reformed smoker for >15 years, current smoker and lifelong non- smoker. Through analysis, it was shown that the six molecules in the signature had significant differences in these four categories (). The risk score and GAPDH displayed significance in follow up treatment (). And it was worth noting that EPAS1, TFRC, and XRCC6 had significant meanings in the primary therapy outcome (). In order to reflect the significance of prognostic signature in immunotherapy, it was verified in melanoma () and renal carcinoma () respectively, and the results showed that the HRGs-built signature had significant significance in immunotherapy (18,19).

Figure 8

Efficacy of the signature in immunotherapy. (A) Smoking correlation analysis of signature factors. (B) Significance of signature factors follow-up treatment. (C) Significance of signature factors primary therapy outcome. (D) Significance of signature factors in immunotherapy in melanoma. (E) Significance of signature factors in immunotherapy in renal cell carcinoma.

Discussion

Studies have found that HRGs play an indispensable role in inflammation, immunity and tumor growth and development. Moreover, it is worth noting that the overall TME infiltration characteristics mediated by the comprehensive action of multiple HRGs have not been fully explored, which is also the original intention of data mining in this study. Clarifying the role of different HRGs modification patterns in the infiltration of TME cells will help us deepen our understanding of the anti-tumor immune response of TME and guide more valid immunotherapy strategies. On the basis of 76 HRGs, we exposed the modification patterns of three different HRGs. These three patterns not only showed significant differences among the clusters, but also displayed significant TME cell infiltration characteristics. Cluster one was qualified by the activation of stromal and carcinogenic pathways, equivalent to immune-excluded phenotype. Cluster two was qualified by competing activation of immune and stromal pathways that was resemble to antagonistic effects, equivalent to immune-desert phenotype. Cluster three was qualified by the activation of immune-related pathways, equivalent to immune-activation phenotype, which presented as a mass infiltration of immune cells in TME. Though the immune-desert phenotype also showed the infiltration of a variety of immune cells, the activation of the matrix pathway lead to mutual antagonism, which weakened the anti-tumor effect of immune cells, so the survival of the immune-desert phenotype was slightly worse than that of immune-activation phenotype. Compatible with the above definition, we discovered that cluster two showed an important matrix activation state, for example, cluster two showed high expression of angiogenesis, EMT, immune checkpoint and Pan-F-TBRS. Interestingly, in this research, the mRNA transcriptome differences among different HRGs modification patterns have been shown to be markedly related to hypoxia and immune-related biological pathways. The DEG were regarded as HRGs, resemble to the clustering results of hypoxia-related gene modification phenotypes. Based on hypoxia characteristic genes, three sets of genome subtypes have been identified, and they were also significantly related to matrix and immune activation, which again proved that hypoxia modification patterns were of great importance for shaping different TME landscapes. Therefore, we constructed a prognostic model of hypoxia-related genes, and by analyzing the clustering and risk score, we could find that cluster one had a higher risk score and cluster three had a lower risk score. This was also consistent with our predicted results in IMvigor210 cohort. This suggested that the prognostic model could be applied to further determine the pattern of TME infiltration, namely the tumor immunophenotype.

In clinical practice, the prognostic signature can be applied to the HRG modification model and the corresponding TME cell immune infiltration model in a single patient, which is helpful to further determine the immunophenotype of the tumor and to fully realize the clinical practice. In addition, we further assessed the clinicopathologic features of the tumor, including tumor stage and mutation burden. More importantly, this research provides new directions for identifying different immunophenotypes and personalized immunotherapy for cancer. Reversing the characteristics of immune infiltration may help to reverse the progression of cancer, which may provide new ideas for the exploration of cancer treatment and mechanisms. Studies have found that the genes in the model play an important role in cancer. The original purpose of this study was to link these eighteen genes associated with hypoxia to predict NSCLC. Adrenomedullin (ADM) works by supplying and amplifying pathologic neovascularization and lymphatics to generate vital signals. Aryl hydrocarbon receptor (AHR) is part of the tobacco-induced carcinogenic effect, and ADM can greatly promote it (20,21). BCL2 interacting killer (BIK) reduces proliferative epithelial cells by releasing calcium from endoplasmic reticulum storage and causing apoptosis, and BIK peptides may have therapeutic potential in airway diseases related to chronic mucosal hypersecretion (22). DNA damage inducible transcript 3 (DDIT3) promotes stress-induced apoptosis and acts as a transcription factor in many diseases (23). As a glycolytic enzyme, ENO1 plays a key role in glucose metabolism and contributes to many cancer tumors. The reduction of phosphorylation of ENO1 leads to the reduction of phosphorylation of PI3K and AKT, which is further related to the decrease of cell proliferation and tumor progression (24,25). NSCLC may have a mechanism of EPAS1 negative feedback regulation (26). FGF3 is a gene belonging to the fibroblast growth factor gene family. It is related to the carcinogenesis of a variety of cancers, including esophageal squamous cell carcinoma, breast cancer and lung cancer (27). Respiratory inhibitor 3-bromopyruvate (3BP) can inhibit respiration by down-regulating the expression of HK-II and GAPDH, thereby significantly reducing the intracellular oxygen consumption rate and reducing the hypoxia of tumors. The high glycolysis rate in cancer cells is a recognized feature of many human tumors, which can provide metabolites that can be used as precursors of anabolic pathways for cancer cells whose metabolites proliferate rapidly. Maintaining a high glycolysis rate depends on the NADH produced by GAPDH on the regeneration of NAD catalyzed by lactate dehydrogenase, because an increase in the ratio of NADH:NAD will inhibit GAPDH (28,29). Macrophage migration inhibitory factor (MIF) is an inflammatory cytokine that plays many roles in inflammation and immunity. MIF is overexpressed in these malignancies in humans and leads to cell cycle disorders, angiogenesis and metastasis (30,31). An inflammatory cytokine that promotes EMT and invasiveness of cancer cells. SLC26A4-AS1 promotes the transcriptional activity of NPTX1 by recruiting NFKB1, thereby producing anti-angiogenic effects on glioma cells (32). The platelet subtype of phosphofructokinase (PFKP) is a rate-limiting enzyme involved in glycolysis. Overexpression of PFKP inhibits apoptosis and promotes the proliferation, migration, invasion and glycolysis of NSCLC cells (33,34). PFK enhanced the role of glycolysis through PPI3K activation and PI3K/AKT-dependent positive feedback regulation. Phosphoglycerate kinase (PGK1), which is both glycolytic enzyme and protein kinase, plays a crucial role in the development of cancer. Inhibition of glycolysis has been an attractive method to treat cancer since evidence suggests that tumor cells rely more on glycolysis than oxidative phosphorylation. Preliminary evidence suggests that inhibition of phosphoglycerate kinase 1 (PGK1) kinase activity will reverse the Warburg effect and cause tumors (35,36). PLAUR induces human LUAD cells to be resistant to gefitinib through the EGFR/P-AKT/survivin signaling pathway (37). SPP1 overexpression is related to the poor prognosis of multiple cancers, while silencing SPP1 inhibits the proliferation, migration and invasion of cancers (38). STC1 is a glycoprotein known to be involved in inflammation and tumor progression. Cancer cells overexpressing STC1 have an inhibitory effect on the migration/infiltration of macrophages. In addition, STC1 is related to calcium/phosphate homeostasis (39). TEK has been identified as a specific receptor for endothelial cells, which plays an important role in the regulation of angiogenesis and remodeling (40). TFRC promotes the proliferation and metastasis of cancer cells by up-regulating the expression of AXIN2 to accelerate the progress of epithelial ovarian cancer (41). Transforming growth factor (TGFA) stimulates cell growth and proliferation, and its overexpression is associated with the survival rate of patients with many tumors (42). PARP6 inhibits the expression of XRCC6 by inducing degradation, thereby affecting the Wnt/β-catenin signaling pathway, thereby helping to inhibit HCC (43). The DNA protein kinase subunit XRCC6 is necessary for the specific response of these outer membrane proteins in macrophages and epithelial cells (44).

The article’s supplementary files as