UPP1 as a Diagnostic Biomarker: Insights from Integrative Bioinformatics and Immune Infiltration Analyses in COPD

Summary

This study presents a reproducible bioinformatics pipeline integrating transcriptomics, machine learning, and immune infiltration analysis to identify potential diagnostic biomarkers and regulatory networks in chronic obstructive pulmonary disease (COPD).

Abstract

COPD is a progressive respiratory disorder characterized by persistent airflow limitation and chronic inflammation, yet the role of N4-acetylcytidine (ac4C) RNA modification in its pathogenesis remains largely unexplored. This study aimed to systematically screen for ac4C-related genes (ac4C-RGs) from a published database and investigate their regulatory networks in COPD, thereby identifying potential biomarkers for further mechanistic studies without assuming a direct regulatory relationship between any specific gene and ac4C modification. Differentially expressed genes (DEGs) were identified from transcriptomic profiles, and weighted gene co-expression network analysis (WGCNA) was applied to uncover key co-expression modules. Cross-analysis among DEGs, significant modules, and ac4C-RGs was conducted. Key genes were screened using LASSO regression, XGBoost, and random forest algorithms, followed by logistic regression‑based diagnostic model construction. Model performance was evaluated by receiver operating characteristic (ROC) curve analysis, area under the curve (AUC) with 95% confidence intervals, calibration curve assessment, and decision curve analysis (DCA). A total of 160 overlapping genes were identified, and six hub genes (PTRF, PRKCDBP, UPP1, TOR3A, FAM168B, and B4GALT2) were consistently selected by all three machine learning algorithms. The diagnostic model demonstrated good discriminative performance, with AUCs of 0.766, 0.759, and 0.723 in the training, internal test, and external validation sets, respectively. Regulatory network analysis suggested potential ceRNA axes and transcription factor interactions, while immune infiltration profiling revealed significant correlations between key genes and multiple immune cell subsets. Drug-gene interaction analysis and molecular docking indicated that fluorouracil, capecitabine, and 5-benzylacyclouridine may exhibit favorable predicted binding affinities with UPP1. In conclusion, PTRF, PRKCDBP, UPP1, TOR3A, FAM168B, and B4GALT2 were identified as potential ac4C-related biomarkers in COPD, potentially involved in immune and metabolic regulation, providing a foundation for future functional investigations and therapeutic exploration.

Introduction

COPD is a chronic and heterogeneous respiratory disorder characterized by progressive airflow limitation resulting from abnormalities in alveolar and airway structures^1,2. A protease-antiprotease imbalance, oxidative stress, chronic inflammation, and cellular senescence constitute the core pathophysiological mechanisms of COPD, leading to structural destruction and functional impairment of lung tissue^3,4. Moreover, COPD is influenced by multiple risk factors, including long-term smoking, environmental pollution, occupational exposure, respiratory infections, and genetic susceptibility^5,6. COPD has imposed a significant burden on the global economy. It is projected to account for 0.111% of global GDP annually from 2020 to 2050⁷. Although current therapeutic strategies, such as bronchodilators, inhaled corticosteroids, pulmonary rehabilitation, and long-term oxygen therapy, can alleviate symptoms, halting disease progression remains challenging. Because of marked clinical heterogeneity, patient outcomes vary widely⁸. Therefore, novel diagnostic biomarkers and prognostic indicators are urgently needed to enhance COPD management and improve patient survival.

RNA modification refers to the chemical alteration of RNA molecules, which can change RNA structure and function to regulate gene expression^9,10. Common RNA modifications include N6-methyladenosine (m6A), pseudouridine (Ψ), 5-methylcytosine (m5C), and ac4C¹¹. The ac4C modification plays a crucial role in preserving mRNA stability and promoting mRNA translation^12,13. NAT10 is the only known eukaryotic enzyme that catalyzes ac4C modification, and its activity is essential for the formation of this modification¹⁴. Studies have shown a strong correlation between oxidative stress, cellular senescence, inflammation, and ac4C modification. For example, CD4+ T lymphocytes in colon tissues from individuals with inflammatory bowel disease (IBD) exhibit markedly elevated levels of NAT10¹⁵. NAT10 enhances ac4C acetylation of the chemokines CCL2 and CXCL1, thereby promoting the infiltration of macrophages and neutrophils and exacerbating inflammatory damage¹⁶. The cellular response to oxidative stress may involve ac4C modification, as evidenced by the marked increase in ac4C levels under oxidative stress. In addition, NAT10 promotes PM_2.5-induced pulmonary fibrosis by stabilizing TGFB1 mRNA via ac4C modification, thereby triggering epithelial-to-mesenchymal transition¹⁷. However, the role of ac4C modification in COPD remains largely unexplored, highlighting the need for further research in this area.

