Integrated Machine Learning Approaches to Explore the Role of Glycosylation-Related Genes in Idiopathic Pulmonary Fibrosis

Summary

This study used integrated transcriptomic analysis, machine learning (LASSO, SVM-RFE, and XGBoost), immune infiltration analysis, molecular docking, and RT-qPCR validation to identify glycosylation-related genes associated with idiopathic pulmonary fibrosis. Ten candidate diagnostic genes were identified, and a diagnostic model with strong performance in lung tissue datasets was developed and validated.

Abstract

Idiopathic pulmonary fibrosis (IPF) is a progressive chronic lung disease with an unclear etiology, and the contribution of glycosylation-related genes (GRGs) to its pathogenesis remains poorly understood. This study focuses on elucidating the potential mechanisms of GRGs in IPF, identifying key biomarkers, and developing a diagnostic model using bioinformatics and machine learning. Transcriptomic data from the GEO database were integrated with GRGs to investigate their role in IPF. Differentially expressed genes (DEGs) were identified and subjected to functional enrichment and protein-protein interaction (PPI) analyses. A machine learning workflow combining LASSO regression, support vector machine-recursive feature elimination (SVM-RFE), and XGBoost was applied to identify key genes and construct a diagnostic model using the GSE150910 training dataset. Model performance was subsequently evaluated in four independent validation datasets. Additional analyses, including gene set enrichment analysis (GSEA), immune infiltration analysis, drug-gene interaction analysis, molecular docking, and RT-qPCR validation using peripheral blood samples from IPF patients and healthy controls, were performed to investigate the biological relevance of the identified genes. A total of 126 glycosylation-related DEGs were identified, and 10 key genes were selected. The diagnostic model achieved area under the curve (AUC) values of 0.963 and 0.814 on lung tissue datasets, and 0.681 and 0.706 on blood datasets. Immune infiltration analysis revealed differences in B-cell abundance between patient subgroups, and RT-qPCR validation confirmed the differential expression of selected genes in clinical samples. These findings provide insight into the potential involvement of GRGs in IPF and support their relevance as candidate diagnostic genes.

Introduction

Idiopathic pulmonary fibrosis (IPF) is a progressive, fibrosing interstitial lung disease characterized by a heterogeneous pathogenesis, which involves genetic predisposition, environmental triggers, and immune dysfunction^1,2. The disease presents significant clinical challenges, primarily due to difficulties in early diagnosis, limited therapeutic efficacy, and a generally poor prognosis. The prevalence of IPF ranges from 20 to 80 patients per 100,000, and without treatment, the median survival time after diagnosis is typically 3–5 years^3,4. Identifying patients at higher risk of rapid disease progression is critical for effective treatment strategies, enabling targeted therapies for those most likely to experience disease decline.

Glycosylation, a co-translational and post-translational modification of proteins, plays a critical role in regulating protein activity and cell biology^5,6. Dysregulated protein glycosylation likely contributes to IPF progression^7,8. In many rheumatic diseases, alterations in protein glycosylation are observed, which are closely associated with inflammatory processes and disease progression⁹. For example, Lyu et al. found that CRLF1 and GLG1 are key players in the pathological processes of osteoarthritis¹⁰. Moreover, glycosylation has been implicated in immune cell activation within pulmonary tissues, further linking it to the pathogenesis of pulmonary diseases¹¹. Given this, dysregulation of protein glycosylation may also contribute to the progression of IPF. However, the role of glycosylation-related genes (GRGs) in IPF remains underexplored, and their potential for diagnosing IPF warrants further investigation.

Herein, we conducted a study integrating machine learning approaches and RNA sequencing data to identify key GRGs associated with IPF. Previous studies in IPF have often relied on conventional differential expression analyses or a single machine learning algorithm, which may limit the robustness and generalizability of candidate biomarkers^12,13. To address these limitations, we employed a multi-step computational workflow that combined differential expression analysis with multiple machine learning algorithms, including LASSO, SVM-RFE, and XGBoost, to improve feature selection reliability. Furthermore, the diagnostic model was evaluated across multiple independent datasets, and key findings were further validated using clinical samples. A diagnostic model was subsequently developed and validated, and immune infiltration characteristics were analyzed. In addition, the immune microenvironment and potential molecular interactions associated with the identified genes were explored. Although this study focused on IPF, the analytical strategy employed here may be applicable to other chronic pulmonary diseases and transcriptomic datasets, providing a reproducible framework for key gene identification and validation. See Figure 1 for the specific process.

