Feature filter multiply labeled peptides

R/clean_DIANN.R

@@ -243,6 +243,12 @@
 #' parenthetical annotation is stripped out of \code{PeptideSequence},
 #' leaving the plain amino acid sequence.
 #'
+#' Peptides with two or more labelable residues are dropped, counting across
+#' all of \code{labeledAminoAcids} combined, since partial labeling is not
+#' supported by the turnover model.  This applies only to the
+#' \code{ModifiedSequence} path; channel-based labeling is not inferred from
+#' sequence content and is left untouched.
+#'
 #' @param dn_input \code{data.table} after column renaming.
 #' @param labeledAminoAcids Character vector of single-letter amino acid codes
 #'   (e.g. \code{c("K")} or \code{c("K", "R")}), or \code{NULL} to skip
@@ -274,6 +280,8 @@
         light_regex <- paste0("\\([^-]+-(?:", aa_pattern, ")-L\\)")
         strip_regex <- paste0("\\([^-]+-(?:", aa_pattern, ")-[HL]\\)")
 
+        dn_input <- .filterMultiplyLabeledPeptides(dn_input, aa_pattern,
+                                                   "\\([^)]*\\)")
         dn_input <- .classifyIsotopeLabelType(dn_input, heavy_regex, light_regex)
         dn_input[, PeptideSequence := gsub(strip_regex, "", PeptideSequence, perl = TRUE)]
     }

R/clean_Spectronaut.R

@@ -179,19 +179,27 @@
 #'
 #' In Spectronaut protein turnover reports, heavy peptides appear in
 #' \code{FG.LabeledSequence} with a bracketed modification, e.g.
-#' \code{_PEPTIDEK[Lys6]_}.  Any sequence that contains
-#' \code{[<heavyLabel>]} is classified as heavy; all others are light.
-#' Sequences that do not have amino acids that can carry the label
-#' are classified as \code{NA}.  For example, if \code{heavyLabels} is 
-#' \code{"Lys6"}, then \code{PEPTIDEZ} is classified as NA since it
-#' has no lysine residues that could be labeled.
+#' \code{_PEPTIDEK[Lys6]_}.  Each entry of \code{heavyLabels} names the
+#' labelable residue and the label together, as \code{<residue>[<label>]}, so
+#' any sequence containing that tag is classified as heavy.  The residue prefix
+#' is what identifies which unlabeled sequences are light: a sequence with no
+#' labeled tag but at least one labelable residue is light, and one with no
+#' labelable residue at all is \code{NA}.  For example, if \code{heavyLabels}
+#' is \code{"K[Lys6]"}, then \code{PEPTIDEK} is classified as light and
+#' \code{PEPTIDEZ} as NA, since it has no lysine residues that could be
+#' labeled.
+#'
+#' Peptides with two or more labelable residues are dropped, counting across
+#' all residues named in \code{heavyLabels} combined, since partial labeling
+#' is not supported by the turnover model.
 #' When \code{heavyLabel} is \code{NULL} the column is left untouched so
 #' that the downstream \code{columns_to_fill} default of \code{"L"} applies,
 #' preserving backwards compatibility.
 #'
 #' @param spec_input `data.table` after column renaming.
-#' @param heavyLabels Character scalar heavy label name (e.g. \code{"Lys6"}),
-#'   or \code{NULL}.
+#' @param heavyLabels Character vector of labelable residue and heavy label
+#'   pairs, written as \code{<residue>[<label>]} (e.g. \code{"K[Lys6]"}), or
+#'   \code{NULL}.
 #' @return `data.table` with \code{IsotopeLabelType} column added or updated.
 #' @keywords internal
 #' @noRd
@@ -208,6 +216,8 @@
         collapse = "|"
     )
 
+    spec_input = .filterMultiplyLabeledPeptides(spec_input, labeled_aa_regex,
+                                                "\\[.*?\\]")
     spec_input = .classifyIsotopeLabelType(spec_input, heavy_regex,
                                             labeled_aa_regex = labeled_aa_regex)
 

R/converters_DIANNtoMSstatsFormat.R

@@ -48,6 +48,17 @@
 #' parenthetical annotation is stripped out of \code{PeptideSequence},
 #' leaving the plain amino acid sequence.
 #'
