Update for SNV/INDEL; options via elipses

Author: nriddiford

NAMESPACE

@@ -1,16 +1,24 @@
 # Generated by roxygen2: do not edit by hand
 
+export(EnrichmentVolcano)
+export(bpSim)
 export(chromDist)
 export(cleanTheme)
+export(dist2motif2)
+export(distOverlay2)
 export(featureEnrichment)
 export(featuresHit)
 export(genTris)
 export(geneEnrichment)
 export(geneHit)
+export(generateData)
 export(getData)
+export(getTriFromFasta)
 export(mutSigs)
 export(mutSpectrum)
+export(plotdistanceOverlay2)
 export(rainfall)
+export(readBed)
 export(samplesPlot)
 export(setCols)
 export(sigTypes)
@@ -19,12 +27,22 @@ export(snvStats)
 export(snvinGene)
 export(triFreq)
 export(tssDist)
-import(BSgenome)
 import(BSgenome.Dmelanogaster.UCSC.dm6)
+import(Biostrings)
 import(GenomicRanges)
 import(RColorBrewer)
+import(colorspace)
+import(cowplot)
 import(data.table)
 import(deconstructSigs)
 import(dplyr)
+import(forcats)
 import(ggplot2)
+import(ggpubr)
+import(plotly)
+import(plyr)
 import(reshape)
+import(scales)
+import(tidyr)
+importFrom(data.table,as.data.table)
+importFrom(data.table,fread)

R/featureEnrichment.R

@@ -10,33 +10,36 @@
 #' @param features File containing total genomic lengths of features [Default 'data/genomic_features.txt']
 #' @param genome_length The total legnth of the genome [Default 118274340 (mappable regions on chroms 2, 3, 4, X & Y for Drosophila melanogastor Dmel6.12)]
 #' @keywords enrichment
