A supervised weighted similarity measure for gene expressions using biological knowledge

ARTICLE IN PRESS GENE-41595; No. of Pages: 11; 4C: Gene xxx (2016) xxx–xxx

Contents lists available at ScienceDirect

Gene journal homepage: www.elsevier.com/locate/gene

Research paper

A supervised weighted similarity measure for gene expressions using biological knowledge Shubhra Sankar Ray a, b , Sampa Misra a, * a b

Machine Intelligence Unit, Indian Statistical Institute, Kolkata 700108, India Center for Soft Computing Research, Indian Statistical Institute, Kolkata 700108, India

A R T I C L E

I N F O

Article history: Received 11 January 2016 Received in revised form 18 August 2016 Accepted 22 September 2016 Available online xxxx Keywords: Gene expression Supervised similarity measure Saccharomyces cerevisiae Gene annotation Computational biology Bioinformatics Pattern recognition

A B S T R A C T A supervised similarity measure for Saccharomyces cerevisiae gene expressions is developed which can capture the gene similarity when multiple types of experimental conditions like cell cycle, heat shock are available for all the genes. The measure is called Weighted Pearson correlation (WPC), where the weights are systematically determined for each type of experiment by maximizing the positive predictive value for gene pairs having Pearson correlation greater than 0.80. The positive predictive value is computed by using the annotation information available from yeast GO-Slim process annotations in Saccharomyces Genome Database (SGD). Genes are then clustered by k-medoid algorithm using the newly computed WPC, and functions of 135 unclassified genes are predicted with a p-value cutoff 10 −5 using Munich Information for Protein Sequences (MIPS) annotations. Out of these genes, functional categories of 55 gene are predicted with pvalue cutoff greater than 10 −10 and reported in this investigation. The superiority of WPC as compared to some existing similarity measures like Pearson correlation and Euclidean distance is demonstrated using positive predictive (PPV) values of gene pairs for different Saccharomyces cerevisiae data sets. The related code is available at http://www.sampa.droppages.com/WPC.html. © 2016 Elsevier B.V. All rights reserved.

1. Introduction The process by which genetic information is used to synthesize gene products is called gene expression (Brenner et al., 1961). These products are usually proteins which perform essential functions such as enzymes, hormones and receptors (Lehninger et al., 1993). Genes that do not code for proteins such as ribosomal RNA, micro RNA or transfer RNA are called non-coding RNAs (Ray and Maiti, 2015; Ray et al., 2013). For understanding the functional similarity among genes, one of the ways is to measure their expression similarity (Eisen et al., 1998) which can be obtained from microarrays. Thousands of gene expressions can be monitored simultaneously by using this powerful technique. Moreover, many researchers conducted various microarray experiments, under different experimental conditions, for the same species (e.g., yeast Saccharomyces cerevisiae) to identify the relationship among genes. These diversified experiments help in finding relationships among genes where genes with similar expression profiles are considered as co-regulated or functionally linked and can be helpful in identifying cellular mechanisms (Lehninger et al., 1993) and gene regulatory networks (Hecker et al., 2009; Lopes et al., 2014). There are various methods for finding the similarity between genes like Pearson correlation (Biedl et al., 2001; Bussemaker et al., 2001), Euclidean distance (Slonim, 2002), and Manhattan distance/Cityblock distance (Krause, 2012; Ray et al., 2007). Mahalanobis Distance (Maesschalck et al., 2000) can also be used for finding the similarity. A similarity measure, called representative distance (RD), to identify cancer samples in gene expression data is described in Yu et al. (2012). RD does not consider all the genes in the microarray data, but only a subset of representative genes and uses the affinity propagation algorithm to detect the subset of representative genes. To identify the relation between genes, a neural network based similarity measure has been suggested in Sawa and Ohno-Machado (2003). The problem of identifying time shifted and inverted gene expression similarities is addressed in Qian et al. (2001). All of the above measures are numerical measures in the sense that they depend on the real values of the expressions. An edge detection based similarity measure for periodic data sets with small sequences is proposed by Filkov et al. (2002). This method looks for local regions in pairs of expression profiles where major changes in expression occur (edges). Balasubramaniyan et al. (2005) introduced a similarity measure to identify locally similar and time-shifted patterns Abbreviations: WPC , Weighted Pearson Correlation; SGD , Saccharomyces Genome Database; MIPS , Munich Information for Protein Sequences; PPV , Positive Predictive Value. * Corresponding author. E-mail addresses: [email protected] (S. Ray), [email protected] (S. Misra).

http://dx.doi.org/10.1016/j.gene.2016.09.033 0378-1119/© 2016 Elsevier B.V. All rights reserved.

Please cite this article as: S. Ray, S. Misra, A supervised weighted similarity measure for gene expressions using biological knowledge, Gene (2016), http://dx.doi.org/10.1016/j.gene.2016.09.033

ARTICLE IN PRESS 2

S. Ray, S. Misra / Gene xxx (2016) xxx–xxx

in gene expression profiles where genes are regarded as similar if they have similar edges. In our previous work (Ray et al., 2007) the distance between two genes is computed using a normalization factor, dependent on the dynamic range of photo multiplier tube and application of it is shown in the gene ordering problem. A huge number of investigations are carried out based on the work of Zhang and Horvath (2005) where gene co-expression networks are identified through weighted Pearson correlation. Those weight values are considered that lead to a network satisfying scale-free topology and ensuring mean connectivity for genes to be high enough so that the network contains enough information for getting functionally enriched clusters (modules) in a later stage. However, functional annotations are not used in the methodology for estimating weights and even the weights are not optimized for a trade-off between a high scale free topology fit and the mean number of connections. Note that, scale-free topology is thought to characterize the evolution of biological systems (Barabási and Albert, 1999) but it is not associated with functional annotations of genes. There are also investigations where weights of Pearson correlation are estimated through higher expression similarity in terms of fitness function in quadratic optimization (Teng and Chan, 2006) and experiment related p-values using gene expression (Engreitz et al., 2010) but again in those works weights are not estimated using functional annotations. Further, these methodologies are mainly designed for gene expressions related with diseases (say normal vs. cancer) where the objective is feature selection (to select a subset of patients/features) and thereby ranking the genes according to their relevance to the disease. In this investigation we deal with gene expressions related with cellular functions where the goal is to design a better similarity measure by assigning higher weights to Pearson correlations for experiment types those show better similarity in terms of functional annotations. The popular measures for finding similarity between genes, using gene expression, are Pearson correlation and Euclidean distance, a statistical measure of linear dependence between random variables. The Pearson correlation is used by Spellman et al. (1998) as a similarity measure between gene expression profiles. Unlike the Euclidean distance, Pearson correlation measures the correlation between two genes. The range lies from −1 to +1. Here 1 indicates that the genes are perfectly correlated, and −1 means that the genes are inversely correlated. Let A = a1 , a2 , .....ai ....., an and B = b1 , b2 , .....bi ....., bn be the expression levels of two genes in terms of log-transformed microarray gene expression data. The Pearson correlation between genes A and B is defined as PA,B =

Cov(A, B) , sA sB

(1)

where Cov(A, B) is the covariance, s a and s b are standard deviation of the gene A and gene B. These are defined in Eqs. (2), (3), and (4) as follows: n   1  Cov(A, B) = (2) Ai − A¯ Bi − B¯ , n i=1

  n  2 1  sA =  Ai − A¯ , and n

(3)

i=1

  n 1  ¯ 2. (Bi − B) sB =  n

(4)

i=1

Here A¯ and B¯ are the means of gene A and gene B, respectively, and n is the total number of expression values for both the genes. Substituting Eqs. (2), (3), and (4) in Eq. (1) we get   n  ¯ Bi − B¯ i=1 Ai − A (5) PA,B = 2  2 .  n  n Ai − A¯ Bi − B¯ i=1

i=1

