R/genomic.R
In circlize: Circular Visualization

Documented in circos.genomicAxis circos.genomicDensity circos.genomicHeatmap circos.genomicIdeogram circos.genomicInitialize circos.genomicLabels circos.genomicLines circos.genomicLink circos.genomicPoints circos.genomicPosTransformLines circos.genomicRainfall circos.genomicRect circos.genomicText circos.genomicTrack circos.genomicTrackPlotRegion circos.initializeCircularGenome circos.initializeWithIdeogram genomicDensity get.current.chromosome getI highlight.chromosome posTransform.default posTransform.text rainfallTransform smartAlign

# == title
# Initialize the circular layout with an ideogram
#
# == param
# -cytoband  A path of the cytoband file or a data frame that already contains cytoband data. By default it is cytoband for hg19.
#            Pass to `read.cytoband`.
# -species Abbreviations of species. e.g. hg19 for human, mm10 for mouse. If this
#          value is specified, the function will download cytoBand.txt.gz from
#          UCSC website automatically. If there is no cytoband for user's species,
#          it will keep on trying to download chromInfo file. Pass to `read.cytoband` or `read.chromInfo`.
# -chromosome.index subset of chromosomes, also used to reorder chromosomes.
# -sort.chr Whether chromosome names should be sorted (first sort by numbers then by letters).
#           If ``chromosome.index`` is set, this argumetn is enforced to ``FALSE``
# -major.by     Increment of major ticks. Pass to `circos.genomicInitialize`.
# -plotType     Which tracks should be drawn. ``ideogram`` for ideogram rectangle, ``axis`` for genomic axis and ``labels`` for chromosome names.
#               If there is no ideogram for specified species, ``ideogram`` will be enforced to be excluded.
#               If it is set to ``NULL``, the function just initialize the plot but draw nothing.
# -track.height Height of the track which contains "axis" and "labels".
# -ideogram.height Height of the ideogram track
# -...    Pass to `circos.genomicInitialize`.
#
# == details
# The function will initialize the circular plot in which each sector corresponds to a chromosome. You can control the order of 
# chromosomes by ``chromosome.index`` or by ``sort.chr``, or by setting a special format of ``cytoband`` (please refer to `read.cytoband` 
# to find out how to control a proper ``cytoband``).
#
# The function finally pass data to `circos.genomicInitialize` to initialize the circular plot.
#
# The style of ideogram is almost fixed, but you can customize it with your self-sefined code. Refer to vignette for demonstration.
#
# == seealso
# https://jokergoo.github.io/circlize_book/book/initialize-genomic-plot.html#initialize-cytoband
#
# == example
# \donttest{
# circos.initializeWithIdeogram()
#
# cytoband.file = system.file(package = "circlize",
#     "extdata", "cytoBand.txt")
# circos.initializeWithIdeogram(cytoband.file)
#
# cytoband.df = read.table(cytoband.file, colClasses = c("character", "numeric",
#     "numeric", "character", "character"), sep = "\t")
# circos.initializeWithIdeogram(cytoband.df)
#
# circos.initializeWithIdeogram(species = "hg18")
#
# circos.initializeWithIdeogram(species = "mm10")
#
# circos.initializeWithIdeogram(chromosome.index = c("chr1", "chr2"))
#
# cytoband = read.table(cytoband.file, colClasses = c("character", "numeric",
#     "numeric", "character", "character"), sep = "\t")
# circos.initializeWithIdeogram(cytoband, sort.chr = FALSE)
#
# cytoband[[1]] = factor(cytoband[[1]], levels = paste0("chr", c(22:1, "X", "Y")))
# circos.initializeWithIdeogram(cytoband, sort.chr = FALSE)
#
# cytoband = read.table(cytoband.file, colClasses = c("character", "numeric",
#     "numeric", "character", "character"), sep = "\t")
# circos.initializeWithIdeogram(cytoband, sort.chr = TRUE)
#
# circos.initializeWithIdeogram(plotType = c("axis", "labels"))
#
# circos.initializeWithIdeogram(plotType = NULL)
#
# circos.par("start.degree" = 90)
# circos.initializeWithIdeogram()
# circos.clear()
#
# circos.par("gap.degree" = rep(c(2, 4), 12))
# circos.initializeWithIdeogram()
# circos.clear()
# }
circos.initializeWithIdeogram = function(
	cytoband = system.file(package = "circlize", "extdata", "cytoBand.txt"), 
	species = NULL, 
	sort.chr = TRUE,
	chromosome.index = usable_chromosomes(species), 
	major.by = NULL,
	plotType = c("ideogram", "axis", "labels"), 
	track.height = NULL, 
	ideogram.height = convert_height(2, "mm"), 
	...) {
	
	if("sector.names" %in% names(list(...))) {
		stop_wrap("Found argumnet `sector.names` is used. Please use the argument `chromosome.index` instead.")
	}

	# proper order will be returned depending on cytoband and sort.chr
	e = try(cytoband <- read.cytoband(cytoband, species = species, sort.chr = sort.chr, chromosome.index = chromosome.index), silent = TRUE)
	if(inherits(e, "try-error") && !is.null(species)) {  # if species is defined
		e2 = try(cytoband <- read.chromInfo(species = species, sort.chr = sort.chr, chromosome.index = chromosome.index), silent = TRUE)
		if(inherits(e2, "try-error")) {
			message(e)
			message(e2)
			stop_wrap("Cannot download either cytoband or chromInfo file from UCSC.")
		} else {
			message_wrap("Downloading cytoBand file from UCSC failed. Use chromInfo file instead. Note ideogram track will be removed from the plot.")
			plotType = setdiff(plotType, "ideogram")

			# because in chromInfo file, there are also many short scaffold
			if(is.null(chromosome.index)) {
	            chromInfo = read.chromInfo(species = species)
	            chr_len = sort(chromInfo$chr.len, decreasing = TRUE)

	            # sometimes there are small scaffold
	            i = which(chr_len[seq_len(length(chr_len)-1)] / chr_len[seq_len(length(chr_len)-1)+1] > 5)[1]
	            if(length(i)) {
	                chromosome = chromInfo$chromosome[chromInfo$chromosome %in% names(chr_len[chr_len >= chr_len[i]])]
	            } else {
	                chromosome = chromInfo$chromosome
	            }
	            cytoband = read.chromInfo(species = species, chromosome.index = chromosome, sort.chr = sort.chr)
	        } 
		}
	} else if(inherits(e, "try-error")) {
		stop(e)
	}
	df = cytoband$df
	chromosome = cytoband$chromosome
	if(is.null(chromosome)) {
		if(is.factor(cytoband[, 1])) {
			chromosome = levels(cytoband$df[, 1])
		} else {
			chromosome = unique(cytoband$df[, 1])
		}
	}

	if(is.null(chromosome.index)) {
		chromosome.index = chromosome
	}

	# here df[[1]] is quite important, should be re-factered
	df[[1]] = factor(as.vector(df[[1]]), levels = chromosome.index)
	
	# sn for sector names, but not for sector index
	sn = unique(as.vector(df[[1]]))

	# we do not need 'chr' prefix if it exits, it holds too much space.
	sn = gsub("chr", "", sn)
	
	o.cell.padding = circos.par("cell.padding")
	circos.par(cell.padding = c(o.cell.padding[1], 0, o.cell.padding[3], 0))
	
	circos.genomicInitialize(df, sector.names = sn, major.by = major.by, plotType = plotType, track.height = track.height, ...)

	if(any(plotType %in% "ideogram")) {
		circos.genomicIdeogram(df, track.height = ideogram.height)
	}
}

# == title
# Add an ideogram track
#
# == param
# -cytoband A data frame or a file path, pass to `read.cytoband`.
# -species Abbreviations of the genome, pass to `read.cytoband`.
# -track.height Height of the ideogram track.
# -track.margin Margins for the track.
#
# == seealso
# https://jokergoo.github.io/circlize_book/book/high-level-genomic-functions.html#ideograms
#
# == author
# Zuguang Gu <z.gu@dkfz.de>
#
# == example
# \donttest{
# circos.initializeWithIdeogram(plotType = c("labels", "axis"))
# circos.track(ylim = c(0, 1))
# circos.genomicIdeogram() # put ideogram as the third track
# }
circos.genomicIdeogram = function(
	cytoband = system.file(package = "circlize", "extdata", "cytoBand.txt"), 
	species = NULL, 
	track.height = mm_h(2),
	track.margin = circos.par("track.margin")) {

	chromosome.index = get.all.sector.index()
	e = try(cytoband <- read.cytoband(cytoband, species = species, chromosome.index = chromosome.index), silent = TRUE)
	if(inherits(e, "try-error")) {
		stop(e)
	}
	df = cytoband$df
	if(all(cytoband.col(df[, 5]) == "#FFFFFF")) {
		warning_wrap("Cannot map colors to cytobands. The ideogram won't be drawn.")
		return(invisible(NULL))
	}
	circos.genomicTrackPlotRegion(df, ylim = c(0, 1), bg.border = NA, track.height = track.height,
		panel.fun = function(region, value, ...) {
			col = cytoband.col(value[[2]])
			circos.genomicRect(region, value, ybottom = 0, ytop = 1, col = col, border = NA, ...)
			xlim = get.cell.meta.data("xlim")
			circos.rect(xlim[1], 0, xlim[2], 1, border = "black")
		}, cell.padding = c(0, 0, 0, 0), track.margin = track.margin
	)
}

# == title
# Initialize circular plot with any genomic data
#
# == param
# -data         A data frame in bed format.
# -sector.names Labels for each sectors which will be drawn along each sector. It will not modify values of sector index.
# -major.by     Increment of major ticks. It is calculated automatically if the value is not set (about every 10 degrees there is a major tick).
# -plotType     If it is not ``NULL``, there will create a new track containing axis and names for sectors.
#               This argument controls which part should be drawn, ``axis`` for genomic axis and ``labels`` for chromosome names
# -tickLabelsStartFromZero Whether axis tick labels start from 0? This will only affect the axis labels while not affect x-values in cells.
# -axis.labels.cex The font size for the axis tick labels.
# -labels.cex   The font size for the labels.
# -track.height If ``PlotType`` is not ``NULL``, height of the annotation track.
# -...          Pass to `circos.initialize`
#
# == details
# The function will initialize circular plot from genomic data. If ``plotType`` is set with value in ``axis`` or ``labels``, there will
# create a new track.
#
# The order of sectors related to data structure of ``data``. If the first column in ``data`` is a factor, the order of sectors
# is ``levels(data[[1]])``; If the first column is just a simple vector, the order of sectors is ``unique(data[[1]]``.
#
# For more details on initializing genomic plot, please refer to the vignettes.
#
# == seealso
# https://jokergoo.github.io/circlize_book/book/initialize-genomic-plot.html#initialize-with-general-genomic-category
#
# == example
# df = read.cytoband()$df
# circos.genomicInitialize(df)
#
# df = data.frame(name = c("TP53", "TP63", "TP73"),
#                 start = c(7565097, 189349205, 3569084),
#                 end = c(7590856, 189615068, 3652765),
#                 stringsAsFactors = FALSE)
# circos.genomicInitialize(df)
# circos.clear()
#
# circos.genomicInitialize(df, tickLabelsStartFromZero = FALSE)
# circos.clear()
#
# circos.genomicInitialize(df, major.by = 5000)
# circos.clear()
#
# circos.genomicInitialize(df, plotType = "labels")
# circos.clear()
#
# circos.genomicInitialize(df, sector.names = c("tp53", "tp63", "tp73"))
# circos.clear()
#
# circos.genomicInitialize(df, sector.names = c("tp53x", "tp63x", "tp73"))
# circos.clear()
#
# df[[1]] = factor(df[[1]], levels = c("TP73", "TP63", "TP53"))
# circos.genomicInitialize(df)
# circos.clear()
#
circos.genomicInitialize = function(
	data, 
	sector.names = NULL, 
	major.by = NULL,
	plotType = c("axis", "labels"), 
	tickLabelsStartFromZero = TRUE,
	axis.labels.cex = 0.4*par("cex"), 
	labels.cex = 0.8*par("cex"), 
	track.height = NULL, 
	...) {

	data = validate_data_frame(data)
	validate_region(data, check_chr = FALSE)
	
	if(is.factor(data[[1]])) {
		fa = levels(data[[1]])
	} else {
		fa = unique(data[[1]])
	}
	
	if(!is.null(sector.names)) {
		if(length(sector.names) != length(fa)) {
			stop_wrap("length of `sector.names` and length of sectors differ.")
		}
	} else {
		sector.names = fa
	}
	
	names(sector.names) = fa
	
	# calculate xlim
	x1 = tapply(data[[2]], data[[1]], min)[fa]
	x2 = tapply(data[[3]], data[[1]], max)[fa]
	
	op = circos.par("cell.padding")
	ow = circos.par("points.overflow.warning")
	circos.par(cell.padding = c(0, 0, 0, 0), points.overflow.warning = FALSE)
	circos.initialize(factor(fa, levels = fa), xlim = cbind(x1, x2), ...)

	if(circos.par$ring) {
		op = c(op[1], 0, op[3], 0)
    	ow = FALSE
	}

	if(is.null(track.height)) {
		if(all(c("axis", "labels") %in% plotType)) {
			track.height = convert_unit_in_canvas_coordinate(1.5, "mm") + strheight("0", cex = axis.labels.cex) + 
				convert_unit_in_canvas_coordinate(0.5, "mm") + strheight("chr", cex = labels.cex)
		} else if("labels" %in% plotType) {
			track.height = strheight("chr", cex = labels.cex)
		} else if("axis" %in% plotType) {
			track.height = convert_unit_in_canvas_coordinate(1.5, "mm") + strheight("0", cex = axis.labels.cex)
		} else {
			track.height = convert_height(3, "mm")
		}
	}

	# axis and chromosome names
	if(any(plotType %in% c("axis", "labels"))) {
		circos.genomicTrackPlotRegion(data, ylim = c(0, 1), bg.border = NA, track.height = track.height,
			panel.fun = function(region, value, ...) {
				sector.index = get.cell.meta.data("sector.index")
				xlim = get.cell.meta.data("xlim")
				
				if(all(c("axis", "labels") %in% plotType)) {
					circos.genomicAxis(h = "bottom", major.by = major.by, tickLabelsStartFromZero = tickLabelsStartFromZero, labels.cex = axis.labels.cex)
					circos.text(mean(xlim), convert_y(1.5, "mm") + convert_y(strheight("chr", cex = axis.labels.cex), "canvas") + convert_y(0.5, "mm"), 
						labels = sector.names[sector.index], cex = labels.cex, adj = c(0.5, 0), niceFacing = TRUE)
				} else if("labels" %in% plotType) {
					circos.text(mean(xlim), 0, labels = sector.names[sector.index], cex = labels.cex, adj = c(0.5, 0), niceFacing = TRUE)
				} else if("axis" %in% plotType) {
					circos.genomicAxis(h = "bottom", major.by = major.by, tickLabelsStartFromZero = tickLabelsStartFromZero,labels.cex = axis.labels.cex)
				}
 			}
		)
	}
	
	circos.par("cell.padding" = op, "points.overflow.warning" = ow)
	return(invisible(NULL))
}


