Data sources and preprocessing (bulk expression, single-cell, proteomics)

We integrated bulk gene-expression datasets, single-cell RNA sequencing (scRNA-seq) datasets, and discovery proteomics to identify COPD-associated biomarkers and pathways. Bulk gene-expression datasets were retrieved from the Gene Expression Omnibus database, including GSE47460, GSE76925, and GSE37768 [10, 11]. GSE47460 contained two platforms, GPL6480 and GPL14550, including healthy controls and COPD samples (healthy controls n = 17 and COPD n = 75 for GPL6480; healthy controls n = 91 and COPD n = 145 for GPL14550). GSE76925 (COPD n = 111, controls n = 40) and GSE37768 (COPD n = 18, controls n = 9) were used as independent validation cohorts. Probe-level preprocessing removed empty probes and those mapping to multiple genes; for genes with multiple probes, the probe with the highest mean expression was retained, and when multiple probes per gene were required, median aggregation was used. Datasets were merged after re-annotation to gene symbols using the official platform annotation files and were batch-corrected and normalized using ComBat (R package sva) [12]. Subsequent analyses used these normalized matrices. scRNA-seq data were obtained from GSE196341 and GSE135893 (12 COPD and 12 control donors in total) [13]. Discovery proteomics were deposited to ProteomeXchange (iProX partner) under PXD068247 [14, 15]. A summary of datasets is provided in Table 1.

Table 1 Datasets included in the studyDifferentially expressed analysis

Differential expression analysis between COPD and control samples was performed using the limma package. When applicable, group and batch/platform information were included in the model design. Differentially expressed genes were defined as genes with |log2 fold change|> 0.585 and Benjamini–Hochberg adjusted P < 0.05, unless otherwise specified. Volcano plots and heatmaps were generated using ggplot2 and pheatmap [16,17,18].

Weighted gene co-expression network analysis (WGCNA)

WGCNA was performed on the training cohort (GSE47460) using the WGCNA R package. The gene expression matrix was filtered to remove genes with standard deviation < 0.45 across samples, log2-transformed (after adding 1) due to high expression values (> 1000), and normalized using normalizeBetweenArrays from the limma package. Sample quality was assessed with goodSamplesGenes (verbose = 3), removing genes with zero variance or > 10% missing values and samples with > 10% missing values. Outlier samples were detected via hierarchical clustering (average linkage on Euclidean distance) and removed using static tree cutting (cutHeight = 1000, minSize = 10), retaining the main cluster. Clinical traits were binarized (Normal = 1/0, Disease = 0/1 for control/COPD samples). A signed co-expression network was constructed using pairwise Pearson correlations raised to the soft-thresholding power β, selected via pickSoftThreshold (powerVector = 1:20; scale-free topology fit, signed R2 ≥ 0.90; mean connectivity criterion). The adjacency matrix (type = "signed") was converted to a topological overlap matrix (TOM) using TOMsimilarity. Gene modules were identified by hierarchical clustering (average linkage) of 1-TOM distances, followed by dynamic tree cutting (cutreeDynamic; deepSplit = 2, minClusterSize = 50, pamRespectsDendro = FALSE). Modules with eigengene dissimilarity > 0.3 were merged (mergeCloseModules, cutHeight = 0.3). Module eigengenes (first principal component) were correlated (Pearson) with COPD status (Disease trait) to identify COPD-associated modules (Bonferroni-adjusted P < 0.05). The module with the strongest association was selected for intersection with differentially expressed genes [19].

Machine-learning feature selection and classification

Batch-corrected expression matrices were used for machine-learning analysis. Only genes shared between the training and validation cohorts were retained, and expression values were z-score standardized before modeling. The phenotype variable was defined as Type. All analyses were performed in R with a fixed random seed of 123. The workflow included feature selection followed by classifier training. Algorithms included penalized regression, support vector machine, random forest, gradient boosting, XGBoost, partial least-squares-based models, Bayesian additive regression trees, and generalized linear models. A no-preselection strategy using all shared genes was also included. Models retaining two or fewer genes were discarded.

The final model was selected based on training-cohort performance and then locked before external validation. Model locking fixed the selected genes, preprocessing workflow, fitted model object, and prediction function. No additional feature selection, refitting, or parameter tuning was performed using the validation cohort. The random forest model was selected as the final classifier and achieved AUCs of 0.996 and 0.834 in the training and validation cohorts, respectively. ROC curves and AUCs were calculated using the pROC package.

Gene set enrichment analysis (GSEA) and functional annotation

