Gene prioritization in Type 2 Diabetes using domain interactions and network analysis

Functional relevance of network motif analysis

Network motifs that are the indicators of specific functional modules in cellular networks [12, 13] were used to predict functions for hypothetical HWEc genes. For instance, network motif analysis showed that Q6MZN5 protein (Wv = 0.93, C_i = 0.0) with domain PF05277 forms a four node motif with partners PF01565 (electron transport), PF02913 (electron transport), PF02127 (proteolysis) suggesting its role in energy metabolism. The predicted functions were also cross-validated using the Protfun server by matching with the functional classes assigned by Protfun. Protfun server also predicted the involvement of Q6MZN5 protein in energy metabolism. Energy generation in mitochondria occurs primarily through oxidative phophorylation and as genes involved in oxidative phosphorylation are known to be coordinately downregulated in T2D, Q6MZN5 might attain significance [14].

Performance of the methods

The benchmarking dataset used to evaluate the performance of our methods, comprised of 19 genes known to be associated with T2D from GWA studies [1, 15] and 353 non-T2D genes lying within 10 Mb regions of 12 chromosomal positions those are never shown to be linked with any disease, as control data set. We observed sensitivity of Wv method to be 0.73 and specificity to be 0.96. For HRC method, sensitivity and specificity were 0.32 and 0.98 respectively. Overall data enrichment of 9 fold and 6.8 fold was observed for the two methods respectively. Using the benchmarking dataset, we compared both Wv and HRC predictions with other known methods of candidate gene prioritization (Table 2). The receiver operating characteristic (ROC) analysis of the two methods, Wv and HRC confirms their better performance (Fig. 2), with only PROSPECTR [16] having a higher sensitivity (0.90) than Wv and HRC. Where, G2D method gave accuracy and sensitivity of 0.73 and 0.53 respectively with the benchmarking dataset.

	Wv	HRC	PROSPECTR	SUSPECTS	G2D	DGP
Accuracy (%)	94.9	94.6	41.4	61.1	72.8	49.8
Sensitivity (%)	73.7	31.6	89.5	52.6	52.6	58.8
Specificity (%)	96	98	38.8	61.5	73.9	49.3

Numbers in bold letters indicate significant values

Further, other gene prioritization methods including PROSPECTR [16], SUSPECTS [17], Disease Gene Prediction (DGP) [18] and G2D [19] were used to prioritize our dataset of 4052 genes. PROSPECTR and SUSPECTS gave a disease gene estimate of 1 in 1.7 (2275 in 3972 and 2254 in 3612 respectively), DGP did 1 in 2.1 (1720 in 3632) and for G2D it was 1 in 7 (524 in 4052). Whereas, Wv method gave disease gene estimate of 1 in 4 (995 in 4052), Wv + HRC provided an estimate of 1 in 9.31 (435 in 4052). Thus, it was striking that PROSPECTR, SUSPECTS and DGP predicted a higher number of genes to be disease candidates compared to our method. Among the other methods, G2D performed the best. While comparing the common candidates predicted by Wv + HRC and other methods we found that 84% (364) were prioritized to be disease candidate by any one of the methods while 16% (71) of the candidates were the unique candidates prioritized by our method (Additional file 2).

Prioritized Candidates by Wv + HRC method

