Bug fixes and addition of geneAnnotationSearch Function

Author: jlmaier12

DESCRIPTION

@@ -17,10 +17,11 @@ Description:
   Automate the detection of gaps and elevations in mapped sequencing 
   read coverage using a 2D pattern-matching algorithm. 'ProActive' detects, characterizes 
   and visualizes read coverage patterns in both genomes and metagenomes. Optionally, 
-  users may provide gene predictions associated with their genome or metagenome 
+  users may provide gene annotations associated with their genome or metagenome 
   in the form of a .gff file. In this case, 'ProActive' will generate an additional 
-  output table containing the gene predictions found within the detected regions of gapped 
-  and elevated read coverage. 
+  output table containing the gene annotations found within the detected regions of gapped 
+  and elevated read coverage. Additionally, users can search for gene 
+  annotations of interest in the output read coverage plots.
 License: GPL-2
 Encoding: UTF-8
 LazyData: true

NAMESPACE

@@ -1,14 +1,17 @@
 # Generated by roxygen2: do not edit by hand
 
-export(ProActive)
+export(ProActiveDetect)
+export(geneAnnotationSearch)
 export(plotProActiveResults)
 import(ggplot2)
+importFrom(dplyr,"%>%")
 importFrom(dplyr,bind_rows)
 importFrom(stats,median)
 importFrom(stats,na.omit)
 importFrom(stringr,regex)
 importFrom(stringr,str_detect)
 importFrom(stringr,str_extract)
 importFrom(stringr,str_extract_all)
+importFrom(stringr,str_which)
 importFrom(utils,capture.output)
 importFrom(utils,write.table)

NEWS.md

@@ -1,3 +1,12 @@
-# ProActive 0.0.1
+## Changes in 0.1.2
++ Changed name of main function from `ProActive()` to `ProActiveDetect()`
++ Added functionality for searching read coverage plots for gene annotations of 
+interest with the geneAnnotationSearch() function
++ Fixed bug with pattern-match start and stop positions on genome/contig chunks
 
-* Initial CRAN submission.
+## Changes in 0.0.2
++ Changes for CRAN re-submission after review
+
+## Changes in 0.0.1
+
++ Initial CRAN submission.

R/GPsInElevGaps.R

@@ -13,10 +13,12 @@
 #' @param chunkContigs TRUE or FALSE, If TRUE and `mode`="metagenome", contigs longer
 #' than the `chunkSize` will be 'chunked' into smaller subsets and pattern-matching
 #' will be performed on each subset. Default is FALSE.
+#' @param chunkSize If `mode`="genome" OR if `mode`="metagenome" and `chunkContigs`=TRUE,
+#' chunk the genome or contigs, respectively, into smaller subsets for pattern-matching.
 #' @importFrom stringr str_detect str_extract_all regex
 #' @importFrom dplyr bind_rows
 #' @keywords internal
