Novel Factors in the Pathogenesis of Psoriasis and Potential Drug Candidates Are Found with Systems Biology Approach

Psoriasis is a multifactorial inflammatory skin disease characterized by increased proliferation of keratinocytes, activation of immune cells and susceptibility to metabolic syndrome. Systems biology approach makes it possible to reveal novel important factors in the pathogenesis of the disease. Protein-protein, protein-DNA, merged (containing both protein-protein and protein-DNA interactions) and chemical-protein interaction networks were constructed consisting of differentially expressed genes (DEG) between lesional and non-lesional skin samples of psoriatic patients and/or the encoded proteins. DEGs were determined by microarray meta-analysis using MetaOMICS package. We used STRING for protein-protein, CisRED for protein-DNA and STITCH for chemical-protein interaction network construction. General network-, cluster- and motif-analysis were carried out in each network. Many DEG-coded proteins (CCNA2, FYN, PIK3R1, CTGF, F3) and transcription factors (AR, TFDP1, MEF2A, MECOM) were identified as central nodes, suggesting their potential role in psoriasis pathogenesis. CCNA2, TFDP1 and MECOM might play role in the hyperproliferation of keratinocytes, whereas FYN may be involved in the disturbed immunity in psoriasis. AR can be an important link between inflammation and insulin resistance, while MEF2A has role in insulin signaling. A controller sub-network was constructed from interlinked positive feedback loops that with the capability to maintain psoriatic lesional phenotype. Analysis of chemical-protein interaction networks detected 34 drugs with previously confirmed disease-modifying effects, 23 drugs with some experimental evidences, and 21 drugs with case reports suggesting their positive or negative effects. In addition, 99 unpublished drug candidates were also found, that might serve future treatments for psoriasis.


Introduction
Psoriasis is a multifactorial inflammatory skin disease. A recent systematic review reported a prevalence from 0% (Taiwan) to 2.1% (Italy) in children and from 0.91% (United States) to 8.5% (Norway) in adults. [1] Genetic predisposition and environmental factors are both important in disease etiology. Several genomewide association studies have been carried out and until now 36 susceptibility loci have been identified. [2] Environmental triggers are also reported such as drugs, smoking, mental stress, skin injury, Streptococcal infection, hormonal changes etc. [3] Psoriasis is an immune-mediated disease. Important immune cells and cytokines have been identified in disease pathogenesis such as IL6, IL17A, TNF etc. [4] Autoimmune basis for chronic inflammation is supposed, although no consistent antigen has been found. Patients with psoriasis have higher risk for metabolic syndrome, and risk increases with disease severity. Both diseases have immunological basis with common cytokines and genetic risk loci like CDKAL1. [5] Keratinocyte hyperproliferation is present in lesional phenotype and is responsible for scale formation. Keratinocyte differentiation markers like keratin 1 and keratin 10 are downregulated and parakeratosis (keratinocytes with nuclei in the stratum granulosum) is also present. [3] Psoriasis is one of the most studied skin diseases. By now more than 34000 hits are available in PubMed for the keyword ''psoriasis'' and the number is increasing. No spontaneous psoriasis-like skin disease is known in animals. Induced mouse models are available which are similar, but not the same as psoriasis in human. [6] Therefore drug discovery is difficult in such models what makes in silico analysis more essential. ''Omics'' data gives the opportunity to examine the disease with systems biology approach.
Stationary changes in gene expression are responsible for fixing phenotypes such as lesional skin areas in psoriasis. Several microarray studies have been carried out to characterize gene expression in healthy and psoriatic skin samples (Table 1). Microarray meta-analysis gives the opportunity to evade biological, regional, and study design-caused variation between studies. [7] Network analysis is a novel and highly developing area of systems biology. Considering gene expression data it is possible to explain alterations in intracellular processes with the analysis of protein-protein and protein-DNA (or gene regulatory) interaction networks. These networks consist of proteins and/or regulated genes as nodes and undirected or directed edges between them. Network centralities like degree or stress are suitable for ranking nodes. Total edge number belong to one node equals its degree in undirected networks. Nodes have in-and out-degrees based on edge directions in directed networks. Degree distribution follows a scale-free power law distribution in biological networks. This fact indicates that highly connected vertices have a large chance of occurring. Nodes with highest degree are called hubs and are essential in network stability. [8] Stress centrality indicates the number of shortest paths (from all shortest paths between any two nodes in the network) passing through the given node thus the capability of a protein for holding together communicating nodes. [9] Interconnecting nodes make up network motifs. Several, such as feed-forward or bifan motif are significantly enriched in biological networks compared to random networks. These elements have important role in network dynamics. [10] We hypothesized that it could be possible to find novel elements of psoriasis pathogenesis with detailed analysis of precisely constructed networks. Network motif enrichment caused by changes in gene expression could have important role in disease development and sustainment. It could be also possible to detect potential drug candidates by analyzing chemical-protein networks. Thus our goal was to construct reliable but yet detailed protein-protein, protein-DNA, merged (containing both proteinprotein and protein-DNA interactions) and chemical-protein interaction networks consisting of differentially expressed genes (DEG) between lesional and non-lesional skin samples and/or the coded proteins. Detailed analysis of these networks could help us to reveal novel players in disease pathomechanism and to identify network motifs and sub-networks with the ability to sustain lesional phenotype.

Microarray Meta-analysis
Six microarray studies examining lesional and non-lesional skin biopsy samples of psoriatic patients were found in Gene Expression Omnibus (GEO) ( Table 1). ''Minimum Information About a Microarray Experiment'' (MIAME) was available for each study. Only non-lesional and lesional samples from affected individuals were used for analysis, samples from healthy people were excluded. Raw.CEL files were downloaded and quality of each sample was assessed with the R package arrayQuality-Metrics. [11] This package defines sample quality with 5 different methods and generates plots for outlier detection. A sample was excluded if it was obviously an outlier in at least 1 measure or had borderline values in at least 2 measures (analysis results are in Dataset S1 compressed file; outliers and argument of exclusion is listed in Table S1). Raw data normalization of remaining samples was carried out with the R package Easy Microarray data Analysis (EMA). [12] GCRMA normalization method was used and probe sets with expression level below 3.5 were discarded. Probe set with the highest interquartile range (IQR) was chosen for common HUGO Gene Nomenclature Committee (HGNC) gene identifiers. Original findings were confirmed with published statistics. For this EMA was used after GCRMA normalization. More DEGs were found in some cases, which might be caused by the pre-filtering process with arrayQualityMetrics (Table S2). The R package MetaQC was used for filtering out low quality studies. [13] The fifty most prevalent gene set were chosen with the software Gene Set Enrichment Analysis (GSEA) and used for external quality control (EQC) score calculation. [14] GSEA was carried out for each study with the following settings: 1000 permutations; minimum set size was 5 and the gene set database was c2.all.4.0.symbols. The resultant study-level p values of a gene set were combined with Fisher's combined probability test. The fifty gene sets with the lowest meta-analysis p value were chosen as input for EQC score calculation. C2.all.4.0.symbols gene set database was chosen as input for consistency quality control (CQCp) value calculation. GSEA input expression matrices contained gene IDs that were present in all studies after EMA filtering. MetaDE package was used to determine DEGs in lesional samples compared to non-lesional ones. [15] DEG p value in individual studies was calculated by two sample T test with unequal variances. Fisher's combined probability test was chosen for meta-analysis statistical method. [16] Fold change of gene expression was given by the ratio between geometrical means of gene expression in lesional and non-lesional samples. [17] Genes with false discovery rate (FDR) less than 0.001 and with fold change higher than 1.5 or less than 21.5 were accepted as DEGs. Construction of protein-protein, protein-DNA and chemical-protein interaction networks STRING database 9.0 was used as resource for protein-protein interactions (PPI). [18] Both directed and undirected networks were created by selecting all interactions between DEG -coded proteins in downloaded raw data. Interaction confidence score cutoff was 900 (''highest confidence'' group) in case of undirected and 800 (containing a part of ''high confidence'' and all ''highest confidence'' interactions) in case of directed interactions. Only directed interactions with ''activation'' or ''ptmod'' actions were used. Chemical-protein interactions between potential drugs, intra-and extracellular compounds and DEG-coded proteins were collected from STITCH database 3.1. [19] The way of interaction confidence score calculation is the same in this database as in STRING thus interactions with the described confidence score cutoff values were selected for network construction. Protein-DNA interaction (PDI) network consisting of DEGs and DEG-coded transcription factors (TF) was created using cis-Regulatory Element Database (CisRED). [20] Regulatory element motifs with pv0:001 were collected from DEG promoter regions. Motifs were coupled with TFs or TF complexes using TRANS-FAC and JASPAR databases. [21,22] Motifs without respective TFs were excluded. Merged DEG-derived network containing PPI and PDI interactions and a network containing only DEG-coded TFs were also generated. Complete PPI, PDI, merged, TF-TF and chemical-protein interaction networks were created for controls using all available interactions in databases with the same statistical threshold as in DEG-derived network construction.

General network analysis, identification of central nodes and motif detection
General network analysis and node centrality value calculation were carried out with NetworkAnalyzer Cytoscape plugin. [23] Isolated nodes and node groups (without connection with the main PPI network) were deleted from graph in order to evade false results. Curve fitting on node degree and stress value distributions was done with MATLAB Curve Fitting Tool (MATLAB R2012b, The Mathworks Inc., Natick, MA). Curve of power law distribution was assessed with Trust-Region algorithm. Goodness of fitting was assessed by R-square and corrected R-square values which prove power law distribution of these node centralities ( Table 2). As power law distribution is asymmetric with a long tail, nodes with centralities above average cannot be assessed using arithmetic mean. A variable with a power-law distribution has a probability P k ð Þ of taking a value k following the function P k ð Þ*Ck {c , where C is constant. First moment (mean value) of a power-law distributed quantity equals: Second moment (variance) of a power-law distributed quantity equals: The sum of first and second moment (mean value and variance) was used as cutoff for centralities with distribution exponent cw3. Expression of variance becomes infinite, when cƒ3, thus only first moment (mean value) was used as cutoff for centralities with distribution exponent2vcv3. [24] Expression of mean value becomes infinite, if cƒ2. In this case weighted mean was used to assess cutoff with the following formula: As bidirectional connections are available in undirected PPI network, stress centrality is independent from edge directions thus both degree and stress had to be above cutoff for central protein selection. As directed networks contain unidirectional interactions, low stress values (i. e. low number of shortest paths cross through the node) can be caused by the dominance of incoming (in-degree) or outgoing (out-degree) interactions. Important nodes with high in-degree or out-degree can still have low stress centrality thus either out-degree or in-degree or stress had to be above cutoff in directed PPI network. As TFs have mainly outgoing interactions, out-degree was used for TF prioritization. Similarly to PPI networks degree and stress had to be above cutoff in undirected chemical -protein interaction network. Drugs with more targets in DEG-derived PPI-networks may have bigger disease modifying effect thus out-degree had to be above cutoff in directed chemical -protein interaction network for drug prioritization ( Table 2).
NetMODE software was used for network motif statistical analysis. Frequency of 3 or 4 node motifs in DEG-derived and complete control networks were compared with 1000 random graphs. Local constant switching mode was used for edge switching method during random network generation. NetMODE p value indicates the number of random networks in which a motif occurred more often than in the input network, divided by total number of random networks. pv0:05 was used as cutoff. [25] Respective sub-networks of enriched motifs were identified with NetMatch Cytoscape plugin. [26] jActiveModules and Cluster-ONE were used for network module and protein complex detection. ClusterONE analysis was carried out with minimum cluster size of 3 with unweighted edges and default advanced parameters. jActiveModules considers gene expression for module search. Input gene expression values have to be between 0 and 1 so normalized expression values got with EMA were scaled between these numbers. [27,28] Functional description of node groups was done with BinGO (''Biological function'' GO terms were selected, FDR,0.001 was used for term enrichment). [29] Results

Detection of DEGs with microarray meta-analysis
In order to get reliable data about gene expression in lesional psoriatic skin samples microarray meta-analysis was carried out. The study by Johnson-Huang et al. was already excluded after sample quality analysis with arrayQualityMetrics package, because at least two samples from one phenotype group are needed for MetaQC analysis and only one non-lesional sample remained after sample filtering. The overall quality of each study was assessed by MetaQC. [13] The software calculated six quality control (QC) measures then created principal component analysis (PCA) biplot and standardized mean rank summary (SMR) score to help in the identification of problematic studies. It was described by authors, that if a study is on the opposite side of arrows in the PCA biplot and has large SMR scores, it's strongly suggested to be excluded from meta-analysis. In contrary, if a study is on the same side of arrows in the PCA biplot and has small SMR scores, it should be included. All five studies were defined as usable based on quality values (Table 1, Figure 1). DEGs were identified by MetaDE. [15] 2307 upregulated and 3056 downregulated genes were found in lesional skin samples compared to non-lesional ones (Table S3). The relatively high number of DEGs can be the result of filtering out low quality samples, which could increase variance and using lower fold change cutoff values than in original studies. DEGs were used for network construction.

General Network analysis
Undirected and directed PPI networks with DEG -coded proteins, directed PDI networks with DEG -coded TFs and regulated DEGs and merged directed networks containing both PPIs and PDIs were created. A TF-TF network consisting of DEG-coded TFs was also generated. The Cytoscape plugin NetworkAnalyzer calculated main network properties for both DEG-derived and control complete networks (Table 3). DEGderived networks had higher diameter (i. e. the length of the longest shortest path in the network) and average shortest path length than control full networks. This may be caused by the inverse correlation of node degree and fold change. [30] Nodes with lower fold change has higher degree. Genes with fold change under cutoff are filtered out from DEG derived networks (between red lines on Figure 2). The remaining nodes has smaller average degree, therefore connectivity of the network is lower resulting in higher diameter and average shortest path length value.

Determination of hubs in DEG-derived networks
Most important nodes of DEG-derived networks were determined using degree and/or stress centralities ( Table 2, full list of nodes and centralities is in Table S4). Numerous already published psoriasis-associated protein-coding genes were found (Table 4). CCNA2, FYN and PIK3R1 proteins are present in top rated hubs in undirected PPI network and are yet unpublished in association with the disease. CCNA2 have role in mitosis regulation. [31] FYN is important in interferon gamma (IFN gamma) signaling, while PIK3R1 is important in insulin-stimulated glucose uptake. [32,33] FYN could be found in jActiveModules cluster with the 2 nd highest score while PIK3R1 were found in cluster with the 3 rd highest score ( Figure S1, S2). Taking account BinGO results these clusters  DEG derived and control networks has similar attributes, but average shortest path length and network diameter is lower in DEG derived networks, which can be explained by lower connectivity (Figure 2). Values for control networks are in brackets. doi:10.1371/journal.pone.0080751.t003 are responsible for signaling and for immune regulation as well (Table S5). A highly connected chemokine-chemokine receptor cluster was also found with ClusterONE analysis ( Figure S3). Central nodes in directed and undirected PPI networks showed overlap (Table 4). CTGF is in top ranked proteins and yet not associated with psoriasis. CTGF is responsible for fibrosis downstream of TGFb signaling. Downregulation of CTGF by psoriasis-associated cytokines INFc and TNFa is already published. [34] PDI network contained DEG-coded TFs and regulated DEGs as nodes and directed edges pointing from the TFs to the regulated genes. TFs were ranked using out-degree centrality. Androgen receptor (AR) and TFDP1 were the highest ranked nodes. AR is a TF, regulating genes that have immunological functions and role in carbohydrate metabolism. [35,36] TFDP1 controls cell cycle progression and is yet not associated with psoriasis. [37] BinGO analysis of TFDP1-regulated genes prove its central role in cell cycle activation (Table S5). MECOM and MEF2A are TFs above centrality cutoff and yet not associated with psoriasis. MECOM have role in cell proliferation and is associated with chronic myeloid leukemia. [38] MEF2A is responsible for the insulin dependent glucose transporter GLUT4 expression and is downregulated in insulin deficient diabetes mellitus. [39] Motif analysis in DEG-derived networks Motifs consisting of 3 or 4 nodes were analyzed in directed DEG-derived and control networks as well (Table 5, Figure 3). Analysis found motifs which were enriched in directed DEGderived but were absent in control networks or vice versa. Some were already generally described in biological systems like Identifying nodes making up motif no. 924 resulted in the high occurrence of central proteins found before. These proteins were associated with the immune system and carbohydrate metabolism.
Motif 332 is enriched in the TF network of lesional skin. This motif is based on the TFDP1-AR reciprocal regulation. Importance of these TFs is already mentioned. An interesting result of motif analysis is the enrichment of feedback loops containing 3 nodes in merged networks compared to separate ones and the enrichment of motif no. 6356 in DEGderived merged network compared to control. Motif no. 6356 consist of a positive feedback loop and all nodes of the loop are controlled by another separated node like IL1B or AR.

Controller sub-network construction
Both lesional and non-lesional skin areas can be found on patients at the same time. We wanted to highlight nodes which may be important in the ''all or none'' switch in lesional skin areas and sustain this phenotype for a long time. It has been argued that hubs in intracellular regulatory networks are enriched with either positive or negative regulatory links and cause much more positive feedback loops than negative ones. [40] It is also proven that positive feedback loops have fundamental role in maintaining autoimmune and autoinflammatory disease states. [41] Enrichment of motif no. 6356 consisting of a positive feedback loop with all nodes controlled by a separated one also suggests central role of positive feedback loops in lesional skin which may be activated by important central proteins like AR or IL1B. This is published that in biological systems interlinked slow and fast positive feedback loops allow systems to convert graded inputs (like several environmental and genetic factors in a psoriatic individual) into decisive all or none outputs (like lesional skin phenotype). [42] Transcriptional regulation needs time so we hypothesized that slow positive feedback loops may consist of at least one gene regulatory interaction. Fast loops may consist of only PPIs. Transcriptional changes of nodes in these loops may be able to sustain the ''switched on'' state.
In order to find most important slow and fast feedback loops containing 2, 3 or 4 nodes, a merged PPI and PDI network was constructed from proteins with centralities above cutoff value. All feedback loops were identified with NetMatch.   Table 5. doi:10.1371/journal.pone.0080751.g003 nodes as well. These may be key components in the metabolicinflammatory interplay in the pathomechanism of psoriasis. ''Slow'' positive feedback loops containing gene regulatory interactions and ''fast'' loops containing only PPIs were also found. All positive feedback loops had common nodes, thus a merged network was generated containing interlinked slow and fast positive feedback loops (Figure 4). Transcriptional changes of all nodes and influence of all edges supported the sustainment of lesional phenotype in this sub-network. Boolean analysis of the resultant controller network was also performed. Nodes with downregulated expression got value of 0 and nodes with upregulated expression got value of 1. Future state of nodes was set based on interactions ( Table 6). The output boolean values were the same as the input state values which prove the role of the controller network in the sustainment of present (lesional) phenotype. Chemical -protein interaction analysis further prove the importance of controller network.

Analysis of chemical-protein interaction networks
Undirected and directed chemical-protein interaction networks were constructed using STITCH database, which contains interactions between proteins and chemical compounds (internal non-protein substances, drugs and environmental substances). [19] Drugs or potential drugs were filtered out from chemicals and ranked by degree and stress centrality in case of undirected and out degree centrality in case of directed networks (Table S4).
Top ranked drugs were grouped into Anatomical Therapeutic Chemical (ATC) classes (Table 7). [43] KEGG DRUG was used for classification. [44] Results show a big overlap between undirected and directed network analysis. Best rated drugs consisted of retinoic acid, cholecalciferol, costicosteroids, methotrexate, sirolimus and tacrolimus, which can be already found in psoriasis guidelines and large clinical trials have proved their effectiveness. [45] Psoriasis studies are available for numerous potential drugs with high centralities. ''Blood glucose lowering drugs'' are promising drug candidates. The biguanide metformin is associated with reduced psoriasis risk in a population based case control study. [46] Many studies are available about ''Thiazolidinedione'' group. A recent meta-analysis showed significant decrease in Psoriasis Area and Severity Index (PASI) scores compared to placebo in case of pioglitazone and non-significant improvement in PASI 50/70 in case of rosiglitazone. [47] Troglitazone normalized histological features in psoriasis models and the lesional phenotype in a small clinical trial. [48] The ''HMG CoA reductase inhibitor'' drug simvastatin was effective in a pilot study, although atorvastatin in the same class showed only a non-significant improvement in a different study. [49,50] Salicylic acid has antifungal effects and it's used as adjuvant because of its keratolytic effect in the treatment of psoriasis. [51] The ''Antineoplastic agent'' methotrexate is a wellknown medication for psoriasis but several additional drugs in the same class were found in our analysis. Studies are available about 5-fluorouracil for the treatment of dystrophic psoriatic fingernails, but it showed only non-significant improvement. [52] Micellar paclitaxel significantly improved psoriasis in a prospective phase II study. [53] A study reported significant effectiveness of topical caffeine. [54] ''Calcium channel blocker'' nifedipine is found to be inductor of the disease in a case control study. [55] A study in 2005 reported significant PASI score reduction of 49.9% by topical theophylline ointment. [56] Mahonia aquafolium extract -consisting of berberine among others -is not classified into ATC classes, but three clinical trials already indicated improvement of psoriasis with this substance. [57] Multiple studies prove efficacy of the terpenoid triptolide in the treatment of psoriasis. [58] A recent study investigated effect of rifampicin on psoriasis and reported a 50.03% mean PASI reduction. [59] Study about the treatment of psoriasis with curcumin was carried out but reported only low response rate. [60] In an in vitro experiment the ''Lipid modifying agent'' clofibrate, but not bezafibrate reversed UVB-light-mediated expression of psoriasis -related inflammatory cytokines (interleukin-6, interleukin-8). [61] Fluvastatin and pravastatin have the potential to inhibit Th17 cell chemotaxis thus lowering immune cell infiltration of psoriatic skin. [62] Anti-proliferative effect of novel COX2 inhibitors on HaCaT keratinocytes was proven in an in vitro experiment and possible therapeutic use in psoriasis was supposed. However no such experiment was carried out with celecoxib which was the only COX2 inhibitor in best rated drugs. [63] N-acetyl-cysteine attenuated TNF alpha -induced cytokine production in primary human keratinocytes, which suggests its anti-psoriatic potential. [64] The ''Thiazolidinedione'' ciglitazone was never used as a medication, but inhibited keratinocyte proliferation in a dose dependent fashion. [48] Histone -deacetylase inhibitor trichostatin A blocked the conversion of regulatory T cells to IL17 expressing T cells suggesting its beneficial role in treating psoriasis. [65] Tse et al. suppose that antiproliferative effect of arsenic compounds could have positive effects on psoriatic skin. [66] The phosphodiesterase inhibitor rolipram has the ability to block enterotoxin B-mediated induction of skin homing receptor on T lymphocytes and may have the potential to inhibit lymphocytic infiltration of lesional skin. [67] The natural polyphenolic compound rottlerin is a potent inhibitor of NFkB and may have disease modulating effects. [68] Case reports are available about psoriasis induction by clonidine, ''agents acting on the renin-angiotensin system'' like captopril or losartan; the ''protein kinase inhibitor'' and ''antineoplastic agent'' imatinib; diclofenac, olanzapine, fluoxetine and chloroquine. Also case reports are available about the beneficial effects of ritonavir; ''antineoplastic agents'' like cytarabine, doxorubicin, and cysplatin; gefitinib, colchicine, lidocaine and nicotine. [69][70][71][72][73][74][75][76][77][78][79][80][81][82][83]    The 32 effective drugs of ''Studies available'' group in Table 7 were filtered out from STITCH data and target proteins were analyzed. All target proteins got an in-degree value reflecting the number of effective drugs acting on it. The group of proteins forming the controller sub-network was compared with the remaining target proteins. The controller sub-network protein group got significantly higher median value (10 vs. 1) using Mann-Whitney Rank Sum Test than the other one, which prove the importance of the controller sub-network in psoriatic lesions. ( Figure 5) (p,0.001; in-degree has power law distribution, thus Ttest could not be used) Higher median value could be caused by higher original degree centralities of controller network proteins in PPI networks, but only weak relation have been found between original degree centrality and the number of effective drugs acting on a protein, which cannot explain the big difference between the median of two groups (corrected R square value in regression analysis: 0.304) In summary, studies are available for 34 drugs found by our analysis, experimental evidence is available for 24 drugs, case reports suggest beneficial or disease-inductor effect of 21 drugs and 98 unpublished drug candidates for the treatment of psoriasis were also found ( Table 7-8).

Microarray Meta-analysis
Previous meta-analysis of psoriasis microarray studies was carried out by Tian et al. 1120 DEGs were found using 5 studies and 1832 DEGs using 3 studies. [84] We used the same 5 studies, but samples with inadequate quality were excluded from each study using arrayQualityMetrics package. The high number of DEGs (5363) in our study may be surprising, but it can be caused by the lower gene expression fold change cutoff (1.5 and 21.5 instead of 2 and 22). The pre -filtering process of samples can decrease variance and can also increase the number of DEGs. Further analysis of DEGs was carried out with Ingenuity Pathway Analysis (IPA) by Tian et al. IPA uses published references, carry out gene set enrichment analysis and TF detection. We used fundamentally different analysis. We generated PPI networks based on the largest PPI database (STRING) available which not only contain experimentally proven interactions but highly reliable interactions based on prediction algorithms or data mining. PDI network was also generated using not only literally proven interactions but interactions based on high fidelity prediction algorithms. Using lower DEG fold-change cutoff and detailed analysis based on node centrality statistics made it possible to identify proteins yet not associated with the disease but may have remarkable impact on pathogenesis. A chemical -protein interaction network based on STITCH database was also created and disease -modifying drug prediction was also possible with this method.

Keratinocyte hyperproliferation and Psoriasis
Keratinocyte hyperproliferation and inhibition of apoptosis are well-known phenomena in psoriasis. Several proteins have been associated with these mechanisms like BCL2, BAX, NFATC1, PPARd, EGF, mTOR, NF-kB etc. [85][86][87][88] Most of them were in central proteins detected by DEG-derived network analysis. Candidate DEG-coded proteins for hyperproliferation like CCNA2, TFDP1 and MECOM were also found. CCNA2 encodes Cyclin A2, that controls S phase and G2/M transition. Not only cell cycle progression is abnormal in lesional skin, but actin cytoskeleton organization as well. [89] A recent study reported that CCNA2 protein has role in cytoskeletal rearrangements and cell migration as well. [31] Cyclin A2 may take part in hyperproliferation and in aberrant actin cytoskeleton organization in psoriatic skin keratinocytes. TFDP1 encodes DP1 protein which is a dimerization partner of E2F transcription factor. The E2F/ DP1 heterodimers regulate cell cycle via DNA replication control and apoptosis. DP1 has E2F-independent function as well: DP1 can stabilize Wnt-on and Wnt-off states in Wnt/b-catenin signaling and determine differential cell fates. [37] TFDP1regulated genes belong to cell cycle progression as shown by BinGO analysis (Table S5). TFDP1 also has a reciprocal gene expression regulation with AR. This interaction was responsible for motif no. 332 enrichment in psoriasis PDI network compared to complete PDI network. This interaction may connect the hyperproliferation machinery to the merged controller subnetwork.

Immunological-metabolic interplay in psoriasis
Psoriasis is an immune-mediated disease. Some proteins which are published as important factors in pathogenesis were absent from DEGs in our microarray-meta analysis, such as TNF alpha, which is an important target in psoriasis therapy. This could be explained by the fact, that increased TNF alpha in psoriatic plaques can be caused mainly by post-transcriptional mechanisms. [90] Many proteins published in association with the immunopathogenesis of psoriasis were highly ranked hubs in PPI networks: IL1, IL8, TGFB1, SP1, STAT1, STAT3, NFKB1, IRF1 etc. [87,[91][92][93][94][95][96][97] A highly interconnected cluster mainly consisting of upregulated chemokines and chemokine receptors was also found by PPI analysis ( Figure S3). The downregulation of src kinase FYN seems to be a counteracting compensatory mechanism as this protein is important in IFN gamma action, in TNF alpha induced COX2 expression and in adipose tissue -mediated inflammation leading to insulin resistance. These processes are important in the pathomechanism of psoriasis. [32,98,99] These data suggest that the FYN inhibitor KBio2_002303 may have beneficial effects in the treatment of psoriasis. An important node in controller subnetwork is IL8. Although its role in psoriasis pathogenesis is published, no trial has been done with IL8 inhibitors. [100] This is true for CCL2 and IRF1 as well. Our study confirms their basic role in sustainment of lesional phenotype. Both can be found in highly ranked hubs and CCL2 is also essential in controller subnetwork by activating two positive feedback loops related to inflammation.
Psoriasis and metabolic syndrome comorbidity is a well-known phenomenon. There is a complicated interaction between the two diseases mediated by inflammatory cytokines among others. [101] Numerous DEG-coded proteins associated with both diseases could be found in central proteins like PPARG, INS-IGF2, LEP etc. [102][103][104] Others, like PIK3R1, AR and MEF2A may have role in the development of metabolic syndrome in psoriasis. PI3KR1 is important in the development of insulin resistance, it propagates inflammatory response in obese mice and may be an important link between the obesity-inflammation interplay in psoriasis. [33] AR has important effect on insulin signaling and thus insulin resistance. It is published that AR knockout mice exhibit insulin resistance. [35] To our knowledge AR has not yet been associated with psoriasis. However it was found in 1981, that lower serum testosterone level therefore decreased AR activation can be detected in psoriatic patients. [105] AR and PPARG connect inflammation-and metabolism-related hubs in controller network thus modulation of these proteins can be beneficial in psoriatic patients, which was also proven by our drug target analysis ( Figure 5). MEF2A is important for GLUT4 expression on insulin-responsive cells. Expression of MEF2A is downregulated in lesional skin samples which suggests a possible mechanism for insulin resistance in psoriasis. Many drugs, which are already widely used as treatment for psoriasis could be found in highly ranked nodes of chemicalprotein interaction networks such as methotrexate, retinoic acid, corticosteroids, sirolimus and tacrolimus. According to STITCH data all of them act through at least one of the hubs in controller sub-network. Top ranked ATC drug classes target members of controller sub-network as well. Blood glucose-lowering drugs act through PPARG and INS-IGF2 activation, which can be the basis of the positive effects of fibrate and HMG-CoA inhibitors in psoriasis as well. [47] Cardiac stimulants such as adrenergic agents also have high impact on lesional skin's PPI and PDI network, mainly by modulating hubs in controller sub-network. ''Sex hormones and modulators of the genital system'' ATC drug class act on AR. The ''antineoplastic drug'' methotrexate mainly acts through the accumulation of adenosine, but other antineoplastic agents may have their effect on keratinocyte hyperproliferation. [106] Studies or case reports already suggest efficacy of some antineoplastic drugs but several new possible agents were found in our analysis. [53,107,108] Mental stress is a known trigger for psoriasis and connection between the neuroendocrine system and skin immune system has been reported. [3,109] This is not surprising that numerous drugs acting on the CNS are enriched in highly ranked drugs. A lot of other drugs which are either classified in ATC classes or just drug candidates are found like kainic acid, cocaine, the HDAC inhibitor sodium butyrate, the PKC inhibitor bisindolylmaleimide I etc. (Table 7) In summary this is the first time PPI, PDI and chemical-protein interaction networks of psoriatic skin samples has been examined with detailed network analysis. Network-building DEGs were identified with fine-quality microarray meta-analysis of 187 nonlesional and 189 lesional samples. Several proteins were found which are yet not associated with psoriasis but may have high impact on the pathogenesis of the disease. Basic disease controller sub-network was also constructed consisting of central nodes coded by DEGs. Numerous anti-psoriatic drugs and drug candidates were also found acting mainly on these nodes. Figure S1 FYN protein in the jActiveModules cluster with 2 nd highest score. Nodes with blue-shaded color are downregulated and nodes with red-shaded color are upregulated. Color intensity is proportional with fold change. (PNG) Figure S2 PIK3R1 protein in the jActiveModules cluster with 3 rd highest score. Nodes with blue-shaded color are downregulated and nodes with red-shaded color are upregulated. Color intensity is proportional with fold change. (PNG) Figure S3 Chemokine-chemokine receptor cluster found by ClusterONE.

(PNG)
Table S1 Outlier samples in arrayQualityMetrics analysis and explanation of exclusion. Numbers from 1 to 5 indicate the number of method used by the software to assess quality. (XLSX)   Dataset S1 Results of arrayQualityMetrics analysis. Only html data can be found in directories, pdf files were deleted due to size restrictions. (ZIP)