Integrated network analysis of transcriptomic and proteomic data in psoriasis

Skin biopsies

Acquisition of the human tissue was approved by the Vavilov Institute of General Genetics of Russian Academy of Sciences review board and the study was conducted after patient's consent and according to the Declaration of Helsinki Principles.

A total of 6 paired nonlesional and lesional (all were plaque-type) skin biopsies from 3 psoriatic patients were profiled using 2D electrophoresis. All the donors who gave biopsy tissue (both healthy controls and individuals with psoriasis) provided a written informed consent for the tissue to be taken and used in this study. Clinical data for all patients are listed in Additional file 1.

Full-thickness punch biopsies were taken from uninvolved skin (at least 2 cm distant from any psoriatic lesion; 6 mm diameter) and from the involved margin of a psoriatic plaque (6 mm diameter) from every patient.

Sample preparation

Skin biopsies from lesional (n = 3) and non-lesional (n = 3) skin (~50 mg each) were mechanically homogenized on ice with mortar and pestle in solubilization buffer comprised of 7 M urea (Fisher Scientific), 2 M thiourea (MERK), 4% (w/v) CHAPS (Fluka), 0.5% (w/v) Triton X-100 (ICN), 0.5%(w/v) ampholines 3/10 (Bio-Rad), 20 mM Tris base (Fisher Scientific), 1 mM MgCl₂ (Fluka), and 0.2 mM PMSF (Acros). The amount of solubilization buffer taken per sample was seven times tissue wet weight. Samples were then carefully sonicated at 4°C for 15 sec (3 bursts for 5 sec) with Bandelin sonicator at 50% power output. Following homogenization and sonication DTT (Acros) was added to reach the final concentration of 65 mM, and samples left at 4°C for 30 min. Then Na₂EDTA was added (2 mM final) and mixture was incubated for additional 2 h and then centrifuged at 12,000 × g for 10 min to remove insoluble debris. Relative protein concentration was determined with GS-800 Calibrated Densitometer (Bio-Rad) after 1D-SDS-PAGE of 3 μl sample solution loaded on 8 × 9 cm mini-gel stained with Coomassie R-250 after electrophoresis.

Two-dimensional electrophoresis

For the first dimension CA-IEF method was used [20]. Glass tubes 20 × 1.5 mm (Bio-Rad) was filled with gel mixture containing 8.3 M urea (Fisher Scientific), 4% (w/v) acrylamide monomers (Acros), 2% CHAPS (Fluka), 1.6% (w/v) 5/8 and 0.4% (w/v) 3/10 ampholines (Bio-Rad), loaded with 120 μg protein and run 15 min at 200 V, 30 min at 300 V, 16 hrs at 400 V, and 1 h at 800 V. Isoelectrofocusing was performed in Protean II xi cell (Bio-Rad). Before second dimension gels were extruded from tubes and equilibrated in 60 mM Tris-buffer pH 6.8, containing 2% (w/v) SDS, 10% (w/v) glycerol, and 1% (w/v) DTT for 30 min. For protein separation in the second dimension, 9-16% linear gradient slab gels prepared according to the standard protocol [20] were used.

Gel image analysis

Protein spots on 2-DE gels were visualized by silver staining [21] and scanned with resolution 300 dpi. Images were analyzed using Melanie III software (GeneBio, Switzerland). The conventional analysis involved (i) protein spot relative volume (%Vol) determination, which was expressed as the sum of pixel intensities in the certain spot divided to the sum of pixel intensities in all spots on the gel; (ii) gel alignment; and (iii) spot matching. Further, sets of %Vol values for every spot were processed by Student test, thereby testing whether there was a significant variation of the certain protein level between two specified groups.

Protein identification by MALDI-TOF mass-spectrometry

Protein spots were cut out (~3 mm3) from 2-DE gels, destained, and in-gel digested with trypsin. Mass-spectrometry of trypsin digested proteins (spots No. 1,2,3,4) was performed using a Microflex MALDI-TOF mass-spectrometer (Bruker, Germany). Peptide samples (0.2-1 μl) were mixed with an equal volume of 2,5-dihydroxybenzoic acid solution (20 mg/ml; Sigma, USA) in 20% acetonitrile and 0.1% trifluoroacetic acid, and the resulting droplets were dried in air. Mass-spectra were obtained for mass range from 800 to 4000 daltons in reflection mode and calibrated using internal standards (trypsin autolysis peaks, MH+ 1046.54, 2212.10 daltons). Peptide peak lists were formed by the SNAP algorithm (XMass software, Bruker). Proteins were identified using the Mascot database search engine. The search parameters were as follows: mass tolerance 100 ppm, NCBI protein sequence database, Homo sapiens taxon, one missed cleavage, variable modifications by propionamide for cysteines and oxidation for methionines.