However, Pearson correlation always assume linear relationship (Song et al., 2012) and value of this measure is affected by the extreme value of expressions (Ray et al., 2007). In many cases, negative correlation values in gene expression domain are hard to interpret for gene function prediction task. A weighted similarity measure, called weighted correlation (WC), between the aforementioned two genes A and B with a weighting vector w, is suggested in Zhou et al. (2005). It is defined as     wi Ai Bi − (1/ wi ) wi Ai wi Bi WC(A, B) = (6)  .

      2 (wi Ai )2 wi B2i − (1/ (wi )) ( (wi Bi )) wi A2i − (1/ (wi )) In Eq. (6) wi is a weighting factor depended on the number of expressions (samples) in each experiment type and the summations are over all the expressions i.e., from 1 to n. Suppose a data set consists of 80 expressions and among them 60 expressions belong to cell cycle, 13 expressions belong to sporulation and 7 expressions belong to diauxic shift experiments. Then each of the 60 cell cycle samples are assigned a weight of 1/60. Similarly, the weights for sporulation and diauxic shift experiments are assigned as 1/13 and 1/7, respectively. Song et al. (2012) used biweight midcorrelation (BICOR) to find similarity between two genes. This correlation is originally proposed in Wilcox (2012) where it is shown to be more robust to outliers than Pearson correlation. Considering the expression profiles of the aforementioned genes A and B, BICOR is defined as n (A) (B) i=1 (Ai − med(A)) wi (Bi − med(B)) wi BICOR(A, B) = (7)   n (A) 2 n (B) 2 j=1 Aj − med(A))wj k=1 Bk − med(B))wk (A)

where wi

2

= 1 − u2i I(1 − |ui |) and ui =

Ai −med(A) . 9MAD(A)

Please cite this article as: S. Ray, S. Misra, A supervised weighted similarity measure for gene expressions using biological knowledge, Gene (2016), http://dx.doi.org/10.1016/j.gene.2016.09.033

ARTICLE IN PRESS S. Ray, S. Misra / Gene xxx (2016) xxx–xxx

3

Here, med(A) and MAD(A) are the median and median absolute deviation of A, respectively. If 1 − |ui | > 0 then I(1 − |ui |) takes on value 1, otherwise 0. Song et al. (2012) also compared BICOR with Spearman Rank correlation (SRC) (Hardin et al., 2007), mutual information (MI) (Daub et al., 2004) and Pearson correlation. The mutual information between two genes is calculated by discretizing the expression values using the equal width discretization method (Daub et al., 2004). This method partitions the interval between the minimum expression value of a gene and the maximum expression value of the same gene into equal-width bins. Here, the number of bins is calculated by taking the root of the total number of expressions. So, the MI estimation methods involve parameter choice. Further, estimating joint density to calculate mutual information between two gene expression requires larger sample size. Most of the measures discussed so far, except BICOR and WC, are not ideal to handle gene expressions when multiple types of experimental conditions like cell cycle, heat shock, and diauxic shift are available for the same set of genes. The similarity measures like Pearson correlation and Euclidean distance cannot capture the variation of ranges for different experimental conditions. These measures are unsupervised and cannot take advantage of the existing functional annotation information available for many genes. Moreover, gene expression data are generally noisy and not very reliable because many factors affect the measurements during experiments (Yang et al., 2001). In this situation, incorporation of biological knowledge (functional annotation) in the similarity measure should take care of the above-mentioned limitations to some extent. In this regard, we introduce a supervised similarity measure, called Weighted Pearson Correlation (WPC), for gene expression profiles in this investigation. The supervision is incorporated through weights which are determined systematically by using functional annotation of known genes available from Saccharomyces Genome Database (SGD) (Dwight et al., 2002). Our measure (WPC) is compared with other similarity measures in terms of Positive Predictive Values (PPV) at various similarity values and PPVs for top gene pairs. Finally WPC is evaluated in terms of predicting gene function. To the best of our knowledge there is no work where the weights of Pearson correlations, obtained from different gene expression experiments, are estimated using functional annotations (as available in GO or MIPS) for designing a similarity measure. While the works in weighted correlation (Zhou et al., 2005) and biweight midcorrelation (Song et al., 2012) assign weights to each expression using the number of expressions in each experiment and median absolute deviation, respectively, in our proposed WPC the weights are estimated using a novel functional annotation based technique. There are also investigations where different types of cluster/module refinement techniques using functional annotations are used in the final stage but in our work functional annotations are used to update the weights at the similarity level which is the first and basic step for any gene expression analysis method. Some similarity measures based on gene annotation and relationship of gene annotation with gene expression are analyzed in Sevilla et al. (2005). Here, rather than exploring the expression similarities between two genes, various possibilities of GO function similarity (Al-Mubaid and Nagar, 2008) between them are explored using path length of annotation (Nagar and Al-Mubaid, 2008), GSesame annotation (Wang et al., 2007) and the techniques described by Resnik (1995) and Lin (1998). Further, correlation coefficients between semantic similarity and expression similarity are also studied using Pearson correlation. A similar type of study is also conducted in Wang et al. (2004) where it is pointed out that there is a strong relation between GO-driven similarity and expression correlation of genes. In (Mistry and Pavlidis, 2008), term overlap (TO) based on Gene Ontology terms is used as a measure of gene functional similarity. A clustering method called CliXO (Clique Extracted Ontology) is developed in Kramer et al. (2014) which uses existing ontology based similarity for genes. However, our investigation involving WPC is focused on measuring expression similarity and hence comparisons with aforementioned ontology based (or functional) similarity measures are not feasible. The rest of the investigation is organized as follows. In Section 2, data sets, the way of scoring PPVs from similarity, proposed WPC involving estimation of weights to maximize PPV, comparison of WPC with other similarity measures and the utility of integrating various experiments are described. Section 3 deals with the experimental results for WPC in terms of PPV values for top gene pairs and function prediction of classified and unclassified genes. Section 4 presents the cross validation results using different training and test sets. In Section 5 the important observations are discussed. Finally, Section 6 concludes the investigation. 2. Methods We mainly focus on assigning weights to different types of experiments in a data set using SGD annotations. The main steps of our methodology can be summarized as follows: a. Assigning initial weights for different experimental conditions, separately, in gene expression data set. b. Systematically determine the values of weights by maximizing Positive Predictive Values (PPV) for gene pairs with new similarity value above 0.8. Before describing these steps of the proposed methodology in Section 2.3 in details, we provide a brief description of data sets used in this investigation and the way to calculate PPV at a particular similarity value. Finally functions of some unclassified genes are predicted in Section 3.2 by clustering genes using k-medoids algorithm on the new similarity measure. 2.1. Description of data sets All Yeast (Eisen et al., 1998), Yeast Complex (Bar-Joseph et al., 2001), and Cell Cycle (Sherlock et al., 2001) data sets are used in this investigation. Table 1 shows the name of the data sets, the number of genes in each data set, number of gene functional categories in

SGD, name of experiment types, ranges of expression values, number of experiments performed under each type, and total number of experiments performed for a particular data set. In All Yeast data set three types of experiments cell cycle, sporulation, and diauxic shift are performed. In Yeast Complex and Cell Cycle four types of experiments cell cycle, sporulation, shock, and diauxic shift are performed. Brief description of these experiments are as follows:

• Cell cycle (Cheng et al., 2013): The cell cycle experiment involves in growth, replication, and division of cells. It can be divided into interphase and mitosis. First, cells monitor their environment and grow (synthesize RNA and proteins) if conditions are favorable. Then the DNA synthesis takes place in the cells to replicate the chromosomal DNA. Finally, cells prepare for mitosis. In mitosis, processes like dissolution of the nuclear membrane, spindle formation etc. take place. • Sporulation (Neiman, 2011): The sporulation experiment of the budding yeast Saccharomyces cerevisiae involves starvation of cells from nitrogen and carbon which activate coordinated and sequential changes in gene expression leading to the production of haploid and stress-resistant spores. • Heat shock (Castells-Roca et al., 2011): The heat shock experiments deal with the most common stresses like elevated temperature and oxidative stress on the cell. The cellular response