# == title
# Initialize a layout for circular genome
#
# == param
# -name Name of the genome (or the "chromosome name").
# -genome_size Size of the genome
# -plotType Pass to `circos.genomicInitialize`.
# -... All goes to `circos.genomicInitialize`.
#
circos.initializeCircularGenome = function(name, genome_size, plotType = "axis", ...) {
	circos.genomicInitialize(data.frame(name, 0, genome_size), ..., ring = TRUE, plotType = plotType)
} 

# == title
# Add genomic axes
#
# == param
# -h Position of the axes. "top" or "bottom".
# -major.at Major breaks. If ``major.at`` is set, ``major.by`` is ignored.
# -labels labels corresponding to ``major.at``. If ``labels`` is set, ``major.at`` must be set.
# -major.by Increment of major ticks. It is calculated automatically if the value is not set (about every 10 degrees there is a major tick).
# -tickLabelsStartFromZero Whether axis tick labels start from 0? This will only affect the axis labels while not affect x-values in cells.
# -labels.cex The font size for the axis tick labels.
# -sector.index Index for the sector
# -track.index  Index for the track
# -... Other arguments pass to `circos.axis`.
#
# == details
# It assigns proper tick labels under genomic coordinate.
# 
# == seealso
# https://jokergoo.github.io/circlize_book/book/high-level-genomic-functions.html#genomic-axes
#
# == example
# circos.initializeWithIdeogram(chromosome.index = paste0("chr", 1:4), plotType = NULL)
# circos.track(ylim = c(0, 1), panel.fun = function(x, y) circos.genomicAxis())
# circos.track(ylim = c(0, 1), track.height = 0.1)
# circos.track(track.index = get.current.track.index(), panel.fun = function(x, y) {
#     circos.genomicAxis(h = "bottom", direction = "inside")
# })
# circos.clear()
#
circos.genomicAxis = function(
	h = "top", 
	major.at = NULL, 
	labels = NULL,
	major.by = NULL, 
	tickLabelsStartFromZero = TRUE,
	labels.cex = 0.4*par("cex"), 
	sector.index = get.current.sector.index(),
	track.index = get.current.track.index(), 
	...) {

	if(!h %in% c("top", "bottom")) {
		stop_wrap("`h` can only be 'top' or 'bottom'.")
	}

	xlim = get.cell.meta.data("xlim", sector.index = sector.index, track.index = track.index)

	if(!is.null(major.at) && !is.null(labels)) {
		if(length(major.at) != length(labels)) {
			stop_wrap("Length of major.at and labels should be the same.")
		}
	}
	if(is.null(major.at) && !is.null(labels)) {
		stop_wrap("If labels is set, major.at should also be set.")
	}

	if(tickLabelsStartFromZero) {
		offset = xlim[1]
		if(is.null(major.by)) {
			major.by = .default.major.by()
		}

		if(is.null(major.at)) {
			major.at = seq(xlim[1], xlim[2], by = major.by)
			major.at = c(major.at, major.at[length(major.at)] + major.by)
		}
		
		if(is.null(labels)) {
			if(major.by >= 1e6) {
				major.tick.labels = paste((major.at-offset)/1000000, "Mb", sep = "")
			} else if(major.by >= 1e3) {
				major.tick.labels = paste((major.at-offset)/1000, "kb", sep = "")
			} else {
				major.tick.labels = paste((major.at-offset), "bp", sep = "")
			}
		} else {
			major.tick.labels = labels
		}
		
	} else {
		if(is.null(major.by)) {
			major.by = .default.major.by()
		}
		if(is.null(major.at)) {
			major.at = seq(floor(xlim[1]/major.by)*major.by, xlim[2], by = major.by)
			major.at = c(major.at, major.at[length(major.at)] + major.by)
		}

		if(is.null(labels)) {
			if(major.by >= 1e6) {
				major.tick.labels = paste(major.at/1000000, "MB", sep = "")
			} else if(major.by >= 1e3) {
				major.tick.labels = paste(major.at/1000, "KB", sep = "")
			} else {
				major.tick.labels = paste(major.at, "bp", sep = "")
			}
		} else {
			major.tick.labels = labels
		}
	}
	circos.axis(h = h, major.at = major.at, labels = major.tick.labels, labels.cex = labels.cex,
		sector.index = sector.index, track.index = track.index, ...)				
}

# == title
# Create a track for genomic graphics
#
# == param
# -data A bed-file-like data frame or a list of data frames
# -ylim If it is ``NULL``, the value will be calculated from data. If ``stack`` is set to ``TRUE``, this value is ignored.
# -stack whether to plot in a "stack" mode.
# -numeric.column Columns of numeric values in ``data`` that will be used for plotting. 
#                 If ``data`` is a data frame list, ``numeric.column`` should be either length of one or length of ``data``.
#                 If value of ``numeric.column`` is not set, its value will depend on the structure of ``data``.
#                 If ``data`` is a data frame, the default value for ``numeric.column`` is all the numeric column starting from the fourth column.
#                 If ``data`` is a list of data frame, the default value for ``numeric.column`` is a vector which have the same length as ``data``
#                 and the value in default ``numeric.column`` is the index of the first numeric column in corresponding data frame.
# -jitter Numeric. Only works for adding points in ``circos.genomicTrackPlotRegion`` under ``stack`` mode
# -panel.fun Self-defined function which will be applied on each sector. Please not it is different
#            from that in `circos.trackPlotRegion`. In this function, there are two arguments (``region`` and ``value``) plus ``...``.
#            In them, ``region`` is a two-column data frame with start positions and end positions in current genomic category (e.g. chromosome). 
#            ``value`` is a data frame which is derived from ``data`` but excluding the first three columns. Rows in ``value`` correspond to 
#            rows in ``region``. ``...`` is mandatory and is used to pass internal parameters to other functions. The definition of
#            ``value`` will be different according to different input data (data frame or list of data frame) and different settings (stacked or not), 
#            please refer to 'details' section and vignettes to detailed explanation.
# -... Pass to `circos.trackPlotRegion`.
#
# == details
# Similar as `circos.trackPlotRegion`, users can add customized graphics by ``panel.fun``, but the behaviour of ``panel.fun``
# will change depending on users' input data and ``stack`` setting.
#
# When ``data`` is a single data frame, ``region`` in ``panel.fun`` is a data frame containing the second and third column in ``data`` in 'current` genomic category (e.g. current chromosome).
# ``value`` is also a data frame containing columns in ``data`` excluding the first three columns.
#
# When ``data`` is a list containing data frames, ``panel.fun`` will be applied iteratively on each data frame, thus, 
# ``region`` is extracted from the data frame which is in the current iteration. For example, if ``data`` contains two data frames, ``panel.fun``
# will be applied with the first data frame in current chromosome and then applied with the second data frame in the same chromosome.
#
# If ``stack`` is set to ``TRUE``, ``ylim`` will be re-defined. in ``stack`` mode, the y-axis will be splitted into several part
# with equal height and graphics will be drawn on each 'horizontal' lines (y = 1, 2, ...). In this case:
#
# When ``data`` is a single data frame containing one or more numeric columns, each numeric column defined in ``numeric.column`` will be treated as a single unit. 
# ``ylim`` is re-defined to ``c(0.5, n+0.5)`` in which ``n`` is number of numeric columns. ``panel.fun`` will be applied iteratively on each numeric column. In each
# iteration, in ``panel.fun``, ``region`` is still the genomic regions in current genomic category, but ``value`` contains current numeric column plus all non-numeric columns.
# Under ``stack`` mode, in ``panel.fun``, all low-level genomic graphical functions will draw on the 'horizontal line' ``y = i`` in which ``i`` is the index of current numeric column 
# and the value of ``i`` can be obtained by `getI`.
#
# When ``data`` is a list containing data frames, each data frame will be treated as a single unit. The situation is quite similar as described in previous paragraph.
# ``ylim`` is re-defined to ``c(0.5, n+0.5)`` in which ``n`` is number of data frames. ``panel.fun`` will be applied iteratively on each data frame. In each
# iteration, in ``panel.fun``, ``region`` is still the genomic regions in current genomic category, and ``value`` contains columns in current data frame excluding the first three columns.
# Under ``stack`` mode, in ``panel.fun``, all low-level genomic graphical functions will draw on the 'horizontal line' ``y = i`` in which ``i`` is the index of current data frame.
#
# Being different from ``panel.fun`` in `circos.trackPlotRegion`, there should be an additional argument ``...`` in ``panel.fun``. This additional
# argument is used to pass hidden values to low-level graphical functions. So if you are using functions like ``circos.genomicPoints``, you should also
# add ``...`` as an additional argument into ``circos.genomicPoints``.
#
# == seealso
# https://jokergoo.github.io/circlize_book/book/genomic-plotting-region.html and https://jokergoo.github.io/circlize_book/book/modes-of-input.html
#
circos.genomicTrackPlotRegion = function(
	data = NULL, 
	ylim = NULL, 
	stack = FALSE,
	numeric.column = NULL, 
	jitter = 0,
	panel.fun = function(region, value, ...) {NULL}, 
	... ) {
	
	if(is.null(data)) {
		all.sector.index = get.all.sector.index()
		data = data.frame(all.sector.index,
		           rep(0, length(all.sector.index)),
				   rep(0, length(all.sector.index)))
	}

    if(is.function(data) || is.function(ylim) || is.function(stack) || is.function(numeric.column)) {
        stop_wrap("The panel function should be set explicitly with the argument name `panel.fun = ...`.")
    }
	
	# re-define panel.fun
	genomicPanelFun = panel.fun

	# now `data` is either a data frame or a list of data frame
	data = normalizeToDataFrame(data)

	if(is.dataFrameList(data)) {
		for(i in seq_along(data)) {
			data[[i]][[1]] = as.character(data[[i]][[1]])
			if(circos.par$ring) {
				l = data[[i]][, 2] > data[[i]][, 3]
				if(any(l)) {
					data[[i]][l, 2] = data[[i]][l, 2] - diff(get.sector.data()[c("min.value", "max.value")])
				}
			}
			validate_region(data[[i]], check_chr = FALSE)
		}
	} else {
		data[[1]] = as.character(data[[1]])
		if(circos.par$ring) {
			l = data[, 2] > data[, 3]
			if(any(l)) {
				data[l, 2] = data[l, 2] - diff(get.sector.data()[c("min.value", "max.value")])
			}
		}
		validate_region(data, check_chr = FALSE)
	}

	# excluding the first three columns
	if(!is.null(numeric.column)) {
		if(is.numeric(numeric.column)) {
			numeric.column = numeric.column - 3
		}
		if(any(numeric.column <= 0)) {
			stop_wrap("Wrong value in `numeric.column`, they should be larger than 3 or character index.")
		}
	}

	# auto calcualte numeric column
	if(is.dataFrameList(data)) {
		if(!is.null(numeric.column)) {
			if(length(numeric.column) == 1) {
				numeric.column = rep(list(numeric.column), length(data))
			} else if(length(numeric.column) == length(data)) {
				
			} else {
				stop_wrap("Length of `numeric.column` should only be one or length of ``data`` if it is a list of data frames.")
			}
			for(i in seq_along(data)) {
				if(!is.numeric(data[[i]][-(1:3)][, numeric.column[i]])) {
					stop_wrap("Some of your `numeric.column` are not numeric.")
				}
			}
		} else {
			numeric.column = sapply(data, function(gr) {
								nc = which(as.logical(sapply(gr[-(1:3)], is.numeric)))
								if(length(nc) == 0) {
									NA
								} else {
									nc[1]
								}
							})
		}
		
	} else {
		# check numeric.column
		if(is.null(numeric.column)) {
			numeric.column = which(as.logical(sapply(data[-(1:3)], is.numeric)))
		} else {
			if(!all(sapply(data[-(1:3)][, numeric.column, drop = FALSE], is.numeric))) {
				stop_wrap("Some of your `numeric.column` are not numeric.")
			}
		}
	}
	
	args = formals(genomicPanelFun)
	if(!(length(args) == 3 && names(args)[3] == "...")) {
		stop_wrap("The `panel.fun` need a third argument `...` to pass special parameters to graphical functions.")
	}

	if(stack) {
	
		if(is.dataFrameList(data)) {
			n = length(data)
			
			circos.trackPlotRegion(ylim = c(0.5, n + 0.5), panel.fun = function(x, y) {
						
					chr = get.current.chromosome()
					for(i in seq_len(n)) {
						l = data[[i]][[1]] == chr
						df = data[[i]][l, , drop = FALSE]
						if(nrow(df)) {
							.param = new.env()
							assign("i", i, envir = .param)
							assign("stack", TRUE, envir = .param)
							assign("jitter", jitter, envir = .param)
							if(!is.null(numeric.column) && !is.na(numeric.column[i])) {
								assign("numeric.column", numeric.column[i], envir = .param)
							}
							genomicPanelFun(df[2:3], df[-(1:3)], .param = .param)
						}
					}

				}, ...)
		} else {
			n = length(numeric.column)
			non.numeric.column = setdiff(seq_along(data[-(1:3)]), numeric.column)
			
			# if there is no numeric column
			if(n == 0) {
				circos.trackPlotRegion(ylim = c(0.5, 1 + 0.5), panel.fun = function(x, y) {
							
						chr = get.current.chromosome()
						l = data[[1]] == chr
						df = data[l, , drop = FALSE]
						i = 1
						if(nrow(df)) {
							.param = new.env()
							assign("i", i, envir = .param)
							assign("stack", TRUE, envir = .param)
							assign("jitter", jitter, envir = .param)
							genomicPanelFun(df[2:3], df[-(1:3)][non.numeric.column], .param = .param)
						}

					}, ...)
			} else {
				circos.trackPlotRegion(ylim = c(0.5, n + 0.5), panel.fun = function(x, y) {
							
						chr = get.current.chromosome()
						for(i in seq_len(n)) {
							l = data[[1]] == chr
							df = data[l, , drop = FALSE]
							if(nrow(df)) {
								.param = new.env()
								assign("i", i, envir = .param)
								assign("stack", TRUE, envir = .param)
								assign("numeric.column", 1, envir = .param)
								assign("jitter", jitter, envir = .param)
								genomicPanelFun(df[2:3], df[-(1:3)][c(numeric.column[i], non.numeric.column)], .param = .param)
							}
						}

					}, ...)
			}
		}			
		
	} else {
		# auto calculate ylim
		if(is.null(ylim)) {
			if(is.dataFrameList(data)) {
				ylim = range(unlist(lapply(seq_along(data), function(i) {
					gr = data[[i]]
					if(is.na(numeric.column[i])) {
						stop_wrap("There is no numeric column in one of your data frame which calculation of `ylim` depends on. Or you can set `ylim` explicitely.")
					}
					range(unlist(lapply(gr[-(1:3)][ numeric.column[i] ], range, na.rm = TRUE)))
				})))
			} else {
				if(length(numeric.column) == 0) {
					stop_wrap("There is no numeric column in your data frame which calculation of `ylim` depends on. Or you can set `ylim` explicitely.")
				}
				ylim = range(unlist(lapply(data[-(1:3)][numeric.column], range, na.rm = TRUE)))
			}

			if(ylim[1] == ylim[2]) {
				stop_wrap("It seems the data points are all the same. Please explicitly set values to `ylim`.")
			}
		}

		if(is.dataFrameList(data)) {
			circos.trackPlotRegion(ylim = ylim, panel.fun = function(x, y) {
					
					chr = get.current.chromosome()
					for(i in seq_along(data)) {
						l = data[[i]][, 1] == chr
						df = data[[i]][l, , drop = FALSE]
						if(nrow(df)) {
							.param = new.env()
							assign("i", i, envir = .param)
							if(!is.na(numeric.column[i])) {
								assign("numeric.column", numeric.column[i], envir = .param)
							}
							genomicPanelFun(df[2:3], df[-(1:3)], .param = .param)
						}
					}
				}, ...)
		} else {
			circos.trackPlotRegion(ylim = ylim, panel.fun = function(x, y) {
					
					chr = get.current.chromosome()
					df = data[data[[1]] == chr, , drop = FALSE]
					if(nrow(df)) {
						.param = new.env()
						assign("i", 1, envir = .param)
						if(length(numeric.column)) {
							assign("numeric.column", numeric.column, envir = .param)
						}
						genomicPanelFun(df[2:3], df[-(1:3)], .param = .param)
					}
				}, ...)
		}
	}
	
}

