Correlation-based network integration of lung RNA sequencing and DNA methylation data in chronic obstructive pulmonary disease

2.1. Case-control study population

In the present study, we analyzed a COPD case-control dataset derived from the Lung Tissue Research Consortium (LTRC). We restricted our analyses to 792 subjects −446 COPD cases and 346 controls − for which both RNA-seq and DNA methylation data were available. Table 1 and Fig. 1 report, respectively, the spirometry and the clinical and demographic features of all COPD and control samples we considered. Interestingly, this cohort showed reasonable consistency in key clinical and demographic features (Fig. 1). Indeed, COPD and control groups had similar values of age and body mass index (BMI) as well as comparable gender distributions (Fig. 1). Regarding smoking status, most COPD subjects (83 %) were former smokers, with a very low percentage of current smokers (6 %) and about 11 % missing values. For the control group, we found again a higher percentage of former smokers (57 %) than current smokers (5 %), and about 38 % missing values (Fig. 1). COPD subjects had greater average smoking intensity than controls, with more than twice the average number of pack-years of cigarette smoking (Fig. 1).

2.2. Study design

The overall study design is shown in Fig. 2. Specifically, we applied SWIM software [9,10] separately to RNA-seq and DNA methylation data and then integrated the output of SWIM to build a coupled network of the preserved interactions between the two single-omics networks. By performing a differential correlation analysis, the links of the coupled network were categorized into four types (i.e., +/+, −/−, +/−, −/+), depending on the sign of pairwise correlations in the single-omics correlation networks. The +/+ type referred to gene pairs that were positively correlated in both single-omics correlation networks; the −/− type referred to gene pairs that were negatively correlated in both single-omics correlation networks; the ± and −/+ types referred to gene pairs whose sign of the correlation coefficient was opposite in the single-omics correlation networks. Further details of these edge types can be found in the section Edge classes of the coupled network of Supplementary Materials.

2.3. SWIM application on RNA-seq and DNA methylation lung tissue data

Starting from 16,251 genes shared between RNA-seq and DNA methylation data based on gene symbols, the SWIM-based exploratory data analysis identified 2,442 differentially expressed genes (DEGs) and 990 differentially methylated genes (DMGs) in COPD vs. control lung tissue samples. Among the 2,442 DEGs, we found that 1,205 (49.3 %) were upregulated and 1,237 (50.7 %) were downregulated in COPD cases (Fig. 3a, first panel). Among the 990 DMGs, we found that 338 (34.2 %) were relatively hypermethylated and 652 (65.8 %) were relatively hypomethylated in COPD cases (Fig. 3b, first panel). Then, the SWIM-based network analysis identified: i) an RNA-based correlation network that contained 1,856 nodes and 113,235 edges, including 801 party hubs, 717 date hubs, and 117 fight-club hubs (Fig. 3a, second panel); and ii) a methylation-based correlation network that contained 892 nodes and 94,876 edges, including 486 date, 252 party, and 108 fight club hubs (Fig. 3b, second panel). To detect the community structure of each single-omics correlation network, SWIM applied the k-means clustering algorithm determining the optimal number of clusters by the elbow method [33] (Fig. 3a–b, third panels). The RNA-based network consisted of four clusters, varying in size from 1,089 nodes for cluster 1, 263 nodes for cluster 2, 295 nodes for cluster 3, and 209 nodes for cluster 4 (Fig. 4, upper panels); whereas the methylation-based network consisted of three clusters, varying in size from 157 nodes in cluster 1, 260 nodes in cluster 2, and 475 nodes in cluster 3 (Fig. 4, upper panels). Then, to assign a topological role to each node of the single-omics correlation networks based on their inter- and intra-cluster interactions, SWIM drew the heat cartography map for the RNA-seq data (Fig. 3a, fourth panel) and for DNA methylation data (Fig. 3b, fourth panel), where party, date, and fight-club hubs are identified by red, orange, and blue coloring, respectively. SWIM identified 112 switch genes in the RNA-based network (Fig. 4, bottom left panel), and most of them were downregulated in COPD cases (90/112, 80.3 %). In addition, SWIM identified 101 switch genes in the methylation-based network and all of them were hypermethylated in COPD cases (Fig. 5, bottom left panel).

