Functional regression clustering

Next-generation sequencing technology (specifically RNA-sequencing or RNA-seq) allows researchers to accurately measure gene expression for all genes in a biological sample (Lister et al., 2008). Until recently, it was prohibitively expensive to perform RNA-seq experiments at more than a few time points at once. RNA-seq is now widespread and affordable enough to use it to investigate time-sensitive biological processes, such as response to environmental stimuli or the organism’s internal clock (Swift & Coruzzi, 2017). In this context, typical biological questions are to detect genes that are differentially expressed over time, and clustering genes according to their expression time courses. Such clustering efforts have mainly been focussed on finding subgroups of genes sharing common time course patterns (Ezer, Shepherd, Brestovitsky, Dickinson, Cortijo, Charoensawan, Box, Biswas, Jaeger & Wigge, 2017; Ezer, Jung, Lan, Biswas, Gregoire, Box, Charoensawan, Cortijo, Lai, Stöckle, Zubieta, Jaeger & Wigge, 2017).

In addition to single RNA-seq time-course experiments, it is now common to perform multiple time series RNA-seq experiments under multiple different treatments (Ezer, Shepherd, Brestovitsky, Dickinson, Cortijo, Charoensawan, Box, Biswas, Jaeger & Wigge, 2017; Ezer, Jung, Lan, Biswas, Gregoire, Box, Charoensawan, Cortijo, Lai, Stöckle, Zubieta, Jaeger & Wigge, 2017), and methods for analysing such data are not well-established. In this paper, we propose a fundamentally different approach to clustering such multiple gene expression time course, by using the link between them—through a functional regression—as a way of clustering genes. This strategy would be able to group together genes that may have very different temporal expression profiles in different conditions, but whose profiles change in the same way across treatments. We hypothesize that such genes would be part of the same pathways. For instance, two genes might have different gene expression patterns under normal growth conditions, because they are normally regulated by different sets of regulatory proteins. However, both these genes might be regulated by the same stress-associated regulatory protein in response to an environmental stimulus, which causes their gene expression pattern to be perturbed in the same way.

Our approach is based on treating gene expression time courses as curves that are sampled discretely, with measurement error, and falls therefore within the realm of functional data analysis (FDA: Ramsay & Silverman, 2005; Ferraty & Vieu, 2006; Horváth & Kokoszka, 2012; Hsing & Eubank, 2015; Kokoszka & Reimherr, 2017), now a prevalent area of statistics, with applications in numerous fields, such as in neuroimaging (Jiang et al., 2009; Zipunnikov et al., 2011; Ivanescu et al., 2015; Li et al., 2015; Jiang et al., 2016; Petersen et al., 2016; Rügamer et al., 2018; Palma et al., 2020), phonetics (Aston et al., 2010; Hadjipantelis et al., 2012, 2015; Pigoli et al., 2018; Tavakoli et al., 2019), or genomics (Serban & Wasserman, 2005; Yao et al., 2005; Reimherr et al., 2014; Ezer & Keir, 2019). Classical statistical modelling tools have been extended to functional data: the functional linear model (FLM) has received a lot of attention (see the review of Morris, 2015, and references therein), and many generalizations have been proposed, such as those inspired by generalized linear models (Müller & Stadtmüller, 2005), (generalized) additive models (Müller & Yao, 2008; McLean et al., 2014; Scheipl et al., 2016), or non-parametric regression (Ferraty & Vieu, 2006).

Parallel to this, there have been many contributions to the clustering of functional data: generalizations of $K$-means have been proposed, usually after projecting the functional observations onto some finite basis of functions (Abraham et al., 2003; Serban & Wasserman, 2005; Auder & Fischer, 2012; Antoniadis et al., 2013), or onto the first functional principal components (FPC; Peng et al., 2008). Extensions of such functional $K$-means have also been proposed: Chiou & Li (2007) proposed using the subspaces spanned by FPCs as representative of the clusters, instead of using cluster means to define the clusters, and Gattone & Rocci (2012) proposed a functional version of the reduced $K$-means algorithm (De Soete & Carroll, 1994). Delaigle et al. (2019) proposed to choose adaptively the projections onto which the $K$-means algorithm is applied, and showed that the technique can yield asymptotically perfect clustering. Mixture models for functional clustering have also been proposed (Chudova et al., 2004; Petrone et al., 2009; Schmutz et al., 2020).

Combining functional linear models and functional clustering has, however, been far less studied. Such a problem is motivated, for instance, by trying to find subgroups in the data characterized by different relationships between the (scalar or functional) response, and one or several functional covariates. An example of such problem arises in plant genomics, where one is interested in clustering the genes of a plant based on the relationship of its circadian expression in summer, say $Y(t)$, where $t\in[0,48]$ is a time observation in hours i.e. from a process observed during two entire days (and 0 represents midnight of the first day), and its relationship with the circadian expressions in autumn, winter and spring, say $X_{1}(t),X_{2}(t),X_{3}(t),t\in[0,48]$. In this setting, one could consider a cluster-specific functional linear model, such as

