An improved method of Wu's functional clustering

The concept of functional clustering was first systematically introduced by Rongling Wu and his team^[1] to classify periodically expressed genes into different patterns through Fourier series expansions. This approach, particularly named Wu's functional clustering (or Wu's FunClu), was built on a mixture-based likelihood function containing weighted mixture components. Each of these components is assumed to follow a multivariate (longitudinal) normal distribution with a mean vector modeled by biologically meaningful mathematical equations, such as Fourier's periodic equations, and a (co)variance matrix modeled by an autoregressive process. Building on this foundation, Wu and his team^[2] further expanded functional clustering to multiple dimensions of data, showing the important application and uniqueness of this approach to dissect and resolve complex systems.

Although there is a rich body of literature on clustering algorithms^[3−5], Wu's FunClu has proved to be more powerful, versatile, and flexible for data clustering. Unlike conventional approaches that treat observations as static vectors, it explicitly considers the dynamic and temporal structure of the data. Its applications are especially important in analyzing omics datasets, including transcriptomics, proteomics, metabolomics, and single-cell data^[2,6,7], where clustering can reveal hidden biological modules and functional relationships. A recent study by Wu et al.^[2] applied functional clustering to large-scale metabolomics data. Metabolites in healthy and diseased groups were grouped into coherent modules, providing new insights into metabolic reprogramming in inflammatory bowel disease^[2]. These examples highlight the broad potential of functional clustering as a versatile method that links advanced statistical modeling with key challenges in biomedicine and other scientific fields.

Modern transcriptome and multi-omics studies often involve tens of thousands to over hundreds of thousands of features, and the increasing sample size presents a significant challenge for data processing^[8,9]. Traditional functional clustering methods struggle with large datasets because of the high computational demands, high memory requirements, and the risk of misclassification, leading to less effective clustering results. Additionally, time-series data are often limited or unavailable in biological research, which hinders the use of temporal clustering methods. Moreover, multi-omics datasets frequently vary in scale, noise characteristics, and correlation structures, requiring flexible and robust clustering techniques^[10−12]. These factors highlight the need for improved functional clustering approaches that can efficiently handle high-dimensional, nontemporal, and heterogeneous biological data while providing meaningful biological insights.

In this study, we propose an improved framework of Wu's FunClu to leverage its capacity for clustering large-scale, multi-omics data. The framework combines the allometric scaling law, advanced optimization techniques, and enhanced covariance modeling approaches to produce accurate, robust, and biologically meaningful clustering results. This framework is designed to be scalable and adaptable, offering a versatile tool for analyzing complex biological datasets.

The rest of the paper is organized as follows. Section "Methods" details the proposed methodology. Section "Silulation" shows the simulation results, and section "Data analysis" demonstrates an application to real biological data. The last section concludes the paper with a discussion and future directions.

We assume that the dataset consists of $ n $ observations, each represented by an $ m $-dimensional profile $ {y}_{i}=({y}_{i1},\cdots ,{y}_{im}{)}^{T} $. The structure of Wu's FunClu is modeled using a finite Gaussian mixture distribution with $ K $ components as follows:

where, $ {\alpha }_{k} $ is the mixing proportion satisfying $ \sum \nolimits_{k=1}^{K}{\alpha }_{k}=1 $, and $ \mathcal{N}(\cdot ) $ denotes the multivariate Gaussian density with a mean vector $ {\mu }_{ki} $and a covariance matrix ${\Sigma }_{k}$.

Traditional approaches are often constructed under the assumption that sufficient temporal or dynamic data are available. Such data enable the direct modeling of time-dependent trajectories. However, in many practical applications, longitudinal measurements are either infeasible to obtain or require excessive experimental and financial resources. This limitation motivates the need for alternative approaches that can extract dynamic-like information from static high-dimensional data.

The study of complex biological systems at the cellular and molecular levels has increasingly drawn parallels to community ecology^[13−15]. Just as ecology examines the interactions between species and their environment, systems biology aims to model how the information flow between molecular components gives rise to the emergent behavior of the system^[16].

A central pillar of this "crosspollination" is the concept of the ecological niche, which characterizes the response of organisms to their environment across time and space^[17]. Within this niche space, a species' survival and growth are fundamentally governed by resource availability^[18,19]. In natural ecosystems, resource availability is the net result of the interactions between all resident organisms and the environmental supply. Consequently, the "net quality" or productivity of such an ecosystem can be summarized by the ecosystem index, a surrogate measure derived from the total abundance of all organisms that reflects the environmental factors essential for the system's survival^[20].