In the RNA-based network, the switch genes were distributed over three clusters with a particular abundance in cluster 1 (59/112, 52.7 %) and cluster 3 (50/112; 44.6 %) (Fig. 4, bottom left panel). In the DNA methylation network, the switch genes are almost all distributed in cluster 1 with a few switch genes in cluster 3 (Fig. 5, bottom left panel). We focused on those switch genes showing a coherent pattern of regulation (i.e., only upregulated or downregulated) and falling in the same network cluster, i.e., the switch genes of cluster 3 of the RNA-based network (Fig. 4, bottom left panel) and the switch genes of cluster 1 of the methylation-based network (Fig. 5, bottom left panel). Further, we looked for negative nearest neighbors of the selected switch genes of both single-omics correlation networks (Fig. 4, Fig. 5, bottom right panels). Interestingly, cluster 2 of both single-omics correlation networks was composed exclusively of the negative nearest neighbors of the switch genes (Fig. 4, Fig. 5). By performing a functional enrichment analysis, we found that the selected switch genes of the RNA-based network were enriched in gene expression regulation activities related to generic transcription factors and RNA-pol II functionalities (Supplementary Table 1), whereas their negative nearest neighbors were enriched in Immune system, Infectious disease, and Inflammasome related pathways according to the Reactome database [34] (Supplementary Table 1). These findings suggested a putative regulatory mechanism where the downregulation of switch genes could be related to an activation of immune and inflammatory components, which is likely important for COPD pathogenesis [2]. Regarding the list of the selected switch genes of the methylation-based network, we observed no statistically significantly enriched pathways (Supplementary Table 2). Nevertheless, their negative nearest neighbors were enriched in both adaptive and innate immune system-related pathways (Supplementary Table 2).

2.4. Coupled network

The RNA-seq and DNA methylation networks built by SWIM were then integrated to obtain the so-called coupled network, where nodes are the genes both differentially expressed and differentially methylated in COPD cases with respect to controls, and a link between two nodes occurs if they were highly correlated or anti-correlated in both single-omics RNA and methylation-based networks (Fig. 6 and Supplementary Table 3). The coupled network contained 97 nodes and 317 links, which were categorized into four types (i.e., +/+, −/−, +/−, −/+), depending on the sign of the correlation coefficient between each gene pair in the starting single-omics correlation networks (Fig. 6).

We observed a predominance of +/+ edge types in the coupled network (Fig. 6), meaning that among the differentially expressed and differentially methylated genes in COPD lung tissue there exist strong positive correlations preserved between transcriptomics and epigenomics layers. Next, we performed a pathway overrepresentation analysis of the nodes of the coupled network, using the Reactome database [34]. We observed a strikingly significant enrichment in Immune System components, both innate and adaptive (Supplementary Fig. 1). In addition, we also found a significant enrichment in cytokine signaling, including multiple interleukin family signaling pathways (Supplementary Fig. 1) such as the interleukin-2 family signaling pathway that was recently associated with the severity of COPD [35].

2.5. Comparison between single vs multi-omics analysis

We compared the results of the functional enrichment analysis performed on the nodes of the single-omics networks and of the coupled network to reveal the specific enriched pathways for each network (Supplementary Fig. 2). Interestingly, although the number of nodes in the coupled network was about 10 times less than the number of nodes in the single-omics networks, the coupled network nodes were specifically enriched in adaptive immune system, tumor necrosis factors (TNFs), and a variety of interleukin pathways (Supplementary Fig. 2) related to both type I inflammation (IL-1, IL-18, and TNFs) and type II inflammation (IL-5). Moreover, to quantify the specific enrichment in immune-related pathways of the nodes of each network, we computed the immune specificity score (see Methods), and we observed that the coupled network had the highest immune specificity compared to the single-omics networks (Supplementary Fig. 3).