+#' In this path only, peptides carrying more than one labelable residue are
+#' removed, and the number removed is reported.  Residues are counted across
+#' all of \code{labeledAminoAcids} together, so with \code{c("K", "R")} a
+#' peptide containing one lysine and one arginine counts as two and is removed.
+#' Such peptides can be only partially labeled, producing more than the two
+#' mass states (fully light and fully heavy) that the turnover model
+#' represents.  Supporting them is future work; turnover analysis is currently
+#' limited to peptides with exactly one labelable residue.  The channel-based
+#' path above is unaffected, since it does not infer labeling from sequence
+#' content.
+#'
 #' When \code{NULL} (default), protein-turnover mode is disabled and all
 #' peptides receive \code{IsotopeLabelType = "Light"}.
 #' @param quantificationColumn Use 'FragmentQuantCorrected'(default) column for quantified intensities for DIANN 1.8.x.

R/converters_SpectronauttoMSstatsFormat.R

@@ -12,17 +12,30 @@
 #' @param peptideSequenceColumn Name of the Spectronaut column that contains the
 #' peptide sequence.  Defaults to \code{"EG.ModifiedSequence"}. The value is 
 #' standardized internally (dots and spaces removed) before column lookup.
-#' @param heavyLabels Character list identifying the heavy isotope labels as it
-#'   appears inside square brackets in the peptide sequence column, e.g.
-#'   \code{c("Lys6")} matches peptides containing \code{[Lys6]}.  
-#'   \code{c("Lys6", "Arg10")} matches peptides containing either \code{[Lys6]} or \code{[Arg10]}.
-#'   Supports any novel label name reported by Spectronaut (e.g. \code{"Leu6"},
-#'   \code{"Phe10"}).  When provided, peptides are
+#' @param heavyLabels Character vector identifying the labelable residue and the
+#'   heavy isotope label as they appear in the peptide sequence column, written
+#'   as the single-letter amino acid code followed by the label in square
+#'   brackets, e.g. \code{c("K[Lys6]")} matches peptides containing
+#'   \code{K[Lys6]}.  \code{c("K[Lys6]", "R[Arg10]")} matches peptides
+#'   containing either \code{K[Lys6]} or \code{R[Arg10]}.  The residue prefix
+#'   is required: it identifies which residues could carry the label, and so
+#'   which unlabeled peptides are light rather than unlabelable.  Supports any
+#'   novel label name reported by Spectronaut (e.g. \code{"L[Leu6]"},
+#'   \code{"F[Phe10]"}).  When provided, peptides are
 #'   classified as heavy (\code{IsotopeLabelType = "H"}), light
 #'   (\code{IsotopeLabelType = "L"}), or unlabeled
 #'   (\code{IsotopeLabelType = NA}) based on its labeled sequence.  When
 #'   \code{NULL} (default) all peptides receive \code{IsotopeLabelType = "L"}.
 #'   Useful for protein turnover experiments.
+#'
+#'   Peptides carrying more than one labelable residue are removed, and the
+#'   number removed is reported.  Residues are counted across all labels
+#'   supplied together, so with \code{c("K[Lys6]", "R[Arg10]")} a peptide containing
+#'   one lysine and one arginine counts as two and is removed.  Such peptides
+#'   can be only partially labeled, producing more than the two mass states
+#'   (fully light and fully heavy) that the turnover model represents.
+#'   Supporting them is future work; turnover analysis is currently limited to
+#'   peptides with exactly one labelable residue.
 #' @param excludedFromQuantificationFilter Remove rows with F.ExcludedFromQuantification=TRUE Default is TRUE.
 #' @param filter_with_Qvalue FALSE(default) will not perform any filtering. TRUE will filter out the intensities that have greater than qvalue_cutoff in EG.Qvalue column. Those intensities will be replaced with zero and will be considered as censored missing values for imputation purpose.
 #' @param qvalue_cutoff Cutoff for EG.Qvalue. default is 0.01.

R/utils_clean_features.R

@@ -309,6 +309,53 @@
 }
 
 
