R/JunctionTools.R

In mffrank/RNAseqToCustomFASTA: Produce Custom Protein Fasta Databases from RNAseq Data

#' Generate a Junction Summary Report
#' @details Saves a table of the types of the Splice junctions to working directory
#' @param sharedjun GRanges object containing all splice junctions
#' @param splicemax Splice Matrix from annotation folder
#' @param ids Id mapping from annotation folder
#' @param txdb Txdb object from annotation folder
#' @param print boolean, should the result be printed in working directory as pdf
#' or be returned as a list.
#' @export

filterSpacePlot <- function(covfilter_unique = c(1,2,3,4,5,10), share_num = c(1,2), bed_directory, splicemax, txdb, ids){

  combinations <- expand.grid(covfilter_unique, share_num)
  colnames(combinations) <- c("covfilter_unique", "share_num")
  combinations$remaining_junctions <- 0
  juncs_list <- list()
  message("Importing Junctions with different filters \n \n")
  for(i in 1:nrow(combinations)){
    juncs <- getStarJunctions(bed_directory, splicemax = splicemax, txdb = txdb, ids = ids,
                              pattern = "\\.tab", share_num = combinations[i,2], extend = 100, skip = 0,
                              covfilter_unique = combinations[i,1], covfilter_multi = 0, overhang = "max_overhang")
    combinations[i,3] <- length(juncs)
    juncs_list[[i]] <- juncs
  }

  m <- ggplot(combinations, aes(x=covfilter_unique, y = remaining_junctions, colour = share_num)) +
    geom_point()

  return(m)
}

#' Plot the distribution of your splice junction types
#' @param junction_type generated by the function: JunctionType
#' @param pdf boolean, should the result be printed in working directory as pdf
#' or be returned as a list.
#' @param outfile Name of the printed pdf
#' @export

junctionSummaryReport <- function(junction_type, pdf = TRUE, outfile = "Junction_summary_report.pdf"){
  simple_types <- as.factor(junction_type$jun_type)
  levels(simple_types) <- list(Novel = "connect a known exon and a non-exon region",
                               Novel = "connect a known exon and a region overlap with known exon",
                               Novel = "connect a region overlap with known exon and a non-exon region",
                               Novel = "connect two non-exon region",
                               Novel = "connect two regions overlaped with known exons",
                               Gene_fusion = "gene fusion",
                               Known ="known junction",
                               AltSplice_junction = "novel alternative splicing junction")
  junction_type$Junction_type <- simple_types
  junction_type$AA_Impact <- junction_type$Junction_type == "AltSplice_junction"
  m <- ggplot(junction_type) +
    geom_bar(aes(x = Junction_type, fill = AA_Impact), stat = "count") +
    scale_y_log10() +
    theme(axis.text.x = element_text(angle=30, hjust=1, vjust=1))
  if(pdf){
    pdf(file = outfile, width = 6, height=4)
    print(m)
    dev.off()
  }else{
    return(m)
  }

}

#' #' Generate a Junction Summary Report
#' #' @details Saves a table of the types of the Splice junctions to working directory
#' #' @param junction_type generated by the function: JunctionType
#' #' @export
#'
#' junctionSummaryReport <- function(junction_type, save = TRUE){
#'
#'   output <- list(junction_type = table(junction_type[, 'jun_type']))
#'   if(save){
#'     saveRDS(output,"VariantSummary.rda")
#'     # pdf("Variant_summary_report.pdf", width = 6, height = 4)
#'     # for(i in 1:length(output)){
#'     #   print(output[[i]])
#'     # }
#'     # dev.off()
#'   } else{
#'     return(output)
#'   }
#' }


#' Plot the overlap of identified Junctions between replicates
#' @param bed_directory Path to folder with all replicate vcf files
#' @param pattern pattern to match
#' @import data.table
#' @import VennDiagram
#' @export

plotJuncReplicateOverlap <- function(bed_directory, pattern = "\\.tab", outfile = "SJ_Replicate_overlap.tiff"){
  files <- list.files(bed_directory, pattern = pattern)
  juncs <- lapply(paste0(bed_directory,"/",files), function(x) fread(x, sep = "\t"))
  vennlist <- lapply(juncs, function(x) paste(x$V1, x$V2, x$V3, sep ="-"))
  names(vennlist) <- paste0("Rep. ", seq_along(vennlist))
  venn.diagram(vennlist,
               filename = outfile,
               print.mode = "raw", fill = c("red", "green"), alpha = 0.5,
               cex = 2.5, margin = 0.06, cat.cex = 2.5, cat.dist = c(0.06,0.06), lty ="blank")

}