Integrative network analysis suggests prioritised drugs for atopic dermatitis

doi:10.21203/rs.3.rs-3579100/v1

Download PDF

Research Article

Integrative network analysis suggests prioritised drugs for atopic dermatitis

https://doi.org/10.21203/rs.3.rs-3579100/v1

This work is licensed under a CC BY 4.0 License

Journal Publication

published 16 Jan, 2024

Read the published version in Journal of Translational Medicine →

Version 1

posted

You are reading this latest preprint version

Background

Atopic dermatitis (AD) is a prevalent chronic inflammatory skin disease whose pathophysiology involves the interplay between genetic and environmental factors, ultimately leading to dysfunction of the epidermis. While several treatments are effective in symptom management, many existing therapies offer only temporary relief and often come with side effects. For this reason, the formulation of an effective therapeutic plan is challenging and there is a need for more effective and targeted treatments that address the root causes of the condition. Here, we hypothesise that modelling the complexity of the molecular buildup of the atopic dermatitis can be a concrete means to drive drug discovery.

Methods

We preprocessed, harmonised and integrated publicly available transcriptomics datasets of lesional and non-lesional skin from AD patients. We inferred co-expression network models of both AD lesional and non-lesional skin and exploited their interactional properties by integrating them with a priori knowledge in order to extrapolate a robust AD disease module. Pharmacophore-based virtual screening was then utilised to build a tailored library of compounds potentially active for AD.

Results

In this study, we identified a core disease module for AD, pinpointing known and unknown molecular determinants underlying the skin lesions. We identified skin- and immune-cell type signatures expressed by the disease module, and characterised the impaired cellular functions underlying the complex phenotype of atopic dermatitis. Therefore, by investigating the connectivity of genes belonging to the AD module, we prioritised novel putative biomarkers of the disease. Finally, we defined a tailored compound library by characterising the therapeutic potential of drugs targeting genes within the disease module to facilitate and tailor future drug discovery efforts towards novel pharmacological strategies for AD.

Conclusions

Overall, our study reveals a core disease module providing unprecedented information about genetic, transcriptional and pharmacological relationships that foster drug discovery in atopic dermatitis.

Systems Biology

Computational Biology

Drug Discovery, Design, & Development

atopic dermatitis

network analysis

disease module

biomarkers

drug discovery

Atopic dermatitis (AD) is one of the most common chronic inflammatory skin diseases characterised by extended skin lesions leading to significant discomfort and reduced quality of life. Its pathophysiology results from a complex interaction of genetic and environmental factors leading to epidermal dysfunction and cutaneous inflammation, which is driven by an abnormal T-cell activation and an increased production of inflammatory cytokines [1].

Considerable efforts have been made to identify the cutaneous molecular alterations associated with AD. Advancements in omics technologies and the massive amounts of data currently available revealed a large number of deregulated genes and pathways in AD, allowing the discovery of novel biomarkers, disease subtypes/endotypes and therapeutic targets [2–5]. Manipulating and integrating these data into biologically meaningful results is still challenging especially in the light of the high disease and study heterogeneity. For these reasons, an exhaustive characterization of the pathogenetic mechanisms and therapeutic opportunities of AD have not been achieved yet [6–9].

Network science is at the forefront to meet this challenge and its potential has been widely demonstrated [10, 11]. Network modelling of biomolecular processes enables the full integration and exploitation of heterogeneous (multi-)omics datasets, uncovering key molecular mechanisms responsible for the onset of several human diseases and facilitating the design of new and more effective treatments. Although today, several effective and specific treatments with tolerable side-effects are available for patients with AD [12], challenges such as the high disease heterogeneity as well as the effect of multimorbidity are demanding the development of further treatment solutions [13, 14]. To date, the potential of network science has not been yet exploited to facilitate drug discovery in AD.

In this study, we retrieved, curated, analysed and integrated 13 public transcriptomics datasets encompassing 337 samples of lesional and 542 samples of unaffected skin from patients affected by AD. We built transcriptome-wide co-expression network models of lesional and non-lesional AD skin. By combining network modelling of the transcriptomics data with multiple layers of curated a priori knowledge [15] we inferred a robust and interpretable AD disease module, pinpointing dysregulated genes and processes underlying the eczematous lesions. Our approach allowed us not only to describe the relationships among known patterns of gene expression driving the AD phenotype, but also to identify novel determinants of the disease. Moreover, we investigated the therapeutic potential of drugs targeting genes within the newly identified AD module, via the identification of relevant pharmacophores and prioritized drugs and drug targets to build a compound library to foster tailored and pharmacologically meaningful drug discovery predictions.

Overall, our investigation broadens the understanding of AD by offering unprecedented insights into the genetic, transcriptional, and pharmacological relationships relevant to the pathology. These findings will assist the identification of novel drugs and molecular targets which can aid drug discovery research for the development of more effective therapeutic strategies towards clinical manifestations of AD. The whole analytical pipeline implemented in this study is shown in Fig. 1.

Data collection and preprocessing

All the raw transcriptomics data collected and utilized in this manuscript are publicly available in the Gene Expression Omnibus (GEO) repository. The preprocessed data consist of 12 microarray-derived gene expression datasets of both lesional (337 samples) and non-lesional skin (542 samples) from patients with AD. GEO IDs of the collected datasets are reported in Table S1.

The preprocessing procedure was carried out as described in [16]. Differentially expressed genes for each dataset were identified through the use of the eUTOPIA software [17] by comparing the lesional skin samples with the non-lesional ones. For the analysis, eUTOPIA default parameters were used.

