library(clusterExperiment)
system.time(rsecFluidigm<-RSEC(se, isCount = TRUE, 
    reduceMethod="PCA", nReducedDims=10,combineMinSize=3,
    ncores=1, random.seed=176201))

#don't call this routine if you ran the above code.
#it will overwrite the rsecFluidigm you made
library(clusterExperiment)
data("rsecFluidigm")

2.3.1 The output

We can look at the object that was created.

rsecFluidigm

## class: ClusterExperiment 
## dim: 7069 65 
## reducedDimNames: no reduced dims stored 
## filterStats: no valid filtering stats stored 
## -----------
## Primary cluster type: mergeClusters 
## Primary cluster label: mergeClusters 
## Table of clusters (of primary clustering):
## -1 m1 m2 m3 m4 m5 m6 
## 13 24 15  4  3  3  3 
## Total number of clusterings: 38 
## Dendrogram run on 'combineMany' (cluster index: 2)
## -----------
## Workflow progress:
## clusterMany run? Yes 
## combineMany run? Yes 
## makeDendrogram run? Yes 
## mergeClusters run? Yes

The print out tells us about the clustering(s) that were created (namely 38 clusterings) and which steps of the workflow have been called (all of them have because we used the wrapper RSEC that does the whole thing). Recall from our brief description above that RSEC clusters the data many times using different parameters before finding an consensus clustering. All of these intermediate clusterings are saved. Each of these intermediate clusterings used subsampling of the data and sequential clustering of the data to find the clustering, while the different clusterings represent the different parameters that were adjusted.

We can see that rsecFluidigm has a built in (S4) class called a ClusterExperiment object. This is a class built for this package and explained in the section on ClusterExperiment Objects. In the object rsecFluidigm the clusterings are stored along with corresponding information for each clustering. Furthermore, all of the information in the original SummarizedExperiment is retained. The print out also tells us information about the “primaryCluster” of rsecFluidigm. Each ClusterExperiment object has a “primaryCluster”, which is the default cluster that the many functions will use unless specified by the user. We are told that the “primaryCluster” for rsecFluidigm is has the label “mergeClusters” – which is the defaul label given to the last cluster of the RSEC function because the last call of the RSEC function is to mergeClusters.

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. The final clustering created by RSEC is saved as the primary clustering of the object.

head(primaryCluster(rsecFluidigm),20)

##  [1] -1 -1  1  1 -1 -1 -1  2  2  2  1  1 -1  1  3  4 -1  2  3 -1

The clusters are encoded by consecutive integers. Notice that some of the samples are assigned the value of -1. -1 is the value assigned in this package for samples that are not assigned to any cluster. Why certain samples are not clustered depends on the underlying choices of the clustering routine and we won’t get into here until we learn a bit more about RSEC. Another special value is -2 discussed in the section on ClusterExperiment objects

This final result of RSEC is the result of running many clusterings and finding the ensembl consensus between them. All of the these intermediate clusterings are saved in rsecFluidigm object. They can be accessed by the clusterMatrix function, that returns a matrix where the columns are the different clusterings and the rows are samples. We show a subset of this matrix here:

head(clusterMatrix(rsecFluidigm)[,1:4])

##            mergeClusters combineMany k0=4,alpha=0.1 k0=5,alpha=0.1
## SRR1275356            -1          -1             -1             -1
## SRR1275251            -1          -1             -1             -1
## SRR1275287             1           1              1              1
## SRR1275364             1           2             -1             -1
## SRR1275269            -1          -1             -1             -1
## SRR1275263            -1          -1             -1             -1

The “mergeClusters” clustering is the final clustering from RSEC and matches the primary clustering that we saw above. The “combineMany” clustering is the result of the initial ensembl concensus among all of the many clusterings that were run, while “mergeClusters” is the result of merging smaller clusters together that did not show enough signs of differences between clusters. The remaining clusters are the result of changing the parameters, and a couple of such clusterings a shown in the above printout of the cluster matrix.

The column names are the clusterLabels for each clustering and can be accessed (and assigned new values!) via the clusterLabels function.

head(clusterLabels(rsecFluidigm))

## [1] "mergeClusters"  "combineMany"    "k0=4,alpha=0.1" "k0=5,alpha=0.1"
## [5] "k0=6,alpha=0.1" "k0=7,alpha=0.1"

We can see the names of more clusterings, and see that the different parameter values tried in each clustering are saved in the names of the clustering. We can also see the different parameter combinations that went into the consensus clustering by using getClusterManyParams (here only 2 different parameters).

head(getClusterManyParams(rsecFluidigm))

##                clusteringIndex k alpha
## k0=4,alpha=0.1               3 4   0.1
## k0=5,alpha=0.1               4 5   0.1
## k0=6,alpha=0.1               5 6   0.1
## k0=7,alpha=0.1               6 7   0.1
## k0=8,alpha=0.1               7 8   0.1
## k0=9,alpha=0.1               8 9   0.1

2.3.2 Visualizing the output

clusterExperiment also provides many ways to visualize the output of RSEC (or any set of clusterings run in clusterExperiment, as we’ll show below).

Visualizing many clusterings The first such useful visualization is a plot of all of the clusterings together using the plotClusters command. For this visualization, it is useful to change the amount of space on the left of the plot to allow for the labels of the clusterings, so we will reset the mar option in par. We also decrease the axisLine argument that decides the amount of space between the axis and the labels to give more space to the labels (axisLine is passed internally to the line option in axis).

defaultMar<-par("mar")
plotCMar<-c(1.1,8.1,4.1,1.1)
par(mar=plotCMar)
plotClusters(rsecFluidigm,main="Clusters from RSEC", whichClusters="workflow", sampleData=c("Biological_Condition","Published2"), axisLine=-1)

This plot shows the samples in the columns, and different clusterings on the rows. Each sample is color coded based on its clustering for that row, where the colors have been chosen to try to match up clusters across different clusterings that show large overlap. Moreover, the samples have been ordered so that each subsequent clustering (starting at the top and going down) will try to order the samples to keep the clusters together, without rearranging the clustering blocks of the previous clustering/row.

We also added a sampleData argument in our call, indicating that we also want to visualize some information about the samples saved in the colData slot (inherited from our original fluidigm object). We chose the columns “Biological_Condition” and “Published2” from colData, which correspond to the original biological condition of the experiment, and the clusters reported in the original paper, respectively. The data from sampleData (when requested) are always shown at the bottom of the plot.

