Updates for CRAN submission

Author: jlmaier12

.Rbuildignore

@@ -6,3 +6,4 @@
 ^docs$
 ^pkgdown$
 ^\.github$
+^data-raw$

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: ProActive
 Title: 
   Detect Elevations and Gaps in Mapped Sequencing Read Coverage
-Version: 0.0.0.9000
+Version: 0.0.1
 Maintainer: Jessie Maier <jlmaier@ncsu.edu>
 Authors@R: 
   c(person(given = "Jessie",
@@ -17,9 +17,9 @@ Description:
   The `ProActive` R package automatically detects regions of gapped and elevated 
   read coverage using a pattern-matching algorithm. `ProActive` can detect, characterize 
   and visualize read coverage patterns in both genomes and metagenomes. Optionally, 
-  users may provide gene annotations associated with their genome or metagenome 
+  users may provide gene predictions 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 ORFs found within the detected regions of gapped 
+  output table containing the gene predictions found within the detected regions of gapped 
   and elevated read coverage. 
 License: GPL-2
 Encoding: UTF-8

NAMESPACE

@@ -8,6 +8,7 @@ importFrom(stats,median)
 importFrom(stats,na.omit)
 importFrom(stringr,regex)
 importFrom(stringr,str_detect)
+importFrom(stringr,str_extract)
 importFrom(stringr,str_extract_all)
 importFrom(utils,capture.output)
 importFrom(utils,write.table)

NEWS.md

@@ -1,3 +1,3 @@
-# ProActive 0.0.1.9000
+# ProActive 0.0.1
 
 * Initial CRAN submission.

R/FullElevGapFunctions.R

@@ -20,18 +20,23 @@ fullElevGap <- function(pileupSubset, windowSize, minSize, maxSize, elevOrGap) {
   minCovSteps <- (abs((minReadCov + quarterReadCov) - minReadCov)) / 2
   maxCoverages <- seq((minReadCov + (quarterReadCov * 3)), maxReadCov, maxCovSteps)
   minCoverages <- seq(minReadCov, (minReadCov + quarterReadCov), minCovSteps)
-  elevLength <- ifelse((nrow(pileupSubset) - (10000 / windowSize)) > (maxSize / windowSize), maxSize / windowSize, nrow(pileupSubset) - (10000 / windowSize))
+  elevLength <- ifelse((nrow(pileupSubset) - (10000 / windowSize)) > (maxSize / windowSize),
+                       maxSize / windowSize,
+                       nrow(pileupSubset) - (10000 / windowSize))
   nonElev <- nrow(pileupSubset) - (elevLength + (5000 / windowSize))
   maxOrMin <- ifelse(elevOrGap == "Elevation", minReadCov, maxReadCov)
   maxOrMin2 <- ifelse(elevOrGap == "Elevation", maxReadCov, minReadCov)
   pattern <- c(rep(maxOrMin, 5000 / windowSize), rep(maxOrMin2, elevLength), rep(maxOrMin, nonElev))
   bestMatchInfo <- collectBestMatchInfo(pattern, pileupSubset, elevOrGap, "Full")
-  lapply(seq_along(maxCoverages), function(maxCov) {
-    bestMatchInfo <<- fullElevGapShrink(minReadCov, windowSize, maxCoverages[[maxCov]], elevLength, nonElev, bestMatchInfo, pileupSubset, minSize, elevOrGap)
-    lapply(seq_along(minCoverages), function(minCov) {
-      bestMatchInfo <<- fullElevGapShrink(minCoverages[[minCov]], windowSize, maxCoverages[[maxCov]], elevLength, nonElev, bestMatchInfo, pileupSubset, minSize, elevOrGap)
-    })
-  })
+  for (maxCov in seq_along(maxCoverages)){
+    bestMatchInfo <- fullElevGapShrink(minReadCov, windowSize, maxCoverages[[maxCov]],
+                                       elevLength, nonElev, bestMatchInfo, pileupSubset, minSize, elevOrGap)
+    for (minCov in seq_along(minCoverages)) {
+      bestMatchInfo <- fullElevGapShrink(minCoverages[[minCov]], windowSize,
+                                         maxCoverages[[maxCov]], elevLength, nonElev,
+                                         bestMatchInfo, pileupSubset, minSize, elevOrGap)
+    }
+  }
   return(bestMatchInfo)
 }
 

R/GPsInElevGaps.R

@@ -17,13 +17,12 @@
 #' @importFrom dplyr bind_rows
 #' @keywords internal
 GPsInElevGaps <- function(elevGapSummList, windowSize, gffTSV, mode, chunkContigs) {
-  GPlist <- list()
   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))
     gffTSV$geneproduct <- product
   }