A priori information about AD were retrieved from several public repositories. GWAS hits were collected from the GWAS catalog (https://www.ebi.ac.uk/gwas/). Disease-gene associations and drug targets were retrieved from OpenTargets (https://www.opentargets.org/).

Data scaling

All of the collected microarray datasets were combined for cross-platform normalization. In particular, the pamr R package (version 1.56.1) [18] was used to mean-adjust the combined microarray data based on a batch variable representing the different datasets downloaded from GEO.

Network inference and analysis

Two distinct co-expression networks were inferred by using the gene expression profiles of the lesional and non-lesional skin samples from all the included studies and selecting the genes common to all the platforms. The co-expression networks were inferred through the use of the INfORM algorithm [19]. We set up INfORM in order to build a robust consensus network by using the clr [20], aracne [21] and mrnet [22] algorithms with the following correlation and mutual information measures: Pearson correlation, Kendall correlation, Spearman correlation, empirical mutual information, Miller-Madow asymptotic bias corrected empirical estimator, Schurmann-Grassberger estimate of the entropy of a Dirichlet probability distribution and a shrinkage estimate of the entropy of a Dirichlet probability distribution, as implemented in the minet Bioconductor package [23]. In order to carry out a network community detection we used the walktrap algorithm [24], implemented in INfORM. All computations performed on the inferred networks were carried out through the use of the igraph Bioconductor package [25]. The rendering of co-expression networks was performed by employing the gephi software [26]. In this manuscript, we show a reduced representation of the networks (subgraph) in order to facilitate the visualisation.

AD gene rank construction

The ranked compendium constructed in this study encompasses genes whose association with AD was inferred in two ways: 1) data driven association, resulting from the meta-analysis of the collected transcriptomic datasets, and the analysis of network connectivity, including differential centrality analysis and the analysis of the bridge genes; 2) association deriving from a priori knowledge, including genomic variants associated with AD, annotated disease-gene associations and drug-target associations.

Genes deriving from data-driven association

1. Gene expression meta-analysis: Meta-analysis of the transcriptomic datasets was implemented as described in del Giudice et al. [27] according to a consensus strategy previously reported to be suitable for gene signature retrieval from omics data [28, 29].

Briefly, the pipeline is based on the combination of effect-size, p-value based and rank-product methods. Fisher test is the standard meta-analysis approach to combine p-values. Effect-size strategies are considered the gold standard to assess within- and between-study variation, while accounting for small sample sizes. Finally, the Rank-product method allows combining individual studies according to the within-study gene expression analysis. The intermediate individual ranks were obtained by using the “sumlog” function of the metap R package (Dewey, M. metap: meta-analysis of significance values. CRAN (2022)), the “effect_sizes” function of the esc R package (Lüdecke, D. Esc: Effect Size Computation for Meta Analysis. CRAN (2019)), and the “RP-advance” function of the RankProd R package [30], respectively. The final consensus gene rank was generated through the Borda function of the TopKlists R package [31].

2. Differential centrality analysis: for each of the networks, their node betweenness, closeness and degree centralities were calculated with the Python’s NetworkX package (Python 3.6, NetworkX 2.3). The nodes were ranked according to each of the centrality measures. For each of the networks, their nodes’ median rank based on the rankings of the three centrality measures were calculated. To compare the network of the lesional skin with the non-lesional one, the absolute difference between the median ranks of the two networks was calculated and the genes were ranked accordingly.

3. Bridge gene analysis: The identification of the bridge genes was conducted as reported in [15]. In our previous study, we hypothesized that, by studying the connectivity patterns among known disease genes, it is possible to identify additional associated genes. Therefore, in this study, we identified all the genes connecting pairs of known AD-associated genes within each of the networks (lesional and non-lesional, respectively), and hence acting as a bridge.

Genes deriving from a priori knowledge association

Genes mapped to genomic variants associated with AD were retrieved from the GWAS catalog [32], https://www.ebi.ac.uk/gwas/].
AD disease genes were retrieved from the DisGeNet repository [33].
Drug-target associations, retrieved from the OpenTarget database [34].

The compendium was ultimately built by merging genes from both categories and ranking them on the basis of the total number of evidence of association with AD.

Inference of the AD disease module

The AD disease module was inferred by using the DIAMOnD algorithm [35]. DIAMOnD works under the hypothesis that disease modules are not usually represented by highly cohesive subgraphs, differently from functional modules, which often overlap with topological communities. DIAMOnD identifies disease modules by making use of a set of genes that are already known to be associated with the disease of interest. The algorithm assesses whether a certain gene holds more connections to seed genes than expected, calculating a connectivity p-value. In our study, the genes of the compendium holding at least 3 pieces of evidence of association with AD were used as seed genes for the algorithm in both the lesional and the non lesional network models. Genes that were identified by the DIAMOnD algorithm showing a p-value ≤ 0.05 in each network were considered for analysis. We therefore selected genes that were significant in the lesional network but not in the non-lesional one. The AD disease module was eventually defined by merging such a set of genes with the seed genes.

Skin and immune cell type specific gene enrichment

