R/chromplot-Internal.R

In chromPlot: Global visualization tool of genomic data

## Track management functions
## ----- ---------- ---------

build.track <- function(datatab, bycol="Chrom", poscol=c("Mid", "Start"),
start.col="Start", end.col="End", category.col=NULL, binsize=1e6,
no.histogram=FALSE, max.segments=200, check.unit=FALSE, 
category.col2=NULL, plot.range=NULL, orig.range=NULL, chroms=NULL) {
    if(is.null(datatab))
        return(NULL)
    if(nrow(datatab)==0L)
        return(NULL)
        
    # Filter chromosomes
    if(!is.null(chroms)) {
        if(any(!chroms %in% datatab[,bycol]))
            warning("Chromosomes not in the dataset: ", paste(setdiff(chroms, 
                    datatab[,bycol]), collapse=", "))
        datatab <- datatab[datatab[,bycol] %in% chroms,]
    }
        
    groupby <- NULL
    Group  <- match(category.col, colnames(datatab))
    Group  <- Group[!is.na(Group)][1]#

    # Find position column
    poscoli <- match(poscol, colnames(datatab))
    if(sum(is.na(poscoli))==length(poscol)) {
        stop("'", paste(poscol[is.na(poscoli)], collapse=", "),
        "' is not a column of datatab.")
    } else {
       poscol  <- poscol[!is.na(poscoli)][1L]
    }
    
    # Remove rows with missing values for postions columns
    datatab <- datatab[!is.na(datatab[,bycol]) & !is.na(datatab[,poscol]),]
    if(check.unit)
        datatab   <- checkunit(datatab)
    
    maxseglen <- max(datatab[,end.col] - datatab[,start.col] - 1)

    
    if(!is.null(category.col)) {
        if(category.col %in% colnames(datatab)) {
            col.idx <- match(category.col, colnames(datatab))
            colnames(datatab)[col.idx] <- "Category"
        }
    }

    if(!is.null(category.col2)) {    
        if(category.col2 %in% colnames(datatab)) {    
            col.idx2 <- match(category.col2, colnames(datatab))
            colnames(datatab)[col.idx2] <- "Category2"
        }
    }

    if("gieStain" %in% names(datatab)) {
        datatab     <- cytoBands(datatab)
    } 
    
    if((maxseglen < binsize | nrow(datatab) >= max.segments) & !no.histogram){
    
        # build track histogram
        if("Category" %in% colnames(datatab)) {
            groupby=split(datatab$Category, datatab[,bycol])
        }
        track_ls <- list()
        for(chr in sort(unique(datatab[,bycol]))) {
            rows <- datatab[,bycol] %in% chr
            track_ls[[chr]] <- split.bins(datatab[rows,poscol], step=binsize, 
                                            fun=length, group.by=groupby[[chr]])     
        }
        track_ls <- track_ls[!sapply(track_ls, is.null)]
        track_ls <- sort.chrom(track_ls)
        # Adjust values to preset scale if provided
        if(!is.null(plot.range)){
            if(length(plot.range)!=2L) {
                stop("plot.range must ve a vector with 2 values, min and max.")
            }
            if(is.null(orig.range)) {
                if(is.null(dim(track_ls[[1]]$counts))){
                    cmin <- min(sapply(track_ls, function(x) min(x$counts)))
                    cmax <- max(sapply(track_ls, function(x) max(x$counts)))
                } else {
                    cmin <- min(sapply(track_ls,
                                function(x) min(rowSums(x$counts))))
                    cmax <- max(sapply(track_ls, 
                                function(x) max(rowSums(x$counts)))) 
                }
            } else {
                cmin <- orig.range[1]
                cmax <- orig.range[2]
            }
            # forces min count to 0 for histograms
            track_ls <- lapply(track_ls, 
                            function(x){x$counts=scale.vec(x$counts, 0, cmax,
                            plot.range[1], plot.range[2]);x}) 
            attr(track_ls, "count_range") <- c(cmin, cmax)
        }
        attr(track_ls, "is_int") <- FALSE
    }else{
        track_ls <- split(datatab, datatab[,bycol])
        track_ls <- sort.chrom(track_ls)
        attr(track_ls, "is_int") <- TRUE        
    }
         
    return(track_ls)
} 