+#' Count regex matches per element, scoring no-match and \code{NA} as 0.
+#' @param x Character vector to search.
+#' @param pattern Perl-compatible regex.
+#' @return Integer vector of match counts, the same length as \code{x}.
+#' @keywords internal
+#' @noRd
+.countRegexMatches = function(x, pattern) {
+    lengths(regmatches(x, gregexpr(pattern, x, perl = TRUE)))
+}
+
+
+#' Drop peptides carrying more than one labelable residue.
+#'
+#' Such peptides can be partially labeled, which the two-state turnover model
+#' cannot represent, so heavy and light rows are dropped together to keep the
+#' light/heavy ratio unbiased.  The number of peptides removed is logged, since
+#' the exclusion is otherwise invisible to the user.
+#'
+#' @param dt \code{data.table} with a \code{PeptideSequence} column.
+#' @param residue_regex Perl-compatible regex matching one labelable residue.
+#' @param strip_regex Perl-compatible regex matching label and modification
+#'   annotations, removed before counting so that residue letters inside an
+#'   annotation are not counted.
+#' @return \code{dt} with multiply labeled rows removed.
+#' @keywords internal
+#' @noRd
+.filterMultiplyLabeledPeptides = function(dt, residue_regex, strip_regex) {
+    stripped = gsub(strip_regex, "", dt[["PeptideSequence"]], perl = TRUE)
+    n_labelable = .countRegexMatches(stripped, residue_regex)
+    is_multiply_labeled = n_labelable >= 2L
+
+    if (any(is_multiply_labeled)) {
+        # Count distinct peptides on the stripped sequence, so that the heavy
+        # and light forms of one peptide are not reported as two.
+        msg = paste("**", data.table::uniqueN(stripped[is_multiply_labeled]),
+                    "peptide(s) with more than one labelable residue were",
+                    paste0("removed (", sum(is_multiply_labeled), " row(s))."),
+                    "Turnover analysis is currently limited to peptides with",
+                    "exactly one labelable residue.")
+        getOption("MSstatsLog")("INFO", msg)
+        getOption("MSstatsMsg")("INFO", msg)
+    }
+
+    dt[!is_multiply_labeled, ]
+}
+
+
 #' Classify IsotopeLabelType from peptide sequence patterns.
 #'
 #' Shared core logic for protein turnover workflows in both Spectronaut and

inst/tinytest/test_clean_DIANN.R

@@ -99,3 +99,54 @@ expect_equal(result_multi_aa$IsotopeLabelType,
              c("H", "H", "L", "L", NA_character_))
 expect_equal(sort(unique(result_multi_aa$PeptideSequence)),
              c("PEPTIDEAC", "PEPTIDEK", "PEPTIDER"))
+
+# Multiply labeled peptides (2+ labelable residues) are filtered out
+dt_multi_label = data.table::data.table(
+    PeptideSequence = c(
+        "PEPTIDEK(SILAC-K-H)",                    # 1 K, heavy -> kept
+        "PEPTIDEK(SILAC-K-L)",                    # 1 K, light -> kept
+        "PEPK(SILAC-K-H)TIDEK(SILAC-K-H)",        # 2 K, heavy -> dropped
+        "PEPK(SILAC-K-H)TIDEK(SILAC-K-L)",        # 2 K, partial -> dropped
+        "PEPK(SILAC-K-L)TIDEK(SILAC-K-L)",        # 2 K, light -> dropped
+        "PEPTIDEAC"                               # 0 K -> kept as NA
+    )
+)
+result_multi_label = MSstatsConvert:::.assignDIANNIsotopeLabelType(
+    dt_multi_label, labeledAminoAcids = c("K"), has_channel = FALSE
+)
+expect_equal(result_multi_label$PeptideSequence,
+             c("PEPTIDEK", "PEPTIDEK", "PEPTIDEAC"))
+expect_equal(result_multi_label$IsotopeLabelType, c("H", "L", NA_character_))
+
+# Count is taken across all labeled amino acids combined
+dt_kr = data.table::data.table(
+    PeptideSequence = c("PEPTIDEK(SILAC-K-H)",
+                        "PEPK(SILAC-K-H)TIDER(SILAC-R-H)",  # 1 K + 1 R -> dropped
+                        "PEPTIDER(SILAC-R-L)")
+)
+result_kr = MSstatsConvert:::.assignDIANNIsotopeLabelType(
+    dt_kr, labeledAminoAcids = c("K", "R"), has_channel = FALSE
+)
+expect_equal(result_kr$PeptideSequence, c("PEPTIDEK", "PEPTIDER"))
+expect_equal(result_kr$IsotopeLabelType, c("H", "L"))
+
+# Residue letters inside an unrelated modification are not counted
+dt_mod = data.table::data.table(
+    PeptideSequence = c("PEPTS(Kmod)IDEK(SILAC-K-H)")
+)
+result_mod = MSstatsConvert:::.assignDIANNIsotopeLabelType(
+    dt_mod, labeledAminoAcids = c("K"), has_channel = FALSE
+)
+expect_equal(nrow(result_mod), 1L)
+expect_equal(result_mod$IsotopeLabelType, "H")
+
+# The Channel path is exempt: labeling is not inferred from sequence content
+dt_channel_multi = data.table::data.table(
+    PeptideSequence = c("PEPKTIDEK", "PEPKTIDEK", "PEPTIDEK"),
+    Channel = c("H", "L", "H")
+)
+result_channel_multi = MSstatsConvert:::.assignDIANNIsotopeLabelType(
+    dt_channel_multi, labeledAminoAcids = c("K"), has_channel = TRUE
+)
+expect_equal(nrow(result_channel_multi), 3L)
+expect_equal(result_channel_multi$IsotopeLabelType, c("H", "L", "H"))