The enrichment of skin and immune cell type-specific genes in the disease module has been evaluated by utilizing the single-cell RNA-Sequencing derived signatures publicly available at the Human Protein Atlas (https://www.proteinatlas.org/). The skin cell types included in this analysis are the following: adipocytes, eccrine cells, endothelial cells, fibroblasts, hair, keratinocytes, Langerhans cells, macrophages, mast cells, melanocytes, mitotic, outer, plasma cells, sebaceous cells, smooth cells, skin T-cells. The selected pool of immune cell types comprise: T-cells, monocytes, macrophages, granulocytes, B-cells, basal cells, dendritic cells, erythroid cells, Langherans cells, NK-cells, plasma and suprabasal cells.

One tail gene set enrichment analyses (GSEA) were performed through Kolmogorov-Smirnov statistics, as implemented in the stats R package. Overrepresentation tests were performed by using the bc3net CRAN package [36].

AD disease module druggability evaluation

To evaluate the druggability of the AD disease module, we combined drug-target, disease-medication and disease-disease similarity information retrieved from DrugBank [37], STITCH [38], Open Targets [39], Pharos [40], and the GWAS catalog [32]. To integrate disease data from different data sources we annotated their disease names with NCBI MedGen concept IDs as a common vocabulary.

Similarity of AD with other diseases was computed based on genes encompassing genetic variants associated with the disease as annotated in the GWAS catalog. We merged phenotypes representing AD in the GWAS catalog (eczema and atopic dermatitis) by combining their associated genes. Diseases with fewer than 5 genes associated were removed from the analysis. Using the gene lists of 611 diseases, we calculated pairwise disease similarities by computing the Jaccard index, cosine similarity, Sørensen–Dice coefficient, and overlap coefficient, as well as disease distances by the Euclidean and Hamming distance. Similarity metrics were inverted to measures of distance by subtraction from 1. Euclidean and Hamming distances were scaled to a 0–1 range by division by the maximum. We computed the Ipsen-Mikhailov distance among the six distance measures using the netdist function from the R nettools package (https://rdrr.io/cran/nettools/). A consensus distance was computed from the six distance measures by taking the mean in a hierarchical fashion based on hierarchical clustering of Ipsen-Mikhailov distances. Based on the consensus distances, the diseases were ranked from the most similar to the least similar to AD.

Subsequently, we extracted disease-drug-target associations from DrugBank, STITCH, Open Targets and Pharos and kept only the targets i.e., genes belonging to the AD disease module. We prioritised drug candidates based on their therapeutic application for diseases similar to AD.

Virtual screening was performed on an in-house library containing chemical compounds from Drugbank, Pharos and Open Targets.

Building 2D and structures 3D structures for the investigated compounds

The 2D sketcher tool of Maestro molecular modeling suite was used to build the 2D structures from the smiles of the compounds.

The LigPrep module of the Maestro molecular modeling suite was used for ligand preparation. Neutral pH was used in the ligand preparation with OPLS4 force field and determining chiralities from 3D structure options were selected [41, 42].

Identifying pharmacophore hypothesis and virtual screening

The Phase tool of Maestro was used to develop a common pharmacophore hypothesis for the active ligands and to use it for pharmacophore-based virtual screening. The virtual screening was performed on an in-house library containing chemical compounds from Drugbank, Pharos and Open Targets to identify compounds that share the selected pharmacophore and features [43].

ADME properties calculation

To identify the ADME (Absorption, Distribution, Metabolism, Excretion) properties for the ligands, the QikProp tool in Maestro was used. (QikProp, Schrödinger, LLC, New York, NY, 2021)

An AD-specific gene module represents relevant functions within the lesional network model

Complex diseases, including atopic dermatitis, often result from the impairment of complex cascades of interactions between molecular effectors. In this context, the construction of disease network models has been shown to be a valid method in order to disentangle such complex interactions, and identify key players in the pathophenotype of diseases [35]. Here, we inferred two transcriptome-wide gene co-expression networks modelling both the lesional and the non-lesional skin molecular buildup of AD patients, respectively. Since the networks were inferred from all the genes common to all the collected transcriptomics datasets, both networks were composed of 17,903 nodes, while the number of edges was 1,649,569 and 1,924,001 for the lesional and the non-lesional, respectively.

Next, we exploited such networks to define a disease module for AD in order to identify a core set of genes that might play a role in the disease and might function as therapeutic targets. To do so, we compiled a ranked compendium of genes whose association with AD is based on 1) a priori knowledge about their involvement in AD and 2) data-driven analyses of the transcriptomics datasets and co-expression networks (see Methods). Subsequently, we selected 1,234 genes (here named “seed genes”) having at least three supporting evidence of their association with AD as reference genes for the identification of the AD disease module. We carried out this analysis on both the lesional and the non-lesional network by utilising the same set of seed genes. We therefore identified 1,218 genes that are significantly more connected with the seed genes in the lesional network but not in the non-lesional one. In this way, we defined the atopic dermatitis module, encompassing 2,452 genes.

We characterised the functional processes represented by such genes by performing an overrepresentation analysis by Fisher’s exact test (Fig. 2).

Interestingly, 181 genes of the disease module are annotated in the “Signaling by Interleukins” pathway, which is the most overrepresented in our analysis (adj. p-value = 2.78e-37). Similarly, other pathways related to interleukin signalling and extracellular matrix remodelling were found to be overrepresented, including “Interleukin-4 and Interleukin-13 signaling” (adj. p-value = 3.97e-29), “Degradation of the extracellular matrix” (adj. p-value = 1.45e-07), “Interleukin-1 family signaling” (adj. p-value = 1.33e-06), “TNFR2 non-canonical NF-kB pathway” (adj. p-value = 1.3e-08), “Interleukin-10 signaling” (adj. p-value = 1.78e-21), “collagen formation” (adj. p-value = 3.18e-08), “collagen degradation” (adj. p-value = 9.8e-08), among others.

