Microbial community segmentation with R

In my previous post I discussed our recent paper in ISME J, in which we used community structure and flow cytometry data to predict bacterial production.  The insinuation is that if you know community structure, and have the right measure of physiology, you can make a decent prediction of any microbial ecosystem function.  The challenge is that community structure data, which often has hundreds or thousands of dimensions (taxa, OTUs, etc.), is not easily used in straightforward statistical models.   Our workaround is to reduce the community structure data from many dimensions to a single categorical variable represented by a number.  We call this process segmentation.

You could carry out this dimension reduction with pretty much any clustering algorithm; you’re simply grouping samples with like community structure characteristics on the assumption that like communities will have similar ecosystem functions.  We  use the emergent self organizing map (ESOM), a neural network algorithm, because it allows new data to be classified into an existing ESOM.  For example, imagine that you are collecting a continuous time series of microbial community structure data.  You build an ESOM to segment your first few years of data, subsequent samples can be quickly classified into the existing model.  Thus the taxonomic structure, physiological, and ecological characteristics of the segments are stable over time.  There are other benefits to use an ESOM.  One is that with many samples (far more than we had in our study), the ESOM is capable of resolving patterns that many other clustering techniques cannot.

There are many ways to construct an ESOM.  I haven’t tried a Python-based approach, although I’m keen to explore those methods.  For the ISME J paper I used the Kohonen package in R, which has a nice publication that describes some applications and is otherwise reasonably well documented.  To follow this tutorial you can download our abundance table here.  Much of the inspiration, and some of the code for this analysis, follows the (retail) customer segmentation example given here.

For this tutorial you can download a table of the closest estimated genomes and closest completed genomes (analogous to an abundance table) here.  Assuming you’ve downloaded the data into your working directory, fire up Kohonen and build the ESOM.

## Kohonen needs a numeric matrix
edge.norm <- as.matrix(read.csv('community_structure.csv', row.names = 1))

## Load the library
library('kohonen')

## Define a grid.  The bigger the better, but you want many fewer units in the grid
## than samples.  1:5 is a good ballpark, here we are minimal.
som.grid <- somgrid(xdim = 5, ydim=5, topo="hexagonal")

## Now build the ESOM!  It is worth playing with the parameters, though in
## most cases you will want the circular neighborhood and toroidal map structure.
som.model.edges <- som(edge.norm, 
                 grid = som.grid, 
                 rlen = 100,
                 alpha = c(0.05,0.01),
                 keep.data = TRUE,
                 n.hood = "circular",
                 toroidal = T)

Congratulations!  You’ve just constructed your first ESOM.  Pretty easy.  You’ve effectively clustered the samples into the 25 units that define the ESOM.  You can visualize this as such:

plot(som.model.edges, type = 'mapping', pch = 19)

There are the 25 map units, with the toroid split and flattened into 2D.  Each point is a sample (row in the abundance table), positioned in the unit that best reflects its community structure.  I’m not going to go into any depth on the ESOM algorithm, which is quite elegant, but the version implemented in the Kohonen package is based on Euclidean distance.  How well each map unit represents the samples positioned within it is represented by the distance between the map unit and each sample.  This can be visualized with:

plot(som.model.edges, type = 'quality', pch = 19, palette.name = topo.colors)

Units with shorter distances in the plot above are better defined by the samples in those units than units with long distances.  What distance is good enough depends on your data and objectives.

The next piece is trickier because there’s a bit of an art to it.  At this point each sample has been assigned to one of the 25 units in the map.  In theory we could call each map unit a “segment” and stop here.  It’s beneficial however, to do an additional round of clustering on the map units themselves.  Particularly on large maps (which clearly this is not) this will highlight major structural features in the data.  Both k-means and hierarchical clustering work fairly well, anecdotally k-means seems to work better with smaller maps and hierarchical with larger maps, but you should evaluate for your data.  Here we’ll use k-means.  K-means requires that you specify the number of clusters in advance, which is always a fun chicken and egg problem.  To solve it we use the within-clusters sum of squares method:

wss.edges <- (nrow(som.model.edges$codes)-1)*sum(apply(som.model.edges$codes,2,var)) 
for (i in 2:15) {
  wss.edges[i] <- sum(kmeans(som.model.edges$codes, centers=i)$withinss)
}

plot(wss.edges,
     pch = 19,
     ylab = 'Within-clusters sum of squares',
     xlab = 'K')

Here’s where the art comes in.  Squint at the plot and try to decide the inflection point.  I’d call it 8, but you should experiment with whatever number you pick to see if it makes sense downstream.

We can make another plot of the map showing which map units belong to which clusters:

k <- 8
som.cluster.edges <- kmeans(som.model.edges$codes, centers = k)

plot(som.model.edges,
     main = '',
     type = "property",
     property = som.cluster.edges$cluster,
     palette.name = topo.colors)
add.cluster.boundaries(som.model.edges, som.cluster.edges$cluster)

Remember that the real shape of this map is a toroid and not a square.  The colors represent the final “community segmentation”; the samples belong to map units, and the units belong to clusters.  In our paper we termed these clusters “modes” to highlight the fact that there are real ecological properties associated with them, and that (unlike clusters) they support classification.  To get the mode of each sample we need to index the sample-unit assignments against the unit-cluster assignments.  It’s a little weird until you get your head wrapped around it:

som.cluster.edges$cluster[som.model.edges$unit.classif]
[1] 5 7 7 5 2 7 5 3 7 5 2 6 1 1 1 7 5 4 7 7 5 7 7 7 7 7 7 1 4 4 4 4 7 7 7 6 6 6 6 1 1 1 7 5 5 5 1 1 1 5 5 7 7 4 8 7 7 4 7 8
[61] 7 7 7 7 6 5 6 7 7 7 6 4 6 5 4 4 6 2 1 1 1 1 1 4 1 4 4 4

A really important thing to appreciate about these modes is that they are not ordered or continuous.  Mode 4 doesn’t necessarily have more in common with mode 5 say, than with mode 1.  For this reason it is important to treat the modes as factors in any downstream analysis (e.g. in linear modeling).  For our analysis I had a dataframe with bacterial production, chlorophyll concentration, and bacterial abundance, and predicted genomic parameters from paprica.  By saving the mode data as a new variable in the dataframe, and converting the dataframe to a zoo timeseries, it was possible to visualize the occurrence of modes, model the data, and test the pattern of modes for evidence of succession.  Happy segmenting!