To investigate the biological pathways associated with candidate genes, gene set enrichment analysis was performed using Gene Ontology and Kyoto Encyclopedia of Genes and Genomes gene sets from the Molecular Signatures Database. Over-representation analysis was also performed using clusterProfiler. Enrichment results were visualized using enrichplot and related R packages. Pathways with adjusted P < 0.05 were considered statistically significant, unless otherwise specified [20].

Single-cell RNA sequencing analysis

Single-cell RNA-seq data were processed using Seurat. Cells were filtered according to quality-control criteria, including mitochondrial RNA content, ribosomal RNA content, erythrocyte RNA content, and the number of detected genes. Data normalization and integration were performed using standard Seurat workflows, and batch effects were mitigated using Harmony. Highly variable genes were identified and used for principal component analysis and uniform manifold approximation and projection.

Cell clustering was performed using a shared nearest-neighbor graph-based method. Cell types were annotated manually according to canonical marker genes and the top marker genes of each cluster. Candidate gene expression was mapped across major lung cell populations, with particular attention to airway epithelial cell populations. Because COPD and control single-cell samples were derived from public datasets, disease-group comparisons were interpreted cautiously in the context of potential dataset-specific effects.

Single‑cell network‑based virtual knockout analysis

To explore the potential regulatory consequences of CYP1B1 perturbation at the single-cell network level, virtual knockout analysis was performed using the scTenifoldKnk framework. Integrated and batch-corrected single-cell RNA-seq data were analyzed, and airway secretory cells were selected for downstream analysis based on prior cell-type localization results.

Raw UMI count matrices were extracted and used for network construction. Highly variable genes were identified using the variance-stabilizing transformation method in Seurat. The final expression matrix included CYP1B1 and selected highly variable genes. scTenifoldKnk constructed gene-regulatory networks through repeated cell subsampling and ensemble network inference. A virtual knockout of CYP1B1 was simulated by computationally removing CYP1B1-associated regulatory edges from the inferred networks. The perturbed networks were compared with corresponding wild-type networks to estimate changes in regulatory influence on other genes.

Perturbation magnitude, Z score, nominal P value, adjusted P value, and network distance metrics were calculated. Genes with adjusted P < 0.05 were considered significantly perturbed. The knocked-out gene itself was excluded from downstream ranking. These results were interpreted as computationally inferred regulatory perturbations rather than direct experimental evidence of gene silencing.

Artificial neural network (ANN)

A feed-forward artificial neural network was constructed to classify COPD and control samples using selected candidate genes. Gene-expression signatures were binarized according to median expression values. The network consisted of an input layer, one hidden layer, and a binary output layer. Model training was performed using Rprop + optimization with a fixed random seed. Model performance was evaluated using confusion matrices, class-specific accuracy, and ROC-AUC. Feature contribution was estimated using the Garson algorithm. The artificial neural network analysis was considered complementary to conventional machine-learning models.

Cigarette smoke extract (CSE) preparation

Mainstream smoke from Hongta brand filter-tipped cigarettes (Yunnan, China) was drawn by a peristaltic pump and bubbled through 10 mL DMEM per cigarette. The pH was adjusted to 7.2, and the solution was sterilized using a 0.22-μm filter to yield 10% CSE (v/v). working concentrations of CSE were prepared freshly for each experiment [21].

Cell culture and treatments

The 16HBE human bronchial epithelial cell line was a kind gift from Guangzhou Medical University (Guangzhou, China), with the original source of the cell line from ATCC. The 16HBE cells were cultured in DMEM (Sigma-Aldrich, Louis, USA) + 10% FBS (Sigma-Aldrich, Louis, USA) at 37 °C, 5% CO₂. Cells were exposed to 2.5% CSE for 24 h. Unless otherwise stated, vehicle-treated cells served as controls [21].

Cigarette smoke-exposed (CS) mouse model

Specific pathogen-free C57BL/6 J mice (15–20 g, from Guangzhou University of Chinese Medicine) were housed at 25 °C, 12/12 h light–dark, with food and water ad libitum. The institutional ethics committee approved protocols. Mice were randomized into two groups: control (n = 8) and CS exposed group (n = 8). CS exposed mouse model was induced via whole-body exposure in a 60 × 57 × 100 cm chamber to 9 cigarettes twice daily, 6 days/week, for 6 months; each session lasted ~ 2 h with ≥ 4 h interval. Control mice were exposed to room air for the same duration [22]. Before exposure, mice received sterile saline aerosol (30 min, twice daily). Tissues (lung), bronchoalveolar lavage fluid (BALF), and blood were collected at endpoint.