2.6. Comparison with other multi-omics integration methods

We compared the results of our analysis with another multi-omics integration method, i.e., the MOFA (Multi-Omics Factor Analysis) algorithm [31]. We chose MOFA for this comparison because it is a general and widely used framework for the unsupervised integration of multi-omics data. Moreover, in the reference describing this method [31], the authors demonstrated over a range of simulations that MOFA tends to be more accurate and more computationally efficient compared to other existing multi-omics integration methods (e.g., iCluster). Here, we applied MOFA by considering as input the RNA-seq matrix of DEGs and the DNA methylation matrix of DMGs. MOFA identified 10 Factors and, among them, we focused on Factors that explained more than 5 % of the variance of both single-omics data types, i.e., Factor 1 and Factor 2 (Supplementary Fig. 4). We also verified that Factor 1 and Factor 2 did not show strong correlations with any measured demographic or smoking covariates (Supplementary Fig. 4b). Then, we considered the genes with the highest weights (i.e., greater than 75th percentile of their distribution) in the RNA-seq and DNA methylation data both on Factor 1 and Factor 2. We found a total of 28 and 17 genes in common between the two molecular layers for Factor 1 and Factor 2, respectively. Most of the 17 genes extracted from Factor 2 (n = 13/17, 76,5 %) were included in the 28 genes extracted from Factor 1. Of interest, 27 of these 28 genes were included in the coupled network, and they were found to be enriched in immune-related pathways (i.e., Immune System, Cytokine Signaling in Immune system, and Innate Immune System) with an immunity specificity score of 3.5 (Supplementary Fig. 3). The additional genes composing the coupled network were found to be more specifically associated with the Adaptive Immune System, thus leading to a higher immunity specificity score (Supplementary Fig. 3). Taken together, these results suggest that our network-based methodology allowed us both to unveil additional genes with specific roles in the adaptive immune responses and to probe putative molecular relationships among them.

3.1. Inspection of the coupled network

In the coupled network, we observed a prevalent sign preservation of the correlation coefficients between the gene pairs of the coupled network moving between transcriptomics and epigenomics layers (Fig. 6). Specifically, most of the gene pairs were found to belong to the +/+ edge class (Fig. 6), meaning that the two genes were positively correlated both in their gene expression profiles and in their methylation profiles. Genes that are co-methylated are likely to be co-expressed in response to an external stimulus [36] (Supplementary Fig. 5), like smoke exposure, which is a known risk factor for COPD affecting the DNA methylation status of genes in small airway epithelial or other lung cells [14]. It is worth noting that the correlation provides information exclusively on the shape and not on the direction of the expression (up/down) or methylation (hypo/hyper), which should be independently investigated (see section Edge classes of the coupled network of Supplementary Materials). We found that most genes (78/88, 89 %) of the coupled network interacting by the +/+ edge class showed an opposite direction between gene expression and methylation of the promoter: 59/88 (67 %) are up-regulated and hypomethylated in COPD patients with respect to controls, whereas 19/88 (22 %) are down-regulated and hypermethylated in COPD patients with respect to controls (Supplementary Table 4). Conversely, the remaining genes (10/88, 11 %) were found to be downregulated and hypomethylated in COPD with respect to controls (Supplementary Table 4).

We also observed changes in the sign of correlation coefficients between some gene pairs of the coupled network ( ± and −/+ edge types) moving from the DNA to RNA layers (Fig. 6), suggesting the potential existence of molecular mechanisms other than DNA methylation regulating gene expression (Supplementary Fig. 6), e.g., post-transcriptional (miRNA targeting, mRNA degradation in cytoplasm), translational (mRNA translational rate), and post-translational (protein modifications and degradation) mechanisms [37]. In these cases, additional molecular data (e.g., miRNA-seq, CHIP-seq, or protein quantification) could allow investigation of other possible molecular mechanisms responsible for these scenarios.

