tests/testthat/test_digestFastqs_cis.R
In fmicompbio/mutscan: Preprocessing and Analysis of Deep Mutational Scanning Data
test_that("digestFastqs works as expected for cis experiments", {
fqc1 <- system.file("extdata/cisInput_1.fastq.gz", package = "mutscan")
fqc2 <- system.file("extdata/cisInput_2.fastq.gz", package = "mutscan")
## default arguments
Ldef <- list(
fastqForward = fqc1, fastqReverse = fqc2,
mergeForwardReverse = TRUE,
minOverlap = 0, maxOverlap = 0, maxFracMismatchOverlap = 1, greedyOverlap = TRUE,
revComplForward = FALSE, revComplReverse = TRUE,
elementsForward = "SUCV", elementsReverse = "SUCVS",
elementLengthsForward = c(1, 10, 18, 96),
elementLengthsReverse = c(1, 7, 17, 96, -1),
adapterForward = "GGAAGAGCACACGTC",
adapterReverse = "GGAAGAGCGTCGTGT",
primerForward = "",
primerReverse = "",
wildTypeForward = "ACTGATACACTCCAAGCGGAGACAGACCAACTAGAAGATGAGAAGTCTGCTTTGCAGACCGAGATTGCCAACCTGCTGAAGGAGAAGGAAAAACTA",
wildTypeReverse = "",
constantForward = "AACCGGAGGAGGGAGCTG",
constantReverse = "GAGTTCATCCTGGCAGC",
avePhredMinForward = 20.0, avePhredMinReverse = 20.0,
variableNMaxForward = 0, variableNMaxReverse = 0,
umiNMax = 0,
nbrMutatedCodonsMaxForward = 1,
nbrMutatedCodonsMaxReverse = 1,
nbrMutatedBasesMaxForward = -1,
nbrMutatedBasesMaxReverse = -1,
forbiddenMutatedCodonsForward = "NNW",
forbiddenMutatedCodonsReverse = "NNW",
useTreeWTmatch = FALSE,
mutatedPhredMinForward = 0.0, mutatedPhredMinReverse = 0.0,
mutNameDelimiter = ".",
constantMaxDistForward = -1,
constantMaxDistReverse = -1,
umiCollapseMaxDist = 0,
filteredReadsFastqForward = "",
filteredReadsFastqReverse = "",
maxNReads = -1, verbose = FALSE,
nThreads = 1, chunkSize = 1000,
maxReadLength = 1024
)
res <- do.call(digestFastqs, Ldef)
## Specify reverse WT - check that it's ignored
LdefrevWT <- Ldef
LdefrevWT$wildTypeReverse <- "ACGT"
expect_warning(resrevWT <- do.call(digestFastqs, LdefrevWT),
"Ignoring 'wildTypeReverse'")
expect_equal(res$filterSummary, resrevWT$filterSummary)
expect_equal(res$summaryTable, resrevWT$summaryTable)
expect_equal(res$filterSummary$nbrTotal, 1000L)
expect_equal(res$filterSummary$f1_nbrAdapter, 126L)
expect_equal(res$filterSummary$f2_nbrNoPrimer, 0L)
expect_equal(res$filterSummary$f3_nbrReadWrongLength, 0L)
expect_equal(res$filterSummary$f4_nbrNoValidOverlap, 0L)
expect_equal(res$filterSummary$f5_nbrAvgVarQualTooLow, 0L)
expect_equal(res$filterSummary$f6_nbrTooManyNinVar, 44L)
expect_equal(res$filterSummary$f7_nbrTooManyNinUMI, 0L)
expect_equal(res$filterSummary$f8_nbrTooManyBestWTHits, 0L)
expect_equal(res$filterSummary$f9_nbrMutQualTooLow, 0L)
expect_equal(res$filterSummary$f10a_nbrTooManyMutCodons, 581L)
expect_equal(res$filterSummary$f10b_nbrTooManyMutBases, 0L)
expect_equal(res$filterSummary$f11_nbrForbiddenCodons, 82L)
expect_equal(res$filterSummary$f12_nbrTooManyMutConstant, 0L)
expect_equal(res$filterSummary$f13_nbrTooManyBestConstantHits, 0L)
expect_equal(res$filterSummary$nbrRetained, 167)
for (nm in setdiff(names(Ldef), c("forbiddenMutatedCodonsForward", "forbiddenMutatedCodonsReverse", "verbose", "fastqForward", "fastqReverse"))) {
expect_equal(res$parameters[[nm]], Ldef[[nm]], ignore_attr = TRUE)
}
for (nm in c("fastqForward", "fastqReverse")) {
expect_equal(res$parameters[[nm]], normalizePath(Ldef[[nm]], mustWork = FALSE),
ignore_attr = TRUE)
}
expect_equal(sum(res$summaryTable$nbrReads), res$filterSummary$nbrRetained)
expect_equal(sum(res$summaryTable$nbrReads == 2), 11L)
expect_equal(sort(res$summaryTable$mutantName[res$summaryTable$nbrReads == 2]),
sort(c("f.14.AAG", "f.15.ATG", "f.19.CAC", "f.1.ACC", "f.20.AAC", "f.21.GTG",
"f.24.AGC", "f.4.CGC", "f.7.GTG", "f.9.GGC", "f.9.GTC")))
expect_true(all(res$summaryTable$nbrReads == res$summaryTable$nbrUmis))
## Check the number of reads with a given number of mismatches
expect_equal(sum(res$summaryTable$nbrReads[res$summaryTable$nbrMutBases == 0]), 77L)
expect_equal(sum(res$summaryTable$nbrReads[res$summaryTable$nbrMutBases == 1]), 88L)
expect_equal(sum(res$summaryTable$nbrReads[res$summaryTable$nbrMutBases == 2]), 2L)
expect_equal(sum(res$summaryTable$nbrReads[res$summaryTable$nbrMutCodons == 0]), 77L)
expect_equal(sum(res$summaryTable$nbrReads[res$summaryTable$nbrMutCodons == 1]), 90L)
expect_equal(sum(res$summaryTable$nbrReads[res$summaryTable$nbrMutAAs == 0]), 97L)
expect_equal(sum(res$summaryTable$nbrReads[res$summaryTable$nbrMutAAs == 1]), 70L)
## Similarly but for the number of unique sequences
expect_equal(sum(res$summaryTable$nbrMutBases == 0), 1L)
expect_equal(sum(res$summaryTable$nbrMutBases == 1), 64L)
expect_equal(sum(res$summaryTable$nbrMutBases == 2), 2L)
expect_equal(sum(res$summaryTable$nbrMutCodons == 0), 1L)
expect_equal(sum(res$summaryTable$nbrMutCodons == 1), 66L)
expect_equal(sum(res$summaryTable$nbrMutAAs == 0), 14L)
expect_equal(sum(res$summaryTable$nbrMutAAs == 1), 53L)
## Check that mutant naming worked (compare to manual matching)
example_seq <- paste0("ACTGATACACGCCAAGCGGAGACAGACCAACTAGAAGATGAGAAGTCT",
"GCTTTGCAGACCGAGATTGCCAACCTGCTGAAGGAGAAGGAAAAACTA")
expect_equal(res$summaryTable$mutantName[res$summaryTable$sequence == example_seq],
"f.4.CGC")
expect_equal(res$summaryTable$mutantNameAA[res$summaryTable$sequence == example_seq],
"f.4.R")
expect_equal(res$summaryTable$mutantNameBase[res$summaryTable$sequence == example_seq],
"f.11.G")
expect_equal(res$summaryTable$mutantNameCodon[res$summaryTable$sequence == example_seq],
"f.4.CGC")
expect_equal(res$summaryTable$mutantNameBaseHGVS[res$summaryTable$sequence == example_seq],
"f:c.11T>G")
expect_equal(res$summaryTable$mutantNameAAHGVS[res$summaryTable$sequence == example_seq],
"f:p.(Leu4Arg)")
expect_equal(res$summaryTable$sequenceAA,
as.character(Biostrings::translate(
Biostrings::DNAStringSet(res$summaryTable$sequence))))
expect_equal(as.numeric(res$summaryTable$varLengths),
nchar(res$summaryTable$sequence))
expect_equal(sum(grepl("WT", res$summaryTable$mutantNameAA) & res$summaryTable$nbrMutBases > 0), 13)
expect_true(all(grepl("silent", res$summaryTable$mutationTypes[grepl("WT", res$summaryTable$mutantNameAA) & res$summaryTable$nbrMutBases > 0])))
expect_true(all(grepl("nonsynonymous|stop", res$summaryTable$mutationTypes[grepl("\\.[1-9][0-9]*\\.", res$summaryTable$mutantNameAA)])))
expect_true(all(grepl("stop", res$summaryTable$mutationTypes[grepl("TAG|TAA|TGA", res$summaryTable$mutantName)])))
expect_equal(sum(res$errorStatistics$nbrMatchForward + res$errorStatistics$nbrMismatchForward),
nchar(Ldef$constantForward[1]) * res$filterSummary$nbrRetained)
expect_equal(sum(res$errorStatistics$nbrMatchReverse + res$errorStatistics$nbrMismatchReverse),
nchar(Ldef$constantReverse[1]) * res$filterSummary$nbrRetained)
expect_equal(res$errorStatistics$nbrMatchForward[res$errorStatistics$PhredQuality == 14], 51L)
expect_equal(res$errorStatistics$nbrMatchForward[res$errorStatistics$PhredQuality == 22], 19L)
expect_equal(res$errorStatistics$nbrMatchForward[res$errorStatistics$PhredQuality == 27], 79L)
expect_equal(res$errorStatistics$nbrMatchForward[res$errorStatistics$PhredQuality == 33], 163L)
expect_equal(res$errorStatistics$nbrMatchForward[res$errorStatistics$PhredQuality == 37], 2690L)
expect_equal(res$errorStatistics$nbrMismatchForward[res$errorStatistics$PhredQuality == 14], 1L)
expect_equal(res$errorStatistics$nbrMismatchForward[res$errorStatistics$PhredQuality == 27], 1L)
expect_equal(res$errorStatistics$nbrMismatchForward[res$errorStatistics$PhredQuality == 37], 2L)
expect_equal(res$errorStatistics$nbrMatchReverse[res$errorStatistics$PhredQuality == 14], 41L)
expect_equal(res$errorStatistics$nbrMatchReverse[res$errorStatistics$PhredQuality == 27], 56L)
expect_equal(res$errorStatistics$nbrMatchReverse[res$errorStatistics$PhredQuality == 33], 96L)
expect_equal(res$errorStatistics$nbrMatchReverse[res$errorStatistics$PhredQuality == 37], 2620L)
expect_equal(res$errorStatistics$nbrMismatchReverse[res$errorStatistics$PhredQuality == 2], 4L)
expect_equal(res$errorStatistics$nbrMismatchReverse[res$errorStatistics$PhredQuality == 14], 7L)
expect_equal(res$errorStatistics$nbrMismatchReverse[res$errorStatistics$PhredQuality == 33], 1L)
expect_equal(res$errorStatistics$nbrMismatchReverse[res$errorStatistics$PhredQuality == 37], 14L)
## Add an additional constant sequence of the wrong length - should fail
Ldefcs <- Ldef
Ldefcs$constantForward <- c("AACCGGAGGAGGGAGCTG", "AACCGGAGGAGGGAGCT")
expect_error(do.call(digestFastqs, Ldefcs),
"All constant sequences must be of the same length")
})
test_that("digestFastqs works as expected when specifying max nbr of mutated bases rather than codons", {
fqc1 <- system.file("extdata/cisInput_1.fastq.gz", package = "mutscan")
fqc2 <- system.file("extdata/cisInput_2.fastq.gz", package = "mutscan")
## default arguments
Ldef <- list(
fastqForward = fqc1, fastqReverse = fqc2,
mergeForwardReverse = TRUE,
minOverlap = 0, maxOverlap = 0, maxFracMismatchOverlap = 1, greedyOverlap = TRUE,
revComplForward = FALSE, revComplReverse = TRUE,
elementsForward = "SUCV", elementsReverse = "SUCVS",
elementLengthsForward = c(1, 10, 18, 96),
elementLengthsReverse = c(1, 7, 17, 96, -1),
adapterForward = "GGAAGAGCACACGTC",
adapterReverse = "GGAAGAGCGTCGTGT",
primerForward = "",
primerReverse = "",
wildTypeForward = "ACTGATACACTCCAAGCGGAGACAGACCAACTAGAAGATGAGAAGTCTGCTTTGCAGACCGAGATTGCCAACCTGCTGAAGGAGAAGGAAAAACTA",
wildTypeReverse = "",
constantForward = "AACCGGAGGAGGGAGCTG",
constantReverse = "GAGTTCATCCTGGCAGC",
avePhredMinForward = 20.0, avePhredMinReverse = 20.0,
variableNMaxForward = 0, variableNMaxReverse = 0,
umiNMax = 0,
nbrMutatedCodonsMaxForward = -1,
nbrMutatedCodonsMaxReverse = -1,
nbrMutatedBasesMaxForward = 2,
nbrMutatedBasesMaxReverse = 2,
forbiddenMutatedCodonsForward = "NNW",
forbiddenMutatedCodonsReverse = "NNW",
useTreeWTmatch = FALSE,
mutatedPhredMinForward = 0.0, mutatedPhredMinReverse = 0.0,
mutNameDelimiter = ".",
constantMaxDistForward = -1,
constantMaxDistReverse = -1,
umiCollapseMaxDist = 0,
filteredReadsFastqForward = "",
filteredReadsFastqReverse = "",
maxNReads = -1, verbose = FALSE,
nThreads = 1, chunkSize = 1000,
maxReadLength = 125
)
res <- do.call(digestFastqs, Ldef)
expect_equal(res$filterSummary$nbrTotal, 1000L)
expect_equal(res$filterSummary$f1_nbrAdapter, 126L)
expect_equal(res$filterSummary$f2_nbrNoPrimer, 0L)
expect_equal(res$filterSummary$f3_nbrReadWrongLength, 0L)
expect_equal(res$filterSummary$f4_nbrNoValidOverlap, 0L)
expect_equal(res$filterSummary$f5_nbrAvgVarQualTooLow, 0L)
expect_equal(res$filterSummary$f6_nbrTooManyNinVar, 44L)
expect_equal(res$filterSummary$f7_nbrTooManyNinUMI, 0L)
expect_equal(res$filterSummary$f8_nbrTooManyBestWTHits, 0L)
expect_equal(res$filterSummary$f9_nbrMutQualTooLow, 0L)
expect_equal(res$filterSummary$f10a_nbrTooManyMutCodons, 0L)
expect_equal(res$filterSummary$f10b_nbrTooManyMutBases, 416L)
expect_equal(res$filterSummary$f11_nbrForbiddenCodons, 212L)
expect_equal(res$filterSummary$f12_nbrTooManyMutConstant, 0L)
expect_equal(res$filterSummary$f13_nbrTooManyBestConstantHits, 0L)
expect_equal(res$filterSummary$nbrRetained, 202L)
expect_named(res$summaryTable, c("mutantName", "sequence","nbrReads",
"maxNbrReads", "nbrUmis", "nbrMutBases",
"nbrMutCodons", "nbrMutAAs", "varLengths",
"mutantNameBase", "mutantNameCodon",
"mutantNameBaseHGVS", "mutantNameAA",
"mutantNameAAHGVS", "mutationTypes",
"sequenceAA"))
for (nm in setdiff(names(Ldef), c("forbiddenMutatedCodonsForward", "forbiddenMutatedCodonsReverse", "verbose", "fastqForward", "fastqReverse"))) {
expect_equal(res$parameters[[nm]], Ldef[[nm]], ignore_attr = TRUE)
}
for (nm in c("fastqForward", "fastqReverse")) {
expect_equal(res$parameters[[nm]], normalizePath(Ldef[[nm]], mustWork = FALSE),
ignore_attr = TRUE)
}
expect_equal(sum(res$summaryTable$nbrReads), res$filterSummary$nbrRetained)
expect_equal(sum(res$summaryTable$nbrReads == 3), 5L)
expect_equal(sort(res$summaryTable$mutantName[res$summaryTable$nbrReads == 3]),
sort(c("f.12.G", "f.17.T", "f.60.G", "f.85.T", "f.96.G")))
expect_true(all(res$summaryTable$nbrReads == res$summaryTable$nbrUmis))
## Check the number of reads with a given number of mismatches
expect_equal(sum(res$summaryTable$nbrReads[res$summaryTable$nbrMutBases == 0]), 77L)
expect_equal(sum(res$summaryTable$nbrReads[res$summaryTable$nbrMutBases == 1]), 88L)
expect_equal(sum(res$summaryTable$nbrReads[res$summaryTable$nbrMutBases == 2]), 37L)
expect_equal(sum(res$summaryTable$nbrReads[res$summaryTable$nbrMutCodons == 0]), 77L)
expect_equal(sum(res$summaryTable$nbrReads[res$summaryTable$nbrMutCodons == 1]), 90L)
expect_equal(sum(res$summaryTable$nbrReads[res$summaryTable$nbrMutCodons == 2]), 35L)
expect_equal(sum(res$summaryTable$nbrReads[res$summaryTable$nbrMutAAs == 0]), 98L)
expect_equal(sum(res$summaryTable$nbrReads[res$summaryTable$nbrMutAAs == 1]), 86L)
expect_equal(sum(res$summaryTable$nbrReads[res$summaryTable$nbrMutAAs == 2]), 18L)
## Similarly but for the number of unique sequences
expect_equal(sum(res$summaryTable$nbrMutBases == 0), 1L)
expect_equal(sum(res$summaryTable$nbrMutBases == 1), 64L)
expect_equal(sum(res$summaryTable$nbrMutBases == 2), 36L)
expect_equal(sum(res$summaryTable$nbrMutCodons == 0), 1L)
expect_equal(sum(res$summaryTable$nbrMutCodons == 1), 66L)
expect_equal(sum(res$summaryTable$nbrMutCodons == 2), 34L)
expect_equal(sum(res$summaryTable$nbrMutAAs == 0), 15L)
expect_equal(sum(res$summaryTable$nbrMutAAs == 1), 69L)
expect_equal(sum(res$summaryTable$nbrMutAAs == 2), 17L)
## Check that mutant naming worked (compare to manual matching)
example_seq <- paste0("ACTGATACACGCCAAGCGGAGACAGACCAACTAGAAGATGAGAAGTCT",
"GCTTTGCAGACCGAGATTGCCAACCTGCTGAAGGAGAAGGAAAAACTA")
expect_equal(res$summaryTable$mutantName[res$summaryTable$sequence == example_seq],
"f.11.G")
expect_equal(res$summaryTable$mutantNameAA[res$summaryTable$sequence == example_seq],
"f.4.R")
expect_equal(res$summaryTable$mutantNameBase[res$summaryTable$sequence == example_seq],
"f.11.G")
expect_equal(res$summaryTable$mutantNameCodon[res$summaryTable$sequence == example_seq],
"f.4.CGC")
expect_equal(res$summaryTable$mutantNameBaseHGVS[res$summaryTable$sequence == example_seq],
"f:c.11T>G")
expect_equal(res$summaryTable$mutantNameAAHGVS[res$summaryTable$sequence == example_seq],
"f:p.(Leu4Arg)")
expect_equal(sum(res$errorStatistics$nbrMatchForward + res$errorStatistics$nbrMismatchForward),
nchar(Ldef$constantForward[1]) * res$filterSummary$nbrRetained)
expect_equal(sum(res$errorStatistics$nbrMatchReverse + res$errorStatistics$nbrMismatchReverse),
nchar(Ldef$constantReverse[1]) * res$filterSummary$nbrRetained)
expect_equal(res$errorStatistics$nbrMatchForward[res$errorStatistics$PhredQuality == 14], 61L)
expect_equal(res$errorStatistics$nbrMatchForward[res$errorStatistics$PhredQuality == 22], 22L)
expect_equal(res$errorStatistics$nbrMatchForward[res$errorStatistics$PhredQuality == 27], 97L)
expect_equal(res$errorStatistics$nbrMatchForward[res$errorStatistics$PhredQuality == 33], 203L)
expect_equal(res$errorStatistics$nbrMatchForward[res$errorStatistics$PhredQuality == 37], 3249L)
expect_equal(res$errorStatistics$nbrMismatchForward[res$errorStatistics$PhredQuality == 14], 1L)
expect_equal(res$errorStatistics$nbrMismatchForward[res$errorStatistics$PhredQuality == 27], 1L)
expect_equal(res$errorStatistics$nbrMismatchForward[res$errorStatistics$PhredQuality == 37], 2L)
expect_equal(res$errorStatistics$nbrMatchReverse[res$errorStatistics$PhredQuality == 14], 42L)
expect_equal(res$errorStatistics$nbrMatchReverse[res$errorStatistics$PhredQuality == 27], 63L)
expect_equal(res$errorStatistics$nbrMatchReverse[res$errorStatistics$PhredQuality == 33], 117L)
expect_equal(res$errorStatistics$nbrMatchReverse[res$errorStatistics$PhredQuality == 37], 3184L)
expect_equal(res$errorStatistics$nbrMismatchReverse[res$errorStatistics$PhredQuality == 2], 5L)
expect_equal(res$errorStatistics$nbrMismatchReverse[res$errorStatistics$PhredQuality == 14], 7L)
expect_equal(res$errorStatistics$nbrMismatchReverse[res$errorStatistics$PhredQuality == 33], 2L)
expect_equal(res$errorStatistics$nbrMismatchReverse[res$errorStatistics$PhredQuality == 37], 14L)
})
test_that("digestFastqs gives the same results regardless of how WT matching is done", {
fqc1 <- system.file("extdata/cisInput_1.fastq.gz", package = "mutscan")
fqc2 <- system.file("extdata/cisInput_2.fastq.gz", package = "mutscan")
## default arguments
Ldef <- list(
fastqForward = fqc1, fastqReverse = fqc2,
mergeForwardReverse = TRUE,
minOverlap = 0, maxOverlap = 0, maxFracMismatchOverlap = 1, greedyOverlap = TRUE,
revComplForward = FALSE, revComplReverse = TRUE,
elementsForward = "SUCV", elementsReverse = "SUCVS",
elementLengthsForward = c(1, 10, 18, 96),
elementLengthsReverse = c(1, 7, 17, 96, -1),
adapterForward = "GGAAGAGCACACGTC",
adapterReverse = "GGAAGAGCGTCGTGT",
primerForward = "",
primerReverse = "",
wildTypeForward = c(wt1 = "ACTGATACACTCCAAGCGGAGACAGACCAACTAGAAGATGAGAAGTCTGCTTTGCAGACCGAGATTGCCAACCTGCTGAAGGAGAAGGAAAAACTA",
wt2 = "CCTGATACACTCCAAGCGGAGACAGACCAACTAGAAGATGAGAAGTCTGCTTTGCAGACCGAGATTGCCAACCTGCTGAAGGAGAAGGAAAAACTA",
wt3 = "ACTGATACACTCCAAGCGGAGACTGACCAACTAGAAGATGAGAAGTCTTCTTTGCAGACCGAGATTGCCAACGTGCTGAAGGAGAAGGAAAAACTA"),
wildTypeReverse = "",
constantForward = "AACCGGAGGAGGGAGCTG",
constantReverse = "GAGTTCATCCTGGCAGC",
avePhredMinForward = 20.0, avePhredMinReverse = 20.0,
variableNMaxForward = 0, variableNMaxReverse = 0,
umiNMax = 0,
nbrMutatedCodonsMaxForward = -1,
nbrMutatedCodonsMaxReverse = -1,
nbrMutatedBasesMaxForward = 2,
nbrMutatedBasesMaxReverse = 2,
forbiddenMutatedCodonsForward = "NNW",
forbiddenMutatedCodonsReverse = "NNW",
useTreeWTmatch = FALSE,
mutatedPhredMinForward = 0.0, mutatedPhredMinReverse = 0.0,
mutNameDelimiter = ".",
constantMaxDistForward = -1,
constantMaxDistReverse = -1,
umiCollapseMaxDist = 0,
filteredReadsFastqForward = "",
filteredReadsFastqReverse = "",
maxNReads = -1, verbose = FALSE,
nThreads = 1, chunkSize = 1000,
maxReadLength = 1024
)
res1 <- do.call(digestFastqs, Ldef)
Ldef$useTreeWTmatch <- TRUE
res2 <- do.call(digestFastqs, Ldef)
expect_equal(res1$filterSummary$nbrTotal, res2$filterSummary$nbrTotal)
expect_equal(res1$filterSummary$f1_nbrAdapter, res2$filterSummary$f1_nbrAdapter)
expect_equal(res1$filterSummary$f2_nbrNoPrimer, res2$filterSummary$f2_nbrNoPrimer)
expect_equal(res1$filterSummary$f3_nbrReadWrongLength, res2$filterSummary$f3_nbrReadWrongLength)
expect_equal(res1$filterSummary$f4_nbrNoValidOverlap, res2$filterSummary$f4_nbrNoValidOverlap)
expect_equal(res1$filterSummary$f5_nbrAvgVarQualTooLow, res2$filterSummary$f5_nbrAvgVarQualTooLow)
expect_equal(res1$filterSummary$f6_nbrTooManyNinVar, res2$filterSummary$f6_nbrTooManyNinVar)
expect_equal(res1$filterSummary$f7_nbrTooManyNinUMI, res2$filterSummary$f7_nbrTooManyNinUMI)
expect_equal(res1$filterSummary$f8_nbrTooManyBestWTHits, res2$filterSummary$f8_nbrTooManyBestWTHits)
expect_equal(res1$filterSummary$f9_nbrMutQualTooLow, res2$filterSummary$f9_nbrMutQualTooLow)
expect_equal(res1$filterSummary$f10a_nbrTooManyMutCodons, res2$filterSummary$f10a_nbrTooManyMutCodons)
expect_equal(res1$filterSummary$f10b_nbrTooManyMutBases, res2$filterSummary$f10b_nbrTooManyMutBases)
expect_equal(res1$filterSummary$f11_nbrForbiddenCodons, res2$filterSummary$f11_nbrForbiddenCodons)
expect_equal(res1$filterSummary$f12_nbrTooManyMutConstant, res2$filterSummary$f12_nbrTooManyMutConstant)
expect_equal(res1$filterSummary$f13_nbrTooManyBestConstantHits, res2$filterSummary$f13_nbrTooManyBestConstantHits)
expect_equal(res1$filterSummary$nbrRetained, res2$filterSummary$nbrRetained)
expect_equal(res1$summaryTable$nbrReads, res2$summaryTable$nbrReads)
})
fmicompbio/mutscan documentation built on March 30, 2024, 9:13 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
test_that("digestFastqs works as expected for cis experiments", {
fqc1 <- system.file("extdata/cisInput_1.fastq.gz", package = "mutscan")
fqc2 <- system.file("extdata/cisInput_2.fastq.gz", package = "mutscan")
## default arguments
Ldef <- list(
fastqForward = fqc1, fastqReverse = fqc2,
mergeForwardReverse = TRUE,
minOverlap = 0, maxOverlap = 0, maxFracMismatchOverlap = 1, greedyOverlap = TRUE,
revComplForward = FALSE, revComplReverse = TRUE,
elementsForward = "SUCV", elementsReverse = "SUCVS",
elementLengthsForward = c(1, 10, 18, 96),
elementLengthsReverse = c(1, 7, 17, 96, -1),
adapterForward = "GGAAGAGCACACGTC",
adapterReverse = "GGAAGAGCGTCGTGT",
primerForward = "",
primerReverse = "",
wildTypeForward = "ACTGATACACTCCAAGCGGAGACAGACCAACTAGAAGATGAGAAGTCTGCTTTGCAGACCGAGATTGCCAACCTGCTGAAGGAGAAGGAAAAACTA",
wildTypeReverse = "",
constantForward = "AACCGGAGGAGGGAGCTG",
constantReverse = "GAGTTCATCCTGGCAGC",
avePhredMinForward = 20.0, avePhredMinReverse = 20.0,
variableNMaxForward = 0, variableNMaxReverse = 0,
umiNMax = 0,
nbrMutatedCodonsMaxForward = 1,
nbrMutatedCodonsMaxReverse = 1,
nbrMutatedBasesMaxForward = -1,
nbrMutatedBasesMaxReverse = -1,
forbiddenMutatedCodonsForward = "NNW",
forbiddenMutatedCodonsReverse = "NNW",
useTreeWTmatch = FALSE,
mutatedPhredMinForward = 0.0, mutatedPhredMinReverse = 0.0,
mutNameDelimiter = ".",
constantMaxDistForward = -1,
constantMaxDistReverse = -1,
umiCollapseMaxDist = 0,
filteredReadsFastqForward = "",
filteredReadsFastqReverse = "",
maxNReads = -1, verbose = FALSE,
nThreads = 1, chunkSize = 1000,
maxReadLength = 1024
)
res <- do.call(digestFastqs, Ldef)
## Specify reverse WT - check that it's ignored
LdefrevWT <- Ldef
LdefrevWT$wildTypeReverse <- "ACGT"
expect_warning(resrevWT <- do.call(digestFastqs, LdefrevWT),
"Ignoring 'wildTypeReverse'")
expect_equal(res$filterSummary, resrevWT$filterSummary)
expect_equal(res$summaryTable, resrevWT$summaryTable)
expect_equal(res$filterSummary$nbrTotal, 1000L)
expect_equal(res$filterSummary$f1_nbrAdapter, 126L)
expect_equal(res$filterSummary$f2_nbrNoPrimer, 0L)
expect_equal(res$filterSummary$f3_nbrReadWrongLength, 0L)
expect_equal(res$filterSummary$f4_nbrNoValidOverlap, 0L)
expect_equal(res$filterSummary$f5_nbrAvgVarQualTooLow, 0L)
expect_equal(res$filterSummary$f6_nbrTooManyNinVar, 44L)
expect_equal(res$filterSummary$f7_nbrTooManyNinUMI, 0L)
expect_equal(res$filterSummary$f8_nbrTooManyBestWTHits, 0L)
expect_equal(res$filterSummary$f9_nbrMutQualTooLow, 0L)
expect_equal(res$filterSummary$f10a_nbrTooManyMutCodons, 581L)
expect_equal(res$filterSummary$f10b_nbrTooManyMutBases, 0L)
expect_equal(res$filterSummary$f11_nbrForbiddenCodons, 82L)
expect_equal(res$filterSummary$f12_nbrTooManyMutConstant, 0L)
expect_equal(res$filterSummary$f13_nbrTooManyBestConstantHits, 0L)
expect_equal(res$filterSummary$nbrRetained, 167)
for (nm in setdiff(names(Ldef), c("forbiddenMutatedCodonsForward", "forbiddenMutatedCodonsReverse", "verbose", "fastqForward", "fastqReverse"))) {
expect_equal(res$parameters[[nm]], Ldef[[nm]], ignore_attr = TRUE)
}
for (nm in c("fastqForward", "fastqReverse")) {
expect_equal(res$parameters[[nm]], normalizePath(Ldef[[nm]], mustWork = FALSE),
ignore_attr = TRUE)
}
expect_equal(sum(res$summaryTable$nbrReads), res$filterSummary$nbrRetained)
expect_equal(sum(res$summaryTable$nbrReads == 2), 11L)
expect_equal(sort(res$summaryTable$mutantName[res$summaryTable$nbrReads == 2]),
sort(c("f.14.AAG", "f.15.ATG", "f.19.CAC", "f.1.ACC", "f.20.AAC", "f.21.GTG",
"f.24.AGC", "f.4.CGC", "f.7.GTG", "f.9.GGC", "f.9.GTC")))
expect_true(all(res$summaryTable$nbrReads == res$summaryTable$nbrUmis))
## Check the number of reads with a given number of mismatches
expect_equal(sum(res$summaryTable$nbrReads[res$summaryTable$nbrMutBases == 0]), 77L)
expect_equal(sum(res$summaryTable$nbrReads[res$summaryTable$nbrMutBases == 1]), 88L)
expect_equal(sum(res$summaryTable$nbrReads[res$summaryTable$nbrMutBases == 2]), 2L)
expect_equal(sum(res$summaryTable$nbrReads[res$summaryTable$nbrMutCodons == 0]), 77L)
expect_equal(sum(res$summaryTable$nbrReads[res$summaryTable$nbrMutCodons == 1]), 90L)
expect_equal(sum(res$summaryTable$nbrReads[res$summaryTable$nbrMutAAs == 0]), 97L)
expect_equal(sum(res$summaryTable$nbrReads[res$summaryTable$nbrMutAAs == 1]), 70L)
## Similarly but for the number of unique sequences
expect_equal(sum(res$summaryTable$nbrMutBases == 0), 1L)
expect_equal(sum(res$summaryTable$nbrMutBases == 1), 64L)
expect_equal(sum(res$summaryTable$nbrMutBases == 2), 2L)
expect_equal(sum(res$summaryTable$nbrMutCodons == 0), 1L)
expect_equal(sum(res$summaryTable$nbrMutCodons == 1), 66L)
expect_equal(sum(res$summaryTable$nbrMutAAs == 0), 14L)
expect_equal(sum(res$summaryTable$nbrMutAAs == 1), 53L)
## Check that mutant naming worked (compare to manual matching)
example_seq <- paste0("ACTGATACACGCCAAGCGGAGACAGACCAACTAGAAGATGAGAAGTCT",
"GCTTTGCAGACCGAGATTGCCAACCTGCTGAAGGAGAAGGAAAAACTA")
expect_equal(res$summaryTable$mutantName[res$summaryTable$sequence == example_seq],
"f.4.CGC")
expect_equal(res$summaryTable$mutantNameAA[res$summaryTable$sequence == example_seq],
"f.4.R")
expect_equal(res$summaryTable$mutantNameBase[res$summaryTable$sequence == example_seq],
"f.11.G")
expect_equal(res$summaryTable$mutantNameCodon[res$summaryTable$sequence == example_seq],
"f.4.CGC")
expect_equal(res$summaryTable$mutantNameBaseHGVS[res$summaryTable$sequence == example_seq],
"f:c.11T>G")
expect_equal(res$summaryTable$mutantNameAAHGVS[res$summaryTable$sequence == example_seq],
"f:p.(Leu4Arg)")
expect_equal(res$summaryTable$sequenceAA,
as.character(Biostrings::translate(
Biostrings::DNAStringSet(res$summaryTable$sequence))))
expect_equal(as.numeric(res$summaryTable$varLengths),
nchar(res$summaryTable$sequence))
expect_equal(sum(grepl("WT", res$summaryTable$mutantNameAA) & res$summaryTable$nbrMutBases > 0), 13)
expect_true(all(grepl("silent", res$summaryTable$mutationTypes[grepl("WT", res$summaryTable$mutantNameAA) & res$summaryTable$nbrMutBases > 0])))
expect_true(all(grepl("nonsynonymous|stop", res$summaryTable$mutationTypes[grepl("\\.[1-9][0-9]*\\.", res$summaryTable$mutantNameAA)])))
expect_true(all(grepl("stop", res$summaryTable$mutationTypes[grepl("TAG|TAA|TGA", res$summaryTable$mutantName)])))
expect_equal(sum(res$errorStatistics$nbrMatchForward + res$errorStatistics$nbrMismatchForward),
nchar(Ldef$constantForward[1]) * res$filterSummary$nbrRetained)
expect_equal(sum(res$errorStatistics$nbrMatchReverse + res$errorStatistics$nbrMismatchReverse),
nchar(Ldef$constantReverse[1]) * res$filterSummary$nbrRetained)
expect_equal(res$errorStatistics$nbrMatchForward[res$errorStatistics$PhredQuality == 14], 51L)
expect_equal(res$errorStatistics$nbrMatchForward[res$errorStatistics$PhredQuality == 22], 19L)
expect_equal(res$errorStatistics$nbrMatchForward[res$errorStatistics$PhredQuality == 27], 79L)
expect_equal(res$errorStatistics$nbrMatchForward[res$errorStatistics$PhredQuality == 33], 163L)
expect_equal(res$errorStatistics$nbrMatchForward[res$errorStatistics$PhredQuality == 37], 2690L)
expect_equal(res$errorStatistics$nbrMismatchForward[res$errorStatistics$PhredQuality == 14], 1L)
expect_equal(res$errorStatistics$nbrMismatchForward[res$errorStatistics$PhredQuality == 27], 1L)
expect_equal(res$errorStatistics$nbrMismatchForward[res$errorStatistics$PhredQuality == 37], 2L)
expect_equal(res$errorStatistics$nbrMatchReverse[res$errorStatistics$PhredQuality == 14], 41L)
expect_equal(res$errorStatistics$nbrMatchReverse[res$errorStatistics$PhredQuality == 27], 56L)
expect_equal(res$errorStatistics$nbrMatchReverse[res$errorStatistics$PhredQuality == 33], 96L)
expect_equal(res$errorStatistics$nbrMatchReverse[res$errorStatistics$PhredQuality == 37], 2620L)
expect_equal(res$errorStatistics$nbrMismatchReverse[res$errorStatistics$PhredQuality == 2], 4L)
expect_equal(res$errorStatistics$nbrMismatchReverse[res$errorStatistics$PhredQuality == 14], 7L)
expect_equal(res$errorStatistics$nbrMismatchReverse[res$errorStatistics$PhredQuality == 33], 1L)
expect_equal(res$errorStatistics$nbrMismatchReverse[res$errorStatistics$PhredQuality == 37], 14L)
## Add an additional constant sequence of the wrong length - should fail
Ldefcs <- Ldef
Ldefcs$constantForward <- c("AACCGGAGGAGGGAGCTG", "AACCGGAGGAGGGAGCT")
expect_error(do.call(digestFastqs, Ldefcs),
"All constant sequences must be of the same length")
})
test_that("digestFastqs works as expected when specifying max nbr of mutated bases rather than codons", {
fqc1 <- system.file("extdata/cisInput_1.fastq.gz", package = "mutscan")
fqc2 <- system.file("extdata/cisInput_2.fastq.gz", package = "mutscan")
## default arguments
Ldef <- list(
fastqForward = fqc1, fastqReverse = fqc2,
mergeForwardReverse = TRUE,
minOverlap = 0, maxOverlap = 0, maxFracMismatchOverlap = 1, greedyOverlap = TRUE,
revComplForward = FALSE, revComplReverse = TRUE,
elementsForward = "SUCV", elementsReverse = "SUCVS",
elementLengthsForward = c(1, 10, 18, 96),
elementLengthsReverse = c(1, 7, 17, 96, -1),
adapterForward = "GGAAGAGCACACGTC",
adapterReverse = "GGAAGAGCGTCGTGT",
primerForward = "",
primerReverse = "",
wildTypeForward = "ACTGATACACTCCAAGCGGAGACAGACCAACTAGAAGATGAGAAGTCTGCTTTGCAGACCGAGATTGCCAACCTGCTGAAGGAGAAGGAAAAACTA",
wildTypeReverse = "",
constantForward = "AACCGGAGGAGGGAGCTG",
constantReverse = "GAGTTCATCCTGGCAGC",
avePhredMinForward = 20.0, avePhredMinReverse = 20.0,
variableNMaxForward = 0, variableNMaxReverse = 0,
umiNMax = 0,
nbrMutatedCodonsMaxForward = -1,
nbrMutatedCodonsMaxReverse = -1,
nbrMutatedBasesMaxForward = 2,
nbrMutatedBasesMaxReverse = 2,
forbiddenMutatedCodonsForward = "NNW",
forbiddenMutatedCodonsReverse = "NNW",
useTreeWTmatch = FALSE,
mutatedPhredMinForward = 0.0, mutatedPhredMinReverse = 0.0,
mutNameDelimiter = ".",
constantMaxDistForward = -1,
constantMaxDistReverse = -1,
umiCollapseMaxDist = 0,
filteredReadsFastqForward = "",
filteredReadsFastqReverse = "",
maxNReads = -1, verbose = FALSE,
nThreads = 1, chunkSize = 1000,
maxReadLength = 125
)
res <- do.call(digestFastqs, Ldef)
expect_equal(res$filterSummary$nbrTotal, 1000L)
expect_equal(res$filterSummary$f1_nbrAdapter, 126L)
expect_equal(res$filterSummary$f2_nbrNoPrimer, 0L)
expect_equal(res$filterSummary$f3_nbrReadWrongLength, 0L)
expect_equal(res$filterSummary$f4_nbrNoValidOverlap, 0L)
expect_equal(res$filterSummary$f5_nbrAvgVarQualTooLow, 0L)
expect_equal(res$filterSummary$f6_nbrTooManyNinVar, 44L)
expect_equal(res$filterSummary$f7_nbrTooManyNinUMI, 0L)
expect_equal(res$filterSummary$f8_nbrTooManyBestWTHits, 0L)
expect_equal(res$filterSummary$f9_nbrMutQualTooLow, 0L)
expect_equal(res$filterSummary$f10a_nbrTooManyMutCodons, 0L)
expect_equal(res$filterSummary$f10b_nbrTooManyMutBases, 416L)
expect_equal(res$filterSummary$f11_nbrForbiddenCodons, 212L)
expect_equal(res$filterSummary$f12_nbrTooManyMutConstant, 0L)
expect_equal(res$filterSummary$f13_nbrTooManyBestConstantHits, 0L)
expect_equal(res$filterSummary$nbrRetained, 202L)
expect_named(res$summaryTable, c("mutantName", "sequence","nbrReads",
"maxNbrReads", "nbrUmis", "nbrMutBases",
"nbrMutCodons", "nbrMutAAs", "varLengths",
"mutantNameBase", "mutantNameCodon",
"mutantNameBaseHGVS", "mutantNameAA",
"mutantNameAAHGVS", "mutationTypes",
"sequenceAA"))
for (nm in setdiff(names(Ldef), c("forbiddenMutatedCodonsForward", "forbiddenMutatedCodonsReverse", "verbose", "fastqForward", "fastqReverse"))) {
expect_equal(res$parameters[[nm]], Ldef[[nm]], ignore_attr = TRUE)
}
for (nm in c("fastqForward", "fastqReverse")) {
expect_equal(res$parameters[[nm]], normalizePath(Ldef[[nm]], mustWork = FALSE),
ignore_attr = TRUE)
}
expect_equal(sum(res$summaryTable$nbrReads), res$filterSummary$nbrRetained)
expect_equal(sum(res$summaryTable$nbrReads == 3), 5L)
expect_equal(sort(res$summaryTable$mutantName[res$summaryTable$nbrReads == 3]),
sort(c("f.12.G", "f.17.T", "f.60.G", "f.85.T", "f.96.G")))
expect_true(all(res$summaryTable$nbrReads == res$summaryTable$nbrUmis))
## Check the number of reads with a given number of mismatches
expect_equal(sum(res$summaryTable$nbrReads[res$summaryTable$nbrMutBases == 0]), 77L)
expect_equal(sum(res$summaryTable$nbrReads[res$summaryTable$nbrMutBases == 1]), 88L)
expect_equal(sum(res$summaryTable$nbrReads[res$summaryTable$nbrMutBases == 2]), 37L)
expect_equal(sum(res$summaryTable$nbrReads[res$summaryTable$nbrMutCodons == 0]), 77L)
expect_equal(sum(res$summaryTable$nbrReads[res$summaryTable$nbrMutCodons == 1]), 90L)
expect_equal(sum(res$summaryTable$nbrReads[res$summaryTable$nbrMutCodons == 2]), 35L)
expect_equal(sum(res$summaryTable$nbrReads[res$summaryTable$nbrMutAAs == 0]), 98L)
expect_equal(sum(res$summaryTable$nbrReads[res$summaryTable$nbrMutAAs == 1]), 86L)
expect_equal(sum(res$summaryTable$nbrReads[res$summaryTable$nbrMutAAs == 2]), 18L)
## Similarly but for the number of unique sequences
expect_equal(sum(res$summaryTable$nbrMutBases == 0), 1L)
expect_equal(sum(res$summaryTable$nbrMutBases == 1), 64L)
expect_equal(sum(res$summaryTable$nbrMutBases == 2), 36L)
expect_equal(sum(res$summaryTable$nbrMutCodons == 0), 1L)
expect_equal(sum(res$summaryTable$nbrMutCodons == 1), 66L)
expect_equal(sum(res$summaryTable$nbrMutCodons == 2), 34L)
expect_equal(sum(res$summaryTable$nbrMutAAs == 0), 15L)
expect_equal(sum(res$summaryTable$nbrMutAAs == 1), 69L)
expect_equal(sum(res$summaryTable$nbrMutAAs == 2), 17L)
## Check that mutant naming worked (compare to manual matching)
example_seq <- paste0("ACTGATACACGCCAAGCGGAGACAGACCAACTAGAAGATGAGAAGTCT",
"GCTTTGCAGACCGAGATTGCCAACCTGCTGAAGGAGAAGGAAAAACTA")
expect_equal(res$summaryTable$mutantName[res$summaryTable$sequence == example_seq],
"f.11.G")
expect_equal(res$summaryTable$mutantNameAA[res$summaryTable$sequence == example_seq],
"f.4.R")
expect_equal(res$summaryTable$mutantNameBase[res$summaryTable$sequence == example_seq],
"f.11.G")
expect_equal(res$summaryTable$mutantNameCodon[res$summaryTable$sequence == example_seq],
"f.4.CGC")
expect_equal(res$summaryTable$mutantNameBaseHGVS[res$summaryTable$sequence == example_seq],
"f:c.11T>G")
expect_equal(res$summaryTable$mutantNameAAHGVS[res$summaryTable$sequence == example_seq],
"f:p.(Leu4Arg)")
expect_equal(sum(res$errorStatistics$nbrMatchForward + res$errorStatistics$nbrMismatchForward),
nchar(Ldef$constantForward[1]) * res$filterSummary$nbrRetained)
expect_equal(sum(res$errorStatistics$nbrMatchReverse + res$errorStatistics$nbrMismatchReverse),
nchar(Ldef$constantReverse[1]) * res$filterSummary$nbrRetained)
expect_equal(res$errorStatistics$nbrMatchForward[res$errorStatistics$PhredQuality == 14], 61L)
expect_equal(res$errorStatistics$nbrMatchForward[res$errorStatistics$PhredQuality == 22], 22L)
expect_equal(res$errorStatistics$nbrMatchForward[res$errorStatistics$PhredQuality == 27], 97L)
expect_equal(res$errorStatistics$nbrMatchForward[res$errorStatistics$PhredQuality == 33], 203L)
expect_equal(res$errorStatistics$nbrMatchForward[res$errorStatistics$PhredQuality == 37], 3249L)
expect_equal(res$errorStatistics$nbrMismatchForward[res$errorStatistics$PhredQuality == 14], 1L)
expect_equal(res$errorStatistics$nbrMismatchForward[res$errorStatistics$PhredQuality == 27], 1L)
expect_equal(res$errorStatistics$nbrMismatchForward[res$errorStatistics$PhredQuality == 37], 2L)
expect_equal(res$errorStatistics$nbrMatchReverse[res$errorStatistics$PhredQuality == 14], 42L)
expect_equal(res$errorStatistics$nbrMatchReverse[res$errorStatistics$PhredQuality == 27], 63L)
expect_equal(res$errorStatistics$nbrMatchReverse[res$errorStatistics$PhredQuality == 33], 117L)
expect_equal(res$errorStatistics$nbrMatchReverse[res$errorStatistics$PhredQuality == 37], 3184L)
expect_equal(res$errorStatistics$nbrMismatchReverse[res$errorStatistics$PhredQuality == 2], 5L)
expect_equal(res$errorStatistics$nbrMismatchReverse[res$errorStatistics$PhredQuality == 14], 7L)
expect_equal(res$errorStatistics$nbrMismatchReverse[res$errorStatistics$PhredQuality == 33], 2L)
expect_equal(res$errorStatistics$nbrMismatchReverse[res$errorStatistics$PhredQuality == 37], 14L)
})
test_that("digestFastqs gives the same results regardless of how WT matching is done", {
fqc1 <- system.file("extdata/cisInput_1.fastq.gz", package = "mutscan")
fqc2 <- system.file("extdata/cisInput_2.fastq.gz", package = "mutscan")
## default arguments
Ldef <- list(
fastqForward = fqc1, fastqReverse = fqc2,
mergeForwardReverse = TRUE,
minOverlap = 0, maxOverlap = 0, maxFracMismatchOverlap = 1, greedyOverlap = TRUE,
revComplForward = FALSE, revComplReverse = TRUE,
elementsForward = "SUCV", elementsReverse = "SUCVS",
elementLengthsForward = c(1, 10, 18, 96),
elementLengthsReverse = c(1, 7, 17, 96, -1),
adapterForward = "GGAAGAGCACACGTC",
adapterReverse = "GGAAGAGCGTCGTGT",
primerForward = "",
primerReverse = "",
wildTypeForward = c(wt1 = "ACTGATACACTCCAAGCGGAGACAGACCAACTAGAAGATGAGAAGTCTGCTTTGCAGACCGAGATTGCCAACCTGCTGAAGGAGAAGGAAAAACTA",
wt2 = "CCTGATACACTCCAAGCGGAGACAGACCAACTAGAAGATGAGAAGTCTGCTTTGCAGACCGAGATTGCCAACCTGCTGAAGGAGAAGGAAAAACTA",
wt3 = "ACTGATACACTCCAAGCGGAGACTGACCAACTAGAAGATGAGAAGTCTTCTTTGCAGACCGAGATTGCCAACGTGCTGAAGGAGAAGGAAAAACTA"),
wildTypeReverse = "",
constantForward = "AACCGGAGGAGGGAGCTG",
constantReverse = "GAGTTCATCCTGGCAGC",
avePhredMinForward = 20.0, avePhredMinReverse = 20.0,
variableNMaxForward = 0, variableNMaxReverse = 0,
umiNMax = 0,
nbrMutatedCodonsMaxForward = -1,
nbrMutatedCodonsMaxReverse = -1,
nbrMutatedBasesMaxForward = 2,
nbrMutatedBasesMaxReverse = 2,
forbiddenMutatedCodonsForward = "NNW",
forbiddenMutatedCodonsReverse = "NNW",
useTreeWTmatch = FALSE,
mutatedPhredMinForward = 0.0, mutatedPhredMinReverse = 0.0,
mutNameDelimiter = ".",
constantMaxDistForward = -1,
constantMaxDistReverse = -1,
umiCollapseMaxDist = 0,
filteredReadsFastqForward = "",
filteredReadsFastqReverse = "",
maxNReads = -1, verbose = FALSE,
nThreads = 1, chunkSize = 1000,
maxReadLength = 1024
)
res1 <- do.call(digestFastqs, Ldef)
Ldef$useTreeWTmatch <- TRUE
res2 <- do.call(digestFastqs, Ldef)
expect_equal(res1$filterSummary$nbrTotal, res2$filterSummary$nbrTotal)
expect_equal(res1$filterSummary$f1_nbrAdapter, res2$filterSummary$f1_nbrAdapter)
expect_equal(res1$filterSummary$f2_nbrNoPrimer, res2$filterSummary$f2_nbrNoPrimer)
expect_equal(res1$filterSummary$f3_nbrReadWrongLength, res2$filterSummary$f3_nbrReadWrongLength)
expect_equal(res1$filterSummary$f4_nbrNoValidOverlap, res2$filterSummary$f4_nbrNoValidOverlap)
expect_equal(res1$filterSummary$f5_nbrAvgVarQualTooLow, res2$filterSummary$f5_nbrAvgVarQualTooLow)
expect_equal(res1$filterSummary$f6_nbrTooManyNinVar, res2$filterSummary$f6_nbrTooManyNinVar)
expect_equal(res1$filterSummary$f7_nbrTooManyNinUMI, res2$filterSummary$f7_nbrTooManyNinUMI)
expect_equal(res1$filterSummary$f8_nbrTooManyBestWTHits, res2$filterSummary$f8_nbrTooManyBestWTHits)
expect_equal(res1$filterSummary$f9_nbrMutQualTooLow, res2$filterSummary$f9_nbrMutQualTooLow)
expect_equal(res1$filterSummary$f10a_nbrTooManyMutCodons, res2$filterSummary$f10a_nbrTooManyMutCodons)
expect_equal(res1$filterSummary$f10b_nbrTooManyMutBases, res2$filterSummary$f10b_nbrTooManyMutBases)
expect_equal(res1$filterSummary$f11_nbrForbiddenCodons, res2$filterSummary$f11_nbrForbiddenCodons)
expect_equal(res1$filterSummary$f12_nbrTooManyMutConstant, res2$filterSummary$f12_nbrTooManyMutConstant)
expect_equal(res1$filterSummary$f13_nbrTooManyBestConstantHits, res2$filterSummary$f13_nbrTooManyBestConstantHits)
expect_equal(res1$filterSummary$nbrRetained, res2$filterSummary$nbrRetained)
expect_equal(res1$summaryTable$nbrReads, res2$summaryTable$nbrReads)
})
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.