-GPsInElevGaps <- function(elevGapSummList, windowSize, gffTSV, mode, chunkContigs) {
+GPsInElevGaps <- function(elevGapSummList, windowSize, gffTSV, mode, chunkContigs, chunkSize) {
   colnames(gffTSV) <- c("seqid", "source", "type", "start", "end", "score", "strand", "phase", "attributes")
   if (TRUE %in% (str_detect(gffTSV[,9], regex('product', ignore_case = T)))){
     product <- str_extract_all(gffTSV [,9], regex("(?<=product=)[\\s\\S]*",ignore_case = T))
@@ -26,13 +28,17 @@ GPsInElevGaps <- function(elevGapSummList, windowSize, gffTSV, mode, chunkContig
     trueRefName <- elevGapSummList[[i]][[8]]
     if(mode == "metagenome"){
       refName <- elevGapSummList[[i]][[8]]
-    if(chunkContigs==TRUE){
+    if(chunkContigs == TRUE){
       refName <- gsub("_chunk_.*", "", refName)
     }
-    gffTSV <- gffTSV[which(gffTSV[, 1] == refName), ]
+      gffTSV <- gffTSV[which(gffTSV[, 1] == refName), ]
+      startPos <- elevGapSummList[[i]][[4]] * windowSize
+      endPos <- elevGapSummList[[i]][[5]] * windowSize
+    } else {
+      chunkNumber <- as.numeric(str_extract(trueRefName, "(?<=\\_)\\d+$")) - 1
+      startPos <- (elevGapSummList[[i]][[4]] * windowSize) + (chunkNumber * chunkSize)
+      endPos <- (elevGapSummList[[i]][[5]] * windowSize) + (chunkNumber * chunkSize)
     }
-    startPos <- elevGapSummList[[i]][[4]] * windowSize
-    endPos <- elevGapSummList[[i]][[5]] * windowSize
     matchRegion <- seq(startPos, endPos, 1)
     GPs <- gffTSV[which(gffTSV[, 5] %in% matchRegion), ]
     GPs$seqid <- rep(trueRefName, nrow(GPs))

R/ProActive-package.R

@@ -1,16 +1,21 @@
 #' @title \bold{ProActive}
 #'
-#' @description Automatic detection of gaps and elevations in mapped sequencing reads using read coverage
-#'   pattern-matching. Pattern-matching results are summarized into a dataframe and users can visualize the
-#'   pattern-matches and associated read coverage patterns. If a .gff file is supplied, ProActive will output
-#'   the ORFs found within the detected regions of elevated and/or gapped read coverage. ProActive works with both
-#'   metagenomes and genomes.
+#' @description `ProActive` automatically detects regions of gapped and elevated read coverage
+#' using a 2D pattern-matching algorithm. `ProActive` detects, characterizes and
+#' visualizes read coverage patterns in both genomes and metagenomes. Optionally,
+#' users may provide gene annotations associated with their genome or metagenome
+#' in the form of a .gff file. In this case, `ProActive` will generate an additional
+#' output table containing the gene annotations found within the detected regions of
+#' gapped and elevated read coverage. Additionally, users can search for gene
+#' annotations of interest in the output read coverage plots.
 #'
-#' @details The two main functions in ProActive are:
+#' @details The three main functions in `ProActive` are:
 #' \enumerate{
-#' \item \code{\link{ProActive}} performs the pattern-matching and characterization of read coverage patterns.
+#' \item \code{\link{ProActiveDetect}} performs the pattern-matching and characterization of read coverage patterns.
 #' \item \code{\link{plotProActiveResults}} plots the results from
-#'  \code{ProActive()}
+#'  \code{ProActiveDetect()}
+#' \item \code{\link{geneAnnotationSearch}} searches classified contigs/chunks for gene annotations that match
+#' user-provided keywords.
 #' }
 #'
 #' @author Jessie Maier \email{[email protected]}

R/ProActive.R R/ProActiveDetect.RR/ProActive.R renamed to R/ProActiveDetect.R

@@ -40,12 +40,12 @@
 #' @return A list containing 6 objects described in the function description.
 #' @export
 #' @examples
-#' metagenome_results <- ProActive(
+#' metagenome_results <- ProActiveDetect(
 #'   pileup = sampleMetagenomePileup,
 #'   mode = "metagenome",
 #'   gffTSV = sampleMetagenomegffTSV
 #' )
-ProActive <- function(pileup, mode, gffTSV, windowSize = 1000, chunkContigs = FALSE,
+ProActiveDetect <- function(pileup, mode, gffTSV, windowSize = 1000, chunkContigs = FALSE,
                       minSize = 10000, maxSize = Inf, minContigLength = 30000,
                       chunkSize = 100000, IncludeNoPatterns = FALSE, verbose = TRUE,
                       saveFilesTo) {
@@ -83,11 +83,11 @@ ProActive <- function(pileup, mode, gffTSV, windowSize = 1000, chunkContigs = FA
   }
   filteredOutContigsDf <- patternMatchSummary[[2]]
   if(verbose){message("Summarizing pattern-matching results")}
-  summaryTable <- classifSumm(pileup, patternMatchSummary[[1]], windowSize, mode)
+  summaryTable <- classifSumm(pileup, patternMatchSummary[[1]], windowSize, mode, chunkSize)
   if (missing(gffTSV) == FALSE) {
     if(verbose){message("Finding gene predictions in elevated or gapped regions of read coverage...")}
     elevGapSummList <- removeNoPatterns(patternMatchSummary[[1]])
-    GPSummTable <- GPsInElevGaps(elevGapSummList, windowSize, gffTSV, mode, chunkContigs)
+    GPSummTable <- GPsInElevGaps(elevGapSummList, windowSize, gffTSV, mode, chunkContigs, chunkSize)
   }
   if(verbose){message("Finalizing output")}
   endTime <- Sys.time()

R/bestMatchListFunctions.R

@@ -8,16 +8,29 @@
 #' 100 bp windows/bins.
 #' @param windowSize The number of basepairs to average read coverage values over.
 #' Options are 100, 200, 500, 1000 ONLY. Default is 1000.
+#' @param chunkSize If `mode`="genome" OR if `mode`="metagenome" and `chunkContigs`=TRUE,
+#' chunk the genome or contigs, respectively, into smaller subsets for pattern-matching.
 #' @param mode Either "genome" or "metagenome"
 #' @keywords internal
-classifSumm <- function(pileup, bestMatchList, windowSize, mode) {
+classifSumm <- function(pileup, bestMatchList, windowSize, mode, chunkSize) {
   if (length(bestMatchList) == 0) {
     stop("No pattern-matches detected")
   }
   refName <- vapply(seq_along(bestMatchList), function(i) {bestMatchList[[i]][[8]]}, character(1))
   elevRatio <- vapply(seq_along(bestMatchList), function(i) {bestMatchList[[i]][[6]]}, numeric(1))
-  startPos <- vapply(seq_along(bestMatchList), function(i) {bestMatchList[[i]][[4]]} * windowSize, numeric(1))
-  endPos <- vapply(seq_along(bestMatchList), function(i) {bestMatchList[[i]][[5]]} * windowSize, numeric(1))
+  startEndPos <- lapply(seq_along(bestMatchList), function(i) {
+    if(mode == "genome"){
+      refName <- bestMatchList[[i]][[8]]
+      chunkNumber <- as.numeric(str_extract(refName, "(?<=\\_)\\d+$")) - 1
+      startPos <- (bestMatchList[[i]][[4]] * windowSize) + (chunkNumber * chunkSize)
+      endPos <- (bestMatchList[[i]][[5]] * windowSize) + (chunkNumber * chunkSize)
+    } else {
+      startPos <- bestMatchList[[i]][[4]] * windowSize
+      endPos <- bestMatchList[[i]][[5]] * windowSize
+    }
+    list(startPos, endPos)})
+  startPos <- vapply(seq_along(startEndPos), function(i){startEndPos[[i]][[1]]}, numeric(1))
+  endPos <- vapply(seq_along(startEndPos), function(i){startEndPos[[i]][[2]]}, numeric(1))
   classification <- vapply(seq_along(bestMatchList), function(i) {bestMatchList[[i]][[7]]}, character(1))
   matchSize <- vapply(seq_along(bestMatchList), function(i) {
     pileupSubset <- pileup[which(pileup[, 1] == bestMatchList[[i]][[8]]), ]

R/chunkingFunctions.R

@@ -39,14 +39,14 @@ contigChunks <- function(pileup, chunkSize) {
         refNameChunk[c(remainIdx:nrow(pileupSubset))] <- rep(paste0(refName, "_chunk_", Z + 1), nrow(contigChunk))
       }
       chunkedPileupSubset <- cbind.data.frame(refNameChunk, pileupSubset$coverage, pileupSubset$position)
-      colnames(chunkedPileupSubset) <- c("refName", "coverage", "position")
       NAIdx <- which(is.na(chunkedPileup[, 1]))[[1]]
       chunkedPileup[c(NAIdx:(NAIdx + nrow(pileupSubset) - 1)), ] <- chunkedPileupSubset
     } else {
       NAIdx <- which(is.na(chunkedPileup[, 1]))[[1]]
       chunkedPileup[c(NAIdx:(NAIdx + nrow(pileupSubset) - 1)), ] <- pileupSubset
     }
   }
+  colnames(chunkedPileup) <- c("refName", "coverage", "position")
   return(chunkedPileup)
 }
 

R/geneAnnotationPlot.R

@@ -0,0 +1,62 @@
+#' Gene annotation plot
+#'
+#' Plot read coverage and location of gene annotations that match the keywords and
+#' search criteria for contig/chunk currently being assessed
+#'
+#' @param geneAnnotSubset Subset of gene annotations to be plotted
+#' @param keywords The key-word(s) used for the search.
+#' @param pileupSubset A subset of the pileup associated with the contig/chunk being assessed
+#' @param colIdx The column index 'gene' or 'product' column
+#' @param startbpRange The basepair at which the search is started if a 'specific' search is used
+#' @param endbpRange The basepair at which the search is ended if a 'specific' search is used
+#' @param elevRatio The maximum/minimum values of the pattern-match
+#' @param pattern The pattern-match information associated with the contig/chunk being assessed
+#' @param windowSize The number of basepairs to average read coverage values over.
+#' @keywords internal
+#' @importFrom stringr str_which
+geneAnnotationPlot <- function(geneAnnotSubset, keywords, pileupSubset,
+                               colIdx, startbpRange, endbpRange, elevRatio,
+                               pattern, windowSize, chunkSize, mode){
+  position <- coverage <- start <- NULL
+  classification <- pattern[[7]]
+  refName <- pileupSubset[1, 1]
+  if(mode == "genome"){
+    chunkNumber <- as.numeric(str_extract(refName, "(?<=\\_)\\d+$")) - 1
+    startPos <- (pattern[[4]] * windowSize) + (chunkNumber * chunkSize)
+    endPos <- (pattern[[5]] * windowSize) + (chunkNumber * chunkSize)
+  } else {
+    startPos <- pattern[[4]] * windowSize
+    endPos <- pattern[[5]] * windowSize
+  }
+  matchIdxs <- str_which(geneAnnotSubset[, colIdx], regex(paste(keywords, collapse = "|"), ignore_case = TRUE))
+  geneAnnotMatches <- geneAnnotSubset[matchIdxs,]
+  geneStartPos <- geneAnnotMatches$start
+  geneAnnotLabels <- paste0("#", c(1:nrow(geneAnnotMatches)), ": ", geneAnnotMatches[, colIdx],
+                            sep = " ", collapse = " \n ")
+  plot <- ggplot(data = pileupSubset, aes(x = position, y = coverage)) +
+    geom_area(fill = "#009E73") +
+    geom_vline(xintercept = geneStartPos, linewidth = 1) +
+    geom_vline(xintercept = c(startPos, endPos), color = "#D55E00", linewidth = 1) +
+    geom_vline(xintercept = c(startbpRange, endbpRange), color = "#D55E00", linewidth = 1, linetype = "dotted") +
+    geom_label(data = geneAnnotMatches, aes(x = start, y = (max(pileupSubset$coverage) / 2),
+                                           label = paste0("#", c(1 : nrow(geneAnnotMatches)))),
+               size = 2.75) +
+    scale_x_continuous(expand = c(0, 0)) +
+    theme(panel.grid.major = element_blank(),
+          panel.grid.minor = element_blank(),
+          panel.background = element_blank(),
+          axis.line = element_line(colour = "black"),
+          text = element_text(size = 15),
+          plot.margin = margin(
+            t = 0,
+            r = 10,
+            b = 0,
+            l = 2
+          )) +
+    labs(title = paste(refName, classification),
+         subtitle = paste("elevation ratio:", round(elevRatio, digits = 4)),
+         x = "Basepair position",
+         caption = geneAnnotLabels,
+         y = "Read coverage")
+  return(plot)
+}