Please cite this article as: S. Ray, S. Misra, A supervised weighted similarity measure for gene expressions using biological knowledge, Gene (2016), http://dx.doi.org/10.1016/j.gene.2016.09.033

ARTICLE IN PRESS 4

S. Ray, S. Misra / Gene xxx (2016) xxx–xxx Table 1 Summary for different microarray data sets. Dataset

All Yeast

Yeast Complex

Cell Cycle

No. of genes Category Experiments performed

6221 SGD (38) Cell cycle (−1.2 to 1.2) 60 Sporulation (−3.0 to 3.0) 13 Diauxic shift (−2.0 to 2.0) 7

Total

80

979 SGD (16) Cell cycle (−1.2 to 1.2) 18+14+15 Sporulation (−3.0 to 3.0) 7+4 Shock (−1.5 to 1.5) 6+4+4 Diauxic shift (−2.0 to 2.0) 7 79

652 SGD (16) Cell cycle (−1.2 to 1.2) 93 Sporulation (−3.0 to 3.0) 9 Shock (−1.5 to 1.5) 56 Diauxic shift (−2.0 to 2.0) 26 184

to heat shock is generally evident from the transcriptional upregulation of genes and encoding of heat shock proteins (HSPs) as a part of the cell’s internal repair mechanism. • Diauxic shift (Galdieri et al., 2010): In diauxic shift experiment, yeast cells are kept in a glucose solution and ideal temperature (say 30 ◦ C) is maintained. While during lag phase the number of cells will not change, in exponential or log phase the number of cells will grow exponentially. Initially it will take the yeast cells some time to realize that they are in a nutrient rich environment. Once realized, they will start rapid cell division to use this rich source of energy. The used data sets of Saccharomyces cerevisiae consist of about 6221, 979, and 652 genes and 80, 79, and 184 experiments. Missing gene expression values of each data set are predicted with LSimpute (Hellem et al., 2004). The genes of the data sets are classified into 38, 16, and 16 groups, respectively, according to yeast GO-Slim process annotations, downloaded on May 2015 from SGD. The names of these groups are mentioned in Table 2. The general information about the genes belonging to different functional annotations for each data set is described in web page links as xls files which can be used for reproducing our results by constructing the annotation vectors through steps as mentioned in our web page http://www.sampa.droppages. com/WPC.html. The links about gene information for different data sets are as follows: a) All Yeast data set: https://drive.google.com/file/d/0BxposWdIam0UUzZFcktWYlJj eTA/view?usp=sharing b) Yeast Complex data set: https://drive.google.com/file/d/0BxposWdIam0UZ2ZYYVJLZ3Z vSmM/view?usp=sharing c) Cell Cycle data set: https://drive.google.com/file/d/0BxposWdIam0UZmN4ZDZqM HF2N00/view?usp=sharing As our ultimate goal is to predict gene functions, rather than location or molecular activity, we have chosen yeast GO-Slim process annotations. Other types of annotations such as yeast GO-Slim component will lead to prediction of locations of genes rather than their involvement in biological process and hence they are not used in our investigation. 2.2. Scoring positive predictive values from similarity In this section, we describe the way of computing positive predictive value (PPV) at a particular similarity value of gene pairs. The

proportion of true positive (TP) gene pairs at various similarity values (for any measure) can be used to compute the PPV, where TP gene pairs are defined as pairs of genes having at least one overlapping yeast GO-Slim process term (Dwight et al., 2002). It may be noted that in yeast GO-Slim process, every gene is annotated in the same level. The PPV (proportion of TP pairs) at a particular similarity value is defined as (Troyanskay et al., 2003)

PPV =

No. of predicted pairs with at least one common GO term . Total no. of predicted pairs (8)

If any of the genes in a particular pair does not belong to any functional category, then that pair is not considered in the total number of predicted pairs. 2.3. New similarity measure As mentioned earlier, the existing measures are still not perfect in reflecting the actual similarity between two gene expressions in the presence of different experimental conditions and cannot utilize the existing biological knowledge in reflecting the functional similarity between genes. In this regard, a supervised measure, called Weighted Pearson Correlation (WPC), is developed where for each type of experiment the weights are determined by maximizing the positive predictive value, using yeast GO-Slim process annotations in Saccharomyces Genome Database (SGD), for all gene pairs having Pearson correlation greater than 0.80. Let us consider, e

e

e

e

e

e

A = a11 , .....ai 1 , a12 , .....ai 2 , ....................., ae1m , .....aeimm and 1

e

e

2

B = b11 , .....bi 1 , b12 , .....bi 2 , ......................, be1m , .....beimm . 1

2

as the expression vectors of the two genes A and B, respectively, over a series of m different types of experiment (e1 , e2 , . . . , em ) consisting of i1 + i2 + . . . + im experiments in total. The new supervised weighted Pearson correlation between A and B is defined as

WPC(A, B) =

m  r=1

  ir   Aj − A¯ Bj − B¯

kr m 

j=1

× .   i  2  2 ir  r   kr ¯ ¯   Aj − A × Bj − B r=1 j=1

(9)

j=1

Please cite this article as: S. Ray, S. Misra, A supervised weighted similarity measure for gene expressions using biological knowledge, Gene (2016), http://dx.doi.org/10.1016/j.gene.2016.09.033

ARTICLE IN PRESS S. Ray, S. Misra / Gene xxx (2016) xxx–xxx

5

Table 2 Summary of different annotation groups in datasets. Dataset

Function name

All Yeast

1. Organelle organization and biogenesis 2. Transport 3. RNA metabolic process 4. DNA metabolic process 5. Transcription 6. Protein modification process 7. Response to stress 8. Cell cycle 9. Response to chemical stimulus 10. Translation 11. Ribosome biogenesis and assembly 12. Vesicle-mediated transport 13. Anatomical structure morphogenesis 14. Lipid metabolic process 15. Cytoskeleton organization and biogenesis 16. Signal transduction 17. Carbohydrate metabolic process 18. Amino acid and derivative metabolic process 19. Cell wall organization and biogenesis 20. Membrane organization and biogenesis 21. Protein catabolic process 22. Generation of precursor metabolites and energy 23. Cofactor metabolic process 24. Meiosis 25. Cellular homeostasis 26. Sporulation 27. Conjugation 28. Cytokinesis 29. Vitamin metabolic process 30. Cellular respiration 31. Cell budding 32. Heterocycle metabolic process 33. Protein folding 34. Aromatic compound metabolic process 35. Pseudohyphal growth 36. Nuclear organization and biogenesis 37. Transposition 38. Electron transport 1. Metabolism 2. Energy 3. Cell cycle and DNA processing 4. Transcription 5. Protein synthesis 6. Protein fate 7. Protein with binding function or cofactor requirement 8. Protein activity regulation 9. Cellular transport, transport facilitation and transport routes 10. Cellular communication/signal transduction mechanism 11. Cell rescue, defense and virulence 12. Interaction with the cellular environment 13. Cell fate 14. Development 15. Biogenesis of cellular components 16. Cell type differentiation 1. Metabolism 2. Energy 3. Transcription 4. Protein synthesis 5. Protein fate 6. Protein with binding function or cofactor requirement 7. Protein activity regulation 8. Cellular transport, transport facilitation and transport routes 9. Cellular communication/signal transduction mechanism 10. Cell rescue, defense and virulence 11. Interaction with the cellular environment 12. Transposable elements, viral and plasmid proteins 13. Cell fate 14. Development 15. Biogenesis of cellular components 16. Cell type differentiation

Yeast Complex

Cell Cycle

The following can be stated about the measure:

⎞

⎛ ⎜ In Eq. (9) each k ⎝say

kr m 

kr

⎟ ⎠ is applied to all possible values of

r=1

