1.1.1 The RSEC Routine

The RSEC algorithm (Resampling-based Sequential Ensemble Clustering) was developed specifically for finding robust clusters in single cell sequencing data. It follows our above suggested workflow, but makes specific choices along the way that we find to be advantageous in clustering large mRNA expression datasets, like those of single cell sequencing. It is implemented in the function RSEC. This vignette serves to show how the pieces of the workflow operate, and use choices that are shown are not those recommended by RSEC. This is both to show the flexibility of the functions and because RSEC is computationally more time-intensive, since it is based on both resampling, and sequential discovery of clusters.

2.3.1 The output

The output of clusterMany is a ClusterExperiment object. This is a class built for this package and explained in the section on ClusterExperiment Objects. In the object ce the clusters are stored, names and colors for each cluster within a clustering are assigned, and other information about the clusterings is recorded. Furthermore, all of the information in the original SummarizedExperiment is retained.

There are many accessor functions that help you get at the information in a ClusterExperiment object and some of the most relevant are described in the section on ClusterExperiment Objects. (ClusterExperiment objects are S4 objects, and are not lists).

For right now we will only mention the most basic such function that retrieves the actual cluster assignments. This is the clusterMatrix function, that returns a matrix where the columns are the different clusterings and the rows are samples. Within a clustering, the clusters are encoded by consecutive integers.

head(clusterMatrix(ce)[,1:3])

##      nVAR=100,k=5 nVAR=500,k=5 nVAR=1000,k=5
## [1,]            1            1             1
## [2,]            2            1             1
## [3,]            3            2             2
## [4,]            1            1             1
## [5,]            4            3             2
## [6,]            1            4             2

Because we have changed clusterLabels above, these new cluster labels are shown here. Notice that some of the samples are assigned the value of -1. -1 has the significance of encoding samples that are not assigned to any cluster. Why certain samples are not clustered depends on the underlying choices of the clustering routine. In this case, the default in clusterMany set the minimum size of a cluster to be 5, which resulted in -1 assignments.

Another special value is -2 discussed in the section on ClusterExperiment objects

4.1.1 Saving the alignment of plotClusters

plotClusters invisibly returns a ClusterExperiment object. In our earlier calls to plotCluster, this would be the same as the input object and so there is no reason to save it. However, the alignment and color assignments created by plotClusters can be requested to be saved via the resetNames, resetColors and resetOrderSamples arguments. If any of these are set to TRUE, then the object returned will be different than those of the input. Specifically, if resetColors=TRUE the colorLegend of the returned object will be changed so that the colors assigned to each cluster will be as were shown in the plot. Similarly, if resetNames=TRUE the names of the clusters will be changed to be integer values, but now the integers will be aligned to try to be the same across clusters (and therefore not consecutive integers, which is why these are saved as names for the clusters and not clusterIds). If resetOrderSamples=TRUE, then the order of the samples shown in the plot will be similarly saved in the slot orderSamples.

Here we make a call to plotClusters, but now ask to reset everything to match this ordering. We’ll save this as a different object so that this is not a permanent change for the rest of the vignette.

par(mar=plotCMar)
ce_temp<-plotClusters(ce,whichClusters="workflow", sampleData=c("Biological_Condition","Cluster2"), 
               main="Clusters from clusterMany, different order",axisLine=-1,resetNames=TRUE,resetColors=TRUE,resetOrderSamples=TRUE)

clusterLegend(ce_temp)[c("mergeClusters","combineMany,final")]

## $mergeClusters
##    clusterIds color     name
## -1 "-1"       "white"   "-1"
## 1  "1"        "#1F78B4" "1" 
## 2  "2"        "#33A02C" "2" 
## 3  "3"        "#E31A1C" "3" 
## 4  "4"        "#FF7F00" "4" 
## 5  "5"        "#6A3D9A" "5" 
## 6  "6"        "#B15928" "6" 
## 
## $`combineMany,final`
##    clusterIds color     name
## -1 "-1"       "white"   "-1"
## 1  "1"        "#1F78B4" "1" 
## 2  "2"        "#33A02C" "2" 
## 3  "3"        "#E31A1C" "3" 
## 4  "4"        "#FF7F00" "4" 
## 5  "5"        "#6A3D9A" "5" 
## 6  "6"        "#B15928" "6"

Now, the clusterLegend slot of the object no longer has the default color/name assignments, but it has names and colors that match across the clusters. Notice, that this means the prefix “m” or “c” that was previously given to distinguish the combineMany result from the mergeClusters result is now gone (the user could manually add them by changing the clusterLegend values).

We can also force plotClusters to use the existing color definitions, rather than create its own. This makes particular sense if you want to have continuity between plots – i.e. be sure that a particular cluster always has a certain color – but would like to do different variations of plotClusters to get a sense of how similar the clusters are.

For example, we set the colors above based on the cluster order from plotClusters where the clusters were ordered according to the workflow. But now we want to plot only the clusters from clusterMany, yet keep the same colors as before so we can compare them. We do this by setting the argument existingColors="all", meaning use all of the existing colors (currently this is the only option available for how to use the existing colors).

par(mar=plotCMar)
par(mfrow=c(1,2))
plotClusters(ce_temp, sampleData=c("Biological_Condition","Cluster2"),
             whichClusters="workflow", existingColors="all",
             main="Workflow Clusters, fix the color of clusters",axisLine=-1)

