Journal Logo

Research Article: Observational Study

Identification of key genes and pathways involved in vitiligo development based on integrated analysis

Lei, Zixian MDa; Yu, Shirong MDb; Ding, Yuan MDb; Liang, Junqin MDb; Halifu, Yilinuer MDb; Xiang, Fang MDb; Zhang, Dezhi MDb; Wang, Hongjuan MDb; Hu, Wen MDb; Li, Tingting MDb; Wang, Yunying MDb; Zou, Xuelian MDb; Zhang, Kunjie MDb; Kang, Xiaojing MD, PhDb,∗

Editor(s): Garcovich., Simone

Author Information
doi: 10.1097/MD.0000000000021297


1 Introduction

Vitiligo is a common acquired, chronic skin condition with the loss of functioning melanocytes, affecting skin, hair or mucosa, or all, characterized by milky white, nonscaly spots with distinct margins.[1] These clinical symptoms seriously affect the mental health and quality of life of individuals.[2] With an estimated prevalence of 0.5% to 1% in the world population,[3] vitiligo is regarded as one of the most important skin conditions causing prejudice and segregation in certain cultures. Both sexes in adults and children seem to be affected by vitiligo equally.[4] However, the etiology and pathogenesis of vitiligo is still obscure and both multifactorial and polygenic. Selective destruction of the pigment-producing cell, namely melanocytes, is core to vitiligo.[5] Autoantibodies against melanocytes were shown to exist in patient serum, and lately have been considered as predictors of disease progression.[6] In addition, several autoimmune comorbidities have been reported in association with vitiligo such as autoimmune thyroid disease,[7] systemic lupus erythematosus,[8] Addison disease[8] and psoriasis,[9] etc. Furthermore, a number of susceptibility loci for vitiligo were detected through genome-wide association studies, lots of which are shared with other autoimmune diseases.[1,7,10]

To regain and maintain the pigment, current therapeutic options for vitiligo are mainly based on 3 aspects: modulating the autoimmune response; reducing oxidative stress in melanocytes; activating regeneration of melanocytes.[11] However, affected individuals require long-term treatment and poor efficacy. So far, vitiligo remains a disorder with restricted therapeutic options and difficult to control its progression effectively for patients, which underlines the unmet need to define the molecular mechanism underlying this disorder and to clear the therapeutic target for new safe and effective therapies to manage this skin condition. Therefore, more intensive research is needed to discover the molecular mechanism, detect novel promising therapeutic targets, and develop efficient therapeutic approaches for vitiligo.

Currently, insights into the pathogenesis of autoimmune diseases at the system level can be provided by researches of systems biology, which can facilitate biological studies. Some datasets of genome-wide gene expression in patients with vitiligo are freely accessible in the public databases, viz. Gene Expression Omnibus.[12,13] Much valuable information in those datasets can be mined and reused by some tools of systems biology. Consequently, weighted gene coexpression network analysis (WGCNA) was performed to categorize the genes expressed aberrantly and significantly into several functional modules,[14,15] expecting to deepen insights into the pathogenesis of vitiligo and seek some novel promising drug targets for the therapy of vitiligo.

2 Materials and methods

2.1 Data resources and principal component analysis (PCA)

