COREclust: a new package for a robust and scalable analysis of complex data

We consider the set of $n$ observations $I=\{X_{1},\cdots,X_{n}\}$ of the studied system. Each observation is described by $p$ quantitative variables, so $X_{i}=(x_{i}^{1},\cdots,x_{i}^{p})$ for $i\in{1,\cdots,n}$. The representative variables in the system will be identified based on a quantitative description of the variable interactions. This requires to define a notion of distance between the observed variables. In this paper, we use standard Pearson correlations but other distance measures could be alternatively used. The interactions are then modelled as an undirected weigted graph where the nodes represent the variables and the edges represent the variable pairs having a non-negligible similarity score. Each edge has also a weight equals to the similarity score between the nodes/variables it connects.

We denote $G(N,E)$ the undirected weigthed graph modeling the observed system. The set $N$ contains $p$ nodes, each of them representing an observed variable, and the set $E$ contains the weighted edges that model the similarity level of the variable pairs.

As mentioned section 1, the representative variables of the studied system are indirectely detected as the centers of CORE-clusters. In our work, this central variable is simply the one that has the lowest average distance to all other variables in the CORE-cluster. The central part of our strategy is then the detection of CORE-clusters of variables based on the graph $G(N,E)$. CORE-clusters are groups of nodes/variables that are densely connected to each other. The notion of $k$-core introduced firstly by Seidman  [20] for unweighted graphs in 1983, formally describe those cohesive groups. To this end, let us consider the degree $deg(n)$ of a node $n\in N$ as the number of edges $e\in E$ incident to itself.

In this paper, we developed the CORE-clustering algorithm based on the notion of CORE-clusters. In order to prevent instability in the process of core detection from one run to the other, the user has to tune the coefficient $\tau$ so that the size of the CORE-clusters to be detected will be in $\left\{\tau,2\tau-1\right\}$. Thus, the $k$-cores to be detected directly depends on this coefficient. The latter controls the level of regularisation when estimating the representative variables: Large values of $\tau$ leads to the detections of large sets of strongly related variables. In that case, the representative variables of these sets is likely to be meaningful even if $n<p$. However each CORE-cluster may also contain several variables that would have ideally been representative. On the contrary, too small values of $\tau$ will be likely to detect all meaningful representative variables but also a non-negligible number of false positive representative variables, especially if the observations are noisy or if $n<p$. Note that another optional parameter is the threshold on the minimal similarity observed between two variable to define an edge in $G(N,E)$. Using it with values stricly higher than $0$ allows to reduce the graph size.

An important characteristic of our algorithm is that it first strongly simplifies $G(N,E)$ using a maximum spanning tree strategy. The resulting graph $G(N,T)$ has the same nodes as $G(N,E)$ and the tree structure $T$ instead of the initial edges $E$, where $T$ is a subset of $E$. All node pairs are then still connected to each other but there is only a single path between each node pair. This strongly sparsifies the graph and makes possible the use of fast and scalable algorithms as described section 3. CORE-clusters detection strategy, which estimates strongly interconnected variable sets also makes sense on $G(N,T)$ as computing a maximum spanning tree aims to find a tree sturcture out of the edges $E$ with the highest possible weights in average.

The very first step of our strategy is to quantify the similarity between the different observed variables. The similarity is first computed and represented as a dense graph $G(N,E)$, where $N$ contains $p$ nodes, each of them representing an observed variable, and $E$ contains $K_{E}=p(p-1)/2$ undirected weighted edges that model a similarity level between all pairs of variables. Various association measures can be used to quantify the similarity between variables $i$ and $j$, $1\leq i,j\leq p$. The most popular ones are Spearman or Pearson correlation measures, different notions of distances, or smooth correlation coefficients. We used the absolute value of Pearson’s correlation coefficient in our package, so the weight between the nodes $i$ and $j$ is the Pearson’s correlation coefficient between the absolute values of the observed variables $i$ and $j$. In this approach, we consider the variables with no connection in between which will have a correlation coefficient set to zero. The algorithmic cost of this estimation using a naive technique is $\mathcal{O}(n^{2}p)$ but it can be easily parallelized using divide and conquer algorithms for reasonably large datasets, as in [18]. For large to very large datasets, correlations should be computed on sparse matrices, using e.g. [7] to make this task computationally tractable.