# == title
# Create a track for genomic graphics
#
# == param
# -... Pass to `circos.genomicTrackPlotRegion`.
#
# == details
# shortcut function of `circos.genomicTrackPlotRegion`.
#
circos.genomicTrack = function(...) {
	circos.genomicTrackPlotRegion(...)
}


# == title
# Which data that ``panel.fun`` is using
#
# == param
# -... Invisible arguments that users do not need to care
#
# == details
# The function should only be put inside ``panel.fun`` when using `circos.genomicTrackPlotRegion`.
#
# If ``stack`` is set to ``TRUE`` in `circos.genomicTrackPlotRegion`, the returned value
# indicates which stack the function will be applied to.
#
# If ``data`` is a list of data frames, the value
# indicates which data frame is being used. Please see the vignette to get a more clear explanation.
getI = function(...) {
	args = list(...)
	if(is.null(args$.param)) {
		stop_wrap("Maybe you should call like `getI(...)`")
	}
	.param = args$.param
	return(.param$i)
}


# == title
# Add points to a plotting region, specifically for genomic graphics
#
# ==param
# -region A data frame contains 2 columns which correspond to start positions and end positions.
# -value  A data frame contains values and other information.
# -numeric.column Which column in ``value`` data frame should be taken as y-value.
#                 If it is not defined, the whole numeric columns in ``value`` will be taken.
# -sector.index Index of sector.
# -track.index Index of track.
# -posTransform Self-defined function to transform genomic positions, see `posTransform.default` for explanation
# -col Color of points. If there is only one numeric column, the length of ``col`` can be either one or number of rows of ``region``.
#      If there are more than one numeric column, the length of ``col`` can be either one or number of numeric columns.
#      Pass to `circos.points`.
# -pch Type of points. Settings are similar as ``col``. Pass to `circos.points`.
# -cex Size of points. Settings are similar as ``col``. Pass to `circos.points`.
# -bg Background colors for points.
# -... Mysterious parameters.
#
# == details
# The function is a low-level graphical function and usually is put in ``panel.fun`` when using `circos.genomicTrack`.
#
# The function behaviours differently from different formats of input, see the examples in 
# the "Examples" Section or go to https://jokergoo.github.io/circlize_book/book/modes-of-input.html for more details.
#
# == example
# circos.par("track.height" = 0.1)
# circos.initializeWithIdeogram(plotType = NULL)
#
# bed = generateRandomBed(nr = 100)
# circos.genomicTrack(bed, panel.fun = function(region, value, ...) {
#     circos.genomicPoints(region, value, pch = 16, cex = 0.5, ...)
# })
#
# circos.genomicTrack(bed, stack = TRUE, panel.fun = function(region, value, ...) {
#     circos.genomicPoints(region, value, pch = 16, cex = 0.5, ...)
#     i = getI(...)
#     cell.xlim = get.cell.meta.data("cell.xlim")
#     circos.lines(cell.xlim, c(i, i), lty = 2, col = "#00000040")
# })
#
# bed1 = generateRandomBed(nr = 100)
# bed2 = generateRandomBed(nr = 100)
# bed_list = list(bed1, bed2)
#
# # data frame list
# circos.genomicTrack(bed_list, panel.fun = function(region, value, ...) {
#     cex = (value[[1]] - min(value[[1]]))/(max(value[[1]]) - min(value[[1]]))
#     i = getI(...)
#     circos.genomicPoints(region, value, cex = cex, pch = 16, col = i, ...)
# })
#
# circos.genomicTrack(bed_list, stack = TRUE,
#     panel.fun = function(region, value, ...) {
#     cex = (value[[1]] - min(value[[1]]))/(max(value[[1]]) - min(value[[1]]))
#     i = getI(...)
#     circos.genomicPoints(region, value, cex = cex, pch = 16, col = i, ...)
#     cell.xlim = get.cell.meta.data("cell.xlim")
#     circos.lines(cell.xlim, c(i, i), lty = 2, col = "#00000040")
# })
#
# bed = generateRandomBed(nr = 100, nc = 4)
# circos.genomicTrack(bed, panel.fun = function(region, value, ...) {
#     cex = (value[[1]] - min(value[[1]]))/(max(value[[1]]) - min(value[[1]]))
#     circos.genomicPoints(region, value, cex = 0.5, pch = 16, col = 1:4, ...)
# })
#
# circos.genomicTrack(bed, stack = TRUE, panel.fun = function(region, value, ...) {
#     cex = (value[[1]] - min(value[[1]]))/(max(value[[1]]) - min(value[[1]]))
#     i = getI(...)
#     circos.genomicPoints(region, value, cex = cex, pch = 16, col = i, ...)
#     cell.xlim = get.cell.meta.data("cell.xlim")
#     circos.lines(cell.xlim, c(i, i), lty = 2, col = "#00000040")
# })
#
# circos.clear()
#
circos.genomicPoints = function(
	region, 
	value, 
	numeric.column = NULL, 
	sector.index = get.cell.meta.data("sector.index"),
    track.index = get.cell.meta.data("track.index"), 
    posTransform = NULL, 
	pch = par("pch"), 
	col = par("col"), 
	cex = par("cex"), 
	bg = par("bg"), 
	...) {
	
	nr = nrow(region)
	if(ncol(region) > 2 && inherits(region[, 1], c("character", "factor"))) {
		region = region[, -1, drop = FALSE]
	}
	
	if(is.atomic(value) && length(value) == 1) {
		value = data.frame(value = rep(value, nr))
	}
	if(is.atomic(value) && length(value) == nr) {
		value = data.frame(value = value)
	}
	if(!is.data.frame(value)) stop_wrap("`value` should be a data frame.")
	
	args = list(...)
	if(!is.null(args$.param)) {
		.param = args$.param
		if(!is.null(.param$stack)) {
			if(.param$stack && is.null(numeric.column)) {
				if(is.null(.param$jitter)) {
					value = data.frame(hline = rep(.param$i, nr))
				} else {
					value = data.frame(hline = rep(.param$i, nr) + (runif(nr) - 0.5)*abs(.param$jitter))
				}
				numeric.column = 1
			}
		} else if(!is.null(.param$numeric.column) && is.null(numeric.column)) {
			numeric.column = .param$numeric.column
		}
	}
	
	if(is.vector(value) && !is.list(value) && length(value) == 1) {
		value = data.frame(value = rep(value, nr))
		numeric.column = 1
	} else if(is.vector(value) && !is.list(value) && length(value) == nr) {
		value = data.frame(value = value)
		numeric.column = 1
	}

	if(ncol(value) == 1) numeric.column = 1
	
	if(!is.null(posTransform)) {
		region = posTransform(region)
	}
	
	if(is.null(numeric.column)) {
		numeric.column = which(as.logical(sapply(value, is.numeric)))
		if(length(numeric.column) == 0) {
			stop_wrap("Cannot find numeric column.")
		}
	}
	
	nc = length(numeric.column)
	pch = .normalizeGraphicalParam(pch, nc, nr, "pch")
	col = .normalizeGraphicalParam(col, nc, nr, "col")
	cex = .normalizeGraphicalParam(cex, nc, nr, "cex")
	bg = .normalizeGraphicalParam(bg, nc, nr, "cex")
	
	if(nc == 1) {
		circos.points( (region[[1]] + region[[2]])/2, value[[ numeric.column ]], 
			pch = pch, col = col, cex = cex, bg = bg,
			sector.index = sector.index, track.index = track.index )
	} else {
		for(i in seq_len(nc)) {
			circos.points( (region[[1]] + region[[2]])/2, value[[ numeric.column[i] ]], 
				pch = pch[i], col = col[i], cex = cex[i], bg = bg[i],
				sector.index = sector.index, track.index = track.index ) 
		}
	}
}