Notice that some samples are white. This indicates that they have the value -1, meaning they were not clustered. In fact, for many clusterings, there is a large amount of white here. This is likely do to the fact that there are only 65 cells here, and the default parameters of RSEC are better suited for a large number of cells, such as seen in more modern single-cell sequencing experiments. The sequential clustering may be problematic for small numbers of cells.

We can use an alternative version of plotClusters called plotClustersWorkflow that will better emphasize the more final clusterings from the ensemble/concensus step and merging steps (it currently does not allow for showing the sampleData as well, however – only clustering results).

par(mar=plotCMar)
plotClustersWorkflow(rsecFluidigm)

Barplots We can examine size distribution of a single clustering with the function plotBarplot. By default, the cluster picked will be the primary cluster, which for the result of clusterMany is rather arbitrarily just the first cluster.

plotBarplot(rsecFluidigm,main="Final Clustering")

We can also pick a particular intermediate clustering, say our intial ensembl clustering before merging.

plotBarplot(rsecFluidigm,whichClusters=c("combineMany" ))

We can also compare two specific clusters with a simple barplot using plotBarplot. Here we compare the “combineMany” and the “mergeClusters” clusterings.

plotBarplot(rsecFluidigm,whichClusters=c("mergeClusters" ,"combineMany"))

Since “combineMany” is a partition of “mergeClusters”, there is perfect subsetting within the clusters of “mergeClusters”.

Co-Clustering We can also visualize the proportion of times samples were together in the individual clusterings (i.e. before the consensus clustering):

plotCoClustering(rsecFluidigm,whichClusters=c("mergeClusters","combineMany"))

Note that this is not the result from any particular subsampling (which was done repeatedly for each clustering), but the proportion of times across the clusterings we ran. The initial consensus clustering in combineMany was made based on these proportions and a particular cutoff of the required proportion of times the samples needed to be together.

Plot of Dendrogram We can visualize how the initial ensembl cluster in combineMany was clustered into a hierarchy and merged to give us the final clustering in mergeClusters:

plotDendrogram(rsecFluidigm,whichClusters=c("combineMany","mergeClusters"),plotType="colorblock",leafType="sample")

As shown in this plot, the individual clusters of the combineMany ensembl clustering were hierarchically clustered (hence the note that the dendrogram was made from the combineMany clustering), and similar sister clusters were merged if there were not enough gene differences between them.

2D plot of clusters Finally, we can plot a 2-dimensional representation of the data with PCA and color code the samples to get a sense of how the data cluster.

plotReducedDims(rsecFluidigm)

We can also look at a higher number of dimensions (or different dimensions) by changing the parameter ‘whichDims’.

plotReducedDims(rsecFluidigm,whichDims=c(1:4))

In this case we can see that higher dimensions show us a greater amount of separation between the clusters than in just 2 dimensions.

par(mar=plotCMar)
plotClusters(ce,main="Clusters from clusterMany", whichClusters="workflow", 
             axisLine=-1)

par(mar=plotCMar)
plotClusters(ce,whichClusters="clusterMany",
               main="Only Clusters from clusterMany",axisLine=-1)

par(mar=plotCMar)
plotClusters(ce,whichClusters="workflow", sampleData=c("Biological_Condition","Published2"), 
               main="Workflow clusters plus other data",axisLine=-1)

5.1.1 Manipulations of colors

In this section we will talk about a number of related ways to manipulate the colors assigned to clusters in conjunction with the plotClusters command.

We will turn off the plotting of the cluster labels to make the plot less cluttered using the option plot=FALSE.

Saving Assignment of colors from 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 into the appropriate slots of the ClusterExperiment object in order to save the color and alignments of samples. This is done 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 (and indeed this is done automatically by mergeClusters so that the combineMany and mergeClusters steps will have aligned color assignments). Similarly, if resetNames=TRUE the names of the clusters will be changed to be integer values, but now those integers will be aligned to try to be the same across clusters (and therefore will not be consecutive integers, which is why these are saved as names for the clusters and not the internal clusterIds). If resetOrderSamples=TRUE, then the order of the samples shown in the plot will be similarly saved in the slot orderSamples.

As an example, we will make a call to plotClusters, but now ask to reset everything to match this clusterAlignment.

First let’s look at what the object’s default colors are for the first two clusterings, accessed by clusterLegend function:

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

## $mergeClusters
##      clusterIds color     name
## [1,] "-1"       "white"   "-1"
## [2,] "1"        "#1F78B4" "M1"
## [3,] "2"        "#33A02C" "M2"
## [4,] "3"        "#FF7F00" "M3"
## [5,] "4"        "#6A3D9A" "M4"
## 
## $`combineMany,final`
##      clusterIds color     name
## [1,] "-1"       "white"   "-1"
## [2,] "1"        "#A6CEE3" "c1"
## [3,] "2"        "#33A02C" "c2"
## [4,] "3"        "#B15928" "c3"
## [5,] "4"        "#1F78B4" "c4"
## [6,] "5"        "#FF7F00" "c5"
## [7,] "6"        "#6A3D9A" "c6"
## [8,] "7"        "#bd18ea" "c7"

The plotClusterLegend creates a quick legend plot of this information for a single clustering:

plotClusterLegend(ce,whichCluster="combineMany,final")

Now, we’ll run plotClusters and save the new colors. We’ll save this as a different object so that this is not a permanent change for the rest of the vignette, though in practice we would usually overwrite the existing ce to save memory on our computer and not have many versions floating around.

ce_temp<-plotClusters(ce,whichClusters="workflow", sampleData=c("Biological_Condition","Published2"), 
               main="Cluster Alignment of Workflow Clusters",clusterLabels=FALSE, axisLine=-1, resetNames=TRUE,resetColors=TRUE,resetOrderSamples=TRUE)

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 values in the clusterLegend). Instead, there are names that “match up” the clusters across the different clusterings.

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

## $mergeClusters
##      clusterIds color     name
## [1,] "-1"       "white"   "-1"
## [2,] "1"        "#E31A1C" "1" 
## [3,] "2"        "#1F78B4" "2" 
## [4,] "3"        "#33A02C" "3" 
## [5,] "4"        "#FF7F00" "4" 
## 
## $`combineMany,final`
##      clusterIds color     name
## [1,] "-1"       "white"   "-1"
## [2,] "1"        "#A6CEE3" "7" 
## [3,] "2"        "#1F78B4" "2" 
## [4,] "3"        "#B15928" "6" 
## [5,] "4"        "#E31A1C" "1" 
## [6,] "5"        "#33A02C" "3" 
## [7,] "6"        "#FF7F00" "4" 
## [8,] "7"        "#bd18ea" "8"