Once $G(N,E)$ is computed, its complexity is strongly reduced to make possible the efficient use of scalable clustering strategies on our data.

The maximum spanning tree of $G(N,E)$ is built by using Kruskal’s algorithm [11]. We denote $G(N,T)$ the resulting graph, where $T$ has a tree-like structure. Details of the algorithm are given in Alg. 1. The algorithmic cost of the sort procedure (line 1) is $\mathcal{O}(K_{E}\log{(K_{E})})$. Then, the for loop (lines 4 to 10) only scans one time the edges. In this for loop, the most demanding procedure is the propagation of label $L(N_{i})$, line 8. Fortunately, the nodes on which the labels are propagated are connex to $N_{j}$ in $G(N,E)$. We can then use a Depth-First Search algorithm [23] for that task, making the average performance of the for loop $\mathcal{O}(K_{E}\log{(p)})$.

In contrast with the other clustering techniques, our CORE-clustering approach detects clusters having an explicitely controled granularity level, and only gathers nodes/variables with a behaviour considered as similar (see section 3.1). Core structures are detected from the maximal spanning tree $G(N,T)$ by gathering iteratively its nodes $N$ in an order that depends on the edge weights in $T$. CORE-structures are identified when a group of gathered nodes has a size in $\{\tau,\cdots,2\tau-1\}$, where $\tau$ is the parameter that controls the granularity level. Algorithm details are given Alg. 2.

To understand how Alg. 2 detects clusters of coherent nodes, we emphasize that the CORE-clustering algorithm is applied on the maximum spanning tree $G(N,T)$ and not on $G(N,E)$. Only the most pertinent relations of each node $N_{n}$ are then taken into account, so the nodes with no pertinent relation will be randomly linked to each other in $T$. As Alg. 2 will first gather these nodes, they will generate many small and meaningless clusters. Parameter $\tau$ should then not be too small to avoid considering these clusters as CORE-clusters. Then, the first pertinent clusters will be established using the edges of $T$ with larger weights, leading to pertinent node groups. Finally, the largest weights of $T$ are treated in the end of Alg. 2 in order to split into several CORE-clusters the nodes related to the most influent nodes/variables.

Remark that the algorithmic structures of Alg. 1 and Alg. 2 are similar. However Alg. 2 runs on $G(N,T)$ and not on $G(N,E)$. The number of edges $K_{T}$ in $G(N,T)$ is much lower than $K_{E}$ since $T$ has a tree structure and not a complete graph structure. It should indeed be slightly higher than $p$ [22], which is much lower than $K_{E}=p(p-1)/2$. Moreover the propagation algorithm (lines 9 and 11) will never propagate then labels on more than $2\tau-1$ nodes. The algorithmic cost of the sort procedure (line 1) is then $\mathcal{O}(K_{T}\log{(K_{T})})$ and the average performance of the for loop $\mathcal{O}(K_{T}\log{(\tau)})$.

We propose an alternative strategy for the CORE-clusters detection algorithms: The edge treatment queue may be ordered by following decreasing edge weights instead of increasing edge weights. The nearest edges will then be first gathered making coherent CORE-clusters as in Alg. 2, although one CORE-cluster may contain several key players. To avoid gathering noisy information, the for loop on the edges (line 6 of Alg. 2) should also stop before meaningless edges are treated. This strategy has a key interest: It can strongly reduce the computational time dedicated to Algs. 1 and 2. By doing so, Alg. 1 and modified Alg. 2 are purely equivalent to Alg. 3.

Again, the structure of this algorithm is very similar to the structure of the Maximum Spanning Tree strategy section 3.1. The algorithmic cost of the sort procedure (line 1) is $\mathcal{O}(K_{E}\log{(K_{E})})$. Then, the loop lines 4 to 13 scans $floor(\gamma\tau)$ edges where $\gamma$ should be sufficiently high to capture the main CORE-clusters but also where $\gamma\tau<<K_{E}$. Labels propagations in this loop (lines 7 and 9) is also limited to $2\tau-1$ nodes. The average performance of the for loop is therefore $\mathcal{O}(\gamma\tau\log{(\tau)})$.

