Table of Contents

cola: A General Framework for Consensus and Hierarchical Partitioning

Author: Zuguang Gu ( z.gu@dkfz.de )

Date: 2019-09-16

Package version: 1.0.1


Introduction

Subgroup classification is a basic task in high-throughput genomic data analysis, especially for gene expression and DNA methylation data analysis. Mostly, unsupervised clustering methods are applied to predict new subgroups or test the consistency with known annotations.

To test the stability of subgroup classifications, consensus clustering or consensus partitioning is always performed. It clusters repeatedly with a randomly sampled subset of data and summarizes the robustness of the clustering, finally it gives a consensus classification of all samples.

Here we present the cola package which provides a general framework for consensus partitioning. It has following advantages:

  1. It modularizes the consensus partitioning processes that various methods can be easily integrated in different steps of the analysis.
  2. It provides rich visualizations for interpreting the results.
  3. It allows running multiple methods at the same time and provides functionalities to compare results in a straightforward way.
  4. It provides a new method to extract features which are more efficient to separate subgroups.
  5. It allows doing partitioning in a hierarchical way to detect subgroups with relatively smaller difference.
  6. It generates detailed HTML reports for the complete analysis.

Following flowchart lists the general steps of consensus partitioning implemented by cola:

The steps are:

  1. Clean the input matrix. The processings are: adjusting outliers, imputing missing values and removing rows with very small variance. This step is optional.
  2. Extract subset of rows with highest scores. Here “scores” are calculated by a certain method. For gene expression analysis or methylation data analysis, \(n\) rows with highest variance are used in most studies, where the “method”, or let's call it “the top-value method” is the variance (by var() or sd()). Note the choice of “the top-value method” can be general. It can be e.g. MAD (median absolute deviation) or any user-defined method.
  3. Scale the rows in the sub-matrix (e.g. gene expression) or not (e.g. methylation data). This step is optional.
  4. Randomly sample a subset of rows from the sub-matrix with probability \(p\) and perform partitioning on the columns of the matrix by a certain partition method, with trying different numbers of subgroups.
  5. Repeat step 4 several times and collect all the partitions.
  6. Perform consensus partitioning analysis and determine the best number of subgroups which gives the most stable subgrouping.
  7. Apply statistical tests to find rows that show significant difference between the predicted subgroups. E.g. to extract subgroup specific genes.
  8. If rows in the matrix can be associated to genes, downstream analysis such as function enrichment analysis can be performed.

All the steps will be explained in detail in following sections.

First of all, we load the cola package.

library(cola)

The input matrix

In most of the analysis tasks, the input matrix contains values of gene expression or DNA methylation (from methylation array or whole genome bisulfite sequencing). If it is an expression matrix, rows correspond to genes, and if it is a methylation matrix, rows correspond to CpG sites. For both types of matrices, columns should always correspond to samples where subgroups are detected. More general, the input matrix can be any type of measurements as long as they represent in a form of matrix, e.g. a matrix where rows are genomic regions and values are the histone modification intensities in the regions, measured from a ChIP-Seq experiment.

In following part of this vignette, we always call the matrix columns as “samples”.

Before doing consensus partition, a simple but important step is to clean the input matrix:

# code is only for demonstration
mat = adjust_matrix(mat)

adjust_matrix() does following preprocessings:

  1. Rows where more than 25% of the samples having NA values are removed;
  2. Use impute::impute.knn() to impute missing data if there is any;
  3. For each row in the matrix, it uses adjust_outlier() (also provided by cola package) to adjust outliers. Values larger than the 95th percentile or less than the 5th percentile are replaced by corresponding percentiles.
  4. Rows with zero variance are removed.
  5. Rows with variance less than 5th percentile of all row variance (after removing rows with zero variance) are removed.

Some of the above steps are optional. For example, methylation matrix does not need to be adjusted for outliers because all the methylation values are already in a fixed data scale (0 ~ 1).

Basic usage

consensus_partition() performs consensus partitioning for a single top-value method and a single partition method. The major arguments for consensus_partition() are:

# code is only for demonstration
res = consensus_partition(mat,
    top_value_method = "MAD",
    top_n = c(1000, 2000, 3000, 4000, 5000),
    partition_method = "kmeans",
    max_k = 6,
    p_sampling = 0.8,
    partition_repeat = 50,
    anno = NULL)

Other arguments can be found in the on-line documentation of consensus_partition().

To get a robust result from consensus partitioning, for a specific top-value method and a specific partition method, consensus_partition() tries different top_n and different number of subgroups. We found that different top_n might give different subgroups for some data sets, thus consensus_partition() pools results from different top_n and aims to give a general consensus subgrouping across different top_n. Later, we will introduce that the bias for top_n on subgrouping can be visualized by membership_heatmap() function.

In most cases, users might not be very sure which top-value method and which partition method are best for their dataset. The helper function run_all_consensus_partition_methods() is a convenient way to try multiple top-value methods and multiple partition methods simultaneously to see which combination of methods gives the best prediction of subgroups.

In run_all_consensus_partition_methods(), most of the arguments are the same as in consensus_partition(), except top_value_method and partition_method now accept a vector of method names and multiple cores can be set by mc.cores argument.

# code is only for demonstration
rl = run_all_consensus_partition_methods(mat, 
    top_value_method = c("sd", "MAD", ...),
    partition_method = c("hclust", "kmeans", ...),
    mc.cores = ...)
cola_report(rl, output_dir = ...)

There are many functions in cola package that can be applied to rl to visualize and compare the results for all combinations of methods simultaneously.

cola_report() function can be applied on rl to generate a HTML report for the complete analysis, with all the plots and tables generated.

Top-value methods

Top-value methods are used to assign scores to matrix rows, later the scores are ordered and only the top \(n\) rows with the highest scores are used for consensus partitioning. The default top-value methods provided in the package are:

all_top_value_methods()
## [1] "sd"  "cv"  "MAD" "ATC"

These top methods are:

These methods can be used in consensus partitioning by providing the name to the top_value_method argument in run_all_consensus_partition_methods() or consensus_partition().

You can register a new top-value method by register_top_value_methods(). The value should be functions. For each function, it should only have one argument which is the matrix for analysis and it must return a vector with scores for rows. In following example, the “max” method uses the row maximum as the row score and we also add the “QCD” (quartile coefficient of dispersion) method as a second method here.

register_top_value_methods(
    max = function(mat) matrixStats::rowMaxs(mat),
    QCD = function(mat) {
        qa = matrixStats::rowQuantiles(mat, probs = c(0.25, 0.75))
        (qa[, 2] - qa[, 1])/(qa[, 2] + qa[, 1])
    })
all_top_value_methods()
## [1] "sd"  "cv"  "MAD" "ATC" "max" "QCD"

By default, the consensus partition functions run_all_consensus_partition_methods() uses all registered top-value methods, but still you can explicitly specify a subset of top-value methods. To remove registered top-value methods, simply use remove_top_value_methods() by providing a vector of names.

remove_top_value_methods(c("max", "QCD"))
all_top_value_methods()
## [1] "sd"  "cv"  "MAD" "ATC"

The ATC method

Choosing the top rows in the matrix is important for the subgroup classification. In most cases, we extract the most variable rows which is defined by row variance. However, sometimes it won't give you meaningful rows which are efficient for subgroup classification. When random noise in the data increases, e.g. for single cell RNASeq data, the most variable genes are too weak to detect any stable subgroups.

If we think reversely, assuming there exist stable subgroups in the data, there must be groups of rows showing similar pattern to support the subgrouping, in other words, rows in the same groups should have high correlations to each other. Thus, if we can get rows that more correlated to others, they are more strong to form a stable subgroup for the samples. According to this thought, we designed the ATC method.

For row \(i\) in a matrix, \(X\) is a vector of the absolute correlation to all other rows, the ATC (ability to correlate to others) for row \(i\) is defined as:

\[ATC_i = 1 - \int_0^1F(x)\]

where \(F(x)\) is the empirical CDF (cumulative distribution function) of \(X\).

In following plot, the line is the CDF curve. ATC is the area above the CDF curve. It can be imagined that when row \(i\) correlates with more other rows, the CDF curve shifts more to the right, thus with higher ATC scores.

plot of chunk unnamed-chunk-9

There can be scenarios when large number of rows are correlated to each other only with very small correlation values. They will gain high ATC value due to the large number of rows (which corresponds to the left part of the read area in above plot that close to \(x = 0\)). To decrease such effect, the ATC definition can be slightly modified to:

\[ATC_i = (1 - \alpha) - \int_{\alpha}^1F(x^{\beta})\]

where now \(ATC_i\) is the red area only on the right of \(x = \alpha\). The coefficient \(\beta\) is the power added to the absolute correlations that it decreases more for the smaller correlations. By Default \(\alpha\) is set to 0 and \(\beta\) is set to 1.

Next we perform a simulation test to show the attributes of ATC method. A matrix with 160 rows, 100 columns with random values are generated as follows:

  1. 100 rows with mean of 0. The covariance matrix is set to 0 and 1 on the diagnal;
  2. 10 rows with mean of 0. The covariance matrix is set to 0.8 and 1 on the diagnal. This simulates high correlations but in a small group of rows;
  3. 50 rows wiht mean of 0. The covariance matrix is set to 0.5 and 1 on the diagnal. This simulates intemediate correlations but in a large group of rows.

The top left figure in following is the heatmap for the random matrix, split by the three groups of rows. In the top right figure, they are ECDF curves of the correlation when calculating ATC scores. The bottom left figure is the ATC scores for all 160 rows and the bottom right figure is the standard deviation for the 160 rows.

plot of chunk unnamed-chunk-10

All the 160 rows have similar variance of 1 that they cannot be distinguished very well by using variance (bottom right figure). As a contrast, the rows with non-zero covariance have higher ATC values (the red and green points), even higher when the number of correlated rows increases (the green points, although the correlation value itself is intermediate, bottom left figure). This shows ATC method can assign higher values for rows which correlate to more other rows.

ATC scores are calculated by ATC() function. By default it uses Pearson correlation. Users can register ATC method with other correlation methods by, e.g.:

# code is only for demonstration
register_top_value_methods(
    ATC_spearman = function(m) ATC(m, method = "spearman"),
    ATC_bicor = function(m) ATC(m, cor_fun = WGCNA::bicor)
)

Partition methods

Partition methods are used to separate samples into \(k\) subgroups where \(k\) is a known parameter for the partition. The default partition methods are:

all_partition_methods()
## [1] "hclust"  "kmeans"  "skmeans" "pam"     "mclust" 
## attr(,"scale_method")
## [1] "z-score" "z-score" "z-score" "z-score" "z-score"

These partition methods are:

Similarly, you can register a new partition method by register_partition_methods(). The value should be functions with two arguments which are the input matrix and number of partitions. There can be a third argument for the function which is ... used for passing more arguments from consensus_partition(). The function should only return a vector of subgroup/class labels or an object that can be imported by clue::cl_membership(). Please note the partition is applied on columns of the matrix and the number of unique levels of subgroup levels which are predicted by the partition method should not exceed \(k\).

Following example registers a partition method which randomly assign subgroup labels to samples:

register_partition_methods(
    random = function(mat, k) {
        sample(letters[1:k], ncol(mat), replace = TRUE)
    }
)

Here the subgroup labels can be in any types (numbers, characters). They only need to be different for different classes. These labels will be re-coded with numeric indices internally (i.e. 1, 2, 3, …).

To remove a partition method, use remove_partition_methods():