build.stat <- function(stat=NULL, statcol=NULL, maxgenecount=5, stat.typ="p", 
scex=1, spty=20, binsize=1e6, margin=margin, statthreshold=NULL,
col.stat="blue", col.stat2="gray57",  statSumm="none", bycol="Chrom",
poscol=c("Mid", "Start"), start.col="Start", end.col="End", id.col="ID",
colors.col="Colors", chroms=NULL) {
    # Sanity checks
    if(!statcol %in% colnames(stat) & !id.col %in% colnames(stat)) {
        stop(paste(statcol, " not a column of stat object."))
    } 
    if(id.col %in% colnames(stat) & statSumm != "none") {
        warning("the statSumm argument should be 'none' if an 'ID'",
        " columns is present. statSum will be ignored for this track.")  
    }
    # Filter chromosomes
    if(!is.null(chroms)) {
        if(any(!chroms %in% stat[,bycol]))
            warning("Chromosomes not in the dataset: ", paste(setdiff(chroms, 
                    stat[,bycol]), collapse=", "))
        stat <- stat[stat[,bycol] %in% chroms,]
    }
    # Find position column
    poscoli <- match(poscol, colnames(stat))
    if(sum(is.na(poscoli))==length(poscol)) {
        stop("'", paste(poscol[is.na(poscoli)], collapse=", "),
             "' is not a column of stat.")
    } else {
        poscol  <- poscol[!is.na(poscoli)][1L]
    }
        
    # Declare variables
    maxstat     <- NULL
    minstat     <- NULL
    Realmaxstat <-NULL
    # Always order stat by position
    stat      <- stat[order(stat[,bycol], stat[,poscol], decreasing=FALSE ) , ] 
        
    if(id.col %in% colnames(stat)) { 
        stat <- stat[!is.na(stat[,poscol]), ]
        stat <- stat[!is.na(stat[,bycol]), ]
        stat[, id.col]<-as.character(stat[,id.col])
        NaID   <- is.na(stat[,id.col])
        NULLID <- is.null(stat[,id.col]) 
        stat[NaID, id.col] <- "" 
        
        if(!statcol %in% colnames(stat)) {
            stat <- cbind(stat, margin)
            colnames(stat)[ncol(stat)] <- statcol
        }
        
        if( colors.col %in% colnames(stat)) {
            stat <-stat[,c(bycol, poscol, end.col, id.col, statcol, colors.col)]
        } else {
            stat[,colors.col] <-"black"
            stat <- stat[,c(bycol, poscol,end.col, id.col,statcol, colors.col)]
        }   
    } 

    if(statSumm== "none" | (id.col %in% colnames(stat))){
        stat_vals<-stat
        colnames(stat_vals)[match(c(poscol,statcol),
                                  colnames(stat_vals))]<-c("mids", "statSumm")
        trackStat     <- split (stat_vals, as.factor(stat_vals[,bycol]))
        maxstat        <- max(stat[,statcol])
        minstat        <- min(stat[,statcol])
        Realminstat    <- minstat           
        Realmaxstat    <- maxstat           
    } else{
        stat_vals  <- data.frame(as.numeric(stat[,poscol]),
                                as.numeric(stat[,statcol]))
        chromstat  <- split(stat_vals, as.factor(stat[,bycol]))   
        trackStat  <- lapply(chromstat, split.bins, step=binsize, fun=statSumm)
        trackStat  <- trackStat[!sapply(trackStat, is.null)]
        maxstat   <- max(unlist(sapply(trackStat, function(x) max(x$statSumm))))
        minstat   <- min(unlist(sapply(trackStat, function(x) min(x$statSumm))))
        Realminstat <- minstat         
        Realmaxstat <- maxstat
    } #else

    if(!is.null(statthreshold)){
        trackStat    <- lapply(trackStat, function(x){
                        x$Colors=ifelse(x$statSumm >= statthreshold,
                                        col.stat, col.stat2);x} )
    } else {
        trackStat    <- lapply(trackStat, function(x){
                        x$Colors=col.stat; x} )
    }    
    sdstat         <- sd(unlist(sapply(trackStat, function(x) x$statSumm)))

    if(sdstat ==0) { # Value was a constant
        trackStat <- lapply(trackStat, function(x){x$statSumm=margin;x})
    } else { # Value varied, must scale
#Use this code to increase observable differences by removing the minimum value.  
         trackStat  <- lapply(trackStat, 
                          function(x){x$statSumm=(x$statSumm-minstat)/
                                      (maxstat-minstat)*
                                      (maxgenecount-margin)+margin;x})
    }
                                
    maxstat  <- max(unlist(sapply(trackStat, function(x) max(x$statSumm))))
    minstat  <- min(unlist(sapply(trackStat, function(x) min(x$statSumm))))

    attr(trackStat, "colScaleStat")  <- col.stat
    attr(trackStat, "is_int" )       <- FALSE
    attr(trackStat, "stats.typ" )    <- stat.typ
    attr(trackStat, "scex" )         <- scex #
    attr(trackStat, "spty" )         <- spty #   
    attr(trackStat, "maxstat" )      <- maxstat #
    attr(trackStat, "Realminstat" )  <- Realminstat #
    attr(trackStat, "Realmaxstat" )  <- Realmaxstat #
    attr(trackStat, "minstat" )      <- minstat
    attr(trackStat, "sdstat")        <- sdstat
     
    return(trackStat)
}



## Plotting functions
## -------- ---------