Manual Assignment of colors A similar setting can be that we have colors we want to manually set for a particular cluster, but we want to have the other clusterings get aligned to the colors of that clustering. We can use plotClusters to assign colors to the other clusterings so that the other clusterings inherit the colors of the cluster of interest.

Let’s manually set the colors for the “mergeClusters” clustering. We’ll again create a (new) ce_temp object so again we don’t overwrite the previous colors for the rest of the vignette. Again, the color information is accessed with the clusterLegend command:

ce_temp<-ce
clusterLegend(ce_temp)[["mergeClusters"]]

##      clusterIds color     name
## [1,] "-1"       "white"   "-1"
## [2,] "1"        "#1F78B4" "M1"
## [3,] "2"        "#33A02C" "M2"
## [4,] "3"        "#FF7F00" "M3"
## [5,] "4"        "#6A3D9A" "M4"

We will just assign new colors to the color column. We can also give them new names too.

clusterLegend(ce_temp)[["mergeClusters"]][,"color"]<-c("white","blue","green","cyan","purple")
clusterLegend(ce_temp)[["mergeClusters"]][,"name"]<-c("Not assigned","Cluster1","Cluster2","Cluster3","Cluster4")

We use the plotClusterLegend to check that we assigned them as we expected:

plotClusterLegend(ce_temp,whichCluster="mergeClusters")

Now we will run the plotClusters alignment plot, but we will direct the alignment to use the cluster colors we just gave for the “mergeClusters” cluster. We do this by using the argument existingColors="firstOnly" – but as the argument option indicates it only works if the “mergeClusters” clustering is the first clustering in the plot.

plotClusters(ce_temp,whichClusters="workflow", sampleData=c("Biological_Condition","Published2"), clusterLabels=FALSE,
               main="Clusters from clusterMany, different order",axisLine=-1,existingColors="firstOnly")

This just created the visualization. We can also save the results of this as we did before with resetColors=TRUE, so that now all of our future clusters make use of this color information. We won’t reset the names, however. Note we can avoid making the plot again with the argument plot=FALSE.

ce_temp<-plotClusters(ce_temp,whichClusters="workflow", sampleData=c("Biological_Condition","Published2"), resetColors=TRUE,
               main="Clusters from clusterMany, different order",axisLine=-1, clusterLabels=FALSE, existingColors="firstOnly",plot=FALSE)

Using only the assigned colors Once we start manually changing the colors, we will sometimes want to align our clusters, but use the existing cluster colors that are saved in the object. We can force plotClusters to use all the existing color assignments, rather than create its own, with the argument existingColors="all". 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 alignment produced in the above plotClusters where the clusterings were ordered according to the workflow and making use of the colors we manually assigned to “mergeClusters”. But now we want to plot only the clusters from clusterMany, yet keep the same colors that we just saved so we can compare them. We do this by setting the argument existingColors="all", meaning use all of the existing colors.

par(mfrow=c(1,2))
plotClusters(ce_temp, sampleData=c("Biological_Condition","Published2"),
             whichClusters="workflow", existingColors="all", clusterLabels=FALSE,
             main="All Workflow Clusters, use existing colors",axisLine=-1)

plotClusters(ce_temp, sampleData=c("Biological_Condition","Published2"), 
             existingColors="all", whichClusters="clusterMany", clusterLabels=FALSE,
             main="clusterMany Clusters, use existing colors",
             axisLine=-1)

We see that while the order of the samples has changed (because I have a different set of clusters to align) the colors assigned to each cluster have stayed the same, so I can more easily compare the plots.

Note that the use of existingColors="firstOnly" and existingColors="all" will not give the same color assignments, even if the colors have been previously aligned to be the same unless its the exact same set of clusterings in the same order. This is because with each set of ordered clusterings, the realignment between clusters changes and so the assignment of colors changes. Here we will make a ClusterAlignment plot for each of these three options of the argument to demonstrate the difference.

par(mfrow=c(2,2))
plotClusters(ce_temp, sampleData=c("Biological_Condition","Published2"), 
             existingColors="all", whichClusters="clusterMany", clusterLabels=FALSE,
             main="clusterMany Clusters, use existing colors",
             axisLine=-1)
plotClusters(ce_temp, sampleData=c("Biological_Condition","Published2"), 
               existingColors="firstOnly", whichClusters="clusterMany", clusterLabels=FALSE,
               main="clusterMany Clusters, use existing of first row only",
               axisLine=-1)
plotClusters(ce_temp, sampleData=c("Biological_Condition","Published2"), 
            existingColors="ignore", whichClusters="clusterMany", clusterLabels=FALSE,
            main="clusterMany Clusters, default\n(ignoring assigned colors)",
            axisLine=-1)

Choosing a different set of colors: plotClusters uses the set of colors defined in the variable massivePalette. This is a set of 484 colors. plotClusters draws from this set of colors sequentially as it goes down the clusterings each time it needs a new color. This is obviously a very large list of colors; the reason for such a long list is that plotClusters will error out if there are not enough colors, so we want to have a large list. If we are running RSEC with many clusterings, and each clustering has many clusters, you can quickly run through many. The first 56 colors are a smaller subset of colors saved as a variable bigPalette, and these have been hand-chosen so that the colors are vibrant and not too similar to each other; the order has also been chosen so that more similar colors are not next to each other. This are the colors that are most frequently seen. massivePalette consists of these colors, plus all of the non-grey colors in colors() that are not already contained in bigPalette

We can examine the colors in bigPalette with the command showPalette:

showPalette()

This plot shows us both the name of the color (on the top) and the index of the color in bigPalette.

We can show a smaller subset or give an arbitrary set of colors to show with showPalette

showPalette(which=1:10)

showPalette(palette())

We can even show the entire massivePalette (notice for list of colors > 100 in length, the index is no longer plotted)

showPalette(massivePalette,cex=0.5)

We can use this information to help change our color choices. For example, suppose I want the unassigned samples to be given the color black instead of white. Then I would like to not use the black in bigPalette to assign to a clustering. I will use the colPalette argument to give a new set of colors that does not include “black”. And I will set the unassigned samples using the option unassignedColor="black" (there is a similar argument missingColor to set the color assigned to those clusters with the identification of “-2”, meaning they were not included in the clustering at all).

