The hypoxia-associated genes in immune infiltration and treatment options of lung adenocarcinoma

Data acquisition of LUAD patients

The analysis process of this study is shown in the flow chart (Fig. 1). The UCSC Xena Browser (https://xenabrowser.net/) was used to download mRNA expression and clinical information for LUAD patients in TCGA, with 503 tumor samples and 35 normal samples. The samples of “primary solid tumors” were selected for follow-up studies. The downloaded RNA-seq dataset was in FPKM format, which reflected the actual expression quantity of the gene when comparing multiple samples. RNA-seq data and clinical information, including 695 tumor samples from three LUAD cohorts from the Gene Expression Omnibus (GEO) (GSE72094, GSE31210, GSE30219) were acquired. To ensure the reliability of the study, with an overall survival (OS) time of ≥30 days were retained and unexpressed genes in ≥50% of the samples were removed to improve data quality. A total of 1198 LUAD patients from four datasets were included. Detailed clinical information for LUAD patient clinical cohorts is summarized in Table S1.

Identification of survival-related HRGs in LUAD

HRGs were acquired from the MsigDB database (https://www.gsea-msigdb.org/gsea/msigdb), and the latest HRGs were saved in an ensemble format. Univariate Cox model regression analysis was performed to screen further the GSE72094 cohort for HRGs associated with patient OS information (P value <0.05).

Identification of hypoxia subtypes of the LUAD cohort based on consensus clustering

As a typical unsupervised clustering method, consensus clustering is a common method for classifying cancer subtypes (Lock & Dunson, 2013). The prognosis-related HRG expression screened by univariate Cox analysis was used to subtype the LUAD cohort in TCGA and three GEO databases, and the most reasonable clustering results were determined by comparison.

Immunological characteristics and clinical features of hypoxia subtypes

Based on human genetic set data, samples from TCGA were subjected to GSVA analysis to obtain pathway scores, and formal difference analysis of pathways between hypoxic subtypes was performed using the limma package (Ritchie et al., 2015) to assess differences in subtype immunological characteristics (Hanzelmann, Castelo & Guinney, 2013). Single-sample gene set enrichment analysis (ssGSEA) was based on 29 immune cell types to quantify immune infiltration in the TCGA cohort using the “GSVA” package (Zhao, Ou & Lin, 2021). Differences in clinical characteristics between the hypoxic subtypes were also evaluated.

HRG prognostic signature development and validation

The GSE72094 dataset was used as a training group based on careful consideration of the number of clinical patients and the sequencing platform. It was tested in three other independent datasets (i.e., TCGA-LUAD, GSE31210 and GSE30219) for the predictive effect of these features to find potential prognostic biomarkers. HRGs related to survival in the training group were further selected and validated by LASSO regression using the “glmnet” package. Multivariate Cox model regression analysis was conducted to filter the HRGs associated with ultimate survival to develop a risk signature for LUAD patients. The characteristics with univariate Cox model results in P < 0.05 were selected for inclusion in the multivariate Cox model analysis. The risk score was calculated using the formula: risk score = ∑coefi * HRGi, where coefi represents the coefficient of i HRG in multivariate Cox model regression analysis, and HRGi denotes the expression level.

LUAD patients were divided into high- and low-risk groups based on the median value of the risk score. Kaplan–Meier (KM) survival analysis and time-ROC curves were used to evaluate the risk model performance in predicting prognosis.

Deterministic curve analysis (DCA) is a method to determine the clinical utility of different predictive models (Kerr et al., 2016). Calibration plots were used to verify the discrimination of the nomograms. The “index.comp” package compared the C-index between different groups.

Tumor mutation burden analysis

Based on somatic mutation data from TCGA, TMB information was obtained for each patient sample in the LUAD-TCGA cohort, and the difference in mutated genes between hypoxic subtypes or risk groups was compared by the R package “maftools”.

Tumor-infiltrating immune cells (TIICs) and immune checkpoint inhibitor analysis

We utilized the MCP-COUNTER (Aran, Hu & Butte, 2017) algorithm to compare immune infiltration differences among risk groups in LUAD-TCGA through the “immunedeconv” package. Additionally, we utilized the “estimate” package to calculate tumor purity scores based on the ratios of immune and stromal cells, ultimately obtaining immune scores, stromal scores, and microenvironment scores.

To evaluate the response to immune checkpoint inhibitors (ICI), we employed the Tumor Immune Dysfunction and Exclusion (TIDE) algorithm, which is capable of predicting the effectiveness of immune checkpoint suppression therapy.

Gene set enrichment analysis

GSEA was performed to reveal KEGG pathways in diverse risk groups of TCGA-LUAD (Subramanian et al., 2005). P < 0.05 and —NES—>1 represent statistical significance. Gene-Concept Network was performed to show complex connections between genes and enriched pathways for visual network presentation by “cnetplot” in R.

Evaluation of candidate targeted drugs

In order to evaluate the predictive capability of the small molecule drug susceptibility risk model, the half-maximal inhibitory concentration (IC50) scores among all candidate targeted drugs were estimated between two risk groups using the “oncoPredict” package. The Connectivity Map (CMap) database (https://clue.io/) was applied to access potential molecule drugs that were highly correlated with both risk groups.

Survival analysis of single HRGs

Prognostic analysis was performed separately for each gene constituting the HRG risk model using sample data from the GSE72094 training set. A minimum P value approach was used to assess the optimal survival cutoff, and patients were divided into high- and low-risk groups based on overall survival.

RT–PCR

Reverse-transcriptase PCR was performed on HRGs. The control gene was GAPDH. The total RNA of A549 cells was extracted by the TRIzol method. The ex-tracted RNA was reverse transcribed into cDNA, followed by PCR amplification. Samples were preamplified at 95 °C for 10 min followed by 40 amplification cycles for 2 s at 95 °F, 20 s at 60 °C, and 10 s at 70 °C. Subsequently, the product DNA was electrophoresed on a 2% agarose gel. The primers used for RT–PCR assays are listed in Table S2. The A549 cell line was obtained from Dr. Jing Liu at the National Center for Nanoscience and Technology, China (NCNST).

Cell transfection analysis

The day before transfection, cells were seeded into 6-well plates, which guar-anteed 70%–90% confluence of cells at the time of transfection. RNAiMAX was used with 20 µmol/L siRNA- or negative control-transfected cells. For one well of the 6-well plate, 4 µL of siRNA and 9 µL of RNAiMAX were dissolved in 125 µL of opti MEM medium, mixed at room temperature, left to stand for 10 min, and then added dropwise to the wells. After 48 h, the transfected cells were collected, and some cells were isolated for total RNA according to a previous method. The knock-down efficiency was measured by RT–PCR with Vazyme’s one-step kit, and the cell phenotypes of the remaining cells were measured, such as proliferation, apoptosis, migration, and invasion. RT–PCR was performed using GAPDH as the internal ref-erence gene, and each group of experiments had three repetitions. siRNA sequences are shown in Table S3.

Cell proliferation analysis

Cell proliferation was assessed by a cell counting kit-8 (CCK-8) Kit. The CCK-8 assay is commonly used to measure cell proliferation, which involves the use of the WST-8 reagent. This reagent comprises of 2-(2-methoxy-4-nitrophenyl)-3-(4-nitrophenyl)-5-(2,4-disulfophenyl)-2H-tetrazolium monosodium salt, which can be reduced by dehydrogenases present in the mitochondria of viable cells. Upon reduction, a water-soluble orange-colored formazan dye is formed. In cell experiments, the Optical Density (OD) value, which indicates the degree of light absorption by the sample, is measured using a spectrophotometer. This value is usually represented as the measured absorbance value. The OD value measured at a wavelength of 450nm in the microplate reader can be used as an indirect measure of cell quantity and proliferation. Briefly, 1 × 103 cells were seeded in each 96-well plate, and initial values were recorded after attachment to the wells and further incubated for 1, 2, 3, 4, and 5 days. Then, 100 µL of CCK-8 reagent was added to each well, and absorbance was detected at 450 nm. There were 6 repli-cates of each set of experiments.

Flow cytometric assay for apoptosis

Apoptosis detection was conducted using an Annexin V-PE/7-AAD Apoptosis Detection Kit. The siRNA-transfected A549 cell lines were cultured for 48 h. Cells harvested with 0.25% trypsin were washed twice more with cold PBS. Cells were stained with annexin v-PE and 7-AAD. The apoptosis rate was detected by flow cytometry. Each set of experiments was repeated three times.

Transwell assay

Cell invasion status was analyzed using 24-well Transwell chambers. After transfection, cells were inoculated into the upper chamber at a density of 1 × 10 $\stackrel{ˆ}{}$5 cells/well using 200 µL of blank medium, which was precoated with Matrigel. The lower chamber was filled with 600 µL of RPMI-1640 containing 20% FBS. Cells were cultured at 37 °C for 24 h, fixed in paraformaldehyde for 15 min, and stained with 0.5% crystal violet for 15 min. Cells on the upper side of the mem-brane were removed with a clean cotton swab, and the number of invading cells was observed with a Nikon Eclipse Ti fluorescence microscope at 20×, and five randomly selected areas were counted under the microscope.

“Scissor” algorithm to identify subpopulations of TAMs associated with the LUAD hypoxic subtype phenotype in single-cell sequencing data

The “Scissor” algorithm presents an innovative approach to identify specific subpopulations from single-cell data that exhibit a high correlation with phenotypic characteristics. This approach is based on integrating phenotypic information obtained from large volumes of RNA-seq data, which allows for a more comprehensive understanding of the cellular heterogeneity within a population (Sun et al., 2022). The algorithm quantified the similarity between bulk and single-cell sequencing data and subsequently optimized the regression model with a sample phenotype correlation matrix to determine subtype classification of single cells. Single-cell sequencing data were obtained from 10 normal lung and 10 fresh lung adenocarcinoma tissue samples taken from clinical patients and processed for single-cell RNA sequencing (Bischoff et al., 2021). TAMs in the single-cell data based on marker genes were screened for subsequent subgroup classification.

Construction of survival-related HRGs in LUAD

A total of 214 HRGs were obtained from the MsigDB database. A total of 197 HRGs for 18,499 genes were then known to be identified in the GSE72094 dataset (n = 386). In the training group (GSE72094), 197 HRGs were used as variables to construct prognostic features and build a risk prediction signature. Univariate COX model regression analysis was performed to find a significant association with OS in LUAD patients. There were 62 survival-related HRGs in the GSE72094 dataset.

Identification of hypoxic subtypes of LUAD cohorts

The four LUAD cohorts were grouped into hypoxic subtypes based on the ex-pression levels of 62 survival-related HRGs using the unsupervised consensus clus-tering method. The analysis yielded that all four LUAD cohorts showed the best dif-ferentiation of clustering groups with K = 2 (Fig. 2A, Figs. S1A, S1C and S1E), with 189 patients belonging to hypoxia-related Cluster 1 and 314 patients in hypox-ia-related Cluster 2 of LUAD. The prognostic analysis of the TCGA-LUAD cohort assessed the survival differences between the two hypoxic subtypes. Kaplan −Meier curves showed significantly higher OS in Cluster 1 than in Cluster 2 (log-rank test, P < 0.0001; Fig. 2B). There were also significant differences in clinical characteristics such as tumor stage, Child grade, and histologic grade between the two groups. Similarly, the other three LUAD cohorts from the GEO database also showed the same prognostic outcome, i.e., Cluster 1 had a significantly higher survival advantage than Cluster 2 (log-rank test, P < 0.0001; Figs. S1B, S1D and S1F). The clustering of hypoxic subtypes based on the “Rtsne” package was subjected to t-SNE dimensionality reduction analysis to further verify that the two groups of hy-poxic subtypes identified based on Consensus Clusterin could be well differentiated (Fig. 2C). Subsequently, a Sankey diagram was created based on two hypoxia-related subtypes of LUAD to assess the distribution of clinical information in the hypoxic subgroups of the TCGA cohort. The proportion of patients with pathological stage 1 and alive was higher in Cluster 1 than in Cluster 2 (Fig. 2D). Prognostic analysis and clinical information analysis based on clustering groupings showed that the expression of HRGs interfered with the development of LUAD tumors.

Immunological characteristics and clinical features of the hypoxic subtypes

Immune infiltration analysis and comparison of the two hypoxic subtypes in the TCGA cohort. GSVA was used to explore the functional heterogeneity of metabolic pathways between the two hypoxic subgroups in LUAD, showing that Cluster 1 was significantly enriched in immune and inflammatory pathways such as complement and coagulation cascades and systemic lupus erythematosus (Fig. 3A), indicating that the enrichment of biological pathways associated with immune acti-vation influencing patients in Cluster 1 had a better prognostic outcome. Boxplot results showed that Cluster 1 had significantly higher ESTIMATE scores, immune scores, and stromal scores than Cluster 2 (Mann–Whitney U test, P < 0.001, Fig. 3B (1-3)) in TCGA cohort. TIDE analysis showed that the TIDE score was significantly lower in Cluster 1 than in Cluster 2 (P = 0.0028; Fig. 3B(4)), indicating that hypoxic Cluster 1 was more receptive to the effects of immunothera-py than Cluster 2. However, TMB analysis between the two hypoxic subgroups showed a higher TMB score in Cluster 1 than Cluster 2 (P < 0.0001; Fig. 3B(5)), predicting a better prognosis for patients with Cluster 1 who received immunotherapy. These studies demonstrated the impact of subtype classification based on the definition of hypoxia on the immune microenvironment and immune cell composition resulting from LUAD.

ssGSEA was used to quantify the differences in the abundances of 28 tumor-infiltrating immune cell types between different hypoxic subgroups. The ex-pression of immune cells such as activated CD4 T cells and myeloid-derived sup-pressor cells was significantly different between the hypoxic Cluster 1 and Cluster 2 subtypes (Fig. 3C). Among them, activated CD4 T cells and central memory CD8 T cells were signifi-cantly enriched in Cluster 1, while mast cells were significantly enriched in Cluster 2. Furthermore, analysis of the clinical characteristics of the patients showed statistically significant differences between the hypoxic subtypes (Fig. 3D). It was shown that the two subtypes differentiated by hypoxic features differ significantly in immune microenvironment and immune cell composition, leading to differences in clinical characteristics.

Role of hypoxia-related genes in predicting the prognosis of lung adenocarcinoma

To further analyze the impact of hypoxia genes on clinical prognosis, prognos-tic models were constructed based on patient mRNA data and clinical information. Sixty-two HRGs were analyzed using LASSO regression, and six key features were obtained (Figs. 4A and 4B). A risk profile of survival-related HRGs consisting of six genes (HK1, PDK3, PFKL, SLC2A1, STC1, and XPNPEP1) was constructed by multivariate Cox model regression analysis (Fig. 4C). Univariate and multivariate Cox model regression analyses were conducted for the six genes that ultimately constituted the risk signature (Table 1). The signatures used to calculate the risk score were as follows: ${\displaystyle \text{risk score}=\left(0.7262630\text{*HK}1\right)+\left(-0.6109597\text{*PDK}3\right)+\left(-0.7375822\text{*PFKL}\right)+\left(0.3149958\text{*SLC}2\mathrm{A}1\right)+\left(0.2331100\text{*STC}1\right)+\left(0.5576957\text{*XPNPEP}1\right).}$

The regression coefficients for the four HRGs, namely, HK1, SLC2A1, STC1, and XPNPEP1, were all positive, which indicated that they are related to poor prognosis. In contrast, PDK3 and PFKL were opposite protective factors in the prognostic signature.

For the training set GSE72094, univariate Cox model regression analysis among the cohort subgroups was performed to assess further the ability of risk characteristics (Fig. 4D). The results showed that risk scores were significantly associated with survival information in most subgroups of GSE72094 patients.

Training and validation of the survival prediction ability of six HRG features in the dataset

According to the median risk score based on the six-HRG signature, the cases in the training dataset GSE72094 were classified into high-risk (n = 193) or low-risk (n = 193) groups. Figure 5 shows the prognostic characteristics of the cohort. The Kaplan–Meier analysis curve showed significantly lower OS in the high-risk group than in the low-risk group (log-rank test, P < 0.0001; Fig. 5A). According to the results of univariate and multivariate Cox model analyses in multi-ple datasets, only the HRG signature was a significant prognostic risk factor (univariate Cox model P value <0.05 & multivariate Cox model P value <0.05), whereas the HRG signature was an independent risk factor for survival prediction (Table S4). Then, a time-dependent ROC analysis was established to evaluate the prognostic ac-curacy of the six-HRG signature. In the GSE72094 dataset, the AUCs for 1-, 3-, and 5-year survival were 0.66, 0.72, and 0.62, respectively (Fig. 5B).

Three datasets were used to validate the survival prediction capability of the six-HRG signature, including GSE30219, GSE31210, and TCGA-LUAD (Figs. 5E–5P). KM analysis was conducted in verification sets, confirming the distinct difference in survival rate between the risk groups (Figs. 5E, 5I and 5M). The results of time-dependent ROC for the three validation cohorts are presented in Figs. 5F, 5J and 5N.

DCA of HRG signatures for 3-year OS in four LUAD clinical cohorts was per-formed to determine the nomograms’ clinical significance. The risk of death resulted in a significant positive net benefit (Figs. 5C, 5G, 5K and 5O). The calibration curves all showed high agreement between the observed and pre-dicted results (Figs. 5D, 5H, 5L and 5P). The C-indices of the 4 cohorts are shown in Table S5.

The heatmap introduced the risk score of six-HRG expression, while the dot plot showed the survival status in the training dataset GSE72094. The results showed that the expression of HK1, SLC2A1, STC1, and XPNPEP1 was higher in death-prone patients with high risk scores (Fig. 6A).

Associations of the risk model with overall survival and clinicopathological characteristics of patients with LUAD

In order to investigate the relationship between the risk model and clinicopathological characteristics of LUAD patients, we conducted a Wilcoxon test on the differences in risk scores among subgroups categorized by age, gender, pathological stage, pathological T stage, pathological N stage, and pathological M stage in the TCGA-LUAD cohort. The results showed that the risk score was significantly correlated with age (P value <0.05), gender (P value <0.001), pathological stage (P value <0.05), pathological T stage (P value <0.01), and pathological N stage (P value <0.05), but not with pathological M stage (Fig. S2).

We further conducted stratified survival analysis on multiple LUAD subtypes categorized by age (<70 and ≥ 70 years), gender (Male and Female), tumor pathological stage (stage I–II and III-IV), pathological T stage (T1-2 and T3-4), pathological N stage (N0 and N1-3), and pathological M stage (M0 and M+), in order to assess the predictive performance of the six-HRG risk model for prognosis (Log-rank test, P value <0.05; Figs. S3A–S3F and S3H–S3L). Results demonstrated that patients in the high-risk group exhibited significantly poorer overall survival than those in the low-risk group across all LUAD subtypes, except for the pathological T1-2 subtype. These findings suggest that the six-HRG risk model is a reliable predictor for LUAD prognosis.

Analysis of tumor mutation burden in high- and low-risk groups in LUAD

In the LUAD-TCGA cohort, mutations in the TP53, TTN, MUC16, CSMD3, LRP1B, RYR2, USH2A, ZFHX4, and KRAS genes were present in the top 10 mutat-ed genes in both the high- and low-risk groups of clinical patients and exceeded 20% in the TCGA patient cohort (Figs. 6B and 6C). This indicates that the frequency of mutated genes is highly similar between the two risk groups.

Estimate of tumor-infiltrating immune cells and ICI analysis

The MCP-counter algorithm was used to assess the expression of various tumor-infiltrating immune cells (TIICs) in the LUAD-TCGA cohort and to compare the differential expression of TIICs between the high- and low-risk groups. The immune infiltration of tumor cells is associated with the prognosis of LUAD (Li et al., 2021). The Wilcoxon rank sum test showed that B cells, cancer-associated fibroblasts (CAFs), myeloid cells, and T cells were significantly different between the two risk groups. CAF infiltration was higher in the high-risk group. In contrast, B cells, myeloid cells, and T cells were more infiltrated in the low-risk group ( P value <0.05, Fig. 7A). Immunoinfiltration analysis was performed in three GEO datasets, including GSE72094 (Fig. S4A), GSE31210 (Fig. S4B), and GSE30219 (Fig. S4C), to evaluate the expression of immune cells in different risk groups. Results showed that CAFs, myeloid cells, and T cells exhibited statistical differences (P value <0.05) between risk groups in all three GEO datasets, and their expression patterns in risk groups were consistent with those in the TCGA-LUAD dataset. B cells showed significantly lower expression levels in the high-risk group compared to the low-risk group in GSE72094 and GSE31210 datasets (P value <0.05).

The xCell algorithm was used to analyze the immune score, stromal score, and tumor microenvironment score of the LUAD-TCGA cohort. Pearson’s correlation analysis was performed between the risk score and TMB, immune score, stromal score, and microenvironment score. A significant positive correlation between risk scores and TMB scores was shown (Fig. 7B). The results of the correlation analysis showed a statistically negative correlation between the risk score and all three immune infiltration scores (Figs. 7C–7E). The analyses suggest that the low-risk group exhibited higher immune activity, which may lead to a better prognosis and further validate the LUAD prognostic model based on HRGs in immunotherapy.

The TIDE score was used to analyze the effects of ICIs on the disease. The re-sults indicated that the high- and low-risk groups were less likely to evade tumor immunization by Wilcoxon test. Statistical analysis showed that the TIDE score was significantly lower in the high-risk group (P value <0.05), indicating that patients in the high-risk group may have a better therapeutic effect after receiving immunotherapy. (Fig. 7F(1)).

GSEA

To explore the potential mechanisms in different prognostic groups, GSEA between risk groups in LUAD-TCGA was performed. Based on the GSEA, 11 path-ways were significantly enriched (Table S6). “GPI-anchor biosynthesis” and “protein export” pathways were expressed in a more activated state in the high-risk group than in the low-risk group. In contrast, “Focal adhesion”, “Ribosome” and “Viral carcinogenesis” were expressed in a relatively inhibited state in the high-risk group. The top five pathways are shown in the ridgeplot and Gene-Concept Network (Fig. 7G).

Evaluation of candidate targeted drugs

Drugs approved for the treatment of NSCLC were collected from the DrugBank database (Sun et al., 2022) as targeted drugs. Using the Wilcoxon statistical analysis, we found that the high-risk group had significantly lower IC50 values for Docetaxel, Erlotinib, Gefitinib, Paclitaxel, and Vinorelbine compared to the low-risk group (P value <0.05) (Fig. 7F(2-6)). Cmap scores of molecule drugs in the A549 cell line were shown through the CMap database. Calculate the CMap scores for the five drugs mentioned above based on their IC50 results. The scores of Paclitaxel and Vinorelbine revealed their significant association and more substantial reversal (Table S7). Results suggest that the risk model is sensitive to targeted therapy.

Survival analysis and pathological atlas on single six-HRGs

Kaplan–Meier analysis of six individual HRGs was performed sequentially in GSE72094 to further explore the mechanism of each specific gene (Figs. 8A–8F). Based on the gene expression of six-HRGs, STC1, PFKL, HK1, PDK3, SLC2A1, and XPNPEP1, there was a significant difference in OS between patients in the two risk groups in the training set (log-rank test, P < 0.0001).

Experimental validation on single six-HRGs

Total RNA was extracted from A549 cells and then subjected to RT −PCR experiments. The amplified products were subjected to agarose gel electrophoresis, and six samples all gave clear bands in the correct region indicated by the marker, confirming that the six predicted six-HRGs (STC1, PFKL, HK1, PDK3, SLC2A1, XPNPEP1) were all expressed in lung adenocarcinoma A549 cells (Fig. 8G).

Correlation analysis was performed in TCGA-LUAD dataset using the gene expression levels of six-HRGs. The results, as depicted in Fig. S5, indicated a significant positive correlation between SLC2A1 and XPNPEP1 and the expression levels of other hypoxia-related genes. Two genes listed as risk factors, SLC2A1 and XPNPEP1, were selected for cell transfection and subsequent phenotypic verification experiments. After cell transfection, total RNA was extracted, and the RT–PCR experiment demonstrated that siRNA against both SLC2A1 and XPNPEP1 genes showed knock-down efficiency after transfection (Fig. 8H). Analysis showed that the OD values of A549 cells with exogenous knock-down of the poor prognosis-associated genes SLC2A1 or XPNPEP1 were significantly lower than those of normal control cells after 5 days of culture (P < 0.001), indicating a significant decrease in the proliferation of A549 cells. It reflected a significant reduction in the cellular activity of A549 cells with knock-down of the above two genes (Figs. 8I and 8J). Flow cytometric analysis showed that both SLC2A1 and XPNPEP1 significantly inhibited apoptosis in LUAD cells (P < 0.001, Figs. 8K and 8L). In addition, transwell assays were performed to investigate the effects of predicted poor prognostic genes on LUAD cell migration and invasion. LUAD cells knocked-down with SLC2A1 or XPNPEP1 showed significantly lower migration and invasion rates than cells stimulated with culture medium (P < 0.001, Figs. 8M and 8N). Knock-down of SLC2A1 or XPNPEP1 substantially suppressed the proliferation and migration of A549 cells, demonstrating that the expression of these genes promotes the proliferation of lung adenocarcinoma tumors. This provides a valid reference for the search for biomarkers in the clinical treatment of LUAD.

“Scissor” algorithm to identify subpopulations of TAMs associated with the LUAD hypoxic subtype phenotype in single-cell sequencing data

Using “Scissor”, an algorithm that combines bulk and single-cell sequencing data, the classification of M1 and M2 cells in TAM cells is distinguished in the single-cell data based on the phenotypic differences in the hypoxic subtypes. The validity and guidance value of the identified hypoxic subtypes were assessed. The marker gene identification method was used to differentiate TAMs, where TSPO was used to define M1 macrophages and CD163 was used to define M2 macrophages (Yuan et al., 2022). UMAP shows the distribution of expression of the two marker genes, TSPO and CD163, in TAMs in single-cell data, thus identifying the M1 and M2 cells (Fig. 9A). Scissor analysis of lung cancer single-cell data was performed using the hypoxic subtype in TCGA-LUAD bulk cell samples. Scissor+ cells were predefined as having a hypoxic Cluster 1 presentation profile, and Scissor cells were predefined as having a hypoxic Cluster 2 presentation profile. UMAP plots were viewed for the distribu-tion of Scissor cells, and out of 12,871 cells from different cell types, Scissor se-lected 351 Scissor+ cells and 377 Scissor- cells, which had a high confidence rela-tionship with the M1 and M2 cell classifications in TAMs (Fig. 9B). The analysis showed that the subpopulation of cells closely associated with the hypoxic Cluster 1 phenotype (with a good prognostic effect), was predominantly distributed among M1 cells in single cells, while the subpopulation of cells closely associated with the hypoxic Cluster 2 phenotype (with a poor prognostic effect) was mainly concentrated in M2 cells. Functional enrichment analysis also confirmed that pathways including glycolysis, glucose metabolism, and hypoxia-related pathways were significantly activated in Scissor+ cells (Fig. 9C).