The RNAseq data (GSE135917 and GSE38792) related to OSAS were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/) (27, 28). Both of their sequencing platforms are GPL6244, total RNA was isolated from human subcutaneous fat. The first group of GSE135917 consisted of 10 OSAS patients and 8 normal people, whereas the second group consisted of 24 OSAS patients who were sampled individually after treatment with a continuous positive airway pressure (CPAP) ventilator. In brief, CPAP was started following the biopsy and after 2 weeks of self-reported CPAP use for more than 4 h per night. The follow-up biopsy was conducted on the other side of the abdomen.

The ReadAffy function from the affy package (version 1.72.0) was utilized to read the raw data included in the cell format files. The data were imported into SIMCA 14.1 (Sartorius, Malmö, Sweden) for analysis, and four components were chosen to construct the PLS-DA model (29, 30).

The raw data are normalized using the Affymetrix platform in R software (version 4.0.5). The RMA functions were then used to process data using the affy package. The Limma package (version 3.50.3) was used to perform DEGs between OSAS and normal groups with the empirical bayesian t-test, and |log2-Fold change (FC)| >0.5 and P < 0.05 were utilized as DEG screening criteria. The heatmap package was used to create a heat map that displayed the top 50 genes with the most significant genes. The ferroptosis-related genes were downloaded from the FerrDb V2 database (http://www.zhounan.org/ferrdb/), which mainly contained gene sets such as “Driver,” “Suppressor,” and “Marker.” A Venn diagram was utilized to illustrate the intersection genes of FerrDb V2 and DEGs (f-DEGs) (31).

The genes in the dataset were ranked by OSAS phenotype and logFC value, and Gene Set Enrichment Analysis (GSEA) was performed on 664 gene sets from WikiPathways. The statistical differences were determined by the normalized enrichment score (|NES| >1), and FDR < 0.25. The dataset was read and subjected to Gene set variation analysis(GSVA) analysis using the GSVA package (version 1.42.0) to acquire a GSVA score for each sample. The Limma package was used to compare the GSVA scores between the OSAS and the normal group with the Bayesian t-test, adj. P < 0.05 and |logFC| > 0.1 represent significant differences.

The f-DEGs were transformed with the org.Hs.eg.db package (version 3.10.0) and then enriched with the cluster profiler package (version 4.2.2) for gene ontology, disease ontology, and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis (32, 33). The q value ≤ 0.05 and gene counts ≥3 were considered significant.

The f-DEGs were analyzed using the STRING database (https://cn.string-db.org/) and a combined interaction score>0.4 was considered statistically significant. Subsequently, the string interactions were visualized using Cytoscape software (version 3.9.0). CytoNCA (2.1.6), a Cytoscape plug-in, was used to filter the density and importance of modules in the PPI networks. Then, the R software screened out the hub genes in the network, based on the following principles: betweenness, closeness, and degree are greater than average.

The THCA gene expression profiles were downloaded from the TCGA database. There were 568 cases of TCGA (510 cases in the tumor group and 58 cases in the normal group). The online database ULCAN (http://ualcan.path.uab.edu) analyzed the gene expression level of HIF1A in cancer tissues and normal tissues of TCGA database samples. Validation of the effect of HIF1A expression on THCA survival through the Kaplan-Meier plotter (http://kmplot.com/analysis/). Transcriptome data for THCA were grouped according to pathology type or gene expression. CIBERSORT was used to compute immune cell composition based on gene expression profiles. R software (including vioplot, ggpubr, and ggExtra packages) was used to analyze the correlation between HIF1A and THCA immune infiltration.

In the GSE135917 dataset, single sample gene set enrichment analysis was used to calculate the per sample infiltration levels of immune cell types. The ggboxplot function in the ggpubr package (version 0.4.0) was used to plot box line plots of immune cells in the two groups of samples. The heatmaps were plotted using the pheatmap package (version 1.0.12), showing the relationship between hub genes and various immune cells.

GSE38792 was utilized to validate the differential expression of hub genes between OSAS patients and healthy volunteers. A variance analysis was performed using the stat_compare_means function in the ggpubr package (version 0.4.0), and the results were presented in box plots (34).

Eight male Sprague-Dawley rats (SPF grade, weight of 230–250 g) were obtained from the central laboratory of the Animal Experimental Center at Renmin Hospital of Wuhan University. The rats were fed in a 12-h alternating day and night environment with corresponding humidity and temperature control. The rats were randomly divided into two groups: (1) Normal control group, in which the rats were kept in a normoxic environment for 8 weeks; (2) Chronic intermittent hypoxia group (CIH), in which the mice were kept in a chronic intermittent hypoxic environment 8 h per day for 8 weeks. The oxygen concentration in the modeling chamber alternated between 30 s of hypoxia (FiO2, 10%) and 60 s of reoxygenation (FiO2, 21%). After 8 weeks, we administered general anesthesia to rats, isolated the abdominal aorta, and obtained arterial blood for arterial blood gas analysis using a blood gas needle.

Total RNA was isolated and purified using RNAiso Plus reagent (9108, Takara, Japan) according to the manufacturer's protocol. Based on the TaqMan probe method, The mRNA expressions of hub genes (HIF1A, ATM, MAPK8, and MAPK14) were detected according to the instructions of probe One-Step qRT-PCR Kit (D7277, Beyotime, China). The reactions were carried out in the cycler under the following conditions: 50°C for 20 min, 95°C for 2 min, 95°C for 15 s and 60°C for 20 s (40 cycles in total). GAPDH was used as a housekeeping gene, and the RT-qPCR primer sequences are listed in Table 1. The relative gene expression level was calculated using the 2– ΔΔCT method (35).

ROS production was detected by in situ staining of fresh lung tissue with a ROS fluorescent probe-dihydroethidium (DHE, D7008, Sigma, 1:500). In this process, sections incubated with stain sections were observed under a fluorescence microscope (IX53, Olympus, Japan). Fluorescein-labeled ROS-positive sections emitted red fluorescence (excitation wavelength 490 nm, emission wavelength 560 nm). At least three sites were selected for each sample. The fluorescence intensity of the cell sections was measured using Image J software.

The lungs were quickly removed after the rats were anesthetized and sacrificed. The saline was rinsed and blotted dry with filter paper. The lung tissue was ground with a cryogenic grinder at 4°C, and the supernatant was centrifuged and assayed according to the manufacturer's instructions for the MDA and GSH kits (A001-3-1/A006-2-1, Nanjing Jiancheng Bioengineering Institute, China). The contents of MDA and GSH were measured in the normal and CIH groups for intra-assay and inter-assay repeated experiments, respectively. The intra- and inter-coefficients of variation were both < 10%.

After the inguinal fat was removed from the rat, the adipose tissue was pressed with oil-absorbing paper (36). The expressions of HIF1A and Atm in rat inguinal adipose tissue were detected using immunohistochemistry. The sections (6 μm) were processed with deparaffinized and rehydration using xylene and different concentration gradients of ethanol. The sections were added to citric acid antigen repair solution (pH 6.0) and then heated in a 95°C water bath for 20 min for antigen repair (37). The primary antibodies used were HIF1A (1:500, K000487P, Solarbio, China), ATM (1:200, K009314P, Solarbio, China). The sections were incubated with HRP-labeled secondary antibodies for 1 h at room temperature. The relative IOD of the immune sections was measured with ImageJ software.

We estimated the area under the ROC curve (AUC) using version 1.8 of the standard pROC tool for the R software. The Pearson correlation coefficient was used to conduct correlation analysis on the hub genes. The CPAP-treated population belongs to the second group in the dataset GSE135917. The paired t-test was performed to compare the variation in hub genes expression before and after the application of a non-invasive ventilator. The DrugBank database was used to search for probe targets (Version 4.2) (38).

GraphPad Prism 8 statistical software (La Jolla, CA, USA) was utilized for statistical analyses. All results were expressed as a mean ± standard deviation (SD) from three independent experiments. The data were tested with the Shapiro-Wilk normality test and the variance homogeneity test. Comparisons were tested using paired or unpaired t-tests and the level of confidence was set at 95% (P < 0.05).

Figure 1A depicts the whole study methodology, including bioinformatics analysis and rat validation model. The animal model of OSAS consists of three components: an intermittent hypoxic chamber, a nitrogen tank, and an oxygen compressor that varies FiO2 between 10 and 21% (Figure 1B).

To evaluate the accuracy of the gene sequencing in the GSE135917 dataset, we employed RNA degradation curve analysis and PLS-DA. The modest slope of each curve shows that mRNA is not degraded explicitly (Supplementary Figure S1A). The 2D and 3D structure of the PLS-DA model is depicted in Supplementary Figures S1B, C. Four components were used to construct the PLS-DA model, which had a cumulative explanatory power of 99.4% and a cumulative predictive power for the dependent variable of 73%. To determine if the PLS-DA model was overfitted, the model was evaluated by holding the X matrix constant and randomly rearranging the Y matrix variables 200 times to generate permutation test results (Supplementary Figure S1D). The results indicate that the normal and OSA groups can be discriminated against without difficulty and that the model is not overfitting.

After GSE135917 Validation, we filtered the DEGs and performed GSEA and GSVA. A total of 23,281 gene expression values were gathered. The OSAS samples contained 1,470 down-regulated genes and 1,674 up-regulated genes, as compared to the normal controls.

The heat map depicts the 50 genes with the greatest differences (positive and negative values) according to logFC values in the dataset (Figure 2A). Ferroptosis, Alzheimer's disease, insulin signaling, mapk signaling pathway, and oxidative stress response were all significantly different between the OSAS and normal groups. The ferroptosis pathway was significantly different between the two groups (NES = 1.8234, P = 0.0005, and FDR = 0.005) (Figure 2B). The WikiPathways were evaluated using gene set variation analysis, and it was determined that there was a significant difference in ferroptosis between the two groups, with elevated sample scores in the OSAS group (Figure 2C). The ferroptosis gene set included 64 genes, and the circular heat map revealed that GCLM, HMGCR, SLC38A1, and CHMP5 were elevated in OSAS (Figure 2D). This grants the theory that ferroptosis may play a significant role in the pathophysiology of OSAS.

To determine whether DEGs are involved in ferroptosis, we retrieved f-DEGs from the FerrDb V2 database and performed functional enrichment on them. The intersection of 543 ferroptosis-related genes from the FerrDb V2 database with DEGs using Venn diagrams generated 74 f-DEGs (Figure 3A). The KEGG pathway database was used to identify 102 enriched pathways, including ferroptosis, endocrine resistance, mitophagy-animal, autophagy-animal, and Kaposi's sarcoma-associated herpesvirus infection (Figure 3B, Supplementary Table S1).

The most enriched pathways in terms of GO terms were those associated with neuron death, regulation of neuron death, cellular response to metal ions, cellular response to chemical stress, and reaction to metallic ions (Figure 3C). On the Metascape website (https://metascape.org/), GO analysis of f-DEG revealed that ferroptosis-related pathways are involved in the regulation of neuronal death (Figure 3D). The f-DEGs were evaluated in the Disease Ontology database, and thyroid carcinoma, papillary thyroid carcinoma, thyroid cancer, and colon cancer were discovered to be enriched disorders (Figure 3E). We discovered that f-DEGs are closely associated with THCA and may play an important role in its onset.

By establishing a PPI network, we were able to identify the hub genes in f-DEGs and establish their relationship. A 71-node, 142-edge PPI network was constructed based on the biological interactions of 74 f-DEGs. The PPI enrichment P-value is 9.9910-16, and the average local clustering coefficient is 0.419. CytoNCA (2.1.6) was applied to analyze f-DEGs, resulting in a network of 53 nodes and 284 edges (Figures 4A, B). Furthermore, the hub genes were discovered using the same screening method: HIF1A, ATM, HSPA5, MAPK14, KRAS, MAPK8, TLR4, and CREB1 (Figure 4C). In the network, HIF1A is the most significant gene. Since f-DEGs are primarily enriched in THCA (Figure 3), we hypothesize that HIF1A may be the principal gene of THCA.

We investigated the association between HIF1A and THCA, which was obtained from the TCGA database. The effect of high HIF1A expression in OSAS on thyroid cancer survival discovered that HIF1A expression in thyroid tissue was also increased, and patients with high HIF1A expression had a lower survival time. As a result, HIF1A overexpression in OSAS may accelerate the occurrence and progression of thyroid cancer (Figures 5A, B). This research examines the relationship between HIF1A, TNM, and thyroid cancer stage to determine the impact of HIF1A on thyroid cancer. HIF1A expression is lowest in normal people and increases with the T and N stages. In the M stage, the 0 stage was highest. HIF1A promotes thyroid carcinoma in patients (Figure 5C). To investigate whether HIF1A promotes THCA, the enrichment pathway was examined. HIF1A is enriched in the ferroptosis, ROS, cancer, and apoptosis pathways of THCA (Figure 5D). The immune milieu plays a crucial role in the development of tumors. HIF1A expression is intimately associated with dendritic cells, macrophages, mast cells, and CD8 T cells in the immunological milieu of THCA (Figures 5E, G). HIF1A causes alterations in the immunological microenvironment of thyroid cancer, which may be associated with ferroptosis of thyroid cells produced by an increase in ROS. OSAS potentially promote the onset and progression of THCA by increasing HIF1A expression.

We performed immune cell analysis on the GSE135917 dataset to understand the role of OSAS on the immune environment of adipocytes. OSAS was associated with an increase in activated CD4 T cells, gamma delta T cells, and regulatory T cells (Supplementary Figure S2A). A correlation heat map was generated based on the immune cell contact relationship (Supplementary Figure S2B). Analysis of the relationship between hub genes and immune cells demonstrated that HIF1A was the most important gene in OSAS immune cells (Supplementary Figures S2C, D). According to Figure 5 and Supplementary Figure S2, we found that HIF1A is critical for the development of THCA and OSAS, which may increase the risk of disease by affecting the immune microenvironment.

We validated the expression of hub genes in another dataset (GSE38792). OSAS had elevated the expression of HIF1A, ATM, HSPA5, MAPK8, MAPK14, TLR4, and CREB1. HIF1A and ATM exhibited statistically significant differences (Figure 6A). HIF1A and Atm mRNA expressions were elevated in rats exposed to CIH, with P-values of 0.037 and 0.002, respectively (Figure 6B). Rearrangement of the OSAS phenotype based on HIF1A expression demonstrated that the ferroptosis pathway remained significantly distinct, NES = 1.4508007, P < 0.05 (Figure 6C). PO2 (CIH group): 51.86 ± 2.97 mmHg, PO2 (Normal group): 69.87 ± 2.45 mmHg. There is a statistical difference between the two groups (P < 0.05).

Figure 7 depicts the values of the oxidative stress markers (ROS, MDA, and GSH) and the hub genes (HIF1A, Atm). CIH exposure increased the expression of ROS in rat lung tissue compared to the control group (Figures 7A, B). MDA was greatly increased in the CIH group compared to the control group, although GSH was dramatically decreased, and the P-values were all significant (P < 0.01) (Figures 7C, D). The expression of HIF1A and Atm proteins in inguinal fat was examined using immunohistochemical techniques in CIH-exposed rats and was significantly elevated in the CIH group with a significant P-value of 0.0049 (Figures 7 E, F). Both Figures 6, 7 demonstrate that HIF1A and ATM play a significant role in the process of OSAS-induced ferroptosis, which is intimately connected to the elevation in ROS levels and the decline in GSH levels.

We performed a diagnostic curve analysis to assess the significance of hub genes (HIF1A and ATM) in the diagnosis of OSAS. The respective areas of HIF1A under the curves (AUC) were 0.838 and 0.80 in both datasets, corresponding to an AUC of 0.838, 0.812 for ATM, respectively (Figures 8A, B). The AUC of the combined HIF1A and ATM model for the diagnosis of OSAS was 0.875 in both datasets, suggesting that the model is beneficial for the identification of OSAS (Figures 8C, D). To assess the interaction analysis between hub genes, Pearson correlation analysis was performed. HIFA and ATM were found to be strongly correlated in both datasets, with correlation coefficients (r) of 0.56, 0.70 and P-values of 0.02, 1.2e-3, respectively (Figures 8E, F).

It was found that there was no significant change in HIF1A and ATM in OSAS patients treated with ventilators, with P-values of 0.74, and 0.31, respectively (Figure 8G). There was no significant change in HIF1A and ATM after short-term non-invasive ventilator (CPAP) therapy, and we investigated pharmacological treatment of OSAS in the following phase.

DrugBank includes the most comprehensive information on medications and their targets, so it provides us with the most important drug targets. Eight drugs targeting these two hub genes were identified from the Drugbank (Table 2). Among these drugs, the modulators are mainly 2-methoxy estradiol, carvedilol, ENMD-1198, and FG-2216, the inducers are hydralazine, and the inhibitors are PX-478 and vadarestat, The stabilizer has caffeine. The approved pharmaceuticals include carvedilol, hydralazine, and caffeine, which are potential OSAS treatments.