plotClusters(ce,whichClusters="workflow", sampleData=c("Biological_Condition","Published2"), unassignedColor="black",colPalette=bigPalette[-grep("black",bigPalette)],
               main="Setting unassigned color",axisLine=-1)

5.1.2 Variant version: plotClustersWorkflow

plotClusters is a generic function for showing any clusterings. In following our workflow, however, there can often be many clusterings from clusterMany and the important clusters from combineMany and mergeMany can get lost from view. The function plotClustersWorkflow implements plotClusters but pulls out specific clusterings designated by the user and puts them apart from the results of clusterMany.

plotClustersWorkflow(ce, axisLine= -1)

By default, the “combineMany” and “mergeClusters” clusterings are highlighted. We can choose other clusterings to be the result, rather than the default ones using the argument whichClusters. We can also sort the clusterings differently so that we sort based on the results of clusterMany first, rather than the highlighted results first and/or choose to have the highlighted clusters on the bottom or the top (independently of the sorting choices). In the following code, we choose to show the different versions of combineMany as our highlighted clusters.

par(mar=plotCMar)
plotClustersWorkflow(ce,whichClusters=c("combineMany,final","combineMany,0.7","combineMany,1"),main="Different choices of proportion in combineMany",sortBy="clusterMany",highlightOnTop=FALSE, axisLine= -1)

There are no limits as to which clusterings can be shown with the whichClusters command, but the non-highlighted clusterings must always be of clusterType “clusterMany”. The user can select a subset or different ordering of the “clusterMany” clusterings with the argument whichClusterMany, which must be a numeric vector giving the indices of the clusters to be chosen.

par(mfrow=c(1,1))
par(mar=defaultMar)
plotHeatmap(ce,main="Heatmap with clusterMany")

5.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.

5.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 
## reducedDimNames: PCA 
## filterStats: var var_combineMany.final 
## -----------
## Primary cluster type: mergeClusters 
## Primary cluster label: mergeClusters 
## Table of clusters (of primary clustering):
## -1 M1 M2 M3 M4 
##  6 13 17 15 14 
## 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","Published2"),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.

5.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.2.4 Setting the breaks

Usually, the breaks that determine the colors of the heatmap are evenly spaced across the range of the data in the entire matrix. When there are a few outlier samples or genes, they can dominate the color and make it impossible to visualize the bulk of the data.

For this reason, the argument breaks in plotHeatmap allows for a value between 0 and 1, to indicate that the range of colors should be chosen as equally spaced between certain quantiles of the data. For example, if breaks=0.99, the range of equally spaced breaks will stop at the top 0.99 quantile of the data and anything above that value gets assigned the single extreme color. If there is negative data in the matrix, then it also will use the lower quantile of the data to stop the range of equally spaced breaks (see ?setBreaks)

Here

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

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

The function setBreaks which is called internally by plotHeatmap is also a stand-alone function that the user can call directly to have greater flexibility in getting breaks for the heatmap. For example it allows the user to specify that the breaks should be symmetric around 0. We also provide some default color spectrum that can be better for different settings or symmetric data around 0 – see ?showHeatmapPalettes

6.1.1 Base clustering algorithms and the ClusterFunction class

This package is meant to be able to use and compare different clustering routines. However, the required input, arguments, etc. of different clustering algorithms varies greatly. We create the ClusterFunction class so that we can take ensure that the necessary information to fit into our workflow is well defined, and otherwise the other details of the algorithm can be ignored. In general, the user will not need to know the details of this class, since they will use built-in functions provided by the package which can be accessed by character values. To see the set of character values that correspond to built in functions,

listBuiltInFunctions()

## [1] "pam"            "clara"          "kmeans"         "hierarchical01"
## [5] "hierarchicalK"  "tight"          "spectral"

If you are interested in implementing your own ClusterFunction object see the documentation of the ClusterFunction class.

There are some important features of any clustering algorithm that are encoded in the ClusterFunction object for which it is important to understand because they affect which algorithms can be used when.

inputType The type of input the algorithm expects, which can be either an \(p x n\) matrix of features, in which case the argument x gives that data, or a \(n x n\) matrix of dissimilarities, in which case the argument diss. Some algorithms can accept either type. To determine the inputType of an algorithm(s),

inputType(c("kmeans","pam","hierarchicalK"))

##        kmeans           pam hierarchicalK 
##           "X"      "either"        "diss"

algorithmType we group together algorithms that cluster based on common strategies that affect how we can use them in our workflow. Currently there are two “types” of algorithms we consider, which we call type “K” and “01”. We can determine the type of a builtin function by the following:

algorithmType(c("kmeans","hierarchicalK","hierarchical01"))

##         kmeans  hierarchicalK hierarchical01 
##            "K"            "K"           "01"

The “K” algorithms are so called because their main parameter requirement is that the user specifies the number of clusters (\(K\)) to be created and require an input of k to the clustering function. Built in ‘K’ algorithms are:

listBuiltInTypeK()

## [1] "pam"           "clara"         "kmeans"        "hierarchicalK"
## [5] "spectral"

The “01” algorithms are so named because the algorithm assumes that the input is a disimilarities between samples and that the similarities encoded in \(D\) are on a scale of 0-1. The clustering functions should 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:

listBuiltInType01()

## [1] "hierarchical01" "tight"

**requiredArgs The different algorithm types correspond to requiring different input types (k versus alpha). This is usually sorted out by clusterMany, which will only dispatch the appropriate one. Clustering functions can also have additional required arguments. See below for more discussion about how these arguments can be passed along to clusterMany or RSEC.

To see all of the required arguments of a function,

requiredArgs(c("hierarchical01","hierarchicalK"))

## $hierarchical01
## [1] "alpha"
## 
## $hierarchicalK
## [1] "k"

6.1.2 Internal clustering procedures

clusterMany iteratively calls a function clusterSingle over the collection of parameters. clusterSingle is the clustering workhorse, and may be used by the user who wants more fine-grained control, see documentation of clusterSingle.

Within each call of clusterSingle, there are three possible steps, depending on the value of subsample and sequential. If these are both false, then just a basic clustering routine is done on the input data (called the “main” clustering). If subsample=TRUE, there is first a step that subsamples and clusters the subsamples to calculate a co-occurance matrix, and that is used as the input for the main clustering step. If sequential=TRUE this process is iterated over and over again to iteratively select the best clusters (see ?seqCluster for a detailed description). Each of these steps has a function that goes with it, but that should not generally be called by the user. However, the documentation of these functions can be useful.