inst/tinytest/test_clean_Spectronaut.R

@@ -73,3 +73,37 @@ expect_equal(result$PeptideSequence,
 dt = make_spec_input(c("_PEPTIDEK_", "_PEPTIDER_"))
 result = MSstatsConvert:::.assignSpectronautIsotopeLabelType(dt, heavyLabels = NULL)
 expect_equal(result, dt)
+
+# Multiply labeled peptides (2+ labelable residues) are filtered out
+dt = make_spec_input(c(
+    "_PEPTIDEK[Lys6]_",         # 1 K, heavy -> kept
+    "_PEPTIDEK_",               # 1 K, light -> kept
+    "_PEPK[Lys6]TIDEK[Lys6]_",  # 2 K, fully heavy -> dropped
+    "_PEPK[Lys6]TIDEK_",        # 2 K, partially labeled -> dropped
+    "_PEPKTIDEK_",              # 2 K, fully light -> dropped
+    "_ACDEGFHI_"                # 0 K -> kept as NA
+))
+result = MSstatsConvert:::.assignSpectronautIsotopeLabelType(
+    dt, heavyLabels = "K[Lys6]")
+expect_equal(result$PeptideSequence,
+             c("_PEPTIDEK_", "_PEPTIDEK_", "_ACDEGFHI_"))
+expect_equal(result$IsotopeLabelType, c("H", "L", NA_character_))
+
+# Count is taken across all labelable residues combined: one K plus one R is
+# doubly labelable when both labels are specified
+dt = make_spec_input(c(
+    "_PEPTIDEK_",   # 1 labelable -> kept
+    "_PEPTIDEKR_",  # 1 K + 1 R -> dropped
+    "_PEPTIDER_"    # 1 labelable -> kept
+))
+result = MSstatsConvert:::.assignSpectronautIsotopeLabelType(
+    dt, heavyLabels = c("K[Lys6]", "R[Arg10]"))
+expect_equal(result$PeptideSequence, c("_PEPTIDEK_", "_PEPTIDER_"))
+expect_equal(result$IsotopeLabelType, c("L", "L"))
+
+# Residue letters inside an unrelated modification tag are not counted
+dt = make_spec_input(c("_S[Kmodification]PEPTIDEK_"))
+result = MSstatsConvert:::.assignSpectronautIsotopeLabelType(
+    dt, heavyLabels = "K[Lys6]")
+expect_equal(nrow(result), 1L)
+expect_equal(result$IsotopeLabelType, "L")

inst/tinytest/test_utils_clean_features.R

@@ -216,6 +216,32 @@ result_multi = MSstatsConvert:::.classifyIsotopeLabelType(
 )
 expect_equal(result_multi$IsotopeLabelType, c("H", "H", "L", "L", NA_character_))
 
+# Test .filterMultiplyLabeledPeptides ----
+# Heavy and light forms of one peptide are reported as a single peptide, and
+# the row count is reported alongside it
+dt_report = data.table::data.table(PeptideSequence = c(
+    "_PEPKTIDEK[Lys6]_", "_PEPKTIDEK[Lys6]_", "_PEPKTIDEK_",  # 1 peptide, 3 rows
+    "_PEPTIDEK_"                                              # kept
+))
+# The appender writes with cat(), so the console output is captured directly
+log_report = capture.output(
+    kept <- MSstatsConvert:::.filterMultiplyLabeledPeptides(
+        dt_report, "K", "\\[.*?\\]")
+)
+expect_true(any(grepl(
+    "1 peptide(s) with more than one labelable residue were removed (3 row(s))",
+    log_report, fixed = TRUE)))
+expect_equal(nrow(kept), 1L)
+
+# Nothing to remove means nothing is reported
+dt_quiet = data.table::data.table(PeptideSequence = c("_PEPTIDEK_", "_ACDEG_"))
+log_quiet = capture.output(
+    all_kept <- MSstatsConvert:::.filterMultiplyLabeledPeptides(
+        dt_quiet, "K", "\\[.*?\\]")
+)
+expect_false(any(grepl("labelable residue", log_quiet, fixed = TRUE)))
+expect_equal(nrow(all_kept), 2L)
+
 # Utility function ----
 expect_equal(MSstatsConvert:::.combine(c("A", "B"), c("A", "B")),
              c("A_A", "B_B"))