The vitiligo datasets GSE65127 and GSE75819 were collected from the NCBI Gene Expression Omnibus (GEO) database ([16] GSE65127 includes a total of 40 samples obtained from 10 patients and 10 healthy volunteers (1 biopsy in matched anatomical areas), comprising tissues from 10 Lesional skins, 10 Peri-Lesional skins, 10 Non-Lesional skins, and 10 healthy volunteers’ skins. GSE75819 contains tissues from 15 Lesional skins and 15 Non-Lesional skins. Series matrix files and clinical information tables were downloaded from the GEO website. PCA is an unsupervised method and data reduction technology, which allows the analysis of major sources of change in multidimensional datasets without introducing inherent deviations. Hence, PCA were performed on disease and health samples to condense the information down to a reduced number of linearly independent variables.

2.2 Differential expression analysis

The “sva” R package was applied to conduct batch normalization of the expression data from the 3 different datasets. Then, the normalized expression matrixes of Non-Lesional, Peri-Lesional, and Lesional differences were analyzed by “limma” R package respectively.[17] The differentially expressed genes (DEGs) were acquired by setting threshold P < .05.

2.3 Coexpression analysis

WGCNA is a method to establish scale-free gene coexpression network.[14] A weighted gene coexpression network of 3 groups of DEGs was constructed using “WGCNA” R package. The soft threshold power of β was set as 9, and weighted adjacency matrix was subsequently generated. Then, hierarchical clustering analysis was carried out. And a weighted adjacency matrix was formed. In addition, the weighted adjacency matrix was transformed into a topological overlap matrix to estimate its connectivity in the network. To evaluate the association between gene coexpression modules and traits, “nlme” R package was adopted to establish a linear mixed model. For each WGCNA module, the gene with kme Pearson correlation (> 0.90)[18] was regarded as a hub gene. Afterward, the receiver operating characteristic (ROC) curves of hub genes were analyzed by “pROC” R package.[19]

2.4 GO function and KEGG pathway enrichment analysis and gene set enrichment analysis (GSEA)

To clarify the possible biological roles of these genes in coexpression networks, the cluster Profiler R package[20] was performed to create gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment paths, results, and plots. P < .05 indicated GO and KEGG pathway with significant difference. Thereafter, GSEA[21] was employed according to the application available at the Broad Institute Gene Set Enrichment Analysis website. For purpose of testing the KEGG pathways gene sets, the activation or inhibition of KEGG pathway was also evaluated by GSVA algorithm.[22]

2.5 Pivot analysis and drug prediction

Pivot is defined as a regulator, including transcription factors, which has a significant regulatory effect on module genes. To explore the pivotal regulators of vitiligo-related module genes, all human transcription factor target data in TRRUST V2 database[23] was utilized as background set for pivot analysis. Setting threshold P < .05. Prediction of drug and drug targets for the vitiligo dysfunction module were also carried out. On the basis of drug target data from Drugbank database,[24] drugs with significant regulatory modules were screened according to hypergeometric test. Setting threshold P < .05.

2.6 RNA expression analysis by quantitative reverse transcription polymerase chain reaction (qRT-PCR)

The whole blood samples of 3 vitiligo patients and 3 healthy controls were collected respectively. Written informed consent was obtained from each subject, and the research protocol was approved by the ethics committee of People's Hospital of Xinjiang Uygur Autonomous Region. The total RNA of whole blood samples were extracted using Trizol reagent (Invitrogen, USA) by following the manufacturer's protocol. Following DNase treatment, RNA was reverse transcribed into complementary DNA with a complementary DNA synthesis kit (Takara, Dalian, China), according to the manufacturer's protocol. qPCR was performed by using SYBR Select Master Mix (Applied Biosystems, USA) and primers (Table 1). The expression of target gene was normalized to GAPDH levels in respective samples as an internal control and calculated using the 2−ΔΔCT method. Students t test was used to compare the differences between 2 groups. P < .05 was considered statistically significant.

Table 1
Table 1:
Primer sequence.

3 Results

3.1 Coexpression of differentially expressed genes in vitiligo

First, PCA was used to analyze the degree of dispersion of different samples. It was found that the 3 groups of disease samples were closer, while the healthy control group samples were farther (Supplemental Digital Content (Figure S1, Then, DEGs were screened between Non-Lesional, Peri-Lesional, Lesional, and healthy skin tissues, respectively, which were considered to be associated with vitiligo. 3132 DEGs were identified in Non-Lesional, 3610 in Peri-Lesional, and 3083 in Lesional (Fig. 1A). Next, WGCNA was performed on the 3 groups of DEGs and a gene coexpression network was constructed. By eliminating the gray module without cooperative expression behavior, 8 modules were obtained, in which the genes had the behavior of coexpression (Fig. 1B–D). Through the correlation study between module and phenotype, it was noteworthy that the highest correlation with Lesional, Peri-Lesional, and Non-Lesional was module brown, blue, and green respectively (Fig. 1E). In addition, gene with the highest connectivity, named hub gene, were identified in each module (Table 2). Hub gene plays a significant role in the module and may function as a key part in the disease. As well, by ROC curve, the ability of hub genes to determine the progression of vitiligo was evaluated (Fig. 2A–H). The Area Under Curve value of NUF2 and HSPB6 increased gradually from Non-Lesional to Lesional, indicating that the more severe the degree of vitiligo, the stronger the evaluation ability of NUF2 and HSPB6. Surprisingly, the diagnostic ability of 7 hub genes was verified in an independent gene expression dataset GSE75819 (Fig. 2I, J). These results implied that most of the DEGs were coexpressed and interacted with each other in the occurrence and development of vitiligo.

Figure 1
Figure 1:
Coexpression of differentially expressed genes in vitiligo was observed. A, The differentially expressed genes between the skin tissue of patients with vitiligo in Lesional, Peri-Lesional, Non-Lesional, and the healthy controls. B, Soft threshold power analysis was used to obtain the scale-free fit index of network topology. C, Hierarchical cluster analysis was conducted to detect coexpression clusters with corresponding color assignments. Different colors represent different modules. D, Heatmap depicts the Topological Overlap Matrix (TOM) of genes selected for weighted coexpression network analysis. Light color represents lower overlap and red represents higher overlap. E, The correlation between module and phenotype. Red reflects positive correlation, blue reflects negative correlation.
Table 2
Table 2:
Hub genes in modules of differentially expressed genes in vitiligo.
Figure 2
Figure 2:
ROC curves of hub genes. A–H, Calculation of Area Under Curve value of hub genes in the modules and ROC curve drawing. Different colors represent different groups. I, J, The seven hub genes validated in dataset GSE75819 and corresponding Area Under Curve values. ROC = receiver operating characteristic.

3.2 Molecular dysregulation mechanism in vitiligo

These module genes may play different roles in regulating the occurrence and development of vitiligo. To identify the biological functions of dysfunctional genes, GO and KEGG enrichment analysis of module genes were performed. Abundant GO terms were collected, including 562 cell composition terms, 842 molecular functional terms, as well as 3916 biological processes terms (Supplemental Digital Content (Table S1, Among them, dysfunctional genes were implicated in macrophage activation involved in immune response, T-helper 17 type immune response, and other related biological functions (Fig. 3A). On the other hand, enrichment analysis of KEGG pathway revealed 157 terms (Supplemental Digital Content (Table S1,, including p53 signaling pathway, chronic myeloid leukemia, etc. (Fig. 3B).These signaling pathways played a vital role in the process of vitiligo. Based on GSVA, dysregulation of these signaling pathways was closely related to vitiligo progression (Fig. 3C). Moreover, the signal pathways dysregulated in vitiligo were also explored through GSEA. The results showed that 6 signal pathways were activated incessantly from Non-Lesional to Peri-Lesional and then to Lesional (Fig. 3D–F). In the dataset GSE65127, 4 of the dysregulated signal pathways were verified by an independent dataset GSE75819 (Fig. 3G, H). These results indicated that the signal pathways such as VEGF and mTOR were continuously dysregulated in the pathogenesis of vitiligo.

Figure 3
Figure 3:
The enrichment analysis of signal pathways of dysregulated gene in vitiligo. A, The enrichment of biological function of modular dysregulated genes. B, The enrichment of signal pathways of modular genes in KEGG. From blue to red, the significance of the nodes gradually increases. C, Thermogram shows the activation state of signal pathways of KEGG in different groups after processing by GSVA. The red node represents upregulation and the blue node represents downregulation. The common signaling pathways in GSEA for Non-Lesional (D), Peri-Lesional (E), and Lesional (F). G, The signaling pathways in GSEA for Lesional in dataset GSE65127. H, The identical signaling pathways of Lesional in dataset GSE65127 and GSE75819. GSEA = gene set enrichment analysis, GSVA = gene set variation analysis, KEGG = Kyoto Encyclopedia of Genes and Genomes.

3.3 Persistent maladjusted genes in the development of vitiligo

In case there was a persistent dysregulation in the signaling pathway, there may be genes with persistent disorders. Consequently, 1197 genes were identified in the 3 groups of DEGs (Fig. 4A). Among them, 269 genes were upregulated constantly from Non-Lesional to Peri-Lesional and then to Lesional, while 82 genes were downregulated (Supplemental Digital Content (Table S2, Hence, we held that these persistent dysregulated genes may be involved in the occurrence and development of vitiligo. Notably, the expression pattern of persistent maladjusted genes in the modules was displayed by thermogram (Fig. 4B). It was obvious that most of the persistent maladjusted genes were clustered in the module blue. Remarkably, HSPB6 was also a gene that continued to be downregulated in this study. In addition, the continuous dysregulated genes in the module had a tendency of up- or down-regulation from Non-Lesional to Lesional. Significantly, the persistent maladjusted genes were found to participate in T-cell receptor signaling pathway, etc. These signaling pathways functioned as a crucial part in the development of vitiligo. The results revealed that there was a constant alteration of gene expression in the course of vitiligo, which may be related to the severity of vitiligo.

Figure 4
Figure 4:
Persistent maladjusted genes in 3 groups of differentially expressed genes. A, Venn diagram of the 3 groups of differentially expressed genes, DEG1 represents Non-Lesional, DEG2 represents Peri-Lesional, and DEG3 represents Lesional. B, Thermogram displays the expression of persistent maladjusted genes in the modules. Red node represents upregulated gene, blue node represents downregulated gene.

3.4 Regulation network of persistent maladjusted genes in vitiligo

Transcription regulation indicates that the level of gene expression alters with the alteration of the transcription rate, playing an essential role in transmission of genetic information accurately and diversely.[25] To this end, transcription regulation of module genes was explored, which showed that 42 transcription factors had significant regulatory effect on the module genes (Supplemental Digital Content (Table S3, By screening the regulators of persistent maladjusted genes, it was found that specificity protein 1 (Sp1) was involved in the development of vitiligo by targeting ITGAV to regulate the thyroid hormone signaling pathway. Alternatively, based on drug prediction, 228 drugs were spotted which may have therapeutic effects on module genes (Supplemental Digital Content (Table S4, Subsequently, focusing on persistent dysregulated genes, transcription factors and drugs were extracted and a map of regulatory network was delineated (Fig. 5A). Furthermore, as expected, the expression of these key genes was verified in GSE75819 dataset (Fig. 5B). The outcomes deciphered that the persistent dysregulated genes affecting the progression of vitiligo were modulated by a number of transcription factors and drugs. Besides, the expressions of hub genes (CACTIN, DCTN1, GPR143, MRPL47, NKTR) in vitiligo patients were verified by qRT-PCR, which were consistent with our findings (Fig. 5C).

Figure 5
Figure 5:
Regulation and verification of key persistent maladjusted genes in vitiligo modules. A, The transcription factors and drugs that regulate the persistent maladjusted genes in the modules, as well as the signal pathway network participated by the target genes. Arrows represent transcription factors, diamonds represent drugs, circles represent target genes, triangles represent modules, and squares represent gene-enriched signaling pathways. B, Data in GSE75819 confirmed the differential expression of key persistent maladjusted genes in vitiligo. The red box represents downregulation while the blue box represents upregulation. C, The qRT-PCR results of hub genes in blood samples of vitiligo patients compared with healthy controls,∗ Indicates P < .05; ∗∗∗ Indicates P < .001.

4 Discussion

The module analysis methods and other related data were applied to study the potential molecular linkages in different process stages of occurrence and development of vitiligo, that was from Non-Lesional, to Peri-Lesional then to Lesional. Under these circumstances, the data from Lesional tissues, Peri-Lesional tissues, and Non-Lesional tissues compared with normal tissues respectively could represent, to a certain extent, the different process stages of occurrence and development of vitligo.

Accordingly, based on WGCNA, the correlation between module and phenotype during the course of vitiligo was identified, including module brown, module blue, and module green, respectively. The module brown was mainly enriched in response to oxygen levels and Wnt pathway. Studies have shown that oxidative stress does play a vital role in causing subsequent autoimmune responses associated with vitiligo.[26,27] Lately, Tang et al[28] demonstrated Vitamin D protected human melanocytes against oxidative damage by activation of Wnt/β-catenin signaling. The module blue was mostly enriched in autophagy and immune response. The damage-associated molecular patterns or autoantigens were generated by stressed melanocytes, which initiated innate immunity and adaptive immunity, activating an inflammatory cascade leading to the dysfunction and death of melanocytes.[29] The vulnerability of melanocytes to oxidative stress may be induced by dysfunctional autophagy.[27] Furthermore, the module green was largely enriched in divalent inorganic cation transport and calcium ion transport. Calreticulin was an omnipresent endoplasmic reticulum protein that regulated translocation of intracellular calcium ions from endoplasmic reticulum lumen to the melanocytes surface under oxidative stress.[30] Besides, calreticulin also prompted pro-inflammatory cytokines relating with immune responses, including interleukin-6 and tumor necrosis factor-α.[30]

Remarkably, for the current study, persistent maladjusted genes from Non-Lesional to Lesional were detected, which may indicate the occurrence and development of vitiligo. Enrichment analysis showed that these genes were significantly enriched in T-cell receptor signaling pathway. Raam et al described in detail that vitiligo was a T-cell-mediated secondary autoimmune disease after treatment of alemtuzumab.[31]

Combined with gene difference analysis and network topology, potential hub genes were identified. Up to present, the relationships between the 8 hub genes identified in this study and vitiligo have not been reported in the literature. They may be potential biomarkers and therapeutic targets for vitiligo in the future. The hub gene of module brown was NUF2 (NUF2 component of NDC80 kinetochore complex), also known as: cell division cycle associated 1. In line with previous findings, Thang et al[32] discovered that knockdown of cell division cycle associated 1 in oral cavity carcinoma cells led to impaired growth and apoptosis in cells. Unfortunately, in the blood samples of vitiligo patients, the expression of NUF2 detected by qRT-PCR was downregulated compared with normal blood samples, which was inconsistent with our findings. However, the expression of NUF2 was verified in an independent dataset GSE75819, in which the sample type was skin tissue. The hub gene of module blue was MRPL47 (mitochondrial ribosomal protein L47). In eukaryotic cells, mitochondria have their own ribosomes (mitoribosomes), which synthesize a small number of proteins, all indispensable for the biogenesis of the oxidative phosphorylation system.[33] The hub gene of module green was GPR143 (G protein-coupled receptor 143), which was targeted to melanosomes in pigment cells. It was reported that mutations of GPR143 were connected with X-linked ocular albinism type 1.[34] Although skin pigmentation was slightly reduced, subcellular abnormalities such as occurrence of giant melanosomes (macromelanosomes), decrease in melanosome number, and changes of melanosome motility were still observed in epidermal melanocytes.[35] The expressions of MRPL47 and GPR143 were confirmed both in the dataset GSE75819 and in the blood samples of vitiligo patients detected by qRT-PCR, which were consistent with our analysis.

Importantly, HSPB6 (heat shock protein family B (small) member 6), as the hub gene of module yellow, was also a persistent maladjusted gene. With mounting evidence accumulated, in response to varied cellular stress or damage, upregulated HSPB6 defended cells against otherwise lethal conditions, such as eliciting smooth muscle relaxation, restraining platelet aggregation, and playing a key regulatory role in cardioprotection and neuroprotection.[36] Although the expression level of HSPB6 proteins in normal skin tissues was low, it played a significant role in the detection data in the current study, suggesting that HSPB6 was also an important gene in the pathogenesis of vitiligo. Besides, the ROC curve displayed that HSPB6 as well as other hub genes had potential diagnostic value in dataset GSE65127. However, compared with normal blood samples, the expression of HSPB6 in the blood samples of vitiligo patients was decreased, but there was no statistical difference. By expanding the sample size, it may provide latent biomarker for future research.

In addition to crucial genes, there were also dysregulated signaling pathways that played a vital role in the course of vitiligo, including p53 signaling pathway. Birlea et al emphasized the active participation involved in defined molecular pathways that modulated the balance of melanocytes and keratinocytes between stemness and differentiation states, including p53 and its downstream effectors regulating melanogenesis and Wnt/β-catenin participating in proliferative, migratory, and differentiation in different pigmentation systems.[37] Lee et al[38] provided that Interferon-inducible T-cell alpha chemoattractant (ITAC) induced the melanocytic migration and hypopigmentation through destabilizing p53 via histone deacetylase 5. Generally, these results were consistent with our integrated analysis.

Moreover, the results presented that transcription factors dominated by SP1 had noteworthy regulatory effects on persistent dysfunctional genes. SP1 is significant to the TATA-less genes, which modulates transcription of manifold target genes related to cell growth, differentiation, apoptosis, and immune response.[39] Previous researches unveiled that SP1 controlled the maximal constitutive activity and basal expression of Ah receptor (AHR) gene by binding to the GC-rich motifs.[40]AHR-129 C > T variant could lead to allele-specific binding of SP1 to AHR promoter and further modify its transcription and downstream effectors in vitiligo.[41] Based on persistent maladjusted genes, drug prediction and drug target information suggested that Celecoxib had important regulatory effects on dysfunction modules. Celecoxib is the first cyclo-oxygenase 2-selective inhibitor, which has been approved for treating inflammation and pain for a long time.[42] Several drugs required further investigation for the therapy or side effects.

In the present study, WGCNA was performed to analyze the DEGs in vitiligo for the first time. In addition, the data acquired from Lesional skins, Peri-Lesional skins, and Non-Lesional skins compared respectively with normal skins may represent the different process stages of occurrence and development of vitligo. However, the main limitation of the present study was the small sample size and there was little literature to support our findings. Regrettably, at this moment, due to the unavailability of tissue samples, it was difficult for us to detect the protein expression of each hub gene. Therefore, we collected blood samples of vitiligo patients for verification of hub genes. However, heterogeneity between different sample types and individuals may lead to inconsistent results, and future larger cohort validation is needed. Taken together, the hub genes and dysregulated signaling pathways as well as transcription factors and drugs may be verified as potential molecular targets and their roles in vitiligo need to elucidate by further studies.

5 Conclusions

In general, based on WGCNA and a series of comprehensive analyses, functional modules were identified and persistent maladjusted genes in the course of vitiligo were also detected. Hub genes of each module were uncovered, which may develop to be diagnostic biomarker or therapeutic target of vitiligo. At least, the present study might provide an integrated and in-depth insight for exploring the underlying mechanism of vitiligo occurrence and development and predicting potential treatment methods and mechanism as well.

Author contributions

Conceptualization: Xiaojing Kang, Zixian Lei, Shirong Yu, Yuan Ding, Junqin Liang.

Methodology: Zixian Lei, Shirong Yu, Yilinuer Halifu.

Software: Zixian Lei, Shirong Yu.

Validation: Xiaojing Kang, Zixian Lei, Shirong Yu, Yuan Ding, Junqin Liang.

Formal analysis: Zixian Lei, Shirong Yu, Yuan Ding, Junqin Liang, Yilinuer Halifu, Fang Xiang, Dezhi Zhang.

Investigation: Zixian Lei, Shirong Yu, Yuan Ding, Yilinuer Halifu, Fang Xiang, Dezhi Zhang.

Resources: Zixian Lei, Hongjuan Wang, Wen Hu, Tingting Li, Yunying Wang, Xuelian Zou, Kunjie Zhang.

Data Curation: Zixian Lei, Yilinuer Halifu, Fang Xiang, Dezhi Zhang, Hongjuan Wang, Wen Hu, Tingting Li.

Writing - Original Draft: Zixian Lei, Shirong Yu.

Writing - Review & Editing: Xiaojing Kang, Zixian Lei.

Visualization: Zixian Lei, Yunying Wang, Xuelian Zou, Kunjie Zhang.

Supervision: Xiaojing Kang, Yuan Ding, Junqin Liang.

Project administration: Xiaojing Kang.

Funding acquisition: Xiaojing Kang.

All authors agreed to be accountable for all aspects of this work.


[1]. Ezzedine K, Eleftheriadou V, Whitton M, et al. Vitiligo. Lancet 2015;386:74–84.
[2]. Elbuluk N, Ezzedine K. Quality of life, burden of disease, co-morbidities, and systemic effects in vitiligo patients. Dermatol Clin 2017;35:117–28.
[3]. Ezzedine K, Lim HW, Suzuki T, et al. Revised classification/nomenclature of vitiligo and related issues: the Vitiligo Global Issues Consensus Conference. Pigment Cell Melanoma Res 2012;25:E1–13.
[4]. Wang KY, Wang KH, Zhang ZP. Health-related quality of life and marital quality of vitiligo patients in China. J Eur Acad Dermatol Venereol 2011;25:429–35.
[5]. Riding RL, Harris JE. The role of memory CD8(+) T cells in vitiligo. J Immunol 2019;203:11–9.
[6]. Jimenez-Brito G, Garza-de-La-Pena E, Perez-Romano B, et al. Serum antibodies to melanocytes in patients with vitiligo are predictors of disease progression. Skinmed 2016;14:17–21.
[7]. Spritz RA. Shared genetic relationships underlying generalized vitiligo and autoimmune thyroid disease. Thyroid 2010;20:745–54.
[8]. Alkhateeb A, Fain PR, Thody A, et al. Epidemiology of vitiligo and associated autoimmune diseases in Caucasian probands and their families. Pigment Cell Res 2003;16:208–14.
[9]. Yen H, Chi CC. Association between psoriasis and vitiligo: a systematic review and meta-analysis. Am J Clin Dermatol 2019;20:31–40.
[10]. Spritz RA. Six decades of vitiligo genetics: genome-wide studies provide insights into autoimmune pathogenesis. J Invest Dermatol 2012;132:268–73.
[11]. Bleuel R, Eberlein B. Therapeutic management of vitiligo. J Dtsch Dermatol Ges 2018;16:1309–13.
[12]. Rung J, Brazma A. Reuse of public genome-wide gene expression data. Nat Rev Genet 2013;14:89–99.
[13]. Sparks R, Lau WW, Tsang JS. Expanding the immunology toolbox: embracing public-data reuse and crowdsourcing. Immunity 2016;45:1191–204.
[14]. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics 2008;9:559.
[15]. Prom-On S, Chanthaphan A, Chan JH, et al. Enhancing biological relevance of a weighted gene co-expression network for functional module identification. J Bioinform Comput Biol 2011;9:111–29.
[16]. Barrett T, Wilhite SE, Ledoux P, et al. NCBI GEO: archive for functional genomics data sets—update. Nucleic Acids Res 2013;41:D991–5.
[17]. Ritchie ME, Phipson B, Wu D, et al. Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 2015;43:e47.
[18]. Pripp AH. [Pearson's or Spearman's correlation coefficients]. Tidsskr Nor Laegeforen 2018;138:
[19]. Robin X, Turck N, Hainard A, et al. pROC: an open-source package for R and S+ to analyze and compare ROC curves. BMC Bioinformatics 2011;12:77.
[20]. Yu G, Wang LG, Han Y, et al. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 2012;16:284–7.
[21]. Subramanian A, Tamayo P, Mootha VK, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A 2005;102:15545–50.
[22]. Hanzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics 2013;14:7.
[23]. Han H, Cho JW, Lee S, et al. TRRUST v2: an expanded reference database of human and mouse transcriptional regulatory interactions. Nucleic Acids Res 2018;46:D380–6.
[24]. Wishart DS, Feunang YD, Guo AC, et al. DrugBank 5.0: a major update to the DrugBank database for 2018. Nucleic Acids Res 2018;46:D1074–82.
[25]. Levine M, Tjian R. Transcription regulation and animal diversity. Nature 2003;424:147–51.
[26]. Xie H, Zhou F, Liu L, et al. Vitiligo: how do oxidative stress-induced autoantigens trigger autoimmunity? J Dermatol Sc 2016;81:3–9.
[27]. He Y, Li S, Zhang W, et al. Dysregulated autophagy increased melanocyte sensitivity to H2O2-induced oxidative stress in vitiligo. Sci Rep 2017;7:42394.
[28]. Tang L, Fang W, Lin J, et al. Vitamin D protects human melanocytes against oxidative damage by activation of Wnt/beta-catenin signaling. Lab Investig A J Tech Methods Pathol 2018;98:1527–37.
[29]. Richmond JM, Frisoli ML, Harris JE. Innate immune mechanisms in vitiligo: danger from within. Curr Opin Immunol 2013;25:676–82.
[30]. Zhang Y, Liu L, Jin L, et al. Oxidative stress-induced calreticulin expression and translocation: new insights into the destruction of melanocytes. J Invest Dermatol 2014;134:183–91.
[31]. Ruck T, Pfeuffer S, Schulte-Mecklenbeck A, et al. Vitiligo after alemtuzumab treatment: secondary autoimmunity is not all about B cells. Neurology 2018;91:e2233–7.
[32]. Thang PM, Takano A, Yoshitake Y, et al. Cell division cycle associated 1 as a novel prognostic biomarker and therapeutic target for oral cancer. Int J Oncol 2016;49:1385–93.
[33]. De Silva D, Tu YT, Amunts A, et al. Mitochondrial ribosome assembly in health and disease. Cell Cycle 2015;14:2226–50.
[34]. d’Addio M, Pizzigoni A, Bassi MT, et al. Defective intracellular transport and processing of OA1 is a major cause of ocular albinism type 1. Hum Mol Genet 2000;9:3011–8.
[35]. De Filippo E, Manga P, Schiedel AC. Identification of novel g protein-coupled receptor 143 ligands as pharmacologic tools for investigating X-linked ocular albinism. Invest Ophthalmol Vis Sci 2017;58:3118–26.
[36]. Li F, Xiao H, Zhou F, et al. Study of HSPB6: insights into the properties of the multifunctional protective agent. Cell Physiol Biochem 2017;44:314–32.
[37]. Birlea SA, Costin GE, Roop DR, et al. Trends in regenerative medicine: repigmentation in vitiligo through melanocyte stem cell mobilization. Med Res Rev 2017;37:907–35.
[38]. Lee E, Choi SY, Bin BH, et al. Interferon-inducible T-cell alpha chemoattractant (ITAC) induces the melanocytic migration and hypopigmentation through destabilizing p53 via histone deacetylase 5: a possible role of ITAC in pigment-related disorders. Br J Dermatol 2017;176:127–37.
[39]. Tan NY, Khachigian LM. Sp1 phosphorylation and its regulation of gene transcription. Mol Cell Biol 2009;29:2483–8.
[40]. Racky J, Schmitz HJ, Kauffmann HM, et al. Single nucleotide polymorphism analysis and functional characterization of the human Ah receptor (AhR) gene promoter. Arch Biochem Biophys 2004;421:91–8.
[41]. Wang X, Li K, Liu L, et al. AHR promoter variant modulates its transcription and downstream effectors by allele-specific AHR-SP1 interaction functioning as a genetic marker for vitiligo. Sci Rep 2015;5:13542.
[42]. Frampton JE, Keating GM. Celecoxib: a review of its use in the management of arthritis and acute pain. Drugs 2007;67:2433–72.

bioinformatics; functional enrichment; hub gene; vitiligo; weighted gene coexpression network analysis

Supplemental Digital Content

Copyright © 2020 the Author(s). Published by Wolters Kluwer Health, Inc.