Add msfragger tsv reader

alphabase/constants/const_files/psm_reader.yaml

@@ -153,6 +153,39 @@ pfind:
     'decoy': ['Target/Decoy', 'Targe/Decoy']
   modification_mapping_type: null  # no mapping required
 
+# https://fragpipe.nesvilab.org/docs/tutorial_fragpipe_outputs.html#psmtsv
+msfragger_psm_tsv:
+  reader_type: msfragger_psm_tsv
+  rt_unit: second
+  column_mapping:
+    'sequence': 'Peptide'
+    'charge': 'Charge'
+    'rt': 'Retention'
+    'score': 'Hyperscore'
+    'proteins': 'Protein'
+    'decoy': 'Is Decoy'
+    'raw_name': 'raw_name'
+    'scan_num': 'scan_num'
+    'precursor_mz': 'Observed M/Z'
+    '_tmp_mods': 'Assigned Modifications'
+    'intensity': 'Intensity'
+  mass_mapped_mods:
+    - 'Oxidation@M'
+    - 'Carbamidomethyl@C'
+    - 'Phospho@S'
+    - 'Phospho@T'
+    - 'Phospho@Y'
+    - 'GlyGly@K'
+    - 'Acetyl@Any_N-term'
+    - 'Glu->pyro-Glu@E^Any_N-term'
+    - 'Gln->pyro-Glu@Q^Any_N-term'
+    - 'TMT6plex@K'
+    - 'TMT6plex@Any_N-term'
+    - 'TMTpro@K'
+    - 'TMTpro@Any_N-term'
+  mod_mass_tol: 0.1
+  modification_mapping_type: 'maxquant'
+
 msfragger_pepxml:
   reader_type: msfragger_pepxml
   rt_unit: second

alphabase/psm_reader/__init__.py

@@ -8,7 +8,7 @@
     "SwathReader",
     "SpectronautReportReader",
     "MaxQuantReader",
-    "MSFragger_PSM_TSV_Reader",
+    "MSFraggerPsmTsvReader",
     "pFindReader",
     "MSFraggerPepXML",
     "MSFraggerPepXMLReader",
@@ -47,9 +47,9 @@
     register_readers as register_mq_readers,
 )
 from alphabase.psm_reader.msfragger_reader import (
-    MSFragger_PSM_TSV_Reader,
     MSFraggerPepXML,
     MSFraggerPepXMLReader,
+    MSFraggerPsmTsvReader,
 )
 from alphabase.psm_reader.msfragger_reader import (
     register_readers as register_fragger_readers,

alphabase/psm_reader/keys.py

@@ -30,6 +30,10 @@ class PsmDfCols(metaclass=ConstantsClass):
         Peptide amino-acid sequence without modification annotations (plain sequence of residues).
     DECOY
         Label indicating whether the matched sequence is a decoy (reverse/shuffled) entry used for FDR estimation.
+    TMP_MODS
+        Intermediate column containing raw modification strings from search engines (e.g., MSFragger's
+        "Assigned Modifications" format: "5S(79.9663), N-term(304.2071)"). This is translated into
+        the standardized MODS and MOD_SITES columns during processing and is not present in the final output.
     MODS
         List of modifications on the peptide in unimod standardization
     SCORE
@@ -105,6 +109,7 @@ class PsmDfCols(metaclass=ConstantsClass):
     MODIFIED_SEQUENCE = "modified_sequence"
     SEQUENCE = "sequence"
     DECOY = "decoy"
+    TMP_MODS = "_tmp_mods"
     MODS = "mods"
     SCORE = "score"
     TO_REMOVE = "to_remove"

alphabase/psm_reader/msfragger_reader.py

@@ -18,10 +18,186 @@
 )
 
 
+class MsFraggerTokens:
+    """String tokens used in MSFragger output formats."""
+
+    MOD_START = "("
+    MOD_STOP = ")"
+    MOD_SEPARATOR = ","
+    N_TERM = "N-term"
+    C_TERM = "C-term"
+
+
 def _is_fragger_decoy(proteins: List[str]) -> bool:
     return all(prot.lower().startswith("rev_") for prot in proteins)
 
 