The GEO database was utilized in this study to identify DEGs between COPD patients and controls. A published list of ac4C-RGs, compiled from multiomics data¹⁸, was then integrated with the DEGs to identify overlapping candidates. Given that ac4C modification is known to influence inflammation and oxidative stress–key processes in COPD–we hypothesized that genes associated with the ac4C regulatory network might be dysregulated in COPD. However, the overlapping genes were not considered direct substrates of NAT10 or genes directly regulated by ac4C; rather, they were considered candidates associated with the ac4C-related network. Key genes were screened using multiple machine learning algorithms, followed by construction and validation of a diagnostic model. Subsequently, relevant pathways involved in COPD pathogenesis were determined, and potential targeted drugs were predicted. Finally, RT-qPCR assays were performed to confirm the expression levels of key genes, providing preliminary insights into COPD pathophysiology and potential avenues for therapeutic exploration.

Protocol

Institutional review board statement

This study was conducted in accordance with the Declaration of Helsinki. The protocol was approved by the Ethics Committee of Shenzhen Luohu Hospital of Traditional Chinese Medicine (approval no. 2024-LHQZYYYXLL-KY-039), and written informed consent was obtained from all participants prior to enrollment. Details of the research tools and materials used in this protocol are provided in the Table of Materials.

Data source and processing

COPD-related gene expression datasets were obtained from the Gene Expression Omnibus (GEO). The GSE54837 dataset was used as the transcriptome dataset, and the GSE112811 dataset served as the validation set (Table 1). The ac4C-RGs were collected from the literature¹⁸. DEGs between the COPD and control groups were identified using the R package limma. DEGs were considered statistically significant if |log₂FC| > 0 and p < 0.05. Volcano plots were generated to visualize the overall distribution of gene expression changes.

Construction of WGCNA

WGCNA was performed on the GSE54837 dataset using R to identify COPD-related modules. Prior to network construction, outlier samples were identified and removed through hierarchical clustering analysis using the hclust function with the average linkage method and a Euclidean distance metric. The optimal soft-thresholding power (β = 10) was selected to achieve a scale-free topology fit index R² ≥ 0.85, balancing scale-free topology and mean connectivity. An adjacency matrix was constructed and transformed into a topological overlap matrix (TOM). Gene modules were identified using the dynamic tree-cutting algorithm (deepSplit = 2, minClusterSize = 50). Modules with eigengene correlations > 0.75 were subsequently merged using the mergeCloseModules function. Module eigengenes were then correlated with clinical traits (COPD status, age, sex, and smoking status) using Pearson correlation coefficients to identify COPD-associated modules for subsequent analysis.

Screening, enrichment analysis, and PPI network analysis of overlapping genes

