Fix makeqset

NEWS.md

@@ -10,6 +10,7 @@ editor_options:
 * Converted `plotPCA` into a method for the `qsea` defined method.
 * `plotPCA` gains a `verbose` option to turn off most of the messages produced.
 * `getSampleTable` is now defined for PCA/UMAP objects.
+* Correctly pass the `fragmentLength` when calling `makeQset` with the `CNVmethod = "MeCap"` option, and fix an issue with hg19 GRanges.  
 
 # mesa 0.5.1
 

R/makeQset.R

@@ -1,6 +1,6 @@
 #' Make the initial qseaSet object from a set of samples
 #'
-#' This function makes the initial qseaSet from a sampleTable that includes the file paths for the bam files and their metadata
+#' This function makes the initial qseaSet from a sampleTable that includes the file paths for the bam files and their metadata. 
 #'
 #' @param sampleTable A data frame with the sample information to be passed to qsea.
 #' @param BSgenome A BSgenome string. See BSgenome::available.genomes() for options.
@@ -10,13 +10,13 @@
 #' @param fragmentLength Average DNA fragment length. Can be set by fragmentType.
 #' @param fragmentSD Standard deviation of the DNA fragment lengths. Can be set by fragmentType.
 #' @param CNVwindowSize What window size (in bp) to use in the calculation of CNV (default 1000000).
-#' @param CNVmethod Which method to use for calculation of CNV. Options include "HMMdefault" (hmmcopy with default parameters) and "None"
-#' @param coverageMethod Whether to use custom methtools method for reading coverage in (set to PairedAndR1s), rather than qsea's (qseaPaired).
-#' @param minMapQual Minimum MAPQ score for a read to be kept. For PairedAndR1s, will keep if either R1 or R2 meet this cutoff in a properly paired read (if MQ tags are set in the bam file with samtools fixmate).
-#' @param minInsertSize For paired reads, only keep them if they are above a minimum length. Can be used for cfDNA size selection. Applies to Input samples as well as MeCap.
-#' @param maxInsertSize For paired reads, only keep them if they are below a maximum length. Can be used for cfDNA size selection. Applies to Input samples as well as MeCap.
-#' @param properPairsOnly Whether to only keep properly paired reads, or to keep high-quality (MAPQ 30+) unpaired R1s as well. Set to TRUE for size selection.
-#' @param minReferenceLength A minimum distance on the genome to keep the read. bwa by default gives 19bp as minimum for a read, which is quite short.
+#' @param CNVmethod Which method to use for calculation of CNV. Options include "HMMdefault" (hmmcopy with default parameters) and "None".
+#' @param coverageMethod Whether to use custom method for reading coverage in (set to PairedAndR1s), rather than qsea's (qseaPaired).
+#' @param minMapQual Minimum MAPQ score for a read to be kept. For coverageMethod = "PairedAndR1s", will keep if either R1 or R2 meet this cutoff in a properly paired read (if MQ tags are set in the bam file with samtools fixmate).
+#' @param minInsertSize For paired reads, only keep them if they are above a minimum length. Only applies to coverageMethod = "PairedAndR1s", and applies to Input samples as well as MeCap.
+#' @param maxInsertSize For paired reads, only keep them if they are below a maximum length. Only applies to coverageMethod = "PairedAndR1s", and applies to Input samples as well as MeCap.
+#' @param properPairsOnly Whether to only keep properly paired reads, or to keep high-quality (MAPQ 30+) unpaired R1s as well. Set to TRUE for size selection. Only applies to coverageMethod = "PairedAndR1s", and applies to Input samples as well as MeCap.
+#' @param minReferenceLength A minimum distance on the genome to keep the read. bwa by default gives 19bp as minimum for a read, which is quite short. Only applies to coverageMethod = "PairedAndR1s", and applies to Input samples as well as MeCap.
 #' @param badRegions A GRanges object containing regions to filter out from the result.
 #' @param hmmCopyGC A data frame containing GC content per bin (each with size `CNVwindowSize`), only for use with hmmcopy.
 #' @param hmmCopyMap A data frame containing Mapability content per bin (each with size `CNVwindowSize`), only for use with hmmcopy.
@@ -132,6 +132,12 @@ makeQset <- function(sampleTable,
                                  Regions = windowsWithoutBlacklist,
                                  window_size = windowSize)
 
+  if(any(genome(getRegions(qseaSet)) != BSgenome)) {
+    regions <- qseaSet %>% getRegions() 
+    genome(regions) <- BSgenome
+    qseaSet@regions <- regions
+  }
+
   # store the extra blacklist file location
   #qseaSet@parameters$badRegions2 <- blacklistBed
 
@@ -249,12 +255,13 @@ makeQset <- function(sampleTable,
     qseaSet <- qsea::addCNV(qseaSet,
                             file_name = "file_name",
                             window_size = CNVwindowSize,
-                            paired = FALSE,
+                            fragment_length = fragmentLength,
+                            paired = TRUE,
                             parallel = parallel,
                             MeDIP = TRUE
     )
 
-    qseaSet <- addMedipsEnrichmentFactors(qseaSet, nCores = ifelse(parallel, BiocParallel::bpworkers(), 1), nonEnrich = TRUE)
+    #qseaSet <- addMedipsEnrichmentFactors(qseaSet, nCores = ifelse(parallel, BiocParallel::bpworkers(), 1), nonEnrich = TRUE)
 
 
   } else if (CNVmethod == "None") {

tests/testthat/test-makeQset.R

@@ -63,7 +63,7 @@ test_that("Making a hg19 qseaSet with qsea coverage method", {
                       windowSize = 300,
                       CNVwindowSize = 1000000,
                       fragmentType = "Sheared",
-                      CNVmethod = "None",
+                      CNVmethod = "MeCap",
                       coverageMethod = "qseaPaired",
                       minMapQual = 10,
                       minInsertSize = 70,
@@ -77,6 +77,7 @@ test_that("Making a hg19 qseaSet with qsea coverage method", {
   expect_equal(testSet %>% qsea::getRegions() %>% width() %>% unique(), 300)
   expect_equal(testSet %>% qsea::getCounts() %>% colSums() %>% unname(), testSet %>% qsea::getLibSize())
   expect_equal(testSet %>% qsea::getRegions() %>% length(), 171015)
+  expect_equal(testSet %>% qsea::getCNV() %>% length(), 51)
 
   expect_true("relH" %in% (testSet %>% addLibraryInformation() %>% qsea::getSampleTable() %>% colnames()))