Refactor data (#86)

* Update minMapQuality doc

* Remove ref genome

* remove need for chr5_distance object

* Remove need for dummy_distance object

* Remove need for dummy_distance object

* Refactor chr5 MD objects to data folder

* Add data files

* Rename read_bam to distance

* Try dynamic examples on mmappr

* Document data

* Try to fix tests

* Update MmapprParam show method

* BAM files in extdata

* Refactor to use extdata BAMs

* Add MMAPPR2data to mmappr and calculateDistance docs

* Use intermediate MDs in docs

* Install MMAPPR2data in Travis script

* Use only dummy BAM in tests

* Fix file validity

* Update tests for new data

* Index every file in BamFileList

* Remove skipDebug

* Redo MD objects

* Update tests for new data

* Remove old test BAM files, update README

* Update docs to match MMAPPR2data change

* Remove .addBamFileIndex from validity check

* Delete extdata

* Bump version

* Update docs

* Fix doc error

* Change system() to system2()

* Change back to system; system2 was failing

* Remove index warning when building param

* Try lazyload=yes

* Travis: only devel

* Try installing ExperimentHub from GitHub

* Remove debug sys.which

* Try using bioc-devel

* Try bioc-devel again

* Try bioc-devel again

* Try another travis build, without sudo

* Remove username from sample objects

* Try possibleDates() again

* Try GitHub BiocManager and sudo=required

* Try ask=false

* Add fake GmapGenomeDirectory to get examples to run

* Update docs with examples

* Add fake VEP for param example

* Create directories for examples

* See if outputMmapprData caused weird cleanEx dev.off problem

* Reinstate outputMmapprData example, comment out prePeak example, which was causing weird failure

* Update docs

* Try comment out prePeak example again (fixed)

* Uncomment prePeak, didn't make a difference

* Don't run both outputMmapprData examples in both places

* Update docs

Author: kjohnsen

.travis.yml

@@ -11,10 +11,17 @@ sudo: required
 
 bioc_packages:
   - BiocCheck
+r_github_packages:
+  - kjohnsen/MMAPPR2data
+  - Bioconductor/BiocManager
 
 before_script:
   # so NAMESPACE is up-to-date
   - R -e 'devtools::document(roclets=c("rd", "collate", "namespace"))'
+  # ensure bioc-devel is being used
+  - R -e 'BiocManager::install(version = "devel", ask=FALSE)'
+  - R -e 'BiocManager::valid()'
+  # run tests here (can't see all output from check)
   - R -e 'devtools::test()'
 
 script:

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: MMAPPR2
 Title: Mutation Mapping Analysis Pipeline for Pooled RNA-Seq
-Version: 0.98.9
+Version: 0.98.10
 Authors@R: c(
     person("Kyle", "Johnsen", email="kjohnsen@byu.edu", role=c("aut")),
     person('Nathaniel', 'Jenkins', role=c('aut')),
@@ -15,10 +15,12 @@ Depends: R (>= 3.5.0)
 License: GPL-3
 Encoding: UTF-8
 LazyData: true
-RoxygenNote: 6.0.1
+LazyLoad: yes
+RoxygenNote: 6.1.0
 Suggests: testthat,
     mockery,
-    roxygen2
+    roxygen2,
+    MMAPPR2data
 Imports: ensemblVEP (>= 1.20.0),
     gmapR,
     Rsamtools,

R/aicc_loess.R

@@ -9,9 +9,7 @@
 #'   of the \code{distance} slot list filled.
 #'   
 #' @examples 
-#' \dontrun{
-#' md <- loessFit(md)
-#' }
+#' postLoessMD <- loessFit(postCalcDistMD)
 #' @export
 loessFit <- function(mmapprData) {
     loessOptResolution <- mmapprData@param@loessOptResolution

R/candidates.R

@@ -14,7 +14,7 @@
 #'
 #' @examples
 #' \dontrun{
-#' md <- generateCandidates(md)
+#' postCandidatesMD <- generateCandidates(postPeakRefMD)
 #' }
 generateCandidates <- function(mmapprData) {
 

R/data.R

@@ -0,0 +1,65 @@
+#' MmapprData object after \code{\link{calculateDistance}} step.
+#'
+#' Produced originally from BAM files of the zy13 mutant
+#' and wild-type pools, aligned to the GRCz11 reference genome
+#' using HISAT2, as available in the \link[MMAPPR2data]{MMAPPR2data}
+#' package.
+#'
+#' @format A \code{\linkS4class{MmapprData}} object.
+#' 
+#' @source \url{https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3613585/}
+"postCalcDistMD"
+
+
+#' MmapprData object after \code{\link{loessFit}} step.
+#'
+#' Produced originally from BAM files of the zy13 mutant
+#' and wild-type pools, aligned to the GRCz11 reference genome
+#' using HISAT2, as available in the \link[MMAPPR2data]{MMAPPR2data}
+#' package.
+#'
+#' @format A \code{\linkS4class{MmapprData}} object.
+#' 
+#' @source \url{https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3613585/}
+"postLoessMD"
+
+
+#' MmapprData object after \code{\link{prePeak}} step.
+#'
+#' Produced originally from BAM files of the zy13 mutant
+#' and wild-type pools, aligned to the GRCz11 reference genome
+#' using HISAT2, as available in the \link[MMAPPR2data]{MMAPPR2data}
+#' package.
+#'
+#' @format A \code{\linkS4class{MmapprData}} object.
+#' 
+#' @source \url{https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3613585/}
+"postPrePeakMD"
+
+
+#' MmapprData object after \code{\link{peakRefinement}} step.
+#'
+#' Produced originally from BAM files of the zy13 mutant
+#' and wild-type pools, aligned to the GRCz11 reference genome
+#' using HISAT2, as available in the \link[MMAPPR2data]{MMAPPR2data}
+#' package.
+#'
+#' @format A \code{\linkS4class{MmapprData}} object.
+#' 
+#' @source \url{https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3613585/}
+"postPeakRefMD"
+
+
+#' MmapprData object after \code{\link{generateCandidates}} step.
+#'
+#' Produced originally from BAM files of the zy13 mutant
+#' and wild-type pools, aligned to the GRCz11 reference genome
+#' using HISAT2, as available in the \link[MMAPPR2data]{MMAPPR2data}
+#' package.
+#'
+#' @format A \code{\linkS4class{MmapprData}} object.
+#' 
+#' @source \url{https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3613585/}
+"postCandidatesMD"
+
+

R/read_bam.R R/distance.RR/read_bam.R renamed to R/distance.R

@@ -1,6 +1,6 @@
 #' Read BAM files and generate Euclidean distance data
 #' 
-#' Initalizes the MMAPPR2 pipeline and precedes the \code{\link{loessFit}}
+#' First step in the MMAPPR2 pipeline. Precedes the \code{\link{loessFit}}
 #' step.
 #'
 #' @param mmapprData The \code{\linkS4class{MmapprData}} object to be analyzed.
@@ -10,14 +10,22 @@
 #' @export
 #'
 #' @examples
-#' \dontrun{
-#' md <- calculateDistance(md)
+#' if (requireNamespace('MMAPPR2data', quietly = TRUE)) {
+#'     ## Ignore these lines:
+#'     MMAPPR2:::.insertFakeVEPintoPath()
+#'     genDir <- gmapR::GmapGenomeDirectory('example', create=TRUE)
+#'     
+#'     mmapprParam <- MmapprParam(refGenome = gmapR::GmapGenome("GRCz11", genDir),
+#'                                wtFiles = MMAPPR2data::zy13wtBam(),
+#'                                mutFiles = MMAPPR2data::zy13mutBam(),
+#'                                species = "danio_rerio")
+#'
+#'     md <- new('MmapprData', param = mmapprParam)
+#'     postCalcDistMD <- calculateDistance(md)
 #' }
 calculateDistance <- function(mmapprData) {
-    if (is.na(Rsamtools::index(wtFiles(param(mmapprData)))))
-        Rsamtools::indexBam(wtFiles(param(mmapprData)))
-    if (is.na(Rsamtools::index(mutFiles(param(mmapprData)))))
-        Rsamtools::indexBam(mutFiles(param(mmapprData)))
+    .indexBamFileList(wtFiles(param(mmapprData)))
+    .indexBamFileList(mutFiles(param(mmapprData)))
 
     chrList <- suppressWarnings(.getFileReadChrList(mmapprData))
 
@@ -28,6 +36,15 @@ calculateDistance <- function(mmapprData) {
 }
 
 
+.indexBamFileList <- function(bfl) {
+    for (i in seq_along(bfl)) {
+        bamFile <- bfl[[i]]
+        if (is.na(Rsamtools::index(bamFile)))
+            Rsamtools::indexBam(bamFile)
+    }
+}
+
+
 .calcDistForChr <- function(chrRange, param){
     startTime <- proc.time()
     tryCatch({
@@ -133,10 +150,12 @@ calculateDistance <- function(mmapprData) {
 
 
 .getFileReadChrList <- function(mmapprData) {suppressWarnings({
-    wtFiles <- mmapprData@param@wtFiles
-    mutFiles <- mmapprData@param@mutFiles
+    bams <- Rsamtools::BamFileList(c(mmapprData@param@wtFiles,
+                                     mmapprData@param@mutFiles))
 
-    chrRanges <- as(GenomeInfoDb::seqinfo(Rsamtools::BamFileList(c(wtFiles, mutFiles))), "GRanges")
+    bamInfo <- Rsamtools::seqinfo(bams)
+    Sys.which('hello-there')
+    chrRanges <- as(bamInfo, "GRanges")
     #cut to standard chromosomes
     chrRanges <- GenomeInfoDb::keepStandardChromosomes(chrRanges, pruning.mode='coarse')
     chrRanges <- GenomeInfoDb::dropSeqlevels(chrRanges, 'chrM', pruning.mode='coarse')

R/main.R

@@ -20,26 +20,34 @@
 #'   desired parameters.
 #'
 #' @return A \code{\linkS4class{MmapprData}} object containing results
-#'   as well as intermediate data.
+#'   and/or intermediate data.
 #' @export
 #'
 #' @examples
-#' \dontrun{
-#' mmapprParam <- MmapprParam(refGenome = GmapGenome("GRCz11"),
-#'                            wtFiles = "wild_type.sorted.bam",
-#'                            mutFiles = "mutant.sorted.bam",
-#'                            species = "danio_rerio")
-#' mmapprData <- mmappr(mmapprParam)
+#' if (requireNamespace('MMAPPR2data', quietly = TRUE)) {
+#'     ## Ignore these lines:
+#'     MMAPPR2:::.insertFakeVEPintoPath()
+#'     genDir <- gmapR::GmapGenomeDirectory('DEFAULT', create=TRUE)
 #' 
-#' ### Alternately, you can navigate the pipeline step by step.
-#' ### This may be helpful for debugging.
-#' md <- new('MmapprData', param = mmapprParam)
-#' md <- calculateDistance(md)
-#' md <- loessFit(md)
-#' md <- prePeak(md)
-#' md <- peakRefinement(md)
-#' md <- generateCandidates(md)
-#' md <- outputMmapprData(md)
+#'     # Specify parameters:
+#'     mmapprParam <- MmapprParam(refGenome = gmapR::GmapGenome("GRCz11", genDir),
+#'                                wtFiles = MMAPPR2data::zy13wtBam(),
+#'                                mutFiles = MMAPPR2data::zy13mutBam(),
+#'                                species = "danio_rerio")
+#'                                
+#'     # Run pipeline:
+#'     mmapprData <- mmappr(mmapprParam)
+#' 
+#'     ### Alternately, you can navigate the pipeline step by step.
+#'     ### This may be helpful for debugging.
+#'     md <- new('MmapprData', param = mmapprParam)
+#'     postCalcDistMD <- calculateDistance(md)
+#'     postLoessMD <- loessFit(postCalcDistMD)
+#'     postPrePeakMD <- prePeak(postLoessMD)
+#'     postPeakRefMD <- peakRefinement(postPrePeakMD)
+#'     \dontrun{postCandidatesMD <- generateCandidates(postPeakRefMD)
+#'     outputMmapprData(postCandidatesMD)
+#'     }
 #' }
 #' 
 #' @seealso \code{\link{calculateDistance}}, \code{\link{loessFit}},
@@ -131,12 +139,10 @@ mmappr <- function(mmapprParam) {
 }
 
 
-.addBamFileIndex <- function(bf) {
-    path <- bf$path
-    index <- paste0(path, ".bai")
-    if (!file.exists(index)){
-        Rsamtools::indexBam(bf)
-    }
-    return(Rsamtools::BamFile(path, index = index))
-}
-
+.insertFakeVEPintoPath <- function() {
+    unlink('/tmp/bin', recursive=TRUE)
+    dir.create('/tmp/bin')
+    file.create('/tmp/bin/vep')
+    system('chmod 777 /tmp/bin/vep')
+    Sys.setenv('PATH'=paste(Sys.getenv('PATH'), '/tmp/bin', sep=':'))
+}

R/output.R

@@ -8,9 +8,8 @@
 #' @export
 #'
 #' @examples
-#'   \dontrun{
-#'     outputMmapprData(md)
-#'   }
+#' dir.create(outputFolder(param(postCandidatesMD))) ## Ignore this line
+#' \dontrun{outputMmapprData(postCandidatesMD)}
 outputMmapprData <- function(mmapprData) {
     stopifnot(class(mmapprData) == "MmapprData")
     oF <- outputFolder(param(mmapprData))

R/peaks.R

@@ -10,9 +10,7 @@
 #' @export
 #'
 #' @examples
-#' \dontrun{
-#' md <- peakRefinement(md)
-#' }
+#' postPeakRefMD <- peakRefinement(postPrePeakMD)
 peakRefinement <- function(mmapprData){
     mmapprData@peaks <- 
         .runFunctionInParallel(mmapprData@peaks,
@@ -102,9 +100,7 @@ peakRefinement <- function(mmapprData){
 #' @export
 #'
 #' @examples
-#' \dontrun{
-#' md <- prePeak(md)
-#' }
+#' postPrePeakMD <- prePeak(postLoessMD)
 prePeak <- function(mmapprData) {
     mmapprData@peaks <- list()
 

R/s4_methods.R

@@ -73,8 +73,15 @@
 #' @export
 #'
 #' @examples
-#' \dontrun{
-#' mp <- MmapprParam(gmapR::GmapGenome('GRCz11'), 'wt.bam', 'mut.bam', 'danio_rerio')
+#' if (requireNamespace('MMAPPR2data', quietly=TRUE)) {
+#'     ## Ignore these lines:
+#'     MMAPPR2:::.insertFakeVEPintoPath()
+#'     genDir <- gmapR::GmapGenomeDirectory('example', create=TRUE)
+#'     
+#'     mmapprParam <- MmapprParam(refGenome = gmapR::GmapGenome("GRCz11", genDir),
+#'                                wtFiles = MMAPPR2data::zy13wtBam(),
+#'                                mutFiles = MMAPPR2data::zy13mutBam(),
+#'                                species = "danio_rerio")
 #' }
 MmapprParam <- function(refGenome, wtFiles, mutFiles, species, vepFlags=NULL,
                         outputFolder=NULL, distancePower=4,
@@ -145,13 +152,7 @@ MmapprParam <- function(refGenome, wtFiles, mutFiles, species, vepFlags=NULL,
         errors <- c(errors, paste0(files, " is not a BamFileList object"))
     for (i in seq_along(files)) {
         file <- files[[i]]
-        if (file.exists(file$path)) {
-            if (length(Rsamtools::index(file)) == 0) {
-                file <- .addBamFileIndex(file)
-                if (length(Rsamtools::index(file)) == 0)
-                    warning(paste0(file$path), " in wtFiles has no index file")
-            }
-        } else {
+        if (!file.exists(file$path)) {
             errors <- c(errors, paste0(file$path, " does not exist"))
         }
     }
@@ -182,7 +183,8 @@ setMethod("show", "MmapprParam", function(object) {
     .customPrint(object@vepFlags, margin)
 
     cat("Other parameters:\n")
-    slotNames <- slotNames("MmapprParam")[5:length(slotNames("MmapprParam"))]
+    slotNames <- slotNames("MmapprParam")[6:length(slotNames("MmapprParam"))]
+    slotNames <- c('species', slotNames)
     slotValues <- sapply(slotNames, function(name) slot(object, name))
     names(slotValues) <- slotNames
     print(slotValues, quote=FALSE)
@@ -255,6 +257,7 @@ setMethod("show", "MmapprData", function(object) {
 #'   peakIntervalWidth peakIntervalWidth<-
 #'   minDepth minDepth<-
 #'   minBaseQuality minBaseQuality<-
+#'   minMapQuality minMapQuality<-
 #'   loessOptResolution loessOptResolution<-
 #'   loessOptCutFactor loessOptCutFactor<-
 #'   naCutoff naCutoff<-
@@ -356,15 +359,15 @@ setMethod("wtFiles<-", "MmapprParam",
           function(obj, value) {
               obj@wtFiles <- Rsamtools::BamFileList(value)
               v <- .validFiles(obj@wtFiles)
-              if (typeof(v) == 'logical') obj else v
+              if (typeof(v) == 'logical') obj else stop(v)
           })
 #' @rdname MmapprParam-functions
 #' @export
 setMethod("mutFiles<-", "MmapprParam",
           function(obj, value) {
               obj@mutFiles <- Rsamtools::BamFileList(value)
               v <- .validFiles(obj@wtFiles)
-              if (typeof(v) == 'logical') obj else v
+              if (typeof(v) == 'logical') obj else stop(v)
           })
 #' @rdname MmapprParam-functions
 #' @export