Protocol

This study was approved by the Institutional Review Board of The First Affiliated Hospital of Bengbu Medical University (Approval Number: 2023YJS162). Written informed consent was obtained from all participants prior to sample collection.

Data collection

The transcriptome data used in this study were obtained from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/; RRID: SCR_005012). The primary training dataset, GSE150910, was generated using the Illumina NovaSeq 6000 platform (GPL24676) and comprised 103 IPF and 103 normal lung tissue samples. To validate the findings, independent datasets GSE24206, GSE110147, GSE93606, and GSE38958 were employed. Detailed information on each dataset is provided in Table 1. Additionally, a total of 636 GRGs were selected from a previously published study¹⁴.

Differential expression analysis and functional characterization of glycosylation-related DEGs

Differential expression analysis between IPF and normal lung tissue samples from the GSE150910 dataset was conducted using the R package DESeq2 (RRID: SCR_015687). Genes with an adjusted P-value (padj) < 0.05 and |log2FoldChange| > 0.5 were considered differentially expressed genes (DEGs). Glycosylation-related DEGs (GR-DEGs) were identified by intersecting the DEGs with a predefined set of 636 GRGs. To further investigate the biological roles of these genes, Gene Ontology (GO) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were performed to elucidate their functional roles and pathway involvement. A protein-protein interaction (PPI) network was generated using the STRING database (RRID: SCR_005223)¹⁵ with an interaction confidence score threshold of > 0.7, to elucidate the molecular interactions and potential regulatory mechanisms of the GR-DEGs in IPF.

Screening of key genes and construction of a diagnostic model 

To screen key genes for IPF from GR-DEGs, we employed multiple machine learning algorithms. Initially, LASSO regression (RRID: SCR_003418) was fitted with binary logistic regression (family = “binomial”) and the optimal penalty parameter λ was selected via 10-fold cross-validation (nfold = 10). The final selection results were the features with non-zero coefficients corresponding to λ_(min) (0.01700442). For SVM-RFE, the rfe function from the R package caret was used. Recursive feature elimination was performed via 10-fold cross-validation (method = “cv”, number = 10), progressively filtering features from 1 to 126, with accuracy used to determine the optimal feature subset. In XGBoost, the objective function was set to binary logistic regression (objective = “binary: logistic”), with the evaluation metric set to log loss (eval_metric = “logloss”), the number of iterations (nrounds) set to 100, and the learning rate (eta) set to 0.1. The top 20 genes were selected based on their feature importance scores (Gain). By intersecting the results from these methods, a refined set of key genes was identified. Based on this gene set, an XGBoost diagnostic model was constructed using the training dataset (GSE150910), and its predictive performance was evaluated using receiver operating characteristic (ROC) analysis on external validation datasets (GSE110147, GSE24206, GSE93606, and GSE38958). Additionally, a nomogram was developed to visualize the contribution of each selected gene to disease probability, and the model’s clinical utility was further assessed through calibration curves and decision curve analysis (DCA).

Exploration of the biological pathways of key genes

To explore the biological context of the key genes identified by machine learning. Gene Set Enrichment Analysis (GSEA)¹⁶ was performed based on the gene lists from Molecular Signatures Database (MSigDB) (RRID: SCR_016863)¹⁷and the pathways were screened for NES > 1. The top enriched pathways were visualized using the enrichplot function.

Exploring biological function and immune landscape differences in IPF subtypes based on key gene scores

Based on the expression profiles of the identified key genes, single-sample Gene Set Enrichment Analysis (ssGSEA) scores were calculated and used to stratify IPF patients into high- and low-scoring groups based on the median score. Differential expression analysis was performed between the two groups, followed by GSEA (RRID: SCR_003199)¹⁸ to conduct GO Biological Processes (GOBP) and KEGG pathway enrichment analyses on the DEGs.

Analysis of immune cell infiltration and key gene expression differences