remove_partition_methods("random")

The built-in hclust method only uses Euclidean distance with “complete” clustering method, It is easy to define another hclust method:

# code is only for demonstration
register_partition_methods(
    hclust_cor = function(mat, k) cutree(hclust(as.dist(1 - cor(mat))), k)
)

Following code registers SOM and NMF partition methods:

# code is only for demonstration
library(kohonen)
register_partition_methods(
    SOM = function(mat, k, ...) {
        kr = floor(sqrt(ncol(mat)))
        somfit = som(t(mat), grid = somgrid(kr, kr, "hexagonal"), ...)
        m = somfit$codes[[1]]
        m = m[seq_len(nrow(m)) %in% somfit$unit.classif, ]
        cl = cutree(hclust(dist(m)), k)
        group = numeric(ncol(mat))
        for(cl_unique in unique(cl)) {
            ind = as.numeric(gsub("V", "", names(cl)[which(cl == cl_unique)]))
            l = somfit$unit.classif %in% ind
            group[l] = cl_unique
        }
        group
    }
)
library(NMF)
register_partition_methods(
    NMF = function(mat, k, ...) {
        fit = nmf(mat, rank = k, ...)
        apply(fit@fit@H, 2, which.max)
    }, scale_method = "max-min"
)

For these two methods, users can simply use register_SOM() and register_NMF() functions in cola package.

# code is only for demonstration
register_SOM()
register_NMF()

In the code above, there is an additional argument scale_method for register_partition_method(). scale_method controls how to scale the matrix rows before partitioning if scaling is turned on consensus_partition() or run_all_consensus_partition_methods(). There are three possible values:

The skmeans method

The skmeans method (the spherical k-means clustering) is powerful to detect subgroups where samples in a same subgroup show strong correlations. skmeans clustering uses cosine similarity and projects data points onto a unit hyper-sphere. As we have tested for many datasets, skmeans is very efficient to detect stable subgroups.

Consensus partitioning

For a given number of top rows \(n_i\), the corresponding matrix with top rows denoted as \(M_i\), a subset of rows with probability of \(p\) are randomly sampled from \(M_i\) and a certain partition method is applied on it, generating a partition \(P_a\). In most of cases, we have no prior knowledge of which \(n_i\) gives better results, thus, cola allows to try multiple \(n_i\) and pool partitions from all \(n_i\) together to find a consensus subgrouping, which also lets rows more on the top of the ranked list give higher weight for determining the final subgroups. Let's assume top rows are tried for \(n_1\), \(n_2\), …, \(n_m\) and the randomly sampling is performed for \(N_s\) times, then, for a given number of subgroups \(k\) for trying, the total number of partitions is \(N_P = m*N_s\).

Consensus matrix

The consensus matrix measures how consistently two samples are in a same subgroup and it can be used to visualize or analysis the stability of the subgrouping. The value \(c_{ij}\) in the consensus matrix is the probability of sample \(i\) and sample \(j\) in a same subgroup in all \(N_P\) partitions. It is calculated as:

\[c_{ij} = \sum_a^{N_p}I(s_{ia}, s_{ja})/N_P\]

where \(s_{ia}\) is the subgroup label for sample \(i\) in partition \(a\) and \(I()\) is the indicator function there \(I(x = y) = 1\) and \(I(x \neq y) = 0\).

Assuming there are stable subgroups in a dataset, which means, for any pair of samples \(i\) and \(j\) both in a same subgroup, \(c_{ij}\) is close to 1, and for any pairs in different subgroups, \(c_{ij}\) is close to 0, if the consensus matrix is visualized as a heatmap, samples in the same subgroup will be represented as a block in the diagonal of the heatmap.

Following two heatmaps visualize two consensus matrices. The left one shows less stability of subgrouping than the right one, while for the right one, there are three very obvious blocks in the diagnal that in each block, the corresponding samples are very likely to be in a same subgroup.

In following two heatmaps, the right one corresponds to a more stable subgrouping.

plot of chunk unnamed-chunk-18

Consensus subgroup labels

As long as we have a list of \(N_P\) partitions for a given subgroup number \(k\), we need to find a consensus partition based on all \(N_P\) partitions.Internally, cola package uses the clue package to construct the “partition ensemble” and predict the consensus subgroups. The “SE” method from clue::cl_consensus() function (please check the on-line documentation of this function) are used to calculate the consensus subgroup labels. The consensus subgroups are labels by integers (i.e. 1, 2, 3, …).

Adjust subgroup labels

The subgroup labels are assigned with numeric indices, however, in each partition, the assignment of the labels can be random, e.g. one same subgroup can be assigned with 1 in one partition, while in the other partition, it can be 2, but they are all identical for the sense of subgrouping. E.g. following partitions are all identical:

1 1 1 1 1 1 1 2 2 2 2 2 2
2 2 2 2 2 2 2 1 1 1 1 1 1
a a a a a a a b b b b b b

The subgroups are identical if switching the subgroup labels. This subgroup label adjustment is called the linear sum assignment problem, which is solved by the solve_LSAP() function in clue package. The aim is to find a mapping \(m()\) between two sets of labels to maximize \(\sum_i I(s_{1i}, m(s_{2i}))\) where \(s_1\) is the first label set and \(s_2\) is the second label set.

In following example, if the mapping is 1 -> 2, 2 -> 1, the second partition in following

1 1 1 1 1 1 1 2 2 2 2 2 # in partition 1
2 2 2 2 2 1 1 1 1 1 1 1 # in partition 2

is adjusted to

1 1 1 1 1 1 1 2 2 2 2 2
1 1 1 1 1 2 2 2 2 2 2 2 # switch 1 <-> 2

For the subgroups predicted by clue::cl_consensus(), the labels are additionally adjusted by the mean distance in each subgroup (calculated from the scaled data matrix), which means, the subgroup with label 1 always has the smallest mean intra-group distance.