-  lapply(seq_along(elevGapSummList), function(i) {
+  GPlist <- lapply(seq_along(elevGapSummList), function(i) {
     trueRefName <- elevGapSummList[[i]][[8]]
     if(mode == "metagenome"){
       refName <- elevGapSummList[[i]][[8]]
@@ -41,8 +40,9 @@ GPsInElevGaps <- function(elevGapSummList, windowSize, gffTSV, mode, chunkContig
       return(NULL)
     }
     GPs$Classification <- elevGapSummList[[i]][[7]]
-    GPlist[[i]] <<- GPs
+    GPs
   })
+  GPlist <- (GPlist[!vapply(GPlist, is.null, logical(1))])
   GPsummTable <- bind_rows(GPlist)
   return(as.data.frame(GPsummTable))
 }

R/ProActive.R

@@ -28,34 +28,27 @@
 #' on. Default is 25000.
 #' @param chunkSize If `mode`="genome" OR if `mode`="metagenome" and `chunkContigs`=TRUE,
 #' chunk the genome or contigs, respectively, into smaller subsets for pattern-matching.
-#' `chunkSize` determines the size (in bp) of each 'chunk'. Default is 50000.
+#' `chunkSize` determines the size (in bp) of each 'chunk'. Default is 100000.
 #' @param IncludeNoPatterns TRUE or FALSE, If TRUE the noPattern pattern-matches will
 #' be included in the ProActive PatternMatches output list. If you would like to visualize
 #' the noPattern pattern-matches in `plotProActiveResults()`, this should be set to TRUE.
+#' @param verbose TRUE or FALSE. Print progress messages to console. Default is TRUE.
 #' @param saveFilesTo Optional, Provide a path to the directory you wish to save
 #' output to. A folder will be made within the provided directory to store
 #' results.
 #' @importFrom utils capture.output write.table
 #' @return A list containing 6 objects described in the function description.
 #' @export
 #' @examples
-#' ## Metagenome mode with gffTSV
 #' metagenome_results <- ProActive(
 #'   pileup = sampleMetagenomePileup,
 #'   mode = "metagenome",
 #'   gffTSV = sampleMetagenomegffTSV
 #' )
-#'
-#' ## Genome mode without gffTSV
-#' genome_results <- ProActive(
-#'   pileup = exampleGenomePileupSubset,
-#'   mode = "genome"
-#' )
-#'
-#' ##gffTSV is optional!
 ProActive <- function(pileup, mode, gffTSV, windowSize = 1000, chunkContigs = FALSE,
                       minSize = 10000, maxSize = Inf, minContigLength = 30000,
-                      chunkSize = 100000, IncludeNoPatterns = FALSE, saveFilesTo) {
+                      chunkSize = 100000, IncludeNoPatterns = FALSE, verbose = TRUE,
+                      saveFilesTo) {
   ## error catching
   if ((chunkSize %% 100) > 0) {
     stop("chunkSize must be divisible by 100")
@@ -73,55 +66,58 @@ ProActive <- function(pileup, mode, gffTSV, windowSize = 1000, chunkContigs = FA
     stop("Pileup file MUST have a windowSize/binsize of 100!")
   }
   startTime <- Sys.time()
-  message("Preparing input file for pattern-matching...")
+  if(verbose){message("Preparing input file for pattern-matching...")}
   pileup <- pileupFormatter(pileup, mode)
   if (mode == "genome") {
     pileup <- genomeChunks(pileup, chunkSize)
   }
-  if (chunkContigs == TRUE) {
+  if (mode == "metagenome" & chunkContigs == TRUE) {
     pileup <- contigChunks(pileup, chunkSize)
   }
-  message("Starting pattern-matching...")
-  patternMatchSummary <- patternMatcher(pileup, windowSize, minSize, maxSize, mode, minContigLength)
-  if (IncludeNoPatterns == TRUE) {
+  if(verbose){message("Starting pattern-matching...")}
+  patternMatchSummary <- patternMatcher(pileup, windowSize, minSize, maxSize, mode, minContigLength, verbose)
+  if (IncludeNoPatterns) {
     classifList <- patternMatchSummary[[1]]
   } else {
     classifList <- removeNoPatterns(patternMatchSummary[[1]])
   }
   filteredOutContigsDf <- patternMatchSummary[[2]]
-  message("Summarizing pattern-matching results")
+  if(verbose){message("Summarizing pattern-matching results")}
   summaryTable <- classifSumm(pileup, patternMatchSummary[[1]], windowSize, mode)
   if (missing(gffTSV) == FALSE) {
-    message("Finding gene predictions in elevated or gapped regions of read coverage...")
+    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)
   }
-  message("Finalizing output")
+  if(verbose){message("Finalizing output")}
   endTime <- Sys.time()
   duration <- difftime(endTime, startTime)
-  message("Execution time: ", round(duration[[1]], 2), units(duration))
-  message(
+  if(verbose){message("Execution time: ", round(duration[[1]], 2), units(duration))}
+  if(verbose){message(
     length(which(
       filteredOutContigsDf[, 2] == "Low read cov"
     )),
     " contigs were filtered out based on low read coverage"
-  )
-  message(
+  )}
+  if(verbose){message(
     length(which(
       filteredOutContigsDf[, 2] == "Too Short"
     )),
     " contigs were filtered out based on length (< minContigLength)"
-  )
+  )}
   arguments <- list(windowSize, mode, chunkSize, chunkContigs)
   cleanSummaryTable <- summaryTable[-which(summaryTable[,2]=="NoPattern"),]
   finalSummaryList <- list(summaryTable, cleanSummaryTable, classifList, filteredOutContigsDf, arguments)
   names(finalSummaryList) <- c("SummaryTable", "CleanSummaryTable", "PatternMatches", "FilteredOut", "Arguments")
   if (missing(gffTSV) == FALSE) {
     finalSummaryList <- c(finalSummaryList, list(GPSummTable))
-    names(finalSummaryList)[6] <- "GenePredictTable"
+    names(finalSummaryList)[6] <- "GeneAnnotTable"
   }
   table <- (table(summaryTable[, 2]))
-  message(paste0(capture.output(table), collapse = "\n"))
+  if(verbose){message(paste0(capture.output(table), collapse = "\n"))}
+  if(mode == "genome" || (mode == "metagenome" & chunkContigs == TRUE)){
+    linkChunks(classifList, pileup, windowSize, mode, verbose)
+  }
   if (missing(saveFilesTo) == FALSE) {
     ifelse(!dir.exists(paths = paste0(saveFilesTo, "\\ProActiveOutput")),
       dir.create(paste0(saveFilesTo, "\\ProActiveOutput")),
@@ -135,7 +131,7 @@ ProActive <- function(pileup, mode, gffTSV, windowSize = 1000, chunkContigs = FA
         GPSummTable,
         file = paste0(
           saveFilesTo,
-          "\\ProActiveOutput\\ProActiveGenePredictstable.csv"
+          "\\ProActiveOutput\\ProActiveGeneAnnotsTable.csv"
         ),
         sep = ",",
         row.names = FALSE

R/bestMatchListFunctions.R

@@ -11,25 +11,19 @@
 #' @param mode Either "genome" or "metagenome"
 #' @keywords internal
 classifSumm <- function(pileup, bestMatchList, windowSize, mode) {
-  refName <- rep(NA, length(bestMatchList))
-  elevRatio <- rep(NA, length(bestMatchList))
-  startPos <- rep(NA, length(bestMatchList))
-  endPos <- rep(NA, length(bestMatchList))
-  classification <- rep(NA, length(bestMatchList))
-  matchSize <- rep(NA, length(bestMatchList))
   if (length(bestMatchList) == 0) {
     stop("No pattern-matches detected")
   }
-  lapply(seq_along(bestMatchList), function(i) {
-    refName[i] <<- bestMatchList[[i]][[8]]
+  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))
+  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]]), ]
     pileupSubset <- changewindowSize(pileupSubset, windowSize, mode)
-    elevRatio[i] <<- bestMatchList[[i]][[6]]
-    classification[i] <<- bestMatchList[[i]][[7]]
-    startPos[i] <<- pileupSubset[bestMatchList[[i]][[4]], 3]
-    endPos[i] <<- pileupSubset[bestMatchList[[i]][[5]], 3]
-    matchSize[i] <<- (length(seq(pileupSubset[bestMatchList[[i]][[4]], 3], pileupSubset[bestMatchList[[i]][[5]], 3], windowSize)) - 1) * windowSize
-  })
+    (length(seq(pileupSubset[bestMatchList[[i]][[4]], 3], pileupSubset[bestMatchList[[i]][[5]], 3], windowSize)) - 1) * windowSize},
+    numeric(1))
   classifSummTable <- cbind.data.frame(refName, classification, elevRatio, startPos, endPos, matchSize)
   return(classifSummTable)
 }
@@ -43,18 +37,14 @@ classifSumm <- function(pileup, bestMatchList, windowSize, mode) {
 #' all contigs/chunks classified by `ProActive()` pattern-matching
 #' @keywords internal
 removeNoPatterns <- function(bestMatchList) {
-  newBestMatchList <- list()
-  length(bestMatchList)
-  X <- 1
-  lapply(seq_along(bestMatchList), function(i) {
+  newBestMatchList <- lapply(seq_along(bestMatchList), function(i) {
     bestMatchInfo <- bestMatchList[[i]]
     classification <- bestMatchInfo[[7]]
     if (classification == "NoPattern") {
       return(NULL)
     } else {
-      newBestMatchList[[X]] <<- bestMatchInfo
-      X <<- X + 1
+      bestMatchInfo
     }
   })
-  return(newBestMatchList)
+  return(return(newBestMatchList[!vapply(newBestMatchList, is.null, logical(1))]))
 }