Once a CORE-cluster is identified, we use a straighforward strategy to select its central variable: the distance between all pairs of variables in each CORE-cluster is computed using a Dijkstra’s algorithm [6], [28] in $G(N,E)$. The central variable is then the one that has the lowest average distance to all other variables in the CORE-cluster. As the CORE-clusters have less than $2\tau$ nodes, the algorithmic cost of this procedure is $\mathcal{O}(\tau^{2})$ times the number of detected CORE-clusters, which should remain low even for large datasets.

The main C++ function of COREclust called by R is principalFunction. Its first inputs are the name of the file containing the similarity level between the variables (arg_InputFile), the minimal size of the cores (arg_MinCoreSize) and the coefficient $\gamma$ (arg_Gamma). Note that $\gamma$ is only used when performing the greedy algorithm. The next input is the name of the file containing the output variable labels (arg_OutputFile). The last one is finally a coefficient $1$ or $2$ that controls whether the standard or the greedy algorithm will be performed (greedyOrNot). It is worth mentionning that the input csv file is a correlation matrix generated by R before executing the C++ code. Then the structure of principalFunction is as described in Alg. 4.

Remark that the C++ code heavily depends on the class CompactGraph, that is summarized in file CoreClust.h. This class contains all variables to store undirected weighted graphs and node Labels. It also contains different methods (i.e. class functions) to read undirected weighted graphs in a list of weighted edges or a squared matrix, to write a graph and its labels, as well as to perform the CORE-clustering strategies described in Alg. 2 and Alg. 3. As mentioned section 2 the scalability of our algorithms directly depend on the algorithm that sorts the graph edges depending of their weights. To make this step computationally efficient, we have used in the SortByIncreasingWeight and SortByDecreasingWeight methods of CompactGraph the reference C++ function std::sort. For $N$ elements, this function performs $\mathcal{O}(N\log{(N)})$ comparisons in average until C++11, and $\mathcal{O}(N\log{(N)})$ comparisons after C++11.

User interactions are made through a R script algo_clust.R that executes the C++ code of Alg. 4. It first takes as inputs the name of a Rdata or a csv file (arg_InputFile). The remaining arguments are the minimal size of the cores (arg_MinCoreSize), the name of the file containing the output variable labels (arg_OutputFile) and optionally the $\gamma$ coefficient (arg_Gamma) that controls the edge number to be scanned to capture the main cores. If it is not called as argument, $\gamma$ is set to zero and the fifth input argument of the principalFunction is set to 1 (CORE-clustering algorithm will be applied on submit data). Otherwise, greedyOrNot is set to 2.

Note that the script algo_clust.R contains a R function corr_Mat that computes the similarity matrix of the observed data. It takes as arguments the name of the file containing the dataset Obs_file_name and a coefficient FileType equals to $1$ if Obs_file_name is a csv file and $2$ otherwise. In the first case, the absolute value of Pearson’s correlation is computed using the function cor of package stat ([24]) and the associated adjacency matrix is built. In the second case, the input matrix that contains the observations is of type sparseMatrix (out of the package Matrix [13]). Then, the absolute value of Pearson’s correlation is computed on its coefficients using the function corSparse of package qlcMatrix [7].

It also includes the function Rcpp::sourceCpp from the package Rcpp that parses the specified C++ file or source code. A shared library is then built and its exported functions and Rcpp modules are made available in the specified environment.

In the R file, the similarity matrix is first computed and saved in the csv file corrMat.csv that will be loaded by the principalFunction of the C++ code.

Many complex networks are known to be scale free, meaning that the network has a little amount of highly connected nodes (hubs) and many poorly connected nodes. Simulating such networks can be seen as a process through which we first generate the profile of leader actors (hubs) and then the profile of the remaining ones around those hubs.

We consider here that the simulated individual profiles are divided in $K$ different clusters of size $n_{C_{1}}$,$n_{C_{2}}$,…,$n_{C_{K}}$. A simulated expression data set $X\in\mathbb{R}^{N\times p}$ is then composed of $p=n_{C1}+n_{C_{2}}+...+n_{C_{K}}$ actor profiles generated for the $K$ clusters and has the structure:

Strategy for simulating actor profiles in a single cluster of size $n_{C}$ is:

1. 1.

   Generate the profile $\mathbf{x}^{(1)}=(x_{1}^{(1)},x_{2}^{(1)},...,x_{N}^{(1)})^{\prime}$ of one variable from a normal distribution, $x_{i}^{(1)}\sim\mathcal{N}(0,1),\forall i=1,2,...,N$. This profile defines the leader or hub profile.

2. 2.

   Choose a minimum correlation $r_{min}$ and a maximum correlation $r_{max}$ between the leader actor and other actors in the cluster.

3. 3.

   Generate $n_{C}-1$ profiles such that the correlation of the j-th profile $\mathbf{x}^{(j)}$ with the leader actor profile $\mathbf{x}^{(1)}$ is forall $j=2,3,...,n_{C}$ close to

   $$r_{j}=r_{min}+\Delta_{r}\left(1-\frac{j}{n_{C}}\right),$$

   where $\Delta_{r}=r_{max}-r_{min}$. The $r_{2}$ correlation is close to $r_{max}$ and the $n_{C}$-th correlation is equal to $r_{min}+\Delta_{r}=r_{max}$. The j-th profile is also generated by adding a Gaussian noise term to the leader profile, to allow a correlation with the leader profile close to the required correlation $r_{j}$

   $$x_{i}^{(j)}=x_{i}^{(1)}+\sqrt{\left(\frac{1}{r_{j}^{2}}-1\right)}\epsilon_{i}^{(j)},\forall i=1,2,...,N$$

   where $\epsilon_{i}^{(j)}\sim\mathcal{N}(0,\alpha)$.

From the matrix including the profile of each actor of the system $S$, the coefficients $s_{ij}\in[0,1]$ of the $p\times p$ matrix of Pearson correlation will be computed and finally represent the edges $E$ of the simulated network $G(N,E)$.

Clustering results can be validated using two scores: The external validation compares the clustering result with a ground truth result. Internal validation also evaluates the clustering quality using the result itself plus the graph properties, and no ground truth data. It can therefore also be used on real-life data. The latter will be introduced in 6.

Standard and greedy algorithms were applied on four different simulations repeated $30$ times, to study their performances and robustness. On the basis of the built dataset in section 4.1, each simulation was generated with well-defined clusters inside which correlations between actor profiles decrease linearly from $r_{max}=1$ to $r_{min}=0.5$. Let simulation be the name of the function modeling the dataset introduced above.

1. a)

   The first simulation matches the generation of the standard dataset of size $100\times p$ defined in section 4.1 with $K=2$ and $K=5$ clusters of size $(2,2)$ and $(7,\cdots,7)$;

2. b)

   The second one consisted in creating noisy data and testing the CoreClustering algorithms on it. To this end, we simulated the standard dataset, and gradually added to each variable an increasing gaussian noise term drawn from a normal distribution $\mathcal{N}(0,\alpha)$ depending on a coefficient $\alpha=\left\{0.25,1.5\right\}$. Examples of simulated matrices with this noise term are shown in Fig. 1;

3. c)

   The third simulation aims at changing the minimal size of core structures detected i.e the level of granularity by the algorithms along the process taking values from $(20,30,40,50)$ and checking their accuracy. The number of clusters was fixed at $K=5$ with sizes randomly drawn from the set $\left\{50,60\right\}$;

4. d)

   The last simulation consists in varying the size of the sample from 5 to 30 and setting $K=5$ and the minimal size of core structures to $20$.

For all the simulated datasets, CoreClustering and CoreClustering_Greedy were applied on their matrices of correlation and then, the cores detected were compared to the real $K$ cluster groups thanks to the quality measurement. The results are shown in Fig. 2 in the form of boxplots. The first four figures (a) and (b) show that the standard CoreClustering is slightely more robust to noisy dataset than the CoreClustering_Greedy . The opposite effect can be observed when the sample size of the dataset decreases ((d)). Moreover, it should be noticed that increasing the number of clusters in the dataset slightly declines the accuracy of the detected cores. However, the two algorithms are really effective when varying the granularity coefficient ((c)) which differentiate it from other clustering techniques mentioned in the introduction. Nevertheless, it must be noted that the CoreClustering_Greedy is really more effective in terms of processing time than the other one.

