A Novel Algorithm for the Precise Calculation of the Maximal Information Coefficient

Measuring associations is an important scientific task. A novel measurement method maximal information coefficient (MIC) was proposed to identify a broad class of associations. As foreseen by its authors, MIC implementation algorithm ApproxMaxMI is not always convergent to real MIC values. An algorithm called SG (Simulated annealing and Genetic) was developed to facilitate the optimal calculation of MIC, and the convergence of SG was proved based on Markov theory. When run on fruit fly data set including 1,000,000 pairs of gene expression profiles, the mean squared difference between SG and the exhaustive algorithm is 0.00075499, compared with 0.1834 in the case of ApproxMaxMI. The software SGMIC and its manual are freely available at http://lxy.depart.hebust.edu.cn/SGMIC/SGMIC.htm.

A ll kinds of relationships determine the development of things [1][2][3] . Relationships and associations should therefore be identified and measured to explore the rules of development. A typical example is measuring the relationships between genes by determining the associations between their expression profiles 4,5 . Many methods have been developed to measure associations through calculation of correlation coefficients, such as Pearson's, Spearman's, mutual information 6,7 , CorGC 8 , and maximal correlation 9 . Recently, Reshef et al. 10 proposed a novel correlation measurement ''maximal information coefficient'' (MIC), and gave a 1-D dynamic programming algorithm, ApproxMaxMI, to calculate MIC. MIC does not rely on the distributional assumptions of measured data and could identify a broad class of associations compared with previous studies. The MIC of two vectors x and y is defined as follows.
MIC~maxfI(x,y)=log 2 minfn x ,n y gg, where I(x,y)~H(x)zH(y){H(x,y) p(x i y j )log 2 1 p(x i y j ) ð1Þ n x ?n y ,B(n), B(n) 5 n 0.6 . In calculating MIC for gene expression profile vectors x and y, n is the number of data points of gene expression profiles, and n x , n y is the number of bins of the partition of the x-and y-axis 10 , respectively. After MIC and its algorithm were published, many applications and discussions appeared [11][12][13][14][15] . Reshef et al. also foresaw the possible disadvantage of the algorithm ''ApproxMaxMI'', and suggested it should be replaced in the future if a method efficiently finds solutions that are closer to optimal or even optimal is developed 10 . Here, we conducted an initial attempt in this direction.

Results
ApproxMaxMI does not optimize the partition of the y-axis. ApproxMaxMI first fixes an equipartition of n data points with horizontal lines, and then calculates MIC values by moving vertical lines to optimize the x-axis partition. However, the partition of the y-axis should be optimized simultaneously instead of being fixed as an equipartition. MIC value of a pair of gene expression profiles (YAL001C:YAL020C) of yeast 16 was 0.30732 according to ApproxMaxMI, but the MIC value optimized by exhaustive algorithm was 0.51582, this value can be obtained directly from the partition scheme in Fig. 1A. Similarly, the MIC value of YAL001C:YAL039C was 0.28519 according to ApproxMaxMI, but the MIC value optimized by exhaustive algorithm was 0.42340 which can be obtained directly from the partition scheme in Fig. 1B. Such significant differences can also be identified from the gene expression profiles of fruit fly (Fig. S1-S4; Droso174figure.zip at http://lxy.depart.hebust.edu.cn/ SGMIC/SGMIC.htm). From the observations, as foreseen by Reshef et al., the MIC values calculated by ApproxMaxMI were not always convergent to optimum MIC values for the following simple proposition: Proposition 1.1: Fixing an equipartition of n data points by horizontal lines is neither a sufficient nor a necessary condition to obtain MIC 5 max{I(x, y)/log 2 min{n x , n y }}. Proposition 1.1 is proven in Supplementary Section 2.1. In this sense, ApproxMaxMI does not achieve the equivalent transformation from a 2-D search to a 1-D search. Theorem 1 and Proposition 6.12 of Reshef et al 10 are related to algorithm convergence of ApproxMaxMI, however they only discuss the convergence to 0 of the MIC value of independent X and Y as the number of data points nR' rather than the convergence to the global optimum MIC value. The convergence of ApproxMaxMI cannot be proven by the approaches because of the following proposition: Proposition 1.  It's worth noting that this proposition guarantees that SG can reach the state corresponding to the optimum MIC values in finite steps from any given initial state. In fact, SG can converge to the global optimum of MIC values shown as follows. Proposition 1.5: SG is convergent. Proposition 1.6: SG is equivalent to the exhaustive algorithm with a sufficient number of iterations.
The equivalence of SG and the exhaustive algorithm is shown in Table S1 of Supplementary Section 3.1.
To set a gold standard for these algorithms, we provide an exhaustive algorithm (Algorithm 3 of Supplementary Section 2.2) to calculate MIC values. The flow chart of SG (Fig. S5) Fig. 1A and 1B). We also compared SG with ApproxMaxMI on a larger data set of three species. From the 4,498,500 relationships mentioned above in the 3,000 genes of yeast 16  For random clouds at sample size n, the larger the value of n is, the more individuals (or chromosomes) and longer running time are needed to show the advantage of SG method over ApproxMaxMI. Specifically, for n , 30 points, 20 individuals are usually enough; for 30 , n , 50, 100 individuals are needed; for 50 , n , 100, we need 1000 individuals; while for 100 , n , 200, at least 10000 individuals are required. These data and running time are available in file populationchangewithN.xlsx (http://lxy.depart.hebust.edu.cn/SGMIC/ SGMIC.htm). For 200 , n , 500, perhaps over 100000 individuals are needed. In this sense, for gene expression profile data, which usually have less than 50 time points, 100 individuals are enough. Because computing MIC value by SG for 500 random points need over 100000 individuals and over 12 hours, we can only calculate a SGMIC value for an example in points500example.xlsx (http://lxy. depart.hebust.edu.cn/SGMIC/SGMIC.htm), and the MIC value by SG with 100000 individuals is 0.150593, while the MIC by ApproxMaxMI is 0.14799.
For non-random relationships, the efficiency of SG is good, as shown in next section, where the majority of MIC values by SG were near 1 when the relationships had no noise, although only 20 individuals are used in SG.

Discussion
Properties of MIC Values by SG. Because SG is a precise method in calculating MIC values, we redescribed the main properties of MIC values for four classes of main relationships, namely, function relationship, non-function relationship, function with noise, and non-function with noise. Fourteen relationships without noise were drawn by us in Fig. 3, and their formulas are presented in Table S4 of Supplementary Section 4. Function relations are represented by trigonometric functions (e.g., sin and tan), power functions (e.g., y 5 x 2 ), exponential functions (e.g., y 5 10 x ), inverse proportion functions (e.g., y~1 x ), and composite functions.
Meanwhile, some non-function relationships, such as the taijitu (symbol for the yin-yang principle, which originated from Yi Jing of ancient China), galaxy figure, heart-shaped line, and polygonal line, were also selected. To investigate the effect of noise on relationmeasure algorithms, three relationships with increasing noise for each of the 14 relationships were drawn by us in Fig. 3.  However, when some relationships (e.g., Fig. 3 A-E) were involved, SG and maximal correlation (ACE) algorithm exhibited a much better performance than that of ApproxMaxMI. Second, the MIC values by SG strictly decreased with increasing noise, a finding indicating that SG has a strong ability to distinguish noise from real signals. However, a strict decreasing trend was not observed with other algorithms. For example, the maximal correlation (ACE) of d1, d2, d3, and d4 was 0.33508, 0.64533, 0.55032, and 0.22122, respectively.
The Application in Predicting Yeast Protein Interaction. To describe the performance of SG more intuitively, we use the Data Repository Yeast Genetic Interactions (DRYGIN for short) 23 , which are derived from large-scale Synthetic Genetic Array (SGA) genetic analysis and the Genetic interaction score (e) can represent the genetic interaction quite accurately. The file sgadata_costanzo2009_ intermediateCutoff_101120.txt.gz contains 76406 genetic interactions with an intermediate cutoff applied (jej . 0.08, p-value , 0.05). The MIC value of 96.46% (73702 out of 76406) genetic interaction by SG is larger than that by ApproxMaxMI, though the maximum MIC values by both methods are the same, i.e. 0.9986. We found 125 interactions (Fig. 4) with high SG score but low ApproxMaxMI, these interactions are not included in the DRYGIN but appear in CCSB interactome 24 . For example, the ApproxMaxMI values of interactions YDR382W and YDL082W, YDR447C and YDL083C, YGL076C and YBR048W, YGR034W and YGL076C are 0.54236, 0.56215, 0.52781 and 0.57648 respectively, while the SG values are all 0.932112. The four interactions are all verified by CCSB interactome but are not included in DRYGIN, which shows the advantage of SG over ApproxMaxMI.

Methods
Datasets. The gene expression profile datasets used are from transcriptome of yeast 16 19,20 . In each thread, the number of chromosomes is 20 by default. Each chromosome consists of genes, which are the abscissa (xgene) of vertical partition lines and the ordinate (ygene) of horizontal partition lines, and these genes form an x-by-y grid. We crossed over an xgene only with another xgene and a ygene only with another ygene during the genetic crossover step. After the mutation, crossover, and simulation annealing, we calculated the fitness of each chromosome and kept the optimum solution. In the annealing, the reproduction operator was derived from Metropolis criteria. The fitness value, championed for over 30 generations, was considered as the MIC value from SG.
Exhaustive Algorithm. Inputting the expression profiles of a pair of genes, we transformed them into a consecutive integer series, kept the order relation of the original x and y coordinates, and denoted the maximum number of the series as n. The computation loop for the x and y coordinates was 1 to 2 n-1 -1. For both x and y coordinates, we transformed the number of the current loop into binary number. Using the binary number, we determined whether to insert a partition line between two space-adjacent expression points (1 indicates a partition line and 0 means no partition line). MIC values were saved for each loop out of (2 n-1 -1) 2 loops, and the maximum value was considered as the MIC value by the exhaustive algorithm.
Mutation Procedure. The mutation with self-adapting parameters is used, and the mutation frequency is defined as: pm~p where f max is the largest fitness of current population, f ave is the average fitness of current population, f is the individual fitness. Parameters pm1 . pm2 with 0 , pm1,pm2 , 1 are mutation parameters. First, an individual (or chromosome) is assigned a random number in (0,1): if the random number is no less than pm, then the individual will mutate; skip the step otherwise. Here, an individual (or chromosome) is a set of vertical or horizontal lines constituting the bins of the partition of the x-and y-axis 10 . Mutating an individual is changing the positions of some vertical or horizontal lines along the x-or y-axis. The first and last horizontal (or vertical) lines are only allowed to move between its original position and the next line, and the other lines are allowed to move between its left and right neighbor lines freely.
P-value of SG algorithm. We use a non-parameter test method: Wilcoxon Rank-sum Test to calculate the p value, the Rank-sum Test is capable of determining if two gene profiles are independent. The p-value can be seen from file Pvalue.xlsx in http://lxy. depart.hebust.edu.cn/SGMIC/SGMIC.htm.
The optimization of B(n). The setting n x ?n y , B(n), B(n)5n 0.6 is mainly based on the following consideration. Based on the Theorems 1 and 2 in Section 6.2 of Reshef et al. 10 , we know that setting B(n) , n can avoid inflated MIC scores with large B(n). Then we search for an optimal value of B(n), above which statistically independent data receive scores bounded away from zero as sample size grows, and below which they receive scores approaching zero. The Supplementary Section 7 Fig. S7 shows the alpha 5 0.6 meets the criterion.