'get_aa_transition_matrix' returns either BLOSUM62 or (the untested) …

…FLU. Use EpitopePrediction as the IC50 prediction tool in tests

R/are_detected.R

@@ -4,16 +4,20 @@
 #' @export
 are_detected <- function(
   protein_sequences,
-  mhc_haplotype = "HLA-A*02-01",
-  percentile = get_ic50_percentile_binder()
+  mhc_haplotype,
+  peptide_length,
+  percentile = get_ic50_percentile_binder(),
+  ic50_prediction_tool
 ) {
   bbbq::check_mhc_haplotype_name(mhc_haplotype)
   results <- rep(NA, length(protein_sequences))
   for (i in seq_along(protein_sequences)) {
     results[i] <- bbbq::is_detected(
       protein_sequence = protein_sequences[i],
       mhc_haplotype = mhc_haplotype,
-      percentile = percentile
+      peptide_length = peptide_length,
+      percentile = percentile,
+      ic50_prediction_tool = ic50_prediction_tool
     )
   }
   results

R/create_consensus_topology_conservation.R

@@ -10,17 +10,19 @@
 #' @inheritParams default_params_doc
 #' @export
 create_consensus_topology_conservation <- function( # nolint indeed a long function name
-  protein_sequences
+  protein_sequences,
+  transition_matrix = get_aa_transition_matrix("BLOSUM62")
 ) {
   protein_sequences_aass <- Biostrings::AAStringSet(protein_sequences)
 
   sink("/dev/null")
   protein_alignment <- msa::msa(protein_sequences_aass)
   sink()
 
-  utils::data(BLOSUM62)
   conservation_scores <- msa::msaConservationScore(
-    protein_alignment, BLOSUM62)
+    protein_alignment,
+    transition_matrix
+  )
   t <- tibble::tibble(
     aa = names(conservation_scores),
     score = as.numeric(conservation_scores),

R/get_aa_transition_matrix.R

@@ -15,7 +15,14 @@
 #'   amino acid name
 #' @author Richèl J.C. Bilderbeek
 #' @export
-get_aa_transition_matrix <- function() {
+get_aa_transition_matrix <- function(
+  transition_matrix_name = "FLU"
+) {
+  if (transition_matrix_name == "BLOSUM62") {
+    library(Biostrings, warn.conflicts = FALSE, quietly = TRUE)
+    utils::data(BLOSUM62)
+    return(BLOSUM62)
+  }
   t <- as.matrix(
     readr::read_csv(
       system.file("extdata", "flu_transitions.csv", package = "bbbq"),

R/is_detected.R

@@ -1,30 +1,29 @@
 #' Determine if the haplotype has at least one binder to the peptide
 #' @inheritParams default_params_doc
-#' @param fragment_length the length of the peptide fragments
-#'   the peptide will be cut into
 #' @author Richèl J.C. Bilderbeek
 #' @export
 is_detected <- function(
   protein_sequence,
   mhc_haplotype,
-  fragment_length = 9,
-  percentile = get_ic50_percentile_binder()
+  peptide_length,
+  percentile = get_ic50_percentile_binder(),
+  ic50_prediction_tool
 ) {
   bbbq::check_mhc_haplotype_name(mhc_haplotype)
 
 
-  ic50_threshold <- mhcnpreds::get_ic50_threshold(
-    peptide_length = fragment_length,
-    mhc_haplotype = mhcnuggetsr::to_mhcnuggets_name(mhc_haplotype),
-    percentile = percentile
+  ic50_threshold <- bbbq::get_ic50_threshold(
+    peptide_length = peptide_length,
+    haplotype = mhc_haplotype,
+    percentile = percentile,
+    ic50_prediction_tool = ic50_prediction_tool
   )
 
-  results <- mhcnuggetsr::predict_ic50s(
-    mhcnuggets_options = mhcnuggetsr::create_mhcnuggets_options(
-      mhc = mhcnuggetsr::to_mhcnuggets_name(mhc_haplotype)
-    ),
+  results <- bbbq::predict_ic50s(
     protein_sequence = protein_sequence,
-    peptide_length = fragment_length
+    peptide_length = peptide_length,
+    haplotype = mhc_haplotype,
+    ic50_prediction_tool = ic50_prediction_tool
   )
   sum(results$ic50 < ic50_threshold) > 0
 }

tests/testthat/test-are_detected.R

@@ -1,11 +1,14 @@
 test_that("use", {
-  skip("Not now")
   expect_equal(
     are_detected(
-      c(
+      protein_sequences = c(
         "VVIILTIAGNILVIMAVSLE",
         "AAAARTIAGRILVIMARSAA"
-      )
+      ),
+      mhc_haplotype = bbbq::get_mhc1_haplotypes()[1],
+      peptide_length = 9,
+      percentile = 0.02,
+      ic50_prediction_tool = "EpitopePrediction"
     ),
     c(TRUE, FALSE)
   )

tests/testthat/test-create_consensus_topology_conservation.R

@@ -10,6 +10,7 @@ test_that("use, simple data", {
     "IMPRESSIVELYFLILAWAYFANSWALKWEETMARKETTRIVIALLYNAILIDENTIFY",
     "IMPRESSIVELYFLIYAWAYFANSWALKSWEETMARKETRIVIALLYNAILIDENTIFY"
   )
+  library(msa)
   t <- create_consensus_topology_conservation(protein_sequences)
   expect_true(tibble::is_tibble(t))
   expect_equal(60, nrow(t))

tests/testthat/test-get_aa_transition_matrix.R

@@ -1,4 +1,8 @@
 test_that("use", {
+  # Follow the BioStrings BLOSUM62 matrix
+  #utils::data(BLOSUM62);
+  #expect_equal(rownames(BLOSUM62), rownames(m))
+
   m <- get_aa_transition_matrix()
   expect_equal(20, nrow(m))
   expect_equal(20, ncol(m))

tests/testthat/test-is_detected.R

@@ -1,18 +1,21 @@
 test_that("use", {
 
-  if (!mhcnuggetsr::is_mhcnuggets_installed()) return()
-
   expect_true(
     is_detected(
       "VVIILTIAGN",
-      mhc_haplotype = mhcnuggetsr::get_mhc_1_haplotypes()[4],
-      percentile = 0.5
+      mhc_haplotype = bbbq::get_mhc1_haplotypes()[1],
+      peptide_length = 9,
+      percentile = 0.5,
+      ic50_prediction_tool = "EpitopePrediction"
     )
   )
   expect_false(
     is_detected(
       "AAAARTIAGR",
-      mhc_haplotype = mhcnuggetsr::get_mhc_1_haplotypes()[20]
+      mhc_haplotype = bbbq::get_mhc1_haplotypes()[1],
+      peptide_length = 9,
+      percentile = 0.01,
+      ic50_prediction_tool = "EpitopePrediction"
     )
   )
 })