Feature filter multiply labeled peptides

R/clean_DIANN.R

@@ -221,32 +221,24 @@
 
 #' Assign IsotopeLabelType for DIANN protein turnover workflows.
 #'
-#' Dispatches to one of two classification paths depending on \code{has_channel}:
+#' Two paths, chosen by \code{has_channel}.  With a \code{Channel} column,
+#' \code{Channel} maps straight to \code{IsotopeLabelType} (\code{"H"} and
+#' \code{"L"} pass through, anything else becomes \code{NA}) and is then
+#' dropped; \code{labeledAminoAcids} acts only as the opt-in flag, and
+#' sequences are not inspected.
 #'
-#' \strong{Channel-based path} (\code{has_channel = TRUE}): \code{Channel}
-#' values are mapped directly to \code{IsotopeLabelType} (\code{"H"} →
-#' \code{"H"}, \code{"L"} → \code{"L"}, anything else → \code{NA}), and the
-#' \code{Channel} column is then dropped.  \code{labeledAminoAcids} acts solely
-#' as the opt-in flag that enables this path; the amino acid codes are
-#' \strong{not} used to validate or filter \code{ModifiedSequence}.
-#'
-#' \strong{ModifiedSequence-parsing path} (\code{has_channel = FALSE}):
-#' \code{PeptideSequence} (the retained \code{ModifiedSequence}) is scanned
-#' for isotope-labeled amino acids, which appear in parentheses immediately
-#' after the labeled residue, in the form \code{(<label>-<aminoAcid>-H)} for
-#' the heavy label or \code{(<label>-<aminoAcid>-L)} for the light label.
-#' For example, \code{K(label-K-H)} marks a heavy-labeled lysine (\code{K}).
-#' Here \code{<aminoAcid>} is one of the single-letter codes in
-#' \code{labeledAminoAcids}, and \code{<label>} is the label name (e.g.
-#' \code{SILAC} or \code{label}). Sequences with no such parenthetical
-#' tags are assigned \code{IsotopeLabelType = NA}. Once classified, the
-#' parenthetical annotation is stripped out of \code{PeptideSequence},
-#' leaving the plain amino acid sequence.
+#' Without one, the label is read from \code{PeptideSequence} (the retained
+#' \code{ModifiedSequence}), where it follows the labeled residue in
+#' parentheses as \code{(<label>-<aminoAcid>-H)} or
+#' \code{(<label>-<aminoAcid>-L)}; \code{K(label-K-H)} is a heavy-labeled
+#' lysine.  Untagged sequences become \code{NA}, the tags are then stripped
+#' from \code{PeptideSequence}, and peptides with two or more labelable
+#' residues (counted across all of \code{labeledAminoAcids}) are dropped, since
+#' partial labeling breaks the two-state turnover model.
 #'
 #' @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
-#'   protein-turnover mode entirely (backwards-compatible default).
+#'   (e.g. \code{c("K", "R")}), or \code{NULL} to skip protein-turnover mode.
 #' @param has_channel Logical; \code{TRUE} when the raw input contained a
 #'   \code{Channel} column that was retained through
 #'   \code{.cleanDIANNSelectRequiredColumns}.
@@ -274,6 +266,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

@@ -171,27 +171,23 @@
 }
 
 
-#' Assign IsotopeLabelType based on heavy label detection.
+#' Assign IsotopeLabelType from Spectronaut heavy label tags.
 #'
-#' When \code{heavyLabel} is provided, each row is classified as heavy
-#' (\code{"H"}), light (\code{"L"}), or unlabeled (\code{NA}) by inspecting
-#' the labeled sequence column for the presence of the label tag.
+#' In Spectronaut turnover reports a labeled residue carries its label in
+#' brackets, e.g. \code{_PEPTIDEK[Lys6]_}.  Each \code{heavyLabels} entry pairs
+#' the residue with the label (\code{"K[Lys6]"}), so the bare residue marks a
+#' peptide that could have been labeled: with \code{"K[Lys6]"},
+#' \code{PEPTIDEK[Lys6]} is heavy, \code{PEPTIDEK} is light, and
+#' \code{PEPTIDEZ} is \code{NA}.  Peptides with two or more labelable residues
+#' (counted across all labels) are dropped, since partial labeling breaks the
+#' two-state turnover model.
 #'