The manifestation of AD lesions arises from a dysregulated cutaneous immune response and an impairment of the skin barrier, which is normally composed by a wide range of cell types contributing to its homeostasis. AD lesional skin is characterised by a chronic inflammatory milieu due to several cellular dysfunctions including a disturbed epidermal barrier integrity, increased immune cell infiltration, abnormal secretion from sebaceous cells and eccrine gland cells. Therefore, we further characterised the genes within the AD module, by considering the known patterns of expression in multiple skin and immune cell types. Our analysis supports the functional annotation reported above and identifies cell types that underlie relevant functions. We, then, exploited publicly available gene expression signatures from several skin cell types from the Human Protein Atlas database by performing a Gene Set Enrichment Analysis (GSEA) to assess the enrichment of cell type-specific genes within the AD disease module (Fig. 3). Our analysis highlighted a significant enrichment of genes specifically expressed in keratinocytes (different stages of differentiation, adj. p-value = 1.64e-12), eccrine gland cells (adj. p-value = 3.29e-15), skin fibroblasts (adj. p-value = 1.64e-12), skin endothelial cells (adj. p-value = 1.91e-10), skin adipocytes (adj. p-value = 4.30e-07), skin mitotic cells (adj. p-value = 9.56e-03), and granular keratinocytes (adj. p-value = 2.72e-02). Similarly, since the dysregulation of immune cells and their crosstalk plays a crucial role in the development and maintenance of atopic dermatitis, we also assessed the presence of cell signatures related to the immune compartment within the AD disease module. As a result, we obtained a significant enrichment for the following immune cells: macrophages (adj. p-value = 1.04e-21), monocytes (adj. p-value = 3.2e-25), T-cells (adj. p-value = 1.04e-21), granulocytes (adj. p-value = 5.2e-19), Langerhans cells (adj. p-value = 7.89e-13), dendritic cells (adj. p-value = 6.71e-14), B-cells (adj. p-value = 2.23e-14), NK-cells (adj. p-value = 2.69e-09). The entire list of skin and immune cell types included in this analysis is described in the Methods section.

Co-expression patterns underlying AD lesional skin drives biomarker discovery

We here propose that investigating connectivity patterns within the disease module could serve as an alternative and effective method for biomarker prioritisation. Traditionally, biomarkers have been understood as disrupted expression profiles of single genes occurring in skin lesions. However, network medicine suggests that specific connectivity patterns involving interactions between two or more genes could be utilized as biomarkers for these lesions.

Therefore, to fully exploit the informative power of the AD disease module, we characterised the edges (representing co-expression relationships between two genes) connecting genes within the module.

We evaluated the relevance of each gene within the AD module by collecting multiple data-driven aspects (transcriptome meta-analysis, differential centrality, and bridge gene analysis) as well as prior-knowledge (GWAS, drug targets, disease genes) (Fig. 1A). We used data-driven and prior evidence to rank the edges connecting genes within the module (Fig. 4).

Based on these criteria, TNFRSF8–IL13, IL22–IL13 and IL4R–IL22 are the top 3 ranked edges suggesting the high relevance of these genes and their co-expression for AD. The top 100 edges connecting relevant genes are reported in Supplementary file 2. Next, we divided the edges in two sets: the first group includes edges connecting pairs of seed genes (used to infer the AD disease module), while in the second, edges for which only one of the genes is a seed gene were included. In this way, we characterised functional processes involving only seed genes, and non-seed genes within the disease module via direct connection to seed genes. The set of edges connecting seed genes with each other enriched several signalling by interleukins pathways, extracellular matrix organisation and degradation, and cell cycle related pathways, such as mitotic G1/S transition and several pathways related with collagen deposition and modification, such as Collagen formation, collagen degradation, collagen biosynthesis and modifying enzymes, collagen chain trimerization. On the other hand, the set of edges connecting seed genes with other genes enriched Processing of capped intron-containing pre-mRNA, mRNA splicing pathways and Rho-GTPase cycle. This analysis allows us to enrich the information contained in the module by interactional properties of data-driven and prior knowledge-driven genes contained into it.

Definition of a custom library of compounds for drug discovery in AD

Subsequently, we constructed a custom library of compounds to foster drug discovery for the definition of efficient therapeutic compounds. We extracted drug-diseases and drug-target relationships from public repositories, including DrugBank [37], STITCH [38], Open Targets [39], Pharos [40]. Notably, we retrieved 285 drugs targeting 635 genes of the AD disease module. To provide further relevance to the drugs within the module, we prioritised compounds that are used to treat diseases that share a similar genetic background with AD. Therefore, we computed a disease-disease similarity and prioritised diseases showing the highest similarity with AD based on genes with genomic variants. We extracted drug-disease relationships and obtained a list of disease-drug-target relationships (Fig. 5A) by merging them with the drug-target relationships. We ranked the list based on the AD-similarity of the disease in the disease-drug-target triplets including AD itself (i.e., drug-target entities linked to AD would always have the highest rank). Afterwards, we filtered such a rank on the drug target genes belonging to the AD disease module (Fig. 5B).

The 285 drugs included in the rank are employed in the treatment of 66 diseases. The most similar disease to AD was autoimmune hepatitis, followed by neuromyelitis optica, asthma, nephrotic syndrome and allergic rhinitis. Drugs used to treat AD such as prednisolone, cyclosporine, tacrolimus, and dupilumab ranked at the top. Moreover, drugs used to treat diseases similar to AD ranked at lower positions. Abatacept and mycophenolate mofetil are used to treat autoimmune hepatitis, the disease with the highest similarity to AD based on genetic susceptibility (Supplementary file 3). We here hypothesised that potentially effective drugs for AD share relevant structural properties, namely pharmacophores.

Therefore, we extrapolated relevant pharmacophores of such drugs in order to identify the minimal active substructure of such compounds (Supplementary file 1 - Table S2).

In order to extrapolate more accurate pharmacophores, we selected 6 features responsible for desired effects.

Among all the drugs included in the disease-drug-target rank computed in the previous step, 19 share the highest ranked pharmacophore hypothesis including prednisolone, tacrolimus, and pimecrolimus, which are approved for treatment of AD. We further tested the hypothesis that additional drugs, not initially included in our list, could be proposed for AD treatment if relevant pharmacophores are included in their molecular structure. To this end, we performed a virtual screening, which consists in screening large libraries of compounds in search of specific structural properties of interest.