Analogously, we treat a sample as a gene ecosystem. In this framework, individual genes are equivalent to organisms, and their expression levels reflect their utilization of the "resources" provided by the cellular environment. Following the notion of the ecosystem index used in ecology and agriculture^[21,22], we define the habitat index $ {E}_{i} $ for each sample $ i $ as

where, $ {y}_{ij} $ represents the expression level of gene $ j $ in sample $ i $.

In complex living systems, the variation of a specific component in response to the state of the whole system is rarely linear. Instead, these relationships are often observed to obey the allometric scaling law, a fundamental biological principle typically described by a power equation^[23−25].

The term "allometry" was first coined by Huxley and Teissier^[26] to investigate the phenomenon of relative growth, specifically how the proportions of an organism's parts change in relation to its total size. Since then, the scope of allometry has been broadly expanded to refer to almost any covarying biological measurements where the part's behavior is dictated by the scale of the whole.

As mentioned above, by conceptualizing each sample as a complex ecosystem of interconnected entities, we use the habitat index—defined as the cumulative abundance of all entities within a specific spatiotemporal context—as a robust proxy for assessing an ecosystem. The observed spectrum of habitat indices captures the systemic variability of gene ecosystems, likely driven by intricate internal interactions and external environmental influences. This framework effectively formalizes the part–whole relationship between individual entities' abundance and the collective habitat index across diverse samples. This methodological approach has already gained significant traction in the biological sciences. Its applications range from monitoring the functional stability of terrestrial soil microbiomes to deciphering the progressive decoupling of gene interactions during biological aging^[27,28].

Translating this principle into the clustering model, the mean function of the $ k $-th cluster is expressed as a scaling function of the habitat index

where, $ {a}_{k} $ and $ {b}_{k} $ are cluster-specific parameters. This representation allows the model to reveal scaling patterns between local observations and the global structure of the dataset, thus providing a biologically meaningful yet statistically flexible formulation for functional clustering.

The specification of the covariance matrix ${\Sigma }_{k}$ is critical for clustering performance. High-dimensional biological data often exhibit complex dependencies, and using an unrestricted covariance matrix may result in excessive free parameters, high computational cost, and numerical instability. To address this, we use a structured covariance model based on an autoregressive formulation.

Let $ {E}_{i} $ denote the expression index of sample $ i $. The variance of sample $ i $ and the covariance between samples $ i $ and $ j $ are given^[29−31] by

where, $ {\nu }^{2} $ represents the measurement error variance and $ \phi \in (-1{,}1) $ is a correlation parameter controlling the strength of dependence. With $ d $ samples, the covariance matrix $\Sigma$ is of dimension $ d\times d $ and admits the form

Given a dataset $ X=\{{x}_{1},{x}_{2},\cdots ,{x}_{n}\} $, the log-likelihood function can be expressed as follows:

Introducing latent variables that indicate cluster membership, the Expetation-Maximization (EM) auxiliary function can be written as follows:

with the posterior responsibility

To estimate the model parameters, we use a hybrid EM algorithm characterized by a dual-optimization strategy in the M-step. Specifically, for parameters with explicit gradients (e.g., mixing proportions), we derive closed-form analytic solutions to ensure computational efficiency. For parameters lacking such analytic properties, we incorporate numerical optimization routines to iteratively maximize the auxiliary function. This hybrid approach balances analytical precision with numerical flexibility, ensuring robust convergence in complex parameter spaces.

The mixing proportions admit the following analytic solution:

However, in practice, it is possible that some clusters receive nearly zero assignments, leading to degeneracy (i.e., "empty clusters"). To mitigate this, we apply a smoothing transformation using the softmax function as follows:

where, $ \gamma $ is a temperature parameter controlling the degree of smoothing. This adjustment ensures that all components retain nonzero mixing weights, thereby preventing the collapse of the mixture model and improving stability.

For the mean function parameters $ ({a}_{k},{b}_{k}) $ and covariance parameters $ ({\phi }_{k},{v}_{k}) $, closed-form solutions are not available. We therefore adopt a gradient-based numerical optimization strategy. Specifically, the Adam algorithm is used to maximize $ Q $ with respect to these parameters^[32].

This mixed update strategy preserves analytic efficiency for simple parameters while leveraging modern optimization for complex components, thereby improving the convergence stability and scalability in high-dimensional applications.