Following subgroup stratification based on ssGSEA scores, immune infiltration differences between the high- and low-score groups were evaluated. First, we calculated the relative abundances of 22 immune cell types in the samples using the CIBERSORT algorithm (RRID: SCR_016955)¹⁹ in combination with the LM22 feature matrix. Specifically, the deconv_tme function in the R package IOBR (parameters: method = “cibersort”, arrays = FALSE, perm = 200) was used for the calculations, and box plots were generated using the ggpubr package (RRID: SCR_021139) to assess differences in immune cell infiltration between the high-score and low-score groups. Additionally, the gsva function in the R package GSVA (using the ssGSEA method) was used to calculate enrichment scores for 28 immune cell types. These scores were then normalized using Min-Max scaling to map them to the [0, 1] interval, facilitating comparisons across cell types. Finally, Wilcoxon rank-sum tests were performed to compare the expression of key genes between normal samples and IPF patients in the GSE150910 and GSE110147 datasets, providing a comprehensive analysis of immune cell infiltration and gene expression differences across IPF subgroups.

Validation of key genes in IPF patients using RT-qPCR analysis

To validate the diagnostic relevance of the identified genes, six genes with the highest importance scores from the XGBoost algorithm were selected for expression level validation in IPF patients and healthy controls using reverse transcription quantitative PCR (RT-qPCR). A total of 20 blood samples, including 9 from IPF patients and 11 from healthy individuals, were collected from the First Affiliated Hospital of Bengbu Medical University. Total RNA was extracted from the blood samples, and RNA concentration was measured using a multifunctional microplate reader. RNA quality was assessed prior to downstream analyses. Genomic DNA was removed during reverse transcription, and primer sequences used for RT-qPCR are listed in Table 2. Primer specificity was verified by melt curve analysis. GAPDH was used as the internal reference gene. Relative gene expression levels were calculated using the 2^-ΔΔCt method. This validation step provides preliminary experimental support for the differential expression and potential diagnostic relevance of the identified genes in IPF.

Statistical analysis

Data were analyzed in R, and the Wilcoxon test was used to detect differences between the two groups. The GSEA, GO, and KEGG enrichment analysis were performed using the R package clusterProfiler (RRID: SCR_016884). A p-value < 0.05 was considered significant unless otherwise specified.

Results

Functional analysis and interaction network analysis of GR-DEGs  

To investigate gene involvement in IPF, differential expression analysis of GSE150910 was conducted, identifying 4,216 DEGs (Figure 2A). By intersecting these DEGs with 636 GRGs, 126 disease-associated glycosylation genes were identified (Figure 2B). Furthermore, to evaluate the potential biological functions of intersection genes, we c...

Discussion

IPF is a chronic, progressive disease of unknown etiology that significantly reduces patients' quality of life and leads to earlier mortality⁴. Despite advances in therapeutic strategies, the median survival time remains low post-diagnosis, and the pathogenesis of IPF is still not fully understood²². Therefore, elucidating the molecular mechanisms underlying IPF is both urgent and critical. Previous studies have established a link between deregulated glycosylation and p...

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
Caret	R Foundation for Statistical Computing	Version 6.0.94	SVM-RFE
ClusterProfiler	R Foundation for Statistical Computing	Version 4.12.6	GSEA, GO and KEGG enrichment analysis
DESeq2	R Foundation for Statistical Computing	Version 1.44.0	Differential expression analysis
Fast First-Strand cDNA Synthesis Mix	Goonie Biotech Co., Ltd.	500-101	cDNA synthesis for RT-qPCR
FastTaq qPCR SYBR Green Master Mix	Gooniebio	500-102	Quantitative PCR amplification
Fluorescent Quantitative PCR System	Bio-Rad Laboratories, Inc.	CFX96	Quantitative PCR amplification
GEO Database	NCBI	https://www.ncbi.nlm.nih.gov/geo/	Transcriptome data source
ggpubr	R Foundation for Statistical Computing	Version 0.6.0	Box plots generation
GSVA	R Foundation for Statistical Computing	Version 1.52.3	ssGSEA
MSigDB	Broad Institute	https://www.gsea-msigdb.org/	GSEA analysis
PCR System	Bio-Rad Laboratories, Inc.	PTC Tempo	cDNA synthesis for RT-qPCR
Primers	Tsingke Biotechnology Co., Ltd.	N/A	Custom synthesis of qPCR primers
R software	R Foundation for Statistical Computing	Version 4.4.1	Statistical analysis and visualization
RNA Isolater Total RNA Extraction Reagent	Vazyme Biotech Co., Ltd.	R401-01	Total RNA extraction from whole blood
STRING Database	STRING Consortium	https://string-db.org/	PPI network construction
Synergy HTX Multi-Mode Reader	BioTek Instruments, Inc.	Synergy HTX	RNA concentration measurement