As a result, we identified 1,505 compounds that share the selected pharmacophore and other relevant structural features (Supplementary file 4). Corticosteroids are the most represented class of drugs at the top of the rank, such as fluocinonide and halopredone acetate. NCX1022, a NO-releasing derivative of hydrocortisone, placed at the 3rd position of the rank, has been investigated in a murine model of AD and showed anti-inflammatory effects [44]. Given the relevance of this investigational compound in the library, we assessed whether such a compound shares molecular targets with drugs that are approved for inflammation management. Interestingly, our prediction shows that NCX1022 shares 21 molecular targets with abrocitinib and upadacitinib.

Desoximetasone, which is a glucocorticoid used for the treatment of inflammatory conditions, skin allergies and dermatoses, also ranked 12th among the drugs.

We finally computed the Adsorption, Distribution, Metabolism and Excretion (ADME) values of the identified library of compounds (Supplementary file 5, Supplementary file 1 - Figure S1).

The ADME analysis reveals that NCX1022 shows higher metabolism rate compared with the other drugs. Moreover, NCX1022 shows medium oral absorption, while the other drugs from the top 10 showed high oral absorption. These findings show that the top 10 ranked drugs have medium to high ADME values, which is a predictive factor for the bioavailability and for the maximum therapeutic efficacy.

Network models provide a powerful framework for integrating heterogeneous datasets and are at the forefront of multi-omics data analysis and interpretation, enabling the discovery of hidden patterns, relationships, and insights that would otherwise be difficult to detect using traditional approaches. By leveraging the interconnectedness and interdependencies within a co-expression network, such as nodes representing genes and edges representing relationships between them, network models offer a holistic view of complex systems, such as human diseases.

In this study, we provide a network-based modelling approach integrating the biggest number of cutaneous transcriptomic datasets of AD, to date. We computed gene co-expression networks of lesional and non-lesional skin, both consisting of 17,903 genes. Notably, the lesional network contained substantially fewer edges (1,649,569) than the non-lesional network (1,924,001) highlighting disrupted patterns of co-expression in AD skin lesions and indicating dysfunctional gene regulatory mechanisms.

In recent years, the advancement of analysis strategies of disease-related networks, led to the identification of so-called “disease modules”. Disease modules are subsets of genes highly interconnected within a disease network from which we can gain insights into the pathophysiology of the disease under study. Genes belonging to a disease module often work together in specific pathways or functional units, shedding light on how impaired molecular mechanisms contribute to the disease development or progression [45, 46]. By combining data-driven evidence with prior knowledge, we here identified a robust disease module of the AD skin lesion encompassing 2,452 genes.

As this community of genes related to the atopic lesion was derived by several features including gene deregulation, gene-gene co-expression relationships, and genetic variation associated with the disease, the computation of the disease module overperforms the computation of differentially expressed genes as usually identified and used for downstream analysis in single-dataset transcriptomic studies [35].

While gene sets derived from differential expression analysis alone usually reveal a pleiotropy of pathways related to inflammation and innate immune responses, the disease module identified in this study showed the strongest enrichment for pathways targeted by new state-of-the-art treatment solutions for AD such as biologics targeting IL4-IL13 signalling (dupilumab and tralokinumab), and drugs targeting janus kinases 1,2 and 3 (barticitinib, upadacitinib and abrocitinib). However, the lack of efficacy in some patients highlights the need for other effective molecules, and underlines the importance of developing targeted medicine tailored to patients characteristics. Further enriched pathways reveal the relevance of TNF, NF-kB, and chemokine signalling as well as a strong signature of epidermal barrier dysfunction. Here we also identified a group of genes that is overrepresented in a high number of known dysregulated pathways in the AD lesion. Such genes, including MMP9, MMP3 and COL1A2 among others, might represent key players in the molecular impairment underlying the aberrant immune response and the dysfunctional barrier as characteristic for AD. In fact, keratinocytes express MMP9 as an important factor in the maintenance of the epithelial barrier function [47].

The analysis of the disease module revealed transcriptomic signatures from abundant cell types in the skin such as keratinocytes, eccrine gland cells, fibroblasts, and endothelial cells as well as immune cells such as lymphocytes and innate immune cells indicating that a complex dysregulated multi-cellular system sustains the phenotype under study. Understanding such complex multi-cellular interactions is crucial for advancing disease modelling and therapeutic interventions. In this perspective, it emerges the need for more sophisticated disease modelling procedures that incorporate the interactions and crosstalk among different cell types, which can lead to more precise and effective treatments.

Importantly, the characterisation of the gene connectivity within the disease module leveraged fundamental insights in the gene deregulation underlying the skin lesion. In fact, while relationships between known AD-associated genes are overrepresented in pathways that are typically disrupted in immune-mediated dermatological diseases, such as interleukin signalling, extracellular matrix organisation, collagen formation and degradation, connection between AD-associated and non-associated genes underlined mRNA processing indicating that the gene regulation machinery is tightly intertwined with known AD biomarkers such as interleukins and chemokines to sustain the AD pathological phenotype.

In this study, we investigated the interactional properties of genes belonging to the disease module, so as to give further insights into its functional characterisation. We demonstrated that, ranking the edges based on data-driven and prior knowledge-driven evidence of association with AD, is a valid approach to fully exploit network medicine principles to prioritise putative biomarkers of disease.

In fact, we hypothesise that exploiting properties of connectivity patterns within the disease network model can be an alternative and efficient tool for biomarker identification. The paradigm is that, while traditionally a biomarker is intended as a disrupted expression profile of a single gene specifically occurring in the skin lesion, network medicine suggests that certain patterns of connectivity involving two or more interacting genes could be exploited as biomarkers of lesions.

