R/snvinGene.R

In nriddiford/mutationProfiles: Collection of functions to explore SNV data (e.g. Mutect2/Freebayes) parsed by get_tris.pl

#' snvinGene
#'
#' Plot all snvs found in a given gene
#' @description Plot all snvs found in a given gene.
#' A 'gene_lengths' file must be provided with the following fields (cols 1..6 required)
#' gene length chrom    start      end      tss scaling_factor
#' This can be generated using the script 'script/genomic_features.pl' and a genome .gtf file
#' @param gene_lengths File containing all genes and their lengths (as generated by 'script/genomefeatures.pl') [Default 'data/gene_lengths.txt']
#' @param gene2plot Name of the gene to plot
#' @import ggplot2 dplyr
#' @keywords gene
#' @export

snvinGene <- function(..., snv_data=NULL, gene_lengths="data/gene_lengths.txt", gene2plot='kuz', annotated=TRUE, col_by_status=TRUE, write=FALSE){
  if(missing(snv_data)){
    snv_data<-getData(...)
  }

  gene_lengths <- read.delim(gene_lengths, header = T)

  region <- gene_lengths %>%
    dplyr::filter(gene == gene2plot) %>%
    droplevels()

  gene_length <-(region$end-region$start)
  wStart<-(region$start - gene_length/10)
  wEnd<-(region$end + gene_length/10)
  wChrom<-as.character(region$chrom)
  wTss<-suppressWarnings(as.numeric(levels(region$tss))[region$tss])

  snv_data <- snv_data %>%
    dplyr::filter(chrom == wChrom & pos >= wStart & pos <= wEnd)

  if(nrow(snv_data) == 0){
    stop(paste("There are no snvs in", gene2plot, "- Exiting", "\n"))
  }

  snv_data$colour_var <- snv_data$feature

  if(annotated){
    snv_data$colour_var <- snv_data$variant_type
    if(col_by_status)
      snv_data$colour_var <- snv_data$status
  }

  p <- ggplot(snv_data)
  p <- p + geom_point(aes(pos/1000000, sample, colour = colour_var, alpha = af, size = 1.5), position=position_jitter(width=0, height=0.1))
  p <- p + guides(size = FALSE, sample = FALSE)
  p <- p + cleanTheme() +
    theme(axis.title.y=element_blank(),
          panel.grid.major.y = element_line(color="grey80", size = 0.5, linetype = "dotted"),
          axis.text.y = element_text(size = 30)
    )
  p <- p + scale_x_continuous("Mbs", expand = c(0,0), breaks = seq(round(wStart/1000000, digits = 2),round(wEnd/1000000, digits = 2),by=0.05), limits=c(wStart/1000000, wEnd/1000000))
  p <- p + annotate("rect", xmin=region$start/1000000, xmax=region$end/1000000, ymin=0, ymax=0.3, alpha=.2, fill="skyblue")
  p <- p + geom_vline(xintercept = wTss/1000000, colour="red", alpha=.7, linetype="solid")

  p <- p + geom_segment(aes(x = wTss/1000000, y = 0, xend= wTss/1000000, yend = 0.1), colour="red")
  middle<-((wEnd/1000000+wStart/1000000)/2)
  p <- p + annotate("text", x = middle, y = 0.15, label=gene2plot, size=6)
  p <- p + ggtitle(paste("Chromosome:", wChrom))
  if(write){
    hit_gene<-paste(gene2plot, "_hits.pdf", sep='')
    cat("Writing file", hit_gene, "\n")
    ggsave(paste("plots/", hit_gene, sep=""), width = 10, height = 10)
  }
  p
}