plot.track <- function(track_ls, which.track, margin, side=1, lwd=3,
     sort_segs=TRUE, col_segs=NULL, ylims=NULL, stack=FALSE, form=NULL, bin=1e6,
     cex=NULL, trckmar=0, legchr=NULL, title=NULL, overlap.ok=FALSE) {
    
    if(is.null(track_ls))
        return(NULL)
    
    track   <- track_ls[[which.track]]
    
    if("ID" %in% colnames(track_ls[[which.track]]) ){ 
        if(any(track$ID != "")) {
            linelength <- strwidth("xx", units="user", cex=cex)
            labtrack <- track[track$ID != "",]
            unstacked <- unstacked(x= 1, y=labtrack$mids, trackId= labtrack, 
            cex=cex)
            if("Colors" %in% colnames(track)){
                 lab_colors <- as.character(track$Colors)
            } else {
                lab_colors <- col_segs
            }
            segments (side*labtrack$statSumm, labtrack$mids,
                        side*(labtrack$statSumm+linelength),
                        unstacked$y, col="black", lwd=0.7, lend=1)
            text(side*(labtrack$statSumm+linelength), unstacked$y, 
                unstacked$ID, col=unstacked$Colors, adj=0.5-(side*0.5),
                lwd=lwd, cex=cex)
        }
    }

    typ     <- attr(track_ls, "stats.typ")
    scex    <- attr(track_ls, "scex")
    spty    <- attr(track_ls, "spty")
    flag    <- 0 

    if(attr(track_ls, "is_int")) {
        if("Colors" %in% colnames(track_ls[[1]])){
             col_segs   <- sort(unique(unlist(sapply(track_ls, "[", "Colors"))))
        }
        if("Category" %in% colnames(track_ls[[1]])) {
            all_categories1 <- unlist(sapply(track_ls, "[", "Category"))
            slevs1 <- sort(unique(all_categories1))
            nlevs1 <- length(slevs1)
            if(length(col_segs)<nlevs1) {
                 stop("Not enough colors provided. There are ", nlevs1, 
                      " categories of elements and ", length(col_segs),
                      " colors.")
             }
        }
        if("Category2" %in% colnames(track_ls[[1]])){
            seg_form        <- c(15, 16, 17, 18, 25, 20:30)
            all_categories2 <- unlist(sapply(track_ls, "[", "Category2"))
            slevs2 <- sort(unique(all_categories2))
            nlevs2 <- length(slevs2)
            if(length(seg_form)<nlevs2) {
                stop("Not enough forms provided. There are ", nlevs2,
                    " categories of elements and ", length(seg_form),
                    " forms.")
            }
        }
    }
    
    if(!is.null(track) ) {
        if(attr(track_ls, "is_int")) {
            if("Colors" %in% colnames(track)){
                seg_colors <- as.character(track$Colors)            }
            if("Category" %in% colnames(track)) {
                seg_colors <- col_segs[match(track$Category, slevs1)]
            } else if(!"Colors" %in% colnames(track)) {
                seg_colors <- rep(col_segs[1L], nrow(track))     
            }
            if("Category2" %in% colnames(track)){
                form <- seg_form[match(track$Category2, slevs2)]
            }
        }
        if(!is.null(typ)) {   
            if("Colors" %in% names(track)) {
                seg_colors <- track$Colors
            } else {
                seg_colors <- rep(col_segs[1L], length(track$mids))     
            }
            flag=1
            switch(typ,
            l=lines(side * (track$statSumm), track$mids,
                    col=seg_colors, type="l"),
            p=lines(side * (track$statSumm), track$mids, 
                    col=seg_colors, type="p", n=NULL, pch=spty, cex=scex),
            hst=attr(track_ls, "is_int")<-FALSE,
            seg=attr(track_ls, "is_int")<-FALSE)     
        }

        if(attr(track_ls, "is_int")) {
            if(overlap.ok) {
                segments(margin*side, track[,"Start"],margin*side,track[,"End"],
                         col=seg_colors, lwd=lwd, lend=1)
            } else {
                if(!stack){
                    plot.lodclust(as.matrix(track[,c("Start", "End")]),
                                  col=seg_colors, lwd=lwd, ymax=margin*side,
                                  stack_ints=FALSE, form=form)                  
                } else {
                    plot.segments(track, seg_colors, startfrom=side*margin, 
                                  lwd=lwd, sort.segs=sort_segs, form=form)
                }
            }
        } else{
            if(flag==0){
            cat("Small segments, I will plot histograms for this track:\n")
                if(trckmar>margin) {
                    track$counts <- scale.vec(track$counts, margin, 1,trckmar,1)
                    margin <- trckmar
                }
                if(!is.null(dim(track$counts))) {
                     plot.stackedHist(track, margin, ylims=NULL, col=col_segs[1], 
                                      stepsize=bin, side=side)
                } else {
                     chromhist(side*margin, side*(track$counts),
                     track$mids, ylims=NULL, col=col_segs[1], stepsize=bin)#bin
                }
            }   
        }
    }    

    # Plot a legend for is_int with Category
    if(!is.null(legchr)) if(!is.na(legchr)) {
        if("Category"%in%colnames(track_ls[[1]])&attr(track_ls, "is_int") & 
            !is.null(legchr)){
            if(which.track == legchr) {
                lcols  <- col_segs[1:nlevs1]
                legend("bottom",legend=slevs1,ncol=max(1,floor(sqrt(nlevs1))-1), 
                       lwd=3, col=lcols, lty=1, xpd=NA,
                cex=cex, bg="white", title=title)
            }
        }

        if("Category2"%in%colnames(track_ls[[1]])&attr(track_ls, "is_int") &
            !is.null(legchr)){
            if(which.track == legchr) {
                lform  <- seg_form[1:nlevs2]
                legend("right",
                legend=slevs2, ncol=max(1,floor(sqrt(nlevs2))-1),
                lwd=3, col="black", lty=1, xpd=NA, pch=lform,
                cex=cex, bg="white", title=paste(title, "2nd Category"))
            }
        }
    }
} # plot.track