# == title
# Add lines to a plotting region, specifically for genomic graphics
#
# == param
# -region A data frame contains 2 column which correspond to start positions and end positions.
# -value  A data frame contains values and other information.
# -numeric.column Which column in ``value`` data frame should be taken as y-value.
#                 If it is not defined, the whole numeric columns in ``value`` will be taken.
# -sector.index Index of sector.
# -track.index Index of track.
# -posTransform Self-defined function to transform genomic positions, see `posTransform.default` for explaination.
# -col col of lines/areas. If there are more than one numeric column, the length of ``col`` can be either one or number of numeric columns.
#      If there is only one numeric column and type is either ``segment`` or ``h``, 
#      the length of ``col`` can be either one or number of rows of ``region``.
#      pass to `circos.lines`
# -lwd Settings are similar as ``col``. Pass to `circos.lines`.
# -lty Settings are similar as ``col``. Pass to `circos.lines`.
# -type There is an additional option ``segment`` which plot segment lines from start position to end position. Settings are similar as ``col``. Pass to `circos.lines`. 
# -area Settings are similar as ``col``. Pass to `circos.lines`.
# -area.baseline Deprecated, use ``baseline`` instead.
# -baseline Settings are similar as ``col``. Pass to `circos.lines`.
# -border Settings are similar as ``col``. Pass to `circos.lines`.
# -pt.col Settings are similar as ``col``. Pass to `circos.lines`.
# -cex Settings are similar as ``col``. Pass to `circos.lines`.
# -pch Settings are similar as ``col``. Pass to `circos.lines`.
# -... Mysterious parameters.
#
# == details
# The function is a low-level graphical function and usually is put in ``panel.fun`` when using `circos.genomicTrack`.
#
# The function behaviours differently from different formats of input, see the examples in 
# the "Examples" Section or go to https://jokergoo.github.io/circlize_book/book/modes-of-input.html for more details.
#
# == examples
# \donttest{
# ### test bed
# circos.par("track.height" = 0.1)
# circos.initializeWithIdeogram(plotType = NULL)
#
# bed = generateRandomBed(nr = 100)
# circos.genomicTrack(bed, panel.fun = function(region, value, ...) {
#     circos.genomicLines(region, value, type = "l", ...)
# })
#
# bed1 = generateRandomBed(nr = 100)
# bed2 = generateRandomBed(nr = 100)
# bed_list = list(bed1, bed2)
#
# circos.genomicTrack(bed_list, panel.fun = function(region, value, ...) {
#     i = getI(...)
#     circos.genomicLines(region, value, col = i, ...)
# })
#
# circos.genomicTrack(bed_list, stack = TRUE, 
#     panel.fun = function(region, value, ...) {
#     i = getI(...)
#     circos.genomicLines(region, value, col = i, ...)
# })
#
# bed = generateRandomBed(nr = 100, nc = 4)
# circos.genomicTrack(bed, panel.fun = function(region, value, ...) {
#     circos.genomicLines(region, value, col = 1:4, ...)
# })
#
# circos.genomicTrack(bed, stack = TRUE, panel.fun = function(region, value, ...) {
#     i = getI(...)
#     circos.genomicLines(region, value, col = i, ...)
# })
#
# bed = generateRandomBed(nr = 100)
# circos.genomicTrack(bed, panel.fun = function(region, value, ...) {
#     circos.genomicLines(region, value, type = "segment", lwd = 2, ...)
# })
#
# circos.clear()
# }
circos.genomicLines = function(
	region, 
	value, 
	numeric.column = NULL, 
	sector.index = get.current.sector.index(),
    track.index = get.current.track.index(), 
    posTransform = NULL, 
	col = ifelse(area, "grey", "black"), 
	lwd = par("lwd"),
    lty = par("lty"), 
    type = "l",
    area = FALSE, 
    area.baseline = NULL, 
    border = "black", 
    baseline = "bottom",
    pt.col = par("col"), 
    cex = par("cex"), 
    pch = par("pch"), 
    ...) {
	
	if(!is.null(area.baseline)) {
		baseline = area.baseline
		warning_wrap("`area.baseline` is deprecated, please use `baseline` instead.")
	}
	nr = nrow(region)
	if(ncol(region) > 2 && inherits(region[, 1], c("character", "factor"))) {
		region = region[, -1, drop = FALSE]
	}

	if(is.atomic(value) && length(value) == 1) {
		value = data.frame(value = rep(value, nr))
	}
	if(is.atomic(value) && length(value) == nr) {
		value = data.frame(value = value)
	}
	if(!is.data.frame(value)) stop_wrap("`value` should be a data frame.")
	
	args = list(...)
	if(!is.null(args$.param)) {
		.param = args$.param
		if(!is.null(.param$stack)) {
			if(.param$stack && is.null(numeric.column)) {
				value = data.frame(hline = rep(.param$i, nr))
				numeric.column = 1
				type = rep("segment", length(type))
			}
		} else if(!is.null(.param$numeric.column) && is.null(numeric.column)) {
			numeric.column = .param$numeric.column
		}
	}
	
	if(is.vector(value) && !is.list(value) && length(value) == 1) {
		value = data.frame(value = rep(value, nr))
		numeric.column = 1
	}
	if(is.vector(value) && !is.list(value) && length(value) == nr) {
		value = data.frame(value = value)
		numeric.column = 1
	}
	if(ncol(value) == 1) numeric.column = 1

	if(!is.null(posTransform)) {
		region = posTransform(region)
	}

	if(is.null(numeric.column)) {
		numeric.column = which(as.logical(sapply(value, is.numeric)))
		if(length(numeric.column) == 0) {
			stop_wrap("Cannot find numeric column.")
		}
	}

	nc = length(numeric.column)
	
	if(all(type %in% c("h", "segment"))) {
		col = .normalizeGraphicalParam(col, nc, nr, "col")
		lwd = .normalizeGraphicalParam(lwd, nc, nr, "col")
		lty = .normalizeGraphicalParam(lty, nc, nr, "col")
	} else {
		col = .normalizeGraphicalParam(col, nc, 1, "col")
		lwd = .normalizeGraphicalParam(lwd, nc, 1, "col")
		lty = .normalizeGraphicalParam(lty, nc, 1, "col")
	}
	pt.col = .normalizeGraphicalParam(pt.col, nc, 1, "col")
	cex = .normalizeGraphicalParam(cex, nc, 1, "col")
	pch = .normalizeGraphicalParam(pch, nc, 1, "col")
	type = .normalizeGraphicalParam(type, nc, 1, "col")
	area = .normalizeGraphicalParam(area, nc, 1, "col")
	baseline = .normalizeGraphicalParam(baseline, nc, 1, "col")
	border = .normalizeGraphicalParam(border, nc, 1, "col")

	if(!is.null(args$hline)) {
		# for(i in seq_len(nr)) {
		# 	circos.lines( c(region[i, 1], region[i, 2]), c(value[i, numeric.column], value[i, numeric.column]), 
		# 		col = col, lwd = lwd, lty = lty, type = "l",
		# 		sector.index = sector.index, track.index = track.index )
		# }
		circos.segments( region[, 1], value[, numeric.column], region[, 2], value[, numeric.column], 
				col = col, lwd = lwd, lty = lty, 
				sector.index = sector.index, track.index = track.index )
	} else if(nc == 1) {
		if(type == "segment") {
			# for(i in seq_len(nr)) {
			# 	circos.lines( c(region[i, 1], region[i, 2]), c(value[i, numeric.column], value[i, numeric.column]), 
			# 		col = col[i], lwd = lwd[i], lty = lty[i], type = "l",
			# 		sector.index = sector.index, track.index = track.index )
			# }
			circos.segments( region[, 1], value[, numeric.column], region[, 2], value[, numeric.column], 
					col = col, lwd = lwd, lty = lty, 
					sector.index = sector.index, track.index = track.index )
		} else {
			circos.lines( (region[[1]] + region[[2]])/2, value[[ numeric.column ]], 
				col = col, lwd = lwd, lty = lty, type = type, 
				area = area, baseline = baseline, 
				border = border, pt.col = pt.col, cex = cex, pch = pch,
				sector.index = sector.index, track.index = track.index )
		}
	} else {
		for(i in seq_len(nc)) {
			if(type[i] == "segment") {
				# for(k in seq_len(nr)) {
				# 	circos.lines( c(region[k, 1], region[k, 2]), c(value[k, numeric.column[i] ], value[k, numeric.column[i] ]), 
				# 		col = col[i], lwd = lwd[i], lty = lty[i], type = "l",
				# 		sector.index = sector.index, track.index = track.index )
				# }
				circos.segments( region[, 1], value[, numeric.column[i] ], region[, 2], value[, numeric.column[i] ], 
						col = col[i], lwd = lwd[i], lty = lty[i], type = "l",
						sector.index = sector.index, track.index = track.index )
			} else {
				circos.lines( (region[[1]] + region[[2]])/2, value[[ numeric.column[i] ]], 
					col = col[i], lwd = lwd[i], lty = lty[i], type = type[i], 
					area = area[i], baseline = baseline[i], 
					border = border[i], pt.col = pt.col[i], cex = cex[i], pch = pch[i],
					sector.index = sector.index, track.index = track.index )
			}
		}
	}
}	

# == title
# Draw rectangle-like grid, specifically for genomic graphics
#
# == param
# -region A data frame contains 2 column which correspond to start positions and end positions.
# -value  A data frame contains values and other information.
# -ytop A vector or a single value indicating top position of rectangles.
# -ybottom A vector or a single value indicating bottom position of rectangles.
# -ytop.column If ``ytop`` is in ``value``, the index of the column.
# -ybottom.column If ``ybottom`` is in ``value``, the index of the column.
# -sector.index Index of sector.
# -track.index Index of track.
# -posTransform Self-defined function to transform genomic positions, see `posTransform.default` for explaination.
# -col The length of ``col`` can be either one or number of rows of ``region``. Pass to `circos.rect`.
# -border Settings are similar as ``col``. Pass to `circos.rect`.
# -lty Settings are similar as ``col``. Pass to `circos.rect`.
# -... Mysterious parameters.
#
# == details
# The function is a low-level graphical function and usually is put in ``panel.fun`` when using `circos.genomicTrack`.
#
# The function behaviours differently from different formats of input, see the examples in 
# the "Examples" Section or go to https://jokergoo.github.io/circlize_book/book/modes-of-input.html for more details.
#
# == example
# \donttest{
# circos.par("track.height" = 0.1, cell.padding = c(0, 0, 0, 0))
# circos.initializeWithIdeogram(plotType = NULL)
#
# bed1 = generateRandomBed(nr = 100)
# bed2 = generateRandomBed(nr = 100)
# bed_list = list(bed1, bed2)
# f = colorRamp2(breaks = c(-1, 0, 1), colors = c("green", "black", "red"))
# circos.genomicTrack(bed_list, stack = TRUE,
#     panel.fun = function(region, value, ...) {
#  
#     circos.genomicRect(region, value, col = f(value[[1]]), 
#         border = NA, ...)
#     i = getI(...)
#     cell.xlim = get.cell.meta.data("cell.xlim")
#     circos.lines(cell.xlim, c(i, i), lty = 2, col = "#000000")
# })
#
# circos.genomicTrack(bed_list, ylim = c(0, 3),
#     panel.fun = function(region, value, ...) {
#     i = getI(...)
#     circos.genomicRect(region, value, ytop = i+0.4, ybottom = i-0.4, col = f(value[[1]]), 
#         border = NA, ...)
#   
#     cell.xlim = get.cell.meta.data("cell.xlim")
#     circos.lines(cell.xlim, c(i, i), lty = 2, col = "#000000")
# })
#
# circos.genomicTrack(bed1, panel.fun = function(region, value, ...) {
#     circos.genomicRect(region, value, col = "red", border = NA, ...)
#
# })
#
# circos.genomicTrack(bed_list, panel.fun = function(region, value, ...) {
#     i = getI(...)
#     circos.genomicRect(region, value, col = i, border = NA, ...)
#
# })
#
# circos.clear()
# }
circos.genomicRect = function(
	region, 
	value = NULL, 
	ytop = NULL, 
	ybottom = NULL, 
	ytop.column = NULL, 
	ybottom.column = NULL,
	sector.index = get.current.sector.index(),
    track.index = get.current.track.index(), 
    posTransform = NULL, 
	col = NA, 
	border = "black", 
	lty = par("lty"), 
	...) {
	
	if(ncol(region) > 2 && inherits(region[, 1], c("character", "factor"))) {
		region = region[, -1, drop = FALSE]
	}
	nr = nrow(region)
	
	args = list(...)
	if(!is.null(args$.param)) {
		.param = args$.param
		if(!is.null(.param$stack)) {
			if(.param$stack) {
				if(is.null(ytop)) ytop = .param$i + 0.5
				if(is.null(ybottom)) ybottom = .param$i - 0.5
			}
		}
	}
	
	if(is.atomic(value) && length(value) == 1) {
		value = data.frame(value = rep(value, nr))
	}
	if(is.atomic(value) && length(value) == nr) {
		value = data.frame(value = value)
	}
	
	# 1. check ytop and ybottom
	# 2. check ytop.colum and ybottom.column
	
	if(!is.null(ytop)) {
		if(length(ytop) == 1) {
			ytop = rep(ytop, nr)
		}
		value = cbind(value, ytop)
		ytop.column = ncol(value)
	}
	
	if(!is.null(ybottom)) {
		if(length(ybottom) == 1) {
			ybottom = rep(ybottom, nr)
		}
		value = cbind(value, ybottom)
		ybottom.column = ncol(value)
	}
	if(is.matrix(value)) value = as.data.frame(value)
	if(!is.data.frame(value)) stop_wrap("`value` should be a data frame.")
	
	ylim = get.cell.meta.data("ylim", sector.index = sector.index, track.index = track.index)
	if(is.null(ybottom.column) && is.null(ytop.column)) {
		# if no ybottom and ytop column are set, the rect will draw along whole ylim
		value = cbind(value, rep(ylim[2], nr), rep(ylim[1], nr))
		ytop.column = ncol(value) - 1
		ybottom.column = ncol(value)
	} else if(is.null(ybottom.column)) {
		value = cbind(value, rep(ylim[1], nr))
		ybottom.column = ncol(value)
	} else if(is.null(ytop.column)) {
		value = cbind(value, rep(ylim[2], nr))
		ytop.column = ncol(value)
	}

	if(!is.null(posTransform)) {
		region = posTransform(region)
	}
	
	if(length(ytop.column) > 1) {
		stop_wrap("Only one ytop columns is allowed.")
	}
	if(length(ybottom.column) > 1) {
		stop_wrap("Only one ybottom columns is allowed.")
	}

	col = .normalizeGraphicalParam(col, 1, nr, "col")
	border = .normalizeGraphicalParam(border, 1, nr, "border")
	lty = .normalizeGraphicalParam(lty, 1, nr, "lty")
	
	# for(i in seq_len(nr)) {
	# 	circos.rect(region[i, 1], value[i, ybottom.column], region[i, 2], value[i, ytop.column],
	# 	            sector.index = sector.index, track.index = track.index,
	# 				col = col[i], border = border[i], lwd = lwd[i], lty = lty[i])
	# }
	circos.rect(region[, 1], value[, ybottom.column], region[, 2], value[, ytop.column],
		            sector.index = sector.index, track.index = track.index,
					col = col, border = border, lty = lty)

}