Available literature suggests that candidates with high Wv and low clustering coefficient values, as prioritized by our method, show considerable relevance to T2D. For instance, the highest ranked candidate GPC1 (Wv-1.0, C_i-0.0) is shown to influence FGF2 signaling pathway [20]. The angiogenic growth factor, FGF2 levels has been found to be associated with cardiovascular events in T2D [21]. LPIN2 and LPIN3 (Wv-0.98, C_i-0.0) are the members of lipin family known to serve as an enzyme for triacylglycerol synthesis as well as transcriptional coactivator in the regulation of lipid metabolism genes [22]. CDK5R2 (Wv-0.98, C_i-0.0) forms the complex with CDK5, a serine/threonine protein kinase that plays crucial role in physiological functions such as glucose stimulated insulin secretion in pancreatic cells [23]. CX3CL1 (Wv-0.97, C_i-0.0), a CX3C chemokine has been shown to have specific role in initiation and progression of atherosclerotic vascular disease [24]. VEGF has already been shown to be associated with T2D complications [25]. PDGFB elevated levels and its induction by PKC activation has been shown to be involved in pathogenesis of diabetic retinopathy [26]. Mutation in AKT2 (R274H) enzyme has been shown to result in autosomal dominant inheritance of severe insulin resistance and T2D [27]. CAMK2 kinase found to be expressed in the pancreatic β-cell has been shown to affect insulin secretion [28]. Fetuin-A has been shown to be an important modulator of insulin resistance [29]. CETP plays an important role in the regulation of HDL metabolism and is shown to be associated with dyslipidemia in GWA studies [30]. DMPK, a serine/threonine protein kinase has been shown to be a positive modulator of insulin action [31]. NDUFA5, NDUFS7, SCO1 and FMO4 all are involved in oxidative phosphorylation and have been shown to be important in diabetes [14]. MOGAT2 was found to be associated with triacylglycerol synthesis and has a role in diet-induced obesity [32]. Interestingly, we found Bactericidal/permeability-increasing protein (BPI) as a highly ranked candidate and it has been suggested over production of BPI could be linked to insulin sensitivity and glucose tolerance [33]. A variation in Cathepsin S (CTSS), a cysteine protease has been shown to be associated human metabolic risk factors for cardiovascular diseases [34]. The other ranked candidate PDE3A in rat has been found to have antilipolytic action of insulin in adipocytes [35].

Selection of dataset

We searched the available literature for genome wide linkage scans in T2D and selected 108 microsatellite markers located on 64 different chromosomal regions with LOD scores ≥ 2.0 to select regions of profound statistical relevance with T2D (Additional file 3). All genes coding for known and unknown proteins within 5 Mb upstream and downstream of the selected markers (10 Mb region) were extracted from Ensembl database (v38) [42]. We ensured a non-redundant dataset comprising of 5441 candidate genes by screening for them in International Protein Index (IPI) [43].

Domain assignment and prediction of partner proteins

The domain definitions used in this study were obtained from Pfam database [44]. The domains were assigned by scanning libraries of Hidden Markov models (HMMs) against the protein sequences resulting in a dataset of 4052 candidates which alone were considered for further analysis. Interacting partner domains for each parent domain (domains contained in the positional candidates) were determined using InterDom database with a confidence score of interactions ≥10 [41]. We utilized 26,058 domain interactions for Wv calculations from InterDom database (Version 2.0) after applying the cut-off of ≥10 for confidence score. These interactions were found to be evident either in DIP, BIND or PDB databases [45–47]. After obtaining the domain partners, all the proteins harboring them were extracted using SwissPfam [44].

Calculation of Weight value using disease specific properties

Binary score of 1 or 0 was assigned to the partner domains on the basis of their presence or absence respectively in the following four classes:

Weight Value (Wv) scoring

Let P_i and P_j be two interacting proteins where P_i is the positional candidate that is being considered for prioritization and P_j is the interacting partner. This interaction is assumed to be a result of interaction between their domains D_i and D_j respectively. Let D_i be the parent domain and D_j the partner, where j = 1, 2 ... N, as more than one domain can interact with D_i.

The weight value for D_i was calculated as follows:

(1)

Where,

a_j = D_ja/T_a

b_j = D_jb/T_b

c_j = P_jc/T_c

d_j = P_jd/T_d

N = Total number of partner domains (D_j interacting with parent domain D_i)

Where,

D_ja and D_jb indicate the presence/absence of D_j in classes a and b respectively.

P_jc and P_jd represent the presence/absence of P_j in classes c and d respectively. Note that presence/absence is scored by a binary system.

However, it is more likely that more than one protein can contain D_j and all the proteins harboring D_j needs to be evaluated for categories c and d.

T_a, T_b, T_c and T_d are the total number of events in the classes 'a' through'd' which are 60, 15, 64 and 1746 respectively