plot.stackedHist <- function(track=NULL, margin=NULL, ylims=NULL, col=NULL, 
stepsize=1e6, side=1) {
    data <- t(track[[2]]) # assuming that first elements is always mids

    for(i in 1:nrow(data)){
        if(i==1){
            x1 <- margin
            x2 <- as.numeric(data[i,])    
        } else{
            x1 <- x2
            x2 <- x1 + as.numeric(data[i,]) - margin
        }
        chromhist(side*x1, side*x2, track$mids, col=col[i], stepsize=stepsize,
                ylims=ylims)
    }
} # plot.stackedHist

plot.int <- function(intervals, intclusters, col, lwd=3, ymax, stack_ints=TRUE, 
form) {
    intervals <- as.matrix(intervals)
    if(ncol(intervals)==1)
        intervals <- t(intervals)
    for(i in 1:length(intclusters)){
        incluster=intclusters[[i]]
        plot.lodclust(intervals[incluster,], col=col[incluster], lwd, ymax,
        stack_ints=stack_ints, form=form[incluster])
    }
} 
chromhist <- function(xleft, xright, y, ylims=NULL, col="red", stepsize=1) {
    ytop <- y - stepsize/2
    ybot <- y + stepsize/2 
    if(!is.null(ylims)) {
        ytop <- unlist(sapply(ytop, max, ylims[1]))
        ybot <- unlist(sapply(ybot, min, ylims[2]))
    }
    rect(xleft, ybot, xright, ytop, col=col, border=NA)
}


draw.scale <-function(y, minval, maxval, lwd, col, minlab=minlab, maxlab=maxval,
title="Counts", cex, ...) { 
    lht      <- sign(diff(par("usr")[3:4])) * strheight(title)

    tickpos  <- pretty_ticks(minval, maxval, 2, ...)
    ticklab  <- pretty_ticks(minlab, maxlab, 2, ...)
    
    ticklab <- format.bignum(ticklab, unit="", sep="")
    tickpos <- tickpos - tickpos[length(tickpos)] *0.5 # move it 50% to the left
    
    lines(c(tickpos[1], tickpos[length(tickpos)]), rep(y+lht*2,2), col=col,
    lwd=lwd, lend=1)
    
    segments(tickpos, y+lht*1, tickpos, y+lht*3, col="black", lwd=1)
    text(tickpos[1] + (tickpos[length(tickpos)]-tickpos[1])*.7, y+lht*0,
    title, adj=.5, font=1, cex=cex) 
    text(tickpos, y+lht*4, ticklab, cex=cex) # Bottom side of scale
} 

plot.segments <- function(segobj, colvec, startfrom, lwd=3, sort.segs=TRUE,
form=NULL) {
    intsmat <- segobj[, c("Start", "End")]
    intslen <- intsmat[,2]-intsmat[,1]
    clusts  <- find.intclusters(intsmat, sort.by.size=sort.segs)
    plot.int(intsmat, clusts, colvec, lwd=lwd, startfrom, stack_ints=TRUE,
            form=form)
}

plot.lodclust<-function(intmat, col, lwd=3, ymax, stack_ints=TRUE, form=NULL){
    intmat  <- matrix(intmat, ncol=2)
    lwdusr  <- pt2usr(lwd,  axisunit="x", lpi=96)
    onestep <- lwdusr * 1.2 * sign(ymax)
    y       <- ymax
    idx     <- 1L:nrow(intmat)
    assigned_levs <- idx
    for(i in 1L:nrow(intmat)) {
        if(i > 1L) {
            if(stack_ints) {
                used_lev<-assigned_levs[sapply(idx,function(x) overlap(intmat[i,],
                intmat[x,]) & x!=i)]
                clear <- min(setdiff(idx, used_lev))
                assigned_levs[i] <- clear
            } else {
                clear <- i
            }   
            y <- ymax + onestep * (clear-1L)
        }
    
        plot.lodint(intmat[i,], col=col[i], lwd=lwd, y=y,form=form[i])
    }
}