# == title
# Draw text in a cell, specifically for genomic graphics
#
# == param
# -region A data frame contains 2 column which correspond to start positions and end positions.
# -value  A data frame contains values and other information.
# -y A vector or a single value indicating position of text.
# -labels Labels of text corresponding to each genomic positions.
# -labels.column If labels are in ``value``, index of column in ``value``.
# -numeric.column Which column in ``value`` data frame should be taken as y-value.
#                 If it is not defined, only the first numeric columns in ``value`` will be taken.
# -sector.index Index of sector.
# -track.index Index of track.
# -posTransform Self-defined function to transform genomic positions, see `posTransform.default` for explanation.
# -facing Passing to `circos.text`. Settings are similar as ``col``.
# -niceFacing   Should the facing of text be adjusted to fit human eyes?
# -direction Deprecated, use ``facing`` instead. 
# -adj Pass to `circos.text`. Settings are similar as ``col``.
# -cex Pass to `circos.text`. Settings are similar as ``col``.
# -col Pass to `circos.text`. The length of ``col`` can be either one or number of rows of ``region``.
# -font Pass to `circos.text`. Settings are similar as ``col``.
# -padding pass to ``posTransform`` if it is set as `posTransform.text`.
# -extend pass to ``posTransform`` if it is set as `posTransform.text`.
# -align_to pass to ``posTransform`` if it is set as `posTransform.text`.
# -... Mysterious parameters.
#
# == details
# The function is a low-level graphical function and usually is put in ``panel.fun`` when using `circos.genomicTrack`.
#
# == example
# circos.par("track.height" = 0.1, cell.padding = c(0, 0, 0, 0))
# circos.initializeWithIdeogram(plotType = NULL)
#
# bed = generateRandomBed(nr = 20)
#
# circos.genomicTrack(bed, ylim = c(0, 1), panel.fun = function(region, value, ...) {
#     circos.genomicText(region, value, y = 0.5, labels = "text", ...)
# })
#
# bed = cbind(bed, sample(letters, nrow(bed), replace = TRUE))
# circos.genomicTrack(bed, panel.fun = function(region, value, ...) {
#     circos.genomicText(region, value, labels.column = 2, ...)
# })
#
# circos.clear()
circos.genomicText = function(
	region, 
	value = NULL, 
	y = NULL, 
	labels = NULL, 
	labels.column = NULL,
	numeric.column = NULL, 
	sector.index = get.current.sector.index(), 
	track.index = get.current.track.index(), 
	posTransform = NULL, 
	direction = NULL, 
	facing = "inside", 
	niceFacing = FALSE,
	adj = par("adj"), 
	cex = 1, 
	col = "black", 
	font = par("font"), 
	padding = 0,
	extend = 0, 
	align_to = "region", 
	...) {
	
	if(!is.null(direction)) {
		facing = direction
		warning_wrap("`direction` is deprecated, please use `facing` instead.")
	}
	
	if(ncol(region) > 2 && inherits(region[, 1], c("character", "factor"))) {
		region = region[, -1, drop = FALSE]
	}
	nr = nrow(region)
	
	if(is.vector(value) && !is.list(value) && length(value) == 1) {
		value = data.frame(value = rep(value, nr))
		numeric.column = 1
	}
	if(is.vector(value) && !is.list(value) && length(value) == nr) {
		value = data.frame(value = value)
		numeric.column = 1
	}
	
	args = list(...)
	if(!is.null(args$.param)) {
		.param = args$.param
		if(!is.null(.param$stack)) {
			if(.param$stack && is.null(numeric.column)) {
				value = data.frame(hline = rep(.param$i, nr))
				numeric.column = 1
			}
		} else if(!is.null(.param$numeric.column) && is.null(numeric.column)) {
			numeric.column = .param$numeric.column
		}
	}
	
	if(!is.null(y)) {
		if(length(y) == 1) {
			y = rep(y, nr)
		}
		value = cbind(value, y)
		numeric.column = ncol(value)
	}
	if(is.matrix(value)) value = as.data.frame(value)
	if(!is.data.frame(value)) stop_wrap("`value` should be a data frame.")

	if(is.null(labels) && is.null(labels.column)) {
		stop_wrap("You should either specify `labels` or `labels.column`.")
	}
	
	if(!is.null(labels)) {
		if(is.vector(labels) && !is.list(labels) && length(labels) == 1) {
			value = cbind(value, labels = rep(labels, nr))
			labels.column = ncol(value)
		}
		if(is.vector(labels) && !is.list(labels) && length(labels) == nr) {
			value = cbind(value, labels = labels)
			labels.column = ncol(value)
		}
	}

	if(is.null(numeric.column)) {
		numeric.column = which(as.logical(sapply(value, is.numeric)))
		if(length(numeric.column) == 0) {
			stop_wrap("Cannot find numeric column.")
		}
		numeric.column = numeric.column[1]
	}
	
	if(length(numeric.column) > 1) {
		stop_wrap("You can only have one numeric column.")
	}
	
	if(!is.null(posTransform)) {
	
		# check settings when it is text-specific transformation
		if(identical(posTransform, posTransform.text)) {
			if(! facing %in% c("clockwise", "reverse.clockwise")) {
				stop_wrap("Only support `facing` in c('clockwise', 'reverse.clockwise') if `posTransform` is `posTransform.text`.")
			}
			region = posTransform(region, value[[ numeric.column ]], value[[labels.column]], cex, font, padding = padding, extend = extend, align_to = align_to)
		} else {
			region = posTransform(region)
		}
	}
	
	nc = length(numeric.column)

	col = .normalizeGraphicalParam(col, nc, nr, "col")
	cex = .normalizeGraphicalParam(cex, nc, nr, "cex")
	font = .normalizeGraphicalParam(font, nc, nr, "font")

	circos.text( (region[[1]] + region[[2]])/2, value[[ numeric.column ]], value[[labels.column]],
		facing = facing, niceFacing = niceFacing, adj = adj, cex = cex, col = col, font = font,
		sector.index = sector.index, track.index = track.index)

}


# == title
# Add links from two sets of genomic positions
#
# == param
# -region1 A data frame in bed format.
# -region2 A data frame in bed format.
# -rou Pass to `circos.link`.
# -rou1 Pass to `circos.link`.
# -rou2 Pass to `circos.link`.
# -col Pass to `circos.link`, length can be either one or nrow of ``region1``.
# -lwd Pass to `circos.link`, length can be either one or nrow of ``region1``.
# -lty Pass to `circos.link`, length can be either one or nrow of ``region1``.
# -border Pass to `circos.link`, length can be either one or nrow of ``region1``.
# -... Pass to `circos.link`.
#
# == details
# Of course, number of rows should be same in ``region1`` and ``region2``.
#
# If you want to have more controls on links, please use `circos.link` directly.
#
# == seealso
# https://jokergoo.github.io/circlize_book/book/genomic-plotting-region.html#genomic-links
#
# == example
# \donttest{
# set.seed(123)
#
# bed1 = generateRandomBed(nr = 100)
# bed1 = bed1[sample(nrow(bed1), 20), ]
# bed2 = generateRandomBed(nr = 100)
# bed2 = bed2[sample(nrow(bed2), 20), ]
# circos.par("track.height" = 0.1, cell.padding = c(0, 0, 0, 0))
# circos.initializeWithIdeogram()
#
# circos.genomicLink(bed1, bed2, col = sample(1:5, 20, replace = TRUE), border = NA)
# circos.clear()
# }
circos.genomicLink = function(
	region1, 
	region2, 
	rou = get_most_inside_radius(), 
	rou1 = rou, 
	rou2 = rou,
    col = "black", 
    lwd = par("lwd"), 
    lty = par("lty"), 
    border = col, 
    ...) {

	if(circos.par$ring) {
		l = region1[, 2] > region1[, 3]
		if(any(l)) {
			region1[l, 2] = region1[l, 2] - diff(get.sector.data()[c("min.value", "max.value")])
		}
		l = region2[, 2] > region2[, 3]
		if(any(l)) {
			region2[l, 2] = region2[l, 2] - diff(get.sector.data()[c("min.value", "max.value")])
		}
	}
	
	region1 = validate_data_frame(region1)
	region2 = validate_data_frame(region2)

	if(ncol(region1) == 2) {
		region1[, 3] = region1[, 2]
	}
	if(ncol(region2) == 2) {
		region2[, 3] = region2[, 2]
	}

	region1 = normalizeToDataFrame(region1, sort = FALSE)
	region2 = normalizeToDataFrame(region2, sort = FALSE)
	
	if(is.dataFrameList(region1)) {
		stop_wrap("`region1` can not be a region list.")
	}
	
	if(is.dataFrameList(region1)) {
		stop_wrap("`region1` can not be a region list.")
	}
	
	if(nrow(region1) != nrow(region2)) {
		stop_wrap("nrow of `region1` and `region2` differ. Please check the chromosome column and make sure all the chromosomes are in the circular layout.")
	}
	
	nr = nrow(region1)
	
	rou1 = .normalizeGraphicalParam(rou1, 1, nr, "rou")
	rou2 = .normalizeGraphicalParam(rou2, 1, nr, "rou")
	col = .normalizeGraphicalParam(col, 1, nr, "col")
	lwd = .normalizeGraphicalParam(lwd, 1, nr, "lwd")
	lty = .normalizeGraphicalParam(lty, 1, nr, "lty")
	border = .normalizeGraphicalParam(border, 1, nr, "border")

	for(i in seq_len(nr)) {
		if(region1[i, 2] == region1[i, 3]) {
			point1 = region1[i, 2]
		} else {
			point1 = c(region1[i, 2], region1[i, 3])
		}
		if(region2[i, 2] == region2[i, 3]) {
			point2 = region2[i, 2]
		} else {
			point2 = c(region2[i, 2], region2[i, 3])
		}
		circos.link(region1[i, 1], point1,
		            region2[i, 1], point2,
					rou1 = rou1[i], rou2 = rou2[i], col = col[i], lwd = lwd[i],
					lty = lty[i], border = border[i], ...)
	}
}


# == title
# Add genomic position transformation lines between tracks
#
# == param
# -data A data frame containing genomic data.
# -track.height Height of the track.
# -posTransform Genomic position transformation function, see `posTransform.default` for an example.
# -horizontalLine Whether to draw horizontal lines which indicate region width .
# -track.margin Margin of tracks.
# -direction Type of the transformation. ``inside`` means position transformed track are located inside 
#       and ``outside`` means position transformed track are located outside.
# -col Color of lines, can be length of one or nrow of ``data``.
# -lwd Width of lines.
# -lty Style of lines.
# -... Pass to `circos.trackPlotRegion`.
#
# == details
# There is one representative situation when such position transformation needs to be applied. 
# For example, there are two sets of regions in a chromosome in which regions in one set regions are
# quite densely to each other and regions in other set are far from others. Heatmap or text is going
# to be drawn on the next track. If there is no position transformation, heatmap or text for those
# dense regions would be overlapped and hard to identify, also ugly to visualize. Thus, a way
# to transform original positions to new positions would help for the visualization. 
circos.genomicPosTransformLines = function(
	data, 
	track.height = 0.1, 
	posTransform = NULL, 
	horizontalLine = c("none", "top", "bottom", "both"), 
	track.margin = c(0, 0),
	direction = c("inside", "outside"), 
	col = "black", 
	lwd = par("lwd"),
    lty = par("lty"), 
    ...) {
	
	horizontalLine = match.arg(horizontalLine)[1]

	data = validate_data_frame(data)
	data = normalizeToDataFrame(data)
	
	if(is.dataFrameList(data)) {
		stop_wrap("`data` can not be list of regions.")
	}
	
	nr = nrow(data)
	if(length(col) == 1) {
		col = rep(col, nr)
	}
	if(length(lwd) == 1) {
		lwd = rep(lwd, nr)
	}
	if(length(lty) == 1) {
		lty = rep(lty, nr)
	}
	
	o.track.margin = circos.par("track.margin")
	circos.par(track.margin = track.margin)
	
	if(direction[1] == "default") direction = "outside"
	if(direction[1] == "reverse") direction = "inside"
	direction = match.arg(direction)[1]
	
	if(direction == "inside") {
		circos.genomicTrackPlotRegion(data, ylim = c(0, 1), bg.border = NA, track.height = track.height, panel.fun = function(region, value, ...) {
			chr = get.current.chromosome()
			l = data[[1]] == chr
			if(!is.null(posTransform)) {
				if(is.function(posTransform)) {
					args = as.list(posTransform)
					if(length(args) == 2) {
						region_new = posTransform(region)
					} else if(length(args) == 3) {
						region_new = posTransform(region, value)
					}
				}
			} else {
				region_new  = region
			}
			
			# for(i in seq_len(nrow(region))) {
			# 	if(horizontalLine == "both" || horizontalLine == "top") {
			# 		circos.lines(c(region[i, 1], region[i, 2]), c(1, 1), col = col[l][i], lwd = lwd[l][i], lty = lty[l][i])
			# 	}
			# 	if(horizontalLine == "both" || horizontalLine == "bottom") {
			# 		circos.lines(c(region[i, 1], region[i, 2]), c(0, 0), col = col[l][i], lwd = lwd[l][i], lty = lty[l][i])
			# 	}
			# 	mid = (region[i, 1] + region[i, 2])/2
			# 	mid_new = (region_new[i, 1] + region_new[i, 2])/2
			# 	circos.lines(c(mid, mid, mid_new, mid_new), c(1, 2/3, 1/3, 0), col = col[l][i], lwd = lwd[l][i], lty = lty[l][i])
			# }
			nr = nrow(region)
			if(horizontalLine == "both" || horizontalLine == "top") {
				circos.segments(region[, 1], rep(1, nr), region[, 2], rep(1, nr), col = col[l], lwd = lwd[l], lty = lty[l])
			}
			if(horizontalLine == "both" || horizontalLine == "bottom") {
				circos.segments(region[, 1], rep(0, nr), region[, 2], rep(0, nr), col = col[l], lwd = lwd[l], lty = lty[l])
			}
			mid = (region[, 1] + region[, 2])/2
			mid_new = (region_new[, 1] + region_new[, 2])/2
			circos.segments(mid, rep(1, nr), mid, rep(2/3, nr), col = col[l], lwd = lwd[l], lty = lty[l])
			circos.segments(mid, rep(2/3, nr), mid_new, rep(1/3, nr), col = col[l], lwd = lwd[l], lty = lty[l])
			circos.segments(mid_new, rep(1/3, nr), mid_new, rep(0, nr), col = col[l], lwd = lwd[l], lty = lty[l])
		}, ...)
	} else {
		circos.genomicTrackPlotRegion(data, ylim = c(0, 1), bg.border = NA, track.height = track.height, panel.fun = function(region, value, ...) {
			chr = get.current.chromosome()
			l = data[[1]] == chr
			region_subset = data[l, , drop = FALSE]
			if(!is.null(posTransform)) {
				if(is.function(posTransform)) {
					args = as.list(posTransform)
					if(length(args) == 2) {
						region_new = posTransform(region)
					} else if(length(args) == 3) {
						region_new = posTransform(region, value)
					}
				}
			} else {
				region_new  = region
			}
			
			# for(i in seq_len(nrow(region))) {
			# 	if(horizontalLine == "both" || horizontalLine == "bottom") {
			# 		circos.lines(c(region[i, 1], region[i, 2]), c(0, 0), col = col[l][i], lwd = lwd[l][i], lty = lty[l][i])
			# 	}
			# 	if(horizontalLine == "both" || horizontalLine == "top") {
			# 		circos.lines(c(region[i, 1], region[i, 2]), c(1, 1), col = col[l][i], lwd = lwd[l][i], lty = lty[l][i])
			# 	}
			# 	mid = (region[i, 1] + region[i, 2])/2
			# 	mid_new = (region_new[i, 1] + region_new[i, 2])/2
			# 	circos.lines(c(mid, mid, mid_new, mid_new), c(0, 1/3, 2/3, 1), col = col[l][i], lwd = lwd[l][i], lty = lty[l][i])
			# }
			nr = nrow(region)
			if(horizontalLine == "both" || horizontalLine == "bottom") {
				circos.segments(region[, 1], rep(0, nr), region[, 2], rep(0, nr), col = col[l], lwd = lwd[l], lty = lty[l])
			}
			if(horizontalLine == "both" || horizontalLine == "top") {
				circos.segments(region[, 1], rep(1, nr), region[, 2], rep(1, nr), col = col[l], lwd = lwd[l], lty = lty[l])
			}
			mid = (region[, 1] + region[, 2])/2
			mid_new = (region_new[, 1] + region_new[, 2])/2
			circos.segments(mid, rep(0, nr), mid, rep(1/3, nr), col = col[l], lwd = lwd[l], lty = lty[l])
			circos.segments(mid, rep(1/3, nr), mid_new, rep(2/3, nr), col = col[l], lwd = lwd[l], lty = lty[l])
			circos.segments(mid_new, rep(2/3, nr), mid_new, rep(1, nr), col = col[l], lwd = lwd[l], lty = lty[l])
		}, ...)
	}
	
	circos.par(track.margin = o.track.margin)
}

