Add option to return all lrs

duke-malaria-collaboratory · Sep 7, 2023 · 10815bd · 10815bd
1 parent 9a2d61d
commit 10815bd
Show file tree

Hide file tree

Showing 4 changed files with 72 additions and 2 deletions.
diff --git a/R/bistro.R b/R/bistro.R
@@ -53,6 +53,8 @@
 #' @param time_limit Time limit in minutes to run the
 #'   [euroformix::contLikSearch()] function on 1 bloodmeal-human pair. Default:
 #'   3
+#' @param return_lrs A boolean indicating whether or not to return log10LRs for
+#'   all bloodmeal-human pairs. Default: FALSE
 #'
 #' @return Tibble with matches for bloodmeal-human pairs including the columns
 #'   listed below. Note that if multiple matches are found for a bloodmeal,
@@ -66,6 +68,12 @@
 #' * `log10_lr`: If match, log10 likelihood ratio (NA otherwise)
 #' * `notes`: Why the bloodmeal does or doesn't have a match
 #'
+#'   If `return_lrs = TRUE`, then a named list of length 2 is returned:
+#' * matches - the tibble described above
+#' * lrs - log10LRs for each bloodmeal-human pair including some
+#' of the columns described above and an additional column:
+#'  `efm_noc`, which is the number of contributors used as input
+#'   into euroformix, which is `min(est_noc, 3)`.
 #'
 #' @export
 #'
@@ -90,7 +98,8 @@ bistro <-
            difftol = 1,
            threads = 4,
            seed = 1,
-           time_limit = 3) {
+           time_limit = 3,
+           return_lrs = FALSE) {
     check_bistro_inputs(
       bloodmeal_profiles,
       human_profiles,
@@ -152,5 +161,10 @@ bistro <-
 
     matches <- identify_matches(log10_lrs, bloodmeal_ids)
 
+    if(return_lrs){
+      matches <- list(matches = matches,
+                      lrs = log10_lrs)
+    }
+
     return(matches)
   }
diff --git a/man/bistro.Rd b/man/bistro.Rd
diff --git a/tests/testthat/_snaps/bistro.md b/tests/testthat/_snaps/bistro.md
@@ -40,3 +40,38 @@
         <chr>              <int>   <dbl> <chr> <chr>       <dbl> <chr>           <dbl>
       1 evid1                 17       2 yes   P1           7.05 passed al~        6.5
 
+---
+
+    Code
+      bistro(bm_evid1, human_profiles, pop_allele_freqs = pop_allele_freqs, kit = "ESX17",
+        peak_thresh = 200, return_lrs = TRUE)
+    Message <simpleMessage>
+      1/17 markers in kit but not in pop_allele_freqs: AMEL
+      Formatting bloodmeal profiles
+      Removing 4 peaks under the threshold of 200 RFU.
+      Formatting human profiles
+      Calculating log10LRs
+      # bloodmeal ids: 1
+      # human ids: 3
+      Bloodmeal id 1/1
+      Human id 1/3
+      Human id 2/3
+      Human id 3/3
+      Identifying matches
+    Output
+      $matches
+      # A tibble: 2 x 8
+        bloodmeal_id locus_count est_noc match human_id log10_lr notes      thresh_low
+        <chr>              <int>   <dbl> <chr> <chr>       <dbl> <chr>           <dbl>
+      1 evid1                 17       2 yes   P1           21.8 passed al~        9.5
+      2 evid1                 17       2 yes   P2           10.3 passed al~        9.5
+      
+      $lrs
+      # A tibble: 3 x 7
+        bloodmeal_id human_id               locus_count est_noc efm_noc log10_lr notes
+        <chr>        <chr>                        <int>   <dbl>   <dbl>    <dbl> <lgl>
+      1 evid1        00-JP0001-14_20142342~          17       2       2    -9.62 NA   
+      2 evid1        P1                              17       2       2    21.8  NA   
+      3 evid1        P2                              17       2       2    10.3  NA   
+      
+
diff --git a/tests/testthat/test-bistro.R b/tests/testthat/test-bistro.R
@@ -31,4 +31,14 @@ test_that("bistro works", {
     ),
     "If `calc_allele_freqs = FALSE`, then `pop_allele_freqs` is required."
   )
+
+  expect_snapshot(bistro(
+    bm_evid1,
+    human_profiles,
+    pop_allele_freqs = pop_allele_freqs,
+    kit = "ESX17",
+    peak_thresh = 200,
+    return_lrs = TRUE
+  ))
+
 })