Haematopoetic stem cells from mouse subjects.

This code produces standard R graphics, rather than the RGL graphics as published, since the former are now the default for 3D KDEs in the ks >= 1.12.0 package.

library(ks)
library(colorspace)
library(RColorBrewer)
library(alphashape3d)

## colour functions
pt.col <- function(pos=1, alpha=1, ...) {return(paste0(brewer.pal(12, "Paired"), format(as.hexmode(round(alpha*255,0)), width=2))[pos])}
seq.col2 <- function(n, alpha=1, ...) {c("transparent",tail(rev(diverge.col(2*n-1, alpha=alpha, ...)), n-1))}
seq.col3d <- function(n, alpha=1, ...){rev(sequential_hcl(n=n, h=290, power=1.1, c.=c(80,70), l=c(30,50), alpha=alpha, ...))}

## auxiliary function
plotashape <- function(x, col=1, alpha=1, ...)
{
    triangles <- x$triang
    iAlpha <- 1
    tr <- t(triangles[triangles[,8+iAlpha]==2 | triangles[,8+iAlpha]==3, c("tr1","tr2","tr3")])
    if (length(tr) != 0) 
        plot3D::triangle3D(x$x[tr, c(1,2,3)], col=col, alpha=alpha, ...)     
}

par(oma=c(0,0,0,0)+0.1, mgp=c(1.8,0.5,0), mar=c(0,2,0,0.3)+0.1, cex.axis=1, cex.lab=1)
xlab <- "CD45.1"
ylab <- "CD45.2"
zlab <- "Ly65Mac1"
xlim <- c(0, 1000)
ylim <- c(0, 1000)
zlim <- c(0, 1000)
	
## prepare data
data(hsct)
hsct6 <- hsct[hsct$subject==6, c(1,4,2)]   ## unsuccessful graft
hsct5 <- hsct[hsct$subject==5, c(1,4,2)]   ## unsuccessful graft
hsct9 <- hsct[hsct$subject==9, c(1,4,2)]   ## successful graft
hsct12 <- hsct[hsct$subject==12, c(1,4,2)] ## successful graft

## remove dead cells with zero fluorescence
hsct5 <- as.matrix(hsct5[apply(hsct5>0, 1, all),])
hsct6 <- as.matrix(hsct6[apply(hsct6>0, 1, all),])
hsct9 <- as.matrix(hsct9[apply(hsct9>0, 1, all),])
hsct12 <- as.matrix(hsct12[apply(hsct12>0, 1, all),])	  

## Fig 1.4
## scatter plot
fhat.hsct <- kde(x=hsct12)

plot(fhat.hsct, cont=0, drawpoints=TRUE, col.pt=pt.col(4, alpha=0.1), cex=0.5, pch=16, xlim=xlim, ylim=ylim, zlim=zlim, xlab=xlab, ylab=ylab, zlab=zlab, phi=15, theta=-40)

## Fig 2.12a
## kernel density estimate with normal scale bandwidth 
fhat.hsct.ns <- kde(x=hsct12,  H=Hns(hsct12))

plot(fhat.hsct.ns, col.fun=seq.col3d, xlim=xlim, ylim=ylim, zlim=zlim, xlab=xlab, ylab=ylab, zlab=zlab, phi=15, theta=-40)

## Fig 2.12b
## kernel density estimate with UCV bandwidth 
fhat.hsct.ucv <- kde(x=hsct12,  H=Hucv(jitter(hsct12)))

plot(fhat.hsct.ucv, col.fun=seq.col3d, xlim=xlim, ylim=ylim, zlim=zlim, xlab=xlab, ylab=ylab, zlab=zlab, phi=15, theta=-40)

# Fig 2.12c
## kernel density estimate with plug-in bandwidth
lot(fhat.hsct, col.fun=seq.col3d, xlim=xlim, ylim=ylim, zlim=zlim, xlab=xlab, ylab=ylab, zlab=zlab, phi=15, theta=-40)

## Fig 2.12d
## kernel density estimate with SCV bandwidth
fhat.hsct.scv <- kde(x=hsct12,  H=Hscv(hsct12))

plot(fhat.hsct.scv, col.fun=seq.col3d, xlim=xlim, ylim=ylim, zlim=zlim, xlab=xlab, ylab=ylab, zlab=zlab, phi=15, theta=-40)

## Fig 5.8a
## kernel density curvature estimate with normal scale bandwidth 
fhat2.hsct.ns <- kdde(x=hsct12, deriv.order=2, H=Hns(hsct12, deriv.order=2))
fhat2.hsct.curv.ns <- kcurv(fhat2.hsct.ns)

plot(fhat2.hsct.curv.ns, col.fun=seq.col2, xlim=xlim, ylim=ylim, zlim=zlim, xlab=xlab, ylab=ylab, zlab=zlab, phi=15, theta=-40)

