diff --git a/README.md b/README.md index 4eb2b8f..c0ac51f 100644 --- a/README.md +++ b/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. diff --git a/protgraph/cli.py b/protgraph/cli.py index 63f4d18..bbf2bbb 100644 --- a/protgraph/cli.py +++ b/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 '=', " + "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", diff --git a/protgraph/export/exporters.py b/protgraph/export/exporters.py index 131410d..da36ed6 100644 --- a/protgraph/export/exporters.py +++ b/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"), ] diff --git a/protgraph/export/peptides/pep_proforma.py b/protgraph/export/peptides/pep_proforma.py new file mode 100644 index 0000000..1196ea8 --- /dev/null +++ b/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 diff --git a/protgraph/protgraph.py b/protgraph/protgraph.py index 6547f15..a1db883 100644 --- a/protgraph/protgraph.py +++ b/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), ] diff --git a/setup.py b/setup.py index 79ae5ed..49ab675 100644 --- a/setup.py +++ b/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, diff --git a/tests/fixtures/merged.txt b/tests/fixtures/merged.txt new file mode 100644 index 0000000..0da3446 --- /dev/null +++ b/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 +// diff --git a/tests/fixtures/minimal.txt b/tests/fixtures/minimal.txt new file mode 100644 index 0000000..d224351 --- /dev/null +++ b/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 +// diff --git a/tests/functional_test.py b/tests/functional_test.py index d27829c..d80c3fe 100644 --- a/tests/functional_test.py +++ b/tests/functional_test.py @@ -214,3 +214,101 @@ def test_issue13(self): def test_issue41(self): args = protgraph.parse_args(["-n", "1", "-epepfasta", os.path.join(self.examples_path, "P49782.txt")]) protgraph.prot_graph(**args) + + def _run_proforma(self, extra_args, out): + """ ProForma export of tests/fixtures/minimal.txt (MKCTMSAK, VARIANT T4M). """ + fixture = os.path.abspath(os.path.join(os.path.dirname(__file__), "fixtures", "minimal.txt")) + args = protgraph.parse_args( + ["-epepproforma", "--pep_proforma_out", out, "--pep_hops", "10", "-d", "trypsin"] + + extra_args + self.procs_num + [fixture] + ) + protgraph.prot_graph(**args) + with open(out) as f: + return [ln.rstrip("\n").split("\t") for ln in f if ln.strip()] + + def test_export_pep_proforma_smoke(self): + args = protgraph.parse_args(["-epepproforma", "--pep_hops", "2"] + self.procs_num + self.example_files) + protgraph.prot_graph(**args) + + def test_export_pep_proforma_pins_mod_placement(self): + import tempfile + with tempfile.TemporaryDirectory() as tmp: + lines = self._run_proforma( + ["-ft", "VARIANT", "-fm", "C:57.021464", "-vm", "M:15.994915", + "--pep_proforma_mod_names", "15.994915=UNIMOD:35"], + os.path.join(tmp, "p.tsv"), + ) + assert all(len(x) == 5 for x in lines) + rows = {tuple(x) for x in lines} + proformas = {x[1] for x in lines} + # exact placement with coordinates: the fixed mod sits on the C + assert ("X9TEST", "C[+57.021464]TMSAK", "3", "8", "0",) in rows + # a variable mod on a VARIANT-introduced residue (no reference + # position!) is placed on that residue, not dropped + assert "C[+57.021464]M[UNIMOD:35]MSAK" in proformas + # variable mods also export their unmodified counterparts + assert {"MK", "M[UNIMOD:35]K"} <= proformas + # the fixed mod is applied to every C, always signed + assert all("C[+57.021464]" in p for p in proformas if "C" in p) + + def test_export_pep_proforma_terminal_mods_and_overwrite(self): + import tempfile + with tempfile.TemporaryDirectory() as tmp: + out = os.path.join(tmp, "p.tsv") + lines = self._run_proforma(["-fm", "NPEPTERM:42.010565", "-vm", "CPEPTERM:79.966"], out) + proformas = {x[1] for x in lines} + proformas.discard("proforma") + assert all(p.startswith("[+42.010565]-") for p in proformas) + # variable C-terminal: modified and unmodified paths both exported + assert any(p.endswith("-[+79.966]") for p in proformas) + assert any(not p.endswith("-[+79.966]") for p in proformas) + # a rerun overwrites instead of appending stale results + again = self._run_proforma(["-fm", "NPEPTERM:42.010565", "-vm", "CPEPTERM:79.966"], out) + assert len(again) == len(lines) + + def test_export_pep_proforma_places_mods_inside_merged_nodes_exactly_once(self): + import re + import tempfile + fixture = os.path.abspath(os.path.join(os.path.dirname(__file__), "fixtures", "merged.txt")) + with tempfile.TemporaryDirectory() as tmp: + out = os.path.join(tmp, "p.tsv") + args = protgraph.parse_args( + [ + "-epepproforma", "--pep_proforma_out", out, "--pep_hops", "10", + "-d", "trypsin", "--pep_miscleavages", "2", + "-fm", "I:23.123", "-fm", "C:57.021464", + ] + self.procs_num + [fixture] + ) + protgraph.prot_graph(**args) + with open(out) as f: + next(f) # header + proformas = {ln.rstrip("\n").split("\t")[1] for ln in f if ln.strip()} + # a mod on an interior residue of a merged node sits on that residue + assert "EI[+23.123]NK" in proformas + # a miscleaved peptide spanning two nodes that EACH carry the same + # fixed mod gets one tag per residue, not the pooled duplicates + assert "AC[+57.021464]KAC[+57.021464]K" in proformas + for p in proformas: + assert "[+57.021464][+57.021464]" not in p and "[+23.123][+23.123]" not in p + assert not re.search(r"\[[^]]*\[", p), p # no tag nested inside a tag + + def test_export_pep_proforma_qualifiers_is_one_column(self): + import tempfile + fixture = os.path.abspath(os.path.join(os.path.dirname(__file__), "fixtures", "merged.txt")) + with tempfile.TemporaryDirectory() as tmp: + out = os.path.join(tmp, "p.tsv") + args = protgraph.parse_args( + [ + "-epepproforma", "--pep_proforma_out", out, "--pep_hops", "10", + "-d", "trypsin", "--pep_miscleavages", "2", + "-fm", "C:57.021464", "--pep_proforma_write_qualifiers", + ] + self.procs_num + [fixture] + ) + protgraph.prot_graph(**args) + with open(out) as f: + widths = {len(ln.rstrip("\n").split("\t")) for ln in f if ln.strip()} + f.seek(0) + next(f) + quals = [ln.rstrip("\n").split("\t")[5] for ln in f if ln.strip()] + assert widths == {6} # header and every row: exactly one qualifiers column + assert any("FIXMOD" in q for q in quals)