Merge pull request #1118 from Proteobench/add-mods-to-entrapment

add modification support in entrapment

Author: rodvrees

docs/available-modules/active-modules/13-entrapment-ion-dia-astral.md

@@ -29,7 +29,7 @@ Each workflow submission is classified as:
 
 > **Use the pre-digested entrapment FASTA — do not enable in-silico digestion.** The entrapment FASTA already contains peptide sequences (not full proteins). Enabling digestion in your search engine will lead to identified peptides without matched entrapments, which makes FDP calculation less accurate.
 
-> **Do not add any variable modifications.**
+> **Only methionine oxidation is supported as a variable modification.** The entrapment mapping enumerates every combination of methionine oxidation and cysteine carbamidomethylation. Enabling any other modification (variable, or an additional fixed modification) causes metric calculation to fail with an error, because the target/entrapment pairing cannot be resolved for it.
 
 ## Data set
 
@@ -53,6 +53,10 @@ ProteoBench reads the search engine output, maps runs to samples, and classifies
 
 The three FDP estimates are computed from the resulting set and compared to the reported FDR threshold (inferred from the output file).
 
+### Modification support
+
+Precursor identifications are matched to the entrapment mapping file using their exact modified sequence, not the stripped peptide. The mapping file enumerates every combination of variable methionine oxidation and fixed cysteine carbamidomethylation, so an identification built only from those two modifications is paired correctly, down to the specific set of oxidized methionines. Any other modification present in the reported sequence is not recognized, and metric calculation fails with an error rather than silently ignoring it.
+
 ## How to use
 
 ### Suggested parameters
@@ -67,7 +71,7 @@ The module currently accepts DIA-NN, FragPipe, FragPipe with DIA-NN quantificati
 | Spectral library | Predicted from entrapment FASTA |
 | Digestion | **None** (FASTA is pre-digested) |
 | Fixed modifications | Carbamidomethylation (C) |
-| Variable modifications | **None** |
+| Variable modifications | Methionine oxidation (M), optional. No other variable modification is supported. |
 | Precursor charge range | 1–5 |
 | Precursor m/z range | 400–1000 |
 | Fragment m/z range | 100–1800 |
@@ -93,35 +97,22 @@ You will receive a link to a GitHub pull request. Save it — it contains your r
 1. Import the raw `.raw` files.
 2. Add the entrapment FASTA. Do not enable "Contaminants" — contaminants are already included in the FASTA.
 3. **Disable in-silico digestion.** The FASTA is pre-digested; use '--cut ' in the additional parameter fields to disable enzymatic cleavage.
-4. Enable library-free search / FASTA-based library generation (activates deep-learning prediction of spectra, RTs, and IMs).
-5. Do not set verbosity / log level higher than 1, otherwise parameter parsing will fail.
-6. Upload `report.tsv` or `report.parquet` for metric calculation, and `report.log.txt` for public submission.
+4. Variable modifications are limited to methionine oxidation (optional); disable all other variable modifications. Carbamidomethylation (C) can be used as a fixed modification.
+5. Enable library-free search / FASTA-based library generation (activates deep-learning prediction of spectra, RTs, and IMs).
+6. Do not set verbosity / log level higher than 1, otherwise parameter parsing will fail.
+7. Upload `report.tsv` or `report.parquet` for metric calculation, and `report.log.txt` for public submission.
 
 ### FragPipe with DIA-NN quantification
 