if gene $i$ belongs to group $k\in\{1,\ldots,K\}$, where $\operatorname{\mathbb{E}}\varepsilon_{i}=0$. The group memberships are of course unknown, but finding clusters based on model (1) is of interest, and gives promising results when applied to gene expressions timecourses of Arabidopsis halleri specimen, see Section 5.1.

To the best of our knowledge, only a handful of papers have considered clustering functional data using functional models similar to (1) Yao et al. (2011) looked at scalar-on-function regression modelling (i.e. the case $Y(t)=Y\in\mathbbm{R}$ in (1)) with one functional covariate, and used FPCA for regularization. Wang et al. (2016) considered the concurrent functional linear model $Y(t)=\sum_{j=1}^{p}\beta_{jk}(t)X_{j}(t)+\varepsilon_{k}(t)$ if the observation comes from group $k$, and used a mixture of Gaussian processes to fit the model using an EM-type algorithm. In this paper, we present a multivariate functional clustering method based on a cluster-specific functional linear model like (1), called Functional Regression Clustering (FRECL). It clusters data of the form $\{(Y_{i},X_{i1},\ldots,X_{ip}):i=1,\ldots,n\}$, where $Y_{i},X_{i1},\ldots,X_{ip}$ are curves, and allowing $p>1$. Our proposed method has been developed with the application to gene expression time courses in mind, but it can also be used in other applications, such as the multiple sclerosis applications of Ivanescu et al. (2015). Indeed, FRECL is versatile enough to be easily applied to cluster any data set in which multiple functional data sets are being compared.

The paper is organized as follows. In Section 2, we describe and present our FRECL model and method for finding the clusters and the regression surfaces $\beta_{jk}(t,s)$. In Section 3, we provide a description of our motivating application for FRECL. In Section 4, we provide an extensive simulation study of the performance of our method in data that resembles that in our motivating example, and compare it to existing (functional) clustering methods. In Section 5, we identify clusters in an expression time course data set utilising FRECL. Section 6 concludes with a discussion. The full code of our method is made available for reproducing our results, and for further applications, see https://gitlab.com/sconde778/frmm_rpackagefunctions.git

Given multivariate functional observations $\{(Y_{i},X_{i1},\ldots,X_{ip}):i=1,\ldots,m\}$, where $Y_{i},X_{i1},\ldots,X_{ip}\in{L^{2}(\mathcal{T})}$, where ${L^{2}(\mathcal{T})}=\{f:\mathcal{T}\to\mathbbm{R}:\lVert f\rVert<\infty\}$ is the Hilbert space of square integrable functions defined over the compact interval $\mathcal{T}\subset\mathbbm{R}$, with norm $\lVert f\rVert=(\int_{\mathcal{T}}f^{2}(t)dt)^{1/2}$, the goal of our paper is to cluster these observations into groups according to the relation between the responses $Y_{i}\in{L^{2}(\mathcal{T})}$, and the predictors $X_{i1},\ldots,X_{ip}\in{L^{2}(\mathcal{T})}$. While we assume that the functional responses and predictors are all defined on the same domain $\mathcal{T}$, our methods can easily be extended to settings with distinct domains for the response and each predictors.

Let $\mathcal{P}_{mK}$ denote the set of partitions of $\{1,\ldots,m\}$ into $K>1$ disjoint sets (which we shall call clusters), with elements of the form $P=\{C_{1},\ldots,C_{K}\}$. Notice that we allow partitions with empty clusters, which implies that $\mathcal{P}_{m,K-1}\subset\mathcal{P}_{m,K}$ provided we identify sets of the form $\{A_{1},\ldots,A_{L},\emptyset\}$ with $\{A_{1},\ldots,A_{L}\}$, where $A_{1},\ldots,A_{L}$ are nonempty sets, $L\geq 1$. Let $\mathbbm{1}_{A}$ the indicator function of the set $A$, defined by $\mathbbm{1}_{A}(x)=1$ if $x\in A$, and $\mathbbm{1}_{A}(x)=0$ otherwise. We assume that the observations $\{(Y_{i},X_{i1},\ldots,X_{ip})\}$ come from the following model,

where $P^{*}=\{C_{1}^{*},\ldots,C_{K}^{*}\}\in\mathcal{P}_{mK}$ is a fixed partition, $\varepsilon_{i}(t)$ is a functional error term, $\beta_{0k}\in{L^{2}(\mathcal{T})}$, and $\beta_{jk}\in L^{2}(\mathcal{T}\times\mathcal{T})$, $j=1,\ldots,p$, $k=1,\ldots,K$. In other words, the same functional linear model links $Y_{i}$ and $(X_{i1},\ldots,X_{ip})$ within each cluster $i\in C^{*}_{k}$, but the functional parameters are (possibly) distinct across the clusters $C^{*}_{1},\ldots,C^{*}_{K}$. The goal is to find the unknown partition $P^{*}$.

