diff --git a/NEWS.md b/NEWS.md index 98a849e..cbb4693 100644 --- a/NEWS.md +++ b/NEWS.md @@ -10,12 +10,12 @@ 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. -* Fixed error when `plotRegionsHeatmap` was given more than one region that overlapped one window. - ### BUG FIXES -* Correct the message produced by `addMedipsEnrichmentFactors` (thanks @daonslog for reporting). +* Fixed error when `plotRegionsHeatmap` was given more than one region that overlapped one window. +* Correct the message produced by `addMedipsEnrichmentFactors`. * `makeQset`, `renameSamples` and `renameQsetNames` will no longer accept sample names that are not valid column names in R without quotation. +* Correctly pass the `fragmentLength` when calling `makeQset` with the `CNVmethod = "MeCap"` option, and fix an issue with hg19 GRanges. # mesa 0.5.1 diff --git a/R/makeQset.R b/R/makeQset.R index 1e88d5c..fdcff51 100644 --- a/R/makeQset.R +++ b/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. @@ -142,6 +142,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 @@ -259,12 +265,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") { diff --git a/man/makeQset.Rd b/man/makeQset.Rd index 0e28c3d..715281a 100644 --- a/man/makeQset.Rd +++ b/man/makeQset.Rd @@ -43,21 +43,21 @@ makeQset( \item{fragmentSD}{Standard deviation of the DNA fragment lengths. Can be set by fragmentType.} -\item{CNVmethod}{Which method to use for calculation of CNV. Options include "HMMdefault" (hmmcopy with default parameters) and "None"} +\item{CNVmethod}{Which method to use for calculation of CNV. Options include "HMMdefault" (hmmcopy with default parameters) and "None".} -\item{coverageMethod}{Whether to use custom methtools method for reading coverage in (set to PairedAndR1s), rather than qsea's (qseaPaired).} +\item{coverageMethod}{Whether to use custom method for reading coverage in (set to PairedAndR1s), rather than qsea's (qseaPaired).} -\item{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).} +\item{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).} -\item{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.} +\item{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.} -\item{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.} +\item{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.} -\item{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.} +\item{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.} \item{badRegions}{A GRanges object containing regions to filter out from the result.} -\item{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.} +\item{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.} \item{hmmCopyGC}{A data frame containing GC content per bin (each with size \code{CNVwindowSize}), only for use with hmmcopy.} @@ -69,5 +69,5 @@ makeQset( A qseaSet object, containing all the information required. } \description{ -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. } diff --git a/tests/testthat/test-makeQset.R b/tests/testthat/test-makeQset.R index f74d9f0..7d56a6c 100644 --- a/tests/testthat/test-makeQset.R +++ b/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()))