+class MSFraggerModificationTranslation:
+    """Translate MSFragger PSM.TSV modifications to alphabase format."""
+
+    def __init__(
+        self,
+        mass_mapped_mods: List[str],
+        mod_mass_tol: float = 0.1,
+    ):
+        """Initialize MSFragger modification translator.
+
+        Parameters
+        ----------
+        mass_mapped_mods : List[str]
+            List of modification names to match against (e.g., ['Phospho@S', 'Oxidation@M'])
+        mod_mass_tol : float
+            Mass tolerance for matching modifications in Daltons. Default: 0.1
+
+        """
+        self._mass_mapped_mods = mass_mapped_mods
+        self._mod_mass_tol = mod_mass_tol
+
+    def __call__(self, psm_df: pd.DataFrame) -> pd.DataFrame:
+        """Translate modifications from MSFragger assigned modifications.
+
+        Parameters
+        ----------
+        psm_df : pd.DataFrame
+            DataFrame with PsmDfCols.TMP_MODS column containing raw assigned modifications strings
+
+        Returns
+        -------
+        pd.DataFrame
+            The input DataFrame with 'mods' and 'mod_sites' columns added
+
+        """
+        mods_list = []
+        sites_list = []
+
+        for _, row in psm_df.iterrows():
+            assigned_mods = row.get(PsmDfCols.TMP_MODS, "")
+            mods, sites = self._parse_assigned_modifications(assigned_mods)
+            mods_list.append(mods)
+            sites_list.append(sites)
+
+        psm_df[PsmDfCols.MODS] = mods_list
+        psm_df[PsmDfCols.MOD_SITES] = sites_list
+
+        return psm_df
+
+    def _parse_assigned_modifications(self, assigned_mods: str) -> Tuple[str, str]:
+        """Parse MSFragger Assigned Modifications string.
+
+        Directly maps mass shifts to modification names without conversion.
+
+        Parameters
+        ----------
+        assigned_mods : str
+            MSFragger format: "5S(79.9663), 7S(79.9663), N-term(304.2071)"
+
+        Returns
+        -------
+        tuple
+            (mods_str, sites_str) where mods and sites are semicolon-separated.
+            Example: ("Phospho@S;Phospho@S;TMT@Any_N-term", "5;7;0")
+
+        """
+        if not assigned_mods:
+            return "", ""
+
+        mod_entries = [
+            m.strip() for m in assigned_mods.split(MsFraggerTokens.MOD_SEPARATOR)
+        ]
+        mods = []
+        sites = []
+
+        for entry in mod_entries:
+            if not entry:
+                continue
+
+            if entry.startswith(MsFraggerTokens.N_TERM):
+                mass_shift = float(
+                    entry.split(MsFraggerTokens.MOD_START)[1].rstrip(
+                        MsFraggerTokens.MOD_STOP
+                    )
+                )
+                mod_name = self._match_mod_by_mass(mass_shift, "Any_N-term")
+                mods.append(mod_name)
+                sites.append("0")
+            elif entry.startswith(MsFraggerTokens.C_TERM):
+                mass_shift = float(
+                    entry.split(MsFraggerTokens.MOD_START)[1].rstrip(
+                        MsFraggerTokens.MOD_STOP
+                    )
+                )
+                mod_name = self._match_mod_by_mass(mass_shift, "Any_C-term")
+                mods.append(mod_name)
+                sites.append("-1")
+            else:
+                parts = entry.split(MsFraggerTokens.MOD_START)
+                if len(parts) != 2:  # noqa: PLR2004
+                    continue
+
+                position_and_amino_acid, mass_string = parts
+                mass_shift = float(mass_string.rstrip(MsFraggerTokens.MOD_STOP))
+
+                position = ""
+                amino_acid = ""
+                for char in position_and_amino_acid:
+                    if char.isdigit():
+                        position += char
+                    else:
+                        amino_acid += char
+
+                if position and amino_acid:
+                    mod_name = self._match_mod_by_mass(mass_shift, amino_acid)
+                    mods.append(mod_name)
+                    sites.append(position)
+
+        return ";".join(mods), ";".join(sites)
+
+    def _match_mod_by_mass(self, mass_shift: float, aa_or_term: str) -> str:
+        """Match mass shift to modification name.
+
+        Parameters
+        ----------
+        mass_shift : float
+            Mass shift in Daltons (e.g., 79.9663 for phosphorylation)
+        aa_or_term : str
+            Amino acid single letter code or terminal (Any_N-term, Any_C-term)
+
+        Returns
+        -------
+        str
+            Modification name in alphabase format (e.g., 'Phospho@S')
+
+        Raises
+        ------
+        ValueError
+            If no matching modification is found
+
+        """
+        for mod_name in self._mass_mapped_mods:
+            if mod_name not in MOD_MASS:
+                continue
+
+            mod_mass = MOD_MASS[mod_name]
+
+            if abs(mass_shift - mod_mass) < self._mod_mass_tol:
+                mod_base, mod_site = mod_name.split("@")
+
+                if aa_or_term in ["Any_N-term", "Any_C-term"]:
+                    if mod_site == aa_or_term or mod_site.endswith(aa_or_term):
+                        return mod_name
+                else:
+                    if mod_site == aa_or_term:
+                        return mod_name
+                    if mod_site in ["X", "Any"]:
+                        return f"{mod_base}@{aa_or_term}"
+
+        raise ValueError(
+            f"Unknown modification: mass_shift={mass_shift:.4f} at {aa_or_term}. "
+            f"Add the modification to 'mass_mapped_mods' in psm_reader.yaml or extend "
+            f"the reader's _mass_mapped_mods list before importing."
+        )
+
+
 def _get_mods_from_masses(  # noqa: PLR0912, C901 too many branches, too complex TODO: refactor
     sequence: str,
     msf_aa_mods: List[str],
@@ -82,15 +258,71 @@ def _get_mods_from_masses(  # noqa: PLR0912, C901 too many branches, too complex
     )
 
 
-class MSFragger_PSM_TSV_Reader(PSMReaderBase):  # noqa: N801 name should use CapWords convention TODO: refactor
+class MSFraggerPsmTsvReader(PSMReaderBase):
     """Reader for MSFragger's psm.tsv file."""
 
+    _reader_type = "msfragger_psm_tsv"
+
     def __init__(
         self,
+        *,
+        column_mapping: Optional[dict] = None,
+        modification_mapping: Optional[dict] = None,
+        fdr: float = 0.01,
+        keep_decoy: bool = False,
+        rt_unit: Optional[str] = None,
         **kwargs,
     ):
-        """Constructor."""
-        raise NotImplementedError("MSFragger_PSM_TSV_Reader for psm.tsv")
+        """Initialize MSFragger PSM TSV reader.
+
+        See PSMReaderBase documentation for parameters.
+
+        """
+        super().__init__(
+            column_mapping=column_mapping,
+            modification_mapping=modification_mapping,
+            fdr=fdr,
+            keep_decoy=keep_decoy,
+            rt_unit=rt_unit,
+            **kwargs,
+        )
+        self._mass_mapped_mods = psm_reader_yaml.get(self._reader_type, {}).get(
+            "mass_mapped_mods", []
+        )
+        self._mod_mass_tol = psm_reader_yaml.get(self._reader_type, {}).get(
+            "mod_mass_tol", 0.1
+        )
+
+    def _translate_modifications(self) -> None:
+        """No-op: modification translation is handled in _load_modifications."""
+
+    def _load_file(self, filename: str) -> pd.DataFrame:
+        """Load MSFragger PSM TSV file."""
+        return pd.read_csv(filename, sep="\t", keep_default_na=False)
+
+    def _pre_process(self, df: pd.DataFrame) -> pd.DataFrame:
+        """MSFragger PSM TSV preprocessing."""
+        df.fillna("", inplace=True)
+        df[[PsmDfCols.RAW_NAME, PsmDfCols.SCAN_NUM]] = (
+            df["Spectrum"].str.split(".").apply(lambda x: pd.Series([x[0], int(x[1])]))
+        )
+        return df
+
+    def _translate_decoy(self) -> None:
+        self._psm_df[PsmDfCols.DECOY] = (
+            self._psm_df[PsmDfCols.DECOY] == "true"
+        ).astype(np.int8)
+
+    def _translate_score(self) -> None:
+        """MSFragger Hyperscore is already in correct format (larger = better)."""
+
+    def _load_modifications(self, origin_df: pd.DataFrame) -> None:  # noqa: ARG002
+        """Parse modifications from PsmDfCols.TMP_MODS column (mapped from 'Assigned Modifications')."""
+        translator = MSFraggerModificationTranslation(
+            mass_mapped_mods=self._mass_mapped_mods,
+            mod_mass_tol=self._mod_mass_tol,
+        )
+        self._psm_df = translator(self._psm_df)
 
 
 class MSFraggerPepXMLReader(PSMReaderBase):
@@ -140,7 +372,7 @@ def __init__(  # noqa: PLR0913, D417 # too many arguments in function definition
         self._mod_mass_tol = psm_reader_yaml["msfragger_pepxml"]["mod_mass_tol"]
 
     def _translate_modifications(self) -> None:
-        pass
+        """No-op: modification translation is handled in _load_modifications."""
 
     def _load_file(self, filename: str) -> pd.DataFrame:
         """Load a MsFragger output file to a DataFrame."""
@@ -221,6 +453,6 @@ def __init__(self, *args, **kwargs):
 
 def register_readers() -> None:
     """Register MSFragger readers."""
-    psm_reader_provider.register_reader("msfragger_psm_tsv", MSFragger_PSM_TSV_Reader)
-    psm_reader_provider.register_reader("msfragger", MSFragger_PSM_TSV_Reader)
+    psm_reader_provider.register_reader("msfragger_psm_tsv", MSFraggerPsmTsvReader)
+    psm_reader_provider.register_reader("msfragger", MSFraggerPsmTsvReader)
     psm_reader_provider.register_reader("msfragger_pepxml", MSFraggerPepXMLReader)

tests/integration/test_psm_readers.py

@@ -20,6 +20,7 @@
     AlphaDiaReaderTsv,
     DiannReader,
     MaxQuantReader,
+    MSFraggerPsmTsvReader,
     SageReaderTSV,
     SpectronautReader,
     SpectronautReportReader,
@@ -193,6 +194,13 @@
 _EVLHLLR_	2	21.775629	0.837	EVLHLLR	440.274169715862	A0AVF1	TTC26	b	705.4294	0.02848413	1	6	noloss
 """
 
+TEST_DATA_MSFRAGGER_PSM_TSV = """Spectrum	Spectrum File	Peptide	Modified Peptide	Extended Peptide	Prev AA	Next AA	Peptide Length	Charge	Retention	Observed Mass	Calibrated Observed Mass	Observed M/Z	Calibrated Observed M/Z	Calculated Peptide Mass	Calculated M/Z	Delta Mass	Expectation	Hyperscore	Nextscore	Probability	Qvalue	Number of Enzymatic Termini	Number of Missed Cleavages	Protein Start	Protein End	Intensity	Assigned Modifications	Observed Modifications	Purity	Is Decoy	Protein	Protein ID	Entry Name	Gene
+file1.01234.01234.3	interact-file1.pep.xml	LENPNNK	LENP[+80]NNK	K.LENPNNK.S	K	S	7	3	228.78	1438.80	1438.79	479.9420	479.9396	1438.79	479.9383	0.01	0.04	16.18	15.64	0.99	0.0001	2	0	93	104	1000.0	7K(304.2071), N-term(304.2071)	0.0	false	tr|Q06AA7|Q06AA7_PIG	Q06AA7	Q06AA7_PIG	TRA2B
+file1.01235.01235.2	interact-file1.pep.xml	PEPTIDE	PEPTIDE	R.PEPTIDE.K	R	K	7	2	120.50	750.35	750.34	375.6750	375.6700	750.34	375.6725	0.01	0.08	12.50	10.20	0.85	0.002	2	0	45	51	800.0		0.0	false	tr|P12345|P12345_HUMAN	P12345	P12345_HUMAN	GENE1
+file2.01236.01236.3	interact-file2.pep.xml	SEQUENCE	SEQUENCE	K.SEQUENCE.R	K	R	8	3	150.75	900.42	900.41	300.8067	300.8033	900.41	300.8050	0.01	0.02	18.90	16.50	0.99	0.0000	2	0	120	127	1500.0		0.0	false	sp|P98765|P98765_MOUSE	P98765	P98765_MOUSE	GENE2
+file2.01237.01237.2	interact-file2.pep.xml	TESTM	TESTM[+16]	R.TESTM.A	R	A	5	2	95.30	600.28	600.27	300.6383	300.6350	600.27	300.6375	0.01	0.10	10.20	8.50	0.75	0.005	2	0	88	92	500.0	5M(15.9949)	0.0	false	tr|Q11111|Q11111_RAT	Q11111	Q11111_RAT	GENE3
+"""
+
 TEST_DATA_ALPHADIA = """base_width_mobility	base_width_rt	rt_observed	mobility_observed	mono_ms1_intensity	top_ms1_intensity	sum_ms1_intensity	weighted_ms1_intensity	weighted_mass_deviation	weighted_mass_error	mz_observed	mono_ms1_height	top_ms1_height	sum_ms1_height	weighted_ms1_height	isotope_intensity_correlation	isotope_height_correlation	n_observations	intensity_correlation	height_correlation	intensity_fraction	height_fraction	intensity_fraction_weighted	height_fraction_weighted	mean_observation_score	sum_b_ion_intensity	sum_y_ion_intensity	diff_b_y_ion_intensity	f_masked	fragment_scan_correlation	template_scan_correlation	fragment_frame_correlation	top3_frame_correlation	template_frame_correlation	top3_b_ion_correlation	n_b_ions	top3_y_ion_correlation	n_y_ions	cycle_fwhm	mobility_fwhm	delta_frame_peak	top_3_ms2_mass_error	mean_ms2_mass_error	n_overlapping	mean_overlapping_intensity	mean_overlapping_mass_error	precursor_idx	rank	scan_center	frame_start	elution_group_idx	frame_stop	frame_center	scan_stop	score	scan_start	i_3	genes	i_1	decoy	mods	channel	rt_library	mobility_library	i_0	flat_frag_start_idx	mod_sites	i_2	flat_frag_stop_idx	mz_library	rt_calibrated	charge	mz_calibrated	sequence	proteins	delta_rt	n_K	n_R	n_P	_decoy	proba	qval	_candidate_idx	valid	candidate_idx	run	mod_seq_hash	mod_seq_charge_hash	pg_master	pg	pg_qval	intensity
 0.000000	38.684296	468.905731	0.000001	117961664.000000	117961664.000000	276840416.000000	95120688.000000	-0.197439	0.197439	781.368469	2849874.750000	2849874.750000	6788314.000000	2273957.500000	0.980586	0.999165	2.000000	0.992594	0.992708	1.000000	1.000000	1.000000	1.000000	0.999540	12.188647	15.848632	-3.659985	1.000000	0.000000	0.000000	0.798432	0.932817	0.908939	0.903936	1.000000	0.898678	11.000000	16.475826	0.000000	-0.950898	0.200608	0.123120	0.000000	0.000000	15.000000	138252	0	0	86989	144294	94514	90601	1	165.561066	0	0.067197	ACTN4_HUMAN;ACTN1_HUMAN	0.340936	0	Carbamidomethyl@C	0	649.565613	0.000000	0.416578	798982	10	0.175290	798994	781.369385	503.363037	2	781.368652	ELPPDQAEYCIAR	O43707;P12814	-34.457306	0	1	2	0.000000	0.000000	0.000000	138252	True	138252	20231017_OA2_TiHe_ADIAMA_HeLa_200ng_Evo011_21min_F-40_05	356312820745498424	356312820745498426	O43707	O43707;P12814	0.000000	262684832.000000
 0.000000	29.428162	799.346985	0.000001	172471648.000000	172471648.000000	372443296.000000	132807928.000000	-0.149760	0.149760	780.903625	6140298.500000	6140298.500000	13225659.000000	4718308.000000	0.999529	0.999349	2.000000	0.839249	0.848623	1.000000	1.000000	1.000000	1.000000	0.993265	12.386384	15.716099	-3.329715	1.000000	0.000000	0.000000	0.915426	0.954805	0.895610	0.949662	1.000000	0.924100	11.000000	16.276836	0.000000	4.326714	0.438759	0.270384	0.000000	0.000000	15.000000	70824	0	0	152005	77362	157724	154714	1	179.995453	0	0.054696	MDHM_HUMAN	0.353536	0		0	934.109009	0.000000	0.428233	395635		0.163536	395647	780.904480	801.701782	2	780.903748	VDFPQDQLTALTGR	P40926	-2.354797	0	1	1	0.000000	0.000000	0.000000	70824	True	70824	20231017_OA2_TiHe_ADIAMA_HeLa_200ng_Evo011_21min_F-40_05	11614389154130684956	11614389154130684958	P40926	P40926	0.000000	475047904.000000
@@ -493,3 +501,13 @@ def test_msfragger_speclib_reader() -> None:
     reader.import_file(input_data)
 
     _assert_reference_df_equal(reader.psm_df, "msfragger_speclib")
+
+
+def test_msfragger_psm_tsv_reader() -> None:
+    """Test the MSFragger PSM TSV reader."""
+    input_data = StringIO(TEST_DATA_MSFRAGGER_PSM_TSV)
+
+    reader = MSFraggerPsmTsvReader()
+    reader.import_file(input_data)
+
+    _assert_reference_df_equal(reader.psm_df, "msfragger_psm_tsv")