Introduction

Comparing samples defined over a single dimension is a straightforward task that relies on standard, well established methods. Meanwhile distance between samples in high dimensional space remains a largely unexplored field. Available solutions rely on multivariate normal distributions, a condition that is both difficult to check and overlooking key behaviors in biological samples, where populations of interest often correspond to a small proportion (<1%) of all the points that have been measured.

We have developed hilbertSimilarity to address the problem of sample similarity in mass cytometry where samples are measured over up to 100 dimensions, and where each sample contains a slightly different number of points (or cells). Our method first transforms each sample into a probability vector, by dividing each dimension into a fixed number of bins and by associating each cell to a specific multidimensional cube. The proportion of cells in each hypercube is characteristic of a given sample. To characterize an compare samples we use measures from Information Theory, since their interpretation in terms of similarity and complexity is straightforward.

To demonstrate the power of hilbertSimilarity we applied the method to a subset of the bodenmiller et al. dataset, comparing the effect of different stimulations and identifying groups of cells that are significantly affected by different treatments.

Compared to other methods, hilbertSimilarity does not rely on expert-driven gating, or require any hypothesis about the number of populations that are present in the data. This makes it a powerful tool to quickly assess a dataset before using traditional methods, or when populations a not known a priori.

Installation

hilbertSimilarity can be installed using the following command

devtools::install_github(yannabraham/hilbertSimilarity)

Once installed the package can be loaded using the standard library command.

library(hilbertSimilarity)

Loading some test data

We first need data from the bodenmiller package:

library(bodenmiller)
data(refPhenoMat)
data(untreatedPhenoMat)
data(refFuncMat)
data(untreatedFuncMat)
data(refAnnots)
refAnnots$Treatment <- 'reference'
data(untreatedAnnots)
fullAnnots <- rbind(refAnnots[,names(untreatedAnnots)],
                    untreatedAnnots)
fullAnnots$Treatment <- factor(fullAnnots$Treatment)
fullAnnots$Treatment <- relevel(fullAnnots$Treatment,'reference')
refMat <- cbind(refPhenoMat,refFuncMat)
untreatedMat <- cbind(untreatedPhenoMat,
                      untreatedFuncMat)
fullMat <- rbind(refMat,untreatedMat)
fullMat <- apply(fullMat,2,function(x) {
  qx <- quantile(x,c(0.005,0.995))
  x[x<qx[1]] <- qx[1]
  x[x>qx[2]] <- qx[2]
  return(x)
})

In this dataset 51936 cells corresponding to 4 have been measured over 23 channels. Cells have been gated into 14 populations. The percentage of each population per treatment is shown below:

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

## # A tibble: 14 x 5
##    Cells          reference `BCR/FcR-XL` `PMA/Ionomycin` vanadate
##    <fct>              <dbl>        <dbl>           <dbl>    <dbl>
##  1 cd14-hladr-         0.14         0.12            0.17     0.43
##  2 cd14-hladrhigh      0.38         0.12            0.03     0.02
##  3 cd14-hladrmid       1.55         0.98            0.37     0.28
##  4 cd14-surf-         11.4         12.1            10.4     28.9 
##  5 cd14+hladr-         0.08         0.11            0.1      1.44
##  6 cd14+hladrhigh      0.13         0.02            0.03     0.02
##  7 cd14+hladrmid       8.69         5.5             2.04     2.45
##  8 cd14+surf-          0.51         0.5             1.47     2.63
##  9 cd4+               25.3         26.7            31.2      7.96
## 10 cd8+               32.1         33.4            35.1     19   
## 11 dendritic           0.37         0.2             0.19     0.36
## 12 igm-                1.47         2.16            1.24     1.42
## 13 igm+                6.77         5.66            6.81     6.25
## 14 nk                 11.1         12.4            10.8     28.8

The smallest cell population, cd14+hladrhigh, corresponds to 0.058% of the total cells.

To demonstrate the use of hilbertSimilarity we will compare analyze the data at the treatment level.

Defining a grid to compare samples

The make.cut function is used to prepare the grid that will be used to process the data matrix. It requires 2 arguments: the number of cuts (or bins), and the minimum number of cells to consider for a density. The latter depends on the technology and on the decision by the scientist running the analysis. For CyTOF we will require a population to correspond to at least 40 cells.

nbins <- 2
cuts <- make.cut(fullMat,
                 n=nbins+1,
                 count.lim=40)

## CD20 
## IgM 
## CD4 
## CD33 
## HLA.DR 
## CD14 
## CD7 
## CD3 
## CD123 
## pStat1 
## pSlp76 
## pBtk 
## pPlcg2 
## pErk 
## pLat 
## pS6 
## pNFkB 
## pp38 
## pStat5 
## pAkt 
## pSHP2 
## pZap70 
## pStat3