-#' 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.
-#' 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.
+#' With \code{NULL} the column is left untouched so the downstream
+#' \code{columns_to_fill} default of \code{"L"} applies.
 #'
 #' @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 \code{<residue>[<label>]} pairs
+#'   (e.g. \code{"K[Lys6]"}), or \code{NULL}.
 #' @return `data.table` with \code{IsotopeLabelType} column added or updated.
 #' @keywords internal
 #' @noRd
@@ -208,6 +204,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

@@ -19,37 +19,28 @@
 #' @param removeFewMeasurements should proteins with few measurements be removed
 #' @param removeOxidationMpeptides should peptides with oxidation be removed
 #' @param removeProtein_with1Feature should proteins with a single feature be removed
-#' @param labeledAminoAcids Character vector of single-letter amino acid codes
-#' that carry the SILAC label in protein turnover experiments, e.g.
-#' \code{c("K")} or \code{c("K", "R")}.  Supplying this vector opts in to
-#' protein-turnover mode; the exact amino acids determine behaviour only in the
-#' \code{ModifiedSequence}-parsing path described below.
+#' @param labeledAminoAcids Character vector of single-letter codes for the
+#' amino acids carrying the SILAC label, e.g. \code{c("K")} or
+#' \code{c("K", "R")}.  Supplying it turns on protein-turnover mode.  Defaults
+#' to \code{NULL}, which marks every peptide
+#' \code{IsotopeLabelType = "Light"}.  How the label is read depends on the
+#' export:
 #'
-#' \strong{Channel-based path} (DIA-NN 2.x exports that include a
-#' \code{Channel} column): when \code{labeledAminoAcids} is non-\code{NULL}
-#' \emph{and} the input contains a \code{Channel} column, \code{Channel} values
-#' are mapped directly to \code{IsotopeLabelType} (\code{"H"} → \code{"H"},
-#' \code{"L"} → \code{"L"}, anything else → \code{NA}).  The amino acid codes
-#' in \code{labeledAminoAcids} are \strong{not} used to validate or filter
-#' \code{ModifiedSequence} in this path.
+#' \strong{DIA-NN 2.x, with a \code{Channel} column.}  \code{Channel} gives the
+#' label directly: \code{"H"} and \code{"L"} pass through, anything else
+#' becomes \code{NA}.  The amino acid codes only switch the mode on here;
+#' sequences are not inspected.
 #'