plot.lodint<-function(lodint, lty=2, col=1, lwd=3, y=NULL, form=1) {
    lowint<-numeric()
    highlim<-numeric()
    if(length(dim(lodint))>1) {
        lowlim<-lodint[1,2]
        highlim<-lodint[nrow(lodint),2]
    }else{
        lowlim<-lodint[1]
        highlim<-lodint[2]
    }
    if(!is.null(form)) {
       line_col <- "darkgray"
    } else {
       line_col <- col
    }
     
    lines(rep(y, 2), c(lowlim, highlim), col=line_col, lwd=lwd, lend=1)
    if(!is.null(form)){
        points(y, mean(c(lowlim, highlim)), col=col, lwd=lwd, lend=1,
        pch=form, cex=1)
    }
     
}

## Auxiliary functions
## --------- ---------

split.bins <-function(dataframe, column=NULL, step, fun=length, group.by=NULL){
## splits data by bins of 'step' size and apply 'fun' to each cell resulting 
## from the combination of levels bin and the group.by column in dataframe
## (if provided)
##
## dataframe    object to split n bins
## column       column to use as variable to create bins. Can 
##                be a integer or a string with colname
## step         bin size. values > 1 are used as absolute
##                size or a fraction of the range in column

    if(is.null(fun))
        return(dataframe)
        
#     if(is.character(fun))
#         fun <- get(fun)
        
    if(is.null(column))
        x<-dataframe
    else
        x<-dataframe[,column] 
        
    if(is.data.frame(x) | is.matrix(x)) {
        if(ncol(x)>1) {
            y=x[,2]
            x=x[,1]
        } else {
            y=x
        }
    }else{
        y=x
    } 
    y<-y[!is.na(y)]
    x<-x[!is.na(x)] 
    if(identical(length(x), 0L))
        return(NULL) 
    span<-diff(range(x, na.rm=TRUE))   
    if(step<=1)
        step<-span*step  
    bin<-ceiling(x/step)  
    mids<-round((bin*step+(bin-1)*step)/2)

    if(is.null(group.by)) {
        byfactor <- list(mids)
    } else {
        byfactor <- list(mids, group.by)
    }

    out<-tapply(y, byfactor, fun)
    if(!is.null(group.by) & identical(fun, length)) {
        out[is.na(out)] <- 0
    }


    outname<-deparse(substitute(fun))
    if(grepl('length', outname)) outname<-'counts'
    if(is.null(group.by))
        out <- as.numeric(out)
        
    output<-list(mids=sort(unique(mids)), out)
    names(output)[length(output)]=outname
    
    return(output)
} # split.bins

chrplotMar <- function(mar, sides, tracks, lwd, lpi=96, stacked=FALSE) {
## Calculate appropriate margin for each track according to what other tracks
## are plotted on the same chromosomal side. If a histogram and segments are
## plotted on the same side, the histogram start where segments end.
    if(length(sides)!=length(tracks))
        stop("Arguments sides and tracks must be of same length.")
    if(any(!sides %in% c(-1,1)))
        stop("values of sides must be in c(-1, 1)")
        
    out      <- rep(mar, length(sides))
    os       <- pt2usr(lwd,  axisunit="x", lpi=lpi)
    are_drawn<- !sapply(tracks, is.null)
    are_ints <- sapply(tracks, attr, "is_int")
    are_ints <- unlist(sapply(are_ints, function(x) if(is.null(x)) FALSE else x))
    off_set  <- which(are_drawn & !are_ints)
    side_height <- list(`-1`=NULL, `1`=NULL)
    

    ## If plotting int segments and stacking, calculate max cluster size to
    ## adjust margin
    track_height <- rep(1, length(tracks))
    if(sum(are_ints)>0 & stacked) {
      track_height[are_ints] <- sapply(tracks[are_ints], max.cluster.height)
    }
    # Max height per chrom side
    side_height_tmp <- tapply(track_height, sides, max) 
    side_height[names(side_height_tmp)] <- side_height_tmp
    
    # Add the hight of one line each segment in the same side as a hist
    if(length(off_set)>0L) {
        for (i in 1:length(off_set)) {
            if(any(sides[-i] %in% sides[i] & are_ints[-i]))
                out[i] <- mar + os * side_height[[sides[i]*.5+1.5]]
        }
    }
    return(out)
}