# == title
# Calculate and add genomic density track
#
# == param
# -data A bed-file-like data frame or a list of data frames. If the input is a list of data frames.
#     there will be multiple density plot in one same track.
# -ylim.force Whether to force upper bound of ``ylim`` to be 1. Ignored if ``count_by`` is set to ``number``.
# -window.size Pass to `genomicDensity`.
# -overlap Pass to `genomicDensity`.
# -count_by Pass to `genomicDensity`.
# -col  Colors. It should be length of one. If ``data`` is a list of data frames, the length of ``col``
#       can also be the length of the list. If multiple sets of genomic regions are visualized in one
#       single track, you should set the colors with transparency to distinguish them.
# -lwd  Width of lines, the same setting as ``col`` argument.
# -lty  Style of lines, the same setting as ``col`` argument.
# -type Type of lines, see `circos.lines`.
# -area See `circos.lines`.
# -area.baseline Deprecated, use ``baseline`` instead.
# -baseline See `circos.lines`.
# -border See `circos.lines`.
# -... Pass to `circos.trackPlotRegion`.
#
# == details
# This function is a high-level graphical function, and it will create a new track.
#
# If you have multiple sets of genomic regions, you should make sure the density ranges 
# for all sets are similar, or I suggest you should put them into different tracks. One example
# can be found in the "Examples" Section where the density range for ``bed_list[[2]]`` is too high
# compared to the range for ``bed_list[[1]]``, thus, it is better to put the two sets of
# regions into two separate tracks.
#
# == seealso
# https://jokergoo.github.io/circlize_book/book/high-level-genomic-functions.html#genomic-density-and-rainfall-plot
#
# == example
# load(system.file(package = "circlize", "extdata", "DMR.RData"))
#
# # rainfall
# \donttest{
# circos.initializeWithIdeogram(plotType = c("axis", "labels"))
#
# bed_list = list(DMR_hyper, DMR_hypo)
# circos.genomicRainfall(bed_list, pch = 16, cex = 0.4, col = c("#FF000080", "#0000FF80"))
#
# circos.genomicDensity(bed_list[[1]], col = c("#FF000080"), track.height = 0.1)
# circos.genomicDensity(bed_list[[2]], col = c("#0000FF80"), track.height = 0.1)
# circos.clear()
#
# ############ draw the two densities in one track  #############
# circos.initializeWithIdeogram(plotType = c("axis", "labels"))
# circos.genomicDensity(bed_list, col = c("#FF000080", "#0000FF80"), track.height = 0.2)
# circos.clear()
# }
circos.genomicDensity = function(
	data, 
	ylim.force = FALSE, 
	window.size = NULL, 
	overlap = TRUE, 
	count_by = c("percent", "number"),
	col = ifelse(area, "grey", "black"), 
	lwd = par("lwd"), 
	lty = par("lty"), 
	type = "l",
	area = TRUE, 
	area.baseline = NULL, 
	baseline = 0, 
	border = NA, 
	...) {
	
	if(!is.null(area.baseline)) {
		baseline = area.baseline
		warning_wrap("`area.baseline` is deprecated, please use `baseline` instead.")
	}
	
	data = normalizeToDataFrame(data)
	
	if(!is.dataFrameList(data)) {
		data = list(data)
	}
	
	if(length(col) == 1) {
		col = rep(col, length(data))
	}
	if(length(lwd) == 1) {
		lwd = rep(lwd, length(data))
	}
	if(length(lty) == 1) {
		lty = rep(lty, length(data))
	}
	if(length(type) == 1) {
		type = rep(type, length(data))
	}
	if(length(area) == 1) {
		area = rep(area, length(data))
	}
	if(length(baseline) == 1) {
		baseline = rep(baseline, length(data))
	}
	if(length(border) == 1) {
		border = rep(border, length(data))
	}

	s = sapply(get.all.sector.index(), function(si) get.cell.meta.data("xrange", sector.index = si))
	if(is.null(window.size)) {
		window.size = 10^nchar(sum(s))/1000  # around 100 major ticks
		#cat(window.size, "is choosen as the window size.\n")
	}
	
	count_by = match.arg(count_by)[1]
	df = vector("list", length = length(data))
	for(i in seq_along(data)) {
		df[[i]] = genomicDensity(data[[i]], window.size = window.size, overlap = overlap, count_by = count_by)
	}
	if(ylim.force && count_by == "percent") {
		ymax = 1
	} else {
		ymax = max(sapply(df, function(gr) max(gr[[4]])))
	}
	if(length(df) == 1) {
		circos.genomicTrackPlotRegion(df[[1]], ylim = c(0, ymax), panel.fun = function(region, value, ...) {
			circos.genomicLines(region, value, col = col, lwd = lwd, lty = lty, type = type, 
				border = border, area = area, baseline = baseline) 
		}, ...)
	} else {
		circos.genomicTrackPlotRegion(df, ylim = c(0, ymax), panel.fun = function(region, value, ...) {
			i = getI(...)
			circos.genomicLines(region, value, col = col[i], lwd = lwd[i], lty = lty[i], type = type[i], 
				border = border[i], area = area[i], baseline = baseline[i]) 
		}, ...)
	}
}

# == title
# Calculate genomic region density
#
# == param
# -region Genomic positions. It can be a data frame with two
#     columns which are start positions and end positions on a single chromosome.
#     It can also be a bed-format data frame which contains the chromosome column.
# -window.size Window size to calculate genomic density
# -n.window number of windows, if it is specified, ``window.size`` is ignored
# -overlap Whether two neighbouring windows have half overlap
# -count_by How to count the value for each window, ``percent``: percent of the window covered by the input regions; ``number``: number of regions that overlap to the window.
# -chr.len the chromosome length. The value should be named vector
#
# == details
# It calculate the percent of each genomic windows that is covered by the input regions.
#
# == values
# If the input is a two-column data frame, the function returns a data frame with three columns: 
# start position, end position and the overlapping (value depends on the ``count_by`` argument). And if the input is a bed-format
# data frame, there will be an additionally chromosome name column.
#
# == example
# bed = generateRandomBed()
# bed = subset(bed, chr == "chr1")
# head(genomicDensity(bed))
# head(genomicDensity(bed, count_by = "number"))
genomicDensity = function(
	region, 
	window.size = 1e7, 
	n.window = NULL, 
	overlap = TRUE, 
	count_by = c("percent", "number"),
	chr.len = NULL) {
	
	region = validate_data_frame(region)

	if(is.character(region[, 1]) || is.factor(region[, 1])) {
		validate_region(region, check_chr = TRUE)
		region[, 1] = as.vector(region[, 1])
		all_chr = unique(region[, 1])
		if(!is.null(chr.len)) {
			if(length(all_chr) == 1 & length(chr.len) == 1) {
				names(chr.len) = all_chr[1]
			}
		}
		return(do.call("rbind", lapply(all_chr, function(chr) {
			l = region[,1] == chr
			if(is.null(chr.len)) {
				max_rg = NULL
				if(is.circos.initialized()) {
					if(chr %in% get.all.sector.index()) {
						max_rg = get.sector.data(sector.index = chr)["max.data"]
					} 
				}
			} else {
				if(chr %in% names(chr.len)) {
					max_rg = chr.len[chr]
				} else {
					max_rg = NULL
				}
			}
			df = genomicDensity(region[l, 2:3, drop = FALSE], window.size = window.size, overlap = overlap, chr.len = max_rg, count_by = count_by)
			cbind(chr = rep(chr, nrow(df)), df)
		})))
	}
	
	validate_region(region, 1, 2, check_chr = FALSE)
	if(ncol(region) >= 3) {
		if(is.numeric(region[, 1])) {
			if(max(region[, 1]) < 100) {
				return(do.call("rbind", lapply(unique(region[, 1]), function(chr) {
					l = region[, 1] == chr
					df = genomicDensity(region[l, 2:3, drop = FALSE], window.size = window.size, overlap = overlap, count_by = count_by)
					cbind(chr = rep(chr, nrow(df)), df)
				})))
			}
		}
	}

	region = region[, 1:2]
	
	region = sort_region(region)
	region = reduce_region(region)

	# make a segmentation
	if(!is.null(chr.len)) {
		max_pos = max(c(chr.len, max(region[[2]])))
	} else {
		max_pos = max(region[[2]])
	}
	if(overlap) {
		if(missing(n.window)) {
			b = seq(1, max_pos, by = window.size/2)
			s = b[-length(b)]
			s = s[-length(s)]
			e = s + window.size - 1
		} else {
			b = seq(1, max_pos, length.out = 2*n.window - 1)
			s = b[-length(b)]
			s = s[-length(s)]
			e = s + b[3] - b[1] - 1
		}
	} else {
		if(missing(n.window)) {
			b = seq(1, max_pos, by = window.size)
			s = b[-length(b)]
			e = s + window.size - 1
		} else {
			b = seq(1, max_pos, length.out = n.window)
			s = b[-length(b)]
			e = s + b[2] - b[1]	
		}
		
	}

	s = as.integer(s)
	e = as.integer(e)

	y = rep(0, length(s))
	names(y) = paste(s, e, sep = ",")
	
	windows = data.frame(start = s, end = e)
	
	op = overlap_region(windows, region, count_by = count_by)

	res = data.frame(start = s, end = e, value = op)
	return(res)
}

# == title
# Highlight chromosomes
#
# == param
# -... pass to `highlight.sector`
#
# == details
# This is only a shortcut function of `highlight.sector`.
#
highlight.chromosome = function(...) {
	
	highlight.sector(...)
}

continuousIndexSegment = function(x, n = NULL, loop = FALSE) {
	if(length(x) == 1) {
		return(list(x))
	} else {
		k = c(0, which(diff(x) > 1), length(x))
		lt = vector("list", length = length(k) - 1)
		for(i in seq_along(k)[-length(k)]) {
			lt[[i]] = x[(k[i] + 1):(k[i+1])]
		}
		
		if(loop && length(lt) > 1) {
			first = lt[[1]]
			last = lt[[length(lt)]]
			if(first[1] == 1 && last[length(last)] == n) {
				lt[[1]] = c(last, first)
				lt = lt[-length(lt)]
			}
		}
		
		return(lt)
	}
}

# == title
# Get current chromosome name
#
# == details
# The function is same as `get.current.sector.index` and
# should only be put inside ``panel.fun`` when using `circos.genomicTrackPlotRegion`.
get.current.chromosome = function() {
	get.current.sector.index()
}




is.dataFrameList = function(data) {
	is.list(data) && all(sapply(data, is.data.frame))
}


normalizeToDataFrame = function(data, sort = FALSE) {

	all.chr = get.all.sector.index()

	if(is.data.frame(data)) {
		data = as.data.frame(data)
		if(ncol(data) < 3) {
			stop_wrap("Your data frame is less than 3 column!.")
		}
		data = data[data[[1]] %in% all.chr, , drop = FALSE]
		if(sort) {
			data = data[order(data[[1]], data[[2]]), , drop = FALSE]
		}
		return(data)
	} else if(is.list(data) && all(sapply(data, is.data.frame))) {
		df = lapply(data, function(gr) {
			if(ncol(gr) < 3) {
				stop_wrap("Your data frame is less than 3 column!.")
			}
			gr = gr[gr[[1]] %in% all.chr, , drop = FALSE]
			if(sort) {
				gr = gr[order(gr[[1]], gr[[2]]), ]
			}
			return(gr)
		})
		return(df)
	} else if(inherits(df, "GRanges")) {
		df = as.data.frame(df)
		return(df)
	} else {
		stop_wrap("The format of `data` should only be a data frame or a list of data frames.")
	}
}


.normalizeGraphicalParam = function(x, nc, nr, name) {
	if(nc == 1) {
		if(!(length(x) == 1 || length(x) == nr)) {
			stop_wrap("The length of `", name, "` (", length(x), ") should be equal to 1 or the number of your regions (", nr, ").")
		} else if(length(x) == 1) {
			x = rep(x, nr)
		}
	} else {
		if(!(length(x) == 1 || length(x) == nc)) {
			stop_wrap("The length of `", name, "` (", length(x), ") should be equal to 1 or the number of your data column (", nc, ").")
		} else if(length(x) == 1) {
			x = rep(x, nc)
		}
	}
	return(x)
}