For each edge class of the coupled network, we highlighted the most connected nodes, i.e., with an internal degree above the 90th percentile (Supplementary Table 4). Among them, we found GPR132, TNFSF8, BATF, AIM2, RHOH, IL10, EVI2A, TRPM2 and IL2RA for the +/+ edge type; KIAA0408 and SLC15A2 for the −/− edge type; FFAR2 and MGAT5B for the ± edge type; and PRDM11 for the −/+ edge type.

The selected genes falling in the +/+ class were significantly up-regulated and their promoter regions were typically significantly hypomethylated in COPD patients with respect to controls, whereas the selected genes falling in the −/− class were significantly down-regulated and their promoter regions were significantly hypermethylated in COPD patients with respect to controls (Supplementary Table 5). A different scenario was observed for FFAR2, MGAT5B and PRDM11, because their mRNA expression did not have opposite direction with respect to the methylation status of their promoter: FFAR2 was significantly up-regulated and its promoter was significantly hypermethylated in COPD subjects compared to controls, whereas MGAT5B and PRDM11 were significantly down-regulated and their promoters were significantly hypomethylated in COPD subjects compared to controls (Supplementary Table 5). This scenario could be related to the heterogeneous activity of DNA methylation of both silencing and triggering gene expression [38].

3.2. Recognition of genes involved in inflammatory processes and immune system regulation

The activation of both innate and adaptive immune systems, leading to the recruitment and activation of inflammatory cells into the lung parenchyma followed by proinflammatory cytokine proteolytic activation and pyropotosis of alveolar cells, has been suggested to provide important pathomechanisms in COPD [39]. Consistent with these observations, we found that nodes of the coupled network are enriched in genes participating in both innate and adaptive immune systems as well as in cytokine signaling mechanisms (Table 2). Among others, TNFSF8, BATF, IL10, and IL2RA are involved in cytokine signaling, and their impaired activity could be related to an increased airway inflammatory response in COPD; AIM2 and TRPM2 are involved in innate immune system mechanisms such as neutrophil degranulation and inflammasome that can be related to COPD pathomechanisms; FFAR2 is involved in the regulation of the immune state through the gut-lung axis [40]; and PRDM11 and RHOH are not directly involved in immune system canonical pathways, but recent laboratory experiments related their altered function to impaired airway inflammation (Table 2). It is also worth noting that some of the highlighted genes (i.e., KIAA0408, EVI2A and SLC15A2) are not currently known to participate in immune-related pathways or in COPD-related mechanisms. However, they were found to interact with important inflammatory factors (e.g., IL10, IL2RA, and TNFSF8) at both transcriptomics and epigenomics layers (Supplementary Table 3). Identifying novel nodes not previously known to be involved in COPD pathogenesis and finding their interactions with known COPD-related genes are important benefits of network applications. In fact, this result could not have been obtained by applying methods that operate at the level of individual genes disregarding their interactions, such as differential gene expression/methylation analysis or MOFA.

Finally, we compared the coupled network with the two single-omics networks to test whether our method can bring advantages in the investigation of COPD-related pathomechanisms. We demonstrated that our coupled network encompasses genes that were more specifically involved in the adaptive immune system compared to the two single-omics networks (Supplementary Fig. 3).

3.3. Limitations and future directions

In this study, we have identified several computational and biological clues that could contribute to a clearer understanding of the putative mechanisms underlying the development and progression of COPD, which may ultimately lead in the future to better treatment of patients affected by this devastating lung disease. Although our approach provided interesting results, it is not free of limitations.

A crucial point of our analysis is the computation and management of the correlation matrix needed to build each single-omics correlation network. As discussed in the section Identification of highly modulated genes of Supplementary Materials, considering the entire transcriptome is not feasible since it would require a very high computational complexity. Setting a FC threshold to reduce the number of input features is currently an unavoidable and widely used step. Undoubtedly, an open challenge to be faced is to make our method more efficient in handling and analyzing larger datasets.

