diff --git a/R/all_tools.R b/R/all_tools.R index d194ebf..91c738c 100755 --- a/R/all_tools.R +++ b/R/all_tools.R @@ -72,11 +72,11 @@ extractPLAF <- function(plafFileName) { #' refFile <- system.file("extdata", "PG0390-C.test.ref", package = "DEploid.utils") #' altFile <- system.file("extdata", "PG0390-C.test.alt", package = "DEploid.utils") #' PG0390CoverageTxt <- extractCoverageFromTxt(refFile, altFile) -#' #PG0390CoverageTxt.deconv <- dEploid(paste( +#' # PG0390CoverageTxt.deconv <- dEploid(paste( #' # "-ref", refFile, "-alt", altFile, #' # "-plaf", plafFile, "-noPanel" -#' #)) -#' #plotProportions(PG0390CoverageTxt.deconv$Proportions, "PG0390-C proportions") +#' # )) +#' # plotProportions(PG0390CoverageTxt.deconv$Proportions, "PG0390-C proportions") #' plotProportions <- function(proportions, title = "Components", cex.lab = 1, cex.main = 1, cex.axis = 1) { @@ -189,7 +189,7 @@ plotAltVsRef <- function(ref, alt, title = "Alt vs Ref", #' #' # Example 2 #' vcfFile <- system.file("extdata", "PG0390-C.test.vcf.gz", package = "DEploid.utils") -#' PG0390CoverageVcf <- extractCoverageFromVcf(vcfFile) +#' PG0390CoverageVcf <- extractCoverageFromVcf(vcfFile, "PG0390-C") #' obsWSAF <- computeObsWSAF(PG0390CoverageVcf$altCount, PG0390CoverageVcf$refCount) #' histWSAF(obsWSAF) #' myhist <- histWSAF(obsWSAF, FALSE) @@ -244,7 +244,7 @@ histWSAF <- function(obsWSAF, exclusive = TRUE, #' #' # Example 2 #' vcfFile <- system.file("extdata", "PG0390-C.test.vcf.gz", package = "DEploid.utils") -#' PG0390CoverageVcf <- extractCoverageFromVcf(vcfFile) +#' PG0390CoverageVcf <- extractCoverageFromVcf(vcfFile, "PG0390-C") #' obsWSAF <- computeObsWSAF(PG0390CoverageVcf$altCount, PG0390CoverageVcf$refCount) #' plafFile <- system.file("extdata", "labStrains.test.PLAF.txt", package = "DEploid.utils") #' plaf <- extractPLAF(plafFile) @@ -328,7 +328,7 @@ plotObsExpWSAF <- function(obsWSAF, expWSAF, #' #' # Example 2 #' vcfFile <- system.file("extdata", "PG0390-C.test.vcf.gz", package = "DEploid.utils") -#' PG0390CoverageVcf <- extractCoverageFromVcf(vcfFile) +#' PG0390CoverageVcf <- extractCoverageFromVcf(vcfFile, "PG0390-C") #' obsWSAF <- computeObsWSAF(PG0390CoverageVcf$altCount, PG0390CoverageVcf$refCount) #' computeObsWSAF <- function(alt, ref) { diff --git a/man/computeObsWSAF.Rd b/man/computeObsWSAF.Rd index 3bd83d0..5d4ed8e 100644 --- a/man/computeObsWSAF.Rd +++ b/man/computeObsWSAF.Rd @@ -26,7 +26,7 @@ obsWSAF <- computeObsWSAF(PG0390CoverageTxt$altCount, PG0390CoverageTxt$refCount # Example 2 vcfFile <- system.file("extdata", "PG0390-C.test.vcf.gz", package = "DEploid.utils") -PG0390CoverageVcf <- extractCoverageFromVcf(vcfFile) +PG0390CoverageVcf <- extractCoverageFromVcf(vcfFile, "PG0390-C") obsWSAF <- computeObsWSAF(PG0390CoverageVcf$altCount, PG0390CoverageVcf$refCount) } diff --git a/man/histWSAF.Rd b/man/histWSAF.Rd index 88c2976..2b96784 100644 --- a/man/histWSAF.Rd +++ b/man/histWSAF.Rd @@ -43,7 +43,7 @@ myhist <- histWSAF(obsWSAF, FALSE) # Example 2 vcfFile <- system.file("extdata", "PG0390-C.test.vcf.gz", package = "DEploid.utils") -PG0390CoverageVcf <- extractCoverageFromVcf(vcfFile) +PG0390CoverageVcf <- extractCoverageFromVcf(vcfFile, "PG0390-C") obsWSAF <- computeObsWSAF(PG0390CoverageVcf$altCount, PG0390CoverageVcf$refCount) histWSAF(obsWSAF) myhist <- histWSAF(obsWSAF, FALSE) diff --git a/man/plotProportions.Rd b/man/plotProportions.Rd index 79d0b51..bcd4cf6 100644 --- a/man/plotProportions.Rd +++ b/man/plotProportions.Rd @@ -32,10 +32,10 @@ panelFile <- system.file("extdata", "labStrains.test.panel.txt", package = "DEpl refFile <- system.file("extdata", "PG0390-C.test.ref", package = "DEploid.utils") altFile <- system.file("extdata", "PG0390-C.test.alt", package = "DEploid.utils") PG0390CoverageTxt <- extractCoverageFromTxt(refFile, altFile) -#PG0390CoverageTxt.deconv <- dEploid(paste( +# PG0390CoverageTxt.deconv <- dEploid(paste( # "-ref", refFile, "-alt", altFile, # "-plaf", plafFile, "-noPanel" -#)) -#plotProportions(PG0390CoverageTxt.deconv$Proportions, "PG0390-C proportions") +# )) +# plotProportions(PG0390CoverageTxt.deconv$Proportions, "PG0390-C proportions") } diff --git a/man/plotWSAFvsPLAF.Rd b/man/plotWSAFvsPLAF.Rd index a3b0c5f..2da144b 100644 --- a/man/plotWSAFvsPLAF.Rd +++ b/man/plotWSAFvsPLAF.Rd @@ -47,7 +47,7 @@ plotWSAFvsPLAF(plaf, obsWSAF) # Example 2 vcfFile <- system.file("extdata", "PG0390-C.test.vcf.gz", package = "DEploid.utils") -PG0390CoverageVcf <- extractCoverageFromVcf(vcfFile) +PG0390CoverageVcf <- extractCoverageFromVcf(vcfFile, "PG0390-C") obsWSAF <- computeObsWSAF(PG0390CoverageVcf$altCount, PG0390CoverageVcf$refCount) plafFile <- system.file("extdata", "labStrains.test.PLAF.txt", package = "DEploid.utils") plaf <- extractPLAF(plafFile) diff --git a/tests/testthat/test-extractFromVCF.R b/tests/testthat/test-extractFromVCF.R index 4114014..4a93b9b 100644 --- a/tests/testthat/test-extractFromVCF.R +++ b/tests/testthat/test-extractFromVCF.R @@ -1,4 +1,3 @@ - #' @title Extract read counts from VCF #' #' @description Extract read counts from VCF file of a single sample. @@ -15,14 +14,15 @@ #' allele count, alternative allele count. #' #' @examples -#' vcfFile = system.file("extdata", "PG0390-C.test.vcf.gz", package = "DEploid") -#' PG0390 = extractCoverageFromVcf(vcfFile) +#' vcfFile <- system.file("extdata", "PG0390-C.test.vcf.gz", package = "DEploid") +#' PG0390 <- extractCoverageFromVcf(vcfFile) #' old_extractCoverageFromVcf <- function(vcfFileName, ADFieldIndex = 2) { # Assume that AD is the second field h <- function(w) { - if (any(grepl("gzfile connection", w))) + if (any(grepl("gzfile connection", w))) { invokeRestart("muffleWarning") + } } gzf <- gzfile(vcfFileName, open = "rb") @@ -38,9 +38,11 @@ old_extractCoverageFromVcf <- function(vcfFileName, ADFieldIndex = 2) { } close(gzf) - vcf <- read.table(gzfile(vcfFileName), skip = numberOfHeaderLines, - header = T, comment.char = "", stringsAsFactors = FALSE, - check.names = FALSE) + vcf <- read.table(gzfile(vcfFileName), + skip = numberOfHeaderLines, + header = T, comment.char = "", stringsAsFactors = FALSE, + check.names = FALSE + ) sampleName <- names(vcf)[10] @@ -53,22 +55,20 @@ old_extractCoverageFromVcf <- function(vcfFileName, ADFieldIndex = 2) { refCount <- as.numeric(unlist(lapply(tmpCov, `[[`, 1))) altCount <- as.numeric(unlist(lapply(tmpCov, `[[`, 2))) - return(data.frame(CHROM = vcf[, 1], - POS = vcf[, 2], - refCount = refCount, - altCount = altCount)) + return(data.frame( + CHROM = vcf[, 1], + POS = vcf[, 2], + refCount = refCount, + altCount = altCount + )) } -test_that("test R vs cpp vcf extract", - { - vcfFile = system.file("extdata", "PG0390-C.test.vcf.gz", package = "DEploid.utils") - vcf = extractCoverageFromVcf(vcfFile, "PG0390-C") - PG0390 = old_extractCoverageFromVcf(vcfFile) - testthat::expect_equal(names(vcf), names(PG0390)) - testthat::expect_equal(vcf$CHROM, PG0390$CHROM) - testthat::expect_equal(vcf$POS, PG0390$POS) - testthat::expect_equal(vcf$refCount, PG0390$refCount) - testthat::expect_equal(vcf$altCount, PG0390$altCount) - } -) +test_that("test R vs cpp vcf extract", { + PG0390 <- old_extractCoverageFromVcf(vcfFileName) + testthat::expect_equal(names(PG0390CoverageVcf), names(PG0390)) + testthat::expect_equal(PG0390CoverageVcf$CHROM, PG0390$CHROM) + testthat::expect_equal(PG0390CoverageVcf$POS, PG0390$POS) + testthat::expect_equal(PG0390CoverageVcf$refCount, PG0390$refCount) + testthat::expect_equal(PG0390CoverageVcf$altCount, PG0390$altCount) +}) diff --git a/tests/testthat/test-vcfreader.R b/tests/testthat/test-vcfreader.R index 15c3cf7..e8ef65f 100644 --- a/tests/testthat/test-vcfreader.R +++ b/tests/testthat/test-vcfreader.R @@ -1,16 +1,3 @@ -vcfFileName <- system.file("extdata", "PG0390-C.test.vcf.gz", - package = "DEploid.utils") -plafFileName <- system.file("extdata", "labStrains.test.PLAF.txt", - package = "DEploid.utils") -panelFileName <- system.file("extdata", "labStrains.test.panel.txt", - package = "DEploid.utils") -refFileName <- system.file("extdata", "PG0390-C.test.ref", package = "DEploid.utils") -altFileName <- system.file("extdata", "PG0390-C.test.alt", package = "DEploid.utils") - -PG0390CoverageVcf <- extractCoverageFromVcf(vcfFileName, "PG0390-C") -plaf <- extractPLAF(plafFileName) - - test_that("Extracted coverage", { PG0390CoverageTxt <- extractCoverageFromTxt(refFileName, altFileName) expect_that(PG0390CoverageTxt, is_a("data.frame")) @@ -35,8 +22,10 @@ test_that("computeObsWSAF", { test_that("WSAF Related", { - obsWSAF <- computeObsWSAF(PG0390CoverageVcf$altCount, - PG0390CoverageVcf$refCount) + obsWSAF <- computeObsWSAF( + PG0390CoverageVcf$altCount, + PG0390CoverageVcf$refCount + ) potentialOutliers <- c(5, 12, 25, 30, 35, 50) expect_that(histWSAF(obsWSAF), is_a("histogram")) @@ -53,8 +42,10 @@ test_that("WSAF Related", { test_that("plotAltVsRef", { - expect_null(plotAltVsRef(PG0390CoverageVcf$refCount, - PG0390CoverageVcf$altCount)) + expect_null(plotAltVsRef( + PG0390CoverageVcf$refCount, + PG0390CoverageVcf$altCount + )) png("AltVsRef.png") plotAltVsRef(PG0390CoverageVcf$refCount, PG0390CoverageVcf$altCount) dev.off() @@ -64,11 +55,12 @@ test_that("plotAltVsRef", { test_that("plotAltVsRefWithOutliers", { potentialOutliers <- c(1, 10, 20, 30, 40) expect_null(plotAltVsRef(PG0390CoverageVcf$refCount, - PG0390CoverageVcf$altCount, - potentialOutliers = potentialOutliers)) + PG0390CoverageVcf$altCount, + potentialOutliers = potentialOutliers + )) png("AltVsRefOutlier.png") plotAltVsRef(PG0390CoverageVcf$refCount, PG0390CoverageVcf$altCount, - potentialOutliers = potentialOutliers) + potentialOutliers = potentialOutliers + ) dev.off() }) -