While the use of precompiled drug libraries is the state of the art in drug discovery, the principles by which the drugs are included are usually not evaluated. However, a rational inclusion process of drugs in libraries for drug discovery would enable a more effective process with a sensible improvement of the predictions, as we recently demonstrated in a case study on COVID-19 drug repositioning [48].

Since disease modules provide a formal perspective over the complexity of disease mechanisms, network pharmacology approaches can take advantage of this knowledge by considering the dependencies among genes belonging to the disease module. Given the high reliability of disease modules identified by integrating data-driven and a priori knowledge about the disease under study, compound libraries can be, in turn, designed in order to target multiple components belonging to the disease module, aiming to optimize the identification of compounds putatively effective for the disease. The definition of tailored compound libraries based on disease modules can leverage the identification of compounds that can have a broader impact on the disease phenotype and offer potential synergistic effects. Although the employment of network pharmacology approaches to identify candidate drugs and targets is well established, the possibility of using such principles to suggest desired structural properties to be considered in the discovery of drug candidates is less explored. Indeed, while current network pharmacology approaches are centered on molecular biology aspects, here we establish the paradigm that intimately integrates molecular biology with computational chemistry in order to construct tailored compound libraries.

We here demonstrate how the disease module can assist in building a custom compound library for drug discovery in AD. To do so, we prioritised common pharmacophores of drugs employed in the treatment of AD or other similar phenotypes whose target genes belong to the disease module. We identified a substructure shared among the investigated drugs. Afterwards, we performed a virtual screening on multiple high-dimensional repositories and defined a custom drug library to tailor drug discovery efforts towards the treatment of atopic dermatitis. This study is subject to certain limitations due to the limited amount of clinical data available along with the transcriptional profiles obtained from public repositories. The absence of detailed clinical information prevents us from establishing meaningful associations between the disease module and relevant clinical parameters, such as disease severity and elapsed time since initial diagnosis. Furthermore, the lack of information on whether patients underwent active pharmacological therapy precludes us from discerning any potential influence of ongoing or terminated pharmacological treatments on the druggability profile of the disease module, which resulted, eventually, in the tailored compound library.

In conclusion, our investigation significantly enhances our understanding of the molecular basis of AD by providing novel and unprecedented insights into the mechanistic and pharmacological relationships driving the AD phenotype. These results hold promise in identifying novel molecular targets and drug candidates, thereby propelling drug discovery research towards the development of effective therapeutic interventions for AD.

Atopic Dermatitis

ADME

Absorption, Distribution, Metabolism, Excretion

COVID-19

Coronavirus Disease 2019

GEO

Gene Expression Omnibus

GSEA

Gene Set Enrichment Analysis

GWAS

Genome Wide Association Study

IL4

Interleukin 4

IL4R

Interleukin 4 Receptor

IL13

Interleukin 13

IL22

Interleukin 22

NF-kB

Nuclear Factor Kappa B

Nitric Oxide

OPLS4

Optimised potential for liquid simulations

TNF

Tumor Necrosis Factor

TNFRSF8

Tumor Necrosis Factor Receptor Superfamily 8

Ethics approval and consent to participate

Not applicable

Consent for publication

Not applicable

Availability of data and material

Not applicable

Competing interests

The authors declare no conflict of interest.

Funding

This work was supported by the EU Innovative Medicines Initiative 2 (IMI2) Biomap Project [grant agreement number 821511], Academy of Finland [grant agreement number 322761], and the European Research Council (ERC) programme, Consolidator project “ARCHIMEDES” [grant agreement number 101043848].

Acknowledgements

AF and was supported by Tampere Institute for Advanced Study (IAS).

Author information

Authors and Affiliations

Finnish Hub for Development and Validation of Integrated Approaches (FHAIVE), Tampere University, Tampere, Finland
Antonio Federico, Lena Möbus, Zeyad Al-Abdulraheem, Alisa Pavel, Giusy del Giudice & Dario Greco
Tampere Institute for Advanced Studies, Tampere University, Tampere, Finland
Antonio Federico
Institute of Biomedicine, University of Eastern Finland, Kuopio, Finland
Vittorio Fortino
Faculty of Medicine, Human Microbiome Research Program, University of Helsinki, Helsinki, Finland Harri Alenius & Nanna Fyhrquist
Institute of Environmental Medicine (IMM), Karolinska Institutet, Stockholm, Sweden Harri Alenius & Nanna Fyhrquist
Division of Pharmaceutical Biosciences, Faculty of Pharmacy, University of Helsinki, 00100 Helsinki, Finland Antonio Federico, Giusy del Giudice & Dario Greco
Institute of Biotechnology, University of Helsinki, Helsinki, Finland Dario Greco

Contributions

AF and DG designed the study. AF and LM performed the analyses. AP, ZAA and VF contributed to the analyses. AF, LM, AP, ZAA and VF participated in data generation and curation. AF, LM, GdG, NF, HA and DG interpreted the results. AF, LM, GdG, and DG wrote the manuscript. DG supervised the study. All authors read and approved the final manuscript.