Another limitation is that our method is based on correlations that are “associations” and do not necessarily imply “causal” relationships. Indeed, we built a correlation network for each omics data type (i.e., transcriptomics and epigenomics) and then we investigated the correlation changes moving from one omics type to the other. Analyzing changes in these pairwise interactions may contribute useful insights into their biological nature. Certainly, as demonstrated by several recent studies [17], adding information by including other omics layers (e.g., proteomics) may further improve the understanding of the system under study, and thus help in bridging the gap from correlation to causation. However, there are several challenges to be considered when the number of omics data types increases such as the difficulty of retrieving complete multi-omics datasets (i.e., all the omics must be obtained for the same sets of samples), the handling of missing values, and the heterogeneity across the different data modalities [55].

In principle, our method can be leveraged to integrate more than two data types simultaneously. For each single-omics, we first build the corresponding correlation network and then the coupled network considering the common pairwise interactions among all the single-omics correlation networks. In the coupled network, the edges’ weight is the class, given by the sign of the pairwise correlation in each single-omics correlation network (e.g., +/+, −/−, +/−, −/+). Generally, considering N omics data leads to 2^N different edge classes in the coupled network. By definition, nodes in the coupled network are common to all single-omics correlation networks. This implies considering different data types but with the same basic units of analysis (e.g., genes) measured on the same samples (like transcriptome, epigenome, and proteome). Thus, a natural limitation of our method is the need to have a common unit of analysis.

A further limitation of our study is that the analysis of the genes' methylation levels was performed considering a narrow region of the promoter. A future extension of the methodology would be to analyze genes' methylation levels considering different Open Reading Frame (ORF) regions, gene bodies, 3′UTR and 5′UTR. Even if gene promoter methylation is a well-established repressive gene regulatory program, recent discoveries have proved that DNA methylation in other genomic regions may increase gene expression [38,56]. Therefore, an important future perspective would be studying the influences of DNA methylation occurring at these other genomic regions in the context of disease (e.g., COPD), compared with the well-established gene promoter methylation regulatory program. In addition, we did not consider QTL effects and splice variants in the analyses and thus they could be included in future developments. The integration of expression and/or methylation QTL could supply information about genomic location of the putative regulatory loci that could influence our integrated correlation network. Moreover, the inclusion of splice variants in our methodology could enrich the gene expression dataset with a supplemental layer of gene regulation.

5.1. Data collection and processing

Lung tissue samples were provided by the Lung Tissue Research Consortium (LTRC). LTRC participants provided written informed consent for excess lung tissue obtained during clinically indicated thoracic surgery procedures to be used in research studies. LTRC subjects were defined as COPD cases if they had a ratio of forced expiratory volume in 1 s (FEV₁) to forced vital capacity (FVC) < 0.70 together with FEV₁ < 80 % of predicted normal spirometry using Hankinson equations [32], consistent with GOLD grade 2–4 spirometry [57]. LTRC subjects were defined as controls if they had normal spirometry, i.e., FEV₁/FVC ratio ≥0.70 and FEV₁ ≥ 80 %. Subjects with pathological diagnoses of interstitial lung disease were excluded. In this work, we considered in total 792 samples −446 COPD cases and 346 controls − retrieved from the freeze dated 30NOV2022, for which RNA-seq and DNA methylation data were available through the NHBLI Trans-Omics for Precision Medicine (TOPMed) Program. Table 1 and Fig. 1 report, respectively, the spirometry and the clinical and demographic features of all COPD and control samples included in this study. Notably, the smoking status variable includes current smokers (i.e., subjects who were smoking within one month of surgery), former smokers (i.e., subjects who quit smoking at least one month prior to surgery) and not reported (i.e., subjects with smoking history not available).

5.2. RNA-seq data preprocessing and normalization