Initialization of parameters is critical for the stability and convergence of the EM algorithm. In our framework, the mean function parameters are initialized using the results of a k-means clustering applied to the raw data matrix. Specifically, the empirical cluster centroids from k-means are taken as the initial estimates for the functional means, which are then projected onto the parametric form $ {\mu }_{k}(E)={a}_{k}{E}^{{{b}_{k}}} $.

The covariance parameters are initialized using identity-based structures with small perturbations to avoid singularity, whereas the mixing proportions are initialized uniformly, i.e., $ \alpha _{k}^{(0)}=1/K $.

This initialization strategy provides a reasonable approximation to the underlying structure of the data, effectively reducing the risk of poor local optima during subsequent EM iterations.

The optimal number of clusters $ K $ is determined by minimizing the Bayesian information criterion (BIC)^[33]. Specifically, for a given model with the likelihood $ L(X\mid {\hat{\Theta }}_{K}) $ and the parameter dimension $ {d}_{K} $, the BIC is defined as

where, $ n $ is the sample size.

We evaluate the BIC across a grid of candidate cluster numbers and select the value of $ K $ that minimizes the criterion. This approach balances the model's fit and complexity, thereby preventing overfitting and ensuring biologically interpretable results.

Multivariate extension of functional clustering

In many real-world applications, data are collected across multiple conditions or tissues, which requires extending the functional clustering model to a multivariate setting. Suppose the dataset consists of paired vectors $ {\mathbf{x}}_{i}=({x}_{i1},\cdots ,{x}_{ip}) $ and $ {\mathbf{y}}_{i}=({y}_{i1},\cdots ,{y}_{iq}) $, representing gene expression trajectories measured under two (or more) conditions. For cluster $ k $, the mean functions are denoted as $ \mu _{k}^{x}(t)=f(\theta _{k}^{x};t) $ and $ \mu _{k}^{y}(t)=f(\theta _{k}^{y};t) $, where $ f(\cdot ) $ is the allometric growth function and $ \theta _{k}^{x},\theta _{k}^{y} $ are cluster-specific parameters.

Assuming a shared covariance structure $ \Sigma $, the joint distribution of observations is given by

where, $ {\alpha }_{k} $ represents the mixing proportions. The corresponding log-likelihood is

Parameter estimation proceeds via the EM algorithm described in section "Methods", with the possibility of smoothing the mixing proportions to avoid degenerate (empty) clusters. This formulation naturally generalizes to more than two tissues or conditions.

To evaluate the performance of the proposed EM algorithm, we conducted a series of simulation studies based on synthetic datasets. Each dataset consisted of $ n=500 $ samples, each with $ p=50 $ features, partitioned into $ k=20 $ latent clusters.

For each cluster $ c\in \{1,\cdots ,k\} $, the mean vector was specified according to a log-linear function as follows:

where, the intercepts $ {a}_{c} $ were drawn independently from the uniform distribution $ U(0.5{,}5) $, and the slopes $ {b}_{c} $ were sampled from $ U(0.1{,}3) $. This formulation ensures heterogeneous mean structures across clusters while maintaining a consistent log-linear growth pattern.

Conditional on the cluster assignment $ {y}_{i} $, the feature vector for sample $ i $ was

The covariance structure is assumed to be common across all components and follows the type form given in Eq. (4).

We considered two correlation settings, $ \phi =0.6 $ and $ \phi =0.9 $, and for each setting, we varied the noise variance across five levels: $ v\in \{0.2{,}\;0.4{,}\;0.6{,}\;0.8{,}\;1.0\} $. For each combination of parameters, we generated 20 independent replicates to reduce Monte Carlo variability.

Methods compared

We compare four algorithms:

FunClu: The full functional clustering pipeline described in section "Methods". In the variant FunClu, we do not apply softmax smoothing to the updated mixing proportions $ {\alpha }_{k} $.

FunClu smoothing: The same as FunClu except that after the analytic update of $ {\alpha }_{k} $ in the M-step, we apply a softmax smoothing given in Eq. (10).

This renormalization avoids zero weights and shrinks extreme small weights toward a nonzero baseline; in practice, it markedly reduces the occurrence of empty clustersand improves numeric stability when some components receive few or zero responsibilities.

We report the results for both variants to demonstrate the concrete effect of smoothing.

GMM: The standard Gaussian mixture model with full covariance, fitted by the EM. We use the same K (true $ k $) for a fair comparison.

KMeans: The standard k-means on the raw $ X $ matrix, with $ k $ set to the true number of clusters.

For all methods, we evaluate clustering recovery relative to the true labels.

Evaluation metrics

We assess clustering performance using the following metrics.