-#' \strong{ModifiedSequence-parsing path} (DIA-NN 1.x exports without a
-#' \code{Channel} column): when \code{labeledAminoAcids} is non-\code{NULL}
-#' and no \code{Channel} column is present, each \code{ModifiedSequence} is 
-#' scanned for isotope-labeled amino acids, which appear in parentheses 
-#' immediately after the labeled residue, 
-#' in the form \code{(<label>-<aminoAcid>-H)} for
-#' the heavy label or \code{(<label>-<aminoAcid>-L)} for the light label.
-#' For example, \code{K(label-K-H)} marks a heavy-labeled lysine (\code{K}).
-#' Here \code{<aminoAcid>} is one of the single-letter codes in
-#' \code{labeledAminoAcids}, and \code{<label>} is the label name (e.g.
-#' \code{SILAC} or \code{label}). Sequences with no such parenthetical
-#' tags are assigned \code{IsotopeLabelType = NA}. Once classified, the
-#' parenthetical annotation is stripped out of \code{PeptideSequence},
-#' leaving the plain amino acid sequence.
-#'
-#' When \code{NULL} (default), protein-turnover mode is disabled and all
-#' peptides receive \code{IsotopeLabelType = "Light"}.
+#' \strong{DIA-NN 1.x, without a \code{Channel} column.}  The label is read
+#' from \code{ModifiedSequence}, where it follows the labeled residue in
+#' parentheses as \code{(<label>-<aminoAcid>-H)} or
+#' \code{(<label>-<aminoAcid>-L)}; for example \code{K(label-K-H)} is a
+#' heavy-labeled lysine.  \code{<label>} is the label name DIA-NN reported
+#' (e.g. \code{SILAC}).  Untagged sequences become \code{NA}, and the tags are
+#' stripped from \code{PeptideSequence} afterwards.  A peptide with two or more
+#' labelable residues (counted across all of \code{labeledAminoAcids}) may be
+#' only partially labeled, which the two-state turnover model cannot represent,
+#' so it is dropped and the number removed is reported.
 #' @param quantificationColumn Use 'FragmentQuantCorrected'(default) column for quantified intensities for DIANN 1.8.x.
 #' Use 'FragmentQuantRaw' for quantified intensities for DIANN 1.9.x.
 #' Use 'auto' for quantified intensities for DIANN 2.x where each fragment intensity is a separate column, e.g. Fr0Quantity.

R/converters_SpectronauttoMSstatsFormat.R

@@ -12,17 +12,23 @@
 #' @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
-#'   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.
+#' @param heavyLabels Character vector naming each labeled residue and its
+#'   heavy label as they appear in the peptide sequence column: the
+#'   single-letter amino acid code, then the label in square brackets.  For
+#'   example \code{"K[Lys6]"}, or \code{c("K[Lys6]", "R[Arg10]")} for a
+#'   double-label experiment.  Any label name Spectronaut reports is accepted.
+#'   The residue letter is required, because it is what tells MSstats which
+#'   peptides could have carried a label.
+#'
+#'   Peptides carrying the label are marked heavy
+#'   (\code{IsotopeLabelType = "H"}), those that could carry it but do not are
+#'   light (\code{"L"}), and those with no labelable residue are \code{NA}.  A
+#'   peptide with two or more labelable residues (counted across all labels
+#'   supplied) may be only partially labeled, which the two-state turnover
+#'   model cannot represent, so it is dropped and the number removed is
+#'   reported.
+#'
+#'   Defaults to \code{NULL}: turnover mode off, every peptide marked light.
 #' @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,22 +309,63 @@
 }
 
 
+#' 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
-#' DIANN converters.  Each peptide is classified as heavy (\code{"H"}), light
-#' (\code{"L"}), or unlabeled (\code{NA}) by matching regex patterns against
-#' the \code{PeptideSequence} column.
+#' Shared by the Spectronaut and DIANN turnover workflows.  Each peptide is
+#' classified as heavy (\code{"H"}), light (\code{"L"}), or unlabeled
+#' (\code{NA}) by matching regexes against \code{PeptideSequence}.
 #'
-#' Two modes are supported:
-#' \describe{
-#'   \item{Spectronaut mode}{Pass \code{labeled_aa_regex}.  The sequence is
-#'     first stripped of all bracket modifications; light is inferred when the
-#'     bare labeled amino acid is present but the heavy bracket form is absent.}
-#'   \item{DIANN mode}{Pass \code{light_regex}.  Both heavy and light patterns
-#'     are matched directly against the modified sequence; absence of either
-#'     yields \code{NA}.}
-#' }
+#' Spectronaut mode (\code{labeled_aa_regex}) strips bracket modifications
+#' first, then infers light from a bare labelable residue with no heavy tag.
+#' DIANN mode (\code{light_regex}) matches heavy and light patterns directly
+#' against the modified sequence.
 #'
 #' @param dt \code{data.table} with a \code{PeptideSequence} column.
 #' @param heavy_regex Perl-compatible regex matching heavy-labeled sequences.

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")