In particular, arguments to these functions that are not set by clusterMany can be passed via named lists: subsampleArgs, mainClusterArgs, and seqArgs. Some of the arguments to these steps can be varied in clusterMany, but more esoteric ones should be sent to these arguments of clusterMany (and they will be fixed for parameter combinations tried in clusterMany).

Main Clustering Step (mainClustering) The main clustering step described above is done by the function mainClustering. In addition to the basic clustering algorithms on the input data, 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. Many of the arguments to mainClustering are arguments to clusterMany as well so that mainClusterArgs is usually not needed. The main exception would be to send more esoteric arguments to the underlying clustering function called in the main clustering step. The syntax for this would be to give a nested list to the argument mainClusterArgs

clusterMany(x,clusterFunction="hierarchicalK", ... , mainClusterArgs=list(clusterArgs=list(method="single") )) 

Here we change the argument method in the clustering function hclust called by the hierarchicalK function to single.

Subsampling (subsampleClustering) A more significant processing that can be coupled with any clustering algorithm is to continually by subsample the data and cluster the subsampled data. This creates a \(n x n\) matrix \(S\) that is a matrix of co-clustering percentages – how many times two samples co-clustered together over the subsamples (there are slight variations as how this can be calculated, see help pages of subsampleClustering ). This does not itself give a clustering, but the resulting \(S\) matrix can then form the basis for clustering the samples. Specifically, the matrix \(D=1-S\) is then given as input to the main clustering step described above. 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.

Note that the juxtaposition of these two steps (the subsampling and then feeding the results to the main clustering function) implies there actually two different possible clustering algorithms (and sets of corresponding parameters) – one for the clustering on the subsampled data, and one for the clustering of the resulting \(D\) based on the percentage of coClustering of samples. This brings up a restriction on the clustering function in the main clustering step – it needs to be able to handle input that is a dissimilarity (inputType is either diss or either). Furthermore, the user might want to set clustering function and corresponding parameters separately for the two steps. The way that clusterMany handles this is that the main arguments of clusterMany focus on varying the parameters related to the main clustering step (the clustering of \(D\) after subsampling). For this reason, the argument clusterFunction varies the clustering function used by the main clustering step, not the subsampling step. The clustering function of the subsampling step can be specified by the user via subsampleArgs, but in this case it is set for all calls of clusterMany and does not vary. Alternatively, if the user doesn’t specify the clusterFunction in subsampleArgs then the default is to use clusterFunction of the main clustering step along with any required arguments given by the user for that function (there are some cases where using the clusterFunction of the main step is not possible for the subsampling step, in which case the default is to use “pam”).

More generally, since few of the arguments to subsampleClustering are allowed to be varied by the direct arguments to clusterMany, it is also more common to want to change these arguments via the argument subsampleArgs. Examples might be resamp.num (the number of subsamples to draw) or samp.p (the proportion of samples to draw in each subsample) – see ?subsampleClustering for a full documentation of the possible arguments. In addition, there are arguments to be passed to the underlying clustering function; like for mainClustering, these arguments would be a nested list to the argument subsampleArgs.

An example of a syntax that sets the arguments for subsampleClustering would be:

clusterMany(x,..., subsampleArgs=list(resamp.num=100,samp.p=0.5,clusterFunction="hiearchicalK", clusterArgs=list(method="single") )) 

Sequential Detection of Clusters Another complicated addition that can be added to the main clustering step 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 a clustering 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 when the option sequential=TRUE is chosen (see ?seqCluster for documentation of how the iteration is done). Sequential clustering can also be quite computationally expensive, particularly when paired with subsampling to determine \(D\) at each step of the iteration.

Because of the iterative nature of the sequential step, there are many possible parameters (see ?seqCluster). Like subsample clustering, clusterMany does not allow variation of very many of these parameters, but they can be set via passing arguments in a named list to seqArgs. An example of a syntax that sets the arguments for seqCluster would be:

clusterMany(x,..., seqArgs=list( remain.n=10)) 

This code changes the remain.n option of the sequential step, which governs when the sequential step stops because there are not enough samples remaining.

6.1.3 Arguments of clusterMany

Now that we’ve explained the underlying architecture of the clustering provided in the package, and how to set the arguments that can’t be varied, we discuss the parameters that can be varied in clusterMany. (There are a few additional arguments available for clusterMany that govern how clusterMany works, 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 in the main clustering step. Recall if subsample=TRUE is part of the combination, then clusterFunction the method that will be used on the matrix \(D\) created from subsampling the data. Otherwise, clusterFunction is the clustering method that will be used directly on the data.
• ks The argument ‘ks’ is interpreted differently for different choices of the other parameters and can differ from between parameter combinations!. If sequential=TRUE is part of the parameter combination, ks defines the argument k0 of sequential clustering (see ?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 main clustering step (and by default that of the subsampled data), and is only relevant if clusterFunction is of type “K”. 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 the documentation of clusterMany for more).
• reduceMethod These are character strings indicating what choices of dimensionality reduction should be tried. The choices are “PCA”,“var”,“mad”,“abscv”, and/or “none”. “PCA” indicates clustering on the top principal components. “var”,“mad”,and “abscv” indicate clustering on the top most variable features, as determined by either “var”, “mad” or “abscv” per gene. And “none” indicates the whole data set should be used (which is usually going to be computationally intractable). 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 to try as we have seen):
  • nFilterDims
  • nReducedDims
