3D plot genomic tracks

依照大多数人的观点,使用3D表示genomic tracks并不是一个很好的选择,因为会使图象变得复杂难懂,难以把握作者的意图。但是有些时候,为了眩,也不得不做。我试着做了两个版本,都不是很理想。

先准备好数据:

library("mvtnorm")
x <- 1:10
y <- 1:5
dens <- matrix(dmvnorm(expand.grid(x, y), 
                       sigma = rbind(c(3, 2), c(2, 3))), 
               ncol=length(y))
library(reshape2)
z <- melt(dens)
head(z)
##   Var1 Var2        value
## 1    1    1 5.827419e-02
## 2    2    1 3.534508e-02
## 3    3    1 1.176536e-02
## 4    4    1 2.149337e-03
## 5    5    1 2.154900e-04
## 6    6    1 1.185695e-05

首先使用scatterplot3d包来画。它的缺点是不能绕y辆旋转。

library(scatterplot3d)
par(xpd=NA)
s3d <- scatterplot3d(z$Var1, z$Var2, z$value, type="n", 
                     grid=FALSE, box=FALSE, 
                     xlab="coordinate", ylab ="", zlab ="coverage", 
                     angle=55, y.ticklabs=paste0("track", y), 
                     lab=c(5, length(y), 7), lab.z=5)
s3d$polygon3d <- function(x, y = NULL, z = NULL, ...){
    xyz <- xyz.coords(c(0, x, 0), c(y[1], y, y[length(y)]), c(0, z, 0))
    if (angle > 2) {
        temp <- xyz$x
        xyz$x <- xyz$y
        xyz$y <- temp
    }
    y2 <- (xyz$y - y.add)/y.scal
    x <- xyz$x/x.scal + yx.f * y2
    y <- xyz$z/z.scal + yz.f * y2
    mem.par <- par(mar = mar, usr = usr)
    on.exit(par(mem.par))
    polygon(x, y, ...)
}
environment(s3d$polygon3d) <- environment(s3d$points3d)
hiCol <- function(x, alpha="CC"){
  sprintf(paste0("#%s", alpha), 
           paste(as.hexmode(as.numeric(col2rgb(x))), collapse=""))
}
for(i in length(y):1) 
  s3d$polygon3d(1:nrow(dens), rep(i, nrow(dens)), 
                dens[, i], col=hiCol(i+1), border=i+1)

plot of chunk unnamed-chunk-2

下面的代码是使用plot3D来画。它的代码简单多了,但是无法控制坐标文字。

library("plot3D")
#scatter3D(z$Var1, z$Var2, z$value)
colvar <- matrix(rep(y, each=length(x)), nrow=length(x))
ribbon3D(x=x, y=y, z=dens, along = "x", 
         curtain = TRUE, space = 0.8, shade = 0.2, 
         col = 2:6, colvar=colvar, colkey=FALSE, 
         zlab="", xlab="", ylab="", 
         ticktype = "detailed", nticks=4, 
         theta = 30, phi=20)

plot of chunk unnamed-chunk-3

Leave a Reply

  

  

  

%d 博主赞过: