Investigation of multiple sclerosis-related pathways through the integration of genomic and proteomic data

Study groups

All patients were diagnosed at Istanbul University-Cerrahpaşa, Cerrahpaşa School of Medicine, Department of Neurology, according to the McDonald criteria (Polman et al., 2011). Eleven unrelated MS patients were randomly selected among the study group of our previous work comprising of 179 MS patients. The patient group had a heterogeneous disease presentation at the time of their CSF sample collection. Sixty independent age- and gender-matched intrinsic healthy controls were included in the analyses. All individuals in the study group are of Turkish origin. All procedures of the study were in accordance with the Helsinki Declaration of 1964 and its later amendments. The Ethics Committee of Istanbul Technical University approved the study (İTÜ-SM.İNAREK-MBG-1), and each individual of the study group gave a written informed consent form prior to the sample collection.

Genome-wide associations

DNA isolation from blood samples was performed (Roche DNA Isolation Kit for Mammalian Blood), and genotyping for 300.000 SNP markers was applied for each individual on the Illumina HumanCytoSNP-12 array. Illumina GenomeStudio software was used for quality control. All individuals showed sample call rates of more than 98% (98.65–99.63%) and were therefore included in the study. Golden Helix SNP & Variation Suite software was used for identity-by-descent detection, suggesting no relatedness between the individuals. SNP filtering was performed: SNPs on the Y chromosome, with call rates lower than 95%, minor allele frequency lower than 0.01, and in strong linkage disequilibrium (r² > 0.5) were excluded from the study. A total of 129.547 SNPs was included in the analysis. Frequency differences of SNPs between cases and controls were assessed, and genotypic association P-values were obtained.

Pathway enrichment analysis