format.bignum <- function(number, unit="b", sep=" ") {
    if(length(number)>1)
        return(unlist(sapply(number, format.bignum, unit, sep)))
    
    gigas <- number/1e9
    megas <- number/1e6
    kilos <- number/1e3
    gig_u <- paste(sep="", "G", unit)  
    meg_u <- paste(sep="", "M", unit)
    kil_u <- paste(sep="", "K", unit) 
    if(gigas[length(gigas)] == 0)
        return(0)
    if(gigas[length(gigas)] == as.integer(gigas[length(gigas)]))
        return(paste(gigas, gig_u, sep=sep))
    else if(megas[length(megas)] == as.integer(megas[length(megas)]))
        return(paste(megas, meg_u, sep=sep))
    else if(kilos[length(kilos)] == as.integer(kilos[length(kilos)]))
        return(paste(kilos, kil_u, sep=sep))
    else
        return(paste(format(number, big.mark=",", trim=TRUE), unit, sep=sep))
}

max.cluster.size <- function(segslist) {
    ints_by_chrom   <- lapply(segslist, function(x) x[,c("Start", "End")])
    clusts_by_chrom <- lapply(ints_by_chrom, find.intclusters,
                                sort.by.size=FALSE)
    output <- max(sapply(clusts_by_chrom[[1]], length))
    return(output)
}

max.cluster.height <- function(segslist) {
    ints_by_chrom   <- lapply(segslist, function(x) x[,c("Start", "End")])
    clusts_by_chrom <- lapply(ints_by_chrom, find.intclusters, 
                                sort.by.size=FALSE)

    ints_by_clust   <- list()
    for (chr in 1:length(ints_by_chrom)) {
        ints_by_clust[[chr]] <- lapply(clusts_by_chrom[[chr]], 
                                        function(x) ints_by_chrom[[chr]][x,])
    }

    output <- max(unlist(sapply(ints_by_clust,
                                function(x) sapply(x, maxclustsize))))
    return(output)
}
## this is just an approximation to the number of levels needed for plotting
## a clusters
maxclustsize <- function(ints) {
   tot_size <- sum(apply(ints, 1, diff))
   
    min_possible_size <- max(ints[,2])-min(ints[,1])
    max_possible_size <- mean(ints[,2])-mean(ints[,1])
    estimate1 <- tot_size/min_possible_size
    estimate2 <- tot_size/max_possible_size
    estimated_size <- ceiling((estimate1+estimate2)/2)
    return(estimated_size)
}
scale.track <- function(track, stat.col="statSumm", mintarget, maxtarget) {
## Converts scale of values in a track to a new range of values
## (mintarget,maxtarget)
    maxstat     <- max(unlist(sapply(track, function(x) max(x$statSumm))))
    minstat     <- min(unlist(sapply(track, function(x) min(x$statSumm))))
  
    for(i in length(track)) {
        track[[i]] <- scale.vec(track[[i]][,stat.col], minstat, maxstat,
                                mintarget, maxtarget)
    }

    return(track)
}

scale.vec <- function(vec, minstat, maxstat, mintarget, maxtarget) {
    (vec-minstat)/(maxstat-minstat)*(maxtarget-mintarget)+mintarget
}

checkunit <- function(featable, column=c("Start", "End"), minlen=10e6) {
    if(max(featable[,column], na.rm=TRUE)<minlen) { # it is not Mb
            featable[,column] <- featable[,column] * 1e6
    }
return(featable)
}

sort.chrom <- function(chroms) {
    if(is.list(chroms))
        return(sort.chrom.list(chroms))
    chroms <- as.character(chroms)
    dchroms <-sub("chr", "", chroms)  
    arenum <- grep("^[[:blank:]]*[[:digit:]]+[[:blank:]]*$", dchroms)
    aresex <- setdiff((1:length(chroms)), arenum)
    outchroms <- chroms[arenum][order(as.numeric(dchroms[arenum]))]
    outchroms <- c(outchroms, sort(chroms[aresex]))
    return(outchroms)
}

sort.chrom.list <- function(chromList) {
    chroms <- names(chromList)
    return(chromList[sort.chrom(chroms)])
}

pretty.choose <- function(number, num.options) {
    if(missing(num.options))
        num.options <- c(0, 5, 10, 20, 50, seq(100, 1000, by=100),
        seq(2e3, 1e4, by=1e3), seq(2e4, 1e5, by=1e4), seq(2e5, 1e6, by=1e5),
        seq(1e-01, 1, by=1e-01))
    if(length(number)>1)
        return(unique(unlist(sapply(number, pretty.choose, num.options))))
    output <- number
    num.options <- sort(num.options, decreasing=TRUE)
    for(i in 1:length(num.options)) {
        if(number > num.options[i] | i==length(num.options)) {
            output <- num.options[i]
        break
        } else if(number > num.options[i+1]) {
            if((num.options[i] - number) <= (number - num.options[i+1]))
                output <- num.options[i]
            else
                output <- num.options[i+1]
            break
        } else 
    next
    }
    return(output)
}

pretty_ticks <- function(minval, maxval, howmany=3, as.is=FALSE, digits=2, ...){
    
    outticks <- signif(seq(minval, maxval, length=howmany), digits=digits)
    
    if(!as.is) {
        outticks <- pretty.choose(outticks)
    }

    return(outticks)
}

