Add ProForma 2.0 peptide exporter (-epepproforma)

README.md

@@ -395,5 +395,26 @@ $ cat e_coli_compact.fasta | grep "^>" | wc -l
 Instead of directly generating a FASTA-file we first create a database, summarizing same sequences and headers. As a post-processing step, the database-entries are exported into FASTA. The generated FASTA has unique sequence as entries and offers some additionally insights. From the first entries we see peptides shared by exactly 2 proteins. Looking at the difference between the compact and
 non-compact FASTA, we see that ~244 000 entries could be summarized into already included entries in the FASTA. This generated FASTA-file can be used for identification.
 
-**NOTE:** Protein-Graphs can contain large amounts of peptides/proteins. Do a dry run with the flags `-cnp`, `-cnpm` or `cnph` (or all of them) WITHOUT the export functionality first and examine the statistics output if it is feasible to generate a FASTA-file. Without a dry run it may happen that a protein like P04637 (P53 Human) with all possible peptides and variants is exported, which will
+### Exporting peptides as ProForma 2.0
+
+Peptides can also be exported in [ProForma 2.0](https://github.com/HUPO-PSI/ProForma)
+notation (the HUPO-PSI standard for writing peptidoforms/proteoforms) into a
+tab-separated file with the columns accession, start, end, misscleavages and
+proforma:
+
+```shell
+$ protgraph -d trypsin --pep_hops 3 -fm "C:57.021464" -vm "M:15.994915" -epepproforma --pep_proforma_out peptides.proforma.tsv examples/e_coli.dat
+$ grep -m 2 "\[" peptides.proforma.tsv
+P0A884	49	53	1	RC[+57.021464]HLR
+P0A884	1	2	0	M[+15.994915]K
+```
+
+Modifications applied via `-fm`/`-vm` are known to ProtGraph only by their
+delta mass, so they are written as mass-delta tags (`C[+57.021464]`), which is
+valid ProForma. If you want named tags instead, declare a mapping for a delta
+you passed: `--pep_proforma_mod_names "15.994915=UNIMOD:35"` then writes
+`M[UNIMOD:35]`. Terminal modifications (`NPEPTERM` etc.) use the ProForma
+terminal syntax (`[+42.010565]-PEPTIDE`).
+
+**NOTE:** Protein-Graphs can contain large amounts of peptides/proteins. Do a dry run with the flags `-cnp`, `-cnpm` or `cnph` (or all of them) WITHOUT the export functionality first and examine the statistics output if it is feasible to generate a peptide output file (e.g., fasta or proforma). Without a dry run it may happen that a protein like P04637 (P53 Human) with all possible peptides and variants is exported, which will
 very likely take up all your disk space.

protgraph/cli.py

@@ -662,6 +662,58 @@ def add_fasta_peptide_export(group):
     )
 
 