## Fig 5.8b
## kernel density curvature estimate with UCV bandwidth
## warning: Hucv may take hours to execute on hsct12
fhat2.hsct.ucv <- kdde(x=hsct12, deriv.order=2, H=Hucv(jitter(hsct12), deriv.order=2))
fhat2.hsct.curv.ucv <- kcurv(fhat2.hsct.ucv)

plot(fhat2.hsct.curv.ucv, col.fun=seq.col2, xlim=xlim, ylim=ylim, zlim=zlim, xlab=xlab, ylab=ylab, zlab=zlab, phi=15, theta=-40)

## Fig 5.8c
## kernel density curvature estimate with plug-in bandwidth
fhat2.hsct <- kdde(x=hsct12, deriv.order=2)
fhat2.hsct.curv <- kcurv(fhat2.hsct)

plot(fhat2.hsct.curv, col.fun=seq.col2, xlim=xlim, ylim=ylim, zlim=zlim, xlab=xlab, ylab=ylab, zlab=zlab, phi=15, theta=-40)

## Fig 5.8d
## kernel density curvature estimate with SCV bandwidth
fhat2.hsct.scv <- kdde(x=hsct12, deriv.order=2, H=Hscv(hsct12, nstage=1, deriv.order=2))
fhat2.hsct.curv.scv <- kcurv(fhat2.hsct.scv)

plot(fhat2.hsct.curv.scv, col.fun=seq.col2, xlim=xlim, ylim=ylim, zlim=zlim, xlab=xlab, ylab=ylab, zlab=zlab, phi=15, theta=-40)

## Fig 6.2a
## modal regions from kernel density estimate
plot(fhat.hsct, cont=50, col=pt.col(10), alpha=0.3, xlim=xlim, ylim=ylim, zlim=zlim, xlab=xlab, ylab=ylab, zlab=zlab, phi=15, theta=-40)

## Fig 6.2b
## modal regions from kernel density curvature estimate
plot(fhat2.hsct.curv, cont=50, col=pt.col(8), alpha=0.3, xlim=xlim, ylim=ylim, zlim=zlim, xlab=xlab, ylab=ylab, zlab=zlab, phi=15, theta=-40)

## Fig 6.8
## kernel mean shift clustering
ms.hsct <- kms(hsct12, verbose=TRUE)

plot(ms.hsct, display="plot3D", col=pt.col((1:5)*2, alpha=0.2), cex=0.4, pch=16, xlim=xlim, ylim=ylim, zlim=zlim, xlab=xlab, ylab=ylab, zlab=zlab, phi=15, theta=-40)
box3d()	   

## Fig 6.12
## significant modal regions
fs.hsct <- kfs(hsct12, verbose=TRUE)

plot(fs.hsct, splom=FALSE, axes=FALSE, col=pt.col(8), xlim=xlim, ylim=ylim, zlim=zlim, xlab=xlab, ylab=ylab, zlab=zlab, phi=15, theta=-40)
for (i in 5:1) plotashape(alphashape3d::ashape3d(unique(ms.hsct$x[ms.hsct$label==i,]), alpha=100), col=pt.col(2*i), alpha=0.05, add=TRUE)

## Fig 7.1
## signficant density difference test
loc.test.hsct <- kde.local.test(x1=hsct6, x2=hsct12, xmin=c(-50,-50,-50), xmax=c(1200,1200,1200))

plot(loc.test.hsct, col=pt.col(c(10,8)), alphavec=c(0.3,0.3), add.legend=FALSE, xlim=xlim, ylim=ylim, zlim=zlim, xlab=xlab, ylab=ylab, zlab=zlab, phi=15, theta=-40)

## Fig 7.2a
## negative control
## NB: this figure is different from the published one
loc.test.hsct.neg.ctrl <- kde.local.test(x1=hsct6, x2=hsct5, xmin=c(-50,-50,-50), xmax=c(1200,1200,1200))

plot(loc.test.hsct.neg.ctrl, col=pt.col(c(10,8)), alphavec=c(0.3,0.3), add.legend=FALSE, xlim=xlim, ylim=ylim, zlim=zlim, xlab=xlab, ylab=ylab, zlab=zlab, phi=15, theta=-40)

## Fig 7.2b
## positive control
loc.test.hsct.pos.ctrl <- kde.local.test(x1=hsct9, x2=hsct12, xmin=c(-50,-50,-50), xmax=c(1200,1200,1200))

plot(loc.test.hsct.pos.ctrl, col=pt.col(c(10,8)), alphavec=c(0.3,0.3), add.legend=FALSE, xlim=xlim, ylim=ylim, zlim=zlim, xlab=xlab, ylab=ylab, zlab=zlab, phi=15, theta=-40)