Low molecular weight proteins (9-20 kDa, spots No. 5,6,7,8,9,10) were identified using nanoLC-MS/MS mass-spectrometry for higher convenience, regarding lack of cleavage sites in such proteins. Analysis of trypsin digest was performed on electrospray ion trap (XCT Ultra Ion Trap Chip Cube 6330 series, Agilent) equipped with chip-cube head. One μl of each sample was subjected (flow rate 3 μl/min) onto reverse-phase in-chip column (40 nl capacity) for 10 minutes under isocratic buffer (5% acetonitrile in 0.1% TFA). Following sample application peptides were separated by linear gradient (0.3 μl/min) of solution A (0.1% TFA/water) and solution B (80% ACN/0.1% TFA/water) for 60 minutes. Spectrum scanning was repeated three times for each sample of protein gel spots and tissue hydrolysate. Mass spectra of eluted peptide were simultaneously obtained under positive polarity for 425-1325 m/z range both in MS and MS/MS mode, 2.1 kV applied accumulation of 85000 ions for 50 milliseconds, averages on 2 spectra. Mass spectra were processed with Spectrum Mill MS Proteomics workbench software (Agilent). Proteins were identified using SwissProt Human Database in the following parameters: score ≥ 7 for peptide and ≥ 20 for protein, minimum S/N ratio 20, maximum peptide ion charge +4, precursor mass tolerance ± 2.5 Da, product mass tolerance ± 0.7 Da, Proteins identification was accomplished with detection of minimum 3 peaks of the same peptide ion with maximum of 2 missed cleavages.

Microarray data analysis

We used recently published data set [22] from GEO data base (http://www.ncbi.nlm.nih.gov/geo/; accession number GSE14095). We compared 28 pairs of samples (in each pair there was a sample of lesional skin and a sample of healthy skin from the same patient). Values for each sample were normalized by sample median value in order to unify distributions of expression signals. For assessment of differential expression we used paired Welch t-test with FDR correction [23]. Probe set was considered as differentially expressed if its average fold change exceeded 2.5 and FDR corrected p-value was less than 0.01.

Overconnection analysis

All network-based analyses were conducted with MetaCore software suite http://www.genego.com. This software employs a dense and manually curated database of interactions between biological objects and variety of tools for functional analysis of high-throughput data.

We defined a gene as overconnected with the gene set of interest if the corresponding node had more direct interactions with the nodes of interest than it would be expected by chance. Significance of overconnection was estimated using hypergeometric distribution with parameters r - number of interactions between examined node and the list of interest; R - degree of examined node, n -sum of interactions involving genes of interest and N - total number of interactions in the database:

Hidden nodes analysis

In addition to direct interacting objects, we also used objects that may not interact directly with objects of interest but are important upstream regulators of those[24]. The approach is generally the same as described above, but the shortest paths instead of direct links are taken into account. As we were interested in transcriptional regulation, we defined a transcriptional activation shortest path as the preferred shortest path from any object in the MetaCore database to the transcription factor target object from the data set. We added an additional condition to include the uneven number of inhibiting interactions in the path (that's required for the path to have activating effect). If the number of such paths containing examined gene and leading to one of objects of interest were higher than expected by chance, this gene was considered as significant hidden regulator. The significance of a node's importance was estimated using hypergeometric distribution with parameters r - number of shortest paths between containing currently examined gene; R - total number of shortest paths leading to a gene of interest through transcriptional factor, n -total number of transcription activation shortest paths containing examined gene and N - total number of transcription activation shortest paths in the database.

Rank aggregation

Both topology significance approaches produced lists of genes significantly linked to a gene or protein set of interest, ranked by corresponding p-values. To combine results of these two approaches, we used a weighted rank aggregation method described in [25]. Weighted Spearman distance was used as distance measure and the genetic algorithm was employed to select the optimal aggregated list of size 20. This part of work was accomplished in R 2.8.1 http://www.r-project.org.

Network analysis

In addition to topology analysis, we examined overexpressed genes and proteins using various algorithms for selecting connected biologically meaningful subnetworks enriched with objects of interest. Significance of enrichment is estimated using hypergeometric distribution.

We first used an algorithm intended to find regulatory pathways that are presumably activated under pathological conditions. It defines a set of transcription factors that are directly regulating genes of interest and a set of receptors whose ligands are in the list of interest and then constructs series of networks; one for each receptor. Each network contains all shortest paths from a receptor to the selected transcriptional factors and their targets. This approach allows us to reveal the most important areas of regulatory machinery affected under the investigated pathological condition. Networks are sorted by enrichment p-value.

The second applied algorithm used was aimed to define the most influential transcription factors. It considers a transcriptional factor from the data base and gradually expands the subnetwork around it until it reaches a predefined threshold size (we used networks of 50 nodes). Networks are sorted by enrichment p-value.

Differentially abundant proteins

Protein abundance was determined by densitometric quantification of the protein spots on 2D-electophoresis gel (Figure 1; see also Additional file 4) followed by MALDI-TOF mass spectrometry. Total of 10 proteins were over-abundant at least 2-fold in lesional skin compared with uninvolved skin: Keratin 14, Keratin 16, Keratin 17, Squamous cell carcinoma antigen, Squamous cell carcinoma antigen-2, Enolase 1, Superoxide dismutase [Mn], Galectin-7, S100 calcium-binding protein A9 and S100 calcium-binding protein A7 (Table 1). Several of these proteins were previously reported to be over-abundant in psoriatic plaques [26–29].

The proteins belonged to a diverse set of pathways and processes. Thus, keratin 17, keratin 14, and keratin 16 are a member of the keratin gene family. The keratins are intermediate filament proteins responsible for the structural integrity of epithelial cells. SERPINB4 and SERPINB3 are serine protease inhibitor to modulate the host immune response against tumor cells. Enolase 1, more commonly known as alpha-enolase, is a glycolytic enzyme expressed in most tissues, one of the isozymes of enolase. Superoxide dismutase 2 protein (SOD2) transforms toxic superoxide, a byproduct of the mitochondrial electron transport chain, into hydrogen peroxide and diatomic oxygen. Galectins are a family of beta-galactoside-binding proteins implicated in modulating cell-cell and cell-matrix interactions. Differential and in situ hybridizations indicate that this lectin is specifically expressed in keratinocytes. The cellular localization and its striking down-regulation in cultured keratinocytes imply a role in cell-cell and/or cell-matrix interactions necessary for normal growth control. S100A9 and S100A7 proteins are localized in the cytoplasm and/or nucleus of a wide range of cells, and involved in the regulation of a number of cellular processes such as cell cycle progression and differentiation. S100A7 is markedly over-expressed in the skin lesions of psoriatic patients.

We attempted to connect the proteins into a network using a collection of over 300,000 manually curated protein interactions and several variants of "shortest path" algorithms applied in MetaCore suite [30] (Figure 2, see Methods for details). The genes encoding overabundant proteins were found to be regulated by several common transcription factors (TFs) including members of the NF-kB and AP-1 complexes, STAT1, STAT3, c-Myc and SP1. Moreover, the upstream pathways activating these TFs were initiated by the overabundant S100A9 through its receptor RAGE [31] and signal transduction kinases (JAK2, ERK, p38 MAPK). This network also included a positive feedback loop as S100A9 expression was determined to be controlled by NF-kB[32]. The topology of this proteomics-derived network was confirmed by several transcriptomics studies [33–38] which showed overexpression of these TFs in psoriasis lesions. Transiently expressed TFs normally have low protein level and, therefore, usually fail to be detected by proteomics methods.

RAGE receptor is clearly the key regulator on this network and plays the major role in orchestrating observed changes of protein abundance. This protein is abundant in both keratinocytes and leukocytes, though normally its expression is low [39]. RAGE participates in a range of processes in these cell types, including inflammation. It is being investigated as a drug target for treatment of various inflammatory disorders [40]. Thus, we may propose that RAGE can also play significant role in psoriasis.

Differentially expressed genes

We used Affymetrix gene expression data set from the recent study [22] involving 33 psoriasis patients. Originally, more than 1300 probe sets were found to be upregulated in lesions as compared with unlesional skin of the same people. We identified 451 genes overexpressed in lesional skin under more stringent statistical criteria (28 samples of lesional skin were matched with their non-lesional counterparts from the same patients in order to exclude individual expression variations, genes with fold change >2.5 and FDR-adjusted p-value < 0.01 were considered as upregulated). The list of overexpressed genes can be found in Additional file 2. The genes encoding 7 out of 10 proteomic markers were overexpressed, well consistent with proteomics data. Expression of Enolase 1, Keratin 14 and Galectin 7 was not altered.

Common transcription regulation for overexpressed genes and differentially abundant proteins

Despite good consistency between the proteomics and expression datasets, the two orders of magnitude difference in list size make direct correlation analysis difficult. Therefore, we applied interactome methods for the analysis of common upstream regulation of the two datasets at the level of transcription factors. First, we defined the sets of the most influential transcription factors using two recently developed methods of interactome analysis [41] and the "hidden nodes" algorithm [42]. The former method ranks TFs based on their one-step overconnectivity with the dataset of interest compared to randomly expected number of interactions. The latter approach takes into account direct and more distant regulation, calculating the p-values for local subnetworks by an aggregation algorithm [42]. We calculated and ranked the top 20 TFs for each data type and added several TFs identified by network analysis approaches (Table 2). The TFs common for both data types were taken as set of 'important pathological signal transducers' (Figure 3). Noticeably, they closely resemble the set of TFs regulating the protein network on Figure 2.

Transcriptomics				Proteomics
Entrez ID	Gene	Over-connection P-Value	Hidden nodes P-Value	Entrez ID	Gene	Over-connection P-Value	Hidden nodes P-Value
6772	STAT1	4.081E-28	1.027E-22	2736	GLI2	4.318E-03	2.080E-04
3659	IRF1	1.344E-22	3.203E-19	6772	STAT1	2.234E-04	1.293E-04
4790	NFKB1	3.677E-13	2.523E-25	3091	HIF1A	2.345E-03	2.402E-04
5970	RELA	3.726E-14	4.120E-24	6778	STAT6	9.420E-04	2.080E-04
5966	REL	4.844E-17	3.580E-16	7392	USF2	1.507E-04	-
6688	SPI1	7.019E-13	1.245E-16	6774	STAT3	2.542E-05	-
3394	IRF8	7.594E-12	1.903E-16	6776	STAT5A	-	2.820E-05
1051	CEBPB	1.590E-12	1.271E-15	668	FOXL2	1.100E-02	1.388E-03
2908	NR3C1	2.307E-05	1.504E-17	2078	ERG	2.386E-04	1.388E-03
10379	IRF9	8.193E-05	4.037E-17	6777	STAT5B	-	2.820E-05
3660	IRF2	3.794E-12	3.914E-09	7020	TFAP2A	4.806E-04	-
6434	SFRS10	-	4.420E-17	5371	PML	1.177E-03	-
6775	STAT4	1.637E-03	5.842E-16	2735	GLI1	1.053E-02	-
6667	SP1	5.508E-12	2.013E-07	6667	SP1	3.269E-05	-
6777	STAT5B	2.725E-04	2.063E-15	3725	JUN	1.198E-03	-
3725	JUN	2.182E-11	7.521E-11	8462	KLF11	6.710E-03	2.775E-03
6774	STAT3	9.643E-10	2.114E-11	2099	ESR1	1.650E-02	-
7157	TP53	9.352E-10	7.171E-08	10538	BATF	1.765E-02	-
30009	TBX21	5.119E-03	1.594E-14	30009	TBX21	-	1.911E-04
1874	E2F4	1.537E-09	8.737E-07	7421	VDR	-	1.693E-04
Additional TFs identified by subnetwork-based analysis
4609	MYC	2.680E-69	-	4609	MYC	1.120E-07	-
				5966	REL	1.120E-07	-
				5970	RELA	1.120E-07	-
				4790	NFKB1	1.120E-07	-

Some genes were considered significant only by one of two topological approaches (this is evident for proteomics data, where low number of proteins limits capabilities of topological analysis). Missing p-value means that correspondent gene has not passed 0.05 significance threshold and has been listed among top factors only due to low p-value for other topological approach. Only one p-value was determined in the case of network analysis (enrichment p-value for the network built around seed transcription factor). Transcriptional factors common for proteomics and transcriptomics level are in bold text.

Identification of influential receptors

In the next step, we applied "hidden nodes" algorithm to identify the most influential receptors that could trigger maximal possible transcriptional response. In total, we found 226 membrane receptors significantly involved into regulation of 462 differentially expressed genes ('hidden nodes' p-value < 0.05). The complete list of receptors can be found in Additional file 3. Assuming that topological significance alone does not necessarily prove that all receptors are involved in real signaling or are even expressed in the sample; we filtered this list by expression performance. The receptors used were those whose encoding genes or corresponding ligands were overexpressed greater than 2.5 fold. We assumed that the pathways initiated by over-expressed receptors and ligands are more likely to be activated in psoriasis. Here we assumed that expression alterations and protein abundance are at least collinear. An additional criterion was that the candidate receptors had to participate in the same signaling pathways with at least one of the common TFs. No receptor was rejected based on this criterion.

In total, 44 receptors passed the transcription cut-off. Of these 24 receptor genes were overexpressed; 23 had overexpressed ligands and 3 cases had overexpression of both ligands and receptors (IL2RB, IL8RA and CCR5; see Figures 4 and 5 and Additional file 4). Interestingly, for several receptors, more than one ligand was over-expressed (Figure 4). Several receptors are composed of several subunits, only one of which was upregulated (for example, IL-2 receptor has only gamma subunit gene significantly upregulated).