Enable triple SILAC (#51)

NAMESPACE

@@ -24,6 +24,7 @@ export(get_medians)
 export(get_notch_per_protein)
 export(get_parsimony_pep2prot)
 export(get_psm_metrics)
+export(get_psm_silac_mod_regex)
 export(get_ratio)
 export(get_sequence)
 export(make_fasta)

R/get_ratio.R

@@ -7,6 +7,8 @@
 #'
 #' @param data `data.frame`.
 #' @param treatment,control `variable`. Name of numeric column in the data.frame.
+#' @param ratio_colname `character`. Name for the output ratio column
+#' @param missing_colname `character`. Name for the output missing column
 #' @param bind `logical`. Should the resulting ratios be added to data
 #' as a column? Default is `TRUE`.
 #'
@@ -21,7 +23,10 @@
 #'
 #' my_data2 <- get_ratio(my_data, treatment, control)
 #' @export
-get_ratio <- function(data, treatment, control, bind = TRUE) {
+get_ratio <- function(data, treatment, control,
+                      ratio_colname = 'ratio',
+                      missing_colname = 'missing',
+                      bind = TRUE) {
   treatment <- deparse(substitute(treatment))
   control <- deparse(substitute(control))
 
@@ -44,10 +49,13 @@ get_ratio <- function(data, treatment, control, bind = TRUE) {
       c(ratio, missing)
     }
   )
-  result <- list(
-    ratio = as.numeric(ratios[1, ]),
-    missing = as.factor(ratios[2, ])
+
+  result <- list(as.numeric(ratios[1, ]),
+                 as.factor(ratios[2, ])
   )
+
+  names(result) <- c(ratio_colname, missing_colname)
+
   if (bind) result <- cbind(data, result)
   result
 }

R/modifications.R

@@ -48,16 +48,26 @@ psm_to_peptide_style_modifications <- function(psm_style_modifications) {
 #'
 #' @param mod_col `character vector` Modification column from Proteome Discoverer
 #' @param level `character` Either 'psm' or 'peptide'
+#' @param psm_modfication_regexes `character vector` One or more regexes to match the expected SILAC modifications.
+#' See \link[camprotR]{get_psm_silac_mod_regex}
 #' @return `character vector` updated Modifications column
 #' @export
-remove_silac_modifications <- function(mod_col, level = 'psm') {
+remove_silac_modifications <- function(mod_col, level = 'psm',
+                                       psm_modfication_regexes=c(get_psm_silac_mod_regex('R_13C6_15N4'),
+                                                                 get_psm_silac_mod_regex('K_13C6_15N2'))){
   if (!level %in% c('psm', 'peptide')) stop('level must be psm or peptide')
 
   if (level == 'psm') {
+
+    psm_modfication_regex <- sprintf(
+      '(; )?(%s)',
+      paste(psm_modfication_regexes, collapse='|')
+      )
+
     mod_col <- gsub(
       '^(; )', '',
       gsub(
-        '(; )?(K|R)\\d{1,2}\\(Label:13C\\(6\\)15N\\((2|4)\\)\\)',
+        psm_modfication_regex,
         '', mod_col
       )
     )
@@ -73,3 +83,45 @@ remove_silac_modifications <- function(mod_col, level = 'psm') {
 
   return(mod_col)
 }
+
+#' Get a pre-defined regex for a SILAC modification in PSM format
+#'
+#' @description This function returns a regex which can be used to remove SILAC
+#' modifications from the modification column in the PSM-level PD output using
+#' \link[camprotR]{remove_silac_modifications}. Call without any arguments to see
+#' a description of the available modifications.
+#'
+#' @param silac_mod SILAC modification name (call this function without
+#' arguments to see available values).
+#' @return Returns a `character vector`, regex for SILAC modification.
+#' @export
+get_psm_silac_mod_regex <- function(silac_mod){
+
+  regexes <- list('R_13C6_15N4' =
+                    list('desc'='Heavy R (10)',
+                         'regex'='R\\d{1,2}\\(Label:13C\\(6\\)15N\\(4\\)\\)'),
+                  'K_13C6_15N2' =
+                     list('desc'='Heavy K (8)',
+                          'regex'='K\\d{1,2}\\(Label:13C\\(6\\)15N\\(2\\)\\)'),
+                  'R_13C6' =
+                    list('desc'='Median R (6)',
+                         'regex'='R\\d{1,2}\\(Label:13C\\(6\\)\\)'),
+                  'K_2H4' =
+                    list('desc'='Median K (4)',
+                         'regex'='K\\d{1,2}\\(Label:2H\\(4\\)\\)'))
+
+  if(missing(silac_mod)){
+
+    table <- do.call('rbind.data.frame', regexes) %>%
+      data.frame() %>%
+      tibble::rownames_to_column('name')
+
+    return(table)
+  }
+  else{
+    silac_mod = match.arg(silac_mod, choices=names(regexes))
+    return(regexes[[silac_mod]][['regex']])
+  }
+}
+
+

R/silac_psm_seq_int.R

@@ -13,7 +13,8 @@
 #' @param include_interference `logical` Should PSM interference be included too?
 #' @param interference_col `string` Column with interference/co-isolation
 #' @param group_cols `string` Additional feature columns to retain, beyond
-#' Sequence and Modification
+#' @param psm_modfication_regexes `character vector` One or more regexes to match the expected SILAC modifications
+#' Sequence and Modification. See \link[camprotR]{get_psm_silac_mod_regex}
 #' @return `data.frame` indicating which SILAC peptides were MS2 matched,
 #' how many PSMs per isotope, and
 #' (optionally) the maximum interference across all PSMs for the peptide
@@ -24,7 +25,9 @@ silac_psm_seq_int <- function(
   mod_col='Modifications',
   include_interference=FALSE,
   interference_col='Isolation.Interference.in.Percent',
-  group_cols=NULL){
+  group_cols=NULL,
+  psm_modfication_regexes=c(get_psm_silac_mod_regex('R_13C6_15N4'),
+                            get_psm_silac_mod_regex('K_13C6_15N2'))){
 
   message('camprotR::silac_psm_seq_int output has changed.
   Columns indicating whether quantification is from PSM are now prefixed with
@@ -47,7 +50,9 @@ silac_psm_seq_int <- function(
     rowwise() %>%
     filter(is.finite(.data$Precursor.Abundance))
 
-  obj[[mod_col]] <- remove_silac_modifications(obj[[mod_col]])
+  obj[[mod_col]] <- remove_silac_modifications(
+    obj[[mod_col]], psm_modfication_regexes=psm_modfication_regexes)
+
   obj[[mod_col]] <- psm_to_peptide_style_modifications(obj[[mod_col]])
 
   obj[[sequence_col]] <- toupper(obj[[sequence_col]])
@@ -69,6 +74,26 @@ silac_psm_seq_int <- function(
            matched_Light=replace_na(.data$matched_Light, FALSE),
            matched_Heavy=replace_na(.data$matched_Heavy, FALSE))
 
+  # Depending on which SILAC modifications in Quan.Channel column,
+  # the output columns need to be updated
+  if('n_Light' %in% colnames(obj_seq)){
+    obj_seq <- obj_seq  %>%
+      mutate(n_Light = replace_na(.data$n_Light, 0),
+             matched_Light = replace_na(.data$matched_Light, FALSE))
+  }
+
+  if('n_Medium' %in% colnames(obj_seq)){
+    obj_seq <- obj_seq  %>%
+      mutate(n_Medium = replace_na(.data$n_Medium, 0),
+             matched_Medium = replace_na(.data$matched_Medium, FALSE))
+  }
+
+  if('n_Heavy' %in% colnames(obj_seq)){
+    obj_seq <- obj_seq  %>%
+      mutate(n_Heavy = replace_na(.data$n_Heavy, 0),
+             matched_Heavy = replace_na(.data$matched_Heavy, FALSE))
+  }
+
   if(include_interference){
 
     if (!interference_col %in% colnames(obj)) {

tests/testthat/test-modifications.R

@@ -18,6 +18,19 @@ test_that("Peptide style SILAC modifications can be removed", {
     "1xCarbamidomethyl [C9]")
 })
 
+test_that("get_psm_silac_mod_regex works with no args", {
+  df <- get_psm_silac_mod_regex()
+
+  expect_true(inherits(df, "data.frame"))
+  expect_equal(colnames(df), c("name", "desc", "regex"))
+})
+
+test_that("get_psm_silac_mod_regex works with args", {
+  regex <- get_psm_silac_mod_regex("K_2H4")
+
+  expect_equal(regex, 'K\\d{1,2}\\(Label:2H\\(4\\)\\)')
+})
+
 #### Sanity checks -------------------------------------------------------------
 test_that("remove_silac_modifications() errors if level is nonsense", {
   expect_error(remove_silac_modifications(