Pearson correlation obtained from a particular (say, r) experiment type. If there are n genes then there are n × n possible Pearson correlations (corresponding to gene pairs) for an experiment. Example: Let us consider two genes A and B with two types of experiments like cell cycle and diauxic shift, and two different weights k1 = 1 and k2 = 2 are assigned for them. If the Pearson correlations between two genes A and B for these two types of experiments are 0.5 and 0.4 then the weighted Pearson correlation (WPC) between two genes A and B is (1 ∗ 0.5)/(1 + 2) + (2 ∗ 0.4)/(1 + 2) = 0.43 Estimation of weights: In WPC measure, the number of weights is equal to the number of experiments performed for each data set, i.e., if there are three experiments like cell cycle, sporulation, and diauxic shift, then three different weights k1 , k2 , and k3 are assigned for them, respectively. The weights for different experiments in WPC is estimated by the following steps: S1) All the weights are initialized as value 1. S2) Similarity values are calculated for all the gene pairs to identify the pairs having similarity above 0.8. S3) PPV is calculated for those gene pairs. S4) The weighting factors are now varied in steps of 0.1 and the steps from 2 to 3 are repeated to find a combination of weights for which the PPV is maximized. Fig. 1 shows a typical example where PPVs, using yeast GO-Slim process annotations, of All Yeast data vary for different values of weight factors ranging from 0 to 100 in steps of 1. Note that the intermediate steps (0.1) are not shown in the figure for better clarity. The curves show instances where one weight factor is varied and the other weight factors are kept constant. The final values of weights are k1 = 1, k2 = 18.5, and k3 = 3.5, where the PPV is maximized.

Positive Predictive Value−−−−−−−>

1) −1 ≤ WPC(A, B) ≤ 1; 2) WPC(A, B) = 0 if and only if A = B; 3) WPC(A, B) = WPC(B, A) (symmetric).

0.99

0.98

0.97

0.96

0.95

0.94 varying k (weight of cell cycle experiment) k =18.5 k =3.5 1

0.93

2

3

varying k (weight of sporulation experiment) k = 1 k = 3.5 2

1

3

varying k (weight of diauxic shift experiment) k =1 k =18.5 3

0.92

0

20

1

40

60

80

2

100

Value of Weight for All Yeast data source−−−−−−−> Fig. 1. Comparing the values of PPV using WPC for All Yeast data, by varying weights of different experiments. The PPVs are computed for similarity values above 0.8. When a particular weight is varied, the other weights are kept constant at the values shown in the figure.

These values (k1 = 1, k2 = 18.5, and k3 = 3.5) indicate that the Pearson correlations of the second and third type of experiments are getting 18.5 and 3.5 times more weights, respectively, than the first type. These values maintain the balance such that the experiment types with higher number of true positive (functionally similar) gene pairs at higher Pearson correlation values will get more weight and the weight value depends on that number. In other words, the weight values indicate the capability of an experiment type in reflecting the functional gene similarity, in the yeast GO-Slim process annotation space, in accordance with the Pearson correlation. Experiment type, where higher correlation values show more functional annotation similarity, should be assigned higher weight than the other ones. For All Yeast data, high value of k2 for sporulation experiment (second one) indicates that at high values of Pearson correlation the degree of accuracy (in terms of functional similarity) for this experiment is

Please cite this article as: S. Ray, S. Misra, A supervised weighted similarity measure for gene expressions using biological knowledge, Gene (2016), http://dx.doi.org/10.1016/j.gene.2016.09.033

ARTICLE IN PRESS S. Ray, S. Misra / Gene xxx (2016) xxx–xxx

Pearson Correlation Proposed WPC 0.03 Weighted Correlation 0.11 Spearman Rank Correlation Biweight midcorrelation 0.14 Mutual Information Euclidean Distance 0.14 Manhattan Distance Mahalanobis Distance 0.07

1

PPV−−−−−>

0.8

0.6

0.4

1

0.8

0.03

0.2

0

Pearson Correlation Proposed WPC Weighted Correlation Spearman Rank Correlation Biweight midcorrelation Mutual Information Euclidean Distance Manhattan Distance Mahalanobis Distance

0.02

PPV−−−−−>

6

0.06 0.14 0.15 0.13

0.08 0.07 0.05

0.6

0.4

0

0.2

0.4

0.6

0.8

0.2

1

Similarity value−−−−−> 0

Fig. 2. The PPVs versus the similarity values are plotted for different similarity measures for All Yeast data set using yeast GO-Slim process annotations. The vertical lines between WPC and PC show the differences between them in terms of PPV at a particular similarity value.

much better than other type of experiments. For other data sets the estimated weight combinations are as follows: i) Yeast Complex: k1 = 1, k2 = 12.8, k3 = 1.3, and k4 = 1 ii) Cell Cycle: k1 = 1, k2 = 1, k3 = 18.8, and k4 = 7.6

0

0.2

0.4 0.6 Similarity value−−−−−>

0.8

1

Fig. 3. For different similarity measures PPVs versus the similarity values are plotted for Yeast Complex data set using yeast GO-Slim process annotations. The vertical lines between WPC and PC show the differences between them in terms of PPV at a particular similarity value.

(within a cluster) which are used to predict functions in the later stage.

2.4. Mapping similarity value to PPV for different data sources

2.5. Utility of integrating different types of experiment

In this section the performance of WPC is compared with other similarity measures in term of PPVs obtained at various similarity values. The performance of WPC for All Yeast, Yeast Complex, and Cell Cycle data sets are shown in Figs. 2, 3, and 4, respectively, using yeast GO-Slim process annotations. The weights of WPC for various data sets are kept constant at values as mentioned in Section 2.3. Here all the measures are scaled within the range 0 to 1. In case of Euclidean, Manhattan and Mahalanobis distance the values are inverted after normalization. From the figures it is clear that PPVs increase for higher similarity values of gene pairs. From Figs. 2, 3, and 4 it can be observed that the performance of WPC is better than the other similarity measures above similarity values 0.5, 0.3 and 0.3 for All Yeast, Yeast Complex and Cell Cycle data sets, respectively. It can also be observed from these figures that the performance of biweight midcorrelation (BICOR) is slightly better than Pearson correlation for all the data sets. While the curves for Spearman rank correlation (SRC) and weighted correlation (WC) suggest that their performance is slightly superior or comparable to Pearson correlation in the similarity value range 0.5 to 0.9, all other measures are noticeably inferior to them considering all the similarity values. From Fig. 4 it can be noticed that at some similarity values around 0.8, Euclidean distance performs better than all other measures for Cell Cycle data set, which is the only case which shows inferior performance of WPC. In Figs. 2, 3, and 4 some vertical lines are plotted to show the difference between the proposed WPC and PC in terms of PPV at various similarity values. This clearly represents the advantage of the WPC above 0.5 similarity value. For example, considering 0.75 and 0.85 similarities in Fig. 2 (All Yeast data and yeast GO-Slim process annotations), the PPVs for WPC are 0.64 and 0.91, respectively, as compared to 0.5 and 0.8 for Pearson correlation. Note that, genes with highly similar expression values mainly group together

The basic idea of integrating different types of experiments in any domain is to enhance the accuracy of the integrated output as compared to the accuracy of its individual components. In this regard, we have now performed experiments by considering each experiment type separately and the PPV vs. the similarity plots for them are

Pearson Correlation Proposed WPC Weighted Pearson Correlation Spearman Rank Correlation Biweight mid correlation Mutual information Euclidean Distance Manhattan Distance Mahalanobis Distance

1

PPV−−−−−>

0.8

0.6

0.02

0.06

0.08

0.01

0.05

0.4

0.2

0

0

0.2

0.4

0.6

0.8

1

Similarity value−−−−−> Fig. 4. The PPVs versus the similarity values are plotted for different similarity measures for Cell Cycle data set using yeast GO-Slim process annotations. The vertical lines between WPC and PC show the differences between them in terms of PPV at a particular similarity value.

Please cite this article as: S. Ray, S. Misra, A supervised weighted similarity measure for gene expressions using biological knowledge, Gene (2016), http://dx.doi.org/10.1016/j.gene.2016.09.033

ARTICLE IN PRESS S. Ray, S. Misra / Gene xxx (2016) xxx–xxx

0.8