+def add_proforma_peptide_export(group):
+    group.add_argument(
+        "--export_peptide_proforma", "-epepproforma", default=False, action="store_true",
+        help="Set this flag to export peptides as ProForma 2.0 strings (PSI standard) "
+        "into a single tab-separated file (accession, proforma, start, end, misscleavages, qualifiers). "
+        "Modifications from -fm/-vm are written as mass-delta tags (e.g. 'S[+79.966]')."
+    )
+    group.add_argument(
+        "--pep_proforma_out", default=os.path.join(os.getcwd(), "peptides.proforma.tsv"),
+        type=str,
+        help="Set the output file for the ProForma peptide export. "
+        "Default: '${pwd}/peptides.proforma.tsv'. NOTE: This will overwrite existing files."
+    )
+
+    def _parse_mod_name(input: str):
+        if "=" not in input:
+            raise ArgumentTypeError(
+                "A modification name mapping needs the form '<DeltaMass>=<Name>', "
+                "e.g. '79.966=UNIMOD:21'. Received: '{}'".format(input)
+            )
+        delta, name = input.split("=", 1)
+        try:
+            delta = float(delta)
+            finite = delta == delta and abs(delta) != float("inf")
+        except ValueError:
+            finite = False
+        if not finite:
+            raise ArgumentTypeError(
+                "The DeltaMass of a name mapping is not a finite number. Received: '{}'".format(input)
+            )
+        name = name.strip()
+        if not name or any(c in name for c in "[]\t\n\r"):
+            raise ArgumentTypeError(
+                "The Name of a name mapping is empty or contains brackets/whitespace "
+                "that would corrupt the output. Received: '{}'".format(input)
+            )
+        return delta, name
+    group.add_argument(
+        "--pep_proforma_mod_names", type=_parse_mod_name, action="append", default=None,
+        help="Optionally map a delta mass used in -fm/-vm to a named ProForma tag, "
+        "e.g. --pep_proforma_mod_names '79.966=UNIMOD:21' writes 'S[UNIMOD:21]' instead "
+        "of 'S[+79.966]'. Can be provided multiple times. No mass lookup is performed; "
+        "the name is taken as declared."
+    )
+    group.add_argument(
+        "--pep_proforma_write_qualifiers",
+        action="store_true",
+        help="Adds the column 'qualifiers' into the ProForma table for tracing back how the peptide "
+        "was retrieved from the protein entry (e.g., 'VARIANT[22:22, T->L, REFERENCE]))."
+    )
+
+
 def add_trie_peptide_export(group):
     group.add_argument(
         "--export_peptide_trie", "-epeptrie", default=False, action="store_true",

protgraph/export/exporters.py

@@ -79,6 +79,7 @@ def __available_exporters(self, **kwargs):
 
             # Peptide Exporter to local filesystem (no setup required)
             (kwargs["export_peptide_fasta"], "protgraph.export.peptides.pep_fasta", "PepFasta"),
+            (kwargs["export_peptide_proforma"], "protgraph.export.peptides.pep_proforma", "PepProForma"),
             (kwargs["export_peptide_trie"], "protgraph.export.peptides.pep_trie", "PepTrie"),
             (kwargs["export_peptide_sqlite"], "protgraph.export.peptides.pep_sqlite", "PepSQLite"),
         ]

protgraph/export/peptides/pep_proforma.py