plotClusters(ce_temp, sampleData=c("Biological_Condition","Cluster2"), 
             existingColors="all", whichClusters="clusterMany",
             main="clusterMany Clusters, fix the color of clusters",
             axisLine=-1)

4.2.1 Displaying clustering or sample information

Like plotClusters, plotHeatmap has a whichClusters option that behaves similarly to that of plotClusters. In addition to the options “all” and “workflow” that we saw with plotClusters, plotHeatmap also takes the option “none”" (no clusters shown) and “primary” (only the primaryCluster). The user can also request a subset of the clusters by giving specific indices to whichClusters like in plotClusters.

Here we create a heatmap that shows the clusters from the workflow. Notice that we choose only the last 2 – from combineMany and mergeClusters. If we chose all “workflow” clusters it would be too many.

whClusterPlot<-1:2
plotHeatmap(ce,whichClusters=whClusterPlot, annLegend=FALSE)

Notice we also passed the option ‘annLegend=FALSE’ to the underlying aheatmap command (with many clusterings shown, it is often not useful to have a legend for all the clusters because the legend doesn’t fit on the page!). The many detailed commands of aheatmap that are not set internally by plotHeatmap can be passed along as well.

Like plotClusters, plotHeatmap takes an argument sampleData, which refers to columns of the colData of that object and can be included.

4.2.2 Additional options for clustering/ordering samples

We can choose to not cluster the samples, but order the samples by cluster. This time we’ll just show the primary cluster (the mergeCluster result) by setting whichClusters="primaryCluster":

plotHeatmap(ce,clusterSamplesData="primaryCluster",
            whichClusters="primaryCluster",
            main="Heatmap with clusterMany",annLegend=FALSE)

As an improvement upon this, we can cluster the clusters into a dendrogram so that the most similar clusters will be near each other. We already did this before with our call to makeDendrogram. We haven’t done anything to change that, so the dendrogram from that call is still stored in the object. We can check this in the information shown in our object:

show(ce)

## class: ClusterExperiment 
## dim: 7069 65 
## Primary cluster type: mergeClusters 
## Primary cluster label: mergeClusters 
## Table of clusters (of primary clustering):
## -1 m1 m2 m3 m4 m5 m6 
## 12  9  7 15 14  4  4 
## Total number of clusterings: 40 
## Dendrogram run on 'combineMany,final' (cluster index: 2)
## -----------
## Workflow progress:
## clusterMany run? Yes 
## combineMany run? Yes 
## makeDendrogram run? Yes 
## mergeClusters run? Yes

We can see, as we expected, that the dendrogram we made from “combineMany,final” is still the active dendrogram saved in the ce object. Now we will call plotHeatmap choosing to display the “mergeClusters” and “combineMany,final” clustering, and ordering the samples by the dendrogram in ce. We will also display some data about the samples from the colData slot using the sampleData argument. Notice that instead of giving the index, we can also give the clusterLabels of the clusters we want to show.

plotHeatmap(ce,clusterSamplesData="dendrogramValue",
            whichClusters=c("mergeClusters","combineMany"),
            main="Heatmap with clusterMany",
            sampleData=c("Biological_Condition","Cluster2"),annLegend=FALSE)

If there is not a dendrogram stored, plotHeatmap will call makeDendrogram based on the primary cluster (with the default settings of makeDendrogram); calling makeDendrogram on ce is preferred so that the user can control the choices in how it is done (which we will discuss below). For visualization purposes, the dendrogram for the combineMany cluster is preferred to that of the mergeCluster cluster, since “combineMany,final” is just a finer partition of the “mergeClusters” clustering.

4.2.3 Using separate input data for clustering and for visualization

While count data is a common type of data, it is also common that the input data in the SummarizedExperiment object might be normalized data from a normalization package such as RUVSeq. In this case, the clustering and all numerical calculations should be done on the normalized data (which may or may not need a log transform). However, these normalized data might not be on a logical count scale (for example, in RUVSeq, the normalize data are residuals after subtracting out gene-specific batch effects).

In this case, it can be convenient to have the visualization of the data (i.e. the color scale), be based on a count scale that is interpretable, even while the clustering is done based on the normalized data. This is possible by giving a new matrix of values to the argument visualizeData. In this case, the color scale (and clustering of the features) is based on the input visualizeData matrix, but all clustering of the samples is done on the internal data in the ClusterExperiment object.

5.1.1 Overview of the implemented clustering procedures

In the quick start section we picked some simple and familiar clustering options that would run quickly and needed little explanation. However, our workflow generally assumes more complex options and more parameter variations are tried. Before getting into the specific options of clusterMany, let us first describe some of these more complicated setups, since many of the arguments of clusterMany depend on understanding them.

Clustering Algorithms (clusterD): Clustering algorithms generally start with a particular predefined distance or dissimilarity between the samples, usually encoded in a \(n x n\) matrix, \(D\). In our package, we consider only such clustering algorithms. The input could also be similarities, though we will continue to call such a matrix \(D\).

