The canonical SMILES representation of PFOS was obtained from the PubChem database. To systematically identify its potential protein targets, we interrogated three independent toxicological databases: ChEMBL, STITCH, and SwissTargetPrediction, with the organism filter set to Homo sapiens. To ensure the reliability of the identified chemical–protein interactions, specific confidence thresholds were applied: for STITCH, only targets with a “combined score” ≥ 0.4 (medium confidence) were retained; for SwissTargetPrediction, targets with a non-zero probability score were included. This multi-database consensus strategy served as an initial reliability filter, prioritizing targets substantiated across different predictive algorithms. The resulting candidate target lists from these resources were consolidated, and duplicate entries were removed to generate a non-redundant set of PFOS-associated molecular targets.

To ensure a comprehensive compilation of NSCLC-associated genes, we retrieved candidate genes from two complementary databases: the Online Mendelian Inheritance in Man (OMIM), an authoritative resource for curated gene-disease relationships, and GeneCards, an integrative database that aggregates information from multiple sources. This dual-database approach was designed to maximize the coverage of potentially relevant genes. The gene sets from both sources were merged and subjected to identifier standardization. Following the removal of duplicates, a unique, consolidated list of NSCLC-associated genes was established for subsequent intersection analysis with the PFOS targets.

To delineate the biological roles and pathways associated with the shared PFOS-NSCLC targets, we performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses using the DAVID bioinformatics platform (13). The analysis comprehensively covered biological processes (BP), cellular components (CC), and molecular functions (MF). Terms with a false discovery rate (FDR)-adjusted p-value <0.05 were deemed statistically significant.

To identify robust diagnostic biomarkers from the identified targets, we employed two complementary machine learning algorithms targeting the 41 overlapping PFOS-NSCLC genes within the GSE33532 transcriptomic dataset. The machine learning workflow was designed and reported in alignment with the DOME (Data, Optimization, Model, and Evaluation) recommendations to ensure transparency and reproducibility in supervised ML for the life sciences. The least absolute shrinkage and selection operator (LASSO) regression was implemented using the “glmnet” R package. Genes retaining non-zero coefficients in the majority of iterations were identified as LASSO-selected candidates. In parallel, support vector machine-recursive feature elimination (SVM-RFE) was conducted using the “caret” and “mlbench” packages after feature standardization. An SVM model with a radial basis function (RBF) kernel was trained, utilizing 10-fold cross-validation to rank features based on their predictive contribution. Genes consistently selected by both LASSO and SVM-RFE were defined as the final core targets for downstream validation (14), a strategy designed to minimize overfitting and enhance reproducibility.

Gene expression data from the GSE33532 cohort, comprising paired NSCLC and adjacent normal tissues, were obtained from the Gene Expression Omnibus (GEO) (15) database. Differential expression analysis of the core genes was performed using the Wilcoxon matched-pairs signed-rank test, with the Benjamini–Hochberg procedure applied to control the false discovery rate (FDR) for multiple comparisons. The diagnostic performance of each candidate gene was evaluated by constructing receiver operating characteristic (ROC) curves and calculating the area under the curve (AUC).

The composition of immune cell infiltrates in NSCLC and normal tissues was characterized using the CIBERSORT algorithm. This deconvolution method estimates the relative proportions of 22 distinct immune cell types based on the LM22 signature matrix, a validated leukocyte gene signature file. The analysis was performed with 1,000 permutations to assess the statistical significance of the deconvolution results. Only samples with a CIBERSORT p-value <0.05 were considered to have accurate immune composition estimates and were retained for downstream analyses. Differences in the abundance of immune cells between groups were assessed using the Wilcoxon rank-sum test. To account for multiple testing across 22 cell types, p-values were adjusted using the Benjamini–Hochberg (FDR) method, with an FDR < 0.05 considered statistically significant. Furthermore, Spearman’s rank correlation analysis was employed to investigate the associations between the expression levels of the core genes and the inferred immune cell fractions (16). Correlation analyses across gene–cell pairs were also subjected to FDR correction, highlighting only robust associations (FDR < 0.05).