@@ -0,0 +1,189 @@
+from decimal import Decimal
+from math import isfinite
+
+from Bio.SeqFeature import UnknownPosition
+
+from protgraph.export.peptides.pep_fasta import PepFasta
+from protgraph.graph_collapse_edges import Or
+
+
+class PepProForma(PepFasta):
+    """
+    Peptide exporter emitting ProForma 2.0 strings (HUPO-PSI standard notation
+    for proteoforms and peptidoforms, LeDuc et al., J. Proteome Res. 2022).
+
+    Writes one tab-separated line per exported peptide path with columns:
+        accession, proforma, start, end, misscleavages, qualifiers
+
+    ProtGraph knows modifications only by their delta mass (-fm/-vm), so
+    modifications are emitted as mass-delta tags by default ("S[+79.966]"),
+    which is valid ProForma. A user-declared mapping (--pep_proforma_mod_names
+    "79.966=UNIMOD:21") replaces a delta by a named tag ("S[UNIMOD:21]").
+    Terminal modifications (NPEPTERM etc.) use ProForma terminal syntax
+    ("[+42.010565]-PEPTIDE" or "[UNIMOD:1]-PEPTIDE").
+
+    A modification is placed on the residue its edge points at: FIXMOD/VARMOD
+    features are attached to the in-edges of the (cloned) modified node when
+    the graph is generated, so the first residue of the traversed edge's
+    target node is the modified one. This also holds for residues introduced
+    by VARIANT/MUTAGEN/CONFLICT, whose features carry no usable reference
+    location, and for nodes merged after annotation.
+
+    NOTE: as with the FASTA exporter, this exports all possible paths, so make
+    sure the traversal can terminate in forseeable future!
+    """
+
+    def start_up(self, **kwargs):
+        super(PepProForma, self).start_up(**kwargs)
+        self.output_file = kwargs["pep_proforma_out"]
+        self.mod_names = dict(kwargs["pep_proforma_mod_names"] or [])
+        self.warned_or_once = False
+        self.write_qualifiers = kwargs["pep_proforma_write_qualifiers"]
+
+    def export(self, prot_graph, queue):
+        tsv_header = [
+            "accession", "proforma", "start", "end", "misscleavages"
+        ] + (["qualifiers"] if self.write_qualifiers else [])
+
+        queue.put((self.output_file, "\t".join(tsv_header) + "\n", True, "w"))
+        super().export(prot_graph, queue)
+
+    def export_peptides(self, prot_graph, l_path_nodes, l_path_edges, l_peptide, l_miscleavages, queue):
+        entries = ""
+        for nodes, edges, misses in zip(l_path_nodes, l_path_edges, l_miscleavages):
+            # start/end come from the first/last node carrying sequence: with
+            # terminal modifications applied, nodes[1] (or nodes[-2]) can be an
+            # empty helper node (see annotate_ptms) whose position 0 does not
+            # exist in 1-based protein coordinates
+            first = next(n for n in nodes[1:-1] if prot_graph.vs[n]["aminoacid"])
+            last = next(n for n in nodes[-2:0:-1] if prot_graph.vs[n]["aminoacid"])
+            acc = self._get_accession_or_isoform(prot_graph.vs[first])
+            start_pos = self._get_position_or_isoform_position(prot_graph.vs[first])
+            end_pos = self._get_position_or_isoform_position(prot_graph.vs[last], end=True)
+            proforma = self._build_proforma(prot_graph, nodes[1:-1], edges)
+            entries += "\t".join(
+                [acc, proforma, str(start_pos), str(end_pos), str(misses)]
+                # one 'qualifiers' COLUMN, matching the header: the individual
+                # qualifier strings are comma-joined (as in the FASTA headers),
+                # and the field is present, possibly empty, on every row
+                + ([",".join(self._get_qualifiers(prot_graph, edges))]
+                   if self.write_qualifiers else [])
+            ) + "\n"
+
+        # w: ensures to have no stale entries, overwriting the file, if it exists.
+        queue.put((self.output_file, entries, False, "w"))
+
+    def _build_proforma(self, prot_graph, inner_nodes, edges):
+        """ ProForma string for one peptide path. """
+        # Get all modifications by iterating over the edges
+        n_term, c_term, by_residue = [], [], {}
+        for edge_id in edges:
+            edge = prot_graph.es[edge_id]
+            for key, delta, loc_end in self._edge_mods(edge):
+                tag = self.mod_names.get(float(delta)) or self._delta_tag(delta)
+                if key in ("NPEPTERM", "NPROTERM"):
+                    n_term.append(tag)
+                elif key in ("CPEPTERM", "CPROTERM"):
+                    c_term.append(tag)
+                else:
+                    by_residue.setdefault(edge.target, []).append((key, tag, loc_end))
+
+        # Build the final proforma sequence
+        parts = []
+        if n_term:
+            parts.append("".join("[{}]".format(t) for t in n_term) + "-")
+        for ie in inner_nodes:
+            aas = prot_graph.vs[ie]["aminoacid"]
+            node_mods = by_residue.get(ie, ())
+            if not node_mods:
+                parts.append(aas)
+            elif len(aas) == 1:
+                # Simply append UniMod
+                parts.append(aas)
+                parts.extend("[{}]".format(t) for _, t, _ in node_mods)
+            else:
+                # We need to add each PTM on its correct position in the merged
+                # node: e.g. merged node "EIN", modification I->+23.123 should
+                # yield "EI[+23.123]N" and NOT "EIN[+23.123]". Only THIS node's
+                # modifications are placed (pooling across the whole path
+                # duplicated same-letter tags), each exactly once, located via
+                # the feature's position relative to the node start; a feature
+                # without a usable position (e.g. on a variant residue) falls
+                # back to the first matching residue letter.
+                node_pos = self._get_position_or_isoform_position(prot_graph.vs[ie])
+                tags_at = {}
+                for key, tag, loc_end in node_mods:
+                    idx = None
+                    if loc_end is not None and isinstance(node_pos, int):
+                        offset = loc_end - node_pos
+                        if 0 <= offset < len(aas):
+                            idx = offset
+                    if idx is None:
+                        idx = aas.find(key[-1])  # positional keys end in the residue letter
+                        if idx == -1:
+                            idx = 0  # never drop a modification the path carries
+                    tags_at.setdefault(idx, []).append(tag)
+                parts.append("".join(
+                    aa + "".join("[{}]".format(t) for t in tags_at.get(i, ()))
+                    for i, aa in enumerate(aas)
+                ))
+
+        if c_term:
+            parts.append("-" + "".join("[{}]".format(t) for t in c_term))
+        return "".join(parts)
+
+    def _edge_mods(self, edge):
+        """ Deduplicated (key, delta, position) modifications of one traversed edge.
+
+        Or-wrapped qualifiers (collapsed parallel edges) are alternatives, not
+        a conjunction: only modifications present in EVERY branch are certain
+        for this traversal and get exported; branch-specific ones are dropped
+        with a one-time warning rather than fabricated onto every path.
+        """
+        if "qualifiers" not in edge.attributes():
+            return []
+        mods = []
+        for f in edge["qualifiers"] or []:
+            if isinstance(f, Or):
+                branches = [self._collect_mods(branch) for branch in f]
+                common = set.intersection(*map(set, branches)) if branches else set()
+                if any(set(b) - common for b in branches) and not self.warned_or_once:
+                    self.warned_or_once = True
+                    print(
+                        "Warning: modifications differ between collapsed edge alternatives; "
+                        "only modifications shared by all alternatives are exported. "
+                        "Re-run with --no_collapsing_edges for exact per-path modifications."
+                    )
+                mods.extend(m for m in dict.fromkeys(x for b in branches for x in b) if m in common)
+            elif getattr(f, "type", None) in ("FIXMOD", "VARMOD"):
+                mods.append(self._feature_mod(f))
+        return list(dict.fromkeys(mods))
+
+    def _collect_mods(self, qualifier):
+        out = []
+        for f in qualifier or []:
+            if isinstance(f, Or):
+                for branch in f:
+                    out.extend(self._collect_mods(branch))
+            elif getattr(f, "type", None) in ("FIXMOD", "VARMOD"):
+                out.append(self._feature_mod(f))
+        return out
+
+    def _feature_mod(self, feature):
+        key, _, delta = feature.qualifiers["note"].rpartition(":")
+        end = feature.location.end
+        return key, delta, None if isinstance(end, UnknownPosition) else int(end)
+
+    def _delta_tag(self, delta):
+        """ '79.966' -> '+79.966': explicitly signed, plain decimal notation.
+
+        ProForma's mass grammar has no exponent and no non-finite forms, so
+        '1e-05' is rewritten as '+0.00001' and nan/inf are refused.
+        """
+        if not isfinite(float(delta)):
+            raise ValueError(
+                "Cannot express the delta mass '{}' as a ProForma mass tag".format(delta)
+            )
+        if "e" in delta.lower():
+            delta = format(Decimal(delta), "f")
+        return delta if delta.startswith("-") else "+" + delta

