BACKGROUND
based on the assumption that co-expression indicates co-regulation, gene expression data clustering aims to reveal gene groups of similar functions in the biological pathways.
this biological rationale is readily supported by both empirical observations and systematic analysis  <cit> .
in particular, consider gene expression time-course experiments, where the data are made up of tens of thousands of genes, each with measurements taken at either uniformly or unevenly distributed time points often with several replicates.
clustering algorithms provide a good initial investigation into such large-scale datasets, which ultimately leads to biological inference.
an excellent review of current techniques and all subsequent analysis can be found in  <cit> .
various model-based methods have been proposed to accommodate the needs for data mining in such massive datasets.
among them are mixed effects models  <cit>  and auto regressive models  <cit> .
the basic approach of these model-based methods is to fit a finite mixture model to the observed data, assuming that there is an underlying true model/density, and then systemically find the optimal parameters so that the fitted model/density is as close to the true model/density as possible.
it is observed that model-based approaches generally achieve superior performance to many others  <cit> .
however, current methods can be problematic, as they often fail to show how clustering can assist in mining gene expression data.
the maximum likelihood estimator  is one of the most extensively used statistical estimation techniques in the literature.
for a variety of models, likelihood functions  <cit> , especially maximum likelihood, have been used for making inferences about parameters of the underlying probability distribution for a given dataset.
the solution often involves a nonlinear optimization such as quasi-newton methods or, more commonly, expectation-maximization  methods  <cit> .
the problem with the former method is that the quantities are estimated only when they satisfy some constraints, while with the latter method all parameters have to be explicitly specified, so the number of clusters k has to be known a priori, which is not practical in microarray data analysis.
there are many unique features of mle, including its efficiency.
however the practical deficiencies of mle, besides those with its optimization, are the lack of robustness against outliers and its sensitivity to the correctness of model specification.
we discuss in this paper the performance of an appealing alternative, the minimum distance estimator   <cit> , which is less explored in this field.
inspired by the work of  <cit> , we propose to incorporate mde in our algorithm for gene expression time-course analysis.
mde provides robust estimation against noise and outliers, which is of particular importance in gene expression data analysis, where data are often noisy and there are few replicates.
tight clustering has been proposed as a response to the needs for obtaining smaller clusters in genomic signal processing.
it was motivated by the fact that the most informative clusters are very often the tight clusters, usually of size 20– <dig> genes  <cit> .
tight clustering refers to methods that can be built upon an existing partition to obtain core patterns that are more interpretable.
the initial partition can be obtained empirically or by using generic algorithms such as k-means.
as a result, more information can possibly be revealed.
for example, if genes in the same functional category are allocated into different tight clusters, one may pursue the underlying explanation by looking into these clusters.
one possible result of such investigation, for example, is new function discovery.
in this sense, to obtain tight clusters, some genes should be classified as scattered genes, if forcing them into clusters will only disturb biologically relevant patterns.
indeed, the issue of scattered genes has received more attention recently  <cit> .
however, in contrast to the prevailing concept that scattered genes should be treated as outliers and discarded from further study, we prove that some scattered genes can be of biological significance.
current methods for gene expression time-course data rarely deal with scattered genes.
to the best of our knowledge,  <cit>  is the first to address this issue, but it results in heavy computation due to the nature of random resampling.
it was proposed in  <cit>  that outliers can be modelled by adding a poisson process component in the mixture model.
however, this method has not been verified in this field, and it relies on correct model specification.
there has been a lot of research focusing on modelling time-course data by splines and autoregressive models, usually followed by em  <cit> .
in  <cit> , the cubic b-spline, which is a linear combination of b-spline basis functions, is used for fitting gene expression time-course data.
to avoid over-fitting, it is suggested not to fit a curve to every individual gene, but to constrain the spline coefficients of co-expressed genes to have the same covariance matrix.
alternatively, we propose in this work a novel approach to fit our spline model.
splinecluster  <cit>  is an efficient hierarchical clustering program based on a regression model with a marginal likelihood criterion.
starting from singleton clusters, the idea is to merge clusters based on marginal likelihood in each iteration.
it is efficient and straightforward to visualize.
the problem is that it overlooks microarray replication information by using only the mean of all replicates, which leads to loss of information.
as microarry experiments are increasingly performed with replicates, the additional information provided by replicated measurements is a valuable source of variability in terms of effective clustering  <cit> .
the outline of this paper is as follows.
in the second section, we describe the mde framework and demonstrate how its excellent properties inspire a partial regression model for fitting gene expression time-course data.
simulated datasets are designed for fitting by both partial mde and mle, to reveal their inherent differences.
built upon the advantages of mde and partial modelling, a robust partial regression clustering algorithm is proposed for tight clustering which naturally incorporates replication information and allows a set of scattered genes to be left out.
the experimental section is made up of two parts.
first, our proposed partial regression clustering algorithm is applied to a simulated dataset to demonstrate its effectiveness.
secondly, it is compared with some recent work by applying the methods to two well studied real datasets.
the superior performance of our algorithm is found through a carefully organized clustering validation, based on both biological knowledge and statistical indices.
in particular, a gene-ontology   <cit>  driven validation measure is proposed, specifically designed for gene expression clustering.
subsequent analysis of the clustering outcome reveals new knowledge generated by incorporating different biological resources.
this study not only explores the differences between the two estimators and the application of partial modelling, but also provide an excellent example of gene expression data mining through the combination of machine learning and biological knowledge.
owning to space restrictions, some discussions, results and elaborations have been relegated .
RESULTS
minimum distance estimation and partial modelling
minimum distance estimator
given a density function f, its corresponding parameters θ and n samples xi, i =  <dig>   <dig>  ..., n, we aim to find the optimal parameters θ^ to approximate the true parameters θ <dig> by minimizing the integrated squared difference
d, f) = ∫  <dig> dx
which gives
d, f) = ∫ f <dig> dx -  <dig> ∫ f f fdx + ∫ f <dig> dx
the last integral ∫ f <dig> dx is a constant with respect to θ, thus can be ignored.
the second integral can be obtained through kernel density estimation  <cit> .
therefore, the mde criterion simplifies to
θ^=arg⁡min⁡θ
there are many interesting features of mde.
first of all, it comes with the same robustness as all other minimum distance techniques  <cit> .
secondly, mde approximates data by making the residuals as close to normal in distribution as possible  <cit> .
these features will be further explained and illustrated in the experiments.
we will also illustrate derivation of the mde criterion for parameter estimation for our partial regression algorithm.
gaussian mixture model with mde
in principle, the finite mixture model methodology assumes that the probability density function, f, can be modelled as the sum of weighted component densities.
the weights are often constrained to have a sum of  <dig>  it is revealed later that this constraint is not necessary.
more flexible models can be obtained by relieving the system from this constraint.
a weighted gaussian mixture model has the form:
f=∑k=1kwkφ,w1+w2+…wk= <dig>
where φ is the gaussian density function, μ, σ are mean and standard deviation, k is the number of components, and wk, k =  <dig>   <dig>  ..., k are the weight parameters.
however, by relieving the constraint of ∑k=1kwk= <dig> the system can be extended for overlapping clustering inference  <cit>  since the sum of the amount of data being modelled in all clusters can exceed the total amount of data.
later, we will further prove that the amount of modelled data can also be less than the total amount of data.
in all cases, wk indicates the proportion of data points that are allocated in the kth component.
let gk  be the part in eq.
to be minimized for a k-component mixture model, we have
gk=∫f2dx−2n∑i=1nf
on the other hand,
∫φ2dx=∫2dx=12σπ∫12πσ2exp⁡222)dx=12σπ
and from section  <dig>  of  <cit>
∫φ1φ2dx=φ∫φdx=φ
by combining eq.,  and , we have
∫f2dx=∫dx=∑k=1kwk22πσk+∫∑k=1k∑l=1kwkwlφkφldx=∑k=1kwk22πσk+∑k=1k∑l=1kwkwlφ
thus from eq.
and  the distance for the k-component gaussian mixture model can be expressed as:
gk=∑k=1kwk22πσk+∑k=1k∑l=1kwkwlφ−2n∑i=1n∑k=1kwkφ
gk  is a closed-form expression, whose minimization can be performed by a standard nonlinear optimization method.
for example, a one-component model has the following mde criterion:
θ^=arg⁡min⁡θ=arg⁡min⁡θ
to further relieve the system from constraints by the weight parameters, while keeping its weighted-component structure, in the next section the idea of partial modelling is presented.
it originated from the fact that incomplete densities are allowed  <cit> , so the model will be fitted to the most relevant data.
partial mixture model with mde
the weight parameters are of particular importance in a partial mixture model.
they allow the model to estimate the component/components, while their value indicates the proportions of fitted data, so the rest of the data can be treated as scattered genes.
this approach is first described in  <cit>  for outlier detection.
it was suggested to accommodate scattered genes by forcing a large scaling parameter in one of the components in the mixture  <cit> .
however, partial modelling provides a better alternative.
although it is suggested in  <cit>  that the unconstrained mixture model can be applied for clustering, through our experiments it is clear that if the data overlap to a certain degree, all components will converge to the biggest component as a result of model freedom.
moreover, it is not practical to formulate the criterion in the form of eq.
when it comes to implementation.
instead, we solve the problem by taking advantage of the one-component model to formulate our clustering algorithm.
partial regression model
to analyse such high dimensional data as gene expression time-course measurements, a regression model with a cubic b-spline basis is set up in order to account for the inherent time dependence.
the linear regression model is capable of handling either uniformly or unevenly distributed time points, while the nonlinear spline basis helps accommodate the underlying stochastic process in the data.
the advantage of using cubic b-spline lies in that the degree of the polynomials is independent of the number of points and that curve shape is controlled locally.
let y be the variables of interest, consisting of gene expression data replicate matrices modelled as
y = α + xβ + ε
x is the design matrix consisting of a linear combination of cubic spline basis functions.
the error term ε represents the residuals taken as a weighted distribution w·n.
α, β = β <dig>  β <dig>  ..., βm, m depending on the choice of x, are the regression parameters.
as stated before, the useful feature of mde is that it fits data in such a way that the residuals are close to a normal distribution.
therefore our model is
ε = y - α - xβ
therefore, given eq.
and , the one-component pmde fit for this model has the form of
θ^=arg⁡min⁡θ=arg⁡min⁡θ
where θ = {w, α, β <dig>  ...βm, σε} and φ is the density of a normal random variable.
altogether there are m +  <dig> parameters to be estimated.
simulated datasets for pmde fitting
the main feature of our model is that it is able to identify the key component, if any, and a set of outliers, in order to find the core structure.
therefore, a feasible parameter estimator is of paramount importance.
we empirically validate our points about the nature of partial modelling and mde through fitting four simple simulated datasets.
the performance of both pmde and mle with a one-component spline regression model  is compared in terms of data fitting accuracy and robustness.
surprisingly, superior performance was achieved for the pmde fits even on such simple datasets.
all datasets are generated by sine functions, modelling cyclic behavior of genes, which are widely employed in the literature  <cit> .
gaussian noise is added to all data.
the number of knots for both spine models is chosen to be  <dig>  to allow for flexibility in curves while avoiding overfitting.
we begin with simulating the situation when the number of components k  is seriously underestimated as illustrated in figure  <dig>  three components are generated from three sine waves simulating gene expression data of three clusters, each with  <dig> time points.
the components comprise 60%, 20% and 20% of the data, respectively.
the pmde fit is highlighted by the pink line and the mle fit is blue.
pmde locates the major component, while mle is biased to all data.
this is strong evidence that pmde is superior to mle in such a scenario.
the fact that the pmde can find the key component without compromising the others suggests a solution to the vexing problem when the number of components is unknown, which is often the situation in gene expression clustering.
histograms of residuals from both fits are plotted in figure  <dig> and  <dig> to prove that pmde fit the data in such a way that the residuals are close to normal.
more datasets shown in figure 1– are used to compare the performances of pmde and mle in different scenarios.
when there are two components of entirely opposite behaviors, we can see from figure  <dig> that the mle fit is almost flat, while pmde fits the larger component .
the situation where lots of outliers are present is simulated in figure  <dig>  where the major component has 60% of the data and the rest  are generated from three different sine waves.
pmde demonstrates its robustness by capturing the major component, while mle is seriously biased.
however, in the case of two clusters of exactly equal size as shown in figure  <dig>  pmde fails, as it is designed to capture only one component but now cannot decide which one to fit.
this can be solved by using a multi-component model.
from these examples, it is observed that pmde has the ability to handle the relevant fraction of data and distinguish it from outliers, while mle blurs the distinction by accounting for all data.
this is of great value for massive datasets, when the data structure is unclear and lots of outliers are present.
the smoother fits of the proposed pmde than that of mle manifest the fact that the former is more robust against noise.
all these suggest pmde a promising tool for microarray data analysis.
interested readers are referred to additional file  <dig>  section  <dig>  for comparison of the two estimators on theoretical ground.
clustering algorithm
when analyzing gene expression time-course data, special attention needs to be paid to the following issues:
• replicates: it is desirable that the algorithm can naturally incorporate replicate information instead of simply using the mean of all replicates.
• number of clusters: the choice of k is always a problem.
the categorization of supervised and unsupervised schemes are usually determined by how k is defined.
in our unsupervised algorithm, new cluster generation automatically terminates when no new cluster can be found in the data.
• scattered genes: recently, many have proposed allowing a noisy set of genes not being clustered  <cit> .
in microarray experiments, it is generally expected that, because of the nature of data and the existence of high noise levels, many genes could show uncorrelated variations and are unrelated to the biological process under investigation.
forcing these genes into clusters will only introduce more false positives, resulting in distorted clusters and difficulty in interpretation.
apart from the aforementioned issues, like other clustering methods, the proposed algorithm needs a stopping criteria.
in this work, a statistical measure of partition quality, the calinski and harabasz  index  <cit> , is used as formulated in eq..
ch=bss/wss/
where bss() and wss() are the between-cluster and within-cluster distances defined as
bss=12∑l=1k∑xi∉cl,xj∈cld <dig>
wss=12∑l=1k∑xi,xj∈cld <dig>
cl in eq.
and  stands for the lth cluster.
the idea behind the ch measure is to compute the pairwise sum of squared errors  between clusters and compare that to the internal sum of squared errors for each cluster.
in effect, it is a measure of between-cluster dissimilarity over within-cluster dissimilarity.
the optimum clustering outcome should be the one that maximizes the ch index in eq.. the ch index was originally meant for squared euclidean distance.
since the residuals are a natural product of our spline regression model, we use the their absolute value as distance measurement in bss and wss but without the square form.
partial regression clustering algorithm
tight clustering, by definition, builds compact clusters upon an existing partition.
the initial partition, if not available, can be obtained by some empirical knowledge or heuristic clustering methods such as k-means.
given an initial partition, the clustering procedure is formulated as in algorithm  <dig>
in the initialization step of the algorithm, an existing partition of a dataset is provided as input.
the tightness threshold, υ, which controls the tightness and the number of the refined clusters produced by the algorithm as output, is defined as the reciprocal of the weighted mean variance of the clusters of the initial partition.
therefore, the greater the threshold is , the tighter the clusters become and the more clusters are formed.
the weights are determined in proportional to the size of the clusters.
in the main loop, after each new cluster is
algorithm  <dig> partial regression clustering
require: initialization
repeat
<dig>  fit partial regression model to each of the clusters;
<dig>  identify potential outliers according to a tightness threshold υ and discard them from the clusters;
<dig>  for all outliers, fit partial regression model to form a new cluster;
repeat
<dig>  for all genes re-evaluate distances to all existing spline regression models, assign them to the closest one;
<dig>  fit partial regression models to all clusters;
<dig>  calculate ch value based on current partitions;
until the clustering quality measured by ch value fails to improve.
<dig>  take the partition with highest ch value.
until no partial regression model can be fitted to the outliers.
<dig>  label all outliers as scattered genes.
generated, all data points are reassigned in the gene redistribution loop, so resultant clusters should be of reasonable size.
the rationale supporting our design is based on the features of partial modelling and robustness of the mde estimator, which we believe is able to find the relevant components in the data, while not being distracted by outliers.
the residuals, as a natural byproduct of model fitting, can be used as the distance between data points and spline regression models.
in this framework, we use deterministic class assignment during the clustering process.
stochastic relaxation or weighted assignment is regarded as more moderate than deterministic assignment.
however, it is also commonly recognised that stochastic relaxation, such as simulated annealing, does not guarantee convergence.
in fact, the selection of starting temperature or the setting of annealing schedule are often heuristic.
an initial temperature, set too high, leads to high computational cost while an initial temperature, set too low, yields similar result as deterministic relaxation but incurs higher computational cost than deterministic relaxation.
after intensive testing with stochastic and deterministic relaxation on the datasets we used, we observed that deterministic assignment strikes a better balance between computational cost and clustering accuracy.
experiment on simulated dataset
simulated datasets are necessary in evaluating the algorithm performance because the biological meanings of real datasets are very often not clear.
besides, simulated datasets provide more controllable conditions to test an algorithm.
to obtain a meaningful result, the simulated data need to share statistical characteristic with biological data.
a simulated dataset is generated from a model x = αi + βiψ  + ε , where ψ  = sin.
α, β, γ, ω are cluster-specific parameters and are chosen according to the normal distribution with mean equal to  <dig> and standard deviation  <dig>  all pattern details are listed .
ψ models the cyclic behavior of gene expression patterns.
<dig> time points are taken from  <dig> of these models, so i ∈  <dig>   <dig>  ...,  <dig>  j ∈  <dig>   <dig>  ...,  <dig>  the cluster sizes are  <dig>   <dig>   <dig>   <dig>   <dig>   <dig>  to model the noisy environment of microarray experiments, gaussian noise ε is added to all data, together with  <dig> outliers generated by adding large variance gaussian noise to three sine waves.
altogether, the simulated dataset is of size  <dig>  finally, we made some perturbations to induce more ambiguity, such as reducing the amplitude of parts of the patterns.
the clustering results are depicted in supplementary figure  <dig> of additional file  <dig>  the correct partition is achieved, with all ten outliers detected as shown in the seventh plot and the whole dataset plotted in the last one.
experiments on yeast cell cycle  dataset
a clustering method can be evaluated on theoretical grounds by internal or external validation, or both.
for internal validation, a statistical measure is preferred.
our algorithm is first validated via the ch measure in a comparison with splinecluster and mclust, two of the most popular clustering methods in the literature.
on the other hand, a measure of agreement such as the adjusted rand index   <cit>  between the resulting partition and the true partition, if known, is often used as an external validation criterion.
although a lot of evaluations for methods of the same kind are conducted in this way  <cit> , we note that there is currently no ground truth, given our knowledge of the biological structures  <cit> .
recognizing this, we set out to evaluate the performance of our algorithm through systematically finding biologically relevant evidence  <cit> .
the key to interpret a clustering outcome is to recognize the functional relationships among genes within a cluster as well as between clusters.
we first provide a quantitative measure based on the graph structure of gene ontology, then pursue biological validation and inference through go enrichment analysis in an empirically way.
clustering y <dig> dataset
a subset of  <dig> genes in yeast saccharomyces cerevisiae cell cycle  dataset  <cit>  measured at  <dig> time points was previously clustered  <cit>  into five clusters based on the first peak time in the cell cycle: early g <dig>  late g <dig>  s, g <dig> and m phase.
the original partition, as shown in supplementary figure  <dig> of additional file 1], indicates ambiguities between groups.
note that this dataset is chosen not only because it is well-studied in the gene expression clustering literature, but also because of its difficulty in terms of clustering.
the original partition makes use of only partial information of gene expression which partly explains why many clustering algorithms have poor performance .
the biological structure is still unclear, even in such heavily investigated organisms as yeast saccharomyces cerevisiae.
moreover, the average cluster size  is still far larger than desirable for efficient biological inference.
it was recently suggested that clustering based on overall profiles is preferred to the original partition on a different subset from the same dataset  <cit> .
we employ the proposed partial regression clustering algorithm to partition the y <dig> dataset into tight clusters.
by obtaining tighter clusters, we expect to obtain more informative and efficient biological inference.
the tightness threshold υ is set to  <dig> as a result of estimation during the initialization and the number of knots for the spline basis is set experimentally to  <dig> to allow flexibility of the curve without overfitting.
the clustering outcome of our algorithm is plotted in figure  <dig>  genes in the bottom right plot are the scattered genes.
the eight clusters  with scattered genes  are then cross-tabulated with the original partition in table  <dig>  the bottom row indicates the sizes of clusters of our partition and the right-most column shows those of the original partition.
the two partitions agree on many genes but also differ in a interesting way.
our partition reveals neat and easily differentiable patterns.
also, we examined the clustering outcome given by our algorithm and by other algorithms.
the left-most column contains the original partition and the top row has the resulting partition, c1–c <dig> are the eight clusters and sg are the set of scattered gene.
each number in the table except the right-most column and bottom row is the number of genes in both clusters corresponding to its row and column.
first of all, to see the effect of scattered gene detection, three algorithms are compared based on the full dataset .
by controlling a parameter in splinecluster we obtained  <dig> clusters for comparison.
the partitions of original, spinecluster and partial regression analysis are illustrated in heatmaps plotted in figure  <dig> for comparison, where an obvious improvement with respect to class distinction can be seen in the last heatmap.
the tick marks on vertical axis in each heatmap indicate where the clusters are located, while in the last heatmap the last  cluster corresponds to the scattered genes.
the second original cluster which is split into the sixth, seventh, and eighth clusters in the splinecluster partition, and the second and fifth cluster in our partition.
a closer look at the seventh and eighth cluster in the splinecluster partition shows they differ only slightly in the peak values.
however, in microarray data analysis, distinct expression patterns are more interesting than different peak values.
this is one of the reasons we use a spline model in our algorithm to capture biologically relevant information.
consider the third cluster in the splinecluster partition, which is split into the sixth and seventh clusters in our partition.
the two clusters show two entirely different patterns, one shifted from the other.
from these results, it is obvious that because of its ability in scattered gene detection, our algorithm reveals more distinguishable patterns in the data.
the set of scattered genes is listed in supplementary table  <dig> of additional file  <dig> with their annotations.
then we use the  <dig> genes , and again obtained  <dig> clusters for splinecluster.
as there is no biological knowledge input, comparison can first be conducted in a purely statistical manner, by the ch index.
mclust  <cit>  is a widely used mixture model-based clustering method.
it is unsupervised, not only in determining the number of clusters, but also in selecting the type of model that best fit the data.
the r implementation of mclust is used in our experiment.
for the 374-gene dataset it decided on the eee  model and also found  <dig> components.
our algorithm achieves the highest ch value of  <dig> , followed by  <dig>  by mclust and  <dig>  by splinecluster.
gene ontology enrichment analysis
to investigate how genes within a cluster are functionally related, and how clustering helps distinguish such functional groups, we apply gene ontology  enrichment analysis to our clustering outcome.
go terms that are likely to be over-represented in each of the clusters are identified.
these go terms are of interest because they represent the most common functions that the genes in a cluster share.
the probability that a given functional class is over-represented in the gene clusters can be estimated by using the hypergeometric distribution  <cit> .
first, for each cluster, all unique go terms that are associated with the genes in the cluster are identified.
then for each term two statistics are needed: the number of genes in the cluster that are annotated at each term and all known genes annotated at each term.
with this information, the hypergeometric distribution can be applied to identify go terms that are associated to more genes than by chance.
the probability is indicated by the resultant p-values.
using the hypergeometric distribution, suppose there are j genes annotated to a function in a total of g genes in the genome, the p-value of observing h or more genes in a cluster of size b annotated to this function is given by
p=1−∑i=0h−1/
the lower the p-value is, the more unlikely the null hypothesis that the terms appear by chance is true.
in this way, the over-represented terms are found for each cluster.
we propose within-cluster compactness  to measure the functional closeness for genes within one cluster based on the corresponding go relationship graph.
for each cluster cl, l ∈ { <dig>   <dig>  ..., k}, the most over-represented go terms tl = {t <dig>  t <dig>  ..., tnl} are found, together with their corresponding p-values pl = {p <dig>  p <dig>  ..., pnl}.
a go relationship graph gl can be plotted using tl as input, linking to their parents until the root 'biological process' is reached.
this measure aims to encourage deeper graphs with lower p-values while discouraging terms in different subgraphs with low p-values.
for example, the go graph in figure  <dig> has two big subgraphs with their node details and p-values listed in supplementary table  <dig> of additional file  <dig>  the measure should be able to represent the large distance between nodes of different subgroups  and their significance in terms of their p-values.
therefore, we define the go distance between two terms as dij = d × ) × ), where d is the shortest path between two terms in go graph and di·= d × 2) is the distance between a term and the root.
as two terms can share parents via multiple paths, the shortest distance between two terms in a go graph is defined as the shortest path by which the two terms reach a shared parent, the lowest common ancestor .
the sum of such distances for all paired go terms can be used to indicate how closely the terms are related within a cluster.
thus, within-cluster compactness for a cluster cl is defined as
wcc=∑ti∈tldi.∑ti∈tl∑tj∈tl,j≠idij
the sum of wcc for all clusters can then serve as a measure for a clustering outcome in terms of its compactness of cluster representation of biological functions.
five clustering algorithms: partial regression, splinecluster, mclust, hierarchical clustering, and k-means are compared by pooling results, using different p-value cut-offs.
using the notion of false discovery rate   <cit> , adjusted p-values are used in accordance to confidence levels, for example 2% of fdr means accepting all tests with adjusted p-values <  <dig>  as significant.
the performances of different algorithms are relatively consistent , revealing a certain robustness of this measure.
our partial regression algorithm has the highest functional cluster closeness among the five methods, indicating superior performance.
to explain what leads to such different yet consistent wcc scores and how the scores reconcile with biological knowledge, we analyse the functional categories that are statistically over-represented in the clusters.
first, we compared the over-represented terms in the resulting clusters of the proposed algorithms  and splinecluster .
for simplicity, we based the following analysis in the biological process ontology .
as indicated by the lowest p-values in each cluster, all pmde clusters have a statistically significant set of cell cycle related terms , while for sc only six out of eight clusters have such significance.
we observed that from the remaining two clusters of poorer quality , some genes involved in dna replication  were combined into pmde cluster  <dig>  resulting in a tight cluster that has a significantly functional over-representation of dna strand elongation  and other functions in dna replication.
such a high quality cluster is essential for predicting unknown functions of genes such as yhr151c and ynl058c within the cluster.
in addition, good agreement was found between known biological functions and gene clusters found by the proposed algorithm.
many clusters are significantly enriched with distinctive cell cycle relevant functions, indicating a good separation of functional clusters.
for example, cluster  <dig> has an over-representation of dna strand elongation  and cluster  <dig> is enriched with microtubule nucleation and chromosome segregation  which is crucial to chromosome division.
consistent with their biological functions, two clusters involving genes expressed in m and earlier phases reveal patterns of slightly different peak time: cluster  <dig> contains an over-representation of genes involved in dna unwinding during replication  and dna geometric change ; and cluster  <dig> is enriched with cytokinesis that is known to occur after replication and segregation of cellular components.
the two gene patterns are both biologically meaningful and statistically sound.
predictive accuracy
we compared five clustering methods: our partial regression algorithm, splinecluster, mclust  <cit> , hierarchical clustering, k-means, in terms of their predictive accuracy established in  <cit> .
since the underlying biological ground truth is unknown, evaluation of clustering algorithms for gene data cannot be carried out by similarity measures such as ari.
instead, predictive accuracy was proposed to test functional prediction accuracy from clustering.
the rationale is that since clustering is aimed at functional prediction of novel genes, if a cluster has exceptionally high occurrences of a certain gene annotation f , all genes in this cluster can be predicted to be in the functional category f. the ratio of the verified predictions to all prediction made reflects the accuracy of a clustering algorithm.
however, we have to bear in mind that this measure greatly depends on the annotation quality of the dataset under study.
since our results involved a set of scattered genes, we propose as described below a slightly different criterion to the one in  <cit> .
suppose a functional category, fi, has vi genes in a dataset of size n. if there are in total v genes belonging to functional categories f <dig>  f <dig>  ..., fm, the remaining n - v genes are denoted as 'unannotated'.
such grouping and the resulting partition c <dig>  c <dig>  ..., ck of a clustering method can be cross-tabulated to form a table.
let nij,  be the  entry of the table denoting the number of annotated genes, pij be the corresponding p-value, and n·j be the size of cluster cj.
given a threshold δ, for a k-cluster solution, its predictive accuracy a is defined as
a = pv /pc
where pv  is the verified predictions and pc is the predictions calculated by
pv=∑j=1k∑i∈{x|pxj<δ}nijpc=∑j=1k∑i∈{x|pxj<δ}n.j
supplementary table  <dig> of additional file  <dig> lists  <dig> genes in y <dig> dataset that are verified to be cell cycle related to their corresponding cell cycle phase, together with their annotations.
the  <dig> genes along with the remaining  <dig> genes denoted as 'unannotated' can then be cross-tabulated with our partition as in supplementary table  <dig> .
the bottom row of supplementary table  <dig> shows the size of clusters and the set of scattered genes.
all scattered genes are excluded from this evaluation.
by pooling results from various thresholds, we obtain a curve of 'prediction made' versus 'accuracy' for each method in comparison .
as shown in figure  <dig>  the curve for our partial regression method is above the others, indicating higher accuracy in functional group prediction.
scattered genes
another important aspect in our investigation is to study the set of scattered genes.
multiple experiments are conducted with various tightness thresholds, υ, in our partial regression method.
in supplementary table  <dig> of additional file  <dig> the set of scattered genes found in eight runs of our program with various thresholds and their annotations are presented.
their frequencies of appearance in these experiments are shown in the column feq.
.
we noticed that although these thresholds result in different numbers of clusters, the set of scattered genes hardly changes .
such consistency leads one to think about the underlying biological meaning.
as has already been pointed out  <cit> , scattered genes can be those individuals that are not relevant to the biological process under study.
however, we stress here that they can also be of significant interest, as each of them might be a key component of the cell cycle that may affect other components and indeed may be a transcription factor themselves.
therefore, its expression pattern can be uncorrelated to others in the set under study.
alternatively, a scattered gene can represent a gene whose expression is controlled by more transcription factors than the other co-regulated genes within clusters.
moreover, because the set of genes under investigation is usually selected after performing gene ranking, there may be others in the complete list that would cluster with scattered genes.
all these considerations drove us to further investigate this set of scattered genes.
among the scattered genes, five are either not well-understood or unknown for their functions.
only one of them, tip1/ybr067c, is verified to be cell cycle related in phase late g <dig>  scb regulated .
indeed, according to supplementary table  <dig>  one would conclude that all the seven genes in late g <dig>  scb regulated phase to have the same behaviour.
however, when their profiles are plotted as in figure  <dig>  we can see that tip1/ybr067c is uncorrelated to the others, making it an interesting subject for further study.
comparative evaluation on scattered gene detection
to further assess the proposed pmde's strength of scattered gene detection, the proposed algorithm is compared with a recent modification of the mclust, which allows an additional component of homogeneous poisson process for scattered genes/noise  <cit> .
the idea is for each method to filter out scattered genes and then, instead of analysing the scattered genes, compare the quality of the filtered datasets in terms of within-cluster sum of squares wss as defined in eq.. if an algorithm is stronger in outlier filtering, tighter clusters should be found in the filtered dataset, hence a smaller value of wss.
since the number of scattered genes identified by the two methods may vary, when the sets of scattered genes filtered out by different methods are of different sizes, we randomly sample a subset of the same size as the smaller set from the lager one and return the leftovers to the filtered dataset so that the filtered datasets to be investigated/clustered are of the same size.
because the clustering quality may be affected by the returned genes, we repeat the process of the random sampling of scattered genes and the clustering of the filtered dataset  <dig> times, and take the average value of wss to compare against the wss of the clustering result by the other method.
we obtain clustering results with the number of clusters k ranging from  <dig> to  <dig> for y <dig> dataset from both the pmde and the mclust.
the results are plotted in figure  <dig>  we can see that the proposed pmde performs better with large number of clusters, k, but not as good as the mclust with smaller k. however, this does not mean that the mclust outperforms the pmde because the pmde is designed to start with an initial set of clusters and iteratively split the current clusters if the splitting can lead to tighter clusters.
therefore, the clustering results by the pmde with smaller values of k are not "final" but just "provisional"; when compared to the "final" results by the mclust, the performance of the pmde appears to be inferior.
however, when the results by the pmde is more mature as k gets bigger, for example when k is greater than or equal to  <dig> as shown in figure  <dig>  the proposed pmde consistently outperforms mclust.
experiments on yeast galactose dataset
experiments are conducted on the yeast galactose dataset  <cit> , which consists of gene expression measurements in galactose utilization in saccharomyces cerevisiae.
gene expression was measured with  <dig> replicate assays across  <dig> experimental conditions .
a subset of measurements of  <dig> genes whose expression patterns reflect four functional categories in the go listings was chosen and clustered previously  <cit> .
compared with y <dig> dataset, yeast galactose dataset show more distinguishable patterns, which is easier for clustering and leads to more agreeable correlation to its functional interpretation.
for this dataset, our partial regression algorithm takes as input all  <dig> replicates of microarray data, yielding  <dig> clusters with  <dig> scattered genes when the tightness threshold is set to low value.
the four clusters  with scattered genes  are then cross-tabulated with the original partition in table  <dig>  we take  <dig> as cluster number, since it is also in accordance with prior knowledge, and obtain partitions from all five algorithms.
following, the results of wcc measure from five algorithms are plotted in figure  <dig> across different p-value cut-offs.
consistent with previous findings  <cit> , the wcc curves in figure  <dig> show that most of the algorithms performed well on this dataset.
the result from partial regression algorithm excels in both biological and statistical validation.
after the scattered genes are excluded by partial regression, the average of wcc scores across different cut-offs are  <dig> ,  <dig> ,  <dig> ,  <dig> , and  <dig> , for partial regression, spline cluster, hierarchical, k-means, and mclust, respectively.
as a mean of statistical validation, ch measure is applied to the above five algorithms, giving values of  <dig> ,  <dig> ,  <dig> ,  <dig> , and  <dig> , respectively.
the bottom row contains cluster sizes for the original partition and the right-most column contains cluster sizes for the resulting partition.
each number in the table except the right-most column and bottom row is the number of overlapping genes in both clusters corresponding to its row and column.
meanwhile, there are interesting findings from the investigation of scattered genes.
for instance, one gene  belonging to the original cluster o <dig> is classified as a scattered gene.
of the other  <dig> genes in original cluster  <dig>   <dig> are clustered into c <dig>   <dig> in c <dig>  and  <dig> in c <dig> .
the expression data of all of the  <dig> genes are plotted in figure  <dig>  revealing very different expression patterns of the  <dig> genes and the  <dig> genes differentiated by our algorithm.
both ykl152c and ymr125w are up-regulated at the beginning with down regulations for all others.
the resulting cluster c <dig> by partial regression is verified by go, since the  <dig> genes share similar annotations among the  <dig> genes in the original cluster o <dig>  for example they are all annotated to glycolysis  observed from the supplementary table  <dig> of additional file  <dig>
as an important transcription factor, ypr186c is an essential protein that binds the 5s rrna gene through the zinc finger domain and directs assembly of a multi-protein initiation complex for rna polymerase iii.
belonging to the original cluster o <dig>  ypr186c is classified as a scattered gene.
we plot its expression levels together with two other genes that are also annotated to go: <dig> , and found dramatic differences among their patterns in figure  <dig>  since this term is quite specific and it should largely reflect a gene's function, mechanisms behind such diverse behaviours are still unclear and are worth further investigations.
in summary, our algorithm receives highest wcc score.
the validity of its partitions are proved through go analysis.
we expect that its ability of scattered gene prediction will be well sought after.
CONCLUSIONS
the aim of clustering gene profiles is to find possible functional relationships among tens of thousands of genes on a microarray.
we propose that while the models for data fitting should be sensitive enough for discriminating individuals/genes, the estimators should be robust enough against noise and possible outliers.
therefore we focused on the differences between estimators by providing experimental comparisons.
the robustness of the minimum distance estimator makes it stand out in our study.
an immediate advantage is that when it is applied to gene expression clustering, it is capable of locating the key components in an unsupervised manner.
as a result, a set of scattered genes that has low correlations is naturally obtained.
besides the go enrichment analysis for the clusters from two real datasets, inference of the sets of scattered genes was also highlighted in this paper.
the partial mixture model  was known to solve problems for low dimensional data.
in fact, one problem with classical pmm is that it cannot fit data of more than  <dig> data points  <cit> .
this is the first time pmm is extended to use on high dimensional data, since current microarray experiments are having more time points and more replicates.
our contributions include introducing mde and the idea of partial modelling to gene expression research, giving comparisons with the most common estimator in the literature – maximum likelihood, and proposing a novel partial regression clustering algorithm.
our spline regression model captures the inherent time dependencies among data.
the error term is of particular importance as it can pick up the noise.
the fact that pmde estimates parameters so the residuals are as close to normal distribution as possible makes it a powerful tool for modelling the error term.
the tightness of resulting clusters can be controlled by a threshold which in a sense decides the number of clusters.
the effectiveness of the algorithm also depends on the model normality.
when model normality holds approximately, clusters can be found.
often gene expression data are transformed during pre-processing so that normality holds approximately.
although many interactions between genes are known, our knowledge of biological networks is far from complete.
no conclusion can be drawn by merely comparing clustering inference with known measure from the biological literature.
in this case, we aim to validate the algorithm and explain the clustering outcome with the help of various biological resources.
as a highlight of this paper, gene ontology clustering validation was applied to the clustering outcomes of yeast cell cycle dataset and yeast galactose dataset.
from current knowledge, it is proved that these clusters can help separate groups of genes with similar functions, while new information can be learned from exploring the go terms.
first we proposed a novel measure based on graph theory and annotation knowledge as functional compactness indication for clusters.
further, predictive accuracy was utilized to compare the annotation prediction power across several common methods.
both measures confirmed that our proposed method has the best performance.
also, gene annotations reveal new knowledge that can be derived from scattered genes.
a concern about go analysis and annotation is that lots of genes and their functions are still unknown or poorly understood.
it is our hope that through clustering, new understanding can be introduced to genome research.
in summary, the proposed system benefits from the robustness of mde to detect scattered genes, the idea of partial modelling for tight clusters, the spline regression model for capturing the expression curves at either uniformly or unevenly distributed time points, and the use of the design matrix for incorporating replicate information.
the proposed algorithm can be applied over an existing clustering to get tighter clusters.
although pmde demonstrates its effectiveness through comparisons with maximum likelihood method, it also has its limits such as relative inefficiency.
the aim of this paper is not to prove which one is better, but rather to provide analytical examples, discussions and insights.
authors' contributions
yy conceived of the study, proposed the formulae, carried out the implementation and prepared the manuscript.
c-tl revised the formulae and advised on the preparation of the manuscript.
rw advised on the preparation of the manuscript.
all authors have approved the final manuscript.
supplementary material
additional file 1
theoretical comparison between mde and mle.
theoretical comparison between minimum divergence estimator and maximum likelihood estimator.
simulated dataset patterns.
patterns for generating simulated dataset.
supplementary figure  <dig>  the resulting clusters by the partial regression clustering algorithm for the simulated dataset.
supplementary figure  <dig>  the original partition of the yeast y <dig> dataset, bottom right plot is the whole dataset.
click here for file
additional file 2
supplementary table  <dig>  the set of scattered genes for yeast y <dig> dataset.
supplementary table  <dig>  table for the go relationship graph.
supplementary table  <dig>  over-represented go terms by partial regression algorithm.
supplementary table  <dig>  over-represented go terms by splinecluster.
supplementary table  <dig>  verified cell cycle related  genes in y <dig> dataset.
supplementary table  <dig>  cross-tabulation of clustering outcome  with verified gene functional categories for yeast y <dig> dataset.
supplementary table  <dig>  over-represented terms in each original cluster for yeast galactose dataset.
click here for file
acknowledgements
we thank three anonymous referees for suggestions and comments that significantly improved this manuscript.

 based on the assumption that co-expression indicates co-regulation, gene expression data clustering aims to reveal gene groups of similar functions in the biological pathways. this biological rationale is readily supported by both empirical observations and systematic analysis. in particular, consider gene expression time-course experiments, where the data are made up of tens of thousands of genes, each with measurements taken at either uniformly or unevenly distributed time points often with several replicates. clustering algorithms provide a good initial investigation into such large-scale datasets, which ultimately leads to biological inference. an excellent review of current techniques and all subsequent analysis can be found in. various model-based methods have been proposed to accommodate the needs for data mining in such massive datasets. among them are mixed effects models and auto regressive models. the basic approach of these model-based methods is to fit a finite mixture model to the observed data, assuming that there is an underlying true model/density, and then systemically find the optimal parameters so that the fitted model/density is as close to the true model/density as possible. it is observed that model-based approaches generally achieve superior performance to many others. however, current methods can be problematic, as they often fail to show how clustering can assist in mining gene expression data. the maximum likelihood estimator is one of the most extensively used statistical estimation techniques in the literature. for a variety of models, likelihood functions, especially maximum likelihood, have been used for making inferences about parameters of the underlying probability distribution for a given dataset. the solution often involves a nonlinear optimization such as quasi-newton methods or, more commonly, expectation-maximization methods. the problem with the former method is that the quantities are estimated only when they satisfy some constraints, while with the latter method all parameters have to be explicitly specified, so the number of clusters k has to be known a priori, which is not practical in microarray data analysis. there are many unique features of mle, including its efficiency. however the practical deficiencies of mle, besides those with its optimization, are the lack of robustness against outliers and its sensitivity to the correctness of model specification. we discuss in this paper the performance of an appealing alternative, the minimum distance estimator, which is less explored in this field. inspired by the work of, we propose to incorporate mde in our algorithm for gene expression time-course analysis. mde provides robust estimation against noise and outliers, which is of particular importance in gene expression data analysis, where data are often noisy and there are few replicates. tight clustering has been proposed as a response to the needs for obtaining smaller clusters in genomic signal processing. it was motivated by the fact that the most informative clusters are very often the tight clusters, usually of size 20– genes. tight clustering refers to methods that can be built upon an existing partition to obtain core patterns that are more interpretable. the initial partition can be obtained empirically or by using generic algorithms such as k-means. as a result, more information can possibly be revealed. for example, if genes in the same functional category are allocated into different tight clusters, one may pursue the underlying explanation by looking into these clusters. one possible result of such investigation, for example, is new function discovery. in this sense, to obtain tight clusters, some genes should be classified as scattered genes, if forcing them into clusters will only disturb biologically relevant patterns. indeed, the issue of scattered genes has received more attention recently. however, in contrast to the prevailing concept that scattered genes should be treated as outliers and discarded from further study, we prove that some scattered genes can be of biological significance. current methods for gene expression time-course data rarely deal with scattered genes. to the best of our knowledge, is the first to address this issue, but it results in heavy computation due to the nature of random resampling. it was proposed in that outliers can be modelled by adding a poisson process component in the mixture model. however, this method has not been verified in this field, and it relies on correct model specification. there has been a lot of research focusing on modelling time-course data by splines and autoregressive models, usually followed by em. in, the cubic b-spline, which is a linear combination of b-spline basis functions, is used for fitting gene expression time-course data. to avoid over-fitting, it is suggested not to fit a curve to every individual gene, but to constrain the spline coefficients of co-expressed genes to have the same covariance matrix. alternatively, we propose in this work a novel approach to fit our spline model. splinecluster is an efficient hierarchical clustering program based on a regression model with a marginal likelihood criterion. starting from singleton clusters, the idea is to merge clusters based on marginal likelihood in each iteration. it is efficient and straightforward to visualize. the problem is that it overlooks microarray replication information by using only the mean of all replicates, which leads to loss of information. as microarry experiments are increasingly performed with replicates, the additional information provided by replicated measurements is a valuable source of variability in terms of effective clustering. the outline of this paper is as follows. in the second section, we describe the mde framework and demonstrate how its excellent properties inspire a partial regression model for fitting gene expression time-course data. simulated datasets are designed for fitting by both partial mde and mle, to reveal their inherent differences. built upon the advantages of mde and partial modelling, a robust partial regression clustering algorithm is proposed for tight clustering which naturally incorporates replication information and allows a set of scattered genes to be left out. the experimental section is made up of two parts. first, our proposed partial regression clustering algorithm is applied to a simulated dataset to demonstrate its effectiveness. secondly, it is compared with some recent work by applying the methods to two well studied real datasets. the superior performance of our algorithm is found through a carefully organized clustering validation, based on both biological knowledge and statistical indices. in particular, a gene-ontology driven validation measure is proposed, specifically designed for gene expression clustering. subsequent analysis of the clustering outcome reveals new knowledge generated by incorporating different biological resources. this study not only explores the differences between the two estimators and the application of partial modelling, but also provide an excellent example of gene expression data mining through the combination of machine learning and biological knowledge. owning to space restrictions, some discussions, results and elaborations have been relegated.  minimum distance estimation and partial modelling minimum distance estimator given a density function f, its corresponding parameters θ and n samples xi, i =..., n, we aim to find the optimal parameters θ^ to approximate the true parameters θ by minimizing the integrated squared difference d, f) = ∫ dx which gives d, f) = ∫ f dx - ∫ f f fdx + ∫ f dx the last integral ∫ f dx is a constant with respect to θ, thus can be ignored. the second integral can be obtained through kernel density estimation. therefore, the mde criterion simplifies to θ^=arg⁡min⁡θ there are many interesting features of mde. first of all, it comes with the same robustness as all other minimum distance techniques. secondly, mde approximates data by making the residuals as close to normal in distribution as possible. these features will be further explained and illustrated in the experiments. we will also illustrate derivation of the mde criterion for parameter estimation for our partial regression algorithm. gaussian mixture model with mde in principle, the finite mixture model methodology assumes that the probability density function, f, can be modelled as the sum of weighted component densities. the weights are often constrained to have a sum of it is revealed later that this constraint is not necessary. more flexible models can be obtained by relieving the system from this constraint. a weighted gaussian mixture model has the form: f=∑k=1kwkφ,w1+w2+…wk= where φ is the gaussian density function, μ, σ are mean and standard deviation, k is the number of components, and wk, k =..., k are the weight parameters. however, by relieving the constraint of ∑k=1kwk= the system can be extended for overlapping clustering inference since the sum of the amount of data being modelled in all clusters can exceed the total amount of data. later, we will further prove that the amount of modelled data can also be less than the total amount of data. in all cases, wk indicates the proportion of data points that are allocated in the kth component. let gk be the part in eq. to be minimized for a k-component mixture model, we have gk=∫f2dx−2n∑i=1nf on the other hand, ∫φ2dx=∫2dx=12σπ∫12πσ2exp⁡222)dx=12σπ and from section of ∫φ1φ2dx=φ∫φdx=φ by combining eq., and, we have ∫f2dx=∫dx=∑k=1kwk22πσk+∫∑k=1k∑l=1kwkwlφkφldx=∑k=1kwk22πσk+∑k=1k∑l=1kwkwlφ thus from eq. and the distance for the k-component gaussian mixture model can be expressed as: gk=∑k=1kwk22πσk+∑k=1k∑l=1kwkwlφ−2n∑i=1n∑k=1kwkφ gk is a closed-form expression, whose minimization can be performed by a standard nonlinear optimization method. for example, a one-component model has the following mde criterion: θ^=arg⁡min⁡θ=arg⁡min⁡θ to further relieve the system from constraints by the weight parameters, while keeping its weighted-component structure, in the next section the idea of partial modelling is presented. it originated from the fact that incomplete densities are allowed, so the model will be fitted to the most relevant data. partial mixture model with mde the weight parameters are of particular importance in a partial mixture model. they allow the model to estimate the component/components, while their value indicates the proportions of fitted data, so the rest of the data can be treated as scattered genes. this approach is first described in for outlier detection. it was suggested to accommodate scattered genes by forcing a large scaling parameter in one of the components in the mixture. however, partial modelling provides a better alternative. although it is suggested in that the unconstrained mixture model can be applied for clustering, through our experiments it is clear that if the data overlap to a certain degree, all components will converge to the biggest component as a result of model freedom. moreover, it is not practical to formulate the criterion in the form of eq. when it comes to implementation. instead, we solve the problem by taking advantage of the one-component model to formulate our clustering algorithm. partial regression model to analyse such high dimensional data as gene expression time-course measurements, a regression model with a cubic b-spline basis is set up in order to account for the inherent time dependence. the linear regression model is capable of handling either uniformly or unevenly distributed time points, while the nonlinear spline basis helps accommodate the underlying stochastic process in the data. the advantage of using cubic b-spline lies in that the degree of the polynomials is independent of the number of points and that curve shape is controlled locally. let y be the variables of interest, consisting of gene expression data replicate matrices modelled as y = α + xβ + ε x is the design matrix consisting of a linear combination of cubic spline basis functions. the error term ε represents the residuals taken as a weighted distribution w·n. α, β = β β..., βm, m depending on the choice of x, are the regression parameters. as stated before, the useful feature of mde is that it fits data in such a way that the residuals are close to a normal distribution. therefore our model is ε = y - α - xβ therefore, given eq. and, the one-component pmde fit for this model has the form of θ^=arg⁡min⁡θ=arg⁡min⁡θ where θ = {w, α, β...βm, σε} and φ is the density of a normal random variable. altogether there are m + parameters to be estimated. simulated datasets for pmde fitting the main feature of our model is that it is able to identify the key component, if any, and a set of outliers, in order to find the core structure. therefore, a feasible parameter estimator is of paramount importance. we empirically validate our points about the nature of partial modelling and mde through fitting four simple simulated datasets. the performance of both pmde and mle with a one-component spline regression model is compared in terms of data fitting accuracy and robustness. surprisingly, superior performance was achieved for the pmde fits even on such simple datasets. all datasets are generated by sine functions, modelling cyclic behavior of genes, which are widely employed in the literature. gaussian noise is added to all data. the number of knots for both spine models is chosen to be to allow for flexibility in curves while avoiding overfitting. we begin with simulating the situation when the number of components k is seriously underestimated as illustrated in figure three components are generated from three sine waves simulating gene expression data of three clusters, each with time points. the components comprise 60%, 20% and 20% of the data, respectively. the pmde fit is highlighted by the pink line and the mle fit is blue. pmde locates the major component, while mle is biased to all data. this is strong evidence that pmde is superior to mle in such a scenario. the fact that the pmde can find the key component without compromising the others suggests a solution to the vexing problem when the number of components is unknown, which is often the situation in gene expression clustering. histograms of residuals from both fits are plotted in figure and to prove that pmde fit the data in such a way that the residuals are close to normal. more datasets shown in figure 1– are used to compare the performances of pmde and mle in different scenarios. when there are two components of entirely opposite behaviors, we can see from figure that the mle fit is almost flat, while pmde fits the larger component. the situation where lots of outliers are present is simulated in figure where the major component has 60% of the data and the rest are generated from three different sine waves. pmde demonstrates its robustness by capturing the major component, while mle is seriously biased. however, in the case of two clusters of exactly equal size as shown in figure pmde fails, as it is designed to capture only one component but now cannot decide which one to fit. this can be solved by using a multi-component model. from these examples, it is observed that pmde has the ability to handle the relevant fraction of data and distinguish it from outliers, while mle blurs the distinction by accounting for all data. this is of great value for massive datasets, when the data structure is unclear and lots of outliers are present. the smoother fits of the proposed pmde than that of mle manifest the fact that the former is more robust against noise. all these suggest pmde a promising tool for microarray data analysis. interested readers are referred to additional file section for comparison of the two estimators on theoretical ground. clustering algorithm when analyzing gene expression time-course data, special attention needs to be paid to the following issues: • replicates: it is desirable that the algorithm can naturally incorporate replicate information instead of simply using the mean of all replicates. • number of clusters: the choice of k is always a problem. the categorization of supervised and unsupervised schemes are usually determined by how k is defined. in our unsupervised algorithm, new cluster generation automatically terminates when no new cluster can be found in the data. • scattered genes: recently, many have proposed allowing a noisy set of genes not being clustered. in microarray experiments, it is generally expected that, because of the nature of data and the existence of high noise levels, many genes could show uncorrelated variations and are unrelated to the biological process under investigation. forcing these genes into clusters will only introduce more false positives, resulting in distorted clusters and difficulty in interpretation. apart from the aforementioned issues, like other clustering methods, the proposed algorithm needs a stopping criteria. in this work, a statistical measure of partition quality, the calinski and harabasz index, is used as formulated in eq.. ch=bss/wss/ where bss() and wss() are the between-cluster and within-cluster distances defined as bss=12∑l=1k∑xi∉cl,xj∈cld wss=12∑l=1k∑xi,xj∈cld cl in eq. and stands for the lth cluster. the idea behind the ch measure is to compute the pairwise sum of squared errors between clusters and compare that to the internal sum of squared errors for each cluster. in effect, it is a measure of between-cluster dissimilarity over within-cluster dissimilarity. the optimum clustering outcome should be the one that maximizes the ch index in eq.. the ch index was originally meant for squared euclidean distance. since the residuals are a natural product of our spline regression model, we use the their absolute value as distance measurement in bss and wss but without the square form. partial regression clustering algorithm tight clustering, by definition, builds compact clusters upon an existing partition. the initial partition, if not available, can be obtained by some empirical knowledge or heuristic clustering methods such as k-means. given an initial partition, the clustering procedure is formulated as in algorithm in the initialization step of the algorithm, an existing partition of a dataset is provided as input. the tightness threshold, υ, which controls the tightness and the number of the refined clusters produced by the algorithm as output, is defined as the reciprocal of the weighted mean variance of the clusters of the initial partition. therefore, the greater the threshold is, the tighter the clusters become and the more clusters are formed. the weights are determined in proportional to the size of the clusters. in the main loop, after each new cluster is algorithm partial regression clustering require: initialization    repeat        fit partial regression model to each of the clusters;        identify potential outliers according to a tightness threshold υ and discard them from the clusters;        for all outliers, fit partial regression model to form a new cluster;       repeat           for all genes re-evaluate distances to all existing spline regression models, assign them to the closest one;           fit partial regression models to all clusters;           calculate ch value based on current partitions;       until the clustering quality measured by ch value fails to improve.        take the partition with highest ch value.    until no partial regression model can be fitted to the outliers.     label all outliers as scattered genes. generated, all data points are reassigned in the gene redistribution loop, so resultant clusters should be of reasonable size. the rationale supporting our design is based on the features of partial modelling and robustness of the mde estimator, which we believe is able to find the relevant components in the data, while not being distracted by outliers. the residuals, as a natural byproduct of model fitting, can be used as the distance between data points and spline regression models. in this framework, we use deterministic class assignment during the clustering process. stochastic relaxation or weighted assignment is regarded as more moderate than deterministic assignment. however, it is also commonly recognised that stochastic relaxation, such as simulated annealing, does not guarantee convergence. in fact, the selection of starting temperature or the setting of annealing schedule are often heuristic. an initial temperature, set too high, leads to high computational cost while an initial temperature, set too low, yields similar result as deterministic relaxation but incurs higher computational cost than deterministic relaxation. after intensive testing with stochastic and deterministic relaxation on the datasets we used, we observed that deterministic assignment strikes a better balance between computational cost and clustering accuracy. experiment on simulated dataset simulated datasets are necessary in evaluating the algorithm performance because the biological meanings of real datasets are very often not clear. besides, simulated datasets provide more controllable conditions to test an algorithm. to obtain a meaningful result, the simulated data need to share statistical characteristic with biological data. a simulated dataset is generated from a model x = αi + βiψ + ε, where ψ = sin. α, β, γ, ω are cluster-specific parameters and are chosen according to the normal distribution with mean equal to and standard deviation all pattern details are listed. ψ models the cyclic behavior of gene expression patterns. time points are taken from of these models, so i ∈..., j ∈..., the cluster sizes are to model the noisy environment of microarray experiments, gaussian noise ε is added to all data, together with outliers generated by adding large variance gaussian noise to three sine waves. altogether, the simulated dataset is of size finally, we made some perturbations to induce more ambiguity, such as reducing the amplitude of parts of the patterns. the clustering results are depicted in supplementary figure of additional file the correct partition is achieved, with all ten outliers detected as shown in the seventh plot and the whole dataset plotted in the last one. experiments on yeast cell cycle dataset a clustering method can be evaluated on theoretical grounds by internal or external validation, or both. for internal validation, a statistical measure is preferred. our algorithm is first validated via the ch measure in a comparison with splinecluster and mclust, two of the most popular clustering methods in the literature. on the other hand, a measure of agreement such as the adjusted rand index between the resulting partition and the true partition, if known, is often used as an external validation criterion. although a lot of evaluations for methods of the same kind are conducted in this way, we note that there is currently no ground truth, given our knowledge of the biological structures. recognizing this, we set out to evaluate the performance of our algorithm through systematically finding biologically relevant evidence. the key to interpret a clustering outcome is to recognize the functional relationships among genes within a cluster as well as between clusters. we first provide a quantitative measure based on the graph structure of gene ontology, then pursue biological validation and inference through go enrichment analysis in an empirically way. clustering y dataset a subset of genes in yeast saccharomyces cerevisiae cell cycle dataset measured at time points was previously clustered into five clusters based on the first peak time in the cell cycle: early g late g s, g and m phase. the original partition, as shown in supplementary figure of additional file 1], indicates ambiguities between groups. note that this dataset is chosen not only because it is well-studied in the gene expression clustering literature, but also because of its difficulty in terms of clustering. the original partition makes use of only partial information of gene expression which partly explains why many clustering algorithms have poor performance. the biological structure is still unclear, even in such heavily investigated organisms as yeast saccharomyces cerevisiae. moreover, the average cluster size is still far larger than desirable for efficient biological inference. it was recently suggested that clustering based on overall profiles is preferred to the original partition on a different subset from the same dataset. we employ the proposed partial regression clustering algorithm to partition the y dataset into tight clusters. by obtaining tighter clusters, we expect to obtain more informative and efficient biological inference. the tightness threshold υ is set to as a result of estimation during the initialization and the number of knots for the spline basis is set experimentally to to allow flexibility of the curve without overfitting. the clustering outcome of our algorithm is plotted in figure genes in the bottom right plot are the scattered genes. the eight clusters with scattered genes are then cross-tabulated with the original partition in table the bottom row indicates the sizes of clusters of our partition and the right-most column shows those of the original partition. the two partitions agree on many genes but also differ in a interesting way. our partition reveals neat and easily differentiable patterns. also, we examined the clustering outcome given by our algorithm and by other algorithms. the left-most column contains the original partition and the top row has the resulting partition, c1–c are the eight clusters and sg are the set of scattered gene. each number in the table except the right-most column and bottom row is the number of genes in both clusters corresponding to its row and column. first of all, to see the effect of scattered gene detection, three algorithms are compared based on the full dataset. by controlling a parameter in splinecluster we obtained clusters for comparison. the partitions of original, spinecluster and partial regression analysis are illustrated in heatmaps plotted in figure for comparison, where an obvious improvement with respect to class distinction can be seen in the last heatmap. the tick marks on vertical axis in each heatmap indicate where the clusters are located, while in the last heatmap the last cluster corresponds to the scattered genes. the second original cluster which is split into the sixth, seventh, and eighth clusters in the splinecluster partition, and the second and fifth cluster in our partition. a closer look at the seventh and eighth cluster in the splinecluster partition shows they differ only slightly in the peak values. however, in microarray data analysis, distinct expression patterns are more interesting than different peak values. this is one of the reasons we use a spline model in our algorithm to capture biologically relevant information. consider the third cluster in the splinecluster partition, which is split into the sixth and seventh clusters in our partition. the two clusters show two entirely different patterns, one shifted from the other. from these results, it is obvious that because of its ability in scattered gene detection, our algorithm reveals more distinguishable patterns in the data. the set of scattered genes is listed in supplementary table of additional file with their annotations. then we use the genes, and again obtained clusters for splinecluster. as there is no biological knowledge input, comparison can first be conducted in a purely statistical manner, by the ch index. mclust is a widely used mixture model-based clustering method. it is unsupervised, not only in determining the number of clusters, but also in selecting the type of model that best fit the data. the r implementation of mclust is used in our experiment. for the 374-gene dataset it decided on the eee model and also found components. our algorithm achieves the highest ch value of, followed by by mclust and by splinecluster. gene ontology enrichment analysis to investigate how genes within a cluster are functionally related, and how clustering helps distinguish such functional groups, we apply gene ontology enrichment analysis to our clustering outcome. go terms that are likely to be over-represented in each of the clusters are identified. these go terms are of interest because they represent the most common functions that the genes in a cluster share. the probability that a given functional class is over-represented in the gene clusters can be estimated by using the hypergeometric distribution. first, for each cluster, all unique go terms that are associated with the genes in the cluster are identified. then for each term two statistics are needed: the number of genes in the cluster that are annotated at each term and all known genes annotated at each term. with this information, the hypergeometric distribution can be applied to identify go terms that are associated to more genes than by chance. the probability is indicated by the resultant p-values. using the hypergeometric distribution, suppose there are j genes annotated to a function in a total of g genes in the genome, the p-value of observing h or more genes in a cluster of size b annotated to this function is given by p=1−∑i=0h−1/ the lower the p-value is, the more unlikely the null hypothesis that the terms appear by chance is true. in this way, the over-represented terms are found for each cluster. we propose within-cluster compactness to measure the functional closeness for genes within one cluster based on the corresponding go relationship graph. for each cluster cl, l ∈ {..., k}, the most over-represented go terms tl = {t t..., tnl} are found, together with their corresponding p-values pl = {p p..., pnl}. a go relationship graph gl can be plotted using tl as input, linking to their parents until the root'biological process' is reached. this measure aims to encourage deeper graphs with lower p-values while discouraging terms in different subgraphs with low p-values. for example, the go graph in figure has two big subgraphs with their node details and p-values listed in supplementary table of additional file the measure should be able to represent the large distance between nodes of different subgroups and their significance in terms of their p-values. therefore, we define the go distance between two terms as dij = d × ) × ), where d is the shortest path between two terms in go graph and di·= d × 2) is the distance between a term and the root. as two terms can share parents via multiple paths, the shortest distance between two terms in a go graph is defined as the shortest path by which the two terms reach a shared parent, the lowest common ancestor. the sum of such distances for all paired go terms can be used to indicate how closely the terms are related within a cluster. thus, within-cluster compactness for a cluster cl is defined as wcc=∑ti∈tldi.∑ti∈tl∑tj∈tl,j≠idij the sum of wcc for all clusters can then serve as a measure for a clustering outcome in terms of its compactness of cluster representation of biological functions. five clustering algorithms: partial regression, splinecluster, mclust, hierarchical clustering, and k-means are compared by pooling results, using different p-value cut-offs. using the notion of false discovery rate, adjusted p-values are used in accordance to confidence levels, for example 2% of fdr means accepting all tests with adjusted p-values < as significant. the performances of different algorithms are relatively consistent, revealing a certain robustness of this measure. our partial regression algorithm has the highest functional cluster closeness among the five methods, indicating superior performance. to explain what leads to such different yet consistent wcc scores and how the scores reconcile with biological knowledge, we analyse the functional categories that are statistically over-represented in the clusters. first, we compared the over-represented terms in the resulting clusters of the proposed algorithms and splinecluster. for simplicity, we based the following analysis in the biological process ontology. as indicated by the lowest p-values in each cluster, all pmde clusters have a statistically significant set of cell cycle related terms, while for sc only six out of eight clusters have such significance. we observed that from the remaining two clusters of poorer quality, some genes involved in dna replication were combined into pmde cluster resulting in a tight cluster that has a significantly functional over-representation of dna strand elongation and other functions in dna replication. such a high quality cluster is essential for predicting unknown functions of genes such as yhr151c and ynl058c within the cluster. in addition, good agreement was found between known biological functions and gene clusters found by the proposed algorithm. many clusters are significantly enriched with distinctive cell cycle relevant functions, indicating a good separation of functional clusters. for example, cluster has an over-representation of dna strand elongation and cluster is enriched with microtubule nucleation and chromosome segregation which is crucial to chromosome division. consistent with their biological functions, two clusters involving genes expressed in m and earlier phases reveal patterns of slightly different peak time: cluster contains an over-representation of genes involved in dna unwinding during replication and dna geometric change; and cluster is enriched with cytokinesis that is known to occur after replication and segregation of cellular components. the two gene patterns are both biologically meaningful and statistically sound. predictive accuracy we compared five clustering methods: our partial regression algorithm, splinecluster, mclust, hierarchical clustering, k-means, in terms of their predictive accuracy established in. since the underlying biological ground truth is unknown, evaluation of clustering algorithms for gene data cannot be carried out by similarity measures such as ari. instead, predictive accuracy was proposed to test functional prediction accuracy from clustering. the rationale is that since clustering is aimed at functional prediction of novel genes, if a cluster has exceptionally high occurrences of a certain gene annotation f, all genes in this cluster can be predicted to be in the functional category f. the ratio of the verified predictions to all prediction made reflects the accuracy of a clustering algorithm. however, we have to bear in mind that this measure greatly depends on the annotation quality of the dataset under study. since our results involved a set of scattered genes, we propose as described below a slightly different criterion to the one in. suppose a functional category, fi, has vi genes in a dataset of size n. if there are in total v genes belonging to functional categories f f..., fm, the remaining n - v genes are denoted as'unannotated'. such grouping and the resulting partition c c..., ck of a clustering method can be cross-tabulated to form a table. let nij, be the entry of the table denoting the number of annotated genes, pij be the corresponding p-value, and n·j be the size of cluster cj. given a threshold δ, for a k-cluster solution, its predictive accuracy a is defined as a = pv /pc where pv is the verified predictions and pc is the predictions calculated by pv=∑j=1k∑i∈{x|pxj<δ}nijpc=∑j=1k∑i∈{x|pxj<δ}n.j supplementary table of additional file lists genes in y dataset that are verified to be cell cycle related to their corresponding cell cycle phase, together with their annotations. the genes along with the remaining genes denoted as'unannotated' can then be cross-tabulated with our partition as in supplementary table. the bottom row of supplementary table shows the size of clusters and the set of scattered genes. all scattered genes are excluded from this evaluation. by pooling results from various thresholds, we obtain a curve of'prediction made' versus'accuracy' for each method in comparison. as shown in figure the curve for our partial regression method is above the others, indicating higher accuracy in functional group prediction. scattered genes another important aspect in our investigation is to study the set of scattered genes. multiple experiments are conducted with various tightness thresholds, υ, in our partial regression method. in supplementary table of additional file the set of scattered genes found in eight runs of our program with various thresholds and their annotations are presented. their frequencies of appearance in these experiments are shown in the column feq.. we noticed that although these thresholds result in different numbers of clusters, the set of scattered genes hardly changes. such consistency leads one to think about the underlying biological meaning. as has already been pointed out, scattered genes can be those individuals that are not relevant to the biological process under study. however, we stress here that they can also be of significant interest, as each of them might be a key component of the cell cycle that may affect other components and indeed may be a transcription factor themselves. therefore, its expression pattern can be uncorrelated to others in the set under study. alternatively, a scattered gene can represent a gene whose expression is controlled by more transcription factors than the other co-regulated genes within clusters. moreover, because the set of genes under investigation is usually selected after performing gene ranking, there may be others in the complete list that would cluster with scattered genes. all these considerations drove us to further investigate this set of scattered genes. among the scattered genes, five are either not well-understood or unknown for their functions. only one of them, tip1/ybr067c, is verified to be cell cycle related in phase late g scb regulated. indeed, according to supplementary table one would conclude that all the seven genes in late g scb regulated phase to have the same behaviour. however, when their profiles are plotted as in figure we can see that tip1/ybr067c is uncorrelated to the others, making it an interesting subject for further study. comparative evaluation on scattered gene detection to further assess the proposed pmde's strength of scattered gene detection, the proposed algorithm is compared with a recent modification of the mclust, which allows an additional component of homogeneous poisson process for scattered genes/noise. the idea is for each method to filter out scattered genes and then, instead of analysing the scattered genes, compare the quality of the filtered datasets in terms of within-cluster sum of squares wss as defined in eq.. if an algorithm is stronger in outlier filtering, tighter clusters should be found in the filtered dataset, hence a smaller value of wss. since the number of scattered genes identified by the two methods may vary, when the sets of scattered genes filtered out by different methods are of different sizes, we randomly sample a subset of the same size as the smaller set from the lager one and return the leftovers to the filtered dataset so that the filtered datasets to be investigated/clustered are of the same size. because the clustering quality may be affected by the returned genes, we repeat the process of the random sampling of scattered genes and the clustering of the filtered dataset times, and take the average value of wss to compare against the wss of the clustering result by the other method. we obtain clustering results with the number of clusters k ranging from to for y dataset from both the pmde and the mclust. the results are plotted in figure we can see that the proposed pmde performs better with large number of clusters, k, but not as good as the mclust with smaller k. however, this does not mean that the mclust outperforms the pmde because the pmde is designed to start with an initial set of clusters and iteratively split the current clusters if the splitting can lead to tighter clusters. therefore, the clustering results by the pmde with smaller values of k are not "final" but just "provisional"; when compared to the "final" results by the mclust, the performance of the pmde appears to be inferior. however, when the results by the pmde is more mature as k gets bigger, for example when k is greater than or equal to as shown in figure the proposed pmde consistently outperforms mclust. experiments on yeast galactose dataset experiments are conducted on the yeast galactose dataset, which consists of gene expression measurements in galactose utilization in saccharomyces cerevisiae. gene expression was measured with replicate assays across experimental conditions. a subset of measurements of genes whose expression patterns reflect four functional categories in the go listings was chosen and clustered previously. compared with y dataset, yeast galactose dataset show more distinguishable patterns, which is easier for clustering and leads to more agreeable correlation to its functional interpretation. for this dataset, our partial regression algorithm takes as input all replicates of microarray data, yielding clusters with scattered genes when the tightness threshold is set to low value. the four clusters with scattered genes are then cross-tabulated with the original partition in table we take as cluster number, since it is also in accordance with prior knowledge, and obtain partitions from all five algorithms. following, the results of wcc measure from five algorithms are plotted in figure across different p-value cut-offs. consistent with previous findings, the wcc curves in figure show that most of the algorithms performed well on this dataset. the result from partial regression algorithm excels in both biological and statistical validation. after the scattered genes are excluded by partial regression, the average of wcc scores across different cut-offs are,,,, and, for partial regression, spline cluster, hierarchical, k-means, and mclust, respectively. as a mean of statistical validation, ch measure is applied to the above five algorithms, giving values of,,,, and, respectively. the bottom row contains cluster sizes for the original partition and the right-most column contains cluster sizes for the resulting partition. each number in the table except the right-most column and bottom row is the number of overlapping genes in both clusters corresponding to its row and column. meanwhile, there are interesting findings from the investigation of scattered genes. for instance, one gene belonging to the original cluster o is classified as a scattered gene. of the other genes in original cluster are clustered into c in c and in c. the expression data of all of the genes are plotted in figure revealing very different expression patterns of the genes and the genes differentiated by our algorithm. both ykl152c and ymr125w are up-regulated at the beginning with down regulations for all others. the resulting cluster c by partial regression is verified by go, since the genes share similar annotations among the genes in the original cluster o for example they are all annotated to glycolysis observed from the supplementary table of additional file as an important transcription factor, ypr186c is an essential protein that binds the 5s rrna gene through the zinc finger domain and directs assembly of a multi-protein initiation complex for rna polymerase iii. belonging to the original cluster o ypr186c is classified as a scattered gene. we plot its expression levels together with two other genes that are also annotated to go:, and found dramatic differences among their patterns in figure since this term is quite specific and it should largely reflect a gene's function, mechanisms behind such diverse behaviours are still unclear and are worth further investigations. in summary, our algorithm receives highest wcc score. the validity of its partitions are proved through go analysis. we expect that its ability of scattered gene prediction will be well sought after.  the aim of clustering gene profiles is to find possible functional relationships among tens of thousands of genes on a microarray. we propose that while the models for data fitting should be sensitive enough for discriminating individuals/genes, the estimators should be robust enough against noise and possible outliers. therefore we focused on the differences between estimators by providing experimental comparisons. the robustness of the minimum distance estimator makes it stand out in our study. an immediate advantage is that when it is applied to gene expression clustering, it is capable of locating the key components in an unsupervised manner. as a result, a set of scattered genes that has low correlations is naturally obtained. besides the go enrichment analysis for the clusters from two real datasets, inference of the sets of scattered genes was also highlighted in this paper. the partial mixture model was known to solve problems for low dimensional data. in fact, one problem with classical pmm is that it cannot fit data of more than data points. this is the first time pmm is extended to use on high dimensional data, since current microarray experiments are having more time points and more replicates. our contributions include introducing mde and the idea of partial modelling to gene expression research, giving comparisons with the most common estimator in the literature – maximum likelihood, and proposing a novel partial regression clustering algorithm. our spline regression model captures the inherent time dependencies among data. the error term is of particular importance as it can pick up the noise. the fact that pmde estimates parameters so the residuals are as close to normal distribution as possible makes it a powerful tool for modelling the error term. the tightness of resulting clusters can be controlled by a threshold which in a sense decides the number of clusters. the effectiveness of the algorithm also depends on the model normality. when model normality holds approximately, clusters can be found. often gene expression data are transformed during pre-processing so that normality holds approximately. although many interactions between genes are known, our knowledge of biological networks is far from complete. no conclusion can be drawn by merely comparing clustering inference with known measure from the biological literature. in this case, we aim to validate the algorithm and explain the clustering outcome with the help of various biological resources. as a highlight of this paper, gene ontology clustering validation was applied to the clustering outcomes of yeast cell cycle dataset and yeast galactose dataset. from current knowledge, it is proved that these clusters can help separate groups of genes with similar functions, while new information can be learned from exploring the go terms. first we proposed a novel measure based on graph theory and annotation knowledge as functional compactness indication for clusters. further, predictive accuracy was utilized to compare the annotation prediction power across several common methods. both measures confirmed that our proposed method has the best performance. also, gene annotations reveal new knowledge that can be derived from scattered genes. a concern about go analysis and annotation is that lots of genes and their functions are still unknown or poorly understood. it is our hope that through clustering, new understanding can be introduced to genome research. in summary, the proposed system benefits from the robustness of mde to detect scattered genes, the idea of partial modelling for tight clusters, the spline regression model for capturing the expression curves at either uniformly or unevenly distributed time points, and the use of the design matrix for incorporating replicate information. the proposed algorithm can be applied over an existing clustering to get tighter clusters. although pmde demonstrates its effectiveness through comparisons with maximum likelihood method, it also has its limits such as relative inefficiency. the aim of this paper is not to prove which one is better, but rather to provide analytical examples, discussions and insights. authors' contributions yy conceived of the study, proposed the formulae, carried out the implementation and prepared the manuscript. c-tl revised the formulae and advised on the preparation of the manuscript. rw advised on the preparation of the manuscript. all authors have approved the final manuscript. supplementary material additional file 1 theoretical comparison between mde and mle. theoretical comparison between minimum divergence estimator and maximum likelihood estimator. simulated dataset patterns. patterns for generating simulated dataset. supplementary figure the resulting clusters by the partial regression clustering algorithm for the simulated dataset. supplementary figure the original partition of the yeast y dataset, bottom right plot is the whole dataset. click here for file additional file 2 supplementary table the set of scattered genes for yeast y dataset. supplementary table table for the go relationship graph. supplementary table over-represented go terms by partial regression algorithm. supplementary table over-represented go terms by splinecluster. supplementary table verified cell cycle related genes in y dataset. supplementary table cross-tabulation of clustering outcome with verified gene functional categories for yeast y dataset. supplementary table over-represented terms in each original cluster for yeast galactose dataset. click here for file acknowledgements we thank three anonymous referees for suggestions and comments that significantly improved this manuscript.
