mvstat.net MVKSA Book GPS Trajectory Contact
Daily
Temp World
Bank Earth
Quake Stem
Cell Foetal
Cardio Grevillea
  Dumb
Bell Air
Quality

Haematopoetic stem cells from mouse subjects. 
library(ks)
library(colorspace)
library(RColorBrewer)
library(rgl)
library(misc3d)
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, ...))}

## rgl plot parameters
hsct.rgl <-function(xlab, ylab, zlab)
{
    par3d(userMatrix=matrix(c(0.78,0.20,-0.59,0.00,-0.63,0.22,-0.75,0.00,-0.02,0.96,0.29,0,0,0,0,1), nrow=4), zoom=1.1, windowRect=c(300,300,1000, 1000), cex=1.5)
    axes3d(edges=c("x--", "y--"), nticks=5, labels=c(0,200,400,600,800,1000))
    axes3d(edges=c("z-+"), nticks=5, labels=c(0,"",400,"",800,""))
    if (!missing(xlab)) mtext3d(xlab, edge="x--", line=2)
    if (!missing(ylab)) mtext3d(ylab, edge="y--", line=2)
    if (!missing(zlab)) mtext3d(zlab, edge="z-+", line=1)
    aspect3d(c(1,1,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)

clear3d()
plot(fhat.hsct, display="rgl", cont=0, drawpoints=TRUE, col.pt=pt.col(4), alpha=0.3, size=8, axes=FALSE, asp=1, xlim=xlim, ylim=ylim, zlim=zlim, xlab=xlab, ylab=ylab, zlab="")
hsct.rgl(zlab=zlab)

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

clear3d()
plot(fhat.hsct.ns, display="rgl", col.fun=seq.col3d, axes=FALSE, asp=1, xlim=xlim, ylim=ylim, zlim=zlim, xlab="", ylab="", zlab="")
hsct.rgl(xlab=xlab, ylab=ylab, zlab=zlab)

## Fig 2.12b
## kernel density estimate with UCV bandwidth 
## warning: Hucv may take hours to execute on hsct12
fhat.hsct.ucv <- kde(x=hsct12,  H=Hucv(jitter(hsct12)))

clear3d()
plot(fhat.hsct.ucv, display="rgl", col.fun=seq.col3d, xlim=xlim, ylim=ylim, zlim=zlim, xlab="", ylab="", zlab="")
hsct.rgl(xlab=xlab, ylab=ylab, zlab=zlab)

# Fig 2.12c
## kernel density estimate with plug-in bandwidth
clear3d()
plot(fhat.hsct, display="rgl", col.fun=seq.col3d, axes=FALSE, xlim=xlim, ylim=ylim, zlim=zlim, xlab="", ylab="", zlab="")
hsct.rgl(xlab=xlab, ylab=ylab, zlab=zlab)

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

clear3d()
plot(fhat.hsct.scv, display="rgl", col.fun=seq.col3d, axes=FALSE, xlim=xlim, ylim=ylim, zlim=zlim, xlab="", ylab="", zlab="")
hsct.rgl(xlab=xlab, ylab=ylab, zlab=zlab)

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

clear3d()
plot(fhat2.hsct.curv.ns, display="rgl", col.fun=seq.col2, xlim=xlim, ylim=ylim, zlim=zlim, xlab="", ylab="", zlab="", axes=FALSE)
hsct.rgl(xlab=xlab, ylab=ylab, zlab=zlab)

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

clear3d()
plot(fhat2.hsct.curv.ucv, display="rgl", col.fun=seq.col2, axes=FALSE, xlim=xlim, ylim=ylim, zlim=zlim, xlab="", ylab="", zlab="")
hsct.rgl(xlab=xlab, ylab=ylab, zlab=zlab)

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

clear3d()
plot(fhat2.hsct.curv, display="rgl", col.fun=seq.col2, axes=FALSE, xlim=xlim, ylim=ylim, zlim=zlim, xlab="", ylab="", zlab="")
hsct.rgl(xlab=xlab, ylab=ylab, zlab=zlab)

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

clear3d()
plot(fhat2.hsct.curv.scv, display="rgl", col.fun=seq.col2, axes=FALSE, xlim=xlim, ylim=ylim, zlim=zlim, xlab="", ylab="", zlab="")
hsct.rgl(xlab=xlab, ylab=ylab, zlab=zlab)

## Fig 6.2a
## modal regions from kernel density estimate
clear3d()
plot(fhat.hsct, display="rgl", cont=50, color=pt.col(8), axes=FALSE, alphavec=0.4, asp=1, xlim=xlim, ylim=ylim, zlim=zlim, xlab=xlab, ylab=ylab, zlab="")
hsct.rgl(zlab=zlab)

## Fig 6.2b
## modal regions from kernel density curvature estimate
clear3d()
plot(fhat2.hsct.curv, display="rgl", cont=50, color=pt.col(8), axes=FALSE, alphavec=0.4, asp=1, xlim=xlim, ylim=ylim, zlim=zlim, xlab=xlab, ylab=ylab, zlab="")
hsct.rgl(zlab=zlab)

## Fig 6.8
## kernel mean shift clustering
## warning:  warning: kms may take hours to execute on hsct12
ms.hsct <- kms(hsct12, verbose=TRUE)

clear3d()
i <- 1; rgl::plot3d(ms.hsct$x[ms.hsct$label==i,1], ms.hsct$x[ms.hsct$label==i,2], ms.hsct$x[ms.hsct$label==i,3], axes=FALSE, size=8, col=pt.col(i*2), alpha=0.3, asp=1, xlim=xlim, ylim=ylim, zlim=zlim, xlab=xlab, ylab=ylab)
for (i in 2:5) rgl::points3d(ms.hsct$x[ms.hsct$label==i,1], ms.hsct$x[ms.hsct$label==i,2], ms.hsct$x[ms.hsct$label==i,3], size=8, col=pt.col(i*2), alpha=0.3)
hsct.rgl(zlab=zlab)
box3d()	   

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

clear3d()
plot(fs.hsct, display="rgl", axes=FALSE, color=pt.col(8), asp=1, xlim=xlim, ylim=ylim, zlim=zlim, xlab=xlab, ylab=ylab, zlab="")
hsct.rgl(zlab=zlab)
box3d()
for (i in 1:5) plot(alphashape3d::ashape3d(unique(ms.hsct$x[ms.hsct$label==i,]), alpha=100), transparency=0.1, edges=FALSE, lines=FALSE, clear=FALSE, col=pt.col(2*i))

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

clear3d()
plot(loc.test.hsct, display="rgl", add.legend=FALSE, axes=FALSE, col=pt.col(c(10,8)), alpha=0.3, asp=1, xlim=xlim, ylim=ylim, zlim=zlim, xlab="", ylab="", zlab="")
hsct.rgl(xlab=xlab, ylab=ylab, zlab=zlab)

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

clear3d()
plot(loc.test.hsct.neg.ctrl, display="rgl", add.legend=FALSE, axes=FALSE, col=pt.col(c(10,8)), alpha=0.3, asp=1, xlim=xlim, ylim=ylim, zlim=zlim, xlab="", ylab="", zlab="")
hsct.rgl(xlab=xlab, ylab=ylab, zlab=zlab)

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

clear3d()
plot(loc.test.hsct.pos.ctrl, display="rgl", add.legend=FALSE, axes=FALSE, col=pt.col(c(10,8)), alpha=0.3, asp=1, xlim=xlim, ylim=ylim, zlim=zlim, xlab="", ylab="", zlab="")
hsct.rgl(xlab=xlab, ylab=ylab, zlab=zlab)