Langan SM, Irvine AD, Weidinger S. Atopic dermatitis. Lancet. 2020 Aug 1;396(10247):345–60.
Renert-Yuval Y, Thyssen JP, Bissonnette R, Bieber T, Kabashima K, Hijnen D, et al. Biomarkers in atopic dermatitis-a review on behalf of the International Eczema Council. J Allergy Clin Immunol. 2021 Apr;147(4):1174-1190.e1.
Bakker D, de Bruin-Weller M, Drylewicz J, van Wijk F, Thijs J. Biomarkers in atopic dermatitis. J Allergy Clin Immunol. 2023 May;151(5):1163–8.
Tokura Y, Hayano S. Subtypes of atopic dermatitis: From phenotype to endotype. Allergol Int. 2022 Jan;71(1):14–24.
Czarnowicki T, He H, Krueger JG, Guttman-Yassky E. Atopic dermatitis endotypes and implications for targeted therapeutics. J Allergy Clin Immunol. 2019 Jan;143(1):1–11.
Tsoi LC, Rodriguez E, Degenhardt F, Baurecht H, Wehkamp U, Volks N, et al. Atopic Dermatitis Is an IL-13-Dominant Disease with Greater Molecular Heterogeneity Compared to Psoriasis. J Invest Dermatol. 2019 Jul;139(7):1480–9.
Paternoster L, Standl M, Waage J, Baurecht H, Hotze M, Strachan DP, et al. Multi-ancestry genome-wide association study of 21,000 cases and 95,000 controls identifies new risk loci for atopic dermatitis. Nat Genet. 2015 Dec;47(12):1449–56.
Ewald DA, Malajian D, Krueger JG, Workman CT, Wang T, Tian S, et al. Meta-analysis derived atopic dermatitis (MADAD) transcriptome defines a robust AD signature highlighting the involvement of atherosclerosis and lipid metabolism pathways. BMC Med Genomics. 2015 Oct 12;8:60.
Möbus L, Rodriguez E, Harder I, Stölzl D, Boraczynski N, Gerdes S, et al. Atopic dermatitis displays stable and dynamic skin transcriptome signatures. J Allergy Clin Immunol. 2021 Jan;147(1):213–23.
Barabási A-L, Oltvai ZN. Network biology: understanding the cell’s functional organization. Nat Rev Genet. 2004 Feb;5(2):101–13.
Barabási A-L, Gulbahce N, Loscalzo J. Network medicine: a network-based approach to human disease. Nat Rev Genet. 2011 Jan;12(1):56–68.
Facheris P, Jeffery J, Del Duca E, Guttman-Yassky E. The translational revolution in atopic dermatitis: the paradigm shift from pathogenesis to treatment. Cell Mol Immunol. 2023 May;20(5):448–74.
Chovatiya R, Paller AS. JAK inhibitors in the treatment of atopic dermatitis. J Allergy Clin Immunol. 2021 Oct;148(4):927–40.
Skou ST, Mair FS, Fortin M, Guthrie B, Nunes BP, Miranda JJ, et al. Multimorbidity. Nat Rev Dis Primers. 2022 Jul 14;8(1):48.
Federico A, Pavel A, Möbus L, McKean D, Del Giudice G, Fortino V, et al. The integration of large-scale public data and network analysis uncovers molecular characteristics of psoriasis. Hum Genomics. 2022 Nov 28;16(1):62.
Federico A, Hautanen V, Christian N, Kremer A, Serra A, Greco D. Manually curated and harmonised transcriptomics datasets of psoriasis and atopic dermatitis patients. Sci Data. 2020 Oct 13;7(1):343.
Marwah VS, Scala G, Kinaret PAS, Serra A, Alenius H, Fortino V, et al. eUTOPIA: solUTion for Omics data PreprocessIng and Analysis. Source Code Biol Med. 2019 Jan 29;14:1.
T. Hastie, R. Tibshirani, Balasubramanian Narasimhan and Gil Chu. Pam: Prediction Analysis for Microarrays. R package version 1.56.1. https://CRAN.R-project.org/package=pamr. CRAN R-project; 2019.
Marwah VS, Kinaret PAS, Serra A, Scala G, Lauerma A, Fortino V, et al. Inform: inference of network response modules. Bioinformatics. 2018 Jun 15;34(12):2136–8.
Faith JJ, Hayete B, Thaden JT, Mogno I, Wierzbowski J, Cottarel G, et al. Large-scale mapping and validation of Escherichia coli transcriptional regulation from a compendium of expression profiles. PLoS Biol. 2007 Jan;5(1):e8.
Margolin AA, Nemenman I, Basso K, Wiggins C, Stolovitzky G, Dalla Favera R, et al. ARACNE: an algorithm for the reconstruction of gene regulatory networks in a mammalian cellular context. BMC Bioinformatics. 2006 Mar 20;7 Suppl 1(Suppl 1):S7.
Meyer PE, Kontos K, Lafitte F, Bontempi G. Information-theoretic inference of large transcriptional regulatory networks. EURASIP J Bioinform Syst Biol. 2007;79879.
Meyer PE, Lafitte F, Bontempi G. minet: A R/Bioconductor package for inferring large transcriptional networks using mutual information. BMC Bioinformatics. 2008 Oct 29;9:461.
Pons P, Latapy M. Computing communities in large networks using random walks. J Graph Algorithms Appl. 2006;10(2):191–218.
Csardi G, Nepusz T. The igraph software package for complex network research. InterJournalComplex Systems,. 2006;1695.
Bastian M, Heymann S, Jacomy M. Gephi: An Open Source Software for Exploring and Manipulating Networks. Unpublished. 2009;
Del Giudice G, Serra A, Saarimäki LA, Kotsis K, Rouse I, Colibaba SA, et al. An ancestral molecular response to nanomaterial particulates. Nat Nanotechnol. 2023 Aug;18(8):957–66.
Kan M, Shumyatcher M, Diwadkar A, Soliman G, Himes BE. Integration of Transcriptomic Data Identifies Global and Cell-Specific Asthma-Related Gene Expression Signatures. AMIA Annu Symp Proc. 2018 Dec 5;2018:1338–47.
Hong F, Breitling R. A comparison of meta-analysis methods for detecting differentially expressed genes in microarray experiments. Bioinformatics. 2008 Feb 1;24(3):374–82.
Del Carratore F, Jankevics A, Eisinga R, Heskes T, Hong F, Breitling R. RankProd 2.0: a refactored bioconductor package for detecting differentially expressed features in molecular profiling datasets. Bioinformatics. 2017 Sep 1;33(17):2774–5.
Lin S. Space oriented rank-based data integration. Stat Appl Genet Mol Biol. 2010 Apr 9;9:Article20.
Sollis E, Mosaku A, Abid A, Buniello A, Cerezo M, Gil L, et al. The NHGRI-EBI GWAS Catalog: knowledgebase and deposition resource. Nucleic Acids Res. 2023 Jan 6;51(D1):D977–85.
Piñero J, Ramírez-Anguita JM, Saüch-Pitarch J, Ronzano F, Centeno E, Sanz F, et al. The DisGeNET knowledge platform for disease genomics: 2019 update. Nucleic Acids Res. 2020 Jan 8;48(D1):D845–55.
Ochoa D, Hercules A, Carmona M, Suveges D, Gonzalez-Uriarte A, Malangone C, et al. Open Targets Platform: supporting systematic drug-target identification and prioritisation. Nucleic Acids Res. 2021 Jan 8;49(D1):D1302–10.
Ghiassian SD, Menche J, Barabási A-L. A DIseAse MOdule Detection (DIAMOnD) algorithm derived from a systematic analysis of connectivity patterns of disease proteins in the human interactome. PLoS Comput Biol. 2015 Apr 8;11(4):e1004120.
de Matos Simoes R, Emmert-Streib F. Bagging statistical network inference from large-scale gene expression data. PLoS ONE. 2012 Mar 30;7(3):e33624.
Wishart DS, Feunang YD, Guo AC, Lo EJ, Marcu A, Grant JR, et al. DrugBank 5.0: a major update to the DrugBank database for 2018. Nucleic Acids Res. 2018 Jan 4;46(D1):D1074–82.
Kuhn M, von Mering C, Campillos M, Jensen LJ, Bork P. STITCH: interaction networks of chemicals and proteins. Nucleic Acids Res. 2008 Jan;36(Database issue):D684-8.
Ochoa D, Hercules A, Carmona M, Suveges D, Baker J, Malangone C, et al. The next-generation Open Targets Platform: reimagined, redesigned, rebuilt. Nucleic Acids Res. 2023 Jan 6;51(D1):D1353–9.
Kelleher KJ, Sheils TK, Mathias SL, Yang JJ, Metzger VT, Siramshetty VB, et al. Pharos 2023: an integrated resource for the understudied human proteome. Nucleic Acids Res. 2023 Jan 6;51(D1):D1405–16.
Shelley JC, Cholleti A, Frye LL, Greenwood JR, Timlin MR, Uchimaya M. Epik: a software program for pK( a ) prediction and protonation state generation for drug-like molecules. J Comput Aided Mol Des. 2007 Dec;21(12):681–91.
Banks JL, Beard HS, Cao Y, Cho AE, Damm W, Farid R, et al. Integrated modeling program, applied chemical theory (IMPACT). J Comput Chem. 2005 Dec;26(16):1752–80.
Dixon SL, Smondyrev AM, Rao SN. PHASE: a novel approach to pharmacophore modeling and 3D database searching. Chem Biol Drug Des. 2006 May;67(5):370–2.
Hyun E, Bolla M, Steinhoff M, Wallace JL, Soldato PD, Vergnolle N. Anti-inflammatory effects of nitric oxide-releasing hydrocortisone NCX 1022, in a murine model of contact dermatitis. Br J Pharmacol. 2004 Nov;143(5):618–25.
Hartwell LH, Hopfield JJ, Leibler S, Murray AW. From molecular to modular cell biology. Nature. 1999 Dec 2;402(6761 Suppl):C47-52.
Kitsak M, Sharma A, Menche J, Guney E, Ghiassian SD, Loscalzo J, et al. Tissue specificity of human disease module. Sci Rep. 2016 Oct 17;6:35241.
Jeong S, Ledee DR, Gordon GM, Itakura T, Patel N, Martin A, et al. Interaction of clusterin and matrix metalloproteinase-9 and its implication for epithelial homeostasis and inflammation. Am J Pathol. 2012 May;180(5):2028–39.
Serra A, Fratello M, Federico A, Ojha R, Provenzani R, Tasnadi E, et al. Computationally prioritized drugs inhibit SARS-CoV-2 infection and syncytia formation. Brief Bioinformatics. 2022 Jan 17;23(1).