(1) Clustering accuracy (CA): The proportion of correctly classified samples under optimal label matching.

(2) Adjusted Rand index (ARI): This accounts for chance agreement in clustering.

(3) Normalized mutual information (NMI): This measures information-theoretic similarity.

Each experiment was replicated 50 times to reduce variability.

Simulation results

Tables 1 and 2 report the averaged performance (mean ± standard deviation [SD] over the replicates) of four methods (FunClu, FunClu Smoothing, GMM, and KMeans) under the noise levels $ v\in \{0.2{,}\;0.4{,}\;0.6{,}\;0.8{,}\;1.0\} $ for two correlation settings $ \phi =0.6 $ and $ \phi =0.9 $. The main findings are as follows.

First, FunClu generally outperforms the baseline methods in terms of CA, ARI, and NMI. This improvement is especially pronounced at low to moderate noise levels ($ \gamma \leq 0.6 $), highlighting the effectiveness of modeling the mean structure via the allometric scaling law in capturing the underlying log-linear relationships.

Second, the smoothing strategy (FunClu_Smoothing) adds stability. Although FunClu alone often achieves the highest numerical scores, with several parameter settings, it suffers from empty clusters, which undermines the interpretability of the clustering solution. FunClu_Smoothing, by contrast, prevents class loss by regularizing the mixing proportions through a softmax transformation, ensuring that all clusters are retained. Although the smoothing variant sometimes yields slightly lower CA, ARI, or NMI values, its solutions remain more reliable and biologically interpretable.

Third, as expected, performance decreases as the noise variance increases. All methods exhibit declining accuracy and agreement metrics when $ \gamma $ grows from 0.2 to 1.0. However, the decline is less severe for FunClu and FunClu_Smoothing compared with GMM and KMeans, confirming the robustness of our proposed functional clustering framework.

Finally, comparing across correlation strengths ($ v=0.9 $ vs $ \mathrm{v}=0.6 $), we observe that stronger correlation structures generally improve performance across all methods, with a higher CA and ARI at $ \mathrm{v}=0.6 $. Nevertheless, the relative advantage of FunClu_Smoothing remains consistent, underscoring its robustness to different covariance conditions.

Summary

In summary, the simulation results consistently demonstrate the advantages of the proposed FunClu framework. By jointly modeling the cluster-specific mean structure through an allometric scaling law and incorporating SAD1 covariance, FunClu achieves more accurate clustering performance than standard approaches in most parameter settings. However, in several high-noise configurations, FunClu occasionally suffers from empty clusters, which, although accompanied by higher numerical indices, may reduce the interpretability of the clustering solution. The smoothing variant effectively mitigates this issue by regularizing the mixing proportions through a softmax transformation, thereby preventing class loss and ensuring that all identified clusters remain biologically interpretable. Collectively, these findings suggest that FunClu, particularly with smoothing, provides a principled and flexible tool for functional clustering with high-dimensional dependent data, laying a solid foundation for subsequent applications to real biological datasets.

The simulation study is intended to evaluate the proposed method under its target modeling assumptions. Testing robustness to severe deviations from the allometric mean structure is beyond the scope of the present work and will be investigated in future studies.

To further demonstrate the practical utility of our proposed method, we applied it to a transcriptomic dataset of Prunus mume 'Meiren', which was previously reported by Meng et al^[34]. The dataset consists of RNA-seq measurements from three distinct tissues (leaves, exocarps, and mesocarps), providing a comprehensive view of tissue-specific gene expression patterns^[34]. This setting is particularly suitable for applying our multivariate functional clustering model, as it allows us to simultaneously capture both shared and tissue-specific gene expression patterns.

The number of clusters was determined by the BIC. As shown in Fig. 1, the BIC curve reached its minimum at K = 72, which was selected as the optimal number of clusters for subsequent analysis.

Figure 1. 

BIC values under different cluster numbers, with the optimal value at K = 72.

On the basis of this choice, the trivariate clustering identified 72 gene modules, each representing a distinct expression pattern across the three tissues. Several clusters exhibited uniform expression across leaves, fruit peel, and fruit flesh, indicating their involvement in core biological processes shared across tissues. In contrast, many clusters demonstrated tissue-specific expression.

The resulting cluster assignments are illustrated in Fig. 2, where the genes are grouped by their expression profiles across the three tissues.

Figure 2. 

Trivariate clustering results showing 72 gene modules across leaf, fruit peel, and fruit flesh tissues.