The simplest scenario is a simple calculation of \(D\) and then a clustering of \(D\), usually dependent on particular parameters. In this package we try to group together algorithms that cluster \(D\) based on common parameters and operations that can be done. Currently there are two “types” of algorithms we consider, which we call type “K” and “01”. The help page of clusterD documents these choices more fully, but we give an overview here. Many of the parameters that are allowed to vary in clusterMany refer to parameters for the clustering of the \(D\) matrix, either for “K” or “01” type algorithms.

The “K” algorithms are so called because their main parameter requirement is that the user specifies the number of clusters (\(K\)) to be created. They assume the input \(D\) are dissimilarities, and depending on the algorithm, may have additional expectations. “pam” and “kmeans” are examples of such types of algorithms.

The “01” algorithms are so named because the algorithm assumes that the input \(D\) consists of similarities between samples and that the similarities encoded in \(D\) are on a scale of 0-1. They use this fact to make the primary user-specified parameter be not the number of final clusters, but a measure \(\alpha\) of how dissimilar samples in the same cluster can be (on a scale of 0-1). Given \(\alpha\), the algorithm then implements a method to then determine the clusters (so \(\alpha\) implicitly determines \(K\)). These methods rely on the assumption that because the 0-1 scale has special significance, the user will be able to make an determination more easily as to the level of dissimilarity allowed in a true cluster, rather than predetermine the number of clusters \(K\). The current 01 methods are “tight”, “hierarchical01” and “pam”.

Subsampling In addition to the basic clustering algorithms on a matrix \(D\), we also implement many other common cluster processing steps that are relevant to the result of the clustering. We have already seen such an example with dimensionality reduction, where the input \(D\) is determined based on different input data. A more significant processing that we perform is calculation of a dissimilarity matrix \(D\) not by a distance on the vector of data, but by subsampling and clustering of the subsampled data. The resulting \(D\) matrix is a matrix of co-clustering percentages. Each entry is a percentage of subsamples where the two samples shared a clustering over the many subsamplings of the data (there are slight variations as how this can be calculated, see help pages of subsampleClustering ). Note that this implies there actually two different clustering algorithms (and sets of corresponding parameters) – one for the clustering on the subsampled data, and one for the clustering of the resulting \(D\) of the percentage of coClustering of samples. clusterMany focuses on varying parameters related to the clustering of \(D\), and generally assumes that the underlying clustering on the subsampled data is simple (e.g. “pam”). (The underlying clustering machinery in our package is performed by a function clusterSingle which is not discussed in this tutorial, but can be called explicitly for fine-grained control over all of the features ). The subsampling option is computationally expensive, and when coupled with comparing many parameters, does result in a lengthy evaluation of clusterMany. However, we recommend it as one of the most useful methods for getting stable clustering results.

Sequential Detection of Clusters Another complicated addition to the clustering that requires additional explanation is the implementation of sequential clustering. This refers to clustering of the data, then removing the “best” cluster, and then re-clustering the remaining samples, and then continuing this iteration until all samples are clustered (or the algorithm in some other way calls a stop). Such sequential clustering can often be convenient when there is very dominant cluster, for example, that is far away from the other mass of data. Removing samples in these clusters and resampling can sometimes be more productive and result in results more robust to the choice of samples. A particular implementation of such a sequential method, based upon (Tseng and Wong 2005), is implemented in the clusterExperiment package. Because of the iterative nature of this process, there are many possible parameters (see help(seqCluster)). clusterMany does not allow variation of very many of these parameters, but instead just has the choice of running the sequential clustering or not. Sequential clustering can also be quite computationally expensive, particularly when paired with subsampling to determine \(D\) at each step of the iteration.

5.1.2 Arguments of clusterMany

Now that we’ve explained the underlying architecture of the clustering provided in the package, we discuss the parameters that can be varied in clusterMany. There are additional arguments available for clusterMany but right now we focus on only the ones that can be given multiple options. Recall that arguments in clusterMany that take on multiple values mean that the combinations of all the multiple valued arguments will be given as input for a clustering routine.

• sequential This parameter consists of logical values, TRUE and/or FALSE, indicating whether the sequential strategy should be implemented or not.
• subsample This parameter consists of logical values, TRUE and/or FALSE, indicating whether the subsampling strategy for determining \(D\) should be implemented or not.
• clusterFunction The clustering functions to be tried. If subsample=TRUE is part of the combination, then clusterFunction the method that will be used on the co-clustering matrix \(D\) created from subsampling the data (where pam clustering is used on the subsampled data). Otherwise, clusterFunction is the clustering method that will be used on dist(t(x))
• ks The argument ‘ks’ is interpreted differently for different choices of the other parameters. If sequential=TRUE is part of the combination, ks defines the argument k0 of sequential clustering (see help(seqCluster)), which is approximately like the initial starting point for the number of clusters in the sequential process. Otherwise, ks is passed to set \(K\) of both the clustering of subsampled data and the actual clustering of the data (if a \(K\) needs to be set, i.e. a type “K” algorithm). When/if findBestK=TRUE is part of the combination, ks also defines the range of values to search for the best k (see the details in help(clusterMany) for more).
• dimReduce These are character strings indicating what choices of dimensionality reduction should be tried. The choices are “PCA”, indicating clustering on the top principal components, “var”, indicating clustering on the top most variable features, and “none”, indicating the whole data set should be used. If either “PCA” or “var” are chosen, the following parameters indicate the number of such features to be used (and can be a vector of values as we have seen):
  • nVarDims
  • nPCADims