The TOPMed program (University of Washington sequencing center) performed the mRNA sequencing experiments. Additional details on the RNA-sequencing processing steps are reported in Supplementary Materials. In general, RNA-seq data could be affected by technical variation, called batch effects, influencing the downstream biological analysis. Standard normalization methods are not able to address the batch effects, therefore methodologies such as ComBat and ComBat-seq were developed to overcome these limitations [58,59]. In particular, Combat-seq can adjust the RNA-seq data for batch effects by fully modeling the distribution of RNA-seq count data through a negative binomial distribution and preserving its integer nature after the adjustment. In the present study, the rationale behind using Combat-seq is to correct the technical variability of the RNA-seq experiments due to the use of different sample sequencing plates. For this correction, we used sequencing plates, COPD affection status, age, sex, and smoking status as covariates in the model. Values for smoking status that were not reported have been considered as missing values.

Afterwards, a normalized matrix, including autosomal and X/Y genes, was produced by using the EDASeq R package [60]. This package allowed us to perform within-sample and between-sample normalization procedures. We selected the GC content within-sample normalization procedure since it was able to reduce the dependencies between the gene-level read counts and the sequence composition in GC that varies from gene to gene. Indeed, genes with higher GC content sequences tend to be underestimated in the sequencing process. In detail, the GC content within-sample normalization procedure stratifies the genes in equally sized bins based on GC-content and then matches the quantiles of the gene count distributions across bins. Then, we chose Upper Quantile (UQ) for between-sample normalization since it was proven to perform well in differential expression and correlation analysis of transcriptomic data [61,62]. The UQ normalization forces the upper quantile of count distributions to be the same among different samples.

5.3. DNA methylation data preprocessing and normalization

Methylation data using DNA extracted from lung tissue samples was obtained using the Illumina EPIC Array (850K) through TOPMed. The DNA methylation data was normalized by using a functional normalization and regression on correlated probes procedure. In the preprocessing phase, we filtered out probes that failed at detection p-value <0.05 in >25 % samples (n = 738). SNP-under-probe effects (n = 11,679), cross-reactive probes (n = 44,493), and multimodal probes were detected using Hartigan's dip test (n = 9,672). The normalized matrix contained DNA methylation Beta-value of 853,552 CpGs sites of both autosomal and X/Y chromosomes for each subject. The Beta-value, ranging from 0 to 1, indicates the percentage of methylation of CpG sites. We ended with 801,234 probes after filtering.

Here, we considered DNA methylation of the gene promoter region, which is known to be related in many cases to the repression of gene transcription [63]. The EPIC array platform was able to measure the DNA-methylation of two gene promoter regions: 1500 and 200 bp upstream to the genes Transcription Start Site (TSS1500 or TSS200). To obtain a gene promoter DNA-methylation measure, we pick the TSS200 probes' B-values for each gene and averaged them (see section Gene level DNA methylation data processing of Supplementary Materials for the details). The rationale for excluding the TSS1500 probes while selecting the TSS200 probes was related to several reasons: i) non-reliability in aggregating the two kinds of probes’ B-values due to low positive correlations between them in many gene promoter regions (Supplementary Fig. 7); ii) demonstrated higher correlation of the DNA-methylation levels between intraregional CpG sites (e.g., TSS200) [38]; and iii) proximity to the Transcription Starting Site and a higher likelihood to be recognized by DNA-methylation binding complexes, which in turn can influence gene expression [38].

Next, we performed a transformation of the DNA methylation data. Indeed, the averaged Beta-value has been shown to have a bounded range that violates the Gaussian distribution assumptions and showed heterogeneity of variation through both low and high ranges [64]. These aspects violate the assumptions of many statistical tests employed for downstream analyses of differential methylation and pairwise correlation [64]. A measure that can be used to overcome these limitations is the M-value, i.e., the log2 ratio of the intensities of methylated probe versus unmethylated probe. M-values do not have a bounded range and show homogeneity of variation through their entire range, which better fit the assumptions of many statistical tests [64]. Thus, we choose to transform the gene level DNA methylation data from Beta-values into M-values to avoid false positives when applying statistical tests for DMGs and correlation analysis.

Finally, the gene level DNA methylation data was corrected for batch effects related to the technical and clinical variability across sample plates by applying the ComBat function of the sva package [58], using array plates, COPD affection status, age, sex, and smoking status as covariates in the model. Values of the smoking status that were not reported have been considered as missing values.