This subgroup label adjustment is frequently used in cola to help the visualization as well as downstream analysis. E.g. for a specific combination of top-value method and partition method, the subgroup labels for different k are adjusted, and for the subgroups from different top-value methods and partition methods, the subgroup labels are also adjusted to make the label difference from different methods minimal.

Membership matrix

The \(N_P\) partitions are stored as a membership matrix where rows are partitions (grouped by top_n) and the subgroup labels in each partition are adjusted according to the consensus partition. Following heatmap is a visualization of all partitions and correspondence to the consensus partition for \(k = 2\). The p* annotation on top of the heatmap is the probability of being in subgroup \(i\) across all partitions.

In following plot, actually we can see the samples in the middle tend to belong to the green subgroup (group with label 1) for small top_n (e.g. top_n = 1000), while, when top_n increases, they go to the red subgroup (group label 2). Since the final subgroups are summarized from all top_n, the probabilities of the middle samples to be in the either subgroup are close, which can also be observed from the probability annotation (p* annotations). This might indicate these are the subset of samples which are in the intermediate state between group 1 and group 2.

Determine the best number of subgroups

Consensus partitioning is applied with a specific number of subgroups (we term it as \(k\)). Normally, a list of \(k\) are tried to find the best \(k\). cola provides metrics to help to determine the best number of subgroups.

The get_stats() function returns statistics for all metrics mentioned below and select_partition_number() plots the statistics with the number of subgroups (an example is here).

Silhouette score

The silhouette scores measures how close one sample is in its own subgroup compared to the closest neighbouring subgroup. For sample \(i\), the mean distance to every subgroups are calculated, denoted as \(d_1\), \(d_2\), …, \(d_k\) (\(d_k\) is the mean Euclidean distance between sample \(i\) and every sample in subgroup \(s\)). The distance to the subgroup where sample \(i\) stays is denoted as \(d_a\) and the silhouette score is defined as:

\[silhouette_i = 1 - d_a/d_b\]

where \(d_b\) is the minimal distance excluding \(d_a\):

\[d_b = min_{j \neq a}^k d_j\]

Following plot illustrates how silhouette score is calculated for sample x_i.

plot of chunk unnamed-chunk-19

The mean silhouette score from all samples is used to choose the best \(k\) where higher the mean silhouette score, better the \(k\).

PAC score

The PAC score measures the proportion of the ambiguous subgrouping. If the subgrouping is stable, in \(N_P\) partitions, sample \(i\) and sample \(j\), in most of the cases, are either always in a same subgroup, or always in different subgroups, which results in that, in the consensus matrix, the values are, in most cases, close to 1 or 0. Then in the CDF of the consensus matrix, the curve will be very flattened between \(x_1\) and \(x_2\) where \(x_1\) is very close to 0 and \(x_2\) is very close to 1 because there are very few values between \(x_1\) and \(x_2\). Thus, the proportion of sample pairs with consensus values in \((x_1, x_2)\) is called the proportion of the ambiguous clustering, which can be calculated by \(F(x_2) - F(x_1)\).

In following plots, the red line in the left plot corresponds to the first consensus heatmap and the blue line corresponds to the second consensus heatmap that is more stable than the first one.It is quite obvious to see the second consensus heatmap has far less PAC value than the first one.

plot of chunk unnamed-chunk-20

In some cases, \(F(x_1)\) or \(F(x_2)\) changes a lot when \(x_1\) has slight change around 0.1, or \(x_2\) has slight change around 0.9. Thus, to make PAC not so sensitive to the selection of \(x_1\) and \(x_2\), PAC value is calculated by removing 5% samples with lowest silhouette scores.

Smaller the PAC score, better the \(k\).

Concordance

The concordance of partitions to the consensus partition is calculated as, for each partition \(P_a\), the probability that it fits the consensus partition:

\[c_{a} = \frac{1}{N_s}\sum_i^{N_s}I(s_{ia} = sc_i)\]

where \(N_s\) is the number of samples, \(s_{ia}\) is the subgroup label of sample \(i\) in partition \(a\) and \(sc_i\) is the consensus subgroup label for sample \(i\). Note class labels in single partitions have already been adjusted to the consensus partition labels.

The final concordance score is the mean value of \(c_a\). Higher the concordance score, better the \(k\).

Area increased

It is the increased area under CDF of the consensus matrix compared to the previous \(k\).

\[A_k = \int F_k(x) - \int F_{k-1}(x)\]

and when \(k = 2\) or for the minimal \(k\):

\[A_k = \int F_k(x)\]

In follow example, there are five consensus heatmaps corresponding to \(k =2,3,4,5,6\). Note the number of subgroups that can be inferred from the consensus heatmap is not necessary to be exactly the same as \(k\). It can be smaller than \(k\).

plot of chunk unnamed-chunk-21

The corresponding CDF curves and the area increased are:

plot of chunk unnamed-chunk-22

The \(k\) before the elbow is taken as the best \(k\) (in above example it is 3). Basically when \(k\) reaches a stable subgrouping, increasing \(k\) won't change the consensus matrix too much, which results in less change of the difference of area under the CDF curve.

Rand index

In some cases, when number of subgroups changes from \(k-1\) to \(k\), all the statistics imply \(k\) is a better choice than \(k-1\). However, when observing the consensus heatmap, basically it is because a very small set of samples are separated to form a new subgroup. In this case, it is better to still keep \(k-1\) subgroups. In other words, the subgrouping with \(k\) is similar as \(k-1\) and it is not worth to increase \(k\) from \(k-1\). In cola package, there are two metrics: Rand index and Jaccard index to measure the similarity of two partitions for \(k-1\) and \(k\). The two metrics are calculated by clue::cl_agreement(..., method = "Rand") and clue::cl_agreement(..., method = "Jaccard").