Supplementaryfile1.pdf
Supplementary file 1
Supplementaryfile2.xlsx
Supplementary file 2
Supplementaryfile3.xlsx
Supplementary file 3
Supplementaryfile4.xlsx
Supplementary file 4
Supplementaryfile5.xlsx
Supplementary file 5

Download PDF

Journal Publication

published 16 Jan, 2024

Read the published version in Journal of Translational Medicine →

Version 1

posted

You are reading this latest preprint version

Integrative network analysis suggests prioritised drugs for atopic dermatitis

Status:

Journal Publication

Version 1

Abstract

Background

Methods

Results

Conclusions

Figures

Background

Methods

Data collection and preprocessing

Data scaling

Network inference and analysis

AD gene rank construction

Inference of the AD disease module

Skin and immune cell type specific gene enrichment

AD disease module druggability evaluation

Building 2D and structures 3D structures for the investigated compounds

Identifying pharmacophore hypothesis and virtual screening

ADME properties calculation

Results

An AD-specific gene module represents relevant functions within the lesional network model

Co-expression patterns underlying AD lesional skin drives biomarker discovery

Definition of a custom library of compounds for drug discovery in AD

Discussion

Conclusions

Abbreviations

Declarations

References

Supplementary Files

Status:

Journal Publication

Version 1