We observed the hospitalization of 56 patients with COPD to observe the real-life disease burden of COPD on patients and to demonstrate the importance of drug development. The overall scheme of COPD drug target screening was shown in Figure 2 . Briefly, we conducted a Summary-data-based Mendelian Randomization (SMR) analysis, using genetic variants that affect the expression of druggable genes in blood and lung tissue as instrumental variables to estimate the causal relationships between these genetic variations and doctor-diagnosed COPD, FEV1/FVC<0.7, and FEV1. The SMR method utilizes summary-level data from GWAS and QTL studies to test for associations between gene expression levels and complex traits of interest. Two-sample MR estimates the causal effects by separately estimating the associations between instrumental variables and both exposure factors and outcomes, requiring two sets of independent sample data. SMR can be used to prioritize potential genes for further research, while two-sample MR is suitable for causal inference. For SMR results meeting the significance threshold and passing the heterogeneity independent instrument (HEIDI) test, we conducted colocalization analyses to determine if the same SNPs were driving both the outcome and exposure. Furthermore, eQTL two-sample MR and protein SMR analyses were employed to ensure consistency in the results. Finally, we matched existing drugs to the identified druggable genes and conducted a PheWAS analysis to preliminarily evaluate the safety of these genes as potential drug targets (i.e., whether they may cause other diseases).

This study used the big data platform of the Clinical Medical Research Centre of the First Affiliated Hospital of Guangzhou Medical University, and the patient’s clinical data were obliterated for privacy-related information before they were provided to the investigators. The observational study was approved by the Ethics Committee of the First Affiliated Hospital of Guangzhou Medical University (ES-2023-017-01).

Druggable genes were obtained from a review by Finan et al. (Finan et al., 2017). A total of 4317 druggable genes on autosomes were utilized in this study.

The eQTL data for whole blood were obtained from the eQTLGen consortia (www.eQTLgen.org), while the lung-specific eQTL data were obtained from the GTEx v.8 consortia (www.gtexportal.org). Detailed descriptions of the data can be found in the original publications (Aguet et al., 2020; Võsa et al., 2021). The eQTLGen consisted of a total of 25,482 whole blood samples from 37 datasets, with the majority of samples being of European ancestry. Gene expression levels of the samples were analyzed by Illumina, Affymetrix expression arrays, and RNA-seq. The GTEx version 8 data contains lung tissue samples from 515 postmortem donors, the majority of whom were of European ethnicity. The mRNA from each tissue sample was sequenced at a median depth of 82.6 million reads. We Used eQTLGen and GTEx v.8 to identify all conditionally independent eQTLs in whole blood and lung tissue, which covered 2531 of the 4317 druggable genes in at least one tissue. The pQTL data for whole blood were obtained from Ferkingstad et al. study, which contained pQTLs for 4907 plasma proteins (Ferkingstad et al., 2021).

For the SMR analysis, we downloaded the cis-eQTL results from eQTLGen and GTEx v.8 consortia in SMR format in May 2023 (yanglab.westlake.edu.cn/software/smr). In the latest version, every SNP-gene combination with a distance <1Mb was included. To further validate the effect of druggable genes at the protein level, we extracted available cis-pQTLs from the study by Ferkingstad et al., considering those with a significance threshold of P < 5e−8.

For the two-sample MR analysis, the cis-eQTLs, which were defined as significant (5e-8) SNPs located within 5 kb upstream of the starting point or 5 kb downstream of the endpoint of a druggable gene, were extracted from the eQTLGen and GTEx v.8 consortia.

FEV1 quantifies the maximal volume of air expelled in the first second of a forceful exhalation, serving as a critical parameter for evaluating pulmonary function. A reduced FEV1/FVC ratio (<0.7) typically indicates airflow obstruction (spirometry defined COPD). The doctor diagnosed COPD in this study comes from the J-10 disease codes record of the healthcare system. This diagnosis incorporates an assessment based on the patient’s exposure history (such as smoking), symptoms (coughing, expectoration, or dyspnea), and post-bronchodilator pulmonary function tests. Conversely, spirometry-defined COPD (with a pre-bronchodilator FEV1/FVC<0.7) is identified during population health screenings. Participants meeting these conditions may not have any exposure history or symptoms (accounting for more than 70% of cases) (Fang et al., 2018) and did not undergo post-bronchodilator pulmonary function tests, making them a distinctly different cohort from those diagnosed with COPD by a doctor. Consequently, this study separately investigates the druggable genes associated with both doctor-diagnosed COPD and spirometry-defined COPD (pre-bronchodilator FEV1/FVC<0.7). All doctor-diagnosed COPD data were obtained from the FinnGen consortium (version R8, downloaded in May 2023), which comprises 16,410 patients with COPD and 283,589 controls (differentiated according to J-10 disease codes in the medical electronic system). The cohorts for FEV1 and FEV1/FVC<0.7 were derived from the UK Biobank and SpiroMeta consortium. For further details regarding the data, please refer to the original literature and the official website ( Supplementary Table 1 ) (Shrine et al., 2019; Higbee et al., 2021).