Genotypic associations and proteomic data were subjected to pathway enrichment analyses using the functional enrichment tool PANOGA, which reveals functionally relevant pathways by identifying genes within the pathways, incorporating protein-protein interaction (PPI) information, and extracting significant pathways. SNPs with genotypic association P-values lower than 0.05 and differentially-expressed proteins with P-values lower than 0.05 from our previous study were used for PANOGA procedure as follows: For SNP associations, because a SNP might affect more than one gene, each SNP was initially assigned to the gene on which the SNP has the most important functional effect, using the SPOT webserver (Saccone et al., 2010). Functional information of SNPs was obtained utilizing SPOT, F-SNP (Lee & Shatkay, 2007), SNPnexus (Dayem Ullah, Lemoine & Chelala, 2012), and SNPinfo (Xu & Taylor, 2009). Genes and proteins were then mapped onto a PPI network, for which Goh et al.’s human PPI network was used (Goh et al., 2007). Next, the jActive Modules algorithm (Ideker et al., 2002 was applied to identify active subnetworks containing a large number of the disease-affected genes and proteins in the PPI network. Each genotypic association and differential expression P-value was taken into account, and active subnetworks that overlap at most 50% with each other were extracted. In order to evaluate the biological importance of the subnetworks, the number of genes and proteins found in a specific KEGG pathway was compared to the total gene and protein number involved in the corresponding pathway. For this functional enrichment procedure, a two-sided hypergeometric test was used, and the Bonferroni method was applied for multiple testing corrections of the P-values. Significantly associated KEGG pathways to our patient group consisted of significantly enriched (P < 0.05) pathways for at least one of the active subnetworks. For the pathways that were enriched in multiple subnetworks, the one with the lowest P-value was reported. PANOGA was run 10 times for both genotypic associations and differentially-expressed proteins, and the lowest P-values over the 10 iterations were reported. The resulting enriched pathways for both datasets were then merged. If a given pathway was identified in both analyses, the corresponding P-values were merged using Fisher’s combined probability test. If the pathway was identified in only one of the analyses, the corresponding P-value was used as the final P-value.

Clustering of enriched pathways

Using a method previously described by Chen et al. (2014) we clustered the enriched pathways to identify similar groups and establish representative pathways. The clustering approach can be summarized as follows: initially overlap index matrix (OI) that consists of the overlap indices between all of the pairs of enriched pathways was calculated. The overlap index (OI_i,j) of a pair of pathways P_i and P_j (denoting the ith and jth pathways, respectively) was defined as: ${\displaystyle O{I}_{i,j}=\frac{|G_i\cap G_j|}{min\left(\left|G_i\right|,\left|G_j\right|\right)}}$ where G_i is the set of all genes in the ith pathway. In the overlap matrix, each row (o_i) represents the gene overlap profile of a pathway (i.e.,  o_i is the gene overlap profile of the ith pathway). To identify similarity between each pair of pathways Pearson correlation coefficients (R) between each pair of overlap profiles (e.g.,  o_i and o_j) were calculated. These correlation coefficients were then converted to pairwise distances (PD) as: ${\displaystyle PD=1-R.}$ Using PD as the distance metric, hierarchical clustering (average-linkage) of the enriched pathways was performed. Examining the hierarchical clustering dendrogram, the dendrogram was partitioned into clusters at the manually selected PD cut-off value of 0.55. The pathway with the lowest P-value in each cluster was selected as the representative pathway for that cluster.

SNP associations

All 71 individuals and a total of 129.547 SNPs passed the quality check and therefore were included in the final analyses. The Manhattan plot summarizing the genome-wide associations is given in Fig. 1. The most significantly associated SNP was rs7873, located on 3′ UTR of the IGF2 gene, whose significance level was the only value close to a classical stringent P-value cut-off for a GWAS (P = 4.39 × 10⁻⁰⁷). When considering SNPs with a significance level of lower than a looser P-value cut off (P < 10⁻⁰⁵) given the small sample size, 4 more SNPs have high associations with MS in the studied patient group, one of which is located on a gene: rs17187282 (P = 1.17 × 10⁻⁰⁶), rs11688088 (P = 2.8 × 10⁻⁰⁶), rs654188 (P = 5.24 × 10⁻⁰⁶), and rs7092208 (P = 9.66 × 10⁻⁰⁶, intronic variant on the MGMT gene).

For enrichment analysis of the SNP associations, the P-value threshold was set lower than 0.05 (Fig. 1, blue line) in order to prevent the elimination of potentially MS-associated SNPs with falsely high P-values due to the low sample size and limitations of the genome-wide association methodology. Later, elimination of the false-positive SNP associations resulted from preferring not to set a conventional stringent P-value was aimed to be achieved during the identification of MS-related pathways using SNP subnetworks, therefore ruling functionally irrelevant SNPs out. A total of 6.594 SNPs had P-values under the threshold of 0.05 (Table S1), which were included in the subsequent analyses.

MS-related pathways

The combination of enriched pathways for the resulting SNPs and differentially expressed proteins resulted in 151 enriched pathways in total (Table S2). In order to identify pathways with similar content and function, this combined list of enriched pathways was clustered (Fig. 2). The dendrogram was then manually partitioned into biologically relevant clusters, and representative pathways were established. In total, 33 clusters, therefore 33 representative pathways were obtained (Table S3). Nine of the representative pathways emerged from both genomic and proteomic analyses, among which the complement and coagulation cascade (hsa04610) was the most significantly associated pathway in the studied group (P = 6.96 × 10⁻³⁰). Figure 3 shows the complement and coagulation cascades with genes identified through the genomic dataset and differentially expressed proteins identified through the proteomic dataset as a representative pathway for the relationships between the genomic and proteomic findings.

Regulation of actin cytoskeleton (hsa04810) and focal adhesions (hsa04510) were the two systems with the second and third lowest P-values even though there were no significant protein level changes related to these pathways (P = 9.64 × 10⁻²⁷ and P = 3.29 × 10⁻²⁶, respectively). Other 8 shared pathways with high significance levels are as follows: Adherens junctions (hsa04520, P = 5.38 × 10⁻¹⁷), colorectal cancer (hsa05210, P = 1.70 × 10⁻¹⁵), pathogenic Escherichia coli infection (hsa05130, P = 7.02 × 10⁻⁹), prion diseases (hsa05020, P = 2.48 × 10⁻⁵), endometrial cancer (hsa05213, P = 6.00 × 10⁻⁹), non-small cell lung cancer (hsa05223, P = 1.40 × 10⁻⁸), bladder cancer (hsa05219, P = 4.83 × 10⁻⁶), and non-homologous end-joining (hsa03450, P = 4.27 × 10⁻⁵). Table 1 shows the detected SNP associations and differentially expressed proteins in the above-mentioned pathways.

The most significantly MS-associated pathways obtained from the SNP associations and differentially expressed proteins did not overlap when the analyses were performed separately (Table S2). The merged data analysis resulted in a combined list of pathways from both datasets, revealing a different set of pathways with the lowest P-values and emphasizing the most relevant pathways that could not have been prioritized by a straight-forward, non-combinatorial omics approach. The top five pathways that emerged from the proteomic-only, genomic-only, and merged dataset analyses are given in Table 2.