• distFunction These are character values giving functions that provide a distance matrix between the samples, when applied to the data. These functions should be accessible in the global environment (clusterMany applies get to the global environment to access these functions). To make them compatible with the dist function, these functions should assume the samples are in the rows, i.e. they should work when applied to t(assay(ce)). We give an example in the next subsection below.
• minSizes these are integer values determining the minimum size required for a cluster (passed to the clusterD part of clustering).
• alphas These are the \(\alpha\) parameters for “01” clustering techniques; these values are only relevant if one of the clusterFunction values is a “01” clustering algorithm. The values given to alphas should be between 0 and 1, with smaller values indicating greater similarity required between the clusters.
• betas These are the \(\beta\) parameters for sequential clustering; these values are only relevant if sequential=TRUE and determine the level of stability required between changes in the parameters to determine that a cluster is stable.
• findBestK This option is for “K” clustering techniques, and indicates that \(K\) should be chosen automatically as the \(K\) that gives the largest silhouette distance between clusters.
• removeSil A logical value as to whether samples with small silhouette distance to their assigned cluster are “removed”, in the sense that they are not given their original cluster assignment but instead assigned -1. This option is for “K” clustering techniques as a method of removing poorly clustered samples (the “01” techniques used by clusterMany generally do this intrinsically as part of the algorithm).
• silCutoff If removeSil is TRUE, then silCutoff determines the cutoff on silhouette distance for unassigning the sample.

clusterMany tries to have generally simple interface, and for this reason makes choices about what is meant by certain combinations. For example, in combinations where findBestK=TRUE, ks=2:10 is taken to mean that the clustering should find the best \(k\) out of the range of 2-10. However, in other combinations ks might indicate the specific number of clusters, \(k\), that should be found. For parameter combinations that are not what is desired, the user should consider making direct calls to clusterSingle where all of these options (and many more) can be explicitly called.

Other parameters for the clustering are kept fixed. As described above, there are many more possible parameters in play than are considered in clusterMany. These parameters can be set via the arguments clusterDArgs, subsampleArgs and seqArgs. These arguments correspond to the different processes described above (the clustering of \(D\), the creation of \(D\) via subsampling, and the sequential clustering process, respectively). These arguments take a list of arguments that are sent directly to clusterSingle. However, these arguments may be overridden by the interpretation of clusterMany of how different combinations interact; again for complete control direct calls to clusterSingle are necessary.

5.1.3 Example changing the distance function and clustering algorithm

Providing different distance functions is slightly more involved than the other parameters, so we give an example here.

First we define distances that we would like to compare. We are going to define two distances that take values between 0-1 based on different choices of correlation.

corDist<-function(x){(1-cor(t(x),method="pearson"))/2}
spearDist<-function(x){(1-cor(t(x),method="spearman"))/2}

These distances are defined so as to give distance of 0 between samples with correlation 1, and distance of 1 for correlation -1.

We will also compare using different algorithms for clustering. Since we chose distances between 0-1, we can use any algorithm. Currently, clusterMany requires that the distances work with all of the clusterFunction choices given. Since some of the clusterFunction algorithms require a distance matrix between 0-1, this means we can only compare all of the algorithms when the distance is a 0-1 distance. (Future versions will probably try to create a work around so that the algorithm just skips algorithms that don’t match the distance).

Note on 0-1 clustering when subsample=FALSE We would note that the default values \(\alpha\) for the 0-1 clustering were set with the distance \(D\) the result of subsampling or other concensus summary in mind. In generally, subsampling creates a \(D\) matrix with high similarity for many samples who share a cluster (the proportion of times samples are seen together for well clustered samples can easily be in the .8-.95 range, or even exactly 1). For this reason the default \(\alpha\) is 0.1 which requires distances between samples in the 0.1 range or less (i.e. a similarity in the range of 0.9 or more). We show an example of the \(D\) matrix from subsampling; we make use of the clusterSingle which is the workhorse mentioned above that runs a single clustering command directly, which gives the output \(D\) from the sampling in the “coClustering” slot of ce. Note that the result is \(1-p_{ij}\) where \(p_{ij}\) is the proportion of times sample \(i\) and \(j\) clustered together.

ceSub<-clusterSingle(ce,dimReduce="mad",ndims=1000,subsample=TRUE,clusterFunction="hierarchical01",subsampleArgs=list(k=8),clusterLabel="subsamplingCluster",clusterDArgs=list(minSize=5))
plotCoClustering(ceSub,colorScale=rev(seqPal5))

We see even here, the default of \(\alpha=0.1\) was perhaps too conservative since only two clusters came out (with size greater than 5).

The distances based on correlation calculated directly on the data, such as we created above, are often used for clustering expression data. But they are unlikely to have distances as low as seen in subsampling, even for well clustered samples. Here’s a visualization of the correlation distance matrix we defined above (using Spearman’s correlation) on the top 1000 most variable features:

dSp<-spearDist(t(transform(ce,dimReduce="mad",nVarDims=1000)))
plotHeatmap(dSp,isSymmetric=TRUE,colorScale=rev(seqPal5))

We can see that the choice of \(\alpha\) must be much higher (and we are likely to be more sensitive to it).