Using simulated data, we showed that CORE-clustering algorithms performs well to detect core structures in a high dimensional network. This section illustrates the use of Core-Clustering algorithms on the well-known synchronized yeast cell cycle data set of [21], which includes 77 samples under various time during the cell cycle and a total of 6179 genes, of which 1660 genes are retained for this analysis after preprocessing and filtering. The goal of this analysis is then to detect core structures among the correlation patterns in the time series of gene expressions of yeast measured along the cell cycle. From this dataset introduced as a matrix, a measure of similarity between all gene pairs was measured with the absolute value of Pearson’s correlation. A total of about $1.3\times 10^{6}$ weighted edges will then be considered when representing the variable correlations in a graph.

The second application of the two developed algorithms is based on a sparse version of the data yeast mentioned above. From the dataset built, $80\%$ of the information on it is set to zero, making the data sparse. The matrix containing those data was then converted to a sparse matrix type using the package Matrix of R and a measure of similarity suitable for it was quantified using e.g. [7].

Those similarity matrices are then saved in a CSV format and passed as the first argument of the functions algo_clust and algo_clust_greedy. The latter are finally run using the following command line Rscript algo_clust.R in_InputFile in_MinCore_Size in_Outputfile or Rscript algo_clust_greedy.R in_InputFile in_MinCore_Size in_Outputfile, where the in_MinCoreSize and in_Outputfile arguments are fixed by the user.

We are interested here in evaluating the performances of the core detection algorithms using internal criteria. Let us intoduce the following evaluation metrics which will quantify the internal quality, i.e, the minimal property of the density of a core.

Let $W=(w_{ik})_{1\leq i,k\leq p}$ be the correlation matrix between each variable $X^{1},...,X^{p}$. Let $\hat{K}$ the number of cores detected. Each variable $X^{1},\cdots,X^{p}$ is assigned to a unique core $\hat{C_{j}}$ ($j\in[1,\hat{K}]$). The quality measurement of internal connections strength in cores is defined as:

where we suppose that a weight equals to zero means that there is no connection between the two nodes involved. For each node of a given core, this score first determines the maximal strength of the edges that connect this node to its adjacent nodes. Over all those coefficents, the minimal one is then divided by the maximal strength of all edges in the core. As the last coefficient is very closed to one, the final score emphasizes the minimal strength of all vertices in each core between its internal nodes. A score $IC(\hat{K})$ closed to one then means that the algorithm reveals a core with highly connected nodes.

In this context, the standard and greedy Core-Clustering algorithms were applied on the matrices of similarity introduced above with the aim of proving their performances on a real complex dataset. Applied on the original one, they detected respectively $6$ and $25$ core structures whereas $7$ and $37$ were catched when it comes to their sparse version. The minimal size of core structures was set at $30$. A visualization of the core clusters detected by the standard algorithm was made (Figure  3) using Cytoscape²²2http://www.cytoscape.org/ where the transparency on the edges linearly depends on the correlation between the variables. Transparency for correlation 0.5 is $1\%$ and transparency for correlation 1 is $20\%$. The quality of the connections strength in core structures detected in each dataset were then internally evaluated using the quality measurement introduced in section 5. The mean of the scores $IC(\hat{K})$ obtained for each cores when the standard algorithm is applied on the first type of data is in $[0.70,0.84]$ whereas it is in $[0.64,0.87]$ for the greedy algorithm. Thus, the assumption that CORE-clusters are sets of highly connected variables is verified. For the sparse version of the data yeast, the scores $IC(\hat{K})$ are respectively in $[0.69,0.73]$ and in $[0.52,0.77]$. Those low scores could be dued to the limited number of information in the sparse data. Additionnaly, the core structures detected in each dataset were compared using the quality measurement of clusters introduced in Section 5. The scoring was about $S_{j}=0.93$ for the first type of data and $S_{j}=0.95$ for its sparse version. However, in view of the different number of core structures detected by each algorithm, this measurement actually compares the first cores recognized by its greedy version with all of the standard one. Thus, the resultant scores show that the two algorithms are both robust when they are applied on a complex dataset and performs efficiently even on very sparse dataset. The computational times for the standard and greedy version of the clustering algorithm applied on the original dataset were respectively about $319.60$ and $26.86$ seconds on a Core i5 computer with a RAM of $4.00$Go. Then, when they were tested on a sparse dataset, the computational times were about 255.92 and 31.78 seconds.