R/plotSE.R

In PoisonAlien/peakseason: Functions for faster ChIP/ATAC-seq analysis and visualization.

#' Draw hocky-stick plot from super enhancer results.
#' @param enhToGene ENHANCER_TO_GENE.txt file generated by ROSE_GeneMapper.
#' @param markGenes SE genes to highlight
#' @param keepHighestSE If multiple genes with name are found to be in SE list, default is to keep top ranking entry. Default TRUE.
#' @export

plotSE = function(enhToGene, markGenes = NULL, keepHighestSE = TRUE){

  enhToGene = data.table::fread(input = enhToGene)
  colnames(enhToGene) = c("Enahncer", "chr", "start", "end", "n_constituent", "size", "K27Ac_signal", "Input_signal",
                          "enhancerRank", "OVERLAP_GENES", "PROXIMAL_GENES", "Hugo_Symbol", "Rank", "isSuper")
  enhToGene = enhToGene[order(as.numeric(Rank))]
  enhToGene[, ip_vector := as.numeric(K27Ac_signal) - as.numeric(Input_signal)]
  enhToGene$ip_vector = ifelse(test = enhToGene$ip_vector < 0, yes = 0, no = enhToGene$ip_vector)
  enhToGene[, plotRank := rev(Rank)]

  x.cutoff = max(enhToGene[,Rank]) - rev(enhToGene[isSuper == 1][, Rank])[1]
  y.cutoff = rev(enhToGene[isSuper == 1][, ip_vector])[1]

  enhToGene.gg = ggplot(data = enhToGene, aes(x = plotRank, y = ip_vector, color = as.character(isSuper), label = Hugo_Symbol))+geom_line(col = "gray70")+
    geom_hline(yintercept = y.cutoff, col = "gray70")+geom_vline(xintercept = x.cutoff, col = "gray70")+
    geom_point(size = 1)+theme_classic()+
    scale_color_manual(values = c('1' = 'black', '0' = 'gray70'))+theme(legend.position = "none")+xlab("Enhancer Rank")+ylab("Signal")+
    ggtitle(label = paste0("#SE: ", nrow(enhToGene[isSuper == 1])), subtitle = paste0("SE cutoff: ", round(y.cutoff, digits = 2)))

  if(!is.null(markGenes)){
    rnklbl = enhToGene[isSuper == 1][Hugo_Symbol %in% markGenes]
    if(keepHighestSE){
      rnklbl = rnklbl[!duplicated(Hugo_Symbol)]
      rnklbl[, Hugo_Symbol := paste0(Hugo_Symbol, " [", Rank, "]")]
    }
    enhToGene.gg = enhToGene.gg+geom_point(data = rnklbl, col = "red")+
      ggrepel::geom_text_repel(data = rnklbl, size = 3, nudge_x = -6000,force = 6, nudge_y = 7000, segment.colour = 'brown', segment.alpha = 0.4)
  }

  enhToGene.gg
}