show.cut(cuts)

Applying the grid to the data

Given a dataset and a grid, the do.cut function will assign each cell to a particular bin in every dimension.

cutFullMat <- do.cut(fullMat,cuts,type='combined')

Effectively, each cell is now associated to a voxel, and each voxel is enriched for a particular cell type that corresponds to the the unique combination of dimensions and ranges it corresponds to.

To uniquely identify each occupied voxel, one can compute a Hilbert index over the grid. The after the dimensions have been ordered to better capture the

library(entropy)
miFullMat <- matrix(0,nrow=ncol(fullMat),ncol = ncol(fullMat) )
for (i in seq(ncol(fullMat)-1)) {
  for (j in seq(i+1,ncol(fullMat))) {
    cur.tbl <- table(cutFullMat[,i],cutFullMat[,j])
    nent <- 1-mi.empirical(cur.tbl)/entropy.empirical(cur.tbl)
    miFullMat[i,j] <- miFullMat[j,i] <- nent
  }
}
dimnames(miFullMat) <- list(colnames(fullMat),colnames(fullMat))
hcFullMat <- hclust(as.dist(miFullMat))
plot(hcFullMat)

col.order <- hcFullMat$labels[rev(hcFullMat$order)]
hc <- do.hilbert(cutFullMat[,col.order],nbins)

Visualizing the Hilbert curve

treatment.table <- with(fullAnnots,
                        table(hc,Treatment))
treatment.pc <- apply(treatment.table,2,function(x) x/sum(x))

av <- andrewsProjection(t(treatment.pc),breaks=30)

treatment.proj <- t(treatment.pc) %*% t(av$freq)

library(ggplot2)
library(reshape2)

## 
## Attaching package: 'reshape2'

## The following object is masked from 'package:tidyr':
## 
##     smiths

melt(treatment.proj) %>% 
    rename(AndrewsIndex=Var2) %>% 
    mutate(AndrewsIndex=av$i[AndrewsIndex]) %>% 
    ggplot(aes(x=AndrewsIndex,y=value))+
    geom_line(aes(group=Treatment,color=Treatment))+
    theme_light(base_size=16)+
    theme(axis.title.x=element_blank(),
          axis.text.x=element_blank(),
          axis.ticks.x=element_blank(),
          axis.title.y=element_blank(),
          axis.text.y=element_blank(),
          axis.ticks.y=element_blank())

proj <- hilbertProjection(hc,target = 2)

## Projecting the Hilbert curve to 2 dimensions will require 4096 bins

fullProj <- fullAnnots %>% 
    mutate(HilbertIndex=hc) %>% 
    group_by_at(vars(one_of(colnames(fullAnnots),'HilbertIndex'))) %>% 
    count() %>%
    bind_cols(as.data.frame(proj[.$HilbertIndex+1,]))
fullProjCount <- fullProj %>% 
    ungroup() %>% 
    count(HilbertIndex,Treatment) %>% 
    arrange(desc(n))
kable(head(fullProjCount))

fullProj %>% 
    group_by(Treatment) %>% 
    mutate(PC=n/sum(n)) %>% 
    ggplot(aes(x=V1,y=V2))+
    geom_tile(aes(fill=PC),
              width=24,
              height=24)+
    facet_wrap(~Treatment)+
    scale_fill_gradient(low='grey80',high='blue')+
    theme_light(base_size=16)+
    theme(axis.title.x=element_blank(),
          axis.text.x=element_blank(),
          axis.ticks.x=element_blank(),
          axis.title.y=element_blank(),
          axis.text.y=element_blank(),
          axis.ticks.y=element_blank())

Comparing cells using Hilbert index

With each cell now associated to a specific hilbert index, each sample can be described by the percentage of cells from a given sample that corresponds to a particular index. The resulting table can be visualized as a heat map:

heatmap(log10(treatment.pc),
        scale = 'none',
        Colv = NA,
        Rowv = NA,
        labRow = NA,
        col = colorRampPalette(c('white',blues9))(256),
        margin = c(12,1))

In this experiment, each sample corresponds to a specific cell type and treatment: to compute the distance between samples we use a distance derived from information theory, the Jensen-Shannon distance. This is done through the js.dist function. The resulting distance matrix can be used to compute a hierarchical cluster:

treatment.dist <- js.dist(t(treatment.table))
treatment.hc <- hclust(treatment.dist)

When ordering samples using treatment.hc we see patterns emerging:

heatmap(log10(treatment.pc),
        scale = 'none',
        Colv = as.dendrogram(treatment.hc),
        Rowv = NA,
        labRow = NA,
        col = colorRampPalette(c('white',blues9))(256),
        margin = c(12,2))

Samples corresponding to reference and BCR/FcR-XL cluster together, while PMA/Ionomycin and vanadate samples show strong differences with every other sample.