R/VDJ_circos.R

In Platypus: Single-Cell Immune Repertoire and Gene Expression Analysis

#' Plots a Circos diagram from an adjacency matrix. Uses the Circlize chordDiagram function. Is called by VDJ_clonotype_clusters_circos(), VDJ_alpha_beta_Vgene_circos() and VDJ_VJ_usage_circos() functions or works on its own when supplied with an adjacency matrix.
#' @param Adj_matrix Adjacency matrix to be plotted. Rownames and Colnames correspond to genes to be matched and entries determine the weight of the connection between the genes (eg. number of clonotypes expressing these two genes).
#' @param group Named list of genes, with list elements corresponding to group-names, and element names being the gene-names. Is generated by VDJ_VJ_usage and VDJ_alpha_beta_Vgene_circos.
#' @param grid.col Named list of genes, with list elements corresponding to color and element names being gene-names. If not supplied it is generated randomly within the function. Is also generated by VDJ_VJ_usage and VDJ_alpha_beta_Vgene_circos.
#' @param label.threshold Genes are only labeled if the count is larger then the label.threshold. By default all label.threshold = 0 (all genes are labeled).
#' @param c.count.label Boolean, lets the user decide if the gene and count labels should be plotted or not. Default = T.
#' @param c.count.label.size Determines the font size of the gene labels. By default the font size for count labels is 0.6.
#' @param gene.label Boolean, lets the user decide if the gene labels should be plotted or not.
#' @param gene.label.size Determines the font size of the gene labels. By default the labelsize is automatically adjusted to 0.7 for labels with two or less digits, 0.6 for labels between 2 and 6 digits, and 0.4 for all longer labels. A manually defined font size will be the same for all labels!
#' @param platy.theme Allows plotting in the new "pretty" theme or the older "spiky" theme without group labels and radial arrangement of gene.labels. Default = "pretty".
#' @param arr.col Data.frame with three columns where the first two indicate the names of genes, clonotypes or clusters to be connected, and the third corresponds to the color of the arrow. Default set to data.frame(c("dummy.clonotype"), c("dummy.cluster"), c("dummy.color")), so no arrow is drawn.
#' @param arr.direction Either 1 or -1 and determines the direction of the arrow. Default=1.
#' @param axis Option to choose the count axis for each gene. "default", "percent" or "max"  possible. Default: "max".
#' @return Returns a circos plot.
#' @export
#' @examples
#' \dontrun{
#' # manual replotting of Circos plot:
#' VDJ_circos(Adj_matrix =  VDJ_alpha_beta_Vgene_circos_output[[2]][[1]], 
#'   grid.col = VDJ_alpha_beta_Vgene_circos_output[[3]], 
#'   group = VDJ_alpha_beta_Vgene_circos_output[[4]], 
#'   c.count.label.size = 0.4, 
#'   gene.label.size = 0.5, 
#'   arr.col = data.frame(c("TRBV10"),c("TRBJ2-7"), c("black")), 
#'   axis="percent")
#' 
#'}