5.4. SWIM-based single-omics analysis

SWIM is a freely downloadable network‐based analytical tool, developed both in MATLAB [9] and in R language [10] to predict important (switch) genes that are associated with large changes in phenotypes. In the present study, SWIM was separately applied on the normalized and corrected RNA-seq matrix and the normalized and corrected DNA methylation matrix. To run SWIM, initial model parameter settings are required, including decision thresholds, statistics, and graphics details. In line with a reverse engineering approach, this initial configuration can be a posteriori fine-tuned until an optimal definition of the switch genes is reached.

SWIM first computed the differentially expressed genes (DEGs) or the differentially methylated genes (DMGs) between COPD and control subjects. Both in the DEGs and DMGs analysis, the standard significance level of 0.05 was set as threshold for the adjusted p-value (FDR correction). Then, to reduce the number of possible false positive results in the list of switch genes, among DEGs and DMGs, SWIM searched for genes that are highly modulated by setting a threshold on the absolute value of their fold-change (FC). All details for a reasonable choice of this threshold can be found in the section Identification of highly modulated genes of Supplementary Materials.

After the identification of highly modulated genes (Fig. 3a–b, first panels), SWIM built a Gene Co-expression Network by computing (positive and negative) correlations between the profiles of each gene pair in RNA or DNA methylation layers. Specifically, SWIM implements a hard thresholding approach to build the correlation networks, where nodes are DEGs/DMGs, and a link occurs if their expression/methylation profiles are highly correlated or anti-correlated. In general, the hard thresholding approach requires the choice of a given correlation threshold and creates networks where two nodes are linked if their absolute correlation exceeds that threshold. The use of this (hard) threshold in SWIM is necessary to avoid ending up with a fully connected network and to remove spurious relationships, thus focusing only on significant associations between highly correlated nodes.

Regarding the choice of the correlation threshold, it should be selected to guarantee an appropriate balance between the number of edges and the number of connected components of the network: the number of edges should be as small as possible to have a manageable network (pointing towards a higher threshold) and the number of connected components should be as small as possible to preserve the integrity of the network (pointing towards a lower threshold). For the lung RNA-seq and DNA methylation data in LTRC, this balance was achieved by setting the correlation threshold to 0.52 and to 0.53, respectively.

Next, SWIM classified each network hub (i.e., nodes with degree ≥5 [65]) as party, date or fight-club on the basis of the Average Pearson Correlation Coefficient (APCC) between its expression or DNA methylation profile and that of its first nearest neighbors. The date/party/fight-club hub trichotomy is generally mirrored by a trimodal distribution of APCC values (Fig. 2a–b, second panels) and corresponds to three different node roles in the correlation network: date hubs show a positive but relatively low APCC value (i.e., low co-expression/co-methylation with their partners); party hubs show a positive and high APCC value (i.e., high co-expression/co-methylation with their partners); and fight-club hubs show a negative APCC value (i.e., inversely correlated with their partners).

To identify communities in the correlation networks, SWIM first searched for clusters using the k-means algorithm and then evaluated the quality of clusters by minimizing the Sum of the Squared Error (SSE), depending on the distance of each feature to its closest centroid. However, the k-means algorithm needs to define a priori the number of clusters to be detected. To choose the most appropriate number of clusters, SWIM leverages the scree plot approach, where this number (3 in Fig. 3a–b, third panels) is suggested by the elbow's position in the SSE plot computed as a function of the number of clusters [33].

To assign a role to each node, SWIM drew a heat cartography map (Fig. 3a–b, fourth panels) by evaluating two parameters related to their inter- and intra-modular connections: the clusterphobic coefficient, which depends on the links of each node to nodes outside its own cluster; and the within-module degree, which depends on how “well-connected” each node is within its own cluster. Nodes with much more external than internal links have a high value of the clusterphobic coefficient and are called connectors, whereas nodes with a high value of the within-module degree are hubs within their community and are called local hubs.