# == title
# Genomic rainfall plot
#
# == param
# -data A bed-file-like data frame or a list of data frames.
# -mode How to calculate the distance of two neighbouring regions, pass to `rainfallTransform`.
# -ylim ylim for rainfall plot track. If ``normalize_to_width`` is ``FALSE``, the value should correspond to ``log10(dist+1)``,
#       and if ``normalize_to_width`` is ``TRUE``, the value should correspond to ``log2(rel_dist)``.
# -col  Color of points. It should be length of one. If ``data`` is a list, the length of ``col``
#       can also be the length of the list.
# -pch  Style of points.
# -cex  Size of points.
# -normalize_to_width If it is ``TRUE``, the value is the relative distance divided by the width of the region.
# -... Pass to `circos.trackPlotRegion`.
#
# == details
# This is high-level graphical function, which mean, it will create a new track.
#
# Rainfall plot can be used to visualize distribution of regions. On the plot, y-axis
# corresponds to the distance to neighbour regions (log-based). So if there is a drop-down on
# the plot, it means there is a cluster of regions at that area.
#
# On the plot, y-axis are log10-transformed.
#
# == seealso
# https://jokergoo.github.io/circlize_book/book/high-level-genomic-functions.html#genomic-density-and-rainfall-plot
#
# == example
# \donttest{
# load(system.file(package = "circlize", "extdata", "DMR.RData"))
#
# # rainfall
# circos.initializeWithIdeogram(plotType = c("axis", "labels"))
#
# bed_list = list(DMR_hyper, DMR_hypo)
# circos.genomicRainfall(bed_list, pch = 16, cex = 0.4, col = c("#FF000080", "#0000FF80"))
#
# circos.genomicDensity(bed_list[[1]], col = c("#FF000080"), track.height = 0.1)
# circos.genomicDensity(bed_list[[2]], col = c("#0000FF80"), track.height = 0.1)
#
# circos.clear()
# }
circos.genomicRainfall = function(
	data, 
	mode = "min", 
	ylim = NULL, 
	col = "black", 
	pch = par("pch"), 
	cex = par("cex"), 
	normalize_to_width = FALSE, 
	...) {
	
	data = normalizeToDataFrame(data)
	
	if(!is.dataFrameList(data)) {
		data = list(data)
	}
	
	for(i in seq_along(data)) {
		validate_region(data[[i]], check_chr = TRUE)
	}

	if(length(col) == 1) {
		col = rep(col, length(data))
	}
	if(length(pch) == 1) {
		pch = rep(pch, length(data))
	}
	if(length(cex) == 1) {
		cex = rep(cex, length(data))
	}

	if(is.null(ylim)) {
		if(normalize_to_width) {
			ylim = c(-10, 10)
		} else {
			ylim = c(0, 9)
		}
 	}
	
	circos.genomicTrackPlotRegion(data, ylim = ylim, panel.fun = function(region, value, ...) {
		df = rainfallTransform(region, mode = mode, normalize_to_width = normalize_to_width)
		i = getI(...)
		if(normalize_to_width) {
			circos.genomicPoints(df[1:2], log10(df[3]), col = col[i], cex = cex[i], pch = pch[i])
		} else {
			circos.genomicPoints(df[1:2], log10(df[3]+1), col = col[i], cex = cex[i], pch = pch[i])
		}
	}, ...)
	
}

# == title
# Calculate inter-distance of genomic regions
#
# == param
# -region Genomic positions. It can be a data frame with two
#     columns which are start positions and end positions on a single chromosome.
#     It can also be a bed-format data frame which contains the chromosome column.
# -mode How to calculate inter-distance. For a region, there is a distance to the 
#       prevous region and also there is a distance to the next region. ``mode``
#       controls how to merge these two distances into one value.
# -normalize_to_width If it is ``TRUE``, the value is the relative distance divided by the width of the region.
#
# == values
# If the input is a two-column data frame, the function returnes a data frame with three columns: start position, end position and distance.
# And if the input is a bed-format data frame, there will be the chromosome column added.
#
# The row order of the returned data frame is as same as the input one.
#
# == example
# bed = generateRandomBed()
# bed = subset(bed, chr == "chr1")
# head(rainfallTransform(bed))
#
rainfallTransform = function(
	region, 
	mode = c("min", "max", "mean", "left", "right"),
	normalize_to_width = FALSE) {
	
	mode = match.arg(mode)[1]

	region = validate_data_frame(region)
	if(is.character(region[, 1]) || is.factor(region[, 1])) {
		region[, 1] = as.vector(region[, 1])
		region = region[1:3]
		region$dist = NA
		for(chr in unique(region[, 1])) {
			l = region[, 1] == chr
			df = rainfallTransform(region[l, 2:3, drop = FALSE], mode = mode)
			region$dist[l] = df$dist
		}
		return(region)
	}
	if(ncol(region) >= 3) {
		if(is.numeric(region[, 1])) {
			if(max(region[, 1]) < 100) {
				region = as.data.frame(region)
				region[[1]] = as.character(region[[1]])
				rainfallTransform(region, mode = mode)
			}
		}
	}
	
	region_order = order_region(region[1:2])
	region_bk = as.data.frame(region[1:2])
	region = as.data.frame(region[region_order, ])
	n = nrow(region)
	dist = numeric(n)
		
	if(n < 2) {
		region = region_bk
		region$dist = NA
		return(region)
	}
		
	dist[1] = region[2, 1] - region[1, 2]
	dist[n] = region[n, 1] - region[n-1, 2]
			
	if(n > 2) {
		d1 = region[3:n, 1] - region[2:(n-1), 2]
		d2 = region[2:(n-1), 1] - region[1:(n-2), 2]
		if(mode == "min") {
			dist[2:(n-1)] = pmin(d1, d2)
		} else if(mode == "max") {
			dist[2:(n-1)] = pmax(d1, d2)
		} else if(mode == "mean") {
			dist[2:(n-1)] = apply(cbind(d1, d2), 1, mean)
		} else if (mode == "left") {
			dist[2:(n - 1)] = d2
		} else if (mode == "right") {
			dist[2:(n - 1)] = d1
		}
	}
		
	dist = ifelse(dist < 0, 0, dist)

	if(mode == "left") {
		dist[1] = NA
	} else if(mode == "right") {
		dist[n] = NA
	}

	region = region_bk
	region$dist[region_order] = dist

	if(normalize_to_width) {
		region$dist = region$dist/(region[, 3] - region[, 2])
	}
	
	return(region)
}

# == title
# Genomic position transformation function
#
# == param
# -region Genomic positions at a single chromosome. It is a data frame with two
#     columns which are start position and end position.
# -... other arguments
#
# == details
# The default position transformation functions transforms position to be equally distributed
# along the chromosome. If users want to define their own transformation function, the requirement
# is that the returned value should be a data frame with two columns: transformed start position
# and transformed end position. The returned value should have same number of rows as the input one.
#
# For details why need to use position transformation, please refer to `circos.genomicPosTransformLines`.
posTransform.default = function(region, ...) {
	xlim = get.cell.meta.data("xlim")
	segment = seq(xlim[1], xlim[2], length.out = nrow(region) + 1)
	return(data.frame(start = segment[-length(segment)], end = segment[-1]))
}

# == title
# Genomic position transformation function specifically for text
#
# == param
# -region Genomic positions at a single chromosome. It is a data frame with two
#     columns which are start position and end position.
# -y positions of texts
# -labels text labels
# -cex  text size
# -font  text font style
# -sector.index sector index
# -track.index track index
# -padding padding of text
# -extend extend to allow labels to be put in an region which is wider than the current chromosome.
#    The value should be a proportion value and the length is either one or two.
# -... other arguments
#
# == details
# This position transformation function is designed specifically for text.
# Under the transformation, texts will be as close as possible to the original positions.
posTransform.text = function(
	region, 
	y, 
	labels, 
	cex = 1, 
	font = par("font"),
	sector.index = get.cell.meta.data("sector.index"),
	track.index = get.cell.meta.data("track.index"), 
	padding = 0, 
	extend = 0, 
	...) {

	if(length(y) == 1) y = rep(y, nrow(region))
	if(length(labels) == 1) labels = rep(labels, nrow(region))

	if(length(extend) == 1) extend = rep(extend, 2)
	if(length(extend) > 2) extend = extend[1:2]

	od = order(region[, 1] + region[, 2])
	od_back = NULL
	od_back[od] = seq_along(od)
	region = region[od, ,drop = FALSE]
	y = y[od]
	labels = labels[od]

	text_height = strheight(labels, cex = cex, font = font)*(1+padding)
	
	d = circlize( (region[[1]] + region[[2]])/2, y, sector.index = sector.index, track.index = track.index)
	alpha1 = d[, "theta"] + as.degree(atan(text_height/2/d[, "rou"]))
	alpha2 = d[, "theta"] - as.degree(atan(text_height/2/d[, "rou"]))
	x1 = reverse.circlize(alpha1, d[, "rou"], sector.index = sector.index, track.index = track.index)[, "x"]
	x2 = reverse.circlize(alpha2, d[, "rou"], sector.index = sector.index, track.index = track.index)[, "x"]
	
	xlim = get.cell.meta.data("xlim", sector.index = sector.index, track.index = track.index)
	xrange = xlim[2] - xlim[1]
	xlim = c(xlim[1] - extend[1]*xrange, xlim[2] + extend[2]*xrange)

	x1_new = x1
	x2_new = x2
	l = x2 - x1 >= xlim[2] - xlim[1]; x1_new[l] = xlim[1]; x2_new[l] = xlim[2]
	l = x1 < xlim[1]; x1_new[l] = xlim[1]; x2_new[l] = x2[l] + xlim[1] - x1[l]
	l = x2 > xlim[2]; x1_new[l] = x1[l] - (x2[l] - xlim[2]); x2_new[l] = xlim[2]
	df = smartAlign(x1_new, x2_new, xlim = xlim)
	return(df[od_back, ,drop = FALSE])
}


# == title
# Add heatmaps for selected regions
#
# == param
# -bed A data frame in bed format, the matrix should be stored from the fourth column.
# -col Colors for the heatmaps. The value can be a matrix or a color mapping function generated by `colorRamp2`.
# -na_col Color for NA values.
# -numeric.column Column index for the numeric columns. The values can be integer index or character index.
#     By default it takes all numeric columns from the fourth column.
# -border Border of the heatmap grids.
# -border_lwd Line width for borders of heatmap grids.
# -border_lty Line style for borders of heatmap grids.
# -connection_height Height of the connection lines. If it is set to ``NULL``, no connection will be drawn.
#        Use `mm_h`/`cm_h`/`inches_h` to set a height in absolute unit.
# -line_col Color of the connection lines. The value can be a vector.
# -line_lwd Line width of the connection lines.
# -line_lty Line style of the connection lines.
# -heatmap_height Height of the heatmap track
# -side Side of the heatmaps. Is the heatmap facing inside or outside?
# -track.margin Bottom and top margins.
#
# == details
# The function visualizes heatmaps which correspond to a subset of regions in the genome.
# The correspondance between heatmaps and regions are identified by connection lines.
#
# The function actually creates two tracks, one track for the connection lines and one track
# for the heamtaps. The heatmaps always fill the whole track.
#
# == seealso
# https://jokergoo.github.io/circlize_book/book/high-level-genomic-functions.html#genomic-heatmap
#
# == example
# \donttest{
# circos.initializeWithIdeogram(plotType = c("labels", "axis"))
# bed = generateRandomBed(nr = 100, nc = 4)
# col_fun = colorRamp2(c(-1, 0, 1), c("green", "black", "red"))
# circos.genomicHeatmap(bed, col_fun, side = "inside", border = "white")
# circos.genomicHeatmap(bed, col_fun, side = "outside", 
#     line_col = as.numeric(factor(bed[[1]])))
# }
circos.genomicHeatmap = function(
	bed, 
	col, 
	na_col = "grey",
	numeric.column = NULL, 
	border = NA, 
	border_lwd = par("lwd"), 
	border_lty = par("lty"), 
	connection_height = mm_h(5),
	line_col = par("col"), 
	line_lwd = par("lwd"), 
	line_lty = par("lty"),
	heatmap_height = 0.15, 
	side = c("inside", "outside"), 
	track.margin = circos.par("track.margin")) {

	bed = validate_data_frame(bed)
	validate_region(bed, check_chr = TRUE)

	mat = bed[, -(1:3), drop = FALSE]
	if(is.null(numeric.column)) {
		numeric.column = which(apply(mat, 2, "is.numeric"))
		if(length(numeric.column) == 0) {
			stop_wrap("You don't have numeric columns in `bed`.")
		}
	} else {
		if(is.numeric(numeric.column)) {
			numeric.column = numeric.column - 3
		}
	}
	mat = mat[, numeric.column, drop = FALSE]

	mat = as.matrix(mat)
	bed = cbind(bed[, 1:3], mat)

	side = match.arg(side)
	if(missing(col)) {
		col = colorRamp2(seq(min(mat, na.rm = TRUE), max(mat, na.rm = TRUE), length.out = 3), c("blue", "#EEEEEE", "red"))
	}
	if(is.function(col)) {
		col = col(mat)
	}
	col[is.na(col)] = na_col
	if(length(border) == 1) {
		if(is.na(border)) {
			border = col
		}
	}
	if(length(border) == 1) {
		border = rep(border, length(col))
		dim(border) = dim(col)
	}
	if(side == "inside") {
		if(!is.null(connection_height)) {
			circos.genomicPosTransformLines(bed, posTransform = posTransform.default, 
			    horizontalLine = "top", track.height = connection_height,
			    track.margin = c(convert_height(1, "mm"), track.margin[2]), 
			    cell.padding = c(0, 0, 0, 0),
			    col = line_col, lwd = line_lwd, lty = line_lty)
		}
		circos.genomicTrackPlotRegion(bed, stack = TRUE, 
		    panel.fun = function(region, value, ...) {
		    	l = bed[, 1] == CELL_META$sector.index
		    	i = getI(...)
		        circos.genomicRect(region, value, col = col[l, i], 
		            border = border[l, i], lwd = border_lwd, lty = border_lty, 
		            posTransform = posTransform.default, ...)
			}, bg.border = NA, track.height = heatmap_height, track.margin = c(track.margin[1], 0),
			cell.padding = c(0, 0, 0, 0))
	} else {
		circos.genomicTrackPlotRegion(bed, stack = TRUE, 
		    panel.fun = function(region, value, ...) {
		    	l = bed[, 1] == CELL_META$sector.index
		    	i = getI(...)
		        circos.genomicRect(region, value, col = col[l, i], 
		            border = border[l, i], lwd = border_lwd, lty = border_lty, 
		            posTransform = posTransform.default, ...)
			}, bg.border = NA, track.height = heatmap_height, track.margin = c(0, track.margin[2]),
			cell.padding = c(0, 0, 0, 0))
		if(!is.null(connection_height)) {
			circos.genomicPosTransformLines(bed, posTransform = posTransform.default, 
			    direction = "outside", horizontalLine = "bottom", track.height = connection_height,
			    track.margin = c(track.margin[2], convert_height(1, "mm")), 
			    cell.padding = c(0, 0, 0, 0), 
			    col = line_col, lwd = line_lwd, lty = line_lty)
		}
	}
}

