R/lollipop.plot.R

In jianhong/trackViewer: A R/Bioconductor package with web interface for drawing elegant interactive tracks or lollipop plot to facilitate integrated analysis of multi-omics data

#' Lolliplots
#' @description Plot variants and somatic mutations
#' @param SNP.gr A object of \link[GenomicRanges:GRanges-class]{GRanges}, 
#' \link[GenomicRanges:GRangesList-class]{GRangesList}
#' or a list of \link[GenomicRanges:GRanges-class]{GRanges}.
#' All the width of GRanges must be 1.
#' @param features A object of \link[GenomicRanges:GRanges-class]{GRanges}, 
#' \link[GenomicRanges:GRangesList-class]{GRangesList}
#' or a list of \link[GenomicRanges:GRanges-class]{GRanges}. 
#' The metadata 'featureLayerID' are used for drawing features in different layers.
#'  See details in vignette.
#' @param ranges A object of \link[GenomicRanges:GRanges-class]{GRanges} or 
#' \link[GenomicRanges:GRangesList-class]{GRangesList}.
#' @param type character. Could be circle, pie, pin, pie.stack or flag.
#' @param newpage Plot in the new page or not.
#' @param ylab Plot ylab or not. If it is a character vector, 
#' the vector will be used as ylab.
#' @param yaxis Plot yaxis or not.
#' @param xaxis Plot xaxis or not. If it is a numeric vector with length greater than 1, 
#' the vector will be used as the points at which tick-marks are to be drawn. 
#' And the names of the vector will be used to as labels to be placed at the tick 
#' points if it has names. 
#' @param ylab.gp,xaxis.gp,yaxis.gp An object of class gpar for ylab, xaxis or yaxis.
#' @param legend If it is a list with named color vectors, a legend will be added.
#' @param legendPosition The position of legend. Possible positions are 'top',
#' 'right', and 'left'.
#' @param cex cex will control the size of circle.
#' @param dashline.col color for the dashed line.
#' @param jitter jitter the position of nodes or labels.
#' @param rescale logical(1), character(1), numeric vector,
#'   or a dataframe with rescale from and to. 
#' Rescalse the x-axis or not.
#' if dataframe is used, colnames must be from.start, from.end,
#'  to.start, to.end. And the from scale must cover the whole plot region.
#' The rescale parameter can be set as "exon" or "intron" to 
#' emphasize "exon" or "intron" region. The "exon" or "intron"
#' can be followed with an integer e.g. "exon_80", or "intron_99".
#' The integer indicates the total percentage of "exon" or "intron" region.
#' Here "exon" indicates all regions in features. 
#' And "intron" indicates all flank regions of the features.
#' @param label_on_feature Labels of the feature directly on them. 
#' Default FALSE.
#' @param lollipop_style_switch_limit The cutoff value for lollipop style for the 'circle' type.
#' If the max score is greater than this cutoff value, trackViewer will only plot one shape at
#' the highest score. Otherwise trackViewer will draw the shapes like `Tanghulu`.
#' @param ... not used.
#' @return NULL
#' @details 
#' In SNP.gr and features, metadata of the GRanges object will be used to control the 
#' color, fill, border, alpha, shape, height, cex, dashline.col, data source of pie if the type is pie. 
#' And also the controls for labels by name the metadata start as 
#' label.parameter.<properties>, and for node labels by name the metadata start as
#' node.label.<properties>,
#' such as label.parameter.rot, label.parameter.gp. The parameter is used for 
#' \link[grid]{grid.text}.or \link[motifStack]{plotMotifLogoA}.
#' The metadata 'featureLayerID' for features are used 
#' for drawing features in different layers. The metadata 'SNPsideID' for SNP.gr
#' are used for determining the side of lollipops. And the 'SNPsideID' could only
#' be 'top' or 'bottom'.
#' @return NULL
#' @import GenomicRanges
#' @import IRanges
#' @import grid
#' @importFrom scales rescale
#' @importClassesFrom grImport Picture
#' @importFrom grImport readPicture grid.picture
#' @importFrom grDevices rgb
#' @export
#' @examples
#' SNP <- c(10, 100, 105, 108, 400, 410, 420, 600, 700, 805, 840, 1400, 1402)
#' x <- sample.int(100, length(SNP))
#' SNP.gr <- GRanges("chr1", IRanges(SNP, width=1, names=paste0("snp", SNP)), 
#'                   value1=x, value2=100-x)
#' SNP.gr$color <- rep(list(c("red", 'blue')), length(SNP))
#' SNP.gr$border <- sample.int(7, length(SNP), replace=TRUE)
#' features <- GRanges("chr1", IRanges(c(1, 501, 1001), 
#'                                     width=c(120, 500, 405),
#'                                     names=paste0("block", 1:3)),
#'                     color="black",
#'                     fill=c("#FF8833", "#51C6E6", "#DFA32D"),
#'                     height=c(0.1, 0.05, 0.08),
#'                     label.parameter.rot=45)
#' lolliplot(SNP.gr, features, type="pie") 
#'