Pearson Correlation Weighted Pearson Correlation Weighted correlation Spearman Rank Correlation Biweight midcorrelation Mutual information Euclidean Distance Manhattan Distance Mahalanobis Distance

0.9

0.8

0.7

0.6

PPV−−−−−−>

PPV−−−−−>

1

PC for All Yeast data set WPC for All Yeast data set PC for first experiment type (cell cycle) PC for second experiment type (sporulation) PC for third experiment type (diauxic shift)

1

7

0.4

0.2

0.6

0.5

0.4

0.3

0

0

0.2

0.4

0.6

0.8

1

Similarity value−−−−−>

0.2

Fig. 5. The PPVs versus the similarity values are plotted for PC and WPC. While WPC is plotted by considering all the experiment types of All Yeast data, PC is plotted for the same as well as considering individual experiment types i.e., cell cycle, sporulation and diauxic shift.

0.1

0

2

4

6

8

Number of top relations−−−−−>

10 4

x 10

Fig. 6. The PPVs versus the number of top gene pairs are plotted for different similarity measure using All Yeast data set and MIPS annotations.

shown in Fig. 5. As expected, from the figure it can be seen that the PPV vs. similarity curve of the integrated output (WPC) is superior to any of those individual experiments. This shows the utility, in terms of increase in PPV, in integrating different types of gene expression experiments. These curves also indicate that only one type of experiment will not provide similar results with those obtained using integration of multiple type of experiments as in WPC. Moreover, the curve of second experiment (sporulation) lies above the first (cell cycle) and third (diauxic shift) experiment. This indicates that order of weights obtained using WPC are in accordance with the order of importance of the various types of experiments as evident from Fig. 5. Interestingly, not only WPC but also PC with all experiments shows better performance than PC of individual experiment.

pairs are 0.68, 0.65, 0.61, 0.57, 0.66, 0.21, 0.49, 0.50 and 0.34 for proposed WPC, PC, WC, SRC, BICOR, MI, Euclidean distance, Manhattan distance and Mahalanobis distance, respectively. While, for Yeast Complex data those values are 0.73, 0.69, 0.68, 0.67, 0.70, 0.45, 0.63, 0.63 and 0.55 for the same number of gene pairs, for Cell Cycle data those values are 0.55, 0.52, 0.51, 0.52, 0.52, 0.45, 0.51, 0.51 and 0.51.

1

Pearson Correlation Proposed WPC Weighted Pearson Correlation Spearman Rank Correlation Biweight mid Correlation Mutual Information Euclidean Distance Manhattan Distance Mahalanobis Distance

0.9

3. Results

3.1. PPV of top gene pair Performance of top gene pairs, in terms of PPV using MIPS annotations, for All Yeast, Yeast Complex, and Cell Cycle data (curves for WPC are at the top) is shown in Figs. 6, 7, and 8, respectively. As expected, the PPV decreases when the number of top gene pairs increases. It can be observed that the performance of WPC is better than the other similarity measures, in terms of PPV, for All Yeast, Yeast Complex and Cell Cycle data sets, respectively. Moreover, PC, WPC, WC, SRC and BICOR measures are noticeably better than Euclidean distance, Manhattan distance, Mahalanobis distance and MI. For example, using All Yeast data the PPVs of top 50,000 gene

PPV−−−−−−>

0.8

As SGD yeast GO-Slim process annotations are used to find the weights of WPC, Munich Information for Protein Sequences (MIPS) (Mewes et al., 2002) annotations are used to evaluate the performance for top 105 gene pairs in terms of PPV. Further, MIPS annotations are used to predict the function of unclassified genes from clusters obtained by using k-medoids algorithm on WPC. According to MIPS there are 510 functional categories available for yeast Saccharomyces cerevisiae.

0.7

0.6

0.5

0.4

0

2

4

6

8

Number of top relations−−−−−>

10 4

x 10

Fig. 7. Comparisons among weighted Pearson correlation and other measures in terms of PPVs versus the number of top gene pairs of Yeast Complex data set using MIPS annotations.

Please cite this article as: S. Ray, S. Misra, A supervised weighted similarity measure for gene expressions using biological knowledge, Gene (2016), http://dx.doi.org/10.1016/j.gene.2016.09.033

ARTICLE IN PRESS 8

S. Ray, S. Misra / Gene xxx (2016) xxx–xxx

1 Pearson Correlation Proposed WPC Weighted Correlation Spearman Rank Correlation Biweight midcorrelation Mutual Information Euclidean Distance Manhattan Distance Mahalanobis Distance

PPV−−−−−−>

0.9

0.8

0.7

0.6

0.5

0.4

0

2

4

6

8

Number of top relations−−−−−>

10 4

x 10

Fig. 8. The PPVs versus the number of top gene pairs are plotted for different similarity measure using Cell Cycle data set and MIPS annotations.

3.2. Function prediction of unclassified gene based on clustering results Genes are clustered by using k-medoids algorithm on the proposed WPC. The k-medoids algorithm (Dodge, 2012) is a partitive clustering algorithm (breaking the data set into a number of groups) which attempts to minimize the distance between points labeled to be in a cluster and a point designated as a center of that cluster. The k-medoids chooses centers from data points and can work with a similarity/distance matrix computed on data points instead of the co-ordinates of data points. The first step of k-medoids algorithm is to define k cluster centers one for each cluster from data set and then assign data point to the nearest medoid, after assigning all the data points k new medoids are recalculated. The whole process is repeated until no more changes are observed. Once the genes are clustered, the biological function of an unclassified gene is predicted from the functional enrichment of the cluster using MIPS annotations and a user defined p-value cut off. The methodology and the results are as follows: S1) As there are 510 MIPS functional categories, the k-value (cluster member) is chosen as 510. S2) The clusters are then studied for their functional enrichment and are divided in three sets based on p-value.

• In the first set 40 clusters are identified with p-value less than 10 −5 . From these clusters functions of 987 classified and 135 unclassified genes are predicted. • In the second set 16 clusters are identified with p-value less than 10 −10 . Functions of 417 classified and 55 unclassified genes are predicted from these clusters. • In the third set 2 clusters are identified with p-value less than 10 −16 and functions of 59 classified and 2 unclassified genes are predicted. Table 3 summarizes the predicted functions for 55 unclassified genes from 16 clusters in the second set. Note that for each of the predicted functions, the related p-values, number of genes within the cluster and within the genome is shown in the table. The importance of the predicted functions for some of these unclassified genes in Table 3 is discussed in the next section. The functions of 135 unclassified genes in the first set are available at http://www.isical.ac.in/~