-FragPipe workflows that produce a DIA-NN-style report are submitted as **FragPipe (DIA-NN quant)**. ProteoBench parses the precursor identifications from the DIA-NN report and extracts workflow metadata from the FragPipe `.workflow` file.
-
-1. Use the ProteoBench entrapment FASTA as the sequence database. Do not add a second contaminant database.
-2. Configure the search without enzymatic digestion. The FASTA contains pre-digested peptide entries, so in-silico digestion must remain disabled throughout the FragPipe/DIA-NN workflow. MSFragger Protein Digestion settings:
-
-```
-Load Rules:    nocleavage
-Cuts 1:        @
-No cuts 1:     @
-```
-
-To set these via the GUI: MSFragger tab -> Protein Digestion -> Load Rules = "nocleavage"; Cuts 1 = "@"; No cuts 1 = "@".
-3. Keep variable modifications disabled. Carbamidomethylation (C) may be used as the fixed modification.
-4. Use the DIA-NN report generated by FragPipe (`report.tsv` or `report.parquet`) for metric calculation.
-5. Upload the FragPipe `.workflow` file for public submission. Do not upload the DIA-NN log as the parameter file for this workflow type.
+FragPipe is currently not supported as the decoy generation is not compatible with the predigested module fasta (see [here](https://github.com/Nesvilab/FragPipe/issues/2847)). 
 
 ### [AlphaDIA](https://github.com/MannLabs/alphadia)
 
 AlphaDIA submissions are parsed from precursor-level output. The entrapment module currently expects AlphaDIA 2.x-style precursor output.
 
 1. Use the ProteoBench entrapment FASTA and disable additional contaminants.
 2. Configure AlphaDIA for a no-enzyme / pre-digested FASTA search: set "no-cleave" as the enzyme parameter.
-3. Keep variable modifications disabled. Use Carbamidomethylation (C) as the fixed modification if alkylation was applied.
+3. Variable modifications are limited to methionine oxidation; keep all other variable modifications disabled. Use Carbamidomethylation (C) as the fixed modification if alkylation was applied.
 4. Upload `precursors.parquet` for metric calculation.
 5. Upload the AlphaDIA `log.txt` file for public submission.
 

proteobench/io/parsing/io_parse_settings/entrapment/DIA/ion/Astral/parse_settings_diann.toml

@@ -13,7 +13,7 @@
 "isalpha" = true
 "isupper" = true
 "pattern" = "\\(([^()]*)\\)"
-"modification_dict" = {"(unimod:35)" = "Oxidation", "(unimod:1)" = "Acetyl", "(unimod:4)" = "Carbamidomethyl"}
+"modification_dict" = {"(unimod:35)" = "Oxidation", "(unimod:4)" = "Carbamidomethyl"}
 
 [general]
 "contaminant_flag" = "Cont_"

proteobench/io/parsing/io_parse_settings/parse_settings_files.toml

@@ -112,6 +112,4 @@
 
 [entrapment_DIA_ion_Astral]
 "DIA-NN" = "parse_settings_diann.toml"
-"FragPipe (DIA-NN quant)" = "parse_settings_fragpipe_diann_quant.toml"
-"FragPipe" = "parse_settings_fragpipe_DIA.toml"
 "AlphaDIA" = "parse_settings_alphadia.toml"

proteobench/io/parsing/parse_ion.py

@@ -305,13 +305,17 @@ def get_proforma_bracketed(
         if idx in pos_mod_dict:
             if idx == 0:
                 new_seq += f"[{pos_mod_dict[idx]}]-"
-            elif idx == len(stripped_seq):
-                new_seq += f"-[{pos_mod_dict[idx]}]"
             else:
                 new_seq += f"[{pos_mod_dict[idx]}]"
         if not before_aa:
             new_seq += aa
 
+    # A modification on the last residue (before_aa=False, i.e. bracket after the
+    # residue) is recorded at idx == len(stripped_seq), which the loop above never
+    # reaches (idx only ranges 0..len(stripped_seq)-1), so it must be appended here.
+    if len(stripped_seq) in pos_mod_dict:
+        new_seq += f"[{pos_mod_dict[len(stripped_seq)]}]"
+
     return new_seq
 
 

proteobench/io/parsing/parse_settings.py

@@ -1023,6 +1023,43 @@ def add_modification_parser(self, parser: ParseModificationSettings):
         """
         self.modification_parser = parser
 
+    def _process_modifications(self, df: pd.DataFrame) -> pd.DataFrame:
+        """
+        Standardize the modified sequence column to ProForma-style bracket notation.
+
+        Tools report modifications differently (e.g. DIA-NN/FragPipe use
+        ``M(UniMod:35)``, AlphaDIA already reports ``M[Oxidation]``). The
+        entrapment mapping file is keyed on the bracketed, human-readable form
+        (e.g. ``M[Oxidation]``, ``C[Carbamidomethyl]``), so tools with a
+        ``[modifications_parser]`` section have their modified-sequence column
+        rewritten into that same notation before mapping. Tools without a
+        modification parser (i.e. already reporting that notation natively)
+        are left untouched.
+
+        Parameters
+        ----------
+        df : pd.DataFrame
+            Input DataFrame containing the modified sequence column.
+
+        Returns
+        -------
+        pd.DataFrame
+            DataFrame with the modified sequence column standardized.
+        """
+        if self.modification_parser is None:
+            return df
+
+        mp = self.modification_parser
+        df[mp.modifications_parse_column] = df[mp.modifications_parse_column].apply(
+            get_proforma_bracketed,
+            before_aa=mp.modifications_before_aa,
+            isalpha=mp.modifications_isalpha,
+            isupper=mp.modifications_isupper,
+            pattern=mp.modifications_pattern,
+            modification_dict=mp.modifications_mapper,
+        )
+        return df
+
     def convert_to_standard_format(self, df: pd.DataFrame) -> pd.DataFrame:
         """
         Convert a software tool output into a generic format supported by the module.
@@ -1054,7 +1091,7 @@ def convert_to_standard_format(self, df: pd.DataFrame) -> pd.DataFrame:
         # df = self._fix_colnames(df)
         # df = self._mark_contaminants(df)
         # df = self._process_species_information(df)
-        # df = self._process_modifications(df)
+        df = self._process_modifications(df)
         # df_melted = self._handle_data_format(df)
         # df_melted = self._filter_zero_intensities(df_melted)
         return df  # self._format_by_analysis_level(df_melted), replicate_to_raw

proteobench/modules/entrapment/entrapment_base_module.py

@@ -145,17 +145,31 @@ def _apply_mapping(
         Filter unmapped peptides, assign target/entrapment labels, and merge pair index.
 
         Reads ``self.mapping_file`` (TSV with ``sequence``, ``peptide_pair_index``,
-        and ``peptide_type`` columns).  Peptides absent from the mapping are removed;
-        if more than ``max_missing_fraction`` of unique peptides are absent an
-        ``EntrapmentError`` is raised.  The caller receives a DataFrame with two extra
-        columns: ``"Target or Entrapment"`` (from ``peptide_type``) and
-        ``"peptide_pair_index"``.
+        and ``peptide_type`` columns). The mapping file may contain both bare peptide
+        stems and modification-specific variants (methionine oxidation and/or cysteine
+        carbamidomethylation, in any combination) that share their stem's identity but
+        have their own ``sequence`` and ``peptide_pair_index`` values. Peptides absent
+        from the mapping are removed; if more than ``max_missing_fraction`` of unique
+        peptides are absent an ``EntrapmentError`` is raised.  The caller receives a
+        DataFrame with two extra columns: ``"Target or Entrapment"`` (from
+        ``peptide_type``) and ``"peptide_pair_index"``.
+
+        Matching is done exactly against the modified ``"Sequence"`` column - there is
+        no fallback to the stripped ``"Peptide"`` stem. ``"Sequence"`` is expected to
+        already be standardized to ProForma-style bracket notation (e.g.
+        ``"M[Oxidation]"``, ``"C[Carbamidomethyl]"``) by the parser's
+        ``[modifications_parser]`` step (see ``ParseSettingsEntrapment._process_modifications``),
+        which normalizes each tool's native notation (e.g. DIA-NN's ``"M(UniMod:35)"``)
+        to that common form. Any modification other than methionine oxidation
+        (``Oxidation``) and/or cysteine carbamidomethylation (``Carbamidomethyl``) is
+        therefore not supported and raises an ``EntrapmentError`` up front, rather than
+        silently collapsing to a mod-agnostic stem match.
 
         Parameters
         ----------
         standard_format : pd.DataFrame
             Standardised input DataFrame (output of ``convert_to_standard_format``).
-            Must contain a ``"Peptide"`` column.
+            Must contain ``"Peptide"`` and ``"Sequence"`` columns.
         max_missing_fraction : float
             Maximum tolerated fraction of unmatched peptides. Defaults to 0.03.
 
@@ -167,11 +181,27 @@ def _apply_mapping(
         Raises
         ------
         EntrapmentError
-            If the fraction of unmatched peptides exceeds ``max_missing_fraction``.
+            If the fraction of unmatched peptides exceeds ``max_missing_fraction``, if
+            any identified sequence carries a modification other than methionine
+            oxidation or cysteine carbamidomethylation, or if a sequence built only from
+            those two modifications is still absent from the mapping file.
         """
+        import re
+
         from proteobench.exceptions import EntrapmentError
 
+        allowed_modifications = {"carbamidomethyl", "oxidation"}
+        modification_pattern = re.compile(r"\[([^\[\]]*)\]")
+
+        def unsupported_modifications(sequence: str) -> List[str]:
+            return [
+                mod
+                for mod in modification_pattern.findall(sequence)
+                if mod.strip().lower() not in allowed_modifications
+            ]
+
         mapping_df = pd.read_csv(self.mapping_file, sep="\t", index_col=False)
+        mapping_df = mapping_df.drop_duplicates(subset=["sequence"], keep="first")
         all_peptides = set(standard_format["Peptide"])
         missing_peptides = all_peptides - set(mapping_df["sequence"])
         missing_fraction = len(missing_peptides) / len(all_peptides) if all_peptides else 0.0
@@ -195,12 +225,47 @@ def _apply_mapping(
         if missing_peptides:
             df = df[~df["Peptide"].isin(missing_peptides)].reset_index(drop=True)
 
+        sequences_with_unsupported_mods = {
+            sequence: mods
+            for sequence in df["Sequence"].unique()
+            for mods in (unsupported_modifications(sequence),)
+            if mods
+        }
+        if sequences_with_unsupported_mods:
+            n_bad = len(sequences_with_unsupported_mods)
+            examples = ", ".join(
+                f"{seq} ({'/'.join(mods)})" for seq, mods in list(sequences_with_unsupported_mods.items())[:5]
+            )
+            raise EntrapmentError(
+                f"{n_bad} identified sequence(s) carry a modification other than methionine oxidation "
+                f"(Oxidation) and/or cysteine carbamidomethylation (Carbamidomethyl). The entrapment mapping "
+                f"file only supports those two modifications.\n\n"
+                f"First {min(5, n_bad)} examples: {examples}"
+            )
+
+        mapping_cols = mapping_df[["sequence", "peptide_pair_index", "peptide_type"]]
+
         df = df.merge(
-            mapping_df[["sequence", "peptide_pair_index", "peptide_type"]],
+            mapping_cols.rename(columns={"sequence": "_mapped_sequence"}),
             how="left",
-            left_on="Peptide",
-            right_on="sequence",
-        ).drop(columns=["sequence"])
+            left_on="Sequence",
+            right_on="_mapped_sequence",
+        ).drop(columns=["_mapped_sequence"])
+
+        unmatched = df["peptide_pair_index"].isna()
+        if unmatched.any():
+            n_unmatched = int(unmatched.sum())
+            examples = ", ".join(sorted(df.loc[unmatched, "Sequence"].unique())[:5])
+            raise EntrapmentError(
+                f"{n_unmatched} identified sequence(s) use only methionine oxidation and/or cysteine "
+                f"carbamidomethylation, but their exact modified form is absent from the entrapment "
+                f"mapping file. This should not happen and likely indicates a gap in the mapping file "
+                f"generation; please report this as a bug.\n\n"
+                f"First {min(5, n_unmatched)} examples: {examples}"
+            )
+
+        df["peptide_pair_index"] = df["peptide_pair_index"].astype("Int64")
+
         df["Target or Entrapment"] = df["peptide_type"].replace("p_target", "entrapment")
         df = df.drop(columns=["peptide_type"])
 

webinterface/pages/markdown_files/Entrapment/DIA/ion/Astral/introduction_DIA_entrapment_ions_astral.md

@@ -4,4 +4,4 @@ The benchmark uses three technical replicates of a human plasma digest acquired
 
 ProteoBench classifies each identified precursor as target or entrapment, maps entrapment peptides to their paired target peptides, and reports lower-bound, combined, and paired false discovery proportion (FDP) estimates. These estimates are compared with the FDR threshold parsed from the submitted tool output.
 
-The module accepts DIA-NN, FragPipe, FragPipe (DIA-NN quant), and AlphaDIA outputs. For a valid comparison, keep the raw file names unchanged, use the ProteoBench entrapment FASTA, disable in-silico digestion, and do not add variable modifications.
+The module accepts DIA-NN and AlphaDIA outputs. For a valid comparison, keep the raw file names unchanged, use the ProteoBench entrapment FASTA, disable in-silico digestion, and do not add variable modifications other than methionine oxidation. Any other modification present in the identified sequences causes metric calculation to fail with an error.

webinterface/pages/markdown_files/Entrapment/DIA/ion/Astral/result_description.md

@@ -6,5 +6,5 @@ It contains the following columns:
 - **Charge**: precursor charge
 - **Q-Value**: precursor-level FDR value parsed from the submitted output
 - **Target or Entrapment**: target/entrapment classification from the ProteoBench mapping file
-- **peptide_pair_index**: identifier linking each entrapment peptide to its paired target peptide
+- **peptide_pair_index**: identifier linking each entrapment peptide to its paired target peptide. A methionine-oxidized or carbamidomethylated identification has its own index, distinct from the unmodified form of the same peptide.
 - **Lower bound FDP**, **Combined FDP**, and **Paired FDP**: empirical FDP estimates used to classify the submitted workflow