A Venn diagram was generated using the R package ggvenn to identify genes overlapping among the DEGs, MEsalmon module genes, and ac4C-RGs. Functional enrichment analysis of the overlapping genes was performed using Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases with the R package clusterProfiler. Protein-protein interaction (PPI) information was obtained from the STRING database (https://string-db.org/) to analyze protein-level interactions among the overlapping genes. Cytoscape software was used to visualize the resulting PPI network.

Identifying key genes through machine learning

Three machine learning techniques were applied: least absolute shrinkage and selection operator (LASSO) regression, extreme gradient boosting (XGBoost), and random forest (RF). LASSO regression was implemented using the glmnet package with 10-fold cross-validation to determine the optimal penalty parameter λ. The type.measure parameter was set to "deviance", and the family parameter was set to "binomial". The optimal λ was selected using the λ_min criterion, which minimizes cross-validated deviance, yielding 17 genes. XGBoost was performed using the xgboost package with the following hyperparameters: nrounds = 100, max_depth = 6, eta = 0.3, subsample = 0.8, colsample_bytree = 0.8, and eval_metric = "logloss". Feature importance was ranked by the gain metric, and the top 30 genes were selected. Random forest was implemented using the randomForest package with ntree = 200. Feature importance was ranked by mean decrease in Gini, and the top 30 genes were selected. The genes selected by the three machine learning methods were intersected to identify key genes for subsequent analyses.

Building and assessing the logistic regression model for risk prediction

The GSE54837 dataset was randomly divided into a training set (70%) and a testing set (30%). A logistic regression model was constructed on the training set using the glm function from the MASS package, with the expression levels of key genes as input features. Model performance was evaluated using ROC curves generated with the pROC package. The 95% confidence intervals for the AUC were computed via 2,000 bootstrap replicates. Model calibration was assessed using calibration curves generated with 1,000 bootstrap resamples (rms package). DCA was performed using the dca package to evaluate net clinical benefit across a range of threshold probabilities. A nomogram was constructed using the nomogram function from the rms package to facilitate individualized risk estimation.

The regression equation was:

logit(P) = 0.5823 + 0.6010 × UPP1 - 0.6563 × PTRF + 0.3853 × B4GALT2 - 0.3972 × FAM168B + 0.1848 × PRKCDBP - 0.4787 × TOR3A.   (1)

Here, P represents the predicted probability of COPD, and each coefficient represents the contribution of the corresponding gene-expression value to the log odds of COPD.

Expression analysis, GeneMANIA network, and molecular regulatory network

Gene expression levels between the COPD and control groups in the GSE54837 dataset were compared using the Wilcoxon rank-sum test. Box plots were generated using the ggplot2 package to visualize the distribution of expression levels, with median, interquartile range (IQR), and individual data points overlaid. GeneMANIA was used to construct gene networks and predict functional interactions. The search was performed with default parameters: species = Homo sapiens, maximum related genes = 20. The resulting network was downloaded and visualized, with edge colors indicating interaction types. A competitive endogenous RNA (ceRNA) network was constructed to investigate post-transcriptional regulatory mechanisms. miRNAs targeting the six key genes were predicted using two independent databases: DIANA-microT (score ≥ 0.8) and miRanda (score ≥ 140, energy ≤ −20 kcal/mol). The intersection of miRNAs identified by both databases was used to construct miRNA-mRNA pairs. Subsequently, lncRNAs targeting these miRNAs were predicted using the StarBase database. An lncRNA-miRNA-mRNA regulatory network was constructed and visualized using Cytoscape. Transcriptional regulatory relationships were predicted using ChIP-X Enrichment Analysis Version 3 (ChEA3). For each key gene with predicted TFs, the top 10 transcription factors with the highest enrichment scores were selected. A TF-target regulatory network was constructed in Cytoscape.

Gene set enrichment analysis and immune cell infiltration assessment

Gene set enrichment analysis (GSEA) was performed using the clusterProfiler package to investigate the biological functions of each key gene. For each key gene, samples were divided into high and low expression groups based on the median value. Differential expression analysis between the two groups was performed using limma, and the resulting gene list was ranked by the signed log₂ fold-change. GSEA was conducted using the gseGO function for GO biological process terms and the gseKEGG function for KEGG pathways, with the following parameters: minGSSize = 10, maxGSSize = 500, pvalueCutoff = 0.05, and nPerm = 1,000. The relative abundance of 28 immune cell types was estimated using single-sample gene set enrichment analysis (ssGSEA) implemented in the GSVA package. A curated gene set signature matrix comprising marker genes for 28 immune cell types was obtained from previous literature¹⁹. For each sample, the gsva function was applied with method = "ssgsea", ssgsea.norm = TRUE, and kcdf = "Gaussian". Spearman correlation coefficients between the ssGSEA enrichment scores and the expression levels of the six key genes were calculated using the cor.test function. The p values were adjusted for multiple testing using the Benjamini-Hochberg method. The correlation matrix was visualized as a heatmap using the pheatmap package.

Drug prediction, molecular docking, and disease association analysis

Potential therapeutic compounds targeting key genes were identified using the DrugBank database. A “key gene-targeting drug” interaction network was constructed in Cytoscape to visualize predicted drug-gene interactions. Molecular docking was performed using the CB-Dock2 platform to assess binding affinities. The 3D protein structure of human UPP1 was retrieved from the Protein Data Bank (PDB ID: 7B8T). Drug molecular structures (SMILES format) were obtained from PubChem. Docking was performed using the AutoDock Vina engine, and outputs were ranked by binding free energy (ΔG, in kcal/mol). Docking complexes were visualized using PyMOL. Associations between key genes and human diseases related to environmental exposures were investigated using the Comparative Toxicogenomics Database (CTD). Each gene was individually queried, and the top ten most strongly associated diseases were extracted and visualized using radar plots.

RT-qPCR protocol

Peripheral venous blood samples were collected from eight COPD patients and eight healthy controls at Shenzhen Luohu Hospital of Traditional Chinese Medicine. COPD was diagnosed according to the Global Initiative for Chronic Obstructive Lung Disease (GOLD) criteria, defined as a post-bronchodilator FEV1/FVC < 0.70. The control group comprised age- and sex-matched healthy volunteers with no history of respiratory diseases and normal pulmonary function tests (FEV1% predicted ≥ 80% and FEV1/FVC ≥ 0.70). The baseline information of the patients is shown in Table 2. Total RNA was extracted from COPD blood samples using a blood RNA extraction kit. For cDNA synthesis, 500 ng of total RNA was reverse-transcribed using a cDNA synthesis kit with genomic DNA removal following the provided protocol. The resulting cDNA was diluted to 150 ng/μL.

RT-qPCR was performed using a SYBR Green-based qPCR master mix on a real-time PCR system. Each 10 μL reaction contained 5 μL of 2x SYBR Green master mix, 0.5 μL each of forward and reverse primers (10 μM), 1 μL of diluted cDNA (15 ng/μL), and 3 μL of nuclease-free water. The cycling conditions were initial denaturation at 95 °C for 5 min, followed by 40 cycles of 95 °C for 10 s and 60 °C for 30 s, with a final melting curve analysis from 60 °C to 95 °C to verify amplification specificity. All reactions were performed in technical triplicates. β-actin was used as the internal reference gene. Primer efficiency for each target gene was validated using standard curve dilution series and ranged from 90% to 110%. Gene expression levels were normalized to β-actin, and relative expression was calculated using the 2^-ΔΔCt method. Statistical comparisons between COPD and control groups were performed using the Mann-Whitney U test.

Statistical analysis

Network visualizations were created using Cytoscape, and statistical analyses were performed using R software. Unless otherwise indicated, the Mann-Whitney U test was used for nonnormally distributed data, and Student's t-test was used for normally distributed data to compare two groups. A value of p < 0.05 was considered statistically significant.

Results

Identification of overlapping genes, enrichment analysis, and PPI network construction

From the GSE54837 dataset, 3,371 DEGs were identified, including 1,675 upregulated and 1,696 downregulated DEGs. The top 10 genes with the most significant up- and downregulation are shown in Figure 1A. Hierarchical clustering of the GSE54837 data was performed (Supplementary Figure 1A), and a soft-thresholding power of 10 was applied to ensure ...

Discussion

The conserved RNA modification ac4C, which predominantly occurs in messenger RNA (mRNA) and transfer RNA (tRNA), enhances mRNA stability and translation efficiency²⁰. NAT10 is the only known RNA acetyltransferase mediating ac4C modification²¹. Studies have shown that NAT10 is upregulated in lung epithelial cells of COPD patients. Knockdown of NAT10 disrupts mitochondrial function and transcriptomic responses²². Based on integrated multiomics analysis...

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
β-actin primers	Tsingke	N/A	Forward: 5’-CATGTACGTTGCTATCCAGGC-3’ Reverse: 5’-CTCCTTAATGTCACGCACGAT-3’
B4GALT2 primers	Tsingke	N/A	Forward: 5’-GGGCAGACTGCTGATCGAG-3’ Reverse: 5’-CCGGTGTCTAAAGGGGATGAT-3’
CB-Dock2	LabShare	Online	Molecular docking
clusterProfiler	Bioconductor	v4.14.6	Enrichment analysis
Cytoscape	Cytoscape Consortium	v3.8.3	Network visualization
DrugBank	University of Alberta	Online	Drug prediction
FAM168B primers	Tsingke	N/A	Forward: 5’-TCTGGGGTTCCCTATGCAAAT-3’ Reverse: 5’-GTAGGATTCGCTCCAGGATACA-3’
glmnet	CRAN	v4.1	LASSO regression
GSVA	Bioconductor	v1.52.3	ssGSEA analysis
Hifair III 1st Strand cDNA Synthesis Supermix	YEASEN	11141ES	cDNA synthesis
Hieff RTPCR SYBR Green Master Mix	YEASEN	11201ES	qPCR amplification
limma	Bioconductor	v3.54.0	Differential expression
LightCycler 480 II System	Roche	LightCycler 480 II	Real-time PCR
PTRF primers	Tsingke	N/A	Forward: 5’-GGGCCGTAGACCAGATCCA-3’ Reverse: 5’-CTTGCTCACCGTATTGCTCGT-3’
PRKCDBP primers	Tsingke	N/A	Forward: 5’-CACGTTCTGCTCTTCAAGGAG-3’ Reverse: 5’-TGTACCTTCTGCAATCCGGTG-3’
R software	R Foundation	v4.4.2	Statistical computing
randomForest	CRAN	v4.7	Random Forest
RNA isolater MolPure Blood RNA Kit	YEASEN	19241ES50	RNA extraction
STRING database	EMBL	Online	PPI network
TOR3A primers	Tsingke	N/A	Forward: 5’-CCCTTGCTCTGTCGTTCCAC-3’ Reverse: 5’-CCCGTCCCGATACAGGTTC-3’
UPP1 primers	Tsingke	N/A	Forward: 5’-CTGTCAGTCATGGTATGGGCA-3’ Reverse: 5’-GAGCACCGGGCATAGTACA-3’
WGCNA	CRAN	v1.72	Co-expression network
xgboost	CRAN	v1.7	XGBoost algorithm