Obtained Wv was normalized as follows:

(2)

Subsequently, this Wv obtained for D_i is allocated to P_i as final outcome to rank the candidate genes dataset. For protein harboring more than one domain, highest Wv obtained amongst all the domains was assigned to the protein. Protein with Wv > 0.50 were called as HWEs.

Ranking by clustering coefficient- HRC method

Proteins in HWEs were further ranked by determining network measure i.e. clustering coefficient (C_i) in the domain interaction network (HRC). The clustering coefficient (C_i) is 1 (high) when all the neighbors of a protein are linked to each other and small (C_i = 0) if network is locally sparse [12]. Feldman et al reported that the two networks i.e. yeast two-hybrid (Y2H) and GeneWays (GW) used in their work are complementary rather than competing views of the human interactome [9]. Disease genes are known to avoid dense clustering neighborhoods and the average clustering coefficient for disease network observed was 0.015 [9]. Therefore, we have considered the lowest C_i for the peripherality feature of the interaction network.

Therefore, HWEs were further ranked by using C_i < 0.015 hypothesizing that those candidates below this C_i could be probable disease candidates classifying them as HWEc. InterDom domain interactions (with a confidence score of interactions ≥10) were used for calculating the C_i in NetworkAnalyzer [50].

Benchmarking of prioritization

To assess the performance of Wv and HRC methodologies for prioritizing T2D related genes, we constructed a test dataset consisting: (i) already identified genes implicated in etiology of T2D by Genome wide association studies (n = 19) [1, 15] and (ii) proteins encoded by genes located within 5 Mb region either side of 12 markers (n = 353) that are not linked to T2D (non-T2D genes). Further we analyzed this dataset by four other known independent sequence based methods viz. PROSPECTR [16], SUSPECTS [17], DGP [18] and G2D [19]. PROSPECTR uses sequence based features like gene length, protein length and percentage identity of homologs in other species and designates candidates as likely disease candidates if scored equal to or over a threshold of 0.5 [16]. SUSPECTS method is based on the annotation data from GO, InterDom and expression libraries [17]. Since, SUSPECTS does not define the cut off score for selecting disease genes, we have used a weighted-score of ≥10 based on the training dataset of genes used by Tiffin et al [3]. DGP assigns a probability score using a decision tree model based on sequence properties i.e. protein length, phylogenetic extent, degree of conservation and paralogy [18]. A probability score of ≥0.5 is assigned to all the disease proteins in the learning set used by DGP. Thus, candidates with a score of ≥0.5 were considered as prioritized disease genes in our dataset. G2D method scans a human genomic region for genes related to a particular disease based on the phenotype of the disorder or their similarity to an already known related gene [19]. The method performs BLASTX of the region against all the GO annotated in Refseq database and extracts putative genes in the chromosomal region and evaluates their possible relation to the disease. The default E-value (<10e^-10) was used here for extracting the putative disease genes based on the "known gene" protocol of G2D. The "known gene" takes already well-known disease genes to be associated with disease of interest for prediction. Thus, here we have considered the same set of "known gene" in T2D used by the GeneWanderer program disease-gene families dataset [5]. As G2D takes chromosomal region for search, we have provided 10 Mb either side of 108 markers for searching putative genes. Measure of performance of methods was done by estimating receiver operating characteristic (ROC) curve, which plots true positive rate (TPR) versus false positive rate (FPR).

Comparison with GeneWanderer

We also compared our prediction using 19 identified genes in T2D with recent GeneWanderer program. This approach ranks each candidate gene in a genomic interval (identified by linkage analysis) by their interactions to genes known to be involved in the phenotype/disease in the protein-protein interaction network [5]. Here, we have used the random walk method of GeneWanderer prioritization approach.

Enrichment ratio

Enrichment ratio was estimated as the proportion of disease genes predicted by the methods divided by the disease genes within all benchmarking set as:

(3)

Where ΣDG is all disease genes and ΣG is total number of genes considered here.

Network analysis