• 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 standard R function dist, 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 mainClustering 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.
  • 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 of parameters. 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 parameter combinations where findBestK=FALSE the same ks might indicate the specific number of clusters, \(K\), that should be found. To see the parameter choices that will be run, the user can set run=FALSE and the output will be a matrix of the parameter values indicated by the choices of the user. For parameter combinations that are not what is desired, the user should consider making direct calls to clusterSingle where all of these options combinations (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 mainClusterArgs, subsampleArgs and seqArgs. These arguments correspond to the different processes described above (the main clustering step, the creation of \(D\) to be clustered 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.

Argument	Dependencies	Passed to	Argument passed to
ks	sequential=TRUE	seqCluster	k0
-	sequential=FALSE, findBestK=FALSE, clusterFunction of type ‘K’	mainClustering	k
-	sequential=FALSE, findBestK=FALSE, subsample=TRUE	subsampleClustering	k
-	sequential=FALSE, findBestK=TRUE, clusterFunction of type ‘K’	mainClustering	kRange
reduceMethod	none	transform	reduceMethod
nFilterDims	reduceMethod in ‘mad’,‘cv’,‘var’	transform	nFilterDims
nReducedDims	reduceMethod=‘PCA’	transform	nReducedDims
clusterFunction	none	mainClustering	clusterFunction
minSizes	none	mainClustering	minSize
distFunction	subsample=FALSE	mainClustering	distFunction
alphas	clusterFunction of type ‘01’	mainClustering	alpha
findBestK	clusterFunction of type ‘K’	mainClustering	findBestK
removeSil	clusterFunction of type ‘K’	mainClustering	removeSil
silCutoff	clusterFunction of type ‘K’	mainClustering	silCutoff
betas	sequential=TRUE	seqCluster	beta

6.1.4 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. 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 may try to create a work around so that the algorithm just skips algorithms that don’t match the distance). Since the distances we defined are between 0-1, however, we can use any algorithm that takes dissimilarities as input.

Note on 0-1 clustering when subsample=FALSE We would note that the default values of \(\alpha\) in clusterMany and RSEC 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).

To illustrate this point, we show an example of the \(D\) matrix from subsampling. To do this we make use of the clusterSingle which is the workhorse mentioned above that runs a single clustering command directly; it gives the output \(D\) from the sampling in the “coClustering” slot of ce when we set replaceCoCluster=TRUE (and therefore we save it as a separate object, so that it doesn’t write over the existing “coClustering” slot in 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,reduceMethod="mad",nDims=1000,subsample=TRUE,subsampleArgs=list(clusterFunction="pam",clusterArgs=list(k=8)),clusterLabel="subsamplingCluster",mainClusterArgs=list(clusterFunction="hierarchical01",clusterArgs=list(alpha=0.1),minSize=5), replaceCoClustering=TRUE)
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 (at leastwith size greater than 5).

However, the distances based on correlation calculated directly on the data, such as we created above, are also often used for clustering expression data directly (i.e. without the subsampling step). But they are unlikely to have dissimilarities 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(transformData(ce,reduceMethod="mad",nFilterDims=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\).

Since we haven’t yet calculated “mad” on this object, it hasn’t been calculated yet. clusterMany does not let you mix and match between uncalculated and stored filters (or dimensionality reductions), so our first step is to store the mad results. We will save these results as a separate object so as to not disrupt the earlier workflow.

ceDist<-makeFilterStats(ce,filterStats="mad")
ceDist

## class: ClusterExperiment 
## dim: 7069 65 
## reducedDimNames: PCA 
## filterStats: var var_combineMany.final mad 
## -----------
## Primary cluster type: mergeClusters 
## Primary cluster label: mergeClusters 
## Table of clusters (of primary clustering):
## -1 M1 M2 M3 M4 
##  6 13 17 15 14 
## 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

ceDist<-clusterMany(ceDist, 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"),
          reduceMethod=c("mad"),nFilterDims=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 “tight” methods did not give relevant results (no samples were clustered)

6.1.5 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,
                        reduceMethod=c("PCA","var"),
                        nFilterDims=c(100,500,1000),nReducedDims=c(5,15,50),run=FALSE)
dim(checkParam$paramMatrix) #number of rows is the number of clusterings

## [1] 108  14

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).

6.3.1 makeDendrogram and plotDendrogram

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 reduceMethod and the number of dimensions via nDims.

ce<-makeDendrogram(ce,reduceMethod="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 
## reducedDimNames: PCA 
## filterStats: var var_combineMany.final var_combineMany.nVAR 
## -----------
## Primary cluster type: combineMany 
## Primary cluster label: combineMany,nVAR 
## Table of clusters (of primary clustering):
## -1 c1 c2 c3 c4 c5 c6 c7 
## 10  7  3  9 15 13  5  3 
## 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,reduceMethod="var",nDims=500,
                   whichCluster="combineMany,final")

We will visualize the dendrogram with plotDendrogram. The default setting plots the dendrogram where there are color blocks equal to the size of the clusters (i.e number of samples in each cluster).

plotDendrogram(ce,leafType="sample",plotType="colorblock")

We can actually use plotDendrogram to compare clusterings too, like plotClusters using the whichClusters argument to identfy which clusters to show. For example, lets compare our different combineMany results

whCM<-grep("combineMany",clusterTypes(ce))
plotDendrogram(ce,whichClusters=whCM,leafType="sample",plotType="colorblock")

Unlike plotClusters, however, there is no aligning of samples to make samples with the same cluster group together.

Making this clustering current Note that if we look at the clusterType of the “combineMany,final” cluster, it is not of type “combineMany”, but “combineMany.x”. This is because we’ve run additional combineMany steps on this data, and the clustering “combineMany,final” is not the most current one. So to distinguish old iterations from the most recent run, the “.x” is added to the “clusterType” where “x” indicates what iteration it was (it would have added “.x” to the clusterLabel as well as if we had assigned it our own label)

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

## [1] "combineMany.3"

We can choose reset this past call to combineMany to be the current ‘combineMany’ output (which will also set this clustering to be the primaryCluster).

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

## class: ClusterExperiment 
## dim: 7069 65 
## reducedDimNames: PCA 
## filterStats: var var_combineMany.final var_combineMany.nVAR 
## -----------
## Primary cluster type: combineMany 
## Primary cluster label: combineMany,final 
## Table of clusters (of primary clustering):
## -1 c1 c2 c3 c4 c5 c6 c7 
##  6  4  8  5  9 15 14  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

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.

6.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 a given False Discovery Rate threshold of 0.05 (using the Benjamini-Hochberg procedure). However, other more sophisticated 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 
## reducedDimNames: PCA 
## filterStats: var var_combineMany.final var_combineMany.nVAR 
## -----------
## Primary cluster type: combineMany 
## Primary cluster label: combineMany,final 
## Table of clusters (of primary clustering):
## -1 c1 c2 c3 c4 c5 c6 c7 
##  6  4  8  5  9 15 14  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

We see in the summary “Dendrogram run on ‘combineMany,final’”, showing us that this is the clustering that will be used (and also showing us the value of giving our own labels to the results of combineMany if we are going to try different strategies).

We will run mergeClusters with the option mergeMethod="adjP". We will also set plotInfo="adjP" meaning that we would like the mergeClusters command to also produce a plot showing the dendrogram and the estimates from the adjP method for each node. We also set calculateAll=FALSE for illustration purposes, meaning the function will only calculate the estimates for the methods we request, but as we explain below, that is not necessarily the best option if you are going to be trying out different cutoffs.