pt2usr <- function(lwd,  axisunit="x", lpi=96) {
    if(!axisunit %in% c("x", "y"))
        stop("Argument 'axisunit' must be in c('x', 'y').")
    pusr   <- par("usr")
    pinch  <- par("pin")
    lwdin  <- lwd * 1/lpi
    if(axisunit=="x")
        usr.in <- (pusr[2]-pusr[1])/pinch[1]
    else
        usr.in <- (pusr[4]-pusr[3])/pinch[1]
        lwdusr <- lwdin * usr.in
    return(lwdusr)
}

usr2pt <- function(lwdusr, axisunit="x", lpi=96) {
# lwdusr  line width in user coordinates units
# lpi     lines per inch (device dependent)
    if(!axisunit %in% c("x", "y"))
        stop("Argument 'axisunit' must be in c('x', 'y').") 
    pusr   <- par("usr")
    pinch  <- par("pin")  
    if(axisunit == "x")
        usr.in <- (pusr[2]-pusr[1])/pinch[1]
    else
        usr.in <- (pusr[4]-pusr[3])/pinch[1]
    lwdin <- lwdusr / usr.in
    lwdpt <- lwdin * lpi
    return(lwdpt)
}
find.intclusters <- function(intmat, sort.by.size=TRUE){
    indices  <- 1:nrow(intmat)
    clusters <- list()
    i=0
    while(length(indices)>0L) {
        i=i+1
    
        lrgst    <- which.max(intmat[indices,2]-intmat[indices,1])
        members  <- indices[lrgst]
        secindx <- setdiff(indices, members)
        if(length(secindx)>0L) {
            member_mem <- 0
            while(!all(members %in% member_mem)) {
                member_mem <- members
                for(j in secindx) {
                    if(overlap(intmat[members,], intmat[j,])) {
                        members <- unique(append(members, j))
                        secindx <- setdiff(indices, members)            
                    }
                }
            }
        }
        if(sort.by.size)
            members <- members[order(intmat[members,2]-intmat[members,1],
            decreasing=TRUE)]
            clusters[[i]] <- members
            indices=setdiff(indices, members)
    }
    cat("  ", length(clusters), "segment clusters.\n")
    return(clusters)
}
overlap<-function(int1, int2) {
    if(length(int1)==0L | length(int2)==0L)
        return(FALSE)
    if(!is.null(dim(int1)))
        int1 <- c(min(int1[,1]), max(int1[,2]))
    if(!is.null(dim(int2)))
        int2 <- c(min(int2[,1]), max(int2[,2]))    
    int1 <- as.numeric(sort(int1))
    int2 <- as.numeric(sort(int2))
    len1 <- diff(int1) 
    len2 <- diff(int2)
    total_len <- len1 + len2  
    span <- max(int1, int2) - min(int1, int2)
    if(length(total_len)!=0L & length(span)!=0L)
        if(total_len >= span) {
            return(TRUE)
        }
    return(FALSE)
}

checkStructure <- function(data = NULL, rmrnd=FALSE, 
                           columnNames=c("Chrom", "Start", "End"),
                           optional.columns=NULL) {

    if(is.null(data)) return(NULL)

    # GenomicRanges
    if("GRanges" %in% class(data)) {
        data <- data.frame(seqnames(data), start(data), end(data), mcols(data))
        colnames(data)[1:3] <- columnNames
    }

    if(is.character(data)) {
        # I will assume that this is a file name
        if(!file.exists(data) &  !grepl("http", data))
            stop("File ", data, "does not exist.")
            
        Lines<-read.delim(data, sep="\n", header=FALSE, fill=FALSE, 
                            comment.char="#")
        if(grepl("[A|C|G|T]{1,}", Lines[2,]) ){
            cat ("Reading AXT file", "\n")
            data <- transformAlign(Lines)    
        } else {
            cat("Reading file ", data, " assuming bed format...", "\n")
            data <- read.delim(data, sep="\t", header=TRUE, fill=FALSE,
                                comment.char="#")
        }
    } 
        
    name <- deparse(substitute(data))

    if(!is.data.frame(data)){
        # Everything failed....
        stop(name, " could not be read.")
    }

    datacols <- colnames(data)
    
    if(!all(columnNames %in% datacols)){ 
        stop(name, " is missing some column(s): ", 
             paste(setdiff(columnNames, datacols), collapse=", "))
    }
    areNA <- apply(is.na(data[,columnNames]), 1, any)
    if(sum(areNA)>0L) {
        warning(sum(areNA), " rows have missing values for columns ",
                paste(columnNames, collapse=", "), 
                "and will be removed.")
        data <- data[!areNA,]
    }
    
    # Assuming first column name is Chromosome
    data[,columnNames[1]] <- as.character(sub("^chr", "",
                                          as.character(data[,columnNames[1]])))
   
    if(rmrnd) {
        #Remove bad chrom
        bad_chroms <-grepl("_",  data[,columnNames[1]])
        if( length ( bad_chroms )>0L )
            data <- data[!bad_chroms,]
    }
    
    if(!is.null(optional.columns)) {
        for(i in 1: length(optional.columns)) {
            if(optional.columns[i] %in% datacols)
                data[,optional.columns[i]] <- as.character(data[,optional.columns[i]])            
        }
    }

    return(data)
}