Notice to calculate the distance in the above plot, we made use of the transform function applied to our ce object to get the results of dimensionality reduction. The transform function gave us a data matrix back that has been transformed, and also reduced in dimensions, like would be done in our clustering routines. transform has similar parameters as seen in clusterMany,makeDendrogram or clusterSingle and is useful when you want to manually apply something to transformed and/or dimensionality reduced data; and you can be sure you are getting the same matrix of data back that the clustering algorithms are using.

Comparing distance functions with clusterMany Now that we have defined the distances we want to compare in our global environment, we can give these to the argument “distFunction” in clusterMany. They should be given as character strings giving the names of the functions. For computational ease for this vignette, we will just choose the dimensionality reduction to be the top 1000 features based on MAD and set K=8 or \(\alpha=0.45\). We will save these results as a separate object so as to not disrupt the earlier workflow.

ceDist<-clusterMany(ce, k=7:9, alpha=c(0.35,0.4,0.45), clusterFunction=c("tight","hierarchical01","pam","hierarchicalK"),findBestK=FALSE,
               removeSil=c(FALSE),dist=c("corDist","spearDist"),
               isCount=TRUE,dimReduce=c("mad"),nPCADims=1000,run=TRUE)
clusterLabels(ceDist)<-gsub("clusterFunction","alg",clusterLabels(ceDist))
clusterLabels(ceDist)<-gsub("Dist","",clusterLabels(ceDist))
clusterLabels(ceDist)<-gsub("distFunction","dist",clusterLabels(ceDist))
clusterLabels(ceDist)<-gsub("hierarchical","hier",clusterLabels(ceDist))
par(mar=c(1.1,15.1,1.1,1.1))
plotClusters(ceDist,axisLine=-2,sampleData=c("Biological_Condition"))

Notice that using the 01 methods did not give relevant results

5.1.4 Dealing with large numbers of clusterings

A good first check before running clusterMany is to determine how many clusterings you are asking for. clusterMany has some limited internal checks to not do unnecessary duplicates (e.g. removeSil only works with some clusterFunctions so clusterMany would detect that), but generally takes all combinations. This can take a while for more complicated clustering techniques, so it is a good idea to check what you are getting into. You can do this by running clusterMany with run=FALSE.

In the following we consider expanding our original clustering choices to consider individual choices of \(K\) (rather than just findBestK=TRUE).

checkParam<-clusterMany(se, clusterFunction="pam", ks=2:10,
                        removeSil=c(TRUE,FALSE), isCount=TRUE,
                        dimReduce=c("PCA","var"),
                        nVarDims=c(100,500,1000),nPCADims=c(5,15,50),run=FALSE)
dim(checkParam$paramMatrix) #number of rows is the number of clusterings

## [1] 108  12

Each row of the matrix checkParam$paramMatrix is a requested clustering (the columns indicate the value of a possible parameter). Our selections indicate 108 different clusterings (!).

We can set ncores argument to have these clusterings done in parallel. If ncores>1, the parallelization is done via mclapply and should not be done in the Rgui interface (see help pages for mclapply).

5.3.1 makeDendrogram

makeDendrogram creates a hierarchical clustering of the clusters as determined by the primaryCluster of the ClusterExperiment object. In addition to being used for merging clusters, the dendrograms created by makeDendrogram are also useful for ordering the clusters in plotHeatmap as has been shown above.

makeDendrogam performs hierarchical clustering of the cluster medoids (after transformation of the data) and provides a dendrogram that will order the samples according to this clustering of the clusters. The hierarchical ordering of the dendrograms is saved internally in the ClusterExperiment object.

Like clustering, the dendrogram can depend on what features are included from the data. The same options for clustering are available for the hierarchical clustering of the clusters, namely choices of dimensionality reduction via dimReduce and the number of dimensions via ndims.

ce<-makeDendrogram(ce,dimReduce="var",ndims=500)
plotDendrogram(ce)

Notice that the plot of the dendrogram shows the hierarchy of the clusters (and color codes them according to the colors stored in colorLegend slot).

Recall that the most recent clustering made is from our call to combineMany, where we experimented with using on some of the clusterings from clusterMany, so that is our current primaryCluster:

show(ce)

## class: ClusterExperiment 
## dim: 7069 65 
## Primary cluster type: combineMany 
## Primary cluster label: combineMany,nVAR 
## Table of clusters (of primary clustering):
## -1 c1 c2 c3 c4 c5 c6 c7 
## 10  5  5  9 15 12  4  5 
## Total number of clusterings: 41 
## Dendrogram run on 'combineMany,nVAR' (cluster index: 1)
## -----------
## Workflow progress:
## clusterMany run? Yes 
## combineMany run? Yes 
## makeDendrogram run? Yes 
## mergeClusters run? No

This is the clustering from combining only the clusterings from clusterMany that use the top most variable genes. Because it is the primaryCluster, it was the clustering that was used by default to make the dendrogram.

We might prefer to get back to the dendrogram based on our combineMany in quick start (the “combineMany, final” clustering). We’ve lost that dendrogram when we called makeDendrogram again. However, we can rerun makeDendrogram and choose a different clustering from which to make the dendrogram.

ce<-makeDendrogram(ce,dimReduce="var",ndims=500,
                   whichCluster="combineMany,final")
plotDendrogram(ce)