ce<-mergeClusters(ce,mergeMethod="adjP",plotInfo=c("adjP"),calculateAll=FALSE)

The info about the merge is saved in the ce object.

mergeMethod(ce)

## [1] "adjP"

mergeCutoff(ce)

## [1] 0.05

nodeMergeInfo(ce)

##    Node                     Contrast isMerged mergeClusterId Storey PC
## 1 Node1 (X6+X3+X2+X7)/4-(X5+X1+X4)/3    FALSE             NA     NA NA
## 2 Node2              X6-(X3+X2+X7)/3     TRUE              2     NA NA
## 3 Node3                 X5-(X1+X4)/2     TRUE              1     NA NA
## 4 Node4                        X1-X4     TRUE             NA     NA NA
## 5 Node5                 X3-(X2+X7)/2     TRUE             NA     NA NA
## 6 Node6                        X2-X7     TRUE             NA     NA NA
##          adjP locfdr MB JC
## 1 0.076672797     NA NA NA
## 2 0.027302306     NA NA NA
## 3 0.047672938     NA NA NA
## 4 0.006648748     NA NA NA
## 5 0.007638987     NA NA NA
## 6 0.005941434     NA NA NA

Notice that nodeMergeInfo gives for each node the proportion estimated to be differentially expressed at each node (as displayed in the plot that we requested), as well as whether that node was merged together in the mergeClusters call (the isMerged column). Because we set calculateAll=FALSE only the methods needed for our command were calculated (adjP). The others have NA values. The column mergeClusterId tells us which nodes in the tree are now equivalent to a cluster; this is different than the isMerged column, since some nodes can be merged but if their parent nodes were also merged, then that node will not be equivalent to a cluster in the “mergeClusters” clustering.

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 plotInfo="all", all of the estimates of the different methods are displayed simultaneously, while before we only showed the values for the specific mergeMethod we requested.

ce<-mergeClusters(ce,mergeMethod="none",plotInfo="all")

Notice that now if we call nodeMergeInfo, all of the methods now have estimates.

nodeMergeInfo(ce)

##    Node                     Contrast isMerged mergeClusterId Storey PC
## 1 Node1 (X6+X3+X2+X7)/4-(X5+X1+X4)/3    FALSE             NA     NA NA
## 2 Node2              X6-(X3+X2+X7)/3     TRUE              2     NA NA
## 3 Node3                 X5-(X1+X4)/2     TRUE              1     NA NA
## 4 Node4                        X1-X4     TRUE             NA     NA NA
## 5 Node5                 X3-(X2+X7)/2     TRUE             NA     NA NA
## 6 Node6                        X2-X7     TRUE             NA     NA NA
##          adjP locfdr MB JC
## 1 0.076672797     NA NA NA
## 2 0.027302306     NA NA NA
## 3 0.047672938     NA NA NA
## 4 0.006648748     NA NA NA
## 5 0.007638987     NA NA NA
## 6 0.005941434     NA NA NA

This means in any future calls to mergeClusters there will be no more need for calculations of per-gene significance, which will speed up the calls if you just want to change the cutoff (all of the methods used the same input of per-gene p-values, so recalling them each time is computationally inefficient). In practice, the default is calculateAll=TRUE, meaning all methods are calculated unless the user specifically requests otherwise.

Now we can pick a cutoff and rerun mergeClusters. We’ll give it a label to keep it separate from the previous merge clusters run we had made. Note, we can turn off plotting completely by setting plot=FALSE.

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

## class: ClusterExperiment 
## dim: 7069 65 
## reducedDimNames: PCA 
## filterStats: var var_combineMany.final var_combineMany.nVAR 
## -----------
## Primary cluster type: mergeClusters 
## Primary cluster label: mergeClusters,v2 
## Table of clusters (of primary clustering):
## -1 m1 m2 
##  6 28 31 
## 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

Notice that the nodeMergeInfo has changed, since different nodes were merged, but the estimates per node stay the same.

nodeMergeInfo(ce)

##    Node                     Contrast isMerged mergeClusterId     Storey
## 1 Node1 (X6+X3+X2+X7)/4-(X5+X1+X4)/3    FALSE             NA 0.35125195
## 2 Node2              X6-(X3+X2+X7)/3     TRUE              2 0.23864762
## 3 Node3                 X5-(X1+X4)/2     TRUE              1 0.27090112
## 4 Node4                        X1-X4     TRUE             NA 0.06125336
## 5 Node5                 X3-(X2+X7)/2     TRUE             NA 0.16056019
## 6 Node6                        X2-X7     TRUE             NA 0.09124346
##           PC        adjP     locfdr         MB         JC
## 1 0.29637573 0.076672797 0.06649586 0.32607158 0.07278762
## 2 0.19728126 0.027302306 0.04179466 0.21092092 0.03323834
## 3 0.22651055 0.047672938 0.06028542 0.25491583 0.05813407
## 4 0.05401465 0.006648748 0.02459745 0.03522422 0.02403522
## 5 0.11930898 0.007638987 0.00000000 0.18885274 0.02231499
## 6 0.07081773 0.005941434 0.02688885 0.10468242 0.03750815

If we want to rerun mergeClusters with a different method, we can do that instead.

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

## class: ClusterExperiment 
## dim: 7069 65 
## reducedDimNames: PCA 
## filterStats: var var_combineMany.final var_combineMany.nVAR 
## -----------
## Primary cluster type: mergeClusters 
## Primary cluster label: mergeClusters,v3 
## Table of clusters (of primary clustering):
## -1 m1 m2 m3 m4 m5 
##  6 13 12  5 15 14 
## Total number of clusterings: 44 
## Dendrogram run on 'combineMany,final' (cluster index: 6)
## -----------
## Workflow progress:
## clusterMany run? Yes 
## combineMany run? Yes 
## makeDendrogram run? Yes 
## mergeClusters run? Yes

We can use plotDendrogram to compare the results. Notice that plotDendrogram can recreate the above plots that were created in the calls to mergeClusters via the argument mergeInfo (of course, this only works after mergeClusters has actually been called so that the information is saved in the ce object).

par(mar=c(1.1,1.1,6.1,2.1))
plotDendrogram(ce,whichClusters=c("mergeClusters,v3","mergeClusters,v2"),mergeInfo="mergeMethod")

