Gravitation field algorithm and its application in gene cluster
- Ming Zheng†1,
- Gui-xia Liu1Email author,
- Chun-guang Zhou1Email author,
- Yan-chun Liang1Email author and
- Yan Wang†1
© Zheng et al; licensee BioMed Central Ltd. 2010
Received: 19 July 2009
Accepted: 20 September 2010
Published: 20 September 2010
Searching optima is one of the most challenging tasks in clustering genes from available experimental data or given functions. SA, GA, PSO and other similar efficient global optimization methods are used by biotechnologists. All these algorithms are based on the imitation of natural phenomena.
This paper proposes a novel searching optimization algorithm called Gravitation Field Algorithm (GFA) which is derived from the famous astronomy theory Solar Nebular Disk Model (SNDM) of planetary formation. GFA simulates the Gravitation field and outperforms GA and SA in some multimodal functions optimization problem. And GFA also can be used in the forms of unimodal functions. GFA clusters the dataset well from the Gene Expression Omnibus.
The mathematical proof demonstrates that GFA could be convergent in the global optimum by probability 1 in three conditions for one independent variable mass functions. In addition to these results, the fundamental optimization concept in this paper is used to analyze how SA and GA affect the global search and the inherent defects in SA and GA. Some results and source code (in Matlab) are publicly available at http://ccst.jlu.edu.cn/CSBG/GFA.
Two of the most challenging tasks of optimization algorithms are to search the global optimum and to find all local optima of the space of solutions in clustering genes from available experimental data , e.g. the gene expression profiles, or given functions. In view of recent technological developments for large-scale measurements of DNA expression level, these two problems can often be formulated and many methods have been proposed. In particular, the heuristic searches are more promising than other kinds of searching approaches. These approaches include GA (genetic algorithm) , SA (simulated annealing) , PSO (Particle Swarm Optimization)  etc. But some inherent drawbacks, especially the inability to the multi-modal functions optimization, can be found from the traditional heuristic search algorithms above. Each of these concepts allows for several modifications, which multiplies the number of possible models for data analysis we can change the algorithm themselves, to find all the valleys of given functions. But we still have a lot of parameters to consider, as known as the number of valleys and the valley distance etc, and the performances of the modifications are not good enough.
GA is a traditional searching-optimization algorithm, this algorithm can search global optimal solution with probability criterion , but it can't converge to the global optimal solution in the theory . So GA always traps in local optima or genetic drift. Anyway, the run-time of this algorithm is acceptable for most cases of searching-optimization problems.
SA is a generic probabilistic meta-algorithm for global optimization problems, namely locating a good approximation to the global optimum of a given function in a large search space. If we search enough time with SA, it will converge in the global optimal solution with probability 1 . But the biggest drawback of SA is that the running-time of SA is so long that the efficiency is not tolerant to us.
Recently, some other searching algorithms have been proposed and discussed, such as PSO etc. These algorithms can search solution like GA. But actually, most of them will not converge to the global optimal solution either.
In this paper, we propose a new heuristic search approach Gravitation Field Algorithm (GFA). It not only can handle unimodal functions optimization, which traditional heuristic searches can do the same, but also can deal with multimodal functions optimization, which traditional heuristic search algorithms can't do. The experiments of the benchmark functions demonstrate that. The idea of GFA is derived from the modern widely accepted theory of planetary formation-- the Solar Nebular Disk Model (SNDM)  in the astronomy.
The complex astronomy theory can be introduced simply as follows:
Several billions years ago, there wasn't any planet in the Solar System; just dusts rounded the whole world. Then the dusts assembled by their own gravitations. For a long time, the rocks had come out. It was the threshold time of the whole Solar System. From then, the rocks moved fast to assemble together, and the bigger rocks arrested the smaller rocks. And they became larger rocks. Finally, the planets came out, and the rocks around them were absorbed out.
GFA is derived from the point of view of the hypothesis theory described above. To start with, all the solutions, which are the dusts in the algorithm model, are initialized randomly, or based on the prior knowledge; what's more, we assign every dust (solution) a weight, we call it mass, whose values are based on the mass function generated from the space of the problem solution; finally, the GFA begins. The power of attraction, which belongs to a certain dust and exists between every two dusts, pulls other dusts, which have the same influence to other dusts. Hence, the dusts assemble together, and the planets come out in the end--they are the optima. If you want to find global optimal solution, the planets assemble again, and the biggest planets will come out. To give a penalty of that the highest mass dust rules the whole space of solution, we propose a distance which can reduce the effect of gravitation field.
Description of GFA
Before all the works start, we design a mass function on the basis of special problem. The mass function is similar to the fitness function in GA. Both of them are score functions which are the criterions for a special solution. We can make the mass function be proportional to, or inversely proportional to, the extreme values of the problem. It all depends.
The algorithm begins with the initial step. We generate and select n dusts di (i = 1, ..., n) in the mass function domain [xbegin, xend] to build the initial solution set. The positions of the dusts can be random or based on the prior knowledge (i.e. there will be greater probability to exist peak values in some positions). Then we assign a mass value to every dust. The mass values, which are described above, are associated with their positions, and are calculated by the mass function. So when the dusts move as we described below, they have a certain probability to find a new position, in which there is a new mass value that is bigger than the mass value of the centre dust. This initialization approach was used because the extension of the space solution could be considered. And the density of the dusts is proportional to the probability of existence of extreme value.
Strategy of division
Where x is the domain of x axis of the mass function; y is the domain of y axis of the mass function.
The rules of the motion of the dusts and strategy of absorption
The rules of the motion of dusts are the kernel part of GFA. This part of GFA decides which one is the peak value among all the dusts in one group. And this step of the algorithm also can generate new mass values which can have bigger mass values than the mass values of centre dust with certain probability.
Each group contains only one dust, called centre dust, that doesn't move within an epoch (but may move in other epochs). The rest dusts in this group move toward the centre dust within the epoch. Formally, centre (di), where i = 1,..., n, is used to represent whether the dust di is the centre dust or not: centre (di) is set to true if di is the centre dust; centre (di) is set to false otherwise.
Eq. (3) is easy but efficient. disi is the distance between the moving surrounding dust and the centre dust. M is a weight value for distance. 0.618/10 is chosen as M for this paper. Actually, many other weight values were tested in our experiment. But the speed and efficiency of GFA with the Golden Ration (0.618) are best. Maybe the value of it is not coincidence. But the reason, which is beyond our comprehension, will not mention in the paper. When the space of a group is small enough, the surrounding dusts will move with small pace calculated using Eq. (3); when the space of the group is big, surrounding dusts will move with big pace. Big pace will appear in two ways: One appears in the late algorithm. In this period, every dust in each group is quasi-optimal, big pace will not affect the efficiency of the GFA. So we make a big pace for fast convergence. The other one appears that we set a small number of dusts in the initial step. In this condition, maybe we want to make a fast convergence. And big pace could accelerate the convergence.
When all the distances between the surrounding dusts and the centre dust are small enough in one small group, such as smaller than a threshold, the surrounding dusts will be absorbed by the centre dust. Based on these rules of motion, we design the strategy of absorption as follows:
All the surrounding dusts are deleted, but the centre dust will not be. We set centre (di) = false, it represents the small group for the next step. When the number of surrounding dusts is bigger than a threshold, we will absorb all the other surrounding dusts for saving running time.
After the absorption of the small groups is complete, the next epoch begins if the algorithm has not converged to the optimal solution. We will divide the space of the solution again, and compute for searching peak values using the survive dusts.
The most important advantage of GFA is the ability to deal with the multimodal objective functions (i.e. the mass functions). GFA can converge for one independent variable mass functions with probability 1. We give a theorem and its strict mathematical proof as follows.
Where m = 1, 2, ..., n. After moving toward xmax in the group, the dusts become x1', x2', ..., xmax', ..., xn'.
First of all, we will give a description of groups and define two conceptions in the graph theory as follows:
Definition 1. In Two Side (ITS): In the line segment of the x m and x max in the variable space, iff ∃ x m ', such that: f (x m ')≥ f (x max '), then we call x m and x max are ITS.
Definition 2. In One Side (IOS): In the line segment of the x m and x max in the variable space, iff ∃ x m ', so that: f (x m ')≥ f (x max '), then we call x m and x max are IOS.
Now we give a theorem and prove it:
The scale of group is small enough; such that the number of peaks is at most one.
The motion of the surrounding dusts is smoothing.
The number of dusts in one group is big enough.
When xmax = xpeak, the group is convergent obviously, f(xmax) = f(xpeak).
- 2)When the xmax ≠ xpeak, then:
If exists xm, max and xm are ITS, the group will be convergent obviously according to condition 2), as Fig. 4 (a). The peak value calculated in the group by the GFA is the real peak value. We call this kind of condition is Real-Peak Condition (RPC).
If doesn't exist xm, such that the xmax and xm are ITS, the group will not be convergent, the xmax will be the pseudo peak value of the group, as Fig. 4 (b) and Fig. 4 (c). Pseudo peak value exists with probability p. we call this kind of condition is Pseudo-Peak Condition. (PPC)
The probability p can be calculated as follow:
We define the power distance of the center dust (see Description of GFA section) in the space of the solution in the group for one independent variable is ltotal. ltotal is shown as the domain of x axis in Fig. 4 (a), (b) or 4(c). And the distance between the xmax and the nearest valley in the group for one independent variable is lmax. lmax is shown as the distance between xmax and the xv in Fig. 4 (a), (b) or 4(c). f(xv) is the valley which is nearest xmax in all valleys in the group.
So when PPC is on, one group will converge in the real peak, like RPC is on. So every group will converge with probability 1.
When every group converges, the surrounding dusts are absorbed out. The real peak and pseudo peak will be the new dusts in the following step.
When Xmax becomes a dust in a new group of a new epoch, the global peak value could be Xmax in one group, obviously.
Hence in a certain domain of variable x, GFA will converge in the global maximum.
Peaks versus valleys
If you want to minimize f(x), like the experiments below, you can do so by maximizing -f(x), because the point at which the maximum of -f(x) occurs is the same as the point at which the minimum of f(x) occurs .
Computation complexity of GFA
Computing the mass function value once needs processing time Tf. Computing motion of every surrounding dust is Tm. Computing the size relationship of every dust once needs processing time Tg. We set the number of dusts in one group is n, the number of groups is m, and the number of times of moving a dust is s.
The computation complexity of Eq. (11) is o(mns).
The computation complexity of Eq. (14) is o(mns).
If the complexity of mass function and the dimensionality of f(x) are big enough, we can say T = Tf. The number of dusts is m × n, so the computation complexity of GFA is the product of the number of dusts and the number of times of moving of surrounding dusts. That is, the computation complexity is o(mns).
Results and Discussion
To test the efficiency of GFA, we assessed the performance of the GFA by employing a suite of different benchmark mathematical functions (Eqs. (16-20)) and by comparing the performance with GA and SA. For each of the five test functions and the each method, 500 minimization runs were performed and mean squared error, standard deviation and mean gauss error values were calculated.
In the following benchmark functions (Eqs. (16-20)), D donates the number of independent variables, and we defined D = 50. The benchmark functions were selected as following:
Where n is the number of runs, f(xi) is the performance for run i and the fopt(xi) is the function value at the global minimum.
Searching global minimum
Configuration of parameter setting for GFA, GA and SA.
Max. numbers of iterations
Number of polulations
MSE, STD, and MGE of GA, GFA, and SA, Best performance (i.e., lowest error) for each function is highlighted in bold underline letters.
From Table 2, we found that sphere's, Rastrigin's and Griewank's functions could be optimized better by GFA than GA and SA, especially for the sphere's and Rastrigin's function.
Mean numbers of epochs until the minimization threshold was reached and mean number of failures.
mean number of epochs
number of failures
mean number of epochs
number of failures
mean number of epochs
number of failures
GFA was able to outperform GA and SA in four of five benchmark functions in terms of epochs needed and least failures. Only in one of the five functions, GFA did not outperform in both speed and robustness. The minimization of the Ackley's function took slightly more epochs with the GFA (97) than with the GA (92). From the Table 3, we know that GA could do better for optimizing the complex function, such as Ackley's function. But for some simple function, such as sphere's function, GFA could do better.
For more accurate computation, we defined the initial number of dusts 10,000, the max power distance of center dust 5, and the epochs 1,000, which was the same as used in GA and SA (for more information of detail parameters, please see description of GFA).
total running-time of GFA, GA and SA with 500 runs.
Although minimum of certain given function could be optimized to resolve the problem for most situations, minimum is not always what we want only. Sometimes a part of minima is also needed. Bayesian network inferring, for example, is a stochastic approach. So the minimum of Bayesian criterion function is not always fit the realistic network best. It is important to find other minima for the problems as cases described above. A number of lowest minima can be retained in the search.
To run multimodal optimization algorithm, we must track more information than the global optimization. Searching exactly a certain number of minima or even all of the minima in the domain is a challenging work, which is also meaningful for some problems, e.g. gene cluster. In the experiment, top 5 minima were searched with GFA. It's beyond GA's and SA's ability.
GFA settings configurations of multi-minima optimization
Domain of each variable
Max numbers of iterations
number of dusts
Number of groups
Top 5 mean minima value of 500 runs for each benchmark function
Application to gene cluster
GFA could be applied in many science research areas, especially bioinformatics. We used it for gene cluster. The clustering algorithm we used was K-means. The dataset used in this paper was from the experiment of Spellman . And the data which was MIAME compliant  could be downloaded from GEO with accession number GDS38. The number of cases was 7,680, of which the 17 missing values were excluded. We computed SMBS correlation coefficient  and excluded missing values with Matlab. We excluded cases when the respondent was dropped only on analyses involving variables that had missing values. The data with 7663 cases and 16 samples was divided into 20 parts, as known as clusters. So we distributed the cases into 20 groups. And our mission was that make sure the mean distance was smallest or the number of runs exceeded 1,000. To compare with GFA, 3 methods, GA, SA and Cluster 3.0 , were used to test the efficiency of the novel algorithm. K-means was used in Cluster 3.0 with correlation (centered) coefficient. For there four clustering methods, 500 runs were performed. And we got the mean value.
To visualize the result of the gene cluster, we used the free software TreeView  from Eisen's Lab. A part results computed by the novel algorithm and other cluster methods of gene cluster were showed in a picture which is Additional file 1. This picture shows the C8.txt and result.txt with Group = 7, which should be opened with excel software.
But there is no single best criterion for obtaining a partition because no precise and workable definition of cluster exists. Clusters can be of any arbitrary shapes and sizes in a multidimensional pattern space . It is impossible to objectively evaluate how good a specific clustering method is without referring to what the clustering will be used for . So we evaluate the result with F test with only mathematical aspect. And P values of all the 16 samples are less than 0.01. The total results of the gene cluster with GFA could be downloaded from our website.
cis-regulatory elements correspond each method
From Table 7, the efficiency of GFA could be seen. The genes with the same cis-regulatory elements could be clustered better by GFA than Cluster 3.0 and SA. Only in the case of the Ace2 cis element, the GFA with optimized parameters could not outperform the Cluster 3.0. Some efficiency is very obvious, especially Ndd1 and Swi5. But it seams that the efficiency of GA and GFA are the same. The results of this cluster experiment indicate that GFA method does work in gene cluster with finding cis-regulatory element in the same cluster. Other clustering results have the same properties.
The developed algorithm software of single-machine for GFA was implemented with matlab R2008a in the software package GFA, which is a short code file, freely available from our website: http://ccst.jlu.edu.cn/CSBG/GFA. You can implement the multi-machines parallel algorithm for GFA based on the detail in the conclusion section below.
In this paper we propose a generalization searching-optimization algorithm called GFA. This algorithm derives from the SNDM theory, and the efficiency of the algorithm is shown as above. We can summarize them into three parts:
In the form of one independent variable, the GFA will converge with probability 1. It gives us a balance level between time and efficiency. If you want to find exactly extreme value of the mass function, you can disperse more dusts to avoid the condition of the p (see mathematical proof). If you want to find the extreme value fast, you can define a small number of initial dusts.
GFA can find all needed the peaks of the solution. No matter the number of initial dusts is big or small.
The running-time of GFA is very short. The reasons are the strategy of the division and the rules of motion described as above. The space of solution is cut into small groups. It's the decomposition of any complex problem. Even more, it will support us a feasibility of parallel computing. In this view, this algorithm is similar to the parallel genetic algorithm . But this mechanism of GFA can be faster than GA's. We could use a large number of idle computers to calculate a complex problem, and the running-time is inversely proportional to the number of the computers.
This work is supported by the NSFC (60873146, 60973092, 60903097); the National HighTech R&D Program of China (863) (2009AA02Z307); the Ph.D. Programs Foundation of Ministry of Education of China (20090061120094); the project of 200810026 support by Jilin University, "211" and "985" project of Jilin University.
- Fang-xiang Wu: Genetic weighted k-means algorithm for clustering large-scale gene expression data. BMC Bioinformatics. 2008, 9 (Suppl 6): S12-10.1186/1471-2105-9-S6-S12. 10.1186/1471-2105-9-S6-S12View ArticleGoogle Scholar
- James T, Shuba G: Correction: genetic algorithm learning as a robust approach to RNA editing site site prediction. BMC Bioinformatics. 2006, 7: 406- 10.1186/1471-2105-7-406View ArticleGoogle Scholar
- Rui J, Hua Y, Fengzhu S: Searching for interpretable rules for disease mutations: a simulated annealing bump hunting strategy. BMC Bioinformatics. 2006, 7: 417- 10.1186/1471-2105-7-417View ArticleGoogle Scholar
- Michael M, Michael S, Gisbert S: Optimized Particle Swarm Optimization (OPSO) and its application to artificial neural network training. BMC Bioinformatics. 2006, 7: 125- 10.1186/1471-2105-7-125View ArticleGoogle Scholar
- Pier-Luigi L, Santo M, Francesco P: Discovery of cancer vaccination protocols with a genetic algorithm driving an agent based simulator. BMC Bioinformatics. 2006, 7: 352- 10.1186/1471-2105-7-352View ArticleGoogle Scholar
- Rudolph G: Convergence properties of canonical genetic algorithms. IEEE Trans Neural Networks. 1994, 5 (1): 96-101. 10.1109/72.265964PubMedView ArticleGoogle Scholar
- Geman S, Gemana D: Stochastic relaxation, Gibbs distributions and the Bayesian restoration of images. IEEE Trans. on Pattern Analysis and Machine Intelligence, PAMI. 1984, 721-741. 10.1109/TPAMI.1984.4767596.Google Scholar
- Safronov VS: Evolution of the Protoplanetary Cloud and Formation of the Earth and the Planets. 1972, Jerusalem: Israel Program for Scientific Translations,Google Scholar
- Genetic Algorithm and Direct Search Toolbox User's Guide, copyright by the MathWorks, Inc. http://www.mathworks.co.uk/access/helpdesk/help/pdf_doc/gads/gads_tb.pdf
- Lohr SL, Rao JNK: Jackknife estimation of mean squared error of small area predictors in nonlinear mixed models. Biometrika. 2009, 2: 457-10.1093/biomet/asp003. 10.1093/biomet/asp003View ArticleGoogle Scholar
- Petros EM, Philippe C: An EWMA chart for monitoring the process standard deviation when parameters are estimated. Computational statistics & data analysis. 2009, 5: 2653-2664.Google Scholar
- Cody WJ: Rational Chebyshev Approximations for the Error Function. Mathematics of Computation. 1969, 23: 631-637. 10.1090/S0025-5718-1969-0247736-4View ArticleGoogle Scholar
- Wan X, Karniadakis GE: A sharp error estimate for the fast Gauss transform. Journal of computational physics. 2006, 11: 7-12. 10.1016/j.jcp.2006.04.016View ArticleGoogle Scholar
- The website of Matlab tools for GA and SA. http://www.mathworks.com/products/gads/
- Paul TS, Gavin S, Michael QZ, : Comprehensive Identification of Cell Cycle-regulated Genes of the Yeast Saccharomyces cerevisiae by Microarray Hybridization. Molecular Biology of the Cell. 1998, 9: 3273-3297.View ArticleGoogle Scholar
- Tim FR, Philippe RS, Paul TS, : A simple spreadsheet-based, MIAME-supportive format for microarray data: MAGE-TAB. BMC Bioinformatics. 2006, 7: 489- 10.1186/1471-2105-7-489View ArticleGoogle Scholar
- Wang Han, Liu Gui-xia, Zhou Chun-guang, : Measuring the Similarity of Co-regulated Genes by Integrating Quantity and Tendency of Gene Expression Changing. The 2nd International Conference on Bioinformatics and Biomedical Engineering. Shanghai, China: Bioinformatics and Biomedical Engineering. 2008, 1896-1900.Google Scholar
- The website for Cluster 3.0. http://bonsai.ims.u-tokyo.ac.jp/~mdehoon/software/cluster
- The website for TreeView. http://rana.lbl.gov/downloads/TreeView/TreeView_vers_1_60.exe
- Jain AK, Dubes RC: Algorithms for Clustering Data. 1988, Prentice Hall, Englewood Cliffs, NJ,Google Scholar
- D'haeseleer P, Liang S, Somogyi R: Genetic network inference: from co-expression clustering to reverse engineering. Bioinformatics. 2000, 8: 707-726. 10.1093/bioinformatics/16.8.707View ArticleGoogle Scholar
- The website for Li-Hsieh Lin's cis elements. http://www.biomedcentral.com/content/supplementary/1471-2105-6-258-S1.xls
- Lin Li-Hsieh, Lee Hsiao-Ching, Li Wen-Hsiung, : Dynamic modeling of cis-regulatory circuits and gene expression prediction via cross-gene identification. BMC Bioinformatics. 2005, 6: 258- 10.1186/1471-2105-6-258PubMedPubMed CentralView ArticleGoogle Scholar
- Petty CC, League MR: A theoretical investigation of a parallel genetic algorithm. Proc. Of 3rd Int. Conf. On Genetic Algorithm, Morgan Kaufmanm. 1989, 398-405.Google Scholar
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.