lolliplot <- function(SNP.gr, features=NULL, ranges=NULL,
                      type="circle",
                      newpage=TRUE, ylab=TRUE, ylab.gp=gpar(col="black"),
                      yaxis=TRUE, yaxis.gp=gpar(col="black"), 
                      xaxis=TRUE, xaxis.gp=gpar(col="black"), 
                      legend=NULL, legendPosition='top', cex=1, 
                      dashline.col="gray80", 
                      jitter=c("node", "label"), 
                      rescale=FALSE, 
                      label_on_feature=FALSE,
                      lollipop_style_switch_limit=10,
                      ...){
    stopifnot(inherits(SNP.gr, c("GRanges", "GRangesList", "list")))
    stopifnot(inherits(features, c("GRanges", "GRangesList", "list")))
    jitter <- match.arg(jitter)
    legendPosition <- checkLegendPosition(legendPosition)
    rescale.old <- rescale
    xaxis.old <- xaxis
    if(any(type!="circle"&jitter=="label")){
      jitter[which(type!="circle"&jitter=="label")] <- "node"
      warning("if jitter set to label, type must be cirle.")
      message("jitter is set to node.")
    }
    SNP.gr.name <- deparse(substitute(SNP.gr))
    if(is(SNP.gr, "GRanges")){
        SNP.gr <- list(SNP.gr)
        if(length(SNP.gr.name)==length(SNP.gr)) {
          names(SNP.gr) <- SNP.gr.name
        }
    }
    len <- length(SNP.gr)
    for(i in seq.int(len)){
        stopifnot(is(SNP.gr[[i]], "GRanges"))
    }
    
    ## prepare the feature
    if(inherits(features, c("GRangesList", "list"))){
      for(i in seq_along(features)){
        stopifnot("features must be a GRanges or GRangesList object"=
                    is(features[[i]], "GRanges"))
      }
      features <- features[seq.int(len)]
    }else{
      stopifnot("features must be a GRanges or GRangesList object"=
                  is(features, "GRanges"))
      #features.name <- deparse(substitute(features))
      features <- list(features)[seq.int(len)]
    }
    
    
    TYPES <- c("circle", "pie", "pin", "pie.stack", "flag")
    if(any(!type %in% TYPES)){
        stop("Error in match argument: ",
             paste0("'type' should be one of '",  
                    paste(TYPES, collapse="', '"), "'."))
    }
    types <- rep(type, length=len)[seq.int(len)]
    rm(type)
    ############### handle legend ####################
    ## set the legend as a list, 
    ## if all the legend for different tracks is same
    ## set draw legend for last track later
    legend <- handleLegend(legend, len, SNP.gr)
    
    ################ handle ranges #####################
    ## if !missing(ranges) set ranges as feature ranges
    ranges <- handleRanges(ranges, SNP.gr, features, len)
    
    ##cut all SNP.gr by the range
    SNP.gr <- cutSNP(SNP.gr, ranges, len)
    
    dots <- list(...)
    if(length(dots$maxYlimNchar)>0) {
      maxYlimNchar<- dots$maxYlimNchar
    } else {
      if(all(c(TRUE, vapply(yaxis, isTRUE, logical(1L))))){
        maxYlimNchar <- getMaxYlimNchar(SNP.gr, 3.5, types)
      } else {
        if(length(yaxis)>1){
          maxYlimNchar <- max(c(1, nchar(as.character(yaxis))), na.rm = TRUE) +
            4
        }else{
          maxYlimNchar<- 2.5
        }
      }
    }
    
    ################## plot ############################
    ## total height == 1
    height <- 1/sum(lengths(ranges))
    args <- as.list(match.call())
    if(length(args$height0)==0){
      height0 <- 0
    }else{
      height0 <- args$height0
    }
    if(newpage) grid.newpage()
    for(i in seq.int(len)){
      if(length(ranges[[i]])>1){## split into multilayers
        args$newpage <- FALSE
        for(j in rev(seq_along(ranges[[i]]))){
          args$ranges <- ranges[[i]][j]
          args$SNP.gr <- SNP.gr[i]
          args$features <- features[[i]]
          args$type <- types[i]
          args$legend <- legend[[i]]
          args$legendPosition <- legendPosition
          args$height0 <- height0
          args$maxYlimNchar <- maxYlimNchar
          height0 <- do.call(what = lolliplot, args = args)
        }
      }else{## is GRanges with length==1
        type <- match.arg(types[i], TYPES)
        if(type=="pin"){ ## read the pin shape file
          pinpath <- system.file("extdata", "map-pin-red.xml", package="trackViewer")
          pin <- readPicture(pinpath)
        }else{
          pin <- NULL
        }
        ## Here we don't know the real height of each tracks
        vp0 <- viewport(x=.5, y=height0 + height*0.5,
                       width=1, height=height)
        pushViewport(vp0)
        LINEW <- as.numeric(convertX(unit(1, "line"), "npc"))
        LINEH <- as.numeric(convertY(unit(1, "line"), "npc"))
        totalH <- as.numeric(unit(1, "npc"))
        if(LINEH > totalH/20){
          LINEH <- totalH/20
        }
        
        ## GAP the gaps between any elements
        GAP <- .2 * LINEH
        
        SNPs <- SNP.gr[[i]]
        strand(SNPs) <- "*"
        SNPs <- sort(SNPs)
        
        feature <- features[[i]]
        
        ## rescale
        rescale <- rescale.old
        xaxis <- xaxis.old
        if(is.logical(rescale)[1]){
          if(rescale[1]){
            range.tile <- tile(ranges[[i]], n = 5)[[1]]
            if(all(width(range.tile)>2)){
              range.tile.cnt <- countOverlaps(range.tile, SNPs)
              feature.start <- feature.end <- feature
              end(feature.start) <- start(feature.start)
              start(feature.end) <- end(feature.end)
              range.tile.cnt2 <- countOverlaps(range.tile, unique(c(feature.start, feature.end)))
              range.tile.cnt <- range.tile.cnt + range.tile.cnt2
              range.width <- width(ranges[[i]])
              range.tile.width <- log2(range.tile.cnt + 1)
              range.tile.width <- range.tile.width/sum(range.tile.width)
              range.tile.width <- range.width * range.tile.width
              range.tile.width <- cumsum(range.tile.width)
              range.tile.width <- start(ranges[[i]]) + c(0, round(range.tile.width)-1)
              rescale <- data.frame(from.start=start(range.tile), from.end=end(range.tile),
                                    to.start=range.tile.width[-length(range.tile.width)],
                                    to.end=range.tile.width[-1])
              rescale$to.start[-1] <- rescale$to.start[-1] + 1
            }
          }
        }else{
          if(is.numeric(rescale)){## rescale is a percentage, convert it into from.start, from.end, to.start, to.end
            ## reset the scale for exons and introns
            feature.rd <- disjoin(c(feature, ranges[[i]]))
            feature.segment.points <- sort(unique(c(start(feature.rd), end(feature.rd))))
            ## filter the points by viewer port.
            feature.segment.points <- 
              feature.segment.points[feature.segment.points>=start(ranges[[i]]) &
                                       feature.segment.points<=end(ranges[[i]])]
            rescale <- rescale/sum(rescale, na.rm = TRUE)
            rescale[is.na(rescale)] <- 0.01
            if(length(rescale)==length(feature.segment.points)-1){
              rescale.ir <- IRanges(feature.segment.points[-length(feature.segment.points)]+1, 
                                    feature.segment.points[-1])
              start(rescale.ir)[1] <- start(rescale.ir)[1]-1
              rescale.ir.width <- sum(width(rescale.ir))
              rescale.ir.new.width <- cumsum(round(rescale.ir.width*rescale, digits = 0))
              rescale <- data.frame(from.start=start(rescale.ir), 
                                    from.end=end(rescale.ir),
                                    to.start=feature.segment.points[1] + 
                                      c(0, rescale.ir.new.width[-length(rescale.ir.new.width)]),
                                    to.end=feature.segment.points[1] + rescale.ir.new.width)
            }else{
              stop("The length of rescale is not as same as the number of segments (including features and non-features).")
            }
          }else{
            if(is.character(rescale)){
              if(grepl("exon|intron", rescale[1], ignore.case = TRUE)){
                ## reset the scale for exons and introns
                feature.rd <- disjoin(c(reduce(feature), ranges[[i]]), ignore.strand=TRUE)
                feature.rd <- subsetByOverlaps(feature.rd, ranges[[1]], ignore.strand=TRUE)
                feature.segment.exon <- subsetByOverlaps(feature.rd, feature, ignore.strand=TRUE)
                feature.segment.intron <- subsetByOverlaps(feature.rd, feature, invert=TRUE, ignore.strand=TRUE)
                feature.segment.exon$type <- rep("exon", length(feature.segment.exon))
                feature.segment.intron$type <- rep("intron", length(feature.segment.intron))
                feature.rd <- sort(c(feature.segment.exon, feature.segment.intron))
                feature.segment.points <- sort(unique(c(start(feature.segment.exon), end(feature.segment.exon))))
                ## filter the points by viewer port.
                feature.segment.points <- 
                  feature.segment.points[feature.segment.points>=start(ranges[[i]]) &
                                           feature.segment.points<=end(ranges[[i]])]
                ratio.to.full.range <- 9
                if(grepl("exon", rescale[1], ignore.case = TRUE)){#set intron width to 1/10
                  if(grepl("^exon_\\d+$", rescale[1], ignore.case = TRUE)){
                    ratio.to.full.range <- 
                      as.numeric(sub("^exon_(\\d+)$", "\\1", rescale[1], ignore.case = TRUE))
                    ratio.to.full.range <- ratio.to.full.range/(100-ratio.to.full.range)
                  }
                  width.full.range <- sum(width(feature.rd)[feature.rd$type=="exon"])
                  rescale <- width(feature.rd)
                  rescale[feature.rd$type=="intron"] <- 
                    rep(ceiling(width.full.range/ratio.to.full.range/sum(feature.rd$type=="intron")),
                        length(rescale[feature.rd$type=="intron"]))
                }else{#set exon width to 1/10
                  if(grepl("^intron_\\d+$", rescale[1], ignore.case = TRUE)){
                    ratio.to.full.range <- 
                      as.numeric(sub("^intron_(\\d+)$", "\\1", rescale[1], ignore.case = TRUE))
                    ratio.to.full.range <- ratio.to.full.range/(100-ratio.to.full.range)
                  }
                  width.full.range <- sum(width(feature.rd)[feature.rd$type=="intron"])
                  rescale <- width(feature.rd)
                  rescale[feature.rd$type=="exon"] <- 
                    rep(ceiling(width.full.range/ratio.to.full.range/sum(feature.rd$type=="exon")),
                        length(rescale[feature.rd$type=="exon"]))
                }
                rescale <- rescale/sum(rescale)
                if(length(rescale)==length(feature.segment.points)-1){
                  rescale.ir <- IRanges(feature.segment.points[-length(feature.segment.points)]+1, 
                                        feature.segment.points[-1])
                  start(rescale.ir)[1] <- start(rescale.ir)[1]-1
                  rescale.ir.width <- sum(width(rescale.ir))
                  rescale.ir.new.width <- cumsum(round(rescale.ir.width*rescale, digits = 0))
                  rescale <- data.frame(from.start=start(rescale.ir), 
                                        from.end=end(rescale.ir),
                                        to.start=feature.segment.points[1] + 
                                          c(0, rescale.ir.new.width[-length(rescale.ir.new.width)]),
                                        to.end=feature.segment.points[1] + rescale.ir.new.width)
                }else{
                  stop("Something wrong with the auto-scale setting. Please report the bug to https://github.com/jianhong/trackViewer/issues")
                }
              }
            }
          }
        }
        if(is.data.frame(rescale)){
          if(all(c("from.start", "from.end", "to.start", "to.end") %in% colnames(rescale))){
            ## check the from coverage the whole region.
            checkflank <- function(x){
              to <- IRanges(x$to.start, x$to.end)
              xol <- findOverlaps(to, drop.self=TRUE, 
                                  drop.redundant=TRUE,minoverlap=2L)
              if(length(xol)>1){
                stop("There is overlaps of the rescale region for 'to' columns.")
              }
              x <- IRanges(x$from.start, x$from.end)
              xol <- findOverlaps(x, drop.self=TRUE, 
                                  drop.redundant=TRUE,minoverlap=2L)
              if(length(xol)>1){
                stop("There is overlaps of the rescale region for 'from' columns.")
              }
              xgap <- gaps(x, start=min(start(x)), end=max(end(x)))
              if(length(xgap)>0){
                stop("There is gaps of the rescale region for 'from' columns.")
              }
            }
            checkflank(rescale)
            rescale.gr <- function(x){
              if(is(x, "GRanges")){
                x.start <- start(x)
                x.end <- end(x)
                y <- c(x.start, x.end)
                x.cut <- cut(y, breaks=c(rescale$from.start[1], rescale$from.end+1),
                             labels=seq.int(nrow(rescale)), right=FALSE)
                y <- mapply(function(a, b){
                  if(!is.na(b)) {
                    rescale(a, to=c(rescale$to.start[b], rescale$to.end[b]),
                            from=c(rescale$from.start[b], rescale$from.end[b]))
                  }else{
                    a
                  }
                }, y, as.numeric(as.character(x.cut)))
                y <- round(y)
                start(x) <- 1
                end(x) <- y[seq_along(x)+length(x)]
                start(x) <- y[seq_along(x)]
                x
              }else{
                x.cut <- cut(x, breaks=c(rescale$from.start[1], rescale$from.end+1),
                             labels=seq.int(nrow(rescale)), right=FALSE)
                y <- mapply(function(a, b){
                  if(!is.na(b)) {
                    rescale(a, to=c(rescale$to.start[b], rescale$to.end[b]),
                            from=c(rescale$from.start[b], rescale$from.end[b]))
                  }else{
                    a
                  }
                }, x, as.numeric(as.character(x.cut)))
                y <- round(y)
                y
              }
            }
            feature <- rescale.gr(feature)
            SNPs <- rescale.gr(SNPs)
            if(is.logical(xaxis)[1]){
              if(xaxis[1]){
                xaxis <- c(rescale$to.start[1], rescale$to.end)
                names(xaxis) <- c(rescale$from.start[1], rescale$from.end)
              }
            }else{
              xaxis.names <- names(xaxis)
              if(length(xaxis.names)!=length(xaxis)){
                xaxis.names <- as.character(xaxis)
              }
              xaxis <- rescale.gr(xaxis)
              names(xaxis) <- xaxis.names
            }
          }
        }
        
        ## top legend
        if(legendPosition$position=='bottom'){
          ## not supported yet.
          message('not suppport yet.')
        }
        ## convert height to npc number
        feature$height <- convertHeight2NPCnum(feature$height)
        ## multiple transcripts in one gene could be separated by featureLayerID
        feature <- setFeatureLayerID(feature, ranges[[i]])
        feature.splited <- split(feature, feature$featureLayerID)
        
        ## bottomblank, the transcripts legend height
        bottomblank <- plotFeatureLegend(feature, as.numeric(convertY(unit(1, "line"), "npc")),
                                         ranges[[i]], xaxis, xaxis.gp, label_on_feature)
        
        ## get the max score and scoreType
        if(length(SNPs$score)>0){
          SNPs$score <- sapply(SNPs$score, mean) ## fix the bug of score is NumericList
        }
        scoreMax0 <- scoreMax <- 
          if(length(SNPs$score)>0) ceiling(max(c(SNPs$score, 1), na.rm=TRUE)) else 1
        if(type=="pie.stack") scoreMax <- length(unique(SNPs$stack.factor))
        if(!type %in% c("pie", "pie.stack")){
          scoreType <- 
            if(length(SNPs$score)>0) all(floor(SNPs$score)==SNPs$score) else FALSE
          if(length(yaxis)>1 && is.numeric(yaxis)){
            if(length(names(yaxis))!=length(yaxis)){
              names(yaxis) <- yaxis
            }
            scoreMax0 <- max(yaxis, scoreMax0)
            scoreMax <- max(yaxis, scoreMax)
          }
          if(scoreMax>lollipop_style_switch_limit) {
            SNPs$score <- 10*SNPs$score/scoreMax
            scoreMax <- 10*scoreMax0/scoreMax
            scoreType <- FALSE
          }else{
            scoreMax <- scoreMax0
          }
        }else{
          scoreType <- FALSE
        }
        
        ## if the type is caterpillar, there are lollipop in both sides
        ## plot the bottom lollipops first. And push a new viewport
        
        IsCaterpillar <- length(SNPs$SNPsideID) > 0
        if(IsCaterpillar){
          if(any(is.na(SNPs$SNPsideID)) || 
             !all(SNPs$SNPsideID %in% c('top', 'bottom'))){
            warning("Not all SNPsideID is top or bottom")
            IsCaterpillar <- FALSE
          }
        }
        
        if(IsCaterpillar){
          SNPs.top <- SNPs[SNPs$SNPsideID=='top']
          SNPs.bottom <- SNPs[SNPs$SNPsideID=='bottom']
        }else{
          SNPs.top <- SNPs
          SNPs.bottom <- GRanges()
        }
        if(length(SNPs.bottom)<1) IsCaterpillar <- FALSE
        ## viewport of plot region
        if(!IsCaterpillar){
          bottomblank <- bottomblank
        }
        
        ## decide the legend position
        # right side keep 3 lines TODO: figure out the plot region to decide the right side margin
        plotRegionWidth <- unit(1 - (maxYlimNchar+3)*LINEW, 'npc')
        plotRegionX <- unit(maxYlimNchar*LINEW, 'npc')
        if(legendPosition$position %in% c('left', 'right')){
          if(length(legendPosition$width)>0){
            if(is.unit(legendPosition$width)){
              plotRegionWidth <- plotRegionWidth - legendPosition$width
            }else if(is.numeric(legendPosition$width)){
              if(legendPosition$width>0 && legendPosition$width<1){
                legendPosition$width <- unit(legendPosition$width, 'npc')
                plotRegionWidth <- plotRegionWidth - legendPosition$width
              }else{
                stop('legendPosition$width must be a number in the range of 0 ~ 1')
              }
            }else{
              stop('legendPosition$width must be a number of a unit.')
            }
          }else{
            legendPosition$width <- convertX(
              unit(max(legendNchar(legend[[i]]))+3, 'char'),
              'npc')
            plotRegionWidth <- plotRegionWidth - legendPosition$width
          }
          if(legendPosition$position=='left'){
            plotRegionX <- legendPosition$width + unit(2*LINEW, 'npc')
          }else{
            plotRegionX <- unit(maxYlimNchar*LINEW, 'npc')
          }
        }
        vp_track <- viewport(x=plotRegionX, y=bottomblank/2 + .5, 
                             width= plotRegionWidth,
                             height= 1 - bottomblank,
                             xscale=c(start(ranges[[i]]), end(ranges[[i]])),
                             clip="off",
                             just = 'left')
        pushViewport(vp_track)
        ratio.yx <- getYXratio()
        ## plot xaxis
        bottomHeight <- 0
        if(IsCaterpillar){
          ## total height == maxscore + extension + gap + labels
          bottomHeight <- getHeight(SNPs=SNPs.bottom, 
                                    ratio.yx=ratio.yx, 
                                    LINEW=LINEW, 
                                    GAP=GAP, 
                                    cex=cex, 
                                    type=type,
                                    scoreMax=scoreMax,
                                    level="data&labels")
          #bottomHeight <- bottomHeight/len
          vp <- viewport(y=bottomHeight, just="bottom",
                         xscale=c(start(ranges[[i]]), end(ranges[[i]])))
          pushViewport(vp)
          xaxis.gp$col <- "gray"
          plot_grid_xaxis(xaxis, gp=xaxis.gp)
          popViewport()
        }else{
          plot_grid_xaxis(xaxis, gp=xaxis.gp)
        }
        
        ## the baseline, the center of the first transcript
        baseline <- 
          max(c(feature.splited[[1]]$height/2, 
                .0001)) + 0.2 * LINEH
        baselineN <- 
          max(c(feature.splited[[length(feature.splited)]]$height/2, 
                .0001)) + 0.2 * LINEH
        
        ##plot features
        feature.height <- plotFeatures(feature.splited, LINEH, bottomHeight, label_on_feature)
        
        if(length(SNPs.bottom)>0){
          plotLollipops(SNPs.bottom, feature.height, bottomHeight, baselineN, 
                        type, ranges[[i]], yaxis, yaxis.gp, scoreMax, scoreMax0, scoreType, 
                        LINEW, cex, ratio.yx, GAP, pin, dashline.col,
                        side="bottom", jitter=jitter)
        }
        feature.height <- feature.height + 2*GAP
        if(length(SNPs.top)>0){
          plotLollipops(SNPs.top, feature.height, bottomHeight, baseline, 
                        type, ranges[[i]], yaxis, yaxis.gp, scoreMax, scoreMax0, scoreType, 
                        LINEW, cex, ratio.yx, GAP, pin, dashline.col,
                        side="top", jitter=jitter)
        }
        
        ## top legend
        this.height <- getHeight(SNPs.top, 
                                 ratio.yx, LINEW, GAP, cex, type,
                                 scoreMax=scoreMax,
                                 level="data&labels")
        this.height0 <- this.height <- this.height + bottomHeight + feature.height
        if(legendPosition$position=='top'){
          this.height <- plotLegend(legend[[i]], this.height, LINEH)
          if('alpha' %in% names(legend[[i]])){
            legend[[i]]$alpha <- NULL
            if('pch' %in% names(legend[[i]])){
              legend[[i]]$pch <- NA
            }
            plotLegend(legend[[i]], this.height0, LINEH)
          }
        }
        popViewport()## vp_track
        
        this.height <- bottomblank + 
          this.height * (1 - bottomblank)
        
        ## ylab
        vp_ylab <- viewport(x=ifelse(legendPosition$position!='left', 0, 
                                     convertX(legendPosition$width, 'npc',
                                              valueOnly = TRUE)),
                            y=this.height*0.5, 
                            width=plotRegionWidth +
                              unit(2*maxYlimNchar*LINEW, 'npc'),
                            height=this.height,
                            just='left')
        pushViewport(vp_ylab)
        if(is.logical(ylab)){
          if(ylab && length(names(SNP.gr))>0){
            grid.text(names(SNP.gr)[i], x = LINEW, 
                      y = .5, rot = 90, gp=ylab.gp)
          }
        }
        if(is.character(ylab)){
          if(length(ylab)==1) ylab <- rep(ylab, len)
          grid.text(ylab[i], x = LINEW,
                    y = .5, rot = 90, gp=ylab.gp)
        }
        popViewport()#vp_ylab
        
        ## left or right legend
        if(legendPosition$position %in% c('left', 'right')){
          legendRegionWidth <- legendPosition$width
          if(legendPosition$position=='left'){
            legendRegionX <- unit(0, 'npc')
          }else{
            legendRegionX <- unit(1, 'npc') - legendPosition$width
          }
          vp_legend <- viewport(x=legendRegionX, y=.5, 
                                width= legendRegionWidth,
                                height= 1,
                                just = 'left')
          pushViewport(vp_legend)
          plotLegend(legend[[i]], 0, LINEH)
          popViewport()## vp_legend
        }
        
        popViewport()#vp0
        height0 <-  height0 + this.height*height
      }
    }
    return(invisible(height0))
}