HiCAGE: Hi-C Annotation and Graphics Ensemble

#' Generate circos plots from overlap data output
#'
#' @param datatable The GRanges datatable generated by overlap function
#' @param plot.subset Allows for plotting interactions with Hi-C scores only.
#' Set to TRUE if no RNAseq or gene expression file was used to generate the
#' overlap output.
#' @param display.legend Adds legend to circos plot
#' @param hic.legend Title for Hi-C legend
#' @param hic.range A vector representing the low and high values for the legend
#'     of Hi-C values
#' @param rna.legend Title for RNA legend
#' @param rna.range A vector representing the low and high values for the legend
#'     of RNA values
#' @param circos.color Passes colors to chord diagram for each group
#' @param ... Addtional arguments passed to chordDiagram function
#' @importFrom circlize chordDiagram circos.clear circos.rect
#' circos.trackPlotRegion colorRamp2 circos.par get.cell.meta.data
#' circos.text
#' @importFrom dplyr left_join
#' @importFrom stats setNames aggregate
#' @importFrom plotrix color.gradient color.legend
#' @import grDevices
#' @export
#' @return Generates a chord diagram containing circluar tracks with an inner
#' track displaying  average Hi-C interaction score and an inner track with
#' displaying average RNA-seq log(FPKM) values for each segment in the plot
#' @examples
#' hic_chr20 <- system.file("extdata", "hic_chr20.txt", package = "HiCAGE")
#' segment_chr20 <- system.file("extdata", "segment_chr20.bed",
#'     package = "HiCAGE")
#' rna_chr20 <- system.file("extdata", "rna_chr20.tsv", package = "HiCAGE")
#' example <- overlap(hicfile = hic_chr20,
#'          segmentfile = segment_chr20,
#'          rnafile = rna_chr20,
#'          rna.columns = c(1, 7))
#' circleplot(example)