VDJ_circos <- function(Adj_matrix, platy.theme, group, grid.col, label.threshold, axis, c.count.label, arr.col, arr.direction, gene.label.size, gene.label, c.count.label.size){

  #### Initailize Variables ####
  
  df <- NULL
  df1 <- NULL
  ylim <- NULL
  xcenter <- NULL
  ycenter <- NULL
  gene.label.length <- NULL
  sector.index <- NULL
  sector.name <- NULL
  label.cex <- NULL
  sub_group <- NULL
  sector_indices <- NULL
  xplot <- NULL
  by <- NULL
  xrange <- NULL

  
  
  if(missing(Adj_matrix)){print("Error: No Adjacency matrix supplied.")}
  if(missing(axis)){axis <- "max"}
  if(missing(label.threshold)){label.threshold <- 0}
  if(missing(grid.col)){grid.col <- stats::setNames(grDevices::rainbow(length(union(rownames(Adj_matrix), colnames(Adj_matrix)))),sample(union(rownames(Adj_matrix), colnames(Adj_matrix))))}
  if(missing(c.count.label)){c.count.label <- T}
  if(missing(c.count.label.size)){c.count.label.size <- 0.6}
  if(missing(arr.col)){arr.col <- data.frame(c("dummy1"), c("dummy2"), c(""))}
  if(missing(arr.direction)){arr.direction <- 1}
  if(missing(gene.label.size)){gene.label.size <- "undef"}
  if(missing(gene.label)){gene.label <- T}
  if(missing(platy.theme)){platy.theme <- "pretty"}
  
  
  if(platy.theme == "pretty"){
    
  #### Rearrange Adj Matrix ####
  
    df <- reshape2::melt(Adj_matrix)
    df1 <- data.frame(Var = c(rownames(Adj_matrix), colnames(Adj_matrix)), value= c(rowSums(Adj_matrix), colSums(Adj_matrix)))
    df1 <- df1[order(df1$value),] #does not really work if group argument is used at the same time...
    
    gene.label.length <- max(nchar(gsub("^\\D+", "", gsub("TRAV", "", gsub("TRBV", "",  gsub("TRAJ", "", gsub("TRBJ", "", gsub("IGKV", "", gsub("IGLV", "", gsub("IGHV", "", gsub("IGKJ", "", gsub("IGLJ", "", gsub("IGHJ", "", append(dimnames(Adj_matrix)[[1]], dimnames(Adj_matrix)[[2]])))))))))))))) # Determine length of longest gene label

  #### Draw Circos-Plot ####
  
    circlize::circos.clear()
    circlize::circos.par(points.overflow.warning=FALSE)
    
    circlize::chordDiagram(df,
                 link.sort = T,
                 link.decreasing = T,
                 order = df1$Var,
                 group = group,
                 directional = arr.direction, 
                 direction.type = c("arrows"),
                 link.arr.col = arr.col,
                 link.arr.length = 0.2,
                 annotationTrack = c("grid"),
                 annotationTrackHeight = c(0.02),
                 preAllocateTracks = list(
                   list(track.height = circlize::mm_h(0.2), track.margin = c(0.01, 0)),  # track 1: margin around plot
                   list(track.height = circlize::mm_h(0.25), track.margin = c(0.01, 0)),  # track 2: line indicating groups
                   list(track.height = circlize::mm_h(0.8+gene.label.length*0.6), track.margin = c(0.01, 0)),    # track 3: gene label
                   list(track.height = circlize::mm_h(3.5), track.margin = c(0.01, 0))),   # track 4: clone count
                 grid.col = grid.col)
    # circlize::circos.info()
   

  #### Add Gene/Cluster labels to circos plot ####
    if(gene.label){
      circlize::circos.track(track.index = 3, panel.fun = function(x, y) {
        if(circlize::get.cell.meta.data("xrange") > label.threshold){
          ycenter <- circlize::get.cell.meta.data("ycenter")
          xcenter <- circlize::get.cell.meta.data("xcenter")
          sector.index <- circlize::get.cell.meta.data("sector.index")
          
          sector.index <- gsub("TRAV", "", sector.index)
          sector.index <- gsub("TRBV", "", sector.index)
          sector.index <- gsub("TRAJ", "", sector.index)
          sector.index <- gsub("TRBJ", "", sector.index)
          sector.index <- gsub("IGKV", "", sector.index)
          sector.index <- gsub("IGLV", "", sector.index)
          sector.index <- gsub("IGHV", "", sector.index)
          sector.index <- gsub("IGKJ", "", sector.index)
          sector.index <- gsub("IGLJ", "", sector.index)
          sector.index <- gsub("IGHJ", "", sector.index)
          sector.index <- gsub("^\\D+", "", sector.index)
          
          if(gene.label.size=="undef"){
            if(nchar(sector.index)>2){
              label.cex <- 0.6
            }
            if(nchar(sector.index)>6){
              label.cex <- 0.4
            }
            if(nchar(sector.index)<=2){
              label.cex <- 0.7
            }
          }else{label.cex <- gene.label.size}
          
          circlize::circos.text(xcenter, ycenter-circlize::mm_h(7+gene.label.length*0.1), sector.index,
                                facing = "clockwise", niceFacing = TRUE, col="black", font = 2, cex=label.cex) 
        }else{
          # ylim <- circlize::get.cell.meta.data("ylim")
          # xcenter <- circlize::get.cell.meta.data("xcenter")
          # sector.index <- circlize::get.cell.meta.data("sector.index")
          # circlize::circos.text(xcenter, ylim[1], sector.index,
          #                       facing = "clockwise", niceFacing = TRUE, adj = c(0, 0.6), cex=0.5, col="white")
        }
      }, bg.border = NA) # here set bg.border to NA is important
      
    }
    
    
    #### Add Group labels to circos plot ####
    
    for(p in unique(unname(group))) {
      
      sub_group = names(group[group == p])
      sector_indices <- circlize::get.all.sector.index() 
      sub_group <- sub_group[which(sub_group %in% sector_indices)]
      
      circlize::highlight.sector(sector.index = sub_group, 
                                 track.index = 2,
                                 col = "black",
                                 text = p,
                                 font = 2,
                                 cex = 0.8,
                                 text.vjust = -0.5,
                                 facing = "bending.inside",
                                 niceFacing = TRUE)
    }
    
    
    #### Add axis to circos plot ####
    
    if(axis=="basic"){
      for(si in circlize::get.all.sector.index()) {
        circlize::circos.axis(h = "top", labels.cex = c.count.label.size, sector.index = si, track.index = 5)
      }
    }
    if(axis=="percent"){
      circlize::circos.track(track.index = 4, panel.fun = function(x, y) {
        xlim = circlize::get.cell.meta.data("xlim")
        ylim = circlize::get.cell.meta.data("ylim")
        sector.name = circlize::get.cell.meta.data("sector.index")
        xplot = circlize::get.cell.meta.data("xplot")
        
        circlize::circos.lines(xlim, c(min(ylim), min(ylim)), lty = 4) # dotted line
        by = ifelse(abs(xplot[2] - xplot[1]) > 30, 0.2, 0.5)
        for(p in seq(by, 1, by = by)) {
          circlize::circos.text(p*(xlim[2] - xlim[1]) + xlim[1], min(ylim) + 0.1,
                                paste0(p*100, "%"), cex = c.count.label.size, adj = c(0.5, 0), niceFacing = TRUE)
        }
      }, bg.border = NA)
    }
    if(axis=="max"){
      for(si in circlize::get.all.sector.index()) {
        circlize::circos.axis(h = "top", labels =  FALSE, sector.index = si, track.index = 5)
      }
      circlize::circos.track(track.index = 5, panel.fun = function(x, y) {
        xlim = circlize::get.cell.meta.data("xlim")
        ylim = circlize::get.cell.meta.data("ylim")
        sector.name = circlize::get.cell.meta.data("sector.index")
        xrange = circlize::get.cell.meta.data("xrange")
        
        #circos.lines(xlim, c(min(ylim), min(ylim)), lty = 1)
        if(c.count.label){
          if(circlize::get.cell.meta.data("xrange")> label.threshold){
            circlize::circos.text(mean(xlim), max(ylim) + 0.8,
                                  round(xrange,0), cex = c.count.label.size, adj = c(0.5, 0), niceFacing = TRUE)
          }
        }
      }, bg.border = NA)
    }
    
    
  }

  if(platy.theme == "spiky"){
    
    #### Rearrange Adj Matrix ####
    
    df <- reshape2::melt(Adj_matrix)
    df1 <- data.frame(Var = c(rownames(Adj_matrix), colnames(Adj_matrix)), value= c(rowSums(Adj_matrix), colSums(Adj_matrix)))
    df1 <- df1[order(df1$value),] #does not really work if group argument is used at the same time...
    
    gene.label.length <- max(nchar(gsub("^\\D+", "", gsub("TRAV", "", gsub("TRBV", "",  gsub("TRAJ", "", gsub("TRBJ", "", gsub("IGKV", "", gsub("IGLV", "", gsub("IGHV", "", gsub("IGKJ", "", gsub("IGLJ", "", gsub("IGHJ", "", append(dimnames(Adj_matrix)[[1]], dimnames(Adj_matrix)[[2]])))))))))))))) # Determine length of longest gene label
    
    #### Draw Circos-Plot ####
    
    circlize::circos.clear()
    circlize::circos.par(points.overflow.warning=FALSE)
    
    circlize::chordDiagram(df,
                           link.sort = T,
                           link.decreasing = T,
                           order = df1$Var,
                           group = group,
                           directional = arr.direction, 
                           direction.type = c("arrows"),
                           link.arr.col = arr.col,
                           link.arr.length = 0.2,
                           annotationTrack = c("grid"),
                           preAllocateTracks = list(track.height = circlize::mm_h(5), track.margin = c(circlize::mm_h(5), 0)),
                           grid.col = grid.col)
    # circlize::circos.info()
    
    
    #### Add Gene/Cluster labels to circos plot ####
    
    if(gene.label.size == "undef"){
      gene.label.size <- 0.5
    }
    
    if(gene.label){
        circlize::circos.track(track.index = 1, panel.fun = function(x, y) {
          if(circlize::get.cell.meta.data("xrange")> label.threshold){
            ylim <- circlize::get.cell.meta.data("ylim")
            xcenter <- circlize::get.cell.meta.data("xcenter")
            sector.index <- circlize::get.cell.meta.data("sector.index")
            circlize::circos.text(xcenter, ylim[1], sector.index,
                                  facing = "clockwise", niceFacing = TRUE, adj = c(0, 0.6), cex=gene.label.size, col="black")
        }else{
          # ylim <- circlize::get.cell.meta.data("ylim")
          # xcenter <- circlize::get.cell.meta.data("xcenter")
          # sector.index <- circlize::get.cell.meta.data("sector.index")
          # circlize::circos.text(xcenter, ylim[1], sector.index,
          #                       facing = "clockwise", niceFacing = TRUE, adj = c(0, 0.6), cex=0.5, col="white")
        }
      }, bg.border = NA) # here set bg.border to NA is important
      
    
    
    
    #### Add Group labels to circos plot ####
  
      # no group labels
    
    #### Add axis to circos plot ####
    
    if(axis=="basic"){
      for(si in circlize::get.all.sector.index()) {
        circlize::circos.axis(h = "top", labels.cex = c.count.label.size, sector.index = si, track.index = 2)
      }
    }
    if(axis=="percent"){
      circlize::circos.track(track.index = 1, panel.fun = function(x, y) {
        xlim = circlize::get.cell.meta.data("xlim")
        ylim = circlize::get.cell.meta.data("ylim")
        sector.name = circlize::get.cell.meta.data("sector.index")
        xplot = circlize::get.cell.meta.data("xplot")
        
        circlize::circos.lines(xlim, c(min(ylim), min(ylim)), lty = 4) # dotted line
        by = ifelse(abs(xplot[2] - xplot[1]) > 30, 0.2, 0.5)
        for(p in seq(by, 1, by = by)) {
          circlize::circos.text(p*(xlim[2] - xlim[1]) + xlim[1], min(ylim) + 0.1,
                                paste0(p*100, "%"), cex = c.count.label.size, adj = c(0.5, 0), niceFacing = TRUE)
        }
      }, bg.border = NA)
    }
    if(axis=="max"){
      for(si in circlize::get.all.sector.index()) {
        circlize::circos.axis(h = "top", labels =  FALSE, sector.index = si, track.index = 2)
      }
      circlize::circos.track(track.index = 2, panel.fun = function(x, y) {
        xlim = circlize::get.cell.meta.data("xlim")
        ylim = circlize::get.cell.meta.data("ylim")
        sector.name = circlize::get.cell.meta.data("sector.index")
        xrange = circlize::get.cell.meta.data("xrange")
        
        #circos.lines(xlim, c(min(ylim), min(ylim)), lty = 1)
        if(c.count.label){
          if(circlize::get.cell.meta.data("xrange")> label.threshold){
            circlize::circos.text(mean(xlim), max(ylim) + 0.5,
                                  round(xrange,0), cex = c.count.label.size, adj = c(0.5, 0), niceFacing = TRUE)
          }
        }
      }, bg.border = NA)
    }
    
    }
    
  }
  
  #### Return Plot ####
  
  return(plot)
}