our proposed partial regression clustering algorithm scores top in gene ontology driven evaluation, in comparison with four other popular clustering algorithms. to apply minimum distance estimation, a partial mixture model that can naturally incorporate replicate information and allow scattered genes is formulated. for the first time partial mixture model is successfully extended to time-course data analysis. the robustness of our partial regression clustering algorithm proves the suitability of the combination of both partial mixture model and minimum distance estimator in this field. on the other hand, there has been extensive use of maximum likelihood techniques for model parameter estimation. in this paper we show the inherent robustness of the minimum distance estimator that makes it a powerful tool for parameter estimation in model-based time-course clustering. we provide experimental results of simulated data fitting, where the minimum distance estimator demonstrates superior performance to the maximum likelihood estimator. we show that tight clustering not only is capable to generate more profound understanding of the dataset under study well in accordance to established biological knowledge, but also presents interesting new hypotheses during interpretation of clustering results. tight clustering arose recently from a desire to obtain tighter and potentially more informative clusters in gene expression studies. by contrast, the minimum distance estimator has been largely ignored.
 tight clustering arose recently from a desire to obtain tighter and potentially more informative clusters in gene expression studies. scattered genes with relatively loose correlations should be excluded from the clusters. however, in the literature there is little work dedicated to this area of research. on the other hand, there has been extensive use of maximum likelihood techniques for model parameter estimation. by contrast, the minimum distance estimator has been largely ignored.  in this paper we show the inherent robustness of the minimum distance estimator that makes it a powerful tool for parameter estimation in model-based time-course clustering. to apply minimum distance estimation, a partial mixture model that can naturally incorporate replicate information and allow scattered genes is formulated. we provide experimental results of simulated data fitting, where the minimum distance estimator demonstrates superior performance to the maximum likelihood estimator. both biological and statistical validations are conducted on a simulated dataset and two real gene expression datasets. our proposed partial regression clustering algorithm scores top in gene ontology driven evaluation, in comparison with four other popular clustering algorithms.  for the first time partial mixture model is successfully extended to time-course data analysis. the robustness of our partial regression clustering algorithm proves the suitability of the combination of both partial mixture model and minimum distance estimator in this field. we show that tight clustering not only is capable to generate more profound understanding of the dataset under study well in accordance to established biological knowledge, but also presents interesting new hypotheses during interpretation of clustering results. in particular, we provide biological evidences that scattered genes can be relevant and are interesting subjects for study, in contrast to prevailing opinion. 