A two-sample MR framework was adopted to infer a potential causal relationship between genetically predicted core gene expression and NSCLC risk (17). Genetic instruments for EIF4EBP1 expression were derived from cis-expression quantitative trait loci (cis-eQTLs) identified by the eQTLGen Consortium (sample size, n = 31,684). Summary-level association statistics for the outcomes—lung adenocarcinoma (LUAD) and lung squamous cell carcinoma (LUSC)—were sourced from the FinnGen consortium R12 release.

Genetic variants located within a ±100 kb cis-window of the EIF4EBP1 gene locus and significantly associated with its expression (p < 5 × 10⁻⁸) were considered as candidate instruments, with a minor allele frequency threshold of >1%. Linkage disequilibrium (LD) clumping was performed using a European reference panel from the 1,000 Genomes Project (clumping window = 10 Mb, r² threshold <0.01) to ensure instrument independence and reduce residual correlation among instruments. Palindromic SNPs and those unavailable in the outcome datasets were excluded. The strength of each instrumental variable was quantified using the F-statistic, with values >10 indicating a low risk of weak instrument bias. Steiger filtering was applied to verify that the instruments explained more variance in the exposure (EIF4EBP1 expression) than in the outcome (NSCLC).

The primary causal estimate for the analysis involving a single genetic instrument was calculated using the Wald ratio method. For scenarios with multiple valid instruments, the inverse-variance weighted (IVW) method would have been employed as the primary analysis, supplemented by sensitivity analyses including MR-Egger, weighted median, and MR-PRESSO to assess pleiotropy and robustness.

The high-resolution three-dimensional crystal structures of the core target proteins were acquired from the RCSB Protein Data Bank (PDB). The three-dimensional structure of PFOS was downloaded from PubChem in SDF format and converted to Mol2 format using Open Babel software. Protein structures were prepared for docking by removing all crystallographic water molecules and adding polar hydrogen atoms to optimize protonation states and hydrogen bonding networks. Molecular docking simulations were performed using AutoDock Vina. The grid box parameters for the binding site were defined, and all necessary input files were prepared using AutoDockTools (v1.5.7). The resulting protein–ligand complexes were visualized and rendered using PyMOL (v3.2), and the conformation with the most favorable (lowest) binding energy was selected for detailed interaction analysis (18).

An AOP was systematically developed to delineate the mechanistic sequence linking PFOS exposure to NSCLC pathogenesis. Potential key events (KEs) were identified by analyzing target–phenotype and gene–phenotype networks. Genes substantiated as shared targets, differentially expressed, and directly interacting with PFOS in docking simulations were designated as candidate molecular initiating events (MIEs). The final AOP was constructed by integrating our computational findings with established literature, mapping a plausible causal trajectory from PFOS exposure (MIE) through intermediate KEs (EIF4EBP1-centered translational reprogramming and immune microenvironment remodeling) to the adverse outcome of NSCLC.

To delineate the molecular landscape linking PFOS exposure to NSCLC pathogenesis, we employed an integrated network toxicology approach. Screening of the ChEMBL, SwissTargetPrediction, and STITCH databases identified 256 putative PFOS targets in Homo sapiens (Supplementary Table S1). Concurrently, to establish a robust set of disease-related genes, we compiled NSCLC-associated genes from both the OMIM and GeneCards databases. The merged and deduplicated list contained 2,173 non-redundant genes (Supplementary Table S2). Cross-referencing these datasets yielded 41 overlapping targets, representing high-confidence candidate genes potentially involved in PFOS–NSCLC intersections (Figure 1). To ensure that these overlaps were not merely artifacts of broad NSCLC oncogenic signaling, we verified that each of the 41 genes exhibited documented toxicological associations with PFOS within the chemical–protein interaction databases. Unlike generic NSCLC-associated gene sets, this specific cohort was significantly enriched in PPAR signaling and lipid-sensing modules (Section 3.2), which are canonical hallmarks of PFAS-mediated metabolic disruption. This dual-criteria approach—requiring both evidence of chemical interaction and disease relevance—was implemented to differentiate PFOS-specific mediators from general cancer background pathways and to mitigate the risk of overinterpretation.