circleplot <- function(datatable,
                       plot.subset = FALSE,
                       display.legend = TRUE,
                       hic.legend = "Avg. Hi-C Score",
                       hic.range = c(0, 130),
                       rna.legend = "Avg. log(FPKM)",
                       rna.range = c(0, 3.2),
                       circos.color = NULL,
                       ...) {

  if (length(unique(datatable$mark1 == 1))) {
    datatable <- datatable[,c(8:14, 1:7, 15)]
    datatable <- setNames(datatable, c("region1chrom",
                                       "region1start",
                                       "region1end",
                                       "mark1",
                                       "segscore1",
                                       "gene1",
                                       "logFPKM1",
                                       "region2chrom",
                                       "region2start",
                                       "region2end",
                                       "mark2",
                                       "segscore2",
                                       "gene2",
                                       "logFPKM2",
                                       "HiCscore"))
  }

  finaltable <- (table(datatable$mark1, datatable$mark2))

  #Data frame for circos plot
  circosR <- data.frame(finaltable)
  circosR <- setNames(circosR, c("mark1", "mark2", "Freq"))

  #Compute average scores for each group of interactions
  scoretable <- aggregate(HiCscore~mark1 + mark2, data = datatable, FUN = mean)
  colnames(scoretable)[3] <- "Avg.Score"

  #Generate table for Heatmap
  circosR$mark1 <- as.character(circosR$mark1)
  circosR$mark2 <- as.character(circosR$mark2)
  heatmap <- left_join(circosR, scoretable, by = c("mark1", "mark2"))

  if (plot.subset == TRUE) {
  circos.clear()
  par(pty="s")
  circos.par(points.overflow.warning=FALSE,
             track.margin = c(0, 0))

  chordDiagram(circosR,
               circos.par(track.margin = c(0.005, 0.005)),
               grid.col = circos.color,
               annotationTrack = "grid",
               preAllocateTracks = 2,
               ...)

  circos.trackPlotRegion(track.index = 1, panel.fun = function(x, y) {
    xlim = get.cell.meta.data("xlim")
    ylim = get.cell.meta.data("ylim")
    sector.name = get.cell.meta.data("sector.index")
    circos.text(mean(xlim), ylim[1],
                sector.name,
                facing = "clockwise",
                niceFacing = TRUE,
                adj = c(0.3, 0.3),
                cex = 0.8)
  }
  , bg.border = NA)
  #Add heatmap to chord diagram
  col_fun = colorRamp2(hic.range,
                       c("white", "red"))

  tab1 <- heatmap[order(heatmap$mark1,rev(heatmap$mark2)),]
  tab1[is.na(tab1)] <- 0
  tab2 <- heatmap[order(heatmap$mark2,rev(heatmap$mark1)),]
  tab2[is.na(tab2)] <- 0
  tab3 <- tab1[,c(2,1,3,4)]
  colnames(tab3) [1] <- "mark1"
  colnames(tab3) [2] <- "mark2"
  tab3$Freq <- 0
  tab1 <- rbind(tab1, tab3)
  tab2 <- rbind(tab2, tab3)

  circos.trackPlotRegion(ylim = c(0, 3),
                         circos.par("track.height" = 0.01),
                         bg.border = NA,
                         track.index = 2,
                         panel.fun = function(x, y) {
                           end <- 0
                           sector = get.cell.meta.data("sector.index")
                           aux.tab1 <- tab1[tab1$mark1 == sector,]
                           aux.tab2 <- tab2[tab2$mark2 == sector,]
                           col_tab1 = col_fun(aux.tab1$Avg.Score)
                           col_tab2 = col_fun(aux.tab2$Avg.Score)
                           for(i in 1:nrow(aux.tab1)){
                             start <- end
                             end <- end + aux.tab1[i,"Freq"]
                             circos.rect(start, 0, end, 2,
                                         border = col_tab1[i],
                                         col = col_tab1[i],
                                         sector.index = sector,
                                         track.index = 2)
                           }
                           for(i in 1:nrow(aux.tab2)){
                             start <- end
                             end <- end + aux.tab2[i,"Freq"]
                             circos.rect(start, 0, end, 2,
                                         border = col_tab2[i],
                                         col = col_tab2[i],
                                         sector.index = sector,
                                         track.index = 2)
                           }
                           end <- 0
                         })
  }
  else {
  #Draw chord diagram
  FPKM1score <- aggregate(logFPKM1~mark1 + mark2, data = datatable, FUN = mean)
  FPKM2score <- aggregate(logFPKM2~mark1 + mark2, data = datatable, FUN = mean)
  heatmap <- left_join(heatmap, FPKM1score, by = c("mark1", "mark2"))
  heatmap <- left_join(heatmap, FPKM2score, by = c("mark1", "mark2"))
  circos.clear()
  par(pty="s")
  circos.par(points.overflow.warning=FALSE,
             track.margin = c(0, 0))

  chordDiagram(circosR,
               circos.par(track.margin = c(0.005, 0.005)),
               grid.col = circos.color,
               annotationTrack = "grid",
               preAllocateTracks = 3,
               ...)

  circos.trackPlotRegion(track.index = 1, panel.fun = function(x, y) {
    xlim = get.cell.meta.data("xlim")
    ylim = get.cell.meta.data("ylim")
    sector.name = get.cell.meta.data("sector.index")
    circos.text(mean(xlim), ylim[1],
                sector.name,
                facing = "clockwise",
                niceFacing = TRUE,
                adj = c(0.3, 0.3),
                cex = 0.8)
    }
    , bg.border = NA)
  #Add heatmap to chord diagram
  col_fun = colorRamp2(hic.range,
                       c("white", "red"))

  tab1 <- heatmap[order(heatmap$mark1,rev(heatmap$mark2)),]
  tab1[is.na(tab1)] <- 0
  tab2 <- heatmap[order(heatmap$mark2,rev(heatmap$mark1)),]
  tab2[is.na(tab2)] <- 0
  tab3 <- tab1[,c(2,1,3:6)]
  colnames(tab3) [1] <- "mark1"
  colnames(tab3) [2] <- "mark2"
  tab3$Freq <- 0
  tab3$Avg.Score <- hic.range[1]
  tab3$logFPKM1 <- 0
  tab3$logFPKM2 <- 0
  tab1 <- rbind(tab1, tab3)
  tab2 <- rbind(tab2, tab3)

  circos.trackPlotRegion(ylim = c(0, 3),
                         circos.par("track.height" = 0.01),
                         bg.border = NA,
                         track.index = 3,
                         panel.fun = function(x, y) {
                           end <- 0
                           sector = get.cell.meta.data("sector.index")
                           aux.tab1 <- tab1[tab1$mark1 == sector,]
                           aux.tab2 <- tab2[tab2$mark2 == sector,]
                           col_tab1 = col_fun(aux.tab1$Avg.Score)
                           col_tab2 = col_fun(aux.tab2$Avg.Score)
                           for(i in 1:nrow(aux.tab1)){
                             start <- end
                             end <- end + aux.tab1[i,"Freq"]
                             circos.rect(start, 0.01, end, 2,
                                         border = col_tab1[i],
                                         col = col_tab1[i],
                                         sector.index = sector,
                                         track.index = 3)
                           }
                           for(i in 1:nrow(aux.tab2)){
                             start <- end
                             end <- end + aux.tab2[i,"Freq"]
                             circos.rect(start, 0.01, end, 2,
                                         border = col_tab2[i],
                                         col = col_tab2[i],
                                         sector.index = sector,
                                         track.index = 3)
                           }
                           end <- 0
                         })

  col_fun2 = colorRamp2(rna.range,
                        c("white", "black"))

  circos.trackPlotRegion(ylim = c(0, 3),
                         circos.par("track.height" = 0.01),
                         bg.border = NA,
                         track.index = 2,
                         panel.fun = function(x, y) {
                           end <- 0
                           sector = get.cell.meta.data("sector.index")
                           aux.tab1 <- tab1[tab1$mark1 == sector,]
                           aux.tab2 <- tab2[tab2$mark2 == sector,]
                           col_tab1 = col_fun2(aux.tab1$logFPKM1)
                           col_tab2 = col_fun2(aux.tab2$logFPKM2)
                           for(i in 1:nrow(aux.tab1)){
                             start <- end
                             end <- end + aux.tab1[i,"Freq"]
                             circos.rect(start, -0.9, end, 1.3,
                                         border = col_tab1[i],
                                         col = col_tab1[i],
                                         sector.index = sector,
                                         track.index = 2)
                           }
                           for(i in 1:nrow(aux.tab2)){
                             start <- end
                             end <- end + aux.tab2[i,"Freq"]
                             circos.rect(start, -0.9, end, 1.3,
                                         border = col_tab2[i],
                                         col = col_tab2[i],
                                         sector.index = sector,
                                         track.index = 2)
                           }
                           end <- 0
                         })
  }

  hiccolor <- color.gradient(c(1,1), c(1,0), c(1,0), nslices = 200)
  rnacolor <- color.gradient(c(1,0), c(1,0), c(1,0), nslices = 200)

  if (display.legend == TRUE & plot.subset == FALSE) {
    color.legend(xl = -0.98,
                 yb = -1.0,
                 xr = -0.50,
                 yt = -0.91,
                 legend = rna.range,
                 rect.col= rnacolor)
    color.legend(xl = 0.50,
                 yb = -1.0,
                 xr = 0.98,
                 yt = -0.91,
                 legend = hic.range,
                 rect.col= hiccolor)
    text(-0.74, -1.045, labels = rna.legend)
    text(0.74, -1.045, labels = hic.legend)
  }
  if (display.legend == "RNA") {
    color.legend(xl = -0.98,
                 yb = -1.0,
                 xr = -0.50,
                 yt = -0.91,
                 legend = rna.range,
                 rect.col= rnacolor)
    text(-0.74, -1.045, labels = rna.legend)
  }
  if (display.legend == TRUE | plot.subset == TRUE) {
    color.legend(xl = 0.50,
                 yb = -1.0,
                 xr = 0.98,
                 yt = -0.91,
                 legend = hic.range,
                 rect.col= hiccolor)
    text(0.74, -1.045, labels = hic.legend)
  }

}