To further extract biologically relevant modules, we constructed a gene interaction network using the idopNetwork^[35] algorithm (results not shown). On the basis of network topology and biological interpretability, we selected three hub modules (M14, M68, and M71) for detailed functional enrichment analysis, and the results are shown in Fig. 3. Gene Ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were performed using the OmicShare tools^[36−38].

Figure 3. 

GO and KEGG enrichment analysis of selected core modules. (a) GO and KEGG enrichment results for modules M14, M68, and M71. (b−d) Classification of enriched KEGG pathways at Level 2 for modules M14, M68, and M71.

The enrichment patterns observed in the three core modules suggest both convergent and divergent regulatory strategies across Prunus mume 'Meiren' tissues.

Both M14 and M68 are enriched in DNA replication and repair pathways, which are fundamental to cell proliferation and genome stability. The concurrent enrichment in protein synthesis/processing (M14) and protein transport/localization (M68) suggests that these two modules capture different stages of cellular maintenance.

Although these modules are less active in the fruit flesh than in the leaves and peel, they collectively highlight their complementary yet distinct regulatory roles in maintaining tissue function.This observation may reflect tissue-specific downregulation of cell-cycle and protein-turnover processes in fruit flesh during its developmental stage.

Unlike M14 and M68, Module M71 highlights secondary metabolism, specifically phenylalanine biosynthesis. This pathway serves as the metabolic entry point for anthocyanin biosynthesis, a class of pigments responsible for purple and red coloration.

The enrichment of this pathway aligns with the distinct pigmentation phenotype of Prunus mume 'Meiren', in which anthocyanins are distributed throughout the vegetative and reproductive tissues.

The finding that this pathway shows reduced expression in fruit flesh relative to leaf and peel suggests the possible spatial regulation of anthocyanin accumulation, with peel and leaf being the predominant pigment-producing tissues. This agrees with the common biological observation that fruit peel, rather than flesh, contributes more strongly to the intensity of pigmentation in many species.

Together, these modules reflect an interplay between primary processes (DNA/protein metabolism) and secondary metabolism (phenylpropanoid pathway).

The simultaneous detection of DNA maintenance pathways and phenylalanine metabolism may indicate that clustering does not only capture the similarity of expression but also functionally meaningful co-regulation. For example, plant cells undergoing stress or differentiation may coordinately downregulate DNA/protein biosynthesis while upregulating pathways that are critical for phenotypic traits such as pigmentation.

In summary, application of the proposed trivariate functional clustering to Prunus mume transcriptome data successfully identified 72 expression clusters, among which the core modules demonstrated clear functional relevance to tissue differentiation and pigmentation. The integration of clustering with network-based module selection and enrichment analysis underscores the utility of the method in uncovering biologically meaningful regulatory pathways in complex plants.

Wu's FunClu is a powerful statistical algorithm for clustering high-dimensional dynamic data on the basis of their underlying trajectories rather than static data points. In this study, we proposed an improved framework of Wu's FunClu that integrates an allometric scaling-based mean structure, a novel covariance modeling strategy, and a hybrid EM optimization procedure. This differs from conventional functional clustering approaches that rely on periodic Fourier basis functions, which restrict analyses to time-series data. The improved FunClu's incorporation of allometric scaling extends its applicability to static gene expression profiles. Moreover, the improved numerical optimization strategy significantly accelerates the convergence speed compared with conventional algorithms. Through extensive simulations, the improved version of Wu's FunClu demonstrates superior accuracy and robustness compared with classical clustering benchmarks, including Kmeans and GMM. Importantly, the use of softmax-based smoothing effectively prevented class collapse, securing robust performance in high-dimensional biological contexts.

Application to real transcriptomic data from Prunus mume revealed biologically meaningful modules across leaf, fruit peel, and fruit flesh tissues^[34]. Notably, our model identified core gene clusters enriched in fundamental processes such as DNA replication, repair, and protein metabolism, as well as tissue-specific pathways like phenylalanine biosynthesis, which is directly linked to anthocyanin production and the purple phenotype of ornamental cultivars. These findings highlight both the statistical reliability and the biological interpretability of our approach.

Although the proposed framework is built upon an allometric assumption, we note that such a relationship is empirically observed in the majority of omics datasets encountered in our applications. This reflects a common part–whole scaling behavior arising from coordinated biological regulation. Nevertheless, the method is intended for use in conjunction with an initial diagnostic assessment, and its application to data with pronounced nonstationarity or abrupt structural changes warrants further methodological extensions.

Overall, the proposed method provides a versatile statistical tool for clustering high-dimensional omics data, offering both methodological innovation and practical utility in functional genomics and beyond.