workflowClusterTable(ce)

##                Iteration
## Type             0  1  2  3  4  5
##   final          0  0  0  0  0  0
##   mergeClusters  1  0  0  1  1  1
##   combineMany    1  1  1  0  1  0
##   clusterMany   36  0  0  0  0  0

head(workflowClusterDetails(ce),8)

##   index          type iteration             label
## 1     1 mergeClusters         0  mergeClusters,v3
## 2     2 mergeClusters         5  mergeClusters,v2
## 3     3 mergeClusters         4   mergeClusters.4
## 4     4   combineMany         4  combineMany,nVAR
## 5     5 mergeClusters         3   mergeClusters.3
## 6     6   combineMany         0 combineMany,final
## 7     7   combineMany         2   combineMany,0.7
## 8     8   combineMany         1     combineMany,1

6.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.

7.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 on the results of rsec:

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

## [1] 10  9

head(pairsAllTop)

##   IndexInOriginal  Contrast Feature    logFC  AveExpr         t      P.Value
## 1            5745 Cl01-Cl02   STMN2 9.839919 8.440126 13.173344 4.850101e-20
## 2            3854 Cl01-Cl02   NRXN1 7.655286 5.583251 10.637778 6.938049e-16
## 3            5085 Cl01-Cl02    RTN1 7.748156 7.266948  9.802468 1.899925e-14
## 4            2325 Cl01-Cl02   GRIA2 7.920474 5.169489  9.587520 4.499016e-14
## 5            3266 Cl01-Cl02    MEG3 7.728490 5.496095  9.047550 3.984190e-13
## 6            4285 Cl01-Cl02  PLXNA2 7.636054 5.379982  9.011605 4.609957e-13
##      adj.P.Val        B
## 1 3.428536e-16 33.03996
## 2 2.452253e-12 24.65642
## 3 4.476856e-11 21.69298
## 4 7.950887e-11 20.91666
## 5 5.431297e-10 18.94485
## 6 5.431297e-10 18.81259

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(rsecFluidigm) 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. “Cl01-Cl02” means a comparison of cluster 1 and cluster 2 (note that these refer to the cluster ids, not any name they might have). 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).

7.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(rsecFluidigm,contrastType="OneAgainstAll",p.value=0.05,number=NROW(rsecFluidigm))
head(best1vsAll)

##   IndexInOriginal ContrastName      Contrast Feature    logFC  AveExpr
## 1            5745         Cl01 Cl01-(Cl02)/1   STMN2 9.839919 8.440126
## 2            3854         Cl01 Cl01-(Cl02)/1   NRXN1 7.655286 5.583251
## 3            5085         Cl01 Cl01-(Cl02)/1    RTN1 7.748156 7.266948
## 4            2325         Cl01 Cl01-(Cl02)/1   GRIA2 7.920474 5.169489
## 5            3266         Cl01 Cl01-(Cl02)/1    MEG3 7.728490 5.496095
## 6            4285         Cl01 Cl01-(Cl02)/1  PLXNA2 7.636054 5.379982
##           t      P.Value    adj.P.Val        B
## 1 13.173344 4.850101e-20 3.428536e-16 33.03996
## 2 10.637778 6.938049e-16 2.452253e-12 24.65642
## 3  9.802468 1.899925e-14 4.476856e-11 21.69298
## 4  9.587520 4.499016e-14 7.950887e-11 20.91666
## 5  9.047550 3.984190e-13 5.431297e-10 18.94485
## 6  9.011605 4.609957e-13 5.431297e-10 18.81259

Notice that now there is both a “Contrast” and a “ContrastName” column, unlike with the pairs comparison. Like before, “Contrast” gives an explicit definition of what is the comparisons, in the form of “(Cl02+Cl03+Cl04+Cl05+Cl06)/5-Cl01”, meaning the mean of the means of clusters 2-6 is compared to the mean of cluster1. Note that the contrasts here are always written in terms of the internal (numeric) cluster id, with an “Cl” in front of the number and a ‘0’ to make the number 2 digits. “ContrastName” interprets this into a more usable name, namely that this contrast can be easily identified as a test of “Cl01” (cluster 1).

We can plot the contrasts with a heatmap for these results. Here we notice that the color next to the gene group matches the cluster that the contrast matches.

plotContrastHeatmap(rsecFluidigm,signifTable=best1vsAll,nBlankLines=10, whichCluster="primary")

7.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.

bestDendro<-getBestFeatures(rsecFluidigm,contrastType="Dendro",p.value=0.05,number=NROW(rsecFluidigm))
head(bestDendro)

##   IndexInOriginal ContrastName Contrast Feature     logFC  AveExpr
## 1            5745        Node1    X2-X1   STMN2 -9.839919 8.440126
## 2            3854        Node1    X2-X1   NRXN1 -7.655286 5.583251
## 3            5085        Node1    X2-X1    RTN1 -7.748156 7.266948
## 4            2325        Node1    X2-X1   GRIA2 -7.920474 5.169489
## 5            3266        Node1    X2-X1    MEG3 -7.728490 5.496095
## 6            4285        Node1    X2-X1  PLXNA2 -7.636054 5.379982
##            t      P.Value    adj.P.Val        B
## 1 -13.173344 4.850101e-20 3.428536e-16 33.03996
## 2 -10.637778 6.938049e-16 2.452253e-12 24.65642
## 3  -9.802468 1.899925e-14 4.476856e-11 21.69298
## 4  -9.587520 4.499016e-14 7.950887e-11 20.91666
## 5  -9.047550 3.984190e-13 5.431297e-10 18.94485
## 6  -9.011605 4.609957e-13 5.431297e-10 18.81259

Again, there is both a “ContrastName” and “Contrast” column. The “Contrast” column identifies which clusters ids were on each side of the node (and hence commpared) and “ContrastName” is the name of the node, determined internally during makeDendrogram.

levels((bestDendro)$Contrast)

## [1] "X2-X1"

We can look at the results again with plotContrastHeatmap.

plotContrastHeatmap(rsecFluidigm,signifTable=bestDendro,nBlankLines=10)

We can plot the dendrogram to help make sense of which contrasts go with which nodes and choose to show the node names with show.node.label=TRUE (plotDendrogram calls plot.phylo from the ape package and can take as imput those arguments like show.node.label).

plotDendrogram(rsecFluidigm,show.node.label=TRUE,whichClusters=c("combineMany","mergeClusters"),leaf="samples",plotType="colorblock")