Lung function testing

Lung mechanics were evaluated using a FinePointe system (Buxco) according to the manufacturer’s instructions, recording airway resistance/compliance indices [23].

Enzyme-linked immunosorbent assays (ELISAs)

Concentrations of IL-1β (Mouse, proteintech, KE10003), IL-6 (Mouse, proteintech, KE10007), IL-8 (Mouse, Abiowell, EM1592), and TNF-α (Mouse, proteintech, KE10002) in BALF and serum were quantified per the manufacturers’ protocols.

Western blotting (WB)

Total protein was extracted from treated 16HBE cells using RIPA lysis buffer (Beyotime, P0013B) supplemented with protease and phosphatase inhibitor cocktails (Cwbio, CW2200). The lysates were incubated on ice and centrifuged at 12,000×g for 15 min at 4 °C. The supernatants were collected, and protein concentrations were determined using a BCA protein assay kit (Bioswamp, BBCAPCK500) according to the manufacturer’s instructions.

Equal amounts of protein were mixed with loading buffer, denatured at 100 ℃ for 10 min, separated by SDS-PAGE, and transferred onto polyvinylidene difluoride membranes (Millopore, ISEQ00010 and IPVH00010). After blocking with QuickBlock™ Blocking Buffer (Beyotime, P0252) for 20 min at room temperature, the membranes were incubated overnight at 4℃ with primary antibodies against CYP1B1 (diagbio, db14044), GPX4 (Selleck, F1580), 4-HNE (Abclonal, A24456), and SOD1 (Selleck, F5059). β-actin (proteintech, 66009-1-Ig) was used as the internal loading control. After washing with TBST, the membranes were incubated with appropriate horseradish peroxidase-conjugated secondary antibodies for 1 h at room temperature. Protein bands were visualized using an enhanced chemiluminescence detection system and quantified using ImageJ software. The relative expression levels of target proteins were normalized to β-actin. At least three independent biological replicates were included for each group.

Immunohistochemistry (IHC), Histology and staining

Lung tissues were fixed in 4% paraformaldehyde (Beyotime, P0099) for 20 min, processed, and embedded in paraffin. Sections (4 μm) were cut and stained with hematoxylin and eosin (H&E), periodic acid–Schiff (PAS/AB-PAS) for mucus/glycoprotein-rich components, and Picrosirius Red for collagen fibers, following standard protocols.

For IHC, sections were deparaffinized, rehydrated, and subjected to antigen retrieval using citrate buffer (pH 6.0) in a microwave for 10 min. Endogenous peroxidase activity was quenched with 3% hydrogen peroxide for 10 min. Sections were blocked with 5% bovine serum albumin (BSA) for 30 min at room temperature and incubated overnight at 4 °C with primary antibodies (CYP1B1 (diagbio, db14044, 1:200)). After washing with PBS, sections were incubated with horseradish peroxidase (HRP)-conjugated secondary antibodies (1:500, Abcam) for 1 h at room temperature. Visualization was performed using DAB (3,3′-diaminobenzidine) substrate. Sections were counterstained with hematoxylin, dehydrated, and mounted.

Images were captured using a light microscope (Olympus BX53), while the scale bars are shown in the lower left corner. Quantification of staining intensity was performed using ImageJ software. For each group, at least three independent biological replicates were analyzed, and image analysis was performed in a blinded manner.

Molecular docking

To explore therapeutic tractability, we performed molecular docking against CYP1B1. Three-dimensional structures for ligands were retrieved (PubChem) and prepared in AutoDockTools 1.5.7; receptor structure (CYP1B1) was obtained from RCSB PDB (www.rcsb.org) and prepared analogously. The grid center and box size were defined around the active site, and docking was performed using the Lamarckian Genetic Algorithm under default parameters. Rigid docking was performed (AutoDockTools 1.6.7), and the top-scoring pose was retained.

Statistical analysis

Analyses were performed at R 4.4.2. Unless otherwise stated, two-sided tests were used, with α = 0.05. For multiple comparisons, the Benjamini–Hochberg correction was applied. Classification performance was summarized by ROC-AUC (95% CI), PR-AUC, accuracy, F1 score, confusion matrices, and calibration curves. Data are reported as mean ± SD (or median [IQR]) with n indicating biological replicates. Due to the relatively small sample size in animal models, correlations between continuous variables were evaluated using the non-parametric Spearman's rank correlation test to ensure robustness.