R/cicero_functions.R

In smorabit/scicat:

#' plot_gene_connections
#'
#' Plot cicero connections that flank the coordinates of a given gene in a given .gtf
#' @param cur_gene name of the gene you want to plot
#' @param connections cicero connections with some extra columns
#' @param gene_anno dataframe containing gtf info
#' @param bp_extend number of basepairs to plot +- relevant peaks. Also the number of bp +- the gene to flank (may change this)
#' @param ... parameters passed to cicero plot_connections function
#' @keywords stats
#' @export
flanking_connections <- function(cur_gene, connections, gene_anno, bp_extend=500, verbose=FALSE, ...) {

  # look for the gene in the gtf (what do if the gene is not there?)
  temp <- subset(gene_anno, gene_name == cur_gene)
  cur_position <- list(temp$seqid[1], temp$start[1] - bp_extend, temp$end[length(temp$end)] + bp_extend)

  # subset by connections where peak1 or peak2 is in the region of interest
  sub_connections1 <- subset(connections, chr == cur_position[[1]] & peak1_start >= cur_position[[2]] & peak1_end <= cur_position[[3]])
  sub_connections2 <- subset(connections, chr == cur_position[[1]] & peak2_start >= cur_position[[2]] & peak2_end <= cur_position[[3]])
  sub_connections <- rbind(sub_connections1, sub_connections2)

  # find lowest and highest coordinates within sub_connections
  chrom <- as.character(cur_position[[1]])
  hi_coord <- max(c(sub_connections$peak1_end, sub_connections$peak2_end))
  lo_coord <- min(c(sub_connections$peak1_start, sub_connections$peak2_start))

  return(list(chrom, lo_coord, hi_coord))
}