meth_plot.R

#!/usr/bin/env R
#
# @package meth_plotter
# @author Anna Diez
# @author Izaskun Mallona

methmax<-1
maxSamples<-100
maxCpG<-100


meth_plot <- function(dd, intervals, SampleOrder, plotType, position, grup, poslabel){
  # recoverign the sample sorting start
  
  if(length(which(table(position)>1))>0) stop("Repeated positions are not allowed")
  if(length(which(dd<0))>0) stop("Negative values not allowed")
  if(length(which(dd>methmax))>0) stop(paste0("Values greater than ",methmax," are not allowed"))

  
  if(SampleOrder == "as-is"){
    dd <- dd[seq(dim(dd)[1],1),]
    grup <- grup[seq(dim(dd)[1],1)]
  }
  
  if(SampleOrder == "by-group"){
    dd <- dd[order(grup,decreasing=TRUE),]
    grup <- grup[order(grup,decreasing=TRUE)]
  }
  
  if(SampleOrder == "by-meth"){
    mdd <- apply(dd,1,mean,na.rm=TRUE)
    dd <- dd[order(mdd),]
    grup <- grup[order(mdd)]
  }
  
  if(SampleOrder == "clust"){
    hc <- hclust(dist(dd))
    dendro <- rev(as.dendrogram(hc))
    dd <- dd[order.dendrogram(dendro),]
    grup <- grup[order.dendrogram(dendro)]
  }
  
  # recoverign the sample sorting end
  
  #################################################################################################
  #################################################################################################
  if (plotType == "proportional" || plotType == "nonproportional"){
    
    if(dim(dd)[2]>maxCpG) stop(paste0("plot type circles only allows a maximum of ", maxCpG," positions."))
    if(dim(dd)[1]>maxSamples) stop(paste0("plot type circles only allows a maximum of ", maxSamples," samples."))   
    
    #################################################################################################
    #################################################################################################
    
    if(plotType=="proportional"){
      xmax <- max(position)+100
      xmin <- min(position)-100
      ymax <- dim(dd)[1]
      par(mai=par()$mai+c(0.5,1,0,0))
      plot(1,type="n",xlim=c(xmin,xmax), ylim=c(1,ymax), axes=FALSE, xlab="", ylab="")
      axis(1, at=position, labels=paste0(poslabel), las=2, srt=45, lwd=0, line=-1)
      axis(2, c(1:ymax), rownames(dd)[1:(dim(dd)[1])], las=1, lwd=0)
      axis(3, at=c(xmin,xmax), labels=c("",""), lwd.ticks=0)
      axis(3, at=position, labels=paste0(position), srt=45, lwd=1)
      abline(h=c(1:ymax), lty=1, col="dimgray")
      
      palette = rev(brewer.pal((dim(intervals)[1]),"RdBu"))
      
      for(j in 1:dim(dd)[1]){
        i = 1:dim(dd)[2]
        colplot = palette[as.numeric(cut(c(0,1,dd[j,i]),breaks=dim(intervals)[1]))[c(-1,-2)]]
        formplot = ifelse(is.na(dd[j,i]),4,21)
        points(position[i], rep(j,times = length(i)), pch=formplot, col="black", bg=colplot, cex=2)
      }
      #legend("right", c(paste0(intervals[,1],"-",intervals[,2]),"Not available"), fill=NULL, border="white", bty="o", pt.cex=1.5, pch=c(rep(21,dim(intervals)[1]),4), pt.bg=c("#FFFFFF",rev(gray.colors((dim(intervals)[1]-2))),"#000000"),bg="white")
      #box(lwd=2)
      par(mai=par()$mai+c(-0.5,-1,0,0))

    }
    
    #################################################################################################
    #################################################################################################
    
    if(plotType=="nonproportional"){
      xmax <- dim(dd)[2]#+((dim(dd)[2]-1)/2)
      xmin <- 1
      ymax <- dim(dd)[1]
      
      par(mai=par()$mai+c(0.5,1,0,0))
      plot(1, type="n", xlim=c(xmin,xmax), ylim=c(1,ymax), axes=FALSE, xlab="", ylab="")
      axis(1, at=c(1:dim(dd)[2]), labels=paste0(position), las=2, srt=45, lwd=0, line=-1)
      axis(2, at=c(1:ymax), rownames(dd)[1:(dim(dd)[1])], las=1, lwd=0)
      axis(3, at=c(xmin,xmax), labels=c("",""), lwd.ticks=0)
      axis(3, at=c(1:dim(dd)[2]), labels=paste0(position), srt=45, lwd=1)
      abline(h=c(1:ymax), lty=1, col="dimgray")
      
      palette = rev(brewer.pal((dim(intervals)[1]),"RdBu"))
      
      for(j in 1:dim(dd)[1]){
        i = 1:dim(dd)[2]
        colplot = palette[as.numeric(cut(c(0,1,dd[j,i]),breaks=dim(intervals)[1]))[c(-1,-2)]]
        formplot = ifelse(is.na(dd[j,i]),4,21)
        points(i, rep(j, times=length(i)), pch=formplot, col="black", bg=colplot, cex=2)
      }
      #legend("right", c(paste0(intervals[,1],"-",intervals[,2]),"Not available"), fill=NULL, border="white", bty="o", pt.cex=1.5, pch=c(rep(21,dim(intervals)[1]),4), pt.bg=c("#FFFFFF",rev(gray.colors((dim(intervals)[1]-2))),"#000000"),bg="white")
      #box(lwd=2)
      par(mai=par()$mai+c(-0.5,-1,0,0))
      
    }
  }
  
  #################################################################################################
  #################################################################################################
  
  if(plotType == "grid"){
    library(lattice)
    
    if(dim(dd)[1]<=dim(intervals)[1]) legendheight <- 1
    if(dim(dd)[1]>dim(intervals)[1]) legendheight <- dim(intervals)[1]/dim(dd)[1]
    
    colkey = rev(brewer.pal((dim(intervals)[1]),"RdBu"))
    
    cuts <- c(intervals[1,1],intervals[1:(dim(intervals)[1]-1),2])-0.000001
    cutfinal<-intervals[dim(intervals)[1],2]
    print(obj1 <- levelplot(t(dd[,order(probePositions)]), at=c(0,intervals[,2]),col.regions=colkey, xlab="", ylab="",
                            scales=list(y=list(labels=paste(grup)),
                                        tck=c(1,0),
                                        x=list(labels=poslabel, rot=90)),
                            border="black"
                            )
          )
  }
  
}