The GO and KEGG pathway analyses were performed on the 41 shared targets to elucidate their biological significance (Supplementary Tables S3, S4). The GO annotation revealed a predominant enrichment in processes related to intracellular receptor signaling, hormone-mediated responses, and transcriptional regulation (Figure 2A). CC analysis localized these targets primarily to the plasma membrane (particularly the external and apical sides) and within nuclear transcriptional complexes, including RNA polymerase II and histone deacetylase complexes (Figure 2B). MF analysis confirmed their core roles in nuclear receptor activity, ligand-activated transcription factor binding, and histone deacetylase activity (Figure 2C).

The KEGG pathway analysis further underscored significant enrichment in the PPAR signaling pathway, chemical carcinogenesis via receptor activation, and thyroid hormone signaling (Figure 3). Additionally, pathways such as transcriptional misregulation in cancer were also prominently featured. Notably, the mapped genes included key regulators such as CD36, PPARA/G, RXRA/B, VDR, ABCB1, HDAC1-3, and EIF4EBP1. These results collectively suggest that the convergence of PFOS toxicity and NSCLC etiology involves systemic perturbations in lipid-sensing nuclear receptors, epigenetic modulation, and xenobiotic transport mechanisms.

To pinpoint the most critical diagnostic mediators among the overlapping targets, we applied two machine learning algorithms to NSCLC transcriptomic data. LASSO regression, utilizing 10-fold cross-validation to mitigate overfitting, identified five candidate genes (Figures 4A,B). Simultaneously, SVM-RFE selection isolated three key features (Figure 4C). The intersection of these independent selection methods robustly highlighted three core targets: CD36, ABCB1, and EIF4EBP1 (Figure 4D), which were prioritized for subsequent validation.

We first validated the expression profiles of the identified core genes in the GSE33532 discovery cohort. Consistent with the network toxicology screening, EIF4EBP1 was significantly upregulated in NSCLC tissues, whereas CD36 and ABCB1 exhibited marked downregulation (Figure 5A). ROC curve analysis confirmed the robust internal discriminatory performance of the core gene signature, with all markers achieving AUC values >0.90 (Figure 5B). Specifically, EIF4EBP1, a central hub in our integrated analysis, yielded an AUC of 0.936. While these metrics indicate high discriminatory potential within the discovery dataset, they are interpreted as exploratory estimates of diagnostic utility pending independent external validation.

To substantiate the generalizability of these findings, we interrogated the larger, independent TCGA-NSCLC cohorts. Remarkably, the expression patterns observed in the discovery cohort were fully recapitulated in the TCGA datasets: EIF4EBP1 was significantly elevated across both LUAD and LUSC (Figure 6A). Furthermore, survival analysis in the TCGA cohort revealed a striking, subtype-specific prognostic divergence for EIF4EBP1. High expression was significantly associated with diminished overall survival in LUAD (HR = 1.9, 95% CI: 1.2–2.9, p = 0.0052; Figure 6B), yet correlated with a paradoxical favorable prognosis in LUSC (HR = 0.67, p = 0.047; Figure 6C).

This cross-cohort consistency in expression, coupled with the subtype-specific clinical impact, reinforces the biological significance of our prioritized targets. Among the three core genes, EIF4EBP1 was designated as the primary focus for downstream mechanistic modeling. This prioritization was strategically motivated by its consistent upregulation across platforms (suggesting an oncogenic driver role) and its central position in the mTOR-mediated translational control pathway. Furthermore, its unique structural affinity and subtype-specific causal association identified in Mendelian randomization designated EIF4EBP1 as the most plausible molecular initiating event (MIE) candidate within the PFOS–NSCLC axis.

