R/plotMotifStackWithRadialPhylog.R

Defines functions plotMotifStackWithRadialPhylog

Documented in plotMotifStackWithRadialPhylog

###############################################################################
######## radial style stack
######## 
###############################################################################
plotMotifStackWithRadialPhylog <- function (phylog, pfms=NULL,
                                            circle=.75, circle.motif=NA, cleaves=1, cnodes=0,
                                            labels.leaves=names(phylog$leaves), clabel.leaves=1,
                                            labels.nodes=names(phylog$nodes), clabel.nodes=0,
                                            draw.box=FALSE,
                                            col.leaves=rep("black", length(labels.leaves)),
                                            col.leaves.bg=NULL, col.leaves.bg.alpha=1,
                                            col.bg=NULL, col.bg.alpha=1,
                                            col.inner.label.circle=NULL, inner.label.circle.width="default",
                                            col.outer.label.circle=NULL, outer.label.circle.width="default",
                                            clockwise =FALSE, init.angle=if(clockwise) 90 else 0,
                                            angle=360, pfmNameSpliter=";", rcpostfix="(RC)", 
                                            motifScale=c("linear","logarithmic"), ic.scale=TRUE,
                                            plotIndex=FALSE, IndexCol="black", IndexCex=.8,
                                            groupDistance=NA, groupDistanceLineCol="red", 
                                            plotAxis=FALSE, font="Helvetica-Bold", fontsize=12)
{
  if (!inherits(phylog, "phylog"))
    stop("Non convenient data")
  leaves.number <- length(phylog$leaves)
  checkLength <- function(tobechecked){
    !((length(tobechecked)>=leaves.number)||is.null(tobechecked))
  }
  checkNA <- function(tobechecked){
    if(is.null(tobechecked)) return(FALSE)
    return(any(is.na(tobechecked)))
  }
  for(tobechecked in c("col.leaves", "col.leaves.bg", "col.bg", "col.inner.label.circle", "col.outer.label.circle")){
    if(checkLength(eval(as.symbol(tobechecked)))) stop(paste("the length of", tobechecked, "should be same as the length of leaves"))
    if(checkNA(eval(as.symbol(tobechecked)))) stop(paste("contain NA in", tobechecked))
  }
  motifScale <- match.arg(motifScale)
  leaves.names <- names(phylog$leaves)
  nodes.number <- length(phylog$nodes)
  nodes.names <- names(phylog$nodes)
  if (length(labels.leaves) != leaves.number)
    labels.leaves <- names(phylog$leaves)
  if (length(labels.nodes) != nodes.number)
    labels.nodes <- names(phylog$nodes)
  if (circle < 0)
    stop("'circle': non convenient value")
  leaves.car <- gsub("[_]", " ", labels.leaves)
  nodes.car <- gsub("[_]", " ", labels.nodes)
  opar <- par(mar = par("mar"), srt = par("srt"))
  on.exit(par(opar))
  par(mar = c(0.1, 0.1, 0.1, 0.1), mfrow=c(1,1))
  dis <- phylog$droot
  max_Dis <- max(dis)
  axis_pos <- circle
  if(!is.na(groupDistance)){
    groupDistance <- (max_Dis - groupDistance) * circle / max_Dis
  }
  dis <- dis/max_Dis
  rayon <- circle
  dis <- dis * rayon
  dist.leaves <- dis[leaves.names]
  dist.nodes <- dis[nodes.names]
  asp <- c(1, 1)
  if(is.null(pfms)){
    plot.default(0, 0, type = "n", asp = 1, xlab = "", ylab = "",
                 xaxt = "n", yaxt = "n", xlim = c(-2, 2), ylim = c(-2, 2),
                 xaxs = "i", yaxs = "i", frame.plot = FALSE)
  }else{
    pin <- dev.size("in") #pin <- par("pin")
    if (pin[1L] > pin[2L]) asp <- c(pin[2L]/pin[1L], 1)
    else asp <- c(1, pin[1L]/pin[2L])
    plot.default(0, 0, type = "n", asp=1, xlab = "", ylab = "",
                 xaxt = "n", yaxt = "n", xlim = c(-2.5, 2.5), ylim = c(-2.5, 2.5),
                 xaxs = "i", yaxs = "i", frame.plot = FALSE)
  }
  d.rayon <- rayon/(nodes.number - 1)
  twopi <- if (clockwise) -2 * pi else 2 * pi
  alpha <- twopi * angle * (1:leaves.number)/leaves.number/360 + init.angle * pi/180
  names(alpha) <- leaves.names
  x <- dist.leaves * cos(alpha)
  y <- dist.leaves * sin(alpha)
  xcar <- (rayon + d.rayon) * cos(alpha)
  ycar <- (rayon + d.rayon) * sin(alpha)
  
  rayonWidth <- max(unlist(lapply(leaves.names, strwidth, units="user", cex=clabel.leaves)))
  circle.motif <- ifelse(is.na(circle.motif), rayon + d.rayon + rayonWidth, circle.motif)
  ##for logos position
  maxvpwidth <- 0
  if(!is.null(pfms)){
    beta <- alpha * 180 / pi
    vpheight <- 2 * pi * angle * circle.motif / 360 / leaves.number / 5
    vpheight <- vpheight * asp[2L]
    xm <- circle.motif * cos(alpha) * asp[1L] / 5 + 0.5
    ym <- circle.motif * sin(alpha) * asp[2L] / 5 + 0.5
    pfmNamesLen <- sapply(pfms, function(.ele) 
      length(strsplit(.ele@name, pfmNameSpliter)[[1]]))
    if(motifScale=="linear")
      vph <- 2.5*vpheight*pfmNamesLen
    else vph <- 2.5*vpheight*(1+log2(pfmNamesLen+0.0001))
    maxvpwidth <- max(mapply(function(.ele, f) ncol(.ele@mat)*f, pfms, vph))
  }
  if(inner.label.circle.width=="default") inner.label.circle.width <- rayonWidth/10
  if(outer.label.circle.width=="default") outer.label.circle.width <- rayonWidth/10
  
  ratio <- if(!is.null(pfms)) 2.5/(circle.motif+maxvpwidth+outer.label.circle.width) else 1
  if(ratio < 1){
    x <- x * ratio
    y <- y * ratio
    xcar <- xcar * ratio
    ycar <- ycar * ratio
    dis <- dis * ratio
    groupDistance <- groupDistance * ratio
    xm <- (xm-.5) * ratio + .5
    ym <- (ym-.5) * ratio + .5
    vpheight <- vpheight * ratio
    axis_pos <- axis_pos * ratio
  }
  ##for plot background
  if(!is.null(col.bg)) col.bg <- highlightCol(col.bg, col.bg.alpha)
  if(!is.null(col.leaves.bg)) col.leaves.bg <- highlightCol(col.leaves.bg, col.leaves.bg.alpha)
  gamma <- twopi * angle * ((1:(leaves.number+1))-0.5)/leaves.number/360 + init.angle * pi/180
  n <- max(2, floor(200*360/leaves.number))
  plotBgArc <- function(r,bgcol,inr){
    if(ratio < 1){
      r <- r*ratio
      inr <- inr*ratio
    }
    t2xy <- function(rx,t) list(x=rx*cos(t), y=rx*sin(t))
    oldcol <- bgcol[1]
    start <- 1
    icnt <- 1
    for(i in 1:leaves.number){
      oldcol <- bgcol[i]
      if(i==leaves.number || bgcol[i+1]!=oldcol){
        P <- t2xy(r, seq.int(gamma[start], gamma[start+icnt], length.out=n*icnt))
        polygon(c(P$x, 0), c(P$y, 0), border=bgcol[i], col=bgcol[i])
        start <- i+1
        icnt <- 1
      } else {
        icnt <- icnt + 1
      }
    }
    if(inr!=0){
      P <- t2xy(inr, seq.int(0, twopi, length.out=n*leaves.number))
      polygon(c(P$x, 0), c(P$y, 0), border="white", col="white")
    }
  }
  if (clabel.leaves > 0) {
    if(!is.null(col.outer.label.circle)) ##plot outer.label.circle
      plotBgArc(circle.motif+maxvpwidth+outer.label.circle.width, col.outer.label.circle, circle.motif+maxvpwidth)
    if(!is.null(col.inner.label.circle)) #plot inner.label.circle
      plotBgArc(circle.motif, col.inner.label.circle, circle.motif - inner.label.circle.width)
    if(!is.null(col.leaves.bg)) ##plot leaves bg
      plotBgArc(circle.motif - inner.label.circle.width, col.leaves.bg, rayon+d.rayon)
    if(!is.null(col.bg)) ##plot center bg
      plotBgArc(ifelse(mean(dist.leaves)/max(dis) > .9, mean(dist.leaves), max(dis)-d.rayon), col.bg, 0)##plotBgArc(mean(dist.leaves), col.bg, 0)
    for(i in 1:leaves.number) {
      par(srt = alpha[i] * 180/pi)
      text(xcar[i], ycar[i], leaves.car[i], adj = 0, col=col.leaves[i], cex = par("cex") *
             clabel.leaves)
      segments(xcar[i], ycar[i], x[i], y[i], col = grey(0.7))
    }
    
    assign("tmp_motifStack_symbolsCache", list(), envir=.globals)
    if(!is.null(pfms)){
      ##extract names
      for(metaChar in c("\\","$","*","+",".","?","[","]","^","{","}","|","(",")"))
      {
        rcpostfix <- gsub(metaChar,paste("\\",metaChar,sep=""),rcpostfix,fixed=TRUE)
      }
      pfmNames <- lapply(pfms, function(.ele) .ele@name)
      for(i in 1:length(pfmNames)){
        pfmname <- unlist(strsplit(pfmNames[[i]], pfmNameSpliter))
        pfmname <- gsub(paste(rcpostfix,"$",sep=""),"",pfmname)
        pfmIdx <- which(makeLeaveNames(labels.leaves) %in% makeLeaveNames(pfmname))
        if(length(pfmIdx)==0) 
          pfmIdx <- which(makeLeaveNames(names(phylog$leaves)) 
                          %in% makeLeaveNames(pfmname))
        if(length(pfmIdx)>0){
          vph <- ifelse(motifScale=="linear",
                        vpheight*length(pfmname),
                        vpheight*(1+log2(length(pfmname))))
          vpw <- vph * ncol(pfms[[i]]@mat) / 2
          vpd <- sqrt(vph*vph+vpw*vpw) / 2
          angle <- median(beta[pfmIdx])
          if(length(pfmIdx)%%2==1){
            this.pfmIdx <- which(beta[pfmIdx] == angle)[1]
            vpx <- xm[pfmIdx[this.pfmIdx]] + vpd * cos(alpha[pfmIdx[this.pfmIdx]]) * asp[1L]
            vpy <- ym[pfmIdx[this.pfmIdx]] + vpd * sin(alpha[pfmIdx[this.pfmIdx]]) * asp[2L]
            vpx1 <- xm[pfmIdx[this.pfmIdx]] - inner.label.circle.width * cos(alpha[pfmIdx[this.pfmIdx]]) * asp[1L] *1.1/5
            vpy1 <- ym[pfmIdx[this.pfmIdx]] - inner.label.circle.width * sin(alpha[pfmIdx[this.pfmIdx]]) * asp[2L] *1.1/5
          }else{
            this.pfmIdx <- order(abs(beta[pfmIdx] - angle))[1:2]
            vpx <- median(xm[pfmIdx[this.pfmIdx]]) + vpd * cos(median(alpha[pfmIdx[this.pfmIdx]])) * asp[1L]
            vpy <- median(ym[pfmIdx[this.pfmIdx]]) + vpd * sin(median(alpha[pfmIdx[this.pfmIdx]])) * asp[2L]
            vpx1 <- median(xm[pfmIdx[this.pfmIdx]]) - inner.label.circle.width * cos(median(alpha[pfmIdx[this.pfmIdx]])) * asp[1L] *1.1/5
            vpy1 <- median(ym[pfmIdx[this.pfmIdx]]) - inner.label.circle.width * sin(median(alpha[pfmIdx[this.pfmIdx]])) * asp[2L] *1.1/5
          }
          pushViewport(viewport(x=vpx, y=vpy, width=vpw, height=vph, angle=angle))
          if(class(pfms[[i]])!="psam"){
            plotMotifLogoA(pfms[[i]], font=font, ic.scale=ic.scale, fontsize=fontsize)
          }else{
            (pfms[[i]], font=font, fontsize=fontsize, newpage=FALSE)
          }
          popViewport()
          if(plotIndex) {
            grid.text(label=i, x=vpx1, 
                      y=vpy1, 
                      gp=gpar(col=IndexCol, cex=IndexCex), rot=angle, just="right")
          }
        }else{
          warning(paste("No leave named as ", paste(pfmname, collapse=", ")), sep="")
        }
      }
    }
    rm(list="tmp_motifStack_symbolsCache", envir=.globals)
  }
  if (cleaves > 0) {
    for (i in 1:leaves.number) points(x[i], y[i], pch = 21, col=col.leaves[i],
                                      bg = col.leaves[i], cex = par("cex") * cleaves)
  }
  ang <- rep(0, length(dist.nodes))
  names(ang) <- names(dist.nodes)
  ang <- c(alpha, ang)
  for (i in 1:length(phylog$parts)) {
    w <- phylog$parts[[i]]
    but <- names(phylog$parts)[i]
    ang[but] <- mean(ang[w])
    b <- range(ang[w])
    a.seq <- c(seq(b[1], b[2], by = pi/180), b[2])
    lines(dis[but] * cos(a.seq), dis[but] * sin(a.seq), col="#222222", lwd=par("cex") * clabel.leaves)
    x1 <- dis[w] * cos(ang[w])
    y1 <- dis[w] * sin(ang[w])
    x2 <- dis[but] * cos(ang[w])
    y2 <- dis[but] * sin(ang[w])
    segments(x1, y1, x2, y2, col="#222222", lwd=par("cex") * clabel.leaves)
  } 
  if(!is.na(groupDistance)){
    if(length(groupDistanceLineCol)!=length(groupDistance)){
      groupDistanceLineCol <- rep(groupDistanceLineCol, ceiling(length(groupDistance)/length(groupDistanceLineCol)))[1:length(groupDistance)]
    }
    for(i in 1:length(groupDistance)){
      symbols(x=0, y=0, circles=groupDistance[i], fg=groupDistanceLineCol[i], lty=2, inches=FALSE, add=TRUE)
    }
  }
  if (cnodes > 0) {
    for (i in 1:length(phylog$parts)) {
      w <- phylog$parts[[i]]
      but <- names(phylog$parts)[i]
      ang[but] <- mean(ang[w])
      points(dis[but] * cos(ang[but]), dis[but] * sin(ang[but]),
             pch = 21, bg = "white", cex = par("cex") * cnodes)
    }
  }
  
  ## draw distance indix  
  if(plotAxis){
    wd <- if(!is.null(pfms)) 5 else 4
    if(clockwise){
      vp <- viewport(x=0.5, y=0.5, width=axis_pos/wd, height=.1, 
                     xscale=c(0, max_Dis), angle=init.angle, just=c(0, 0))
      pushViewport(vp)
      grid.xaxis(gp=gpar(cex = par("cex") * clabel.leaves, col="lightgray"),
                 main=TRUE)
      popViewport()
    }else{
      vp <- viewport(x=0.5, y=0.5, width=axis_pos/wd, height=.1, 
                     xscale=c(0, max_Dis), angle=init.angle, just=c(0, 1))
      pushViewport(vp)
      grid.xaxis(gp=gpar(cex = par("cex") * clabel.leaves, col="lightgray"),
                 main=FALSE)
      popViewport()
    }
  }
  
  points(0, 0, pch = 21, cex = par("cex") * 2 * clabel.leaves, bg = "red")
  if (clabel.nodes > 0) {
    delta <- strwidth(as.character(length(dist.nodes)), cex = par("cex") *
                        clabel.nodes)
    for (j in 1:length(dist.nodes)) {
      i <- names(dist.nodes)[j]
      par(srt = (ang[i] * 360/2/pi + 90))
      x1 <- dis[i] * cos(ang[i])
      y1 <- dis[i] * sin(ang[i])
      symbols(x1, y1, delta, bg = "white", add = TRUE,
              inches = FALSE)
      text(x1, y1, nodes.car[j], adj = 0.5, cex = par("cex") *
             clabel.nodes)
    }
  }
  
  if (draw.box)
    box()
  return(invisible())
}

Try the motifStack package in your browser

Any scripts or data that you put into this service are public.

motifStack documentation built on Nov. 1, 2018, 4:27 a.m.