Note that the clusterType of this clustering is not “combineMany”, but “combineMany.x”, where “x” indicates what iteration it was:

clusterTypes(ce)[which(clusterLabels(ce)=="combineMany,final")]

## [1] "combineMany.3"

We can choose reset this past call to combineMany to the current ‘combineMany’ output (which will also set this clustering to be the primaryCluster). We don’t need to recall makeDendrogram, since we already used this clustering to make the dendrogram by our explicit call of the argument whichCluster.

ce<-setToCurrent(ce,whichCluster="combineMany,final")
show(ce)

## class: ClusterExperiment 
## dim: 7069 65 
## Primary cluster type: combineMany 
## Primary cluster label: combineMany,final 
## Table of clusters (of primary clustering):
## -1 c1 c2 c3 c4 c5 c6 
## 12  9  7 15 14  4  4 
## Total number of clusterings: 41 
## Dendrogram run on 'combineMany,final' (cluster index: 3)
## -----------
## Workflow progress:
## clusterMany run? Yes 
## combineMany run? Yes 
## makeDendrogram run? Yes 
## mergeClusters run? No

5.3.2 Merging clusters with little differential expression

We then can use this hierarchy of clusters to merge clusters that show little difference in expression. We do this by testing, for each node of the dendrogram, for which features is the mean of the set of clusters to the right split of the node is equal to the mean on the left split. This is done via the getBestFeatures (see section on getBestFeatures), where the type argument is set to “Dendro”.

Starting at the bottom of the tree, those clusters that have the percentage of features with differential expression below a certain value (determined by the argument cutoff) are merged into a larger cluster. This testing of differences and merging continues until the estimated percentage of non-null DE features is above cutoff. This means lower values of cutoff result in less merging of clusters. There are multiple methods of estimation of the percentage of non-null features implemented. The option mergeMethod="adjP" which we showed earlier is the simplest: the proportion found significant by calculating the proportion of DE genes at a given False Discovery Rate threshold (using the Benjamini-Hochberg procedure). However, other methods are also implemented (see the help of mergeClusters).

Notice that mergeClusters will always run based on the clustering that made the currently existing dendrogram. So it is always good to check that it is what we expect.

ce

## class: ClusterExperiment 
## dim: 7069 65 
## Primary cluster type: combineMany 
## Primary cluster label: combineMany,final 
## Table of clusters (of primary clustering):
## -1 c1 c2 c3 c4 c5 c6 
## 12  9  7 15 14  4  4 
## Total number of clusterings: 41 
## Dendrogram run on 'combineMany,final' (cluster index: 3)
## -----------
## Workflow progress:
## clusterMany run? Yes 
## combineMany run? Yes 
## makeDendrogram run? Yes 
## mergeClusters run? No

mergeClusters can also be run without merging the cluster, and simply drawing a plot showing the dendrogram along with the estimates of the percentage of non-null features to aid in deciding a cutoff and method. By setting plotType="all", all of the estimates of the different methods are displayed simultaneously, while before in the QuickStart, we only showed the default values.

mergeClusters(ce,mergeMethod="none",plotType="all")

## Note: Merging will be done on ' combineMany,final ', with clustering index 3

## Warning in locfdr::locfdr(tstats, plot = 0): f(z) misfit = 2.9. Rerun with
## increased df

## Warning in locfdr::locfdr(tstats, plot = 0): f(z) misfit = 1.6. Rerun with
## increased df

Now we can pick a cutoff. We’ll give it a label to keep it separate from the previous run we had made.

ce<-mergeClusters(ce,cutoff=0.05,mergeMethod="adjP",clusterLabel="mergeClusters,v2")

## Note: Merging will be done on ' combineMany,final ', with clustering index 3

## Warning in locfdr::locfdr(tstats, plot = 0): f(z) misfit = 2.9. Rerun with
## increased df

## Warning in locfdr::locfdr(tstats, plot = 0): f(z) misfit = 1.6. Rerun with
## increased df

ce

## class: ClusterExperiment 
## dim: 7069 65 
## Primary cluster type: mergeClusters 
## Primary cluster label: mergeClusters,v2 
## Table of clusters (of primary clustering):
## -1 m1 m2 m3 
## 12 31  7 15 
## Total number of clusterings: 42 
## Dendrogram run on 'combineMany,final' (cluster index: 4)
## -----------
## Workflow progress:
## clusterMany run? Yes 
## combineMany run? Yes 
## makeDendrogram run? Yes 
## mergeClusters run? Yes

If we want to rerun mergeClusters with a different cutoff, we can go ahead and do that (notice in the summary above, that the dendrogram is still there, and associated with the same cluster, so we can rerun mergeClusters and know that it will still pull from the same input clustering from combineMany as we were before)

ce<-mergeClusters(ce,cutoff=0.15,mergeMethod="MB",
                  clusterLabel="mergeClusters,v3")

## Note: Merging will be done on ' combineMany,final ', with clustering index 4

## Warning in locfdr::locfdr(tstats, plot = 0): f(z) misfit = 2.9. Rerun with
## increased df

## Warning in locfdr::locfdr(tstats, plot = 0): f(z) misfit = 1.6. Rerun with
## increased df

ce

## class: ClusterExperiment 
## dim: 7069 65 
## Primary cluster type: mergeClusters 
## Primary cluster label: mergeClusters,v3 
## Table of clusters (of primary clustering):
## -1 m1 m2 m3 m4 m5 m6 
## 12  9  7 15 14  4  4 
## Total number of clusterings: 43 
## Dendrogram run on 'combineMany,final' (cluster index: 5)
## -----------
## Workflow progress:
## clusterMany run? Yes 
## combineMany run? Yes 
## makeDendrogram run? Yes 
## mergeClusters run? Yes

5.4.1 Designate a Final Clustering

A final protected clusterTypes is “final”. This is not created by any method, but can be set to be the clusterType of a clustering by the user (via the clusterTypes command). Any clustering marked final will be considered one of the workflow values for commands like whichClusters="workflow". However, they will NOT be renamed with “.x” or removed if eraseOld=TRUE. This is a way for a user to ‘save’ a clustering as important/final so it will not be changed internally by any method, yet still have it show up with the “workflow” clustering results. There is no limit to the number of such clusters that are so marked, but the utility of doing so will drop if too many such clusters are chosen.

For best functionality, particularly if a user has determined a single final clustering after completing clustering, a user will probably want to set the primaryClusterIndex to be that of the final cluster and rerun makeDendrogram. This will help in plotting and visualizing. The setToFinal command does this.

Here we will demonstrate marking a cluster as final. We go back to our previous mergeClusters that we found with cutoff=0.05 and mark it as our final clustering. First we need to find which cluster it is. We see from our above call to the workflow functions above, that it is clusterType equal to “mergeClusters.4” and label equal to “mergeClusters,v2”. In our call to setToFinal we will decide to change it’s label as well.

ce<-setToFinal(ce,whichCluster="mergeClusters,v2",
               clusterLabel="Final Clustering")
par(mar=plotCMar)
plotClusters(ce,whichClusters="workflow")

Note that because it is labeled as “final” it shows up automatically in “workflow” clusters in our plotClusters plot. It has also been set as our primaryCluster and has the new clusterLabel we gave it in the call to setToFinal.

This didn’t get rid of our undesired mergeClusters result that is most recent. It still shows up as “the” mergeClusters result. This might be undesired. We could remove that “mergeClusters” result with removeClusters. Alternatively, we could manually change the clusterTypes to mergeClusters.x so that it doesn’t show up as current. A cleaner way to do this would have been to first set the desired cluster (“mergeClusters.4”) to the most current iteration with setToCurrent, which would have bumped up the existing mergeClusters result to be no longer current.

6.1.1 All Pairwise

The option type="Pairs", which we saw earlier, performs all pair-wise tests between the clusters for each feature, testing for each pair of clusters whether the mean of the feature is different between the two clusters. Here is the example from above using all pairwise comparisons:

pairsAllTop<-getBestFeatures(ce,contrastType="Pairs",p.value=0.05)
dim(pairsAllTop)

## [1] 150   9

head(pairsAllTop)

##   IndexInOriginal Contrast Feature     logFC   AveExpr         t
## 1            6565    X1-X2     VIM -8.608470  7.366498 -7.374382
## 2            3309    X1-X2    MIAT  7.875230  7.987637  6.479908
## 3             552    X1-X2  BCL11B  9.186117  4.120686  6.468115
## 4            4501    X1-X2 PRTFDC1 -7.635334  1.913974 -6.235150
## 5            2210    X1-X2    GLI3 -7.460028  5.121621 -6.189493
## 6            1284    X1-X2   CXADR  6.796402 10.516764  5.723578
##        P.Value    adj.P.Val         B
## 1 5.322547e-10 3.865591e-07 10.834597
## 2 1.815383e-08 7.103104e-06  8.032926
## 3 1.901249e-08 7.357625e-06  7.995912
## 4 4.724549e-08 1.522698e-05  7.265297
## 5 5.643513e-08 1.765221e-05  7.122301
## 6 3.408406e-07 7.608638e-05  5.669800

Notice that compared to the quick start guide, We didn’t set the parameter number which is passed to topTable, so we can get out at most 10 significant features for each contrast/comparison (because the default value of number in topTable is 10). Similarly, if we didn’t set a value for p.value, topTable would return the top number genes per contrast, regardless of whether they were all significant or not. These are the defaults of topTable, which we purposefully do not modify, but we urge the user to read the documentation of topTable carefully to understand what is being asked for. In the QuickStart, we set number=NROW(ce) to make sure we got all significant genes.

In addition to the columns provided by topTable, the column “Contrast” tells us what pairwise contrast the result is from. “X1-X2” means a comparison of cluster 1 and cluster 2. The column “IndexInOriginal” gives the index of the gene to the original input data matrix, namely assay(ce). The other columns are given by topTable (with the column “Feature” renamed – it is usually “ProbeID” in limma).

6.1.2 One Against All

The choice type="OneAgainsAll" performs a comparison of a cluster against the mean of all of the other clusters.

best1vsAll<-getBestFeatures(ce,contrastType="OneAgainstAll",p.value=0.05)
head(best1vsAll)