protgraph/protgraph.py

@@ -184,6 +184,7 @@ def __call__(self, parser, namespace, values, option_string=None):
         ("citus_peptide_export", cli.add_citus_peptide_export),
         ("sqlite_peptide_export", cli.add_sqlite_peptide_export),
         ("fasta_peptide_export", cli.add_fasta_peptide_export),
+        ("proforma_peptide_export", cli.add_proforma_peptide_export),
         ("trie_peptide_export", cli.add_trie_peptide_export),
         ("gremlin_graph_export", cli.add_gremlin_graph_export),
     ]

setup.py

@@ -11,7 +11,7 @@
 
 setup(
     name='protgraph',
-    version='0.3.12',
+    version='0.3.13',
     author="Dominik Lux",
     description="ProtGraph, a graph generator for proteins.",
     long_description=long_description,

tests/fixtures/merged.txt

@@ -0,0 +1,12 @@
+ID   MERGED_TEST             Reviewed;          11 AA.
+AC   X9MRG0;
+DT   01-JAN-2020, integrated into UniProtKB/Swiss-Prot.
+DT   01-JAN-2020, sequence version 1.
+DT   01-JAN-2020, entry version 1.
+DE   RecName: Full=Merged node modification test;
+OS   Homo sapiens (Human).
+OC   Eukaryota; Metazoa; Chordata.
+OX   NCBI_TaxID=9606;
+SQ   SEQUENCE   11 AA;  1200 MW;  0123456789ABCDEF CRC64;
+     KEINKACKACK
+//

tests/fixtures/minimal.txt

@@ -0,0 +1,16 @@
+ID   TEST_HUMAN              Reviewed;           8 AA.
+AC   X9TEST;
+DT   01-JAN-2020, integrated into UniProtKB/Swiss-Prot.
+DT   01-JAN-2020, sequence version 1.
+DT   01-JAN-2020, entry version 1.
+DE   RecName: Full=Minimal test protein;
+GN   Name=TEST;
+OS   Homo sapiens (Human).
+OC   Eukaryota; Metazoa; Chordata.
+OX   NCBI_TaxID=9606;
+FT   VARIANT         4
+FT                   /note="T -> M (test variant)"
+FT                   /id="VAR_900001"
+SQ   SEQUENCE   8 AA;  920 MW;  0123456789ABCDEF CRC64;
+     MKCTMSAK
+//