# == title
# Add labels to specified genomic regions
#
# == param
# -bed A data frame in bed format.
# -labels A vector of labels corresponding to rows in ``bed``.
# -labels.column If the label column is already in ``bed``, the index for this column in ``bed``.
# -facing fFacing of the labels. The value can only be ``"clockwise"`` or ``"reverse.clockwise"``.
# -niceFacing Whether automatically adjust the facing of the labels.
# -col Color for the labels.
# -cex Size of the labels.
# -font Font of the labels.
# -padding Padding of the labels, the value is the ratio to the height of the label.
# -connection_height Height of the connection track.
# -line_col Color for the connection lines.
# -line_lwd Line width for the connection lines.
# -line_lty Line type for the connectioin lines.
# -labels_height Height of the labels track.
# -side Side of the labels track, is it in the inside of the track where the regions are marked?
# -labels.side Same as ``side``. It will replace ``side`` in the future versions.
# -track.margin Bottom and top margins.
#
# == details
# The function adds labels for the specified regions. The positions of labels are arranged
# so that they are not overlapping to each other.
#
# This function creates two tracks, one for the connection lines and one for the labels.
#
# == seealso
# https://jokergoo.github.io/circlize_book/book/high-level-genomic-functions.html#labels
# 
# == example
# \donttest{
# circos.initializeWithIdeogram()
# bed = generateRandomBed(nr = 50, fun = function(k) sample(letters, k, replace = TRUE))
# bed[1, 4] = "aaaaa"
# circos.genomicLabels(bed, labels.column = 4, side = "inside")
# circos.clear()
#
# circos.initializeWithIdeogram(plotType = NULL)
# circos.genomicLabels(bed, labels.column = 4, side = "outside",
#     col = as.numeric(factor(bed[[1]])), line_col = as.numeric(factor(bed[[1]])))
# circos.genomicIdeogram()
# circos.clear()
# }
circos.genomicLabels = function(
	bed, 
	labels = NULL, 
	labels.column = NULL,
	facing = "clockwise", 
	niceFacing = TRUE,
	col = par("col"), 
	cex = 0.8, 
	font = par("font"), 
	padding = 0.4,
	connection_height = mm_h(5), 
	line_col = par("col"), 
	line_lwd = par("lwd"), 
	line_lty = par("lty"),
	labels_height = min(c(cm_h(1.5), max(strwidth(labels, cex = cex, font = font)))),
	side = c("inside", "outside"), 
	labels.side = side,
	track.margin = circos.par("track.margin")) {

	bed = validate_data_frame(bed)
	validate_region(bed, check_chr = TRUE)
	
	if(is.null(labels)){
		labels = bed[[labels.column]]
	}
	bed[[4]] = as.vector(labels)
	bed[[1]] = factor(as.vector(bed[[1]]), levels = get.all.sector.index())

	od = order(bed[[1]], bed[[2]])
	bed = bed[od, , drop = FALSE]
	bed[[1]] = as.vector(bed[[1]])

	# order `bed`

	if(length(col) == 1) col = rep(col, nrow(bed))
	if(length(cex) == 1) cex = rep(cex, nrow(bed))
	if(length(font) == 1) font = rep(font, nrow(bed))
	if(length(line_col) == 1) line_col = rep(line_col, nrow(bed))
	if(length(line_lwd) == 1) line_lwd = rep(line_lwd, nrow(bed))
	if(length(line_lty) == 1) line_lty = rep(line_lty, nrow(bed))

	col = col[od]
	cex = cex[od]
	font = font[od]
	line_col = line_col[od]
	line_lwd = line_lwd[od]
	line_lty = line_lty[od]

	if(!circos.par$xaxis.clock.wise) {
		all_chr_vec = as.character(bed[, 1])
		bed[, 2] = .CIRCOS.ENV$.SECTOR.DATA[all_chr_vec, "max.value"] - bed[, 2] + .CIRCOS.ENV$.SECTOR.DATA[all_chr_vec, "min.value"]
		bed[, 3] = .CIRCOS.ENV$.SECTOR.DATA[all_chr_vec, "max.value"] - bed[, 3] + .CIRCOS.ENV$.SECTOR.DATA[all_chr_vec, "min.value"]
		bed[, 2:3] = bed[, 3:2]
		circos.par(xaxis.clock.wise = TRUE)
		on.exit(circos.par(xaxis.clock.wise = FALSE))
	}

	chr = get.all.sector.index()[1]
	sector_data = get.sector.data(chr)
	
	## find the largest gap
	all_chr = unique(bed[, 1])
	rho = NULL
	for(cr in all_chr) {
		sub_bed = bed[bed[, 1] == cr, ]
		rho = c(rho, circlize(sub_bed[, 2], y = rep(1, nrow(sub_bed)), sector.index = cr)[, 1])
	}
	if(length(rho) == 0) {
		anchor = ((sector_data["start.degree"] + sector_data["end.degree"])/2 + 180) %% 360
	} else if(length(rho) == 1) {
		anchor = (rho + 180) %% 360
	} else {
		rho = sort(rho)
		rho = c(rho, rho[1] + 360)
		i = which.max(diff(rho))
		anchor = ((rho[i] + rho[i+1])/2) %% 360
	}

	# map all other chromosomes to the first chromosome
	chr_width = sector_data["start.degree"] - sector_data["end.degree"]
	# extend = (360 - chr_width)/chr_width
	# extend = c(0, extend)
	extend = numeric(2)
	s1 = sector_data["start.degree"] %% 360
	s2 = sector_data["end.degree"] %% 360
	# if the anchor in inside the first sector
	if(s1 < s2) { # the first sector go across theta = 0
		s1 = s1 + 360
	} 

	if(anchor >= s2 && anchor <= s1) { # anchor inside sector
		if(s1 - s2 > 180) {
			extend[1] = (abs(s1 - anchor) %% 360)/chr_width
			extend[2] = -(abs(s2 - anchor) %% 360)/chr_width
		} else {
			extend = (360 - chr_width)/chr_width
			extend = c(0, extend)
		}
	} else {
		extend[1] = ((anchor - s1) %% 360)/chr_width  # goes reverse clockwise
		extend[2] = ((s2 - anchor) %% 360)/chr_width  # goes clockwise
	}
	extend = abs(extend)

	if(0) segments(0, 0, cos(anchor/180*pi), sin(anchor/180*pi))

	# start.degree is always larger than end.degree
	new_chr_range = c(sector_data["start.degree"] + chr_width*extend[1], sector_data["end.degree"] - chr_width*extend[2])

	all_chr = unique(bed[, 1])
	bed2 = NULL
	for(cr in all_chr) {
		sub_bed = bed[bed[, 1] == cr, ]
		if(cr != chr) {
			dfx1 = circlize(sub_bed[, 2], y = rep(1, nrow(sub_bed)), sector.index = cr)
			dfx2 = circlize(sub_bed[, 3], y = rep(1, nrow(sub_bed)), sector.index = cr)
			degree_diff = dfx2[, 1] - dfx1[, 1]
			
			l = dfx1[, 1] > new_chr_range[1] | dfx1[, 1] < new_chr_range[2]
			if(any(l)) dfx1[l, 1] = (dfx1[l, 1] - new_chr_range[2]) %% 360 + new_chr_range[2]
			x1 = reverse.circlize(dfx1, sector.index = chr)[, 1]
			
			dfx2[, 1] = dfx1[, 1] + degree_diff
			x2 = reverse.circlize(dfx2, sector.index = chr)[, 1]
			sub_bed[, 2:3] = data.frame(start = x1, end = x2)
		}
		bed2 = rbind(bed2, sub_bed)
	}
	bed2[, 1] = chr

	if(!facing %in% c("clockwise", "reverse.clockwise")) {
		stop_wrap("facing can only be 'clockwise' or `reverse.clockwise`.")
	}

	op = circos.par("points.overflow.warning")
	circos.par("points.overflow.warning" = FALSE)
	labels.side = side = match.arg(side)[1]
	if(labels.side == "inside") {
		# an empty track
		circos.genomicTrackPlotRegion(bed2, ylim = c(0, 1), 
			track.margin = c(convert_height(0.5, "mm"), track.margin[2]), 
		    cell.padding = c(0, 0, 0, 0),
		    track.height = connection_height, bg.border = NA)
		i_track = get.cell.meta.data("track.index")

		# add labels
		circos.genomicTrackPlotRegion(bed2, ylim = c(0, 1), panel.fun = function(region, value, ...) {
			l = bed2[[1]] == CELL_META$sector.index
			circos.genomicText(region, value, y = 1, labels.column = 1, facing = facing, adj = c((facing == "clockwise") + 0, 0.5),
				posTransform = posTransform.text, col = col[l], cex = cex[l], font = font[l], niceFacing = niceFacing,
				padding = padding, extend = extend)
		}, track.height = labels_height, bg.border = NA, track.margin = c(track.margin[1], 0),
		cell.padding = c(0, 0, 0, 0))

		circos.genomicPosTransformLines(bed2, 
			posTransform = function(region, value) {
				l = bed2[[1]] == get.cell.meta.data("sector.index", track.index = i_track + 1)
				posTransform.text(region, y = 1, labels = value[[1]], cex = cex[l], font = font[l], 
					track.index = i_track + 1, padding = padding, extend = extend)
			},
			direction = "inside", horizontalLine = "top", track.index = i_track,
			cell.padding = c(0, 0, 0, 0),
			col = line_col, lwd = line_lwd, lty = line_lty
		)
	} else {
		circos.genomicTrackPlotRegion(bed2, ylim = c(0, 1), panel.fun = function(region, value, ...) {
			l = bed2[[1]] == CELL_META$sector.index
			circos.genomicText(region, value, y = 0, labels.column = 1, facing = facing, adj = c((facing == "reverse.clockwise") + 0, 0.5),
				posTransform = posTransform.text, col = col[l], cex = cex[l], font = font[l], niceFacing = niceFacing,
				padding = padding, extend = extend)
		}, track.height = labels_height, bg.border = NA, track.margin = c(0, track.margin[2]),
		cell.padding = c(0, 0, 0, 0))
		i_track = get.cell.meta.data("track.index")

		circos.genomicPosTransformLines(bed2, 
			posTransform = function(region, value) {
				l = bed2[[1]] == get.cell.meta.data("sector.index", track.index = i_track)
				posTransform.text(region, y = 0, labels = value[[1]], cex = cex[l], font = font[l], 
					track.index = i_track, padding = padding, extend = extend)
			},
			direction = "outside", horizontalLine = "bottom",
			col = line_col, lwd = line_lwd, lty = line_lty, track.height = connection_height,
			track.margin = c(track.margin[1], convert_height(0.5, "mm")), 
		    cell.padding = c(0, 0, 0, 0)
		)
	}
	circos.par("points.overflow.warning" = op)
}

# == title
# Adjust positions of text
#
# == param
# -x1 Position which corresponds to the top of the text.
# -x2 Position which corresponds to the bottom of the text.
# -xlim Ranges on x-axis.
#
# == details
# used internally
smartAlign = function(x1, x2, xlim) {
	
	ncluster.before = -1
	ncluster = length(x1)
	while(ncluster.before != ncluster) {
		ncluster.before = ncluster
		cluster = rep(0, length(x1))
		i_cluster = 1
		cluster[1] = i_cluster
		for(i in seq_along(x1)[-1]) {
			# overlap with previous one
			if(x1[i] <= x2[i-1]) {
				cluster[i] = i_cluster
			} else {
				i_cluster = i_cluster + 1
				cluster[i] = i_cluster
			}
		}
		ncluster = length(unique(cluster))
		
		if(ncluster.before == ncluster) break
		
		# tile intervals in each cluster and re-assign x1 and x2
		new_x1 = numeric(length(x1))
		new_x2 = numeric(length(x2))
		for(i_cluster in unique(cluster)) {
			index = which(cluster == i_cluster)
			len = x2[index] - x1[index]
			total_len = sum(len)
			mid = (min(x1[index]) + max(x2[index]))/2
			if(total_len > xlim[2] - xlim[1]) {
				tp = c(0, cumsum(len)) + (xlim[1] + xlim[2])/2 - total_len/2
			} else if(mid - total_len/2 < xlim[1]) {
				tp = c(0, cumsum(len)) + xlim[1]
			} else if(mid + total_len/2 > xlim[2]) {
				tp = xlim[2] - total_len + c(0, cumsum(len))
			} else {
				tp = c(0, cumsum(len)) + mid - total_len/2
			}
			new_x1[index] = tp[-length(tp)]
			new_x2[index] = tp[-1]
		}
		
		x1 = new_x1
		x2 = new_x2
	}
	
	return(data.frame(start = x1, end = x2))
}