##   IndexInOriginal ContrastName              Contrast Feature     logFC
## 1            6565            1 X1-(X2+X3+X4+X5+X6)/5     VIM -5.569268
## 2            5753            1 X1-(X2+X3+X4+X5+X6)/5   STRAP -5.808820
## 3            2960            1 X1-(X2+X3+X4+X5+X6)/5    LDHA -6.447274
## 4            6090            1 X1-(X2+X3+X4+X5+X6)/5 TMEM258 -4.881551
## 5             552            1 X1-(X2+X3+X4+X5+X6)/5  BCL11B  5.678299
## 6            2350            1 X1-(X2+X3+X4+X5+X6)/5   GSTP1 -5.005114
##    AveExpr         t      P.Value    adj.P.Val        B
## 1 7.366498 -6.372501 2.764171e-08 6.776853e-06 8.270257
## 2 7.975761 -6.358567 2.918907e-08 6.877918e-06 8.223424
## 3 7.200328 -5.436302 1.015365e-06 1.469819e-04 5.159886
## 4 6.563927 -5.423463 1.065741e-06 1.516858e-04 5.117965
## 5 4.120686  5.340430 1.456449e-06 1.936990e-04 4.847522
## 6 8.755384 -5.125553 3.245503e-06 3.781724e-04 4.153399

Notice that now there is both a “Contrast” and a “ContrastName” column. Like before, “Contrast” gives an explicit definition of what is the comparisons, in the form of “(X2+X3+X4+X5+X6)/5-X1”, meaning the mean of the means of clusters 2-6 is compared to the mean of cluster1. “ContrastName” interprets this into a more usable name, namely that this contrast can be easily identified as a test of “cluster 1”.

6.1.3 Dendrogram

The option type="Dendro" is more complex; it assumes that there is a hierarchy of the clusters (created by makeDendrogram and stored in the ClusterExperiment object). Then for each node of the dendrogram, getBestFeatures defines a contrast or comparison of the mean expression between the daughter nodes.

Notice the problem if we try to run the command, however.

bestDendroTest1<-getBestFeatures(ce,contrastType="Dendro",p.value=0.05)

## Error in .local(x, ...): Primary cluster does not match the cluster on which the dendrogram was made. Either replace existing dendrogram with on using the primary cluster (via 'makeDendrogram'), or reset primaryCluster with 'primaryClusterIndex' to be equal to index of 'dendo_index' slot

This is because the dendrogram we have does not correspond to the primaryCluster.

show(ce)

## class: ClusterExperiment 
## dim: 7069 65 
## Primary cluster type: mergeClusters 
## Primary cluster label: mergeClusters,v3 
## Table of clusters (of primary clustering):
## -1 m1 m2 m3 m4 m5 m6 
## 12  9  7 15 14  4  4 
## Total number of clusterings: 43 
## Dendrogram run on 'combineMany,final' (cluster index: 5)
## -----------
## Workflow progress:
## clusterMany run? Yes 
## combineMany run? Yes 
## makeDendrogram run? Yes 
## mergeClusters run? Yes

getBestFeatures only works on the primaryCluster, and has noticed that the dendrogram in our ce object was not made from the primaryCluster. So we need to rerun makeDendrogram with our primaryCluster. It’s regardless generally a good idea to make sure the dendrogram you have is the one you want.

ce<-makeDendrogram(ce,dimReduce="var",ndims=500)
bestDendro<-getBestFeatures(ce,contrastType="Dendro",p.value=0.05)
head(bestDendro)

##   IndexInOriginal ContrastName                  Contrast Feature     logFC
## 1            2210        Node1 (X4+X5+X1+X6)/4-(X2+X3)/2    GLI3 -8.212906
## 2            5745        Node1 (X4+X5+X1+X6)/4-(X2+X3)/2   STMN2  8.524708
## 3            4285        Node1 (X4+X5+X1+X6)/4-(X2+X3)/2  PLXNA2  8.197781
## 4            1465        Node1 (X4+X5+X1+X6)/4-(X2+X3)/2    DLK1 -5.338176
## 5            3886        Node1 (X4+X5+X1+X6)/4-(X2+X3)/2  NUDCD2 -7.065937
## 6            5276        Node1 (X4+X5+X1+X6)/4-(X2+X3)/2   SFRP1 -7.849545
##    AveExpr          t      P.Value    adj.P.Val        B
## 1 5.121621 -11.139064 2.528864e-16 2.234567e-12 25.24467
## 2 8.742343  11.057577 3.418369e-16 2.416445e-12 24.98301
## 3 5.572182  10.071251 1.393355e-14 6.156015e-11 21.73829
## 4 3.227247  -9.162249 4.633852e-13 1.191666e-09 18.63158
## 5 4.109444  -8.995111 8.893410e-13 2.087345e-09 18.04976
## 6 5.533309  -8.979601 9.449010e-13 2.087345e-09 17.99562

Again, there is both a “ContrastName” and “Contrast” column. The “Contrast” column identifies which clusters where on each side of the node (and hence commpared) and “ContrastName” is a name for the node.

levels((bestDendro)$Contrast)

## [1] "(X4+X5+X1+X6)/4-(X2+X3)/2" "X4-(X5+X1+X6)/3"          
## [3] "X2-X3"                     "X5-(X1+X6)/2"             
## [5] "X1-X6"

We can plot the dendrogram showing the node names to help make sense of which contrasts go with which nodes (plotDendrogram calls plot.phylo from the ape package and can take those arguments).

plotDendrogram(ce,show.node.label=TRUE)