For all pairs of samples, denote following symbols (https://en.wikipedia.org/wiki/Rand_index#Definition):

the Rand index which is the percent of pairs of samples that are both in a same cluster or both are not in a same cluster in the partition of \(k\) and \(k-1\).

\[Rand = \frac{a+b}{a+b+c+d}\]

If Rand index is too high, it means the two subgroupings are very similar and it is not sufficient to increase from \(k-1\) to \(k\).

Jaccard index

The Jaccard index is the ratio of pairs of samples that are both in a same subgroup in the partition of \(k\) and \(k-1\) and the pairs of samples are both in a same subgroup in the partition of \(k\) or \(k-1\).

\[Jaccard = \frac{a}{a+c+d}\]

In following plots, when the number of subgroups increases from 3 to 4, there is only one single sample separated from other subgroups to form a new subgroup. The Rand index or the Jaccard index for \(k=4\) is close to 1, which means, the subgroups at \(k=4\) are highly similar as \(k=3\), thus we ignore \(k=4\) and take \(k=3\) as the better subgrouping.

Rule

cola provides a suggest_best_k() function which suggests the best \(k\). It is based on following rules:

suggest_best_k() only gives suggestion on selecting a reasonable \(k\). Users still need to look at the plots (e.g. by select_partition_number() or consensus_heatmap() functions), or even by checking whether the subgrouping gives a reasonable signatures by get_signatures(), to pick a reasonable \(k\) that best explains their study.

Find signatures

As long as there are stable subgroups, we can look for rows which show distinct difference in one subgroup compared to others. They can be called signature genes or signature CpG sites if the corresponding dataset is gene expression data or methylation data.

By default, samples with silhouette scores less than 0.5 are removed. cola provides following methods:

Users can also provide their own method by providing a function with the matrix and subgroup labels as inputs and a vector of FDR as output.

Gene ontology Enrichment

If rows in the matrix can be associated to genes, GO_enrichment() can be applied to perform GO enrichment analysis to the signature genes.

Compare multiple methods

consensus_partition() is the core function for consensus partitioning. But it can only perform analysis with a single top-value method and a single partition method. In most cases, we have no idea of which combination of top-value method and partition method gives better results. Here run_all_consensus_partition_methods() can perform analysis with multiple methods simultaneously:

# code is only for demonstration
rl = run_all_consensus_partition_methods(mat, 
    top_value_method = c("sd", "MAD", ...),
    partition_method = c("hclust", "kmeans", ...),
    mc.cores = ...)

By default it runs analysis for all combinations of top-value methods in all_top_value_methods() and partition methods in all_partition_methods().

cola package provides functions to collect plots from all combinations of methods for straightforward comparisons.

# code is only for demonstration
collect_plots(rl, fun = consensus_heatmap, k = ...)
collect_plots(rl, fun = membership_heatmap, k = ...)
collect_plots(rl, fun = get_signatures, k = ...)

And collect_classes() compares consensus partition from all methods:

# code is only for demonstration
collect_classes(rl, k = ...)

The plots from collect_plots() and collect_classes() can be found here.

Implementation of the package

cola is implemented in an object-oriented way. There are two main classes where ConsensusPartition class contains results for a single top-value method and a single partition method, while ConsensusPartitionList class contains results for multiple top-value methods and multiple partition methods.

In following example code, TCGA_GBM_subgroup.rds is generated by the code demonstrated here (Running cola analysis is always time-consuming, so we use the object that has already be generated.).

download.file("https://jokergoo.github.io/cola_examples/TCGA_GBM/TCGA_GBM_subgroup.rds", 
    destfile = "TCGA_GBM_subgroup.rds", quiet = TRUE)
rl = readRDS("TCGA_GBM_subgroup.rds")
file.remove("TCGA_GBM_subgroup.rds")

Simply entering the variable name gives you the summary of the analysis and a list of functions that can be applied to this object:

rl
## A 'ConsensusPartitionList' object with 24 methods.
##   On a matrix with 11268 rows and 173 columns.
##   Top rows are extracted by 'sd, cv, MAD, ATC' methods.
##   Subgroups are detected by 'hclust, kmeans, skmeans, pam, mclust, NMF' method.
##   Number of partitions are tried for k = 2, 3, 4, 5, 6.
##   Performed in total 30000 partitions by row resampling.
## 
## Following methods can be applied to this 'ConsensusPartitionList' object:
##  [1] "GO_enrichment"         "cola_report"           "collect_classes"       "collect_plots"        
##  [5] "collect_stats"         "colnames"              "get_anno"              "get_anno_col"         
##  [9] "get_classes"           "get_matrix"            "get_membership"        "get_stats"            
## [13] "ncol"                  "nrow"                  "rownames"              "show"                 
## [17] "suggest_best_k"        "test_to_known_factors" "top_rows_heatmap"      "top_rows_overlap"     
## 
## You can get result for a single method by, e.g. object["sd", "hclust"] or object["sd:hclust"]
## or a subset of methods by object[c("sd", "cv")], c("hclust", "kmeans")]

To get results for a single top-value method and partition method, you can subset rl by the name of the combination of the methods.

rl["sd:hclust"]
## A 'ConsensusPartition' object with k = 2, 3, 4, 5, 6.
##   On a matrix with 11268 rows and 173 columns.
##   Top rows (1000, 2000, 3000, 4000, 5000) are extracted by 'sd' method.
##   Subgroups are detected by 'hclust' method.
##   Performed in total 1250 partitions by row resampling.
##   Best k for subgroups seems to be 2.
## 
## Following methods can be applied to this 'ConsensusPartition' object:
##  [1] "GO_enrichment"           "cola_report"             "collect_classes"        
##  [4] "collect_plots"           "collect_stats"           "colnames"               
##  [7] "compare_signatures"      "consensus_heatmap"       "dimension_reduction"    
## [10] "get_anno"                "get_anno_col"            "get_classes"            
## [13] "get_consensus"           "get_matrix"              "get_membership"         
## [16] "get_param"               "get_signatures"          "get_stats"              
## [19] "membership_heatmap"      "ncol"                    "nrow"                   
## [22] "plot_ecdf"               "rownames"                "select_partition_number"
## [25] "show"                    "suggest_best_k"          "test_to_known_factors"

Functions on ConsensusPartitionList class that are important to use are:

Functions on ConsensusPartition class:

Visualizations

cola package provides rich visualizations for the results generated by a single method or multiple methods.

On the ConsensusPartition object

The object which is generated with a single top-value method and a single partition method belongs to the class ConsensusPartition. There are several visualization functions that can be applied to it. select_partition_number() makes several plots to show different statistics along with different \(k\), which helps to determine the “best k”.

res = rl["MAD:kmeans"] # the ConsensusPartition object
select_partition_number(res)

plot of chunk unnamed-chunk-29

The heatmap for the consensus matrix with a certain \(k\):

consensus_heatmap(res, k = 4)

plot of chunk unnamed-chunk-30

The heatmap for the membership matrix with a certain \(k\):

membership_heatmap(res, k =4)

plot of chunk unnamed-chunk-31

The PCA plot with a certain \(k\):

dimension_reduction(res, k = 4)
## use UMAP

plot of chunk unnamed-chunk-32

The heatmap for the signature rows with a certain \(k\). The heatmap is split into two parts by columns. The left heatmap where the barplots on top are in black contains samples with silhouette scores larger than 0.5 and the right heatmap where the barplot son top are in grey contains samples with silhouette scores less than 0.5. Rows are automatically split by k-means.

get_signatures(res, k = 4)

plot of chunk unnamed-chunk-33

collect_classes() which is applied on the ConsensusPartition object visualizes how subgroups are corresponded with increasing \(k\). Same row in all heatmaps corresponds to a same sample.

collect_classes(res)

plot of chunk unnamed-chunk-34

collect_plots() which is applied on the ConsensusPartition object puts all the plots from all \(k\) into one single page.

collect_plots(res)

plot of chunk unnamed-chunk-35

On the ConsensusPartitionList object

run_all_consensus_partition_methods() returns a ConsensusPartitionList object. There are two main functions which can visualize results from all combinations of methods and compare directly.

collect_plots(rl, fun = consensus_heatmap, k = 4)

plot of chunk unnamed-chunk-36

fun can also be membership_heatmap or get_signatures that membership heatmap and signature heatmap for each method will be plotted.

collect_classes() which is applied on the ConsensusPartitionList object plots the partition for each combination of methods and the lightness correspond to the silhouette scores for samples in each method. Rows are clustered by the dissimilarity measurement from clue::cl_dissimilarity(..., method = "comembership"). On top the consensus subgroup is inferred from all methods by taking the mean silhouette scores as weight.

collect_classes(rl, k = 4)

plot of chunk unnamed-chunk-37

collect_stats() helps to compare statistics from multiple methods and multiple metrics.

collect_stats(rl, k = 4)

plot of chunk unnamed-chunk-38

The HTML report

All the content introduced above is mainly for the deep understanding of the package. In real data analysis, users do not need to type that amount of code. cola_report() function wraps all the code and performs the complete analysis automatically. Normally, applying cola analysis, following three lines of code are enough for you.

# code is only for demonstration
mat = adjust_matrix(mat) # for some datasets, you don't need this line.
rl = run_all_consensus_partition_methods(mat, mc.cores = ...)
cola_report(rl, output_dir = ...) # Alles ist da!

A Real-world example

In this example, we use the TCGA GBM microarray data set where four subtypes is predicted. The two files (unifiedScaled.txt and TCGA_unified_CORE_ClaNC840.txt) for use is from here. The web page for the analysis is https://jokergoo.github.io/cola_examples/.

Following code is used to perform analysis of consensus partition.

# code is only for demonstration
mat = read.table("https://jokergoo.github.io/cola_examples/unifiedScaled.txt", 
    header = TRUE, row.names = 1, check.names = FALSE)
mat = as.matrix(mat)

subtype = read.table("https://jokergoo.github.io/cola_examples/TCGA_unified_CORE_ClaNC840.txt", 
    sep = "\t", header = TRUE, check.names = FALSE, stringsAsFactors = FALSE)
subtype = structure(unlist(subtype[1, -(1:2)]), names = colnames(subtype)[-(1:2)])

mat = mat[, names(subtype)]

mat = adjust_matrix(mat)
rl = run_all_consensus_partition_methods(mat, top_n = c(1000, 2000, 3000, 4000), 
    max_k = 6, mc.cores = 4, anno = data.frame(subtype = subtype), 
    anno_col = list(subtype = structure(seq_len(4), names = unique(subtype))))

Simply typing rl gives a summary of the analysis:

rl
## A 'ConsensusPartitionList' object with 24 methods.
##   On a matrix with 11268 rows and 173 columns.
##   Top rows are extracted by 'sd, cv, MAD, ATC' methods.
##   Subgroups are detected by 'hclust, kmeans, skmeans, pam, mclust, NMF' method.
##   Number of partitions are tried for k = 2, 3, 4, 5, 6.
##   Performed in total 30000 partitions by row resampling.
## 
## Following methods can be applied to this 'ConsensusPartitionList' object:
##  [1] "GO_enrichment"         "cola_report"           "collect_classes"       "collect_plots"        
##  [5] "collect_stats"         "colnames"              "get_anno"              "get_anno_col"         
##  [9] "get_classes"           "get_matrix"            "get_membership"        "get_stats"            
## [13] "ncol"                  "nrow"                  "rownames"              "show"                 
## [17] "suggest_best_k"        "test_to_known_factors" "top_rows_heatmap"      "top_rows_overlap"     
## 
## You can get result for a single method by, e.g. object["sd", "hclust"] or object["sd:hclust"]
## or a subset of methods by object[c("sd", "cv")], c("hclust", "kmeans")]

The HTML report for the analysis is available at https://jokergoo.github.io/cola_examples/TCGA_GBM/TCGA_GBM_subgroup_cola_report/cola_report.html.

Hierarchical partitioning

Normal consensus partition methods aim to find \(k\) subgroups at the same time. However, when 1. there are dominant subgroups, or 2. the number of potential subgroups are large, it is difficult to find secondary subgroups which show less difference with normal consensus partition process. To solve this problem, the subgroups can be found in a hierarchical way, where dominant subgroups are found first, and secondary subgroups are detected afterwards.

cola package implements hierarchical partition and the flowchart is as follows.

Generally, at each recursive step, the consensus partition is performed for a small set of \(k\) (because the final larger number of subgroups will be detected hierarchically) and a subset of samples, If the PAC score of the “best k” is less than 0.2, the samples are split into two subgroups where the first group contains samples with largest distance to all other subgroups and the second groups are all other samples. For each of the two groups, if the number of samples is less than 6, the hierarchical partition stops, or it repeatedly performs the hierarchical partition on the subset of samples in the corresponding groups.

Hierarchical partition is performed by hierarchical_partition() function. Note you can only use one single top-value method and a single partition methods here.

In following example, we still use the TCGA GBM microarray datasets. The consensus partition which is summarized from all methods (from run_all_consensus_partition_methods()) are added as an annotation to compare to the subgroups predicted by hierarchical partition.

# code is only for demonstration
library(RColorBrewer)
rh = hierarchical_partition(mat, top_n = c(1000, 2000, 3000, 4000),
    top_value_method = "MAD", partition_method = "kmeans",
    anno = data.frame(subtype = subtype,
        consensus = get_classes(rl, k = 4)$class), 
    anno_col = list(subtype = structure(seq_len(4), names = unique(subtype)),
        consensus = structure(brewer.pal(4, "Set1"), names = 1:4)))

rh is already generated and can be downloaded from here.

download.file("https://jokergoo.github.io/cola_examples/TCGA_GBM/TCGA_GBM_subgroup_hierarchical_partition.rds", 
  destfile = "TCGA_GBM_subgroup_hierarchical_partition.rds", quiet = TRUE)
rh = readRDS("TCGA_GBM_subgroup_hierarchical_partition.rds")
file.remove("TCGA_GBM_subgroup_hierarchical_partition.rds")

Simply typing rh gives the summary of the analysis.

rh
## A 'HierarchicalPartition' object with 'ATC:skmeans' method.
##   On a matrix with 11268 rows and 173 columns.
##   Performed in total 21750 partitions by row resampling.
##   There are 15 groups.
## 
## Hierarchy of the partition:
##   0, 173 cols
##   |-- 02, 53 cols
##   |   |-- 021, 16 cols
##   |   `-- 020, 37 cols
##   |       |-- 0201, 21 cols
##   |       |   |-- 02011, 9 cols
##   |       |   `-- 02010, 12 cols
##   |       `-- 0200, 16 cols
##   `-- 00, 120 cols
##       |-- 001, 43 cols
##       |   |-- 0012, 11 cols
##       |   `-- 0010, 32 cols
##       |       |-- 00103, 9 cols
##       |       `-- 00100, 23 cols
##       |           |-- 001001, 8 cols
##       |           `-- 001000, 15 cols
##       `-- 000, 77 cols
##           |-- 0001, 37 cols
##           |   |-- 00011, 21 cols
##           |   |   |-- 000111, 10 cols
##           |   |   `-- 000110, 11 cols
##           |   `-- 00010, 16 cols
##           |       |-- 000101, 7 cols
##           |       `-- 000100, 9 cols
##           `-- 0000, 40 cols
##               |-- 00001, 22 cols
##               `-- 00000, 18 cols
##                   |-- 000001, 7 cols
##                   `-- 000000, 11 cols
## 
## Following methods can be applied to this 'HierarchicalPartition' object:
##  [1] "GO_enrichment"         "all_leaves"            "all_nodes"             "cola_report"          
##  [5] "collect_classes"       "collect_plots"         "colnames"              "dimension_reduction"  
##  [9] "get_anno"              "get_anno_col"          "get_classes"           "get_matrix"           
## [13] "get_signatures"        "max_depth"             "ncol"                  "nrow"                 
## [17] "rownames"              "show"                  "suggest_best_k"        "test_to_known_factors"
## 
## You can get result for a single node by e.g. object["01"]

cola uses a special way to encode the node in the hierarchy. The length of the node name is the depth of the node in the hierarchy and the substring excluding the last digit is the node name of the parent node. E.g. for the node 0011, the depth is 4 and the parent node is 001. For each node, the last digit represents the index of the subgroup detected in parent node. and 0 represents all other subgroups. E.g. 02 means this node contains samples which belong to subgroup 2 when doing consensus partition at node 0 (subgroup 2 is selected as a node because the intra-group distance is smallest compared to other subgroups), and 00 contains samples which do not belong to subgroup 2 at node 0.

On each node of the partition hierarchy, it is a ConsensusPartition object (note it is not a sub-hierarchy) and it can be retrieved by using node name as index.

rh["00"]
## A 'ConsensusPartition' object with k = 2, 3, 4.
##   On a matrix with 11268 rows and 120 columns.
##   Top rows (1000, 2000, 3000, 4000, 5000) are extracted by 'ATC' method.
##   Subgroups are detected by 'skmeans' method.
##   Performed in total 750 partitions by row resampling.
##   Best k for subgroups seems to be 3.
## 
## Following methods can be applied to this 'ConsensusPartition' object:
##  [1] "GO_enrichment"           "cola_report"             "collect_classes"        
##  [4] "collect_plots"           "collect_stats"           "colnames"               
##  [7] "compare_signatures"      "consensus_heatmap"       "dimension_reduction"    
## [10] "get_anno"                "get_anno_col"            "get_classes"            
## [13] "get_consensus"           "get_matrix"              "get_membership"         
## [16] "get_param"               "get_signatures"          "get_stats"              
## [19] "membership_heatmap"      "ncol"                    "nrow"                   
## [22] "plot_ecdf"               "rownames"                "select_partition_number"
## [25] "show"                    "suggest_best_k"          "test_to_known_factors"

collect_classes() plots the hierarchy of subgroups as well as the annotations that are set before.

collect_classes(rh)

plot of chunk unnamed-chunk-46

In above plot, generally hierarchical partition found similar subgroups as consensus partition, but hierarchical partition additionally found two subgroups for the Mesenchymal subtype samples, which makes totally 5 subgroups. However, if we directly check the 5 subgroups in rl, actually the subgrouping is not stable (see following consensus heatmap). This means the difference between the two subgroups in Mesenchymal is so small that it cannot be distinguished if we take all samples in the analysis.

consensus_heatmap(rl["MAD:kmeans"], k = 5)

plot of chunk unnamed-chunk-47

In the hierarchical partition, the two subgroups of Mesenchymal subtype are under the node 001. If only applying consensus partition on node 001, actually there are two obvious subgroups and there are quite a lot signature genes which are significantly different between the two subgroups. This proves it is very meanningful that there are two secondary subgroups in Mesenchymal subtype.

get_signatures(rh["001"], k = 2)

plot of chunk unnamed-chunk-48

get_sigatures() can also be applied to the rh object to visualize the signature rows found in every node in the partition hierarchy.

get_signatures(rh, verbose = FALSE)

plot of chunk unnamed-chunk-49

Similar as the ConsensusPartitionList object, cola_report() function can also be applied to this HierarchicalPartition object. The full report for rh can be found at https://jokergoo.github.io/cola_examples/TCGA_GBM/TCGA_GBM_subgroup_hierarchical_partition_cola_report/cola_hc.html.

# code is only for demonstration
cola_report(rh, output_dir = "~/your-output-folder/tcga_cola_rh_report")

Following functions can be applied to the HierarchicalPartition class:

Session info

sessionInfo()
## R version 3.6.1 (2019-07-05)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.3 LTS
## 
## Matrix products: default
## BLAS:   /home/biocbuild/bbs-3.9-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.9-bioc/R/lib/libRlapack.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8       
##  [4] LC_COLLATE=C               LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C              
## [10] LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] grid      stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] GetoptLong_0.1.7     mvtnorm_1.0-11       matrixStats_0.55.0   circlize_0.4.8      
## [5] ComplexHeatmap_2.0.0 cola_1.0.1           knitr_1.24           markdown_1.1        
## 
## loaded via a namespace (and not attached):
##  [1] bitops_1.0-6         bit64_0.9-7          RColorBrewer_1.1-2   httr_1.4.1          
##  [5] data.tree_0.7.8      tools_3.6.1          backports_1.1.4      R6_2.4.0            
##  [9] lazyeval_0.2.2       DBI_1.0.0            BiocGenerics_0.30.0  colorspace_1.4-1    
## [13] gridExtra_2.3        tidyselect_0.2.5     bit_1.1-14           compiler_3.6.1      
## [17] Biobase_2.44.0       xml2_1.2.2           influenceR_0.1.0     microbenchmark_1.4-6
## [21] slam_0.1-45          scales_1.0.0         readr_1.3.1          genefilter_1.66.0   
## [25] askpass_1.1          stringr_1.4.0        digest_0.6.20        pkgconfig_2.0.2     
## [29] htmltools_0.3.6      umap_0.2.3.1         highr_0.8            htmlwidgets_1.3     
## [33] rlang_0.4.0          GlobalOptions_0.1.0  rstudioapi_0.10      RSQLite_2.1.2       
## [37] impute_1.58.0        visNetwork_2.0.8     shape_1.4.4          jsonlite_1.6        
## [41] mclust_5.4.5         dendextend_1.12.0    dplyr_0.8.3          rgexf_0.15.3        
## [45] RCurl_1.95-4.12      magrittr_1.5         Matrix_1.2-17        Rcpp_1.0.2          
## [49] munsell_0.5.0        S4Vectors_0.22.1     reticulate_1.13      viridis_0.5.1       
## [53] lifecycle_0.1.0      stringi_1.4.3        blob_1.2.0           parallel_3.6.1      
## [57] crayon_1.3.4         lattice_0.20-38      splines_3.6.1        annotate_1.62.0     
## [61] hms_0.5.1            zeallot_0.1.0        pillar_1.4.2         igraph_1.2.4.1      
## [65] rjson_0.2.20         stats4_3.6.1         XML_3.98-1.20        glue_1.3.1          
## [69] evaluate_0.14        downloader_0.4       png_0.1-7            vctrs_0.2.0         
## [73] gtable_0.3.0         openssl_1.4.1        purrr_0.3.2          tidyr_1.0.0         
## [77] clue_0.3-57          assertthat_0.2.1     ggplot2_3.2.1        xfun_0.9            
## [81] eulerr_5.1.0         xtable_1.8-4         skmeans_0.2-11       RSpectra_0.15-0     
## [85] viridisLite_0.3.0    survival_2.44-1.1    tibble_2.1.3         AnnotationDbi_1.46.1
## [89] memoise_1.1.0        IRanges_2.18.2       cluster_2.1.0        Rook_1.1-1          
## [93] DiagrammeR_1.0.1     brew_1.0-6