DermaGene and VitmiRS: a comprehensive systems analysis of genetic dermatological disorders

Abstract

Background

The interaction of genetic variants and their distribution in the genome is firmly believed to contribute to genetic dermatological disorders. For the convenience of clinicians and researchers to explore such genetic variants, we identified and validated the genetic association of human SNPs with dermatological disorders through manual curation and computational analysis.

Methods

Multiple online resources were investigated for creating a comprehensive list of dermatological disorders followed by identification and manual curation of data description of the SNP–disease relationship. The process for database creation involved extensive review of published literature and relevant articles that probed association of miRNAs and SNPs with vitiligo. Furthermore, computational analysis for deleterious SNPs and polypharmacological studies was performed with an effort to provide potential novel drug candidates for optimal therapeutic interventions.

Results

We established a dedicated database on dermatological disorders, DermaGene, which proffers a user-friendly interface to enable systematic querying and analysis of SNPs associated with dermatological disorders. Based on our disease network analysis, we further extended our present work to construct another comprehensive database, VitmiRS (Vitiligo associated miRNAs and SNPs), which furnishes detailed information on each miRNA and SNP association that renders susceptibility to vitiligo. Also, our systems analysis approach unraveled potential molecular determinants that may pose as novel drug candidates and can be targeted for efficient therapeutic approaches in vitiligo treatment.

Conclusions

Both the databases are freely accessible and will serve as a significant resource providing insights into previously undiscovered disease–gene and disease–SNP relationships. Our analysis unveiled significant findings that we believe may furnish a comprehensive understanding of the biological mechanisms that mediate vitiligo disease pathogenesis, thereby, driving the way towards better therapeutic interventions for disease management.

Background

The importance of systems and personalized medicine in the field of healthcare has grown vastly, especially in the last decade. Recent research has growingly demonstrated that many of the seemingly disparate diseases have a common molecular mechanism and strong association among them. This signifies the deeper interplay of genes, in contrast to one gene-one disease commonality well indicating towards the possibility of the related diseases to occur together in an individual (Yu et al. 2015). The exhibition of comorbidity association suggests that the occurrence of one disease will increase the likelihood of the other thereby contributing to compounded healthcare costs and multiplex clinical management. Also, a higher degree of interaction between many diseases indicates a greater prevalence and mortality rate associated with it. Due to the contribution of disease comorbidities as a burden on societies at large, it becomes crucial for risk stratification and charting treatment course. One of the most powerful approaches in the present times, the Genome-Wide Association Studies (GWAS), serves as the primary pipeline for identifying disease variants and researching association of diseases and their phenotypes (Lewis et al. 2011). The complete sequencing of several genomes has uncovered that the primary provenance of the complexity of human disorders is in parallel with the genetic variations attributable to single nucleotide polymorphisms rather than barely protein–protein interactions. The interaction of genetic variants and the distribution of SNPs in the genome are firmly believed to contribute to the genetic risk of phenotypic variation in human complex diseases (Cordell 2009; Irvine and McLean 2003; Shen et al. 2016). Thus, integrative network analysis utilizing multi-omics data has become increasingly widespread because each approach has intrinsic caveats. A systems-level understanding of the complexities of genetic variations is expected to further enhance disease diagnosis, prognosis, and treatment with improved patient risk stratification and planning of treatments (Park et al. 2009).

The advent of SNP genotyping of dermatological disorders has revolutionized our understanding of the genetics of such disorders (Maruthappu et al. 2014). Relevant information from several published literature concerning the disease–SNP association in dermatological disorders are interspersed which impels the need for a repository that stores SNPs having prime roles in the initiation of dermatological disorders. Such a resource would be valuable for clinicians and researchers making it convenient for them to study the regulatory mechanisms of SNPs that underpin disease pathogenesis and translate their discoveries into clinical practice. We, therefore, aimed to fulfill this void concerning a dedicated genetic dermatological disorder database. DermaGene is the first attempt at providing information pertaining to human genetic dermatological disorders establishing itself as the primary pipeline for identifying disease variants and researching association of dermatological diseases and their phenotypes. The database as of now stores 114 unique diseases, as described by 244 different genes along with 871 associated SNPs and have been made freely accessible online. The present work not only provides information on experimentally validated dermatological disorders associated SNPs curated from a large scale of literature but also offers computationally predicted insights into previously undiscovered disease–gene and disease–SNP relationships.

Our disease association network analysis highlighted vitiligo to exhibit maximum interactions with different genes from our database. Vitiligo is a chronic, acquired depigmentation disorder of the skin resulting in an episodically progressive loss of functional melanocytes causing pigment dilution in the affected areas of the skin (Picardo et al. 2015). Affecting 0.5–1% of the world population (Boisseau-Garsaud et al. 2000; tHowitz et al. 1977) and with a prevalence rate of 0.5–2.5% in India (Handa and Kaur 1999), vitiligo can develop at any age irrespective of the type of skin, gender, race, or geographical location underlying both genetic and non-genetic factors in a complex interactive manner. We, therefore, investigated the role of specific miRNAs and SNPs reported to be associated with vitiligo in the human genome to identify their plausible effect on vitiligo susceptibility using computational platforms. The cataloging of susceptible miRNAs and SNPs is essential for narrowing down the plausible concomitant genetic determinants of vitiligo. Protein–protein interactions are virtually intrinsic for every cellular and regulatory process, and a damaging alteration in such interactions have been deduced to cause and sometimes even accelerate human diseases. The regulation or impediment of a known detrimental protein–protein interaction delineates a principal target for drug discovery. Hence, a systems biology approach was implemented that unveiled significant interconnections and revealed intricate patterns of disease association. Such a network analysis is helpful in studying the gene expressions and analyzing a large set of disease-associated proteins (Barabasi and Oltvai 2004). Furthermore, we prioritized a few proteins in our protein–protein interaction network as pertinent hub proteins which may be targeted for treating vitiligo. We performed computational analysis for harmful SNPs and scrutinized drug–target and drug–similarity interactions with an effort to provide potential novel drug candidates for optimal therapeutic interventions. Subsequently, we manually curated all miRNAs and SNPs and established a vitiligo specific database, VitmiRS (Vitiligo associated miRNAs and SNPs), to provide a comprehensive, user-friendly interface to detailed information on each miRNA and SNP association that renders susceptibility to vitiligo supported with pertinent data. VitmiRS currently houses 41 miRNAs and 134 associated SNPs for 84 genes which have been made online allowing easy retrieval of data. Till date, no systematic efforts have been made for the compilation of such data in a single platform that is specifically designed for dermatological disorders thereby contributing to the uniqueness of both the databases. A comprehensive view of the databases is shown in Fig. 1. We hope they serve as a significant resource that would benefit the scientific community. The present polypharmacogenomic network integration approach in vitiligo disease module is novel to the best of our knowledge that attempts to furnish a comprehensive understanding of the biological mechanisms that mediate vitiligo disease pathogenesis and offering novel interventional drug targets for efficient therapeutic interventions for genetically susceptible individuals. The graphical representation of the workflow is shown in Fig. 2.

Fig. 1

A comprehensive view of DermaGene and VitmiRS database

Fig. 2

Graphical representation of the workflow

Methods

Construction of DermaGene: data sourcing, curation, and reassessment

We initiated by creating a comprehensive list of dermatological disorders by referring to various online and offline resource databases followed by identification and mining of data to make every entry in our database descriptive of the SNP–disease relationship. The three top databases that we referred to included UniProt, DisGeNET (Pinero et al. 2015), and OMIM.

Database creation was initiated by procuring the entire UniProt database, followed by converting the data obtained into the suitable and workable MS Excel format. A significant volume of data downloaded was delimited which required to be segregated into distinct columns, and we applied text to column—delimited to get a succinct view. This first stage download of database was subjected to filters from where we manually selected every possible dermatological disorder from the pool of 74,177 entries. This was achieved by referring to our previously generated comprehensive list of dermatological disorders. At this stage, we had the Gene name and rs ID indicating the SNP and the corresponding disease name. We thus reduced our database to 5000 entries, followed by removal of incomplete or duplicate entries thereby condensing it further to 2800 entries described by three attributes.