selectchr <- function(chrom, data) {
    if( any (!chrom %in% names(data))) 
        stop("Some chromsomes in 'chr' cannot be drwan.") 
        data <- data[names (data) %in% chrom]
    return (data)
}
transformAlign <- function(filename = NULL){
    y      <- nrow(filename) # size file
    row    <- seq(1, y, by=3) # vector with the number of rows that
    ## contains chromosomal position
    data <- filename[row, ]
    data<-as.character(data)
    data <- strsplit(data, " ", fixed=TRUE)
    matrix <- t(sapply(data, unlist)) #list to matrix
    aux    <-matrix[, c(2:5)]
    dataBed <- as.data.frame(aux) 
    colnames(dataBed) <- c("Chrom", "Start", "End", "Group")
    dataBed$Start     <-as.numeric(as.character(dataBed$Start))
    dataBed$End       <-as.numeric(as.character(dataBed$End))
    dataBed$Group     <-as.character(dataBed$Group) 
    dataBed$Chrom     <-sub("^chr([^:digit:]|[M|X|Y|].+)", "\\1", dataBed$Chrom)
       
    return(dataBed)
}
getAnnotBiomaRt <-function(org=NULL, annotBiomaRt=NULL ) {  
    if(!is.null(org) ){ 
        if(annotBiomaRt){

            dataset <-paste(org, "_gene_ensembl", sep="")
            ensembl <- useMart("ENSEMBL_MART_ENSEMBL",dataset=dataset,
                        host="www.ensembl.org")
            annot<-getBM(attributes=c("chromosome_name","start_position", 
            "end_position","ensembl_gene_id", "strand"),
            filters="chromosome_name",
            values=as.list(as.character(c(1:40, paste(1:40, "L", sep=""),
                           paste(1:40, "R", sep=""), "X", "Y")))
            , mart = ensembl)
            
            colnames(annot)[1]<-"Chrom"
            colnames(annot)[2]<-"Start"
            colnames(annot)[3]<-"End"
            colnames(annot)[4]<-"Name"
            colnames(annot)[5]<-"Strand"    
            return(annot)
        } 
    }
}


#For plotting ID from tables
unstacked<-function(x=NULL, y=NULL, trackId=NULL, cex=NULL){

    rest      <- NULL
    trackId$x <- x
    trackId$y <- y
 
    for(i in 1:nrow(trackId)){
        rest <- 0  
        if(i<=nrow(trackId)-1){
            rest<-trackId[i+1,"mids"] - trackId[i,"mids"]
            offset <- trackId[i,"mids"] + abs(strheight(trackId[i,"ID"], cex=cex)*1.05)- trackId[i+1,"mids"]
            if(offset>0){
              trackId[i+1, "y"]<- trackId[i+1, "y"]+offset    
            } 
        }
     
    }
    return(trackId)
}
cytoBands<- function(data=NULL){
    if(!"Colors" %in% names(data)) {
        if("gieStain" %in% names(data)) {
            data        <- data[!is.na(data[ ,"gieStain"]), ]
            type.bands <- match(data$gieStain, c("acen", "gneg", "gpos",
            "gvar", "stalk", "gpos25", "gpos50","gpos66", "gpos75", "gpos100"))
            bandcol<- c("#CCCCCC", "grey96", "#000000", "grey96",  "#CCCCCC", 
                    "grey90", "grey70","#666666", "grey40", "grey0")[type.bands]
            data$Colors<-bandcol
        } else{
                stop("The parameter bands does not contain a color column or",
                " Giemsa stain results (UCSC format).")
          }
    } else{
        if("Group" %in% names(data)) {
            data<-NULL
        }
    }
    return(data)
}

makeCache <- function() {
    cache <- new.env(parent=emptyenv())
    list(get = function(key) cache[[key]],
         set = function(key, value) cache[[key]] <- value,
         load = function(rdaFile) load(rdaFile, cache),
         ls = function() ls(cache))
}

# Save an object in a chache
saveObj <- function(object, oname, cname="chromPlot-cache") {
      chrpltchache <- new.env(parent=emptyenv())
      assign(oname, object, envir=chrpltchache)
      attach(chrpltchache, name=cname)
} # saveObj


# Get object from cache. If not there, returns NULL silently
getObj <- function(oname, cname="chromPlot-cache") {
    if(exists(oname)) {
        out <- get(oname)
        return(out)
    }
    return(NULL)
} # getObj