To explore the immune microenvironment, we deconvoluted the abundance of 22 immune cell subsets using the CIBERSORT algorithm. Compared to normal lung tissue, NSCLC samples displayed a marked shift toward an adaptive immune profile, characterized by significantly increased infiltration of plasma cells, activated CD4 + memory T cells, follicular helper T cells, γδ T cells, and M1 macrophages. Conversely, innate effectors such as monocytes, NK cells, neutrophils, and eosinophils were significantly depleted (all p < 0.05; Figures 7A,B). Co-infiltration analysis identified two distinct, anticorrelated immune modules: a humoral/T-helper module (plasma cells and activated CD4 + memory T cells) and a myeloid-granulocytic module (monocytes and neutrophils), reflecting divergent immune polarization states within the TME (Figure 7C). These immune shifts primarily reflect tumor-associated remodeling; PFOS relevance is inferred from target overlap and EIF4EBP1-associated patterns rather than direct exposure measurements.

We further examined the specific association between the core genes and these immune subsets (Supplementary Table S5). While CD36 and ABCB1 generally correlated with myeloid and resting lymphocyte subsets, EIF4EBP1 exhibited a strong positive correlation with adaptive immune components—specifically plasma cells, activated CD4 + memory T cells, and M1 macrophages—and a negative correlation with monocytes and M2 macrophages (Figure 7D). This distinct pattern suggests that EIF4EBP1 upregulation in the tumor milieu is closely linked to an adaptive-skewed, yet potentially remodeled and dysregulated, immune microenvironment.

We employed two-sample MR to investigate the causal relationship between genetically predicted EIF4EBP1 expression and NSCLC risk using summary statistics from the FinnGen R12 release (Supplementary Table S6). Using rs28565141 as a valid genetic instrument (F-statistic = 28.622), we found that genetically predicted higher expression of EIF4EBP1 was causally associated with an increased risk of LUAD (β_MR = 1.434, SE = 0.635, p = 0.024; OR = 4.196, 95% CI: 1.209–14.565) (Figure 8A). In contrast, no significant causal effect was observed for LUSC (β_MR = −0.130, SE = 0.692, p = 0.851; OR = 0.878, 95% CI: 0.226–3.412) (Figure 8B). These findings provide genetic support consistent with a subtype-specific association, implicating EIF4EBP1 as a candidate contributor to LUAD risk.

To assess the structural basis of the direct PFOS-EIF4EBP1 interaction, we performed molecular docking simulations. PFOS was predicted to occupy a shallow groove on the EIF4EBP1 surface, adopting an extended conformation stabilized by a mixture of hydrophobic and polar residues (Figure 9A).

Detailed analysis of the binding interface revealed a network of stabilizing interactions (Figures 9B, 10). The terminal polar group of PFOS forms critical hydrogen bonds with the side-chain amide of Gln187, the carboxylate of Asp44A, and the amine of Lys88A (H-bond distances 1.9–3.1 Å). Notably, the interaction with Lys88A exhibits partial salt-bridge character. Beyond polar contacts, the hydrophobic fluoroalkyl tail of PFOS engages in extensive van der Waals interactions and parallel π–π stacking with the indole rings of Trp43A and Trp89A. This cooperative network of hydrogen bonds, hydrophobic contacts, and stacking interactions suggests that PFOS can form a stable, non-covalent complex with EIF4EBP1, potentially modulating its functional state.

Integrating the findings from multi-omic analyses, we constructed a putative AOP linking chronic PFOS exposure to NSCLC pathogenesis. In this framework, the molecular initiating event (MIE) involves the direct binding of PFOS to core targets such as EIF4EBP1. This interaction is hypothesized to elicit downstream key events (KEs), potentially culminating in NSCLC initiation/progression and therapy resistance. These molecular and pathway-level changes propagate to cellular responses, specifically the remodeling of the tumor immune microenvironment toward a dysregulated adaptive state. Collectively, this cascade of events culminates in the adverse outcome (AO) of NSCLC initiation, progression, and potential therapy resistance (Figure 11).