@@ -1261,3 +1261,153 @@ test_that("digestFastqs works as expected for trans experiments, if collapsing t
12611261 expect_equal(res $ summaryTable $ varLengths , " 80,16_16,80" )
12621262})
12631263
1264+ test_that(" reads are filtered out if there are too many constant hits" , {
1265+ fqt1 <- system.file(" extdata/transInput_1.fastq.gz" , package = " mutscan" )
1266+ fqt2 <- system.file(" extdata/transInput_2.fastq.gz" , package = " mutscan" )
1267+ # # default arguments
1268+ Ldef <- list (
1269+ fastqForward = fqt1 , fastqReverse = fqt2 ,
1270+ mergeForwardReverse = FALSE ,
1271+ minOverlap = 0 , maxOverlap = 0 , maxFracMismatchOverlap = 0 , greedyOverlap = TRUE ,
1272+ revComplForward = FALSE , revComplReverse = FALSE ,
1273+ elementsForward = " SUCVV" , elementsReverse = " SUCVV" ,
1274+ elementLengthsForward = c(1 , 10 , 18 , 80 , 16 ),
1275+ elementLengthsReverse = c(1 , 8 , 20 , 16 , 80 ),
1276+ adapterForward = " GGAAGAGCACACGTC" ,
1277+ adapterReverse = " GGAAGAGCGTCGTGT" ,
1278+ primerForward = " " ,
1279+ primerReverse = " " ,
1280+ wildTypeForward = " ACTGATACACTCCAAGCGGAGACAGACCAACTAGAAGATGAGAAGTCTGCTTTGCAGACCGAGATTGCCAACCTGCTGAAGGAGAAGGAAAAACTA" ,
1281+ wildTypeReverse = " ATCGCCCGGCTGGAGGAAAAAGTGAAAACCTTGAAAGCTCAGAACTCGGAGCTGGCGTCCACGGCCAACATGCTCAGGGAACAGGTGGCACAGCTT" ,
1282+ constantForward = " AACCGGAGGAGGGAGCTG" ,
1283+ constantReverse = " GAAAAAGGAAGCTGGAGAGA" ,
1284+ avePhredMinForward = 20.0 , avePhredMinReverse = 20.0 ,
1285+ variableNMaxForward = 0 , variableNMaxReverse = 0 ,
1286+ umiNMax = 0 ,
1287+ nbrMutatedCodonsMaxForward = 1 ,
1288+ nbrMutatedCodonsMaxReverse = 1 ,
1289+ nbrMutatedBasesMaxForward = - 1 ,
1290+ nbrMutatedBasesMaxReverse = - 1 ,
1291+ forbiddenMutatedCodonsForward = " NNW" ,
1292+ forbiddenMutatedCodonsReverse = " NNW" ,
1293+ useTreeWTmatch = FALSE ,
1294+ collapseToWTForward = TRUE ,
1295+ collapseToWTReverse = TRUE ,
1296+ mutatedPhredMinForward = 0.0 , mutatedPhredMinReverse = 0.0 ,
1297+ mutNameDelimiter = " ." ,
1298+ constantMaxDistForward = - 1 ,
1299+ constantMaxDistReverse = - 1 ,
1300+ umiCollapseMaxDist = 0 ,
1301+ filteredReadsFastqForward = " " ,
1302+ filteredReadsFastqReverse = " " ,
1303+ maxNReads = - 1 , verbose = FALSE ,
1304+ nThreads = 1 , chunkSize = 1000 ,
1305+ maxReadLength = 1024
1306+ )
1307+
1308+ # # Forward
1309+ Ldef1 <- Ldef ; Ldef1 $ constantForward <- c(c1 = " AACCGGAGGAGGGAGCTA" , c2 = " AACCGGAGGAGGGAGCTC" )
1310+ res1 <- do.call(digestFastqs , Ldef1 )
1311+
1312+ expect_equal(res1 $ filterSummary $ nbrTotal , 1000L )
1313+ expect_equal(res1 $ filterSummary $ f1_nbrAdapter , 314L )
1314+ expect_equal(res1 $ filterSummary $ f2_nbrNoPrimer , 0L )
1315+ expect_equal(res1 $ filterSummary $ f3_nbrReadWrongLength , 0L )
1316+ expect_equal(res1 $ filterSummary $ f4_nbrNoValidOverlap , 0L )
1317+ expect_equal(res1 $ filterSummary $ f5_nbrAvgVarQualTooLow , 7L )
1318+ expect_equal(res1 $ filterSummary $ f6_nbrTooManyNinVar , 0L )
1319+ expect_equal(res1 $ filterSummary $ f7_nbrTooManyNinUMI , 0L )
1320+ expect_equal(res1 $ filterSummary $ f8_nbrTooManyBestWTHits , 0L )
1321+ expect_equal(res1 $ filterSummary $ f9_nbrMutQualTooLow , 0L )
1322+ expect_equal(res1 $ filterSummary $ f10a_nbrTooManyMutCodons , 192L + 95L + 68L + 37L )
1323+ expect_equal(res1 $ filterSummary $ f10b_nbrTooManyMutBases , 0L )
1324+ expect_equal(res1 $ filterSummary $ f11_nbrForbiddenCodons , 6L + 2L )
1325+ expect_equal(res1 $ filterSummary $ f12_nbrTooManyMutConstant , 0L )
1326+ expect_equal(res1 $ filterSummary $ f13_nbrTooManyBestConstantHits , 279L )
1327+ expect_equal(res1 $ filterSummary $ nbrRetained , 0L )
1328+
1329+ # # Reverse
1330+ Ldef1 <- Ldef ; Ldef1 $ constantReverse <- c(c1 = " GAAAAAGGAAGCTGGAGAGG" , c2 = " GAAAAAGGAAGCTGGAGAGT" )
1331+ res1 <- do.call(digestFastqs , Ldef1 )
1332+
1333+ expect_equal(res1 $ filterSummary $ nbrTotal , 1000L )
1334+ expect_equal(res1 $ filterSummary $ f1_nbrAdapter , 314L )
1335+ expect_equal(res1 $ filterSummary $ f2_nbrNoPrimer , 0L )
1336+ expect_equal(res1 $ filterSummary $ f3_nbrReadWrongLength , 0L )
1337+ expect_equal(res1 $ filterSummary $ f4_nbrNoValidOverlap , 0L )
1338+ expect_equal(res1 $ filterSummary $ f5_nbrAvgVarQualTooLow , 7L )
1339+ expect_equal(res1 $ filterSummary $ f6_nbrTooManyNinVar , 0L )
1340+ expect_equal(res1 $ filterSummary $ f7_nbrTooManyNinUMI , 0L )
1341+ expect_equal(res1 $ filterSummary $ f8_nbrTooManyBestWTHits , 0L )
1342+ expect_equal(res1 $ filterSummary $ f9_nbrMutQualTooLow , 0L )
1343+ expect_equal(res1 $ filterSummary $ f10a_nbrTooManyMutCodons , 192L + 95L + 68L + 37L )
1344+ expect_equal(res1 $ filterSummary $ f10b_nbrTooManyMutBases , 0L )
1345+ expect_equal(res1 $ filterSummary $ f11_nbrForbiddenCodons , 6L + 2L )
1346+ expect_equal(res1 $ filterSummary $ f12_nbrTooManyMutConstant , 0L )
1347+ expect_equal(res1 $ filterSummary $ f13_nbrTooManyBestConstantHits , 279L )
1348+ expect_equal(res1 $ filterSummary $ nbrRetained , 0L )
1349+ })
1350+
1351+ test_that(" digestFastqs runs with revComplForward = TRUE" , {
1352+ fqt1 <- system.file(" extdata/transInput_1.fastq.gz" , package = " mutscan" )
1353+ fqt2 <- system.file(" extdata/transInput_2.fastq.gz" , package = " mutscan" )
1354+ # # default arguments
1355+ Ldef <- list (
1356+ fastqForward = fqt1 , fastqReverse = fqt2 ,
1357+ mergeForwardReverse = FALSE ,
1358+ minOverlap = 0 , maxOverlap = 0 , maxFracMismatchOverlap = 0 , greedyOverlap = TRUE ,
1359+ revComplForward = TRUE , revComplReverse = FALSE ,
1360+ elementsForward = " SUCVV" , elementsReverse = " SUCVV" ,
1361+ elementLengthsForward = c(1 , 10 , 18 , 80 , 16 ),
1362+ elementLengthsReverse = c(1 , 8 , 20 , 16 , 80 ),
1363+ adapterForward = " GGAAGAGCACACGTC" ,
1364+ adapterReverse = " GGAAGAGCGTCGTGT" ,
1365+ primerForward = " " ,
1366+ primerReverse = " " ,
1367+ wildTypeForward = " TAGTTTTTCCTTCTCCTTCAGCAGGTTGGCAATCTCGGTCTGCAAAGCAGACTTCTCATCTTCTAGTTGGTCTGTCTCCGCTTGGAGTGTATCAGT" ,
1368+ wildTypeReverse = " ATCGCCCGGCTGGAGGAAAAAGTGAAAACCTTGAAAGCTCAGAACTCGGAGCTGGCGTCCACGGCCAACATGCTCAGGGAACAGGTGGCACAGCTT" ,
1369+ constantForward = " CAGCTCCCTCCTCCGGTT" ,
1370+ constantReverse = " GAAAAAGGAAGCTGGAGAGA" ,
1371+ avePhredMinForward = 20.0 , avePhredMinReverse = 20.0 ,
1372+ variableNMaxForward = 0 , variableNMaxReverse = 0 ,
1373+ umiNMax = 0 ,
1374+ nbrMutatedCodonsMaxForward = 1 ,
1375+ nbrMutatedCodonsMaxReverse = 1 ,
1376+ nbrMutatedBasesMaxForward = - 1 ,
1377+ nbrMutatedBasesMaxReverse = - 1 ,
1378+ forbiddenMutatedCodonsForward = " " ,
1379+ forbiddenMutatedCodonsReverse = " " ,
1380+ useTreeWTmatch = FALSE ,
1381+ collapseToWTForward = TRUE ,
1382+ collapseToWTReverse = TRUE ,
1383+ mutatedPhredMinForward = 0.0 , mutatedPhredMinReverse = 0.0 ,
1384+ mutNameDelimiter = " ." ,
1385+ constantMaxDistForward = - 1 ,
1386+ constantMaxDistReverse = - 1 ,
1387+ umiCollapseMaxDist = 0 ,
1388+ filteredReadsFastqForward = " " ,
1389+ filteredReadsFastqReverse = " " ,
1390+ maxNReads = - 1 , verbose = FALSE ,
1391+ nThreads = 1 , chunkSize = 1000 ,
1392+ maxReadLength = 1024
1393+ )
1394+
1395+ res <- do.call(digestFastqs , Ldef )
1396+
1397+ expect_equal(res $ filterSummary $ nbrTotal , 1000L )
1398+ expect_equal(res $ filterSummary $ f1_nbrAdapter , 314L )
1399+ expect_equal(res $ filterSummary $ f2_nbrNoPrimer , 0L )
1400+ expect_equal(res $ filterSummary $ f3_nbrReadWrongLength , 0L )
1401+ expect_equal(res $ filterSummary $ f4_nbrNoValidOverlap , 0L )
1402+ expect_equal(res $ filterSummary $ f5_nbrAvgVarQualTooLow , 7L )
1403+ expect_equal(res $ filterSummary $ f6_nbrTooManyNinVar , 0L )
1404+ expect_equal(res $ filterSummary $ f7_nbrTooManyNinUMI , 0L )
1405+ expect_equal(res $ filterSummary $ f8_nbrTooManyBestWTHits , 0L )
1406+ expect_equal(res $ filterSummary $ f9_nbrMutQualTooLow , 0L )
1407+ expect_equal(res $ filterSummary $ f10a_nbrTooManyMutCodons , 394L )
1408+ expect_equal(res $ filterSummary $ f10b_nbrTooManyMutBases , 0L )
1409+ expect_equal(res $ filterSummary $ f11_nbrForbiddenCodons , 0L )
1410+ expect_equal(res $ filterSummary $ f12_nbrTooManyMutConstant , 0L )
1411+ expect_equal(res $ filterSummary $ f13_nbrTooManyBestConstantHits , 0L )
1412+ expect_equal(res $ filterSummary $ nbrRetained , 285L )
1413+ })
0 commit comments