shubhra/WPC/unclassifiedprediction135.xls and those of 55 unclassified genes in the second set are available at http://www.isical.ac. in/~shubhra/WPC/unclassifiedprediction55.xls. It will be now prudent to discuss the advantage of WPC over PC in terms of gene function prediction as the computational burden of WPC is higher than PC. By considering Fig. 2 it can be observed that at the similarity values 0.65, 0.8 and 0.85 the differences between WPC and PC in terms of PPV are 0.07, 0.14 and 0.11, respectively. The effect of these differences can be observed in the clustering results as well as in function prediction of the unclassified genes. For example, consider the unclassified gene YER0561C mentioned in row 1 of Table 3. In this cluster corresponding to row 1, there are 2 unclassified genes (YER0561C and YHR049W) and 23 classified genes belonging to function “Ribosome biogenesis”. The average Pearson correlation between the unclassified gene YER0561C and the classified genes YIL052C, YLR344W, YIL018W, YLR062C, YER064C, YKR094C, YLL044W, YCR031C, YGL076C, YER074W, YGR034W, YEL054C, YBR079C, YDL082W, YDR471W, YJL189W, YPL143W, YLR185W, YHR049W, YNL119W, YDR417C, YDR418W and YLR405W is 0.70. The gene YER0561C would have never formed a cluster with the aforementioned genes using k-medoids as the average Pearson correlation value (0.70) might not be good enough to group YER0561C with the classified ones. On the other hand, the proposed WPC assigned higher correlation values between YER0561C and the classified genes, the average of which is 0.87. This leads to grouping of YER0561C with the aforementioned classified genes in the same cluster using kmedoids. This grouping also enabled us to predict the function of unclassified gene YER0561C as “Ribosome biogenesis” with p-value 2.7403 × 10 −16 . Without WPC, it might have not been possible to predict this function. Two other examples are as follows: i) The function of unclassified gene YPR203W, mentioned in row 2 of Table 3, is predicted as “DNA topology” with pvalue 2.7986 × 10 −16 using WPC. This gene is clustered with classified genes YLR467W, YER111C, YGR296W, YER190W, YPR202W, YPR204W, YJL225C, YDL127W, YIL177C, YHL050C, YLL066C, YBL110C, YPR018W, YCL022C, YDR545W, YPL283C, YLR463C, YLR464W, YLR465W, YLR466W, YBR158C, YEL077C and YNL339C, having function “DNA topology”, with average WPC value 0.77 using k-medoids. In contrast, the average Pearson correlation value between this unclassified gene and other classified genes in the cluster is 0.65. So it might have not been possible to predict the function of the gene YPR203W using Pearson correlation. ii) The function of unclassified gene YNL060C, mentioned in row 6 of Table 3, is predicted as “rRNA processing” with p-value 7.1314 × 10 −13 using WPC. This gene is clustered with classified genes YMR049C, YPR110C, YCL053C, YCL054C, YPL093W, YJL069C, YLR419W, YOR001W, YIL096C, YLR138W, YNL075W, YBR266C, YDL121C, YLR073C, YOR091W, YGL078C, YPL044C, YLR223C, YER046C, YBR069W, YPL126W, YMR146C, YDR395W, YMR290C, YJL109C, YER025W, YLL008W, YCLX02C, YNL112W, YMR310C and YLR409C, having function “rRNA processing”, with average WPC value 0.74 using k-medoids. In contrast, the average Pearson correlation value between this unclassified gene and other classified genes in the cluster is 0.6. So it might have not been possible to predict the function of the gene YNL060C using Pearson correlation. Other cases in Table 3 are also similar in nature.

4. Cross validation The purpose of validation using different training and test sets is served to some extent as we have used yeast GO-Slim process annotations for weight determination and a different annotation

Please cite this article as: S. Ray, S. Misra, A supervised weighted similarity measure for gene expressions using biological knowledge, Gene (2016), http://dx.doi.org/10.1016/j.gene.2016.09.033

ARTICLE IN PRESS S. Ray, S. Misra / Gene xxx (2016) xxx–xxx

9

Table 3 Function predictions of top 55 unclassified genes (according to mips classification) at p-value less than 10 −10 . Row no.

Unclassified gene

Functional category

p-Value

1 2 3

YER0561C YHR049W YPR203W YLR464W YPR202W YBL110C YGL081W YHR079BC YER0441C YKL121W YNL269W YNL211C YLR376C YCLX01W YJR030C YJL018W YNL060C YCLX02C YCL053C YIL096C YDL121C YLR073C YOR091W YMR310C YBL101W-A YFL-TYA YFL-TYB YBR246W

Ribosome biogenesis DNA topology Meiosis

2.7403e-16 16 2.7986e-16 13 9.5217e-11 13

284 52 156

rRNA processing Extension/polymerization activity rRNA processing

1.1504e-15 15 9.1865e-11 7 7.1314e-13 13

167 37 167

4.0313e-13 8

34

9.7631e-15 12

79

2.8927e-15 10

48

3.6401e-15 16

167

1.2227e-11 10 2.6407e-16 19

154 221

9.1519e-17 14

126

1.654e-17 19 9.4146e-12 10

453 126

6.8945e-13 11

126

4 5 6 7 8 9 10 11 12 13 14 15 16