The MR analysis is based on three fundamental assumptions: 1) genetic variation must be associated with the exposure; 2) genetic variation must not be associated with any confounding factors; and 3) genetic variation must not have a direct relationship with the outcome (Sanderson et al., 2022). As an extension of the MR concept, SMR employs two-step least squares (2SLS) to estimate the effect of the exposure on the outcome (e.g., doctor diagnosed COPD) using genetic variation that is significantly associated with the exposure (e.g., genes targetable for drug intervention at the transcriptional or protein expression level) (Sanderson et al., 2022). The causal associations were computed as follows:

The odds ratio (OR) estimates of the effect of druggable gene expression on COPD risk are calculated using the formula: OR_{druggable gene-COPD}=exp(β_{druggable gene-COPD}), where OR represents the odds ratio estimate per 1-ln increase in the level of the druggable gene, and exp is the base of the natural logarithm.

SMR analyses were completed using the SMR software for Linux version 1.0.3 with default options (https://yanglab.westlake.edu.cn/software/smr). To account for multiple testing and control for genome-wide type I error rates, we used the Bonferroni correction to adjust the P values (Armstrong, 2014).

We conducted the HEIDI test to assess whether the observed associations were influenced by vertical pleiotropy rather than linkage disequilibrium (LD) of the dependent variable (Zhu et al., 2016). LD estimation was performed using the European pedigree genomes from the 1000 Genomes Project Consortium as a reference (Zheng-Bradley and Flicek, 2017). Associations with a significance level of P < 0.01 in the HEIDI test were excluded as they may be attributed to LD rather than pleiotropy.

An alternative Bayesian test was performed to examine the colocalization of the two traits and estimate the posterior probability of shared variation (Giambartolomei et al., 2014). For each leading SNP in the outcome GWAS database, we retrieved all SNPs within a 100 kb range upstream and downstream. The posterior probability of H4 (PP.H4) was estimated, with the commonly used threshold of PP.H4 > 0.8 for GWAS and QTL associations (Navarro et al., 2021).

To further validate the results of the SMR analysis, the two-sample MR approach was employed to estimate the effect on outcomes by using eQTL of druggable genes as instruments (Li et al., 2023). We assessed the strength of the correlation between instrumental variables and exposure using the F-statistic (>10) (Burgess and Thompson, 2011). In the presence of directional pleiotropy, we utilized the MR-egger estimated effect size as the outcome (Burgess and Thompson, 2017). The study proceeded with the two-sample MR analysis using MR-Egger, weighted median, inverse variance weighting, simple mode, and weighted mode methods. When multiple analytical methods yield similar results, we consider the MR results to be more robust. Cochran’s Q test was used to assess the overall heterogeneity of Wald ratios (Higgins et al., 2003). In addition, the Steiger test was conducted to determine whether the exposure precedes the outcome (https://jean997.github.io/cause) (Hemani et al., 2017).

To further assess the potential horizontal pleiotropy and any possible side effects of potential drug targets, we conducted a PheWAS analysis using the AstraZeneca portal (https://azphewas.com/). The original study utilized data from approximately 15,500 binary phenotypes and around 1,500 continuous phenotypes from about 450,000 participants who underwent exome sequencing, as published by the UK Biobank. The complete methodology can be found in the original article (Wang et al., 2021). We performed multiple corrections and set the threshold at 2E-9 to account for possible false positives.

All DNA positions are based on the human reference genome build hg19 (GRCh37). Data processing was carried out using R software version 4.2 with “TwosampeMR” and “easyMR” packages.

mNGS and hybridization capture-based targeted mNGS were used to detect pathogens. Briefly, about 1 mL of sample was centrifuged at 12,000 g for 5 min to collect the pathogens and human cells. Next, 50 μL of precipitate underwent depletion of host nucleic acid using 1 U of Benzonase (Sigma) and 0.5% Tween 20 (Sigma) and incubated at 37°C for 5 min. Subsequently, the DNA nucleic acid was extracted using a QIAamp UCP pathogen minikit. Qubit 4.0 (Thermo Fisher Scientific, MA, USA) was used to measure extracted DNA concentrations, which were used to construct metagenomics libraries. The inspected and qualified library was sequenced on the Nextseq 550 platform (Illumina, San Diego, USA). For targeted NGS, the constructed library from each sample was used for hybrid capture-based enrichment of microbial probes. To remove adapter and low-quality, low complexity, and short reads of < 35 bp, raw data generated by the sequencing were filtered. The clean reads were blasted against a microbial database downloaded from the National Center for Biotechnology Information. Finally, microbial information at the species level can be obtained.

Recurrent acute exacerbations and hospitalization are crucial characteristics of severe COPD. We observed 56 patients with COPD admitted to the hospital from July 2020 to January 2024, most of whom were hospitalized with a combination of lung infections. During hospitalization, 94.6% of patients received antibiotics, 21.4% underwent mechanical ventilation, and 64.2% faced hospitalization costs over 20,000 RMB. Additionally, 51.8% were hospitalized multiple times during follow-up, with a mortality rate of 5.4% ( Supplementary Table 2 ). The mNGS results suggested that Pseudomonas aeruginosa and Human gammaherpesvirus 4 were the most common bacterial and viral pathogens ( Supplementary Figure 1 ).

We used eQTLs from whole blood and lung tissues to intersect with druggable genes to obtain druggable eQTLs. Finally, cis-eQTLs for 2529 druggable genes from eQTLGen’s whole blood and 1042 druggable genes from GTEx v.8’s lung tissues were extracted, both of which were used for SMR analyses with COPD and lung function as outcomes.

Among the genes identified by SMR/HEIDI and colocalization analysis, we obtained cis-pQTL data for 9 genes from the proteomics summary data of whole blood samples provided by Ferkingstad et al.

Our exposure and outcome samples were both predominantly European, so it can be assumed that the genetic variants used as instrumental variables (cis-eQTLs and cis-pQTLs) had a consistent role between the two different sample sets. The exposure and outcome samples were from completely independent studies, so there was no participant overlap.

We performed SMR analysis using doctor-diagnosed COPD GWAS summary data from the FinnGen consortium. Initially, we identified 3 genes (SLC22A5, C2, and GPC2) whose expression levels in whole blood were found to be correlated with doctor-diagnosed COPD (corrected for P-value by Bonferroni, P_HEIDI≥0.01, see Supplementary Table 3 for detailed data).

We utilized the FEV1/FVC<0.7 GWAS summary data from the UK biobank, and identified 23 genes whose expression levels in whole blood or lung tissues were associated with spirometry-defined COPD (corrected for P-value by Bonferroni, P_HEIDI≥0.01, see Supplementary Table 3 for detailed data).

The SMR analysis using FEV1 GWAS summary data identified 36 genes that were significantly associated with FEV1. Out of these genes, six were also found to be associated with spirometry-defined COPD.

Furthermore, colocalization analysis identified 2, 13, and 16 genes that shared the same SNPs (PP.H4 > 0.8 and/or PH4/PH3 > 2, Supplementary Table 4 ) in the results of doctor diagnosed COPD, FEV1/FVC<0.7 and FEV1, respectively ( Figure 3 ).

We performed eQTL two-sample MR ( Table 1 ; Supplementary Table 5 ) on druggable genes that had passed the colocalization analysis to assess the consistency of direction and statistical significance.

We further performed protein SMR analysis ( Table 1 ; Supplementary Table 5 ). Of the identified genes, 1 (GPC2) was associated with doctor-diagnosed COPD, 10 (TESK2, AKR1A1, LMCD1, ERBB3, PSMA4, MAST2, MMP15, MRC2, BMP4, and LTK) with FEV1/FVC <0.7, and 11 (APH1A, CHI3L2, PTK7, RASGRP3, HSP90AA1, TESK2, ADAM33, COL15A1, LINGO1, LTBR, and SULT1A1) with FEV1 in protein SMR results, without contradicting the expression SMR results. These genes were categorized as high priority.

In addition, we tested the underlying assumptions and sensitivity analyses for two-sample MR ( Supplementary Table 6 ). The F-Statistic for the gene instrumental variables were all greater than 10. The Egger test indicated no evidence of horizontal pleiotropy for 16 genes and suggested that the IVW results could be used as the primary ones. The Steiger test confirmed the correct causal direction.

Among the high-priority drug-targeted genes, MMP15 is associated with Montelukast and COL-3 (NSC-683551), while PSMA4 is linked to MARIZOMIB. These three drugs are small molecule inhibitors, which may reduce the risk of spirometry-defined COPD. On the other hand, the APH1A gene, which codes for the gamma-secretase subunit APH-1A, appears to negatively impact FEV1 (β<0 in MR analysis). Several existing inhibitors against this target (e.g., TARENFLURBIL, BEGACESTAT, SEMAGACESTAT, NIROGACESTAT, RG-4733, and AVAGACESTAT) may counteract this harmful effect. Additionally, high expression of the ERBB3 gene may lower the risk of spirometry-defined COPD. However, various existing antitumor drugs, which are inhibitors of this target, may increase the risk of spirometry-defined COPD ( Table 2 ).

The PheWAS analysis revealed that increased GP9 expression may have adverse effects on platelets in addition to possibly increasing the risk of spirometry-defined COPD. Conversely, increased ERBB3 expression is associated with a reduced risk of spirometry-defined COPD and inflammation suppression. Moreover, increased MRC2 expression affects cardio metabolism and the nervous system ( Figure 4 ).