Finally, switch genes of each single-omics network were defined as a subset of fight-club nodes with the following characteristics: they were network connectors that mainly interact outside their own cluster; they were not local hubs; and they were mainly anti-correlated with their interaction partners.

It is worth to stress again that a good choice of all above-mentioned thresholds should be reflected in an appropriate definition of the switch genes that should guarantee that their properties are verified. If this does not happen, each threshold must be varied within a range that respects the rationale behind its definition.

5.5. Coupled network construction

The coupled network was built by integrating the two single-omics correlation networks for RNA-seq and DNA methylation, thus preserving the nodes and the edges in common between the two single-omics correlation networks. More precisely, the two single-omics correlation networks were integrated by selecting all gene pair interactions (edges) occurring both in RNA and methylation-based networks. This allowed us to construct a new integrated network, named the coupled network, where nodes are genes present in both single-omics correlation networks (and thus differentially expressed and differentially methylated in COPD cases compared to controls); ii) and a link between two genes occurs if a correlation exists between their expression and methylation profiles simultaneously.

5.6. Pathway enrichment analysis

The functional enrichment analysis of a given gene list was performed by querying the Reactome database [34] through the EnrichR web tool [66]. Reactome is a free, peer-reviewed, and large pathway database which provides manually curated molecular information about physiological and pathological processes in humans, including both hereditary and acquired diseases. It systematically associates human proteins to their molecular functions, thus providing a resource for the visualization, interpretation and analysis of pathway knowledge as well as for the discovery of novel functional relationships in data such as gene expression profiles [34]. In this study, we used the EnrichR web tool [66] to compute gene set enrichment in Reactome pathways of a chosen input gene list. For pathway enrichment analysis, all of the genes annotated in the Reactome database were considered as background and the value of 0.05 was set as threshold for the adjusted p-value (FDR correction) associated to each pathway.

5.7. MOFA analysis

MOFA (Multi‐Omics Factor Analysis) is an unsupervised method for the integration of multi-omics data which can be viewed as a generalization of principal component analysis [31]. Indeed, given multiple omics data types on the same or on partially overlapping sets of samples, MOFA infers an interpretable low-dimensional representation in terms of a small number of latent factors that should capture the relevant signal in the input data. It performs a dimensionality reduction from a complex high dimensional matrix (i.e., the matrix of the observations), whose size scales with the number N of the features, to a low dimensional latent space (i.e., the matrix of the factors), whose size scales with the number M of latent factors. M must be much lower than N, but there is no general rule regarding the optimal value for M. In the present study, we set the number of latent factors M equal to 10 and the rationale for choosing this value was to follow the guidelines provided by developers of the MOFA method [31]. Indeed, they stated that the optimal number of factors depends on the aim of the analysis, the dimensions of the assays, and the complexity of the data: in general, if the aim is to identify the major sources of biological variation, like in this study, one would typically choose the top 10 factors; in other tasks, such as imputation of missing values, even small sources of variation can be important and hence models should be trained with a larger number of factors.

Next, the learned factors can be used for a variety of downstream analyses, including visualization, clustering and classification of samples, data imputation, and factor annotation using gene set enrichment analysis. In the present study, we used MOFA2 R package (version 1.8.0) to analyze the RNA-seq and DNA methylation data. The two data matrices were scaled to have a unit variance. The MOFA model training was started with 10 factors. Other parameters were set as default.

5.8. Immune specificity score

In order to quantify the specific involvement in immune-related pathways of a given gene set, we introduced a score named “immune specificity”. Specifically, after a pathway enrichment analysis, the immune specificity was defined as follows:

where $N$ refers to the total number of immune-related pathways in which the input gene lists were enriched (FDR <0.05). This means that the higher the score, the more the selected genes are specifically involved in pathways related to the immune response, which is known to play a crucial role in the development and progression of COPD [[67], [68], [69], [70]]. As immune-related pathways, we considered the root of the pathways associated to the immune system (i.e., Immune System, Cytokine Signaling in Immune system, Adaptive Immune System and Innate Immune System) according the Reactome database [34].