-#' @import dplyr
-#' @return A data frame with FC scores for all genes seen at least n times in snv data
-#' @export 
+#' @import dplyr ggpubr
+#' @return A snv_data frame with FC scores for all genes seen at least n times in snv snv_data
+#' @export
+featureEnrichment <- function(..., snv_data=NULL, features='data/genomic_features.txt', genome_length=118274340, write=FALSE){
 
-featureEnrichment <- function(features='data/genomic_features.txt', genome_length=118274340){
   genome_features<-read.delim(features, header = T)
-  data<-getData()
-  mutCount<-nrow(data)
-  
+
+  if(missing(snv_data)){
+    snv_data<-getData(...)
+  }
+  mutCount<-nrow(snv_data)
+
   # To condense exon counts into "exon"
-  data$feature<-as.factor(gsub("exon_.*", "exon", data$feature))
-  
-  classCount<-table(data$feature)
+  snv_data$feature<-as.factor(gsub("exon_.*", "exon", snv_data$feature))
+
+  classCount<-table(snv_data$feature)
   classLengths<-setNames(as.list(genome_features$length), genome_features$feature)
-  
+
   fun <- function(f) {
     # Calculate the fraction of geneome occupied by each feature
     featureFraction<-classLengths[[f]]/genome_length
-    
-    # How many times should we expect to see this feature hit in our data (given number of obs. and fraction)?
+
+    # How many times should we expect to see this feature hit in our snv_data (given number of obs. and fraction)?
     featureExpect<-(mutCount*featureFraction)
-    
-    # observed/expected 
+
+    # observed/expected
     fc<-classCount[[f]]/featureExpect
-    fc<-round(fc,digits=1)
+    Log2FC<-log2(fc)
     featureExpect<-round(featureExpect,digits=3)
-    
+
     # Binomial test
     if(!is.null(classLengths[[f]])){
       if(classCount[f] >= featureExpect){
@@ -47,17 +50,63 @@ featureEnrichment <- function(features='data/genomic_features.txt', genome_lengt
         stat<-binom.test(x = classCount[f], n = mutCount, p = featureFraction, alternative = "less")
         test<-"depletion"
       }
-      sig_val<-'F'
-      if(stat$p.value <= 0.05){ sig_val<-'T'}
+      sig_val <- ifelse(stat$p.value <= 0.001, "***",
+                        ifelse(stat$p.value <= 0.01, "**",
+                               ifelse(stat$p.value <= 0.05, "*", "")))
+
       p_val<-format.pval(stat$p.value, digits = 3, eps=0.0001)
-      list(feature = f, observed = classCount[f], expected = featureExpect, fc = fc, test = test, sig = sig_val, p_val = p_val)
+      list(feature = f, observed = classCount[f], expected = featureExpect, Log2FC = Log2FC, test = test, sig = sig_val, p_val = p_val)
     }
   }
-  
-  enriched<-lapply(levels(data$feature), fun)
+
+  enriched<-lapply(levels(snv_data$feature), fun)
   enriched<-do.call(rbind, enriched)
   featuresFC<-as.data.frame(enriched)
   # Sort by FC value
-  featuresFC<-arrange(featuresFC,desc(as.integer(fc)))
-  return(featuresFC)
+  featuresFC<-dplyr::arrange(featuresFC,desc(abs(as.numeric(Log2FC))))
+  featuresFC$Log2FC<-round(as.numeric(featuresFC$Log2FC), 1)
+
+  if(write){
+    featuresFC <- filter(featuresFC, observed >= 5)
+    first.step <- lapply(featuresFC, unlist)
+    second.step <- as.data.frame(first.step, stringsAsFactors = F)
+    ggpubr::ggtexttable(second.step, rows = NULL, theme = ttheme("mGreen"))
+    feat_enrichment_table <- paste("feature_enrichment_table.tiff")
+    cat("Writing to file: ", 'plots/', feat_enrichment_table, sep = '')
+    ggsave(paste("plots/", feat_enrichment_table, sep=""), width = 5.5, height = (nrow(featuresFC)/3), dpi=300)
+  } else{
+    return(featuresFC)
+  }
+}
+
+
+featureEnrichmentPlot <- function(write=FALSE) {
+  feature_enrichment<-featureEnrichment()
+
+  feature_enrichment$feature <- as.character(feature_enrichment$feature)
+  feature_enrichment$Log2FC <- as.numeric(feature_enrichment$Log2FC)
+
+  feature_enrichment <- transform(feature_enrichment, feature = reorder(feature, -Log2FC))
+
+  feature_enrichment <- filter(feature_enrichment, observed >= 5)
+
+  # Custom sorting
+  # feature_enrichment$feature <- factor(feature_enrichment$feature, levels=c("intron", "intergenic", "exon", "3UTR", "ncRNA", "5UTR"))
+
+  p<-ggplot(feature_enrichment)
+  p<-p + geom_bar(aes(feature, Log2FC, fill = as.character(test)), stat="identity")
+  p<-p + guides(fill=FALSE)
+  p<-p + ylim(-2,2)
+  p<-p + cleanTheme() +
+    theme(panel.grid.major.y = element_line(color="grey80", size = 0.5, linetype = "dotted"),
+          axis.text.x = element_text(angle = 45, hjust=1),
+          axis.text = element_text(size=20)
+    )
+
+  if(write){
+    feat_plot <- paste("feat_plot.pdf")
+    cat("Writing file", feat_plot, "\n")
+    ggsave(paste("plots/", feat_plot, sep=""), width = 5, height = 10)
+  }
+  p
 }

R/featuresHit.R

@@ -5,30 +5,39 @@
 #' @keywords features
 #' @export
 
+featuresHit <- function(..., write=FALSE){
+  snv_data<-getData(...)
 
-featuresHit <- function(){
-  data<-getData()
-  
   # To condense exon counts into "exon"
-  data$feature<-as.factor(gsub("exon_.*", "exon", data$feature))
-  
+  snv_data$feature<-as.factor(gsub("exon_.*", "exon", snv_data$feature))
+
   # Reoders descending
-  data$feature<-factor(data$feature, levels = names(sort(table(data$feature), decreasing = TRUE)))
-  
-  #cols<-set_cols(data, "feature")
-  
-  p<-ggplot(data)
-  p<-p + geom_bar(aes(feature, fill = feature))
+  snv_data$feature<-factor(snv_data$feature, levels = names(sort(table(snv_data$feature), decreasing = TRUE)))
+
+  snv_data <- snv_data %>%
+    dplyr::group_by(feature) %>%
+    dplyr::add_tally() %>%
+    ungroup() %>%
+    dplyr::filter(n >= 5) %>%
+    droplevels()
+  #cols<-setCols(snv_data, "feature")
+
+  p <- ggplot(snv_data)
+  p <- p + geom_bar(aes(feature, fill = feature))
   #p<-p + cols
-  p<-p + cleanTheme() +
+  p <- p + cleanTheme() +
     theme(axis.title.x=element_blank(),
           panel.grid.major.y = element_line(color="grey80", size = 0.5, linetype = "dotted"))
-  p<-p + scale_x_discrete(expand = c(0.01, 0.01))
-  p<-p + scale_y_continuous(expand = c(0.01, 0.01))
-  
-  features_outfile<-paste("hit_features_count.pdf")
-  cat("Writing file", features_outfile, "\n")
-  ggsave(paste("plots/", features_outfile, sep=""), width = 20, height = 10)
-  
+  p <- p + scale_x_discrete(expand = c(0.01, 0.01))
+  p <- p + scale_y_continuous(expand = c(0.01, 0.01))
+
+  # colour to a pub palette:
+  # p<-p + ggpar(p, palette = 'jco')
+
+  if(write){
+    features_outfile<-paste("hit_features_count.pdf")
+    cat("Writing file", features_outfile, "\n")
+    ggsave(paste("plots/", features_outfile, sep=""), width = 20, height = 10)
+  }
   p
-}
+}