MVKSA | Haematopoetic stem cells
Mvksa Home   Daily Temp   World Bank  
Earth Quake   Stem Cell   Cardio  
Grevillea   Dumbbell   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)