Letting ${\boldsymbol{\beta}}_{k}=(\beta_{0k},\beta_{1k},\ldots,\beta_{pk})\in{L^{2}(\mathcal{T})}\times L^{2}(\mathcal{T}\times\mathcal{T})\times\cdots\times L^{2}(\mathcal{T}\times\mathcal{T})$ be a functional vector, and defining

then (2) can be rewritten as

Another way of expressing model (2) is to say that the model

holds within each cluster of $P^{*}$, with possibly distinct functional parameters ${\boldsymbol{\beta}}$. Note that it is not straightforward to view model (2) as a mixture model since density functions are generally not well defined in a functional context (Delaigle et al., 2010). Notice that (4) is a function on function regression model with multiple functional predictors, and that $\Phi((X_{i1},\ldots,X_{ip}),{\boldsymbol{\beta}})$ can be viewed as the conditional expectation of the functional response $Y$ given $(X_{1},\ldots,X_{p})$ and ${\boldsymbol{\beta}}$, provided $\operatorname{\mathbb{E}}\varepsilon(t)=0,t\in\mathcal{T}$.

Our strategy for finding $P^{*}$ is inspired by the $K$-means algorithm (MacQueen et al., 1967), which we can view as an iterative procedure alternating between a model fitting (the computation of the cluster means given the cluster allocations) and a partition update (assigning each observation in the next iteration to the current closest cluster mean). We therefore propose to estimate $P^{*}$ using an iterative procedure, summarized as follows:

1. 1.

   Pick $P\in\mathcal{P}_{mK}$ with non-empty clusters at random,

2. 2.

   Fit model (4) within each cluster of $P$, thus obtaining $K$ fitted models,

3. 3.

   Reassign the $i$th observation, $i=1,\ldots,m$ to the best fitting model $\hat{k}(i)\in\{1,\ldots,K\}$, thus defining a new partition $P^{+}$,

4. 4.

   Set $P\leftarrow P^{+}$ and repeat steps 2–3 until convergence.

Notice that step 2 requires a fitting method $\mathcal{F}$. We discuss the choice of $\mathcal{F}$ below. Step 2 results in $K$ fitted models of the form (4), with estimates $\hat{\boldsymbol{\beta}}_{k}$. Step 3 requires finding the best fitting model for each observation $(Y_{i},X_{i1},\ldots,X_{ip})$; we choose to do this by computing the norm of its fitted residuals under each model,

Other methods for updating clusters can be chosen, such as choosing a different norm in (5). The full version of our generic FRECL algorithm, using the fitting method $\mathcal{F}$, is given in Algorithm 1.

Because the output of Algorithm 1 depends on the initial random partition, we propose to use consensus clustering (Monti et al., 2003) to produce more consistent results. Consensus clustering consists of running a clustering algorithm multiple times, with different initial partitions, and then aggregating the obtained clusters. We describe it formally in Algorithm 2.

Algorithms 1 and 2 depend on a number of parameters, which we now briefly discuss and make recommendations about. A more detailed discussion of the choice of some of these parameters is deferred to Section 4.

Whilst Schmutz et al. (2018) claim that FunHDDC, a ‘specific’ method for clustering multivariate functional data, also works for univariate functional data, our simulation results show that it is outperformed by the other approaches we consider; even by those developed for non-functional, high dimensional variables such as HDDC. In one of the illustrations included in Bouveyron & Jacques (2011), FunHDDC is outperformed by HDDC, too.

Functional mixture models make possible to study relationships between explanatory and response variables over time allowing for clusters characterized by these relationships. Chamroukhi (2015) presents a clustering method for mixture regression that involves a penalised likelihood, where the penalty is the total entropy. FRECL is an alternative to these settings that does not need to consider a constrained optimization problem, and consequently does not need to estimate the value of the Lagrange multiplier. Yao et al. (2011) propose a functional mixture model implemented with FPCs in order to overcome overfitting with a finite number of observations and infinite-dimensional parameters and as usual selecting a certain number of components. The FPCs are necessarily computed (only) with the explanatory variables. Since the estimation of the slope parameters involves only the same number of selected FPCs, as well, the slope parameter space is restricted, and thus their estimates may be far from the truth. FRECL, clustering based on the functional ‘mixture’ model (2), improves existing methods, although these were not developed specifically for generating processes with an underlying functional mixture model. Consequently, it is a step further in the development of functional data clustering.

Software in the form of R code, together with a sample input data set and complete documentation is available at https://gitlab.com/sconde778/frmm_rpackagefunctions.git.

Conflict of Interest: None declared. The authors thank Ioannis Kosmidis for suggesting appropriate venues for this manuscript. This project was funded by the Alan Turing Institute Research Fellowship under EPSRC Research grant (TU/A/000017) to DE; EPSRC/BBSRC Innovation Fellowship (EP/S001360/1) to DE and SC. ST would like to thank the Isaac Newton Institute for Mathematical Sciences, Cambridge, for support and hospitality during the programme Statistical Scalability where work on this paper was undertaken. This work was supported by EPSRC grant no EP/R014604/1.