Transposable elements, viral and plasmid proteins YDR233C YCLX06C Sugar, glucoside, polyol and carboxylate catabolism YER158C YJR080C YOR205C YGL107C Electron transport and membraneYOR006C YGR206W associated energy conservation YCR016W YBR271W YIL110W rRNA processing YDR161W YDR198C YMR157C Mitochondrion YLR204W YHR059W YML030W Ribosomal proteins YDR493W YKR074W YLR253W YCL056C Proteasomal degradation (ubiquitin/proteasomal pathway) YMR266W Protein synthesis YOR164C YCLX05C Proteasomal degradation (ubiquitin/proteasomal pathway) YIR016W YBR151W YOL098C Proteasomal degradation (ubiquitin/proteasomal pathway

database, MIPS, is used for evaluation (Figs. 6, 7, and 8) using the same weights. Further, to evaluate in a conventional way 2 fold cross validation is used where half of the genes are randomly selected (along with their yeast GO-Slim process annotation profiles) for the weight determination (k1 = 1.8, k2 = 18.1, and k3 = 4.8) process and the remaining half is used for evaluation using their GO-Slim process annotation profiles. The process is repeated 5 times and the results for one of those instances using All Yeast data set are presented in Fig. 9. From the figure it is clear that the performance of WPC is better than PC using even half of the genes. As expected, the curves obtained using half of the genes in the testing process are below than those obtained using all the genes as the number of gene pairs are reduced to one fourth as compared to the number of gene pairs obtained when all the genes are used. The cross validation results for Yeast Complex and Cell Cycle data sets are available in supplementary file at https://drive.google.com/open? id=0BxposWdIam0UNjdWUUVqSWpWMUk.

No. of genes within the genome

Evaluation of WPC is performed in three widely studied ways (Marcotte et al., 1999; Troyanskay et al., 2003) as follows:

• In the first way, the PPVs of WPC are compared with those of related measures at various similarity values using multiple Saccharomyces cerevisiae data set. Above 0.5 similarity value better performance of WPC is clearly visible from Figs. 2, 3, and 4 for All Yeast, Yeast Complex and Cell Cycle data sets, respectively. • In the second way, the PPVs of WPC are compared with those of related measures for top gene pairs ranging from 1 to 105 . Using the same data sets the superior performance of WPC can be observed from Figs. 6, 7, and 8.

WPC using 2 fold cross validation PC using half Data set WPC using full Data set PC using full Data set

1

5. Discussion

0.8

PPV−−−−−>

In this investigation, three yeast data sets namely All Yeast, Yeast Complex, Cell Cycle are used where multiple type of experimental conditions are available. For a data set with single type of experimental condition any one of the existing measures can be used to compute gene similarity. As discussed earlier these measures are unable to compute the exact similarity between genes with multiple types of experiments as the experiment with higher range (e.g., range of cell cycle is −1.2 to 1.2 as compared to those of sporulation experiments which is −3.0 to 3.0) dominates the lower one. Moreover, existing similarity measures cannot utilize the biological knowledge in reflecting the functional similarity between genes. In this regard, we developed a supervised similarity measure, weighted Pearson correlation (WPC), which will systematically take care of the various issues by incorporating the biological knowledge through weights to various experimental conditions. The weights are supervised by maximizing the PPVs of gene pairs above 0.8 similarity value using SGD annotations.

No. of genes within cluster

0.6

0.4

0.2

0

0

0.2

0.4

0.6

0.8

1

Similarity value−−−−−> Fig. 9. Comparing WPC with PC in terms of PPVs versus similarity values using half of the genes in All Yeast data set. The top two curves are obtained using all of the genes as shown in Section 2.4.

Please cite this article as: S. Ray, S. Misra, A supervised weighted similarity measure for gene expressions using biological knowledge, Gene (2016), http://dx.doi.org/10.1016/j.gene.2016.09.033

ARTICLE IN PRESS 10

S. Ray, S. Misra / Gene xxx (2016) xxx–xxx

• Finally, WPC is evaluated in terms of predicting gene function using MIPS annotation. Here, genes are clustered using k-medoids algorithm on WPC measure and function of 135 unclassified and 987 classified genes is predicted at a p-value cutoff 10 −5 on functionally enriched clusters. Function predictions of top 55 unclassified genes at a p-value less than 10 −10 are already shown in Table 3.

Now we will discuss in details about the predicted function of some of those unclassified genes. From Table 3, we find that functional enrichment of the gene YMR266W is in MIPS category “Protein synthesis”. Our examination revealed that the related protein is of unknown function but its over-expression suppresses NaCl sensitivity of sro7 mutant cells by restoring sodium pump (Ena1p) localization to the plasma membrane (Dwight et al., 2002). Further studies (Vijaranakul et al., 1997) indicate that synthesis of some proteins is increased in response to NaCl stress. So our prediction of YMR266W in “Protein synthesis” is a possible one. The gene YHR219W has shown functional enrichment in “DNA topology”. Our analysis predicts that YHR219W is involved in “DNA topology” because the putative protein of unknown function for this gene is similar to helicases (Shiratori et al., 1999). These proteins play important roles in various cellular processes including DNA replication, DNA repair, and maintenance of chromosomal stability. Hence, our prediction for YHR219W as “DNA topology” is a highly likely one. The gene YMR157C has 20 cluster members and 10 out of 20 genes show functional enrichment in MIPS category “Mitochondrion”. Our literature study (Hess et al., 2009) reveals that it is a null mutant, which displays elevated frequency of mitochondrial genome loss and the corresponding protein is detected in purified mitochondria in high-throughput experiments. So our function prediction for YMR157C as “Mitochondrion” is a promising one. The functional category of YIR016W is predicted as “Proteasomal degradation”. According to SGD this putative protein expression is directly regulated by the metabolic and meiotic transcriptional regulator and over-expression of it causes a cell cycle delay or arrest. The proteasomal degradation pathway is essential for many cellular processes, including the cell cycle and the regulation of gene expression (Pfaff, 2002). Hence, our prediction for category “Proteasomal degradation” may be a possible one. WPC can easily be computed for any species where functional annotations of some genes and expression values for different types of experiments are available. First, one has to construct the functional annotation profile for each gene of the new species and then compute the annotation similarity by checking overlaps between genes using their annotation profiles. In our investigation we have performed this annotation profile construction through a computational approach. In this regard we have now designed a webpage at http://www.sampa.droppages.com/WPC.html where the source code for annotation profile construction as well as the main program for calculating the weights of Pearson correlations, obtained from different experiment types, is available. Once the WPC is computed for any species it can be used for gene function prediction for that species as shown in Section 3.2. If new functions for genes are identified for Saccharomyces cerevisiae then one has to update the annotation profiles for them. For example, say the annotation profile, representing the categories “Cellular respiration”, “RNA splicing”, and “Translation”, for the gene ‘YER054C’ is 0 1 0. If after sometime (say 1 year) this gene, in addition to “RNA splicing”, is also found to be involved in the category “Cellular respiration” then one has to update the profile as 1 1 0. On the other hand, if a new category (say “Vacuole organization”) is discovered for that gene then one have to increase the annotation profile length for all the genes, and only for ‘YER054C’ the corresponding category value will be 1. The new profile will be 0 1 0 1 for ‘YER054C’. The various steps for

annotation profile construction of genes can be easily executed by a programmatic approach as shown in our webpage. 6. Conclusion In this study, we developed a new framework for similarity measure called weighted Pearson correlation (WPC) where the weights are estimated using the functional annotations available from SGD. WPC is then used to predict the function of 135 unclassified and 987 classified Saccharomyces cerevisiae genes using a p-value cutoff 10 −5 . The functions are predicted using MIPS annotations as yeast GO-Slim process annotations are used for finding the weights. The performance of WPC is found to be better than the existing Pearson correlation, Euclidean distance, Manhattan distance, weighted correlation, biweight midcorrelation, Spearman rank correlation, Mahalanobis distance, and Mutual information in terms of PPV at various similarity values and also for top 80,000 gene pairs. Further, the utility of integrating different types of experiments through WPC is also established. Our results indicate that WPC is capable of identifying the similarity between gene expression profiles in a better manner in terms of functional annotations. The WPC is found to be useful in predicting gene functions using k-medoids algorithm. Evaluation of WPC on gene function confirmed its potential value and utility for data sets where multiple types of microarray experiments and functional annotations for a high number of genes are available. The WPC can also be computed for any species where functional annotations of some genes and expression values for different types of experiments are available. Even under the same experimental category the weights can be assigned to various time points for finding functionally similar genes in a supervised manner.

References Al-Mubaid, H., Nagar, A., 2008. Comparison of four similarity measures based on go annotations for gene clustering. IEEE Symposium on Computers and Communications. pp. 531–536. Balasubramaniyan, R., Hüllermeierand, E., Weskamp, N., Kämper, J., 2005. Clustering of gene expression data using a local shape-based similarity measure. Bioinformatics 21 (7), 1069–1077. Bar-Joseph, Z., Gifford, D.K., Jaakkola, T.S., 2001. Fast optimal leaf ordering for hierarchical clustering. Bioinformatics 17 (suppl. 1), S22–S29. Barabási, A.L., Albert, R., 1999. Emergence of scaling in random networks. Science 286 (5439), 509–512. Biedl, T., Brejova, B., Demaine, E.D., Hamel, A.M., Vinar, T., 2001. Optimal Arrangement of Leaves in the Tree Representing Hierarchical Clustering of Gene Expression Data. 14. Dept. Computer Sci., Univ., Waterloo. Brenner, S., Jacob, F., Meselson, M., 1961. An unstable intermediate carrying information from genes to ribosomes for protein synthesis. Nature 190 (4776), 576–581. Bussemaker, H.J., Li, H., Siggia, E.D., 2001. Regulatory element detection using correlation with expression. Nat. Genet. 27 (2), 167–174. Castells-Roca, L., Garcíamartínez, J., Moreno, J., Herrero, E., Bellí, G., Péréz-Ortí, J., 2011. Heat shock response in yeast involves changes in both transcription rates and mrna stabilities. PloS one e17272. Cheng, C., Fu, Y., Shen, L., Gerstein, M., 2013. Identification of yeast cell cycle regulated genes based on genomic features. BMC Syst. Biol. 7 (1), 70. Daub, C., Steuer, R., Selbig, J., Kloska, S., 2004. Estimating mutual information using b-spline functions-an improved similarity measure for analysing gene expression data. BMC Bioinf. 5 (1), 1. Dodge, Y., 2012. Statistical Data Analysis Based on the L1-Norm and Related Methods. Birkhäuser. Dwight, S., Harris, M., Dolinski, K., Ball, C., Binkley, G., Christie, K., Fisk, D., IsselTarver, L., Schroeder, M., Sherlock, G., Sethuraman, A., Weng, S., Botstein, D., Cherrya, J., 2002. Saccharomyces genome database (SGD) provides secondary gene annotation using the gene ontology (GO). Nucleic Acids Res. 30 (1), 69–72. Eisen, M., Spellman, P., Brown, P., Botstein, D., 1998. Cluster analysis and display of genome-wide expression patterns. Proc. Natl. Acad. Sci. 95 (25), 14863–14867. Engreitz, J., Morgan, A., Dudley, J., Chen, R., Thathoo, R., Altman, R., Butte, A., 2010. Content-based microarray search using differential expression profiles. BMC Bioinf. 11 (1), 603. Filkov, V., Skiena, S., Zhi, J., 2002. Analysis techniques for microarray time-series data. J. Comput. Biol. 9 (2), 317–330. Galdieri, L., Mehrotra, S., Yu, S., Vancura, A., 2010. Transcriptional regulation in yeast during diauxic shift and stationary phase. Omics: J. Integr. Biol. 14 (6), 629–638.

Please cite this article as: S. Ray, S. Misra, A supervised weighted similarity measure for gene expressions using biological knowledge, Gene (2016), http://dx.doi.org/10.1016/j.gene.2016.09.033

ARTICLE IN PRESS S. Ray, S. Misra / Gene xxx (2016) xxx–xxx Hardin, J., Mitani, A., Hicks, L., VanKoten, B., 2007. A robust measure of correlation between two genes on a microarray. BMC Bioinf. 8 (1), 1. Hecker, M., Lambeck, S., Toepfer, S., Va Someren, E., Guthke, R., 2009. Gene regulatory network inference: data integration in dynamic models - a review. Biosystems 96 (1), 86–103. Hellem, B., Dysvik, B., Jonassen, I., 2004. LSimpute: accurate estimation of missing values in microarray data with least squares methods. Nucleic Acids Res. 32 (3), e34. Clore, Hess, D., Myers, C., Huttenhower, C., Hibbs, M., Hayes, A., Paw, J., Mendoza, R., Luis, B., Nislow, C., Giaever, G., Costanzo, M., Troyanskaya, O., Caudy, A., 2009. Computationally driven, driven, quantitative experiments discover genes required for mitochondrial biogenesis. PLoS Genetics 5 (3). e1000407-e1000407. Kramer, M., Dutkowski, J., Yu, M., Bafna, V., Ideker, T., 2014. Inferring gene ontologies from pairwise similarity data. Bioinformatics 30 (12), i34–i42. Krause, E., 2012. Taxicab Geometry: An Adventure in Non-Euclidean Geometry. Courier Corporation. Lehninger, A., Nelson, D., Cox, M., 1993. Principles of Biochemistry. Worth, New York. Lin, D., 1998. An information-theoretic definition of similarity. ICML, vol. 98. Citeseer., pp. 296–304. Lopes, F., Ray, S., Hashimoto, R., Cesar, R., 2014. Entropic biological score: a cell cycle investigation for GRNs inference. Gene 541 (2), 129–137. Maesschalck, R., Jouan-Rimbaud, D., Massart, D., 2000. The Mahalanobis distance. Chemom. Intell. Lab. Syst. 50 (1), 1–18. Marcotte, E., Pellegrini, M., Thompson, M., Yeates, T., Eisenberg, D., 1999. A combined algorithm for genome-wide prediction of protein function. Nature 402 (6757), 83–86. Mewes, H., Frishman, D., Güldener, U., Mannhaupt, G., Mayer, K., Mokrejs, M., Morgenstern, B., Münsterkötter, M., Rudd, S., Weil, B., 2002. Mips: a database for genomes and protein sequences. Nucleic Acids Res. 30 (1), 31–34. Mistry, M., Pavlidis, P., 2008. Gene ontology term overlap as a measure of gene functional similarity. BMC Bioinf. 9, 1–11. Nagar, A., Al-Mubaid, H., 2008. A new path length measure based on go for gene similarity with evaluation using sgd pathways. In 21st IEEE International Symposium on Computer-Based Medical Systems, CBMS’08. pp. 590–595. Neiman, A., 2011. Sporulation in the budding yeast Saccharomyces cerevisiae. Genetics 189 (3), 737–765. Pfaff, D., 2002. Hormones, Brain and Behavior. Academic Press. Qian, J., Dolled-Filhart, M., Lin, J., Yu, H., Gerstein, M., 2001. Beyond synexpression relationships: local clustering of time-shifted and inverted gene expression profiles identifies new, biologically relevant interactions. J. Mol. Biol. 314 (5), 1053–1066. Ray, S., Bandyopadhyay, S., Pal, S., 2007. Dynamic range-based distance measure for microarray expressions and a fast gene-ordering algorithm. IEEE Trans. Syst. Man Cybern. B Cybern. 37 (3), 742–749. Ray, S., Maiti, S., 2015. Noncoding RNAs and their annotation using metagenomics algorithms. Wiley Interdisciplinary Reviews: Data Mining and Knowledge Discovery. 5(1). pp. 1–20. Ray, S., Pal, J., Pal, S., 2013. Computational approaches for identifying cancer miRNA expressions. Gene Expr. 15 (5-6), 243–253. Resnik, P., 1995. Using information content to evaluate semantic similarity in a taxonomy. arXiv preprint cmp-lg/9511007 Sawa, T., Ohno-Machado, L., 2003. A neural network-based similarity index for clustering DNA microarray data. Comput. Biol. Med. 33 (1), 1–15.

11

Sevilla, J., Segura, V., Podhorski, A., Guruceaga, E., Mato, J., Martinez-Cruz, L., Corrales, F., Rubio, A., 2005. Correlation between gene expression and go semantic similarity. IEEE/ACM Trans. Comput. Biol. Bioinf. 2 (4), 330–338. Sherlock, G., Hernandez-Boussard, T., Kasarskis, A., Binkley, G., Matese, J., Dwight, S., Kaloper, M., Weng, S., Jin, H., Ball, C., Eisen, M., Spellman, P., Brown, P., Botstein, D., Cherry, J., 2001. The Stanford microarray database. Nucleic Acids Res. 29 (1), 152–155. Shiratori, A., Shibata, T., Arisawa, M., Hanaoka, F., Murakami, Y., Eki, T., 1999. Systematic identification, classification, and characterization of the open reading frames which encode novel helicase-related proteins in Saccharomyces cerevisiae by gene disruption and northern analysis. Yeast 15 (3), 219–253. Slonim, D., 2002. From patterns to pathways: gene expression data analysis comes of age. Nat. Genet. 32, 502–508. Song, L., Langfelder, P., Horvath, S., 2012. Comparison of co-expression measures: mutual information, correlation, and model based indices. BMC Bioinf. 13 (1), 1–21. Spellman, P., Sherlock, G., Zhang, M., Iyer, V., Anders, K., Eisen, M., Brown, P., Botstein, D., Futcher, B., 1998. Comprehensive identification of cell cycle-regulated genes of the yeast Saccharomyces cerevisiae by microarray hybridization. Mol. Biol. Cell 9 (12), 3273–3297. Teng, L., Chan, L., 2006. Biclustering gene expression profiles by alternately sorting with weighted correlated coefficient. Inproceedings of the 2006 16th IEEE Signal Processing Society Workshop. Troyanskay, O.G., Dolinski, K., Owen, A.B., Altman, R.B., Botstein, D., 2003. A Bayesian framework for combining heterogeneous data sources for gene function prediction (in Saccharomyces cerevisiae). Proc. Natl. Acad. Sci. 100, 8348–8353. 14. Vijaranakul, U., Nadakavukaren, M., Bayles, D., Wilkinson, B., Jayaswal, R., 1997. Characterization of an NaCl-sensitive Staphylococcus aureus mutant and rescue of the NaCl-sensitive phenotype by glycine betaine but not by other compatible solutes. Appl. Environ. Microbiol. 63 (5), 1889–1897. Wang, H., Azuaje, F., Bodenreider, O., Dopazo, J., 2004. Gene expression correlation and gene ontology-based similarity: an assessment of quantitative relationships. In Computational Intelligence in Bioinformatics and Computational Biology, CIBCB’04. pp. 25–31. Wang, J., Du, Z., Payattakool, R., Philip, S., Chen, C.-F., 2007. A new method to measure the semantic similarity of go terms. Bioinformatics 23 (10), 1274–1281. Wilcox, R., 2012. Introduction to Robust Estimation and Hypothesis Testing. Academic Press. Yang, Y., Dudoit, S., Luu, P., Speed, T., 2001. Normalization for cDNA microarray data. BiOS 2001 The International Symposium on Biomedical Optics,. pp. 141–152. International Society for Optics and Photonics. Yu, Z., You, J., Li, L., Wong, H.-S., Han, G., 2012. Representative distance: a new similarity measure for class discovery from gene expression data. IEEE Trans. NanoBiosci. 11 (4), 341–351. Zhang, B., Horvath, S., 2005. A general framework for weighted gene co-expression network analysis. Stat. Appl. Genet. Mol. Biol. 4 (1). Zhou, Y., Young, J., Santrosyan, A., Chen, K., Yan, S., Winzeler, E., 2005. In silico gene function prediction using ontology-based pattern identification. Bioinformatics 21 (7), 1237–1245.

Please cite this article as: S. Ray, S. Misra, A supervised weighted similarity measure for gene expressions using biological knowledge, Gene (2016), http://dx.doi.org/10.1016/j.gene.2016.09.033