As the next step, we utilized disGeNET, to individually query every single gene (obtained from the UniProt database) to identify the p value and the PubMed ID. This was achieved by manually matching the results obtained with their corresponding rs IDs. Additionally, we populated the database with corresponding PubMed ID, p value, and disease class. Further, OMIM was put to use to support filling of gaps in the database. Extensive research helped us identify the missing values, which we were required to eliminate owing to insufficient data, thus reducing the database to 1116 entries.

Once we were finished with supplementing the database in its creation stage with requisite data, we proceeded to utilize PubMed as our primary source and to access all the research papers linked through their PubMed ID with our database. The work culminated into exhaustive reading of research papers to discern the ethnicities associated with the research conducted. Exhaustive study of research papers enabled us to determine ethnicities and to subsequently complete the database. Once the database was generated, the data was carefully re-assessed to ensure authenticity, through manual analysis of input data with the intention to establish the veracity of information as sourced from the specified databases. Besides, identification of alternative names for the specific diseases allowed for scope expansion of the database.

Each entry in this database therefore describes the genes, the SNPs associated with diseases, the PubMed IDs of texts referred, and the frequency of disease occurrence along with the population ethnicity where it exists. Also, we have taken into consideration the ethnicity of the genotypic variation to identify disease susceptibility in specific populations, the p value for the occurrence of genetic dermatological disorders which enables recognition of how strongly a particular gene influences the mentioned diseases. We thus have included six fields, and the data obtained was subsequently used to chart patterns within multiple genes and diseases. The information sourced aims to provide in-depth knowledge on common dermatological disorders, subsequently provisioning us to analyze the risk of genetic variants on their contribution to disease susceptibility.

Disease classification and network analysis of DermaGene data

Network theory involves the analysis of complex systems and serves as a potent biological tool, thereby, facilitating a better understanding of gene interactions and their role in causing diseases. To perform disease association network analysis, we constructed disease–gene bipartite network (DGN) and disease–SNP bipartite network (DSN) using Cytoscape (Shannon et al. 2003) with each network comprising a pair of disjoint set of nodes. The two sets of nodes are chosen such that there is no intra-connection but capable of exhibiting inter-connection. The nodes communicate with each other via edges, where the edges are the representatives of the nature of the interaction between the pair of nodes. In DGN, all the genes present in our database formed one set of nodes and the other set of nodes were reflected by all the known dermatological disorders in DermaGene. An edge connected the two sets of nodes if a disease–gene association was found. Whereas in DSN, all the diseases in our database formed one set of nodes and the other set of nodes comprised of the associated SNPs. The diseases were linked only if they showed common genetic variants. We also illustrated the correlation between the diseases as a result of a shared gene between the two. The diseases were selected as nodes while the edge attributes represented the genes common between them.

Data collection for vitiligo disease module analysis

To analyze the role of miRNAs associated with vitiligo, the miRNA information was obtained from published literatures available in NCBI PubMed server (https://www.ncbi.nlm.nih.gov/pubmed/) and mined publicly available online databases, namely, HMDD v2.0 (the Human microRNA Disease Database) (Li et al. 2014), Entrez GENE database of NCBI, and miRBase (Griffiths-Jones 2006) for relevant entries. A total of 41 types of miRNAs were found to be associated with vitiligo.

On the other hand, we extracted the information of genes associated with vitiligo and its SNPs from the GENE database and dbSNP of NCBI. The build 141 of NCBI dbSNP database is the latest release containing nearly 44 million validated human SNPs (Sherry et al. 2001). A total of 186 genes were reported to be associated with the disease among which 134 polymorphisms for 84 genes were reported to be positively associated. UniProt IDs of the proteins were also noted.

Identification of miRNA target genes and construction of miRNA–target gene network

The miRNAs identified to be associated with vitiligo were used to find their respective target genes. TargetScanHuman 7.1 (Agarwal et al. 2015) was used to detect targets in the 3′UTR of the protein-coding transcripts by base-pairing rules where predictions with both broadly conserved and poorly conserved sites are provided. These target sites are the conserved sites that match or are complementary to the seed region of the miRNA that ultimately facilitates the binding of miRNA with the mRNA to functionally degrade the mRNA thereby resulting in gene silencing.

A structured network layout explaining network integrity is the core requirement to justify the interaction between miRNA and disease. A miRNA–target gene interaction network was therefore constructed and analyzed to understand the miRNA–target gene relationship and validate the miRNA–disease association. The miRNA–target gene bipartite network consists of two sets of nodes—one set represents the miRNAs, and the other set represents the target genes. Nodes from the two sets were connected if a particular miRNA is associated with a particular target gene. The datasheet prepared which included all the miRNAs associated with vitiligo along with their target genes was used to generate the network in Cytoscape (Shannon et al. 2003).

Construction of vitiligo-associated protein–protein interaction network

Interactions between genes whose expression profiles are correlated with disease pathogenesis may contribute to the progression of the disease (Srivastava et al. 2017). To identify such interactions, we constructed the protein–protein interaction network using STRING (The Search Tool for the Retrieval of Interacting Genes) (Von Mering et al. 2005). It is a comprehensive database capable of providing an overall view of all the known and predicted protein–protein interactions of physical and functional associations. The PPI network generated based on STRING online database was then visualized in Cytoscape whose common feature lies in combing biological interaction networks with relevant large databases into a unified framework.

Functional module and enrichment analysis

Considering the connectivity properties of a network, we resolved to identify the significant clusters or modules enriched in biological processes from the complex bipartite network to extract biologically meaningful interactions. The modular analysis can provide a better insight into the relationship of the interconnected proteins assuming that the highly connected nodes in a network could form a cluster. As we know, cellular processes and functions are modular; therefore, it is more feasible to predict the structural and functional behavior of a particular module than that of an individual gene. Such module analysis has played a significant role in the past in determining disease mechanisms (Mitra et al. 2013). We used Markov Clustering Algorithm (MCL) method among the other clustering algorithms as it is the most widely used unsupervised clustering algorithms for functional module analysis and assigns a fast and reliable scalable method for finding functionally enriched clusters in complex networks. The granularity parameter for MCL clustering was kept at 1.8.

To interpret the biological impetus of the clusters in the network, we further performed functional enrichment analysis of the clustered gene groups using DAVID (The Database for Annotation, Visualization and Integrated Discovery) 6.8 (Huang et al. 2009). Functional enrichment analysis uses statistical approaches to identify the clusters of genes or proteins that are highly expressed or enriched in a large set of genes or proteins which may have an association with disease (Huang et al. 2008). All the genes were mapped into DAVID with the default settings selected as “Homo sapiens” in both the species background as well as the current background for the analysis.

Pathway analysis

To analyze the pathways in which these genes are involved, we conducted a pathway analysis using KEGG (Kyoto Encyclopedia of Genes and Genomes) pathways in DAVID 6.8. KEGG Pathway database is the most comprehensive and widely used database of annotation information. The pathway classification with p ≤ 0.001 and Benjamini–Hochberg FDR ≤ 0.01 was considered to have the most biological significance. Pathway analysis helps to interpret the data in the context of biological interactions and identify related proteins within a pathway.

A single protein may be involved in multiple pathways that are of importance to many biological processes. The biological cause of a disease can be explored by examining the changes in gene expression in a pathway. Also, the same pathway can be targeted for novel drug candidates (Wang et al. 2010). Deciphering the pathways which are explicitly targeted by the essential proteins may provide insight into their regulatory mechanisms. To better analyze the pathway analysis results, we constructed a protein pathway network to elucidate which pathways were eminently targeted by the proteins, specifically, the hub proteins. Such pathways, in turn, can be targeted for therapy and treatment approaches.

SNP analysis

The feasibility of the identification of the nsSNPs that vest susceptibility or resistance to human diseases has been improved with the use of in silico tools. To elucidate of the function of mutations in vitiligo-susceptible genes, we investigated the pathogenic effect of 134 SNPs which were reported to be associated with vitiligo. Among the 134 SNPs, the functional context of 36 nsSNPs was analyzed by employing various computational platforms. We used a combination of computational tools, namely, SIFT (Kumar et al. 2009), PolyPhen 2.0 (Adzhubei et al. 2010), PROVEAN (Choi et al. 2012), SNPs&GO (Magesh and Doss 2014), I-Mutant Suite 3.0 (Capriotti et al. 2008), and PANTHER Evolutionary Analysis of Coding SNP (Mi et al. 2005) to identify the nsSNPs that potentially affect the structure and function of proteins associated with vitiligo.

Additionally, MutPred 2 (Pejaver et al. 2017) was used to interpret the possible molecular cause of disease-inducing amino acid substitutions. Furthermore, NetSurfP (Petersen et al. 2009) was used to analyze the effect of such mutations in the stability of the protein by predicting the solvent accessibility of the substituted residue.

The information from different computational platforms was combined to prioritize the deleterious SNPs to increase the predictive power and accuracy of the results of the in silico techniques.

Construction of drug–target network and drug–similarity network

A single drug can target multiple proteins. To analyze this relationship between drug and protein targets (disease–gene products) and to understand how they intervene therapeutically in disease processes, we constructed a drug–target network (Yildirim et al. 2007). We extracted information about drugs with respect to our drug candidates from DrugBank (Knox et al. 2011), which is a chemo-informatics resource that is updated and maintained with the Food and Drug Administration (FDA) information. We used known FDA drugs (both approved and approved-investigational) concerning to the protein targets to generate the drug–target bipartite network.

To identify which drugs act similarly on the same target, we constructed a drug–similarity network to analyze the interactions of drugs and its action on protein targets. DrugBank provides clinically relevant drug interactions, and the information of the interacting drugs and its interacting mechanism with respect to a particular drug was taken from DrugBank. Such a network can provide insights on drug–drug interaction for a potential drug candidate where an interacting drug can act synergistically or antagonistically with another drug altering the benefit or effectiveness of the drug on disease conditions. An interacting drug of a particular drug can also pose as a potential drug for the protein targeted by that drug.

Protein modeling of TLR2

Finding the 3D structure of proteins is helpful in predicting the impact of SNPs on the structural level and in showing the degrees of alteration. To elucidate the molecular dynamic behavior of the SNP in TLR2 protein, we performed preliminary protein modeling by G23D (Genomic variant to 3D protein data) (Solomon et al. 2016) which is a tool for the conversion of human genomic coordinates to protein structures. G23D allows the mapping of evolutionary related as well as identical protein of genomic variants in a 3D model structure assisting in the feasibility of structural insight. Along with the mutated sites, it also displays the wild-type residue and other functional sites on the modeled 3D protein structure to facilitate better interpretation of the variant.

We further carried out TLR2 protein–protein interaction analysis to analyze its interaction with other proteins which might be influenced by the mutation in TLR2 protein.

Construction of VitmiRS

The primary data in VitmiRS constitutes the collection of miRNAs and SNPs that have been reported to be associated with vitiligo. The process for database creation relied on extensive review of high-quality published journals, relevant articles published in PubMed, and multiple online resources discussing the association of miRNAs and SNPs with vitiligo for the identification and the subsequent extrication of relevant mentions of the fields that we explored. Information from multiple comprehensive and relevant databases, namely, dbSNP (Sherry et al. 2001), UniProt, HMDD v2.0 (Li et al., 2014), miRBase (Griffiths-Jones 2006), and GENE database, was also taken into consideration. Every paper was read manually to eliminate the errors that may arise from machine-generated information. The information included was filtered at the thresholds of p value ≤ 0.05 and carefully re-assessed through manual analysis to ensure authenticity. A user-friendly web interface was developed for the ease of data retrieval. Apache was used as the application server and PHP as the programming language. The presentation layer of VitmiRS was created using XHTML and CSS, while for the backend database, MySQL was used.

Results

DermaGene

We curated a total of 244 unique genes in our database, with each gene exhibiting multiple SNPs and describing 114 diseases. A total of 871 unique SNPs exist that exhibits both pleiotropy and epistasis over dermatological disorders. The present DermaGene version relies on the review of over thousand literature sources to generate a comprehensive source comprising of 1116 tuples described through six attributes. It provides a user-friendly interface for easy retrieval of data to query detailed information on each SNP–disease association and hopes to facilitate analysis of the comorbidity associations asserted between human variation and observed disease conditions.

Disease classification and network analysis of DermaGene data

Disease classification

Analysis of DermaGene database revealed the existence of maximum SNPs accounting for dermatological disorders to lie on chromosome 17. The distribution of SNPs was also analyzed to observe their association with the different disease classes included in our database. We identified that the primary disease class which exists in our database is skin and connective tissue diseases (632). This indicates that the maximum number of SNPs, as in our database accounts for the skin and connective tissue class of diseases, followed by neoplasia (116) and nutritional and metabolic diseases (78). The results are summarized in Fig. 3. Our database reveals the possibility for disease comorbidity existence. We have identified few diseases that share the SNP with other diseases, indicating association among them and the subsequent combinatorial occurrence of the diseases. This shows that in the presence of co-morbid associations, the incidence of one of them in an individual may elevate the possibility of another disease to occur. An additional file (see Additional files 1 and 2) lists the 12 incidences of disease comorbidity arising from shared SNPs and 22 disease–gene associations, identified from DermaGene respectively which may be responsible for comorbid associations.

Fig. 3

SNP distribution and association in DermaGene database. a Chromosomal distribution of SNPs. b Pie chart depicting the distribution of disease class in DermaGene

Disease–gene bipartite network (DGN)

Our analysis of DGN (Fig. 4a) identified vitiligo from our database to be under the influence of a maximum number of different genes since it exhibited maximum disease correlations, therefore, indicating a probable multi-factorial genetic inheritance. From our database, various genes are showed a higher degree of influence compared to other genes which are defined by the number of diseases it influences (Fig. 4b). The genes that exhibit maximum association include IL2RA, PTPN22, and FGFR2 and are shown in Fig. 5. We observed that the gene IL2RA causes three distinct diseases—alopecia areata, urticaria pigmentosa, and vitiligo. At the same time, another gene PTPN22 accounts for alopecia areata, systemic lupus erythematosus, and vitiligo (Fig. 5a). The third gene IKZ4 also contributes to alopecia and vitiligo, which confronts us with possible idea of a genetic association between these diseases. The relationship between alopecia and vitiligo has been validated through previously conducted research studies (Huang et al. 2013; Krishnaram et al. 2013). Similarly, the gene FGFR2 contributes to Apert syndrome, Crouzon syndrome, and Beare–Stevenson cutis gyrata syndrome (Fig. 5b), in spite of the absence of definite literature, indicating epistasis, if at all. This means a probable novel finding which can be useful for doctors and researchers alike.

Fig. 4

Disease–gene bipartite network (DGN). a Vitiligo disease–gene bipartite network. The node size is indicative of the number of associations. b Disease–disease correlation network. The node size is indicative of the number of associations

Fig. 5

Mapping disease comorbidity in disease–gene bipartite network (DGN). a Maximum association of the genes IL2RA and PTPN22 with multiple diseases. b Maximum association of the gene FGFR2 with multiple diseases

Disease–SNP bipartite network (DSN)

From our DSN (Fig. 6), we inferred that most SNPs are associated with single diseases, although we also report the occurrence of few SNPs that exhibit pleiotropic effect since they are linked to multiple diseases. We identified the SNPs rs28931615, rs121912737, rs2476601, rs121913478, rs28931594, rs11540654, rs1042522, rs59856285, and rs60723330 from our DermaGene database to reveal a probable pleiotropy effect defined through association with multiple diseases. Published studies explicitly authenticate our findings for the already existing relationship between acanthosis nigricans and Crouzon syndrome (Mir et al. 2013); acrokeratosis verruciformis and Darier disease (Harman et al. 2016); systemic lupus erythematosus and vitiligo (Nath and Kelly 2001); alopecia and vitiligo (Krishnaram et al. 2013); alopecia and systemic lupus erythematosus (Werth et al. 1992); ichthyosis hystrix with deafness and keratitis–ichthyosis–deafness syndrome (Dalamon et al. 2016); associations of ectodermal dysplasia (ED) and incontinentia pigmenti (IP) caused by IBKG gene, as resulting from complete loss of function in case of IP and X-linked recessive inheritance as in the case of ED (Zonana et al. 2000); Costello syndrome and seborrheic keratosis caused by germline mutation of Ras/MAPK signaling pathway (Hafner and Groesser 2013); and keratoderma palmoplantar which is often accompanied by pachyonychia congenita as described by Smith et al. (2000). There is thus a likelihood of epistatic effect and that the shared SNPs serve as an index of the biological relationship existing between them with subsequent co-manifestation of these diseases.

Fig. 6

Disease–SNP bipartite network (DSN). The triangular nodes represent the SNPs, and circular nodes indicate the disease. The size of circular nodes is proportional to the SNPs that the particular disease is associated with. a The diseases xeroderma pigmentosum and keloid share two SNPs rs1042522 and rs11540654. b Relationship between the diseases Beare–Stevenson cutis gyrata and Brooke–Spiegler syndrome as caused by the SNP rs121913478

We also found that diseases xeroderma pigmentosum and keloid share two SNPs rs1042522 and rs11540654, which served as a revelation on the previously unidentified association between the two (Fig. 6a). Our database sheds light on the relationship between the diseases Beare–Stevenson cutis gyrata and Brooke–Spiegler syndrome as caused by the SNP rs121913478 which exhibits no prior known correlation (Fig. 6b).

We have identified various other genetic associations through network analysis as well, wherein we observe the role of a single gene in causing various dermatological disorders. Published evidence for a probable relationship between alopecia and Marie–Una syndrome can be inferred from a published study conducted on a seven generation British family, although the study is not directly indicative of the particular association as such (Marren et al. 1992). Xeroderma pigmentosum (XP) and Cockayne syndrome (CS), although phenotypically different, share genetic similarities, since the gene ERCC5 is responsible for both the cases and this is probably because the several mutations accounting for UV sensitivity may also account for somatic and brain growth failure, thus resulting in neurological problems. Additionally, an existence of literature for LMNA gene contributing to both Hutchinson–Gilford progeria syndrome and congenital muscular dystrophy establish the authenticity of our findings from the DermaGene database (Barateau et al. 2017). Also, we find that Menkes syndrome is a more severe variant of cutis laxa, the latter being initiated during early middle childhood. Both the diseases arising from a mutation in ATP2A2 gene exhibit loose skin as a phenotypic contribution to the affected individual. Lastly, previous reports show clear progression of chilblain lupus erythematosus to systemic lupus erythematosus (SLE) (Bansal and Goel 2014), with around 20% individuals developing SLE (Hedrich et al. 2008).

We also identified the gene ABCC6 located at chromosome 16 on short arm (p) at position 13.11, to have the maximum SNPs as shown in Fig. 7. This gene contributes to pseudoxanthoma elasticum and hence can play a key in disease prognosis. With a total of 91 SNPs within the gene, ABCC6 gene serves as a disease hotspot, with the existence of multiple SNPs contributing to the same disease. This reveals that ABCC6 gene is a relatively large gene, of size 74,766 bp to allow for a high number of SNPs.

Fig. 7

Mapping ABCC6 gene to pseudoxanthoma elasticum. Multiple SNPs of ABCC6 gene are shown to be contributing to the same disease pseudoxanthoma elasticum

Data collection for vitiligo analysis

We found 41 miRNAs to be associated with vitiligo with 34 upregulated and 7 downregulated in disease conditions. A total of 134 polymorphisms for 84 genes were found to be associated with vitiligo. Among the 134 SNPs, 36 are nsSNPs, 7 are coding synonymous SNPs, and 68 SNPs were found in the non-coding region.

Identification of miRNA target genes and construction of miRNA–target gene network

The respective target genes of the miRNAs were identified to explicate the biological targets of these miRNAs whereby they control gene expression and ultimately regulate the cellular and molecular responses during disease development and progression. The miRNA–target gene bipartite network represented a total of 41 miRNAs and 98 unique target genes consisting of 139 nodes and 220 edges (Fig. 8). The target genes are considered to be connected in the network if they share a common miRNA. To identify the hub miRNAs and target genes associated with vitiligo, among the number of methods available for hub identification, we chose maximum clique centrality (MCC) along with two other topological parameters, namely, betweenness centrality and bottleneck (Chin et al. 2014), and normalized the data to identify the top hubs in the network. In our analysis, we found seven hub miRNAs (hsa-miR-99b, hsa-miR-577, hsa-miR-9, hsa-miR-155, hsa-miR-211, hsa-miR-10a, and hsa-miR-145).

Fig. 8

miRNA–target gene interaction network. The green colored nodes represent the miRNAs, the white colored nodes represent the target genes, and yellow colored nodes represent the hub miRNAs along with the target genes of the highest degree in the network

Protein–protein interaction network

Protein–protein interaction network of vitiligo-associated proteins exhibited significant interconnections between the proteins. It comprised of 71 nodes and 322 edges (Fig. 9). The proteins are considered to be connected in the network if they interact with each other. If each protein is not inclined to interact with another protein or specifically interacts with a single protein, then the bipartite network would be disconnected into many single nodes corresponding to specific or unique proteins with few or no edges between the nodes. Rather, the protein–protein interaction network generated displayed many interactions between the proteins. For prioritizing proteins as hubs, we chose MCC along with betweenness centrality and bottleneck as topological parameters (Chin et al. 2014) and normalized the data to identify the top 15 hubs in the network. These 15 essential (hub) proteins are IL10, IFNG, IL4, CD44, IL1B, CTLA4, GZMB, FOXP3, TNF, IL2RA, CAT, ESR1, TLR2, HLA-A, and GSTP1.

Fig. 9

Protein–protein interaction network. The purple nodes represent the hub proteins in the network

Functional module and enrichment analysis

The functional module analysis of vitiligo-associated proteins revealed four functional modules (Fig. 10). The average size of the clusters was 17.75, and they were ranked by their modularity score of 1. The majority of the proteins were found to form a single large cluster. This implies that these proteins have a biological similarity in their functions.

Fig. 10

Functional module (clusters) network. The shades of green color, from light to dark, represent the decrease in the number of interactions. The genes with the maximum number of interactions show lighter shades in the cluster

We performed functional enrichment analysis of the larger functional module consisting of 64 proteins using DAVID 6.8. Keeping the classification stringency at highest and considering the enrichment score value ≤ 1.3 to be significant, the given set of target genes was classified into nine functionally enriched clusters (see Additional file 3) that involved 15 genes from the given set of 64 genes. All the nine clusters were observed to be primarily associated with the immune system regulatory processes, such as MHC class I/II-like antigen recognition protein, conserved site of immunoglobulin/major histocompatibility complex, conserved site of interleukin-10, positive regulation of JAK-STAT cascade, apoptotic signaling, chemokines, and TNF.

Pathway analysis

According to the threshold of hypergeometric test p ≤ 0.001 and Benjamini–Hochberg FDR ≤ 0.01, the mapped genes were found to be enriched in a total of 30 pathways, which involved a total of 40 genes from the given gene set. The most significant pathway was found to be the allograft rejection pathway with a p value of 5.70E−17 and FDR value of 1.51E−14 involving 13 genes. Among these 30 pathways, almost half of the pathways (15 pathways) were associated with immune system responses and related disorders, and autoimmunity, while 14 pathways were associated with infectious (viral, bacterial, and parasitic) disease pathways that have been reported to impair the proper regulation of the immune system.

In the protein–pathway network, a protein is connected to a pathway if the protein is known to be involved in that particular pathway. The network consisted of 71 nodes and 270 interactions representing 40 proteins and 30 pathways (Fig. 11).

Fig. 11

Protein–pathway network. The red colored nodes represent the proteins, the dark blue colored nodes represent the pathways, and the light blue colored nodes represent the hub proteins

Based on degree value, TNF showed maximum interaction in the network implying its participation in most of the pathways. Among the 30 pathways, TNF was found to be associated with 24 pathways. Similarly, IL1B, IFNG, IL10, and TLR2 were found to be involved in 19, 18, 15, and 14 pathways respectively. It was also observed that TLR2 was specifically associated with the infectious disease pathways that were responsible for deregulating the immune system processes as depicted in the results of KEGG pathway see analysis. An additional file (see Additional file 4) lists the pathways in which the hub proteins were found to be involved. Among the 15 hub proteins identified in the protein–protein interaction network, 11 were found in this network suggesting that these 11 essential proteins are involved in the filtered significant pathways.

SNP analysis

Among the 134 SNPs, 36 (26.9%) are nsSNPs and 7 (5.2%) are SNPs, while 68 (50.8%) SNPs were found in the non-coding region. SNPs in the non-coding region comprises of 44 (32.8%) SNPs in the intronic region, 19 (14.2%) in the near-gene region, and 5 (3.8%) in the mRNA UTR region. The rest 23 (17.1%) among 134 SNPs are intergenic. The distribution of the SNPs is shown in Fig. 12a.

Fig. 12

Distribution of SNPs. a Graph of the distribution of SNPs. b Graph of the distribution of the deleterious and benign nsSNPs as predicted by SIFT, PolyPhen, PROVEAN, I-Mutant Suite, SNPs&GO, and PANTHER

For our analysis, we selected the nsSNPs and UTR-region SNPs since UTRs are central for the post-transcriptional regulation of gene expression, and alterations in the functional UTR region can lead to severe pathology (Conne et al. 2000). The 36 nsSNPs were analyzed by using a combination of SIFT, PolyPhen, PROVEAN, SNPs&GO, I-Mutant Suite, and PANTHER Evolutionary Analysis of Coding SNP tools.

SIFT predicted six nsSNPs (16.7%) to be deleterious with a tolerance score cut off which is ≤ 0.05. Further analysis of the nsSNPs using PolyPhen predicted two nsSNPs to be “probably damaging” and five nsSNPs to be “possibly damaging” with a tolerance cut off score ≥ 0.5. Consequently, seven nsSNPs (19.4%) were characterized as damaging. Analysis using PROVEAN revealed a similar result as that of PolyPhen tool. Based on a tolerance cut off score value ≤ 0.05, it predicted seven nsSNPs (19.4%) to be damaging. Of these seven nsSNPs, one (rs5743708) was reported to be highly damaging with a tolerance score of 0.

To improve the prediction accuracy, we used I-Mutant Suite and SNPs&GO tools for further analysis. The output results in I-Mutant Suite predict the free energy change (DDG) and classifies it into any three of the mentioned classes, namely, largely unstable (DDG < − 0.5 kcal/mol), largely stable (DDG > 0.5 kcal/mol), or neutral (− 0.5 ≤ DDG ≤ 0.5 kcal/mol). The output result of a protein as largely unstable indicates the effect of a deleterious nsSNP. In our I-Mutant Suite results, we found that 11 nsSNPs (30.6%) exhibited a DDG value of less than − 0.5, which indicates that these are largely unstable resulting in disease-associated mutations. Whereas, SNPs&GO predicted 10 nsSNPs (27.8%) to be related to disease-associated mutations by using a tolerance cut off score ≥ 0.5.

A disease-causing mutation in the highly conserved regions of the genome may pose a high risk of that particular mutation to be damaging which is why we decided to carry out SNP conservation analysis. The probability of a mutation to cause a disease increases monotonically with the increase in the degree of site conservation (Vitkup et al. 2003). Conservation analysis by PANTHER Evolutionary Analysis of Coding SNP predicted 14 nsSNPs (38.9%) to be deleterious based on their preservation time. A longer preservation time implies a greater possibility of functional impact on a protein (Mi et al. 2005). Since different in silico tools have a diverse set of alignments and molecular characteristics, the results of the six tools were slightly different. Accordingly, we combined the results of SIFT, PolyPhen, PROVEAN, SNPs&GO, I-Mutant Suite, and PANTHER Evolutionary Analysis of Coding SNP (see Additional file 5) to predict the deleterious nsSNPs common in all the analysis.

Figure 12b shows the distribution of deleterious and benign nsSNPs obtained using SIFT, PolyPhen, PROVEAN, I-Mutant Suite, SNPs&GO and PANTHER.

Of all of the predictions, 16.7%, 19.4%, 19.4%, 30.6%, 27.8%, and 38.9% deleterious nsSNPs were specifically found by SIFT, PolyPhen, PROVEAN, SNPs&GO, I-Mutant Suite and PANTHER respectively. Combining the results of all the six tools, three nsSNPs, namely, rs1801133 (MTHFR), rs5743708 (TLR2), and rs11575993 (SOD2), were predicted to be functionally significant. Tables 1 and 2 presents the deleterious nsSNPs obtained through the SIFT, PolyPhen, PROVEAN, I-Mutant Suite, SNPs&GO, and PANTHER analysis of the vitiligo-associated nsSNPs.

Table 1 nsSNPs found to be deleterious using SIFT, PolyPhen, and PROVEAN

Table 2 nsSNPs found to be deleterious using I-Mutant Suite, SNPs&GO, and PANTHER

MutPred predicted the molecular cause of the nsSNPs to become deleterious in MTHFR (rs1801133), TLR2 (rs5743708), and SOD2 (rs11575993). Analysis of the results showed an interrelation of the SNPs to be damaging with respect to the solvent accessibility of the protein. The type of mutated residue and its position in the sequence affect the stability of the protein which due to mutation decreases with the decrease in solvent accessibility of a residue (Vitkup et al. 2003).

NetSurfP predicted the surface solvent accessibility of amino acids by using the protein FASTA sequence as a query. The solvent accessibility is predicted to be buried or exposed, based on the accessibility of the amino acid residues to the solvent. The reliability of relative surface accessibility is verified in the form of Z-score which highlights the surface prediction reliability. As given in Table 3, the class assignment does not change for the three nsSNPs. Although there were very minimal changes in the relative surface accessibility (RSA) values for the three nsSNPs, a considerable drift in the Z-score was not observed between the wild-type and mutant-type proteins. A decrease in RSA value has been observed in the mutant types of MTHFR and SOD2 while there was an increase in the RSA value of the TLR2 mutant type.

Table 3 Solvent accessibility analysis of the mutated proteins by NetSurfP

Drug–target network and drug–similarity network

Drug–target network

The bipartite network of drug–protein target interaction consisted of 109 nodes and 84 interactions (Fig. 13a). A drug and protein are considered to be connected to each other if the protein is a known target of the drug, giving rise to a drug–target network.

Fig. 13

Polypharmacological analysis. a Drug–target network. The light blue colored nodes represent the hub protein targets, the dark blue colored nodes represent the protein targets, and the pink colored nodes represent the drugs. b Drug–similarity network. The pink colored nodes represent the hub protein targets, the yellow colored nodes represent the protein targets, the green colored nodes represent the drugs, and the blue colored nodes represent the interacting drugs

In our analysis, we found that the most of the drugs targeting a particular protein did not show any interaction with other protein targets in the network except for the two drugs, etanercept and carfilzomib. Etanercept targets both TNF and LTA while carfilzomib targets PSMB8 as well as PSMB9. Also, we found that out of the 15 hub proteins, only 9 protein targets were found to be present in this network. This illustrates the other 6 hub proteins (IL10, IL4, GZMB, FOXP3, TLR2, AND HLA-A) as potential drug candidates for which drug information is currently not available. Also, we found that 4 hub proteins, namely, CD44, CD152, CAT, and GSTP1, were targeted by a single drug. This highlights the imperative need to discover more effective drugs that target these proteins which may play a major role in therapeutics to alleviate disease conditions in patients.

Another notable finding in our analysis was that the drugs which showed a high degree in the network were mostly indicated for the treatment of autoimmune diseases and deregulated immune responses. One of them is etanercept that targets TNF, a major proinflammatory cytokine that affects various aspects of the immune response, and LTA as well. Etanercept is a genetically engineered decoy receptor that consists of the ligand-binding domain of TNFR2 and the Fc component of human IgG1. It competitively binds with high affinity to TNFR2 inhibiting the binding of both TNF-α and TNF-β to the cell surface receptors, consequently, inhibiting inflammation induced melanocyte death. It has been indicated to be clinically used for rheumatoid arthritis, psoriatic arthritis, ankylosing spondylitis, and Crohn’s disease (Nanda and Bathon 2004). However, etanercept has been reported to be less efficient as a monotherapy in vitiligo patients requiring the need of a combinative therapy (Rigopoulos et al. 2007). LTA, on the other hand, is involved in the follicular dendritic cells development and has been observed to induce signals leading to lymphoid neo-organogenesis driving the inflammatory responses in autoimmune diseases like rheumatoid arthritis (Takemura et al. 2001). This suggests that etanercept may suppress lymphoid neo-organogenesis and reduce the proliferation of mature dendritic cells in vitiligo lesions (Wang et al. 2011). The other drug, carfilzomib, is a tetrapeptide epoxyketone-based proteasome inhibitor that targets PSMB9 and PSMB8. Peptides generated from ubiquitin-tagged cytosolic proteins are presented to CTLs by MHC class-I molecules which are degraded by multi-catalytic, cytosolic immune-proteasome complex called LMP2 and LMP7 encoded by PSMB9 and PSMB8 genes respectively (Cresswell et al. 2005). This intrinsic enzymatic activity of immune-proteasomes may be altered by genetic variations which reduce the expression of PSMB8 and PSMB9 in vitiligo PBMCs after IFNG stimulation. This leads to defective proteolytic degradation and accumulation of ubiquitinylated proteins in the epidermis of vitiligo patients leading to ROS production and auto-inflammatory immune responses which may be detrimental for the manifestation of vitiligo (Dani et al. 2017). Carfilzomib targets the catalytic activity of immune-proteasomes and irreversibly inactivates the proteasome, thereby inhibiting aberrant immune function (Miller et al. 2013).

Drug–similarity network

The drug similarity tripartite network of protein targets, drugs, and interacting drugs comprised of 178 nodes and 1331 interactions (Fig. 13b). A drug and its interacting drug are considered to be connected if they share a common protein target. Interacting drug partners of a particular drug may enhance the efficacy of the drug or may even target the same protein. Such interacting drugs may represent themselves as potential drug repositioning candidates. With this concept, we constructed the drug–similarity network that displayed interconnections between drugs and their interacting drug partners.

Apart from the drugs etanercept and carfilzomib targeting more than one protein as shown in the drug–target network, we found two other drugs showing interaction with another protein target in the network. Both diethylstilbestrol and conjugated equine estrogens targeting ESR1, which is one of the hub protein identified in the protein–protein interaction network analysis, were found to be interacting with the drug entacapone targeting COMT. This signifies that the two drugs interacting with entacapone might target COMT which is targeted by entacapone alone. Also, we found that there are no interacting drugs reported for etanercept targeting LTA and our drug similarity network analysis suggests that the interacting drugs for etanercept targeting TNF might as well target LTA. Etanercept has been less efficient as a monotherapy as mentioned earlier; therefore, the interacting drugs for etanercept as shown in the network might catalyze its efficiency when used in combination. Thus, further comprehensive study is required to validate the effectiveness of these drugs in combinative therapy. In addition, GSTP1, also a hub protein, was shown to be targeted by a single drug (clomipramine) in the drug–target network; however, seven interacting drugs were shown to be connected to this drug in this network. These interacting drugs can be further studied to investigate their potential as drug repositioning candidates for vitiligo treatment.

Protein modeling of TLR2

Among the three predicted functionally significant nsSNPs, we performed protein modeling of TLR2 since our drug–target network analysis showed TLR2 as a potential drug target for which drug information is currently not available. The G23D tool generated a full-length protein model for the TLR2 protein based on E-value which was 0.001 by default as generated from the PDB library and ModBase hits. Since we did not perform any molecular dynamics simulation studies for protein structure optimization, therefore, the 3D homology model generated by G23D is a preliminary model implicating the disruptive role of the SNP (rs5743708) on TLR2 protein.

Figure 14a shows the cartoon representation of the modeled protein structure including both the wild-type (green) and the mutant (light blue) amino acid. It also displayed the other variants as reported in ClinVar, COSMIC, and dbSNP for TLR2 which is colored as red, yellow, and dark green respectively in stick representation in the 3D structure.

Fig. 14

Protein modeling. a 3D structure of modeled mutant (R753Q) TLR2 protein. The wild-type residue is represented in green color, and the mutant residue is colored light blue. Other variants are also displayed in this structure with the dbSNP variant colored dark green, COSMIC variant colored orange, and the ClinVar variant is represented in red color. The backbone of the protein is shown in gray color. b Disease mechanism influenced by R753Q TLR2 polymorphism

An A/G mutation (rs5743708) caused the substitution of the amino acid from arginine into glutamine at position 753 (R753Q) of the protein. The mutation was identified to be located within the TIR domain as annotated in UniProt. Since it is a membrane protein, it is exposed on the surface which is in agreement with the previous NetSurfP results. Analysis of the 3D structure suggests that the residue glutamine in position 753 is located in the middle of a helix. There was a difference in the size observed in the wild-type and mutant amino acids with the mutant residue being smaller than the wild-type residue causing conformational changes on the DD loop which might lead to loss of interactions. Also, the wild-type residue is positively charged, but the amino acid substitution leads to a decrease in the positive charge that changes the interaction surface within the TIR domain via altered electrostatic potential. These may affect TLR2 dimerization causing loss of interactions with other molecules or residues affecting the functional activity of the protein (Xiong et al. 2012). The R753Q polymorphism compromises the TLR2/1 or TLR2/6 assembly resulting in deficient tyrosine phosphorylation and impaired recruitment of MyD88. This reduces the phosphorylation of IRAK1 and diminishes the activation of MAPKs and NF-κB resulting in the deficient production of cytokines thus altering TLR2 signaling competence. Therefore, the compromised signaling is due to deficient tyrosine phosphorylation and not due to lower mutant receptor expression. The reduced activation of NF-κB signaling pathway results in melanocyte apoptosis suggesting their decisive role in the increased risk for the development of vitiligo (Karaca et al. 2013; Traks et al. 2015).

The protein–protein interaction analysis of TLR2 showed its interaction with TLR1, TLR6, LY96, MyD88, IRAK4, TRAF6, NFKB1, TIRAP, CD14, and HMGB1. All these proteins were found to have a fundamental role in regulating the innate immune responses. The amino acid substitution of arginine to glutamine at position 753 has been reported to be located in the TIR domain. TLR2 activates NF-κB in combination with either TLR1 or TLR6 by bringing together the TIR domains and triggering tyrosine phosphorylation. The TIR domain of TIRAP binds to the TIR domain of TLR2 which then extends a docking platform for MyD88 recruitment. MyD88 facilitates the recruitment of IRAK4 which activates TRAF6 leading to the translocation of activated NF-κB to the nucleus where it induces target gene expression generating inflammatory responses (Oliveira-Nascimento et al. 2012). Mutations in the TIR domain tend to have more severe impact on signaling than those affecting the extracellular domain (Karaca et al. 2013). Genetic alterations in TLR2 might affect its interaction with both TLR1 and TLR6 consequently affecting NF-κB activation (Brown et al. 2006). It was also found that this mutation significantly reduced NF-κB activation by about 50 and 75% (Merx et al. 2007; Ben-Ali et al. 2011). In the case of vitiligo, altered NF-kB signaling may result in impaired melanogenesis inciting human melanocytes susceptible to TNF induced apoptosis (Shang et al. 2002) (Fig. 14b).

This implicates the possible detrimental effect of this mutation on the interaction of TLR2 with TIRAP, MyD88, IRAK4, TRAF6, and the consequent signal transduction. HMGB1 has been found to upregulate the expression of TLR2 of the NK cells and promote NK cell activation mediating inflammatory responses (Qiu et al. 2014), while CD14 (a co-receptor of TLR2) accelerates the microbial ligand transfer from CD14 to TLR2, resulting in an increased TLR2 signaling (Raby et al. 2013). LY96 have been found to enhance the expression of both TLR2 and TLR4 and vice versa enabling them to generate highly sensitive responses to a broad range of microbial lipopolysaccharide (LPS) structures (Dziarski et al. 2001). Altered expression of TLR2 will invariably affect the interaction of TLR2 with these proteins thereby influencing TLR2 signaling and inflammatory responses resulting in defective immune response to some antigens such as viruses in the case of vitiligo (Karaca et al. 2013).

VitmiRS

Currently, VitmiRS has information on 41 miRNAs and 134 associated SNPs for 84 unique genes. VitmiRS provides a user-friendly interface to detailed information on disease association of miRNAs and their target genes, SNPs and their position and location, the change in alleles and the resulting amino acid change due to such mutational events, and their respective literature reference. The users can search and access information for their respective query in the database browsing through query categories of miRNAs, SNPs, genes, and amino acid change. Further, the database interface also allows for the selection of particular attributes, such as PubMed ID, miRNA sequences, chromosome number, chromosome location, odds ratio, and p value; the geographical location and population in which the association of specific miRNAs and SNPs with vitiligo has been reported. Also, it incorporates the results of SNP analysis, thereby providing information on the pathogenic effect of the each SNPs as well. Thus, the database can be queried individually or in combination with any of the categories as mentioned above according to the need and interest of the user. The information available in VitmiRS may, therefore, be used for various associative studies relating to miRNAs and genetic variants to uncover their role and importance in vitiligo disease susceptibility and pathogenesis.

Discussion

An enhanced understanding of the disease causative genes and their relationship with other genes along with other diseases offers an improved resource for finding the candidate genes that mark the onset of a disease. Several unrelated diseases have been revealed to show common molecular mechanisms with strong association among them (Yu et al. 2015). As a result of such disease associations, the possibility of related diseases to occur together in an individual cannot be ruled out. Analysis of disease-gene bipartite network allowed the identification of several diseases which are under the influence of multiple genes wherein it highlighted vitiligo to be influenced by a multiplex of genes. Additionally, we also identified the genes which exhibit pleiotropy, thus accounting for multiple diseases and may play a crucial role in determining disease comorbidity. We explored some unique relationships with their previous published studies being non-existential which include the genetic association between rosacea and vitiligo, as caused by the gene BTNL2; the association between GM1-gangliosidosis and eczema, both being caused by GLB1 gene; and the association between piebaldism and urticaria pigmentosa as a result of a mutation in the KIT gene. This provides with a new perspective in understanding dermatosis opening new avenues for conducting genetic research. In the DSN analysis, we found shared SNPs between diseases xeroderma pigmentosum and keloid whose association has been previously unidentified. Also, our analysis shows a correlation between the diseases Beare–Stevenson cutis gyrata and Brooke–Spiegler syndrome. Our database, therefore, provides a fresh perspective to the possible association between these diseases, as no published literature source indicating their correlation is available as of now. However, validation of our findings through wet-lab research can further strengthen our findings.

Based on our network analysis results of disease–gene association, we further carried out an integrative systems analysis of vitiligo. Vitiligo is a result of convoluted interactions of biological, environmental, and immunological events; hence, a single concept cannot be attributable to all the conditions of functional melanocyte loss. A large-scale analysis and integration of miRNA–disease associations will offer a platform to investigate the patterns of the miRNAs and its associations with diseases. Molecular signatures of miRNAs as reported in vitiligo patients suggest that these are actively involved and have a significant role in disease pathogenesis. Having a cardinal role in maintaining physiological homeostasis and disease development and progression, miRNAs are significant for melanocyte development and survival (Mansuri et al., 2016). Our miRNA–target gene network analysis revealed seven hub miRNAs, namely, hsa-miR-99b, hsa-miR-577, hsa-miR-9, hsa-miR-155, hsa-miR-211, hsa-miR-10a, and hsa-miR-145 implicating their role in vitiligo pathogenesis. The upregulation of hsa-miR-99b reduces the cytotoxic activity (cytokine effector functions) of NK cells which are crucial for the normal BCR signaling and proliferation of B cells. This causes deregulation of genes involved in B cell maturation and development resulting in the dysfunctioning of the immune system indicating them to be important players in vitiligo immunopathogenesis (Nandagopal et al. 2014; Šahmatova et al. 2016).

TYRP1 is targeted by hsa-miR-577, and its reduced expression as induced by miR-577 leads to increased sensitivity of melanocytes to oxidative stress causing early cell death of vitiligo melanocytes (Manga et al. 2006; Sturm and Duffy 2012). Also, the downregulation of PTPN22 was observed to be influenced by miR-577 which triggers the overexpression of T cells and suppresses anti-apoptotic AKT kinase inducing melanocyte destruction, thereby, rendering susceptibility to autoimmunity in vitiligo patients (Mansuri et al., 2016). Elevated levels of SIRT1 have been reported to protect cells from oxidative stress and inflammatory microenvironment (Han et al. 2008). Increased expression of miR-9 downregulates SIRT1 resulting in melanocytes apoptosis in vitiligo (Saunders et al. 2010). SIRT1 has been shown to regulate stress-activated MAPK pathway via Akt and ASK1 in vitiligo keratinocytes (Becatti et al. 2014). Previous studies have demonstrated the influence of miR-145 on the genes involved in the pigmentation process (Dynoodt et al. 2013). The genes targeted by miR-145 also regulate MAPK pathway along with JNK and TGFB signaling pathway and are related to the functional groups that might indirectly influence cellular processes in vitiligo wherefore they interfere with melanocytes function and viability (Šahmatova et al. 2016) suggesting the role of both miR-9 and miR-145 in the destruction of melanocytes in vitiligo.

TGFBR2 is targeted by miR-211, and downregulated miR-211 increases the expression of TGFBR2 which in turn downregulates MITF. MITF is a known primary regulator of melanocyte development and its survival (Levy et al. 2010); thus, deregulated MITF will considerably affect melanocyte development implying the role of miR-211 in vitiligo pathogenesis. Overexpression of miR-155 was found to modulate the levels of several interferon-regulated genes, such as SOCS1, IFITM1, and IRF1 that inhibits the expression of melanogenesis-associated genes, such as, TYRP1, YWHAE, SDCBP, and SOX10 in melanocytes and particularly YWHAE in keratinocytes. This suggests that upregulated miR-155 is associated with vitiligo pathogenesis which alters interferon signaling as well as targets melanogenesis-associated genes (Šahmatova et al. 2016). Both miR-155 and miR-10a are on their own largely dispensable for regulatory T cell (Treg) function and stability which is responsible for suppressing autoimmune pathology. Inhibition of miR-10a expression leads to reduced FOXP3 expression levels which subsequently decrease the stability of Treg cells (Jeker et al. 2012) resulting in the insufficient suppression of inflammation in autoimmune diseases which could likely happen in vitiligo patients.

Identifying the susceptible genes and their variants which drive the way to the onset of disease is fundamental to unravel their contribution in disease progression. Most of the vitiligo-associated genes are plausible biological candidate genes which when altered are responsible for stimulating melanocyte-specific immune responses. These candidate genes encode immunoregulatory and melanocyte proteins constituting a dense immunoregulatory network, an alteration in which highlights the systems and pathways mediating vitiligo susceptibility (Spritz 2013). Genetic variation alters or damages protein structure disrupting protein–protein interactions which are otherwise essential for regulatory processes constituting the pretext of disease development. Network-based studies of these interacting proteins may impart an insight into disease pathogenesis initiating better diagnosis and the feasibility of personalized treatment for vitiligo patients in the future. Our PPI network analysis results identified 15 hub proteins, namely, IL10, IFNG, IL4, CD44, IL1B, CTLA4, GZMB, FOXP3, TNF, IL2RA, CAT, ESR1, TLR2, HLA-A, and GSTP1 to be associated with vitiligo. The balance between pro- and anti-inflammatory cytokines plays a significant role in the pathogenesis of vitiligo. Higher concentrations of IFNG, a pro-inflammatory cytokine, enhance T cell-melanocyte attachment in the skin initiating T cell-mediated apoptosis of melanocytes in vitiligo. On the other hand, reduced concentrations of IL-10, a potent regulator of anti-inflammatory immune responses, were observed in vitiligo patients (Singh et al. 2012). An increased IL-10 concentration with an increase in the IFNG levels exhibited a positive correlation with disease duration as reported in vitiligo patients (Ala et al. 2015). IL-4, an immunomodulatory cytokine, stimulates B cell proliferation and T cell development that leads to the elevation of baseline IgE levels inducing inflammation (Del Prete et al. 1988). Polymorphisms in the IL4 gene are known to increase its expression increasing the lgE levels thereby implicating its role in autoimmunity-mediated vitiligo susceptibility (Imran et al. 2012). CTLA4 expressed by Tregs is a negative regulator of T cell function and fosters tolerance to self-antigens. Decreased levels of CTLA4 mRNA and deregulated CTLA4 expression due to genetic variations have been found in vitiligo patients (Dwivedi et al. 2011) suggesting its involvement in susceptibility to vitiligo. Also, upregulated CD44 expression in response to naive T cell proliferation as induced by autoimmune melanocyte destruction concomitantly increases T cell development implicating the complex regulation of self-reactive T cells in vitiligo (Byrne et al. 2014).

Increased mRNA levels of IL1B increases SOD levels leading to increased H₂O₂ production as observed in vitiligo patients. Genetic variability in IL1B resulting in altered IL1B transcript levels might be associated with elevated NPY levels in patients with vitiligo whose synthesis is governed by IL1B (Laddha et al. 2014). Increased NPY levels lead to epidermal and dermal hypoxia which might potentiate melanocyte death in vitiligo (Tu et al. 2001). Alterations in CAT have been reported to result in the reduction of the catalase enzyme activity and consequently evoke excess H₂O₂ accumulation in the entire epidermis of vitiligo patients (Casp et al. 2002). Although the genetic mechanisms of estrogen in increased pigment cell activity are not largely known yet, ESR1 expression on human melanocytes has been demonstrated to have specific actions in human pigmentation (Im et al. 2002). Also, genetic variation in ESR1 gene has been reported to show its association with vitiligo (Jin et al. 2004). Additionally, GSTP1 is broadly expressed in defense against oxidative stress wherein they detoxify a variety of electrophilic compounds generated by ROS-induced damaged cells (Nebert and Vasiliou 2004). Altered GSTP1 expression fails to protect cells against chemical toxicity and stress contributing to melanocyte death in vitiligo patients (Dušinská et al. 2001; Liu et al. 2009).

Effector functions of IL2RA and GZMB in the target cell killing by cytotoxic T cells (CTLs) and NK cells activation-induced cell death terminate immune responses and mediate melanocyte killing in vitiligo (Spritz 2010). GZMB also have a role to play in cleaving melanocyte proteins that constitute vitiligo auto-antigens activating auto-antigens that initiate and propagate autoimmunity directed against melanocytes (Darrah and Rosen 2010). FOXP3, the master regulator of Treg cells, have a vital role in maintaining immune balance and its alteration triggers autoimmune diseases including vitiligo (Jahan et al. 2015). TNF downregulates MITF affecting melanocyte development and proliferation, and ultimately affecting melanogenesis. Also, TNF-α downregulates MSHR binding activity and reduces MC1-R expression, both of which are known inducers of melanogenesis (Camara-Lemarroy and Slas-Alanis 2013). TRAIL, a TNF family death receptor, activates caspases and cleaves melanocyte proteins and also promotes dendritic cell-mediated melanocyte death eliciting apoptosis of primary human melanocytes (Larribere et al. 2004). TNF, thus, acts as the central regulatory effector in the immunopathological mechanisms involved in vitiligo.

Among TLRs, TLR2 is fundamental for immune responses against mycobacterial infections, in sensing oxidative stress and cellular necrosis and, also in inducing apoptosis (Petry and Gaspari 2009). It also has the propensity to recognize a wide array of antigens evincing its instrumental role in the evolution of self-reactive diseases (Borrello et al. 2011). Altered expression and signaling due to TLR2 polymorphism have been proposed to be the reason for inadequate immune responses to viral or pathogenic antigens in vitiligo (Karaca et al. 2013; Traks et al. 2015).

The functional enrichment analysis result of the single large cluster consisting of 64 proteins demonstrated the vitiligo-associated genes to be primarily involved in immune response regulation by cytokines and oxidative stress, and apoptotic processes. Oxidative stress in the melanocytes stimulates local inflammatory responses whereby it leads to the activation of innate immune processes as a result of which melanocyte-specific cytotoxic immune responses are evoked in vitiligo patients. Vitiliginous melanocytes show increased surface expression of HLA-A, a class I MHC receptor, which enables it to present multiple autoantigens to T cells destroying skin melanocytes (Hayashi et al. 2016). Also, increased expression levels of HLA class II molecules triggers an increased production of immunostimulatory cytokines that may act as an adjuvant during the presentation of autoantigens (Cavalli et al. 2016), tying together with HLA class I molecules in the development of autoimmunity in vitiligo patients. Also, alteration in the concentration of various pro-inflammatory and anti-inflammatory cytokines such as IL-10, IL-2, TNF, and IFNG has been associated with many autoimmune disorders (Singh et al. 2012).

Apart from exhibiting pathways associated with immune responses and autoimmunity, our pathway analysis results also consisted of pathways corresponding to infectious diseases, particularly viral infections. Several studies have implicated the etiopathogenesis of vitiligo to multiple viral infections as epidermal melanocytes are important targets of viruses (Duvic et al. 1987; Grimes et al. 1996). Also, viral infectious diseases, in most cases, impair the body’s systemic immune response. This explains the reason why the pathways associated with infectious diseases were also observed to be significant in the results. In the protein–pathway network, TNF was found to be involved in 24 out of 30 significant pathways indicating it to be a prime regulator of vitiligo immunopathogenesis (Birol et al. 2006).

Our SNP analysis predicted rs1801133 (MTHFR), rs5743708 (TLR2), and rs11575993 (SOD2) to be functionally significant for vitiligo susceptibility. Previous studies have shown these mutations to be associated with other diseases as well. The A222V variant of MTHFR has been reported to be associated with vascular diseases (Morita et al. 1997) and cancer (Hubner et al. 2007) The R753Q variant of TLR2 has been reported to be associated with colorectal cancer (Boraska Jelavić et al. 2006), atherosclerosis (Bielinski et al. 2011), and tuberculosis (Ogus et al. 2004). However, the L84F variant of SOD2 has only been reported in vitiligo till date (Laddha et al. 2013).

There is no definite cure available for vitiligo, and hence, the various treatment options available merely aim to improve skin appearance by repigmentation or stabilizing depigmentation without the assurance of reoccurrence or extension of depigmentation (Njoo and Westerhof 1999). Our drug–target network analysis revealed novel potential drug candidates which could be explored for improved therapeutics for vitiligo. It is noteworthy that most of the drugs defined for the protein targets, even a few hub protein targets, are not indicated to be used for the treatment of vitiligo as reported in DrugBank. This indicates a pressing need to evaluate these drugs and perform investigational studies to identify new indications and elucidate the efficiency of these drugs for the treatment of vitiligo. This would lead to significant contributions in drug discovery complementing the existing drug pipelines, thereby, improving the quality of life in vitiligo patients. Additionally, some of the hub proteins are targeted by a single drug which can be further examined to contrive better effective drugs to enhance the success rate of treatments. Interacting drug partners of a particular drug targeting a particular protein might either directly target that protein or enhance the efficacy of the drug. With this concept, we constructed the drug–similarity network and found interacting drugs for those proteins which were previously shown to be targeted by a single drug in the drug–target network. These interacting drugs might function as an alternative to the native drug with enhanced efficacy.

Polymorphism in TLR2 (rs5743708) was found to be deleterious in our SNP analysis results indicating its potentiality to induce vitiligo. TLR2 was also found to be one of the hub protein targets for which no drug information is available yet. The R753Q mutation was identified to be located within the TIR domain, an intracellular signaling domain, which compromises the signaling capacity of TIR domain impairing MyD88–TLR2 assembly. This inactivates NF-κB signaling pathway which can invariably influence the regulation of inflammatory processes and can even impair melanogenesis suggesting its role in vitiligo pathogenesis (Karaca et al. 2013; Traks et al. 2015). Our PPI network of TLR2 shows its interacting proteins suggesting that an altered TLR2 might have an impact on its interaction with other proteins essential for many biological functions and signaling processes. This indicates the need to analyze the structural details of the protein and the effect of mutation on its structure and function, and carry out further experimental studies to discover new drugs targeting the mutated TLR2 protein associated with vitiligo.

Conclusions

Investigating the patterns of miRNAs and its associations with diseases and identifying the susceptible genes variants triggering disease onset is significant for unraveling the key factors that underpin disease manifestation. We, therefore, identified vitiligo-associated miRNAs and their targets, and susceptible genes, and carried out a comprehensive network analysis of these data which revealed the association of significant hub miRNAs and proteins with disease susceptibility. We validated their functional role and carried out SNP analysis wherein we identified mutation in TLR2 (R753Q) as deleterious. Our drug–target network and drug–similarity network unveiled novel molecular determinants and drug repositioning candidates for vitiligo. Both our databases DermaGene and VitmiRS provide authentic and relevant information on the essentialities of disease association which may facilitate researchers in identifying dermatological disorder causation. We intend to regularly update both the databases, DermaGene and VitmiRS, respectively, through a bi-annual review to ensure consistency of data and to ensure it is maintained up-to-date. The procedure for regular update would rely on identification of latest entries which may be missing in our database and merging the new data with the existing data. Our approach can provide an insight into the regulatory mechanisms of disease manifestation, thereby, implicating its role in improved therapeutic and diagnostic interventions.