GenoToPep

GenoToPep is a reproducible command-line package for generating peptide candidate universes from metagenome assemblies or predicted-protein FASTA files.

Installation

Python 3.10 or newer is required.

PyPI

python -m pip install genotopep
genotopep --version
genotopep --help

Bioconda

GenoToPep is available from Bioconda. The current public Bioconda release is 0.1.0.1 (noarch, build pyhdfd78af_0):

conda install -c conda-forge -c bioconda genotopep

For an exact reproducible install of the current public release:

conda install -c conda-forge -c bioconda genotopep=0.1.0.1

The local 0.2.0 development tree documented below is not yet published to Bioconda.

Input routes

• Nucleotide assembly: --assembly sample.fna[.gz]
• Predicted proteins: --proteins sample.faa[.gz]
• Predicted proteins with genomic provenance: --proteins sample.faa --gff sample.gff3[.gz]

--assembly and --proteins are mutually exclusive. GFF3 is optional and is only valid with FAA input.

FNA ORF calling

Assembly input is translated with Pyrodigal metagenome mode [1,2]. The default short-ORF threshold is:

min_protein_aa = 10
min_gene_nt = 10 * 3 + 3 = 33 nt
closed = false
max_overlap = 0

Pyrodigal still applies its coding and start-site models: 33 nt is an eligible minimum, not a guarantee that every open reading frame of that length is emitted. Contig-edge partial calls are retained and marked in the ORF provenance output. Terminal stop symbols are not written to protein FASTA.

Generation strategies

Intact

Emits the complete supplied or Pyrodigal-called protein when its length is within the inclusive range.

genotopep generate intact \
  --assembly assembly.fna.gz \
  --min-length 10 --max-length 100 \
  --output results/assembly-intact

genotopep generate intact \
  --proteins proteins.faa.gz \
  --min-length 10 --max-length 100 \
  --output results/proteins-intact

Windows

Emits contiguous subsequences for every inclusive length, using stride 1 by default.

genotopep generate windows \
  --assembly assembly.fna.gz \
  --min-length 10 --max-length 50 --stride 1 \
  --output results/assembly-windows

genotopep generate windows \
  --proteins proteins.faa.gz \
  --min-length 10 --max-length 50 --stride 1 \
  --output results/proteins-windows

Cleavage

Builds a boundary union from selected deterministic rules, creates elementary fragments, and emits every length-valid joining of adjacent fragments.

genotopep generate cleavage \
  --assembly assembly.fna.gz \
  --cleavage-mode all \
  --min-length 10 --max-length 100 \
  --output results/assembly-cleavage

genotopep generate cleavage \
  --proteins proteins.faa.gz \
  --cleavage-mode trypsin cnbr \
  --min-length 10 --max-length 100 \
  --output results/proteins-cleavage

Supported lowercase modes follow simplified deterministic profiles derived from published specificity summaries and ExPASy PeptideCutter conventions [3-5]:

• trypsin: simplified deterministic profile that cuts after K/R unless followed by P; reported context-dependent exceptions and counter-exceptions are not modeled [3,5]
• chymotrypsin: high-specificity profile that cuts after F/Y unless followed by P and after W unless followed by M/P [3,5]
• cnbr: boundary after M; product chemistry, methionine oxidation, and condition-dependent incomplete cleavage are not modeled [5,6]
• all: union of trypsin, chymotrypsin, and CNBr boundaries

Omitting --cleavage-mode resolves to all. Multiple individual modes are accepted. all cannot be combined with an individual mode. The cited sources support residue-specific boundary conventions; the union of selected boundaries and enumeration of every contiguous adjacent-fragment joining are GenoToPep-defined candidate-generation semantics.

Pre-generation estimates

estimate reports the exact occurrence count needed for workload planning without writing peptide outputs or performing sequence deduplication. An occurrence is a location-derived candidate, including duplicate sequences: if the same sequence is generated from multiple proteins or coordinates, each location is counted separately. Therefore unique peptides <= retained occurrences; the retained occurrence count is an exact upper bound on the final exact-sequence unique peptide count. Runtime, occurrence-TSV rows, and candidate-processing work scale primarily with occurrences rather than unique sequences.

The exclusion count is occurrence-level too. Because the exclusion TSV aggregates by source protein, strategy, and reason, one exclusion row can represent many excluded occurrences. counts.exclusions is the excluded-candidate count, whereas counts.exclusion_groups is the number of written aggregate rows.

genotopep estimate windows \
  --proteins proteins.faa.gz \
  --min-length 10 --max-length 50 --stride 1

genotopep estimate intact cleavage \
  --assembly assembly.fna.gz \
  --cleavage-mode trypsin cnbr \
  --min-length 10 --max-length 100

genotopep estimate \
  --assembly assembly.fna.gz \
  --cleavage-mode all \
  --min-length 10 --max-length 100 --stride 1

Omitting strategy selection defaults to all; explicitly specifying strategy all has the same effect. Both expand to intact, windows, and cleavage. all cannot be combined with individual strategy names. For assembly input, Pyrodigal ORF calling is performed to obtain the protein sequences required for exact counting, but ORF and peptide output files are not written. Reported occurrences are exact pre-deduplication counts and an upper bound on unique peptides. The exact unique count is produced by generation and exact-sequence deduplication.

Estimate JSON schema version 1 includes strategy-specific counts, arithmetic totals across the selected strategies, and non-blocking scale warnings. A strategy with at least 10,000,000 candidate events (occurrences + exclusions) receives a large_candidate_universe warning. The threshold does not limit generation, require --force, or predict runtime, memory, compressed size, or temporary-database size. Exhaustive windows grow rapidly: one canonical 300-aa protein produces 11,111 windows for lengths 10-50 at stride 1, so one million such proteins would produce 11,111,000,000 occurrences. Run estimate before large generation jobs and provision compute and storage from representative local measurements.

Changes from 0.1.x

GenoToPep 0.2.0 is an output-schema migration. Manifest schema 2 is not resume-compatible with schema-1 or 0.1.x output directories. Exclusion TSV files now aggregate rows by source and reason while preserving exact excluded-occurrence totals in excluded_occurrences; they no longer enumerate every excluded candidate sequence. Unique peptide FASTA records are ordered lexicographically by exact amino-acid sequence rather than first occurrence. Estimate JSON schema 1 is a separate contract from manifest schema 2. Re-run generation under 0.2.0 when schema-2 outputs are required; do not merge 0.1.x and 0.2.0 partial outputs.

Outputs

For input sample.proteins.faa.gz and strategy windows, the main outputs are:

sample.proteins.windows.peptides.faa.gz
sample.proteins.windows.peptide_occurrences.tsv.gz
sample.proteins.windows.exclusions.tsv.gz
sample.proteins.command.txt
sample.proteins.run-config.yaml
sample.proteins.run.log
sample.proteins.manifest.json
sample.proteins.checksums.sha256

FNA runs additionally write:

sample.orfs.faa.gz
sample.orfs.fna.gz
sample.orfs.gff3.gz
sample.orfs.tsv.gz

The ORF TSV protein_sequence and nucleotide_sequence columns intentionally duplicate the sequences in orfs.faa.gz and orfs.fna.gz, respectively, so each ORF row is self-contained for tabular analysis. GenoToPep 0.2.0 does not provide an option to omit these columns; account for this duplication when planning storage for very large assemblies.

Relationships among output files

• *.peptides.faa.gz contains one record per exact peptide sequence; the FASTA header is the unique peptide_id (>sha256:...).
• *.peptide_occurrences.tsv.gz joins to the peptide FASTA on peptide_id. Many occurrence rows can reference one FASTA record, so this is a many-to-one relationship from occurrences to peptides.
• For FNA input, *.orfs.tsv.gz joins id to *.peptide_occurrences.tsv.gz source_gene_id. FAA input does not produce an ORF TSV; source_gene_id then refers to the supplied protein identifier.
• *.exclusions.tsv.gz has no peptide_id. For FNA input its source_gene_id can join to ORF id, but exclusions and retained occurrence rows are mutually exclusive at the candidate-event level.

The occurrence TSV has one data row per retained candidate occurrence. The exclusion TSV instead aggregates excluded candidates by source protein, strategy, and reason: counts.exclusions equals the sum of excluded_occurrences, while counts.exclusion_groups equals the exclusion TSV data-row count. Occurrence and exclusion TSV rows must not be added as row counts. Use counts.occurrences + counts.exclusions for the complete candidate-event total; aggregated exclusions do not retain each excluded candidate sequence or coordinate.

The peptide FASTA contains one record per exact amino-acid sequence. Its ID is sha256:<digest>. Exact sequence, rather than the digest, is the internal deduplication key. The occurrence TSV retains every source, protein-relative coordinate, genomic coordinate when available, strategy, window, and cleavage occurrence.

Occurrence TSV columns

Interpret coordinate_scope before the genomic-coordinate columns. protein_and_genomic means columns parent_contig through peptide_end_nt_1based are populated from an FNA-called ORF or a successful FAA+GFF3 mapping. protein_only means genomic mapping is unavailable; for example, an empty gene_strand in an FAA run without GFF3 means unknown genomic context, not a strandless feature.

The remaining conditional groups are:

• window_length and window_stride: populated only for windows; otherwise not applicable to that generation strategy.
• cleavage_modes, elementary_fragment_count, and internal_boundaries_joined: populated only for cleavage; otherwise not applicable to that generation strategy.
• orf_caller, orf_caller_version, and min_gene_nt: populated only for FNA input; not applicable to supplied FAA input.

All integer coordinates are inclusive and 1-based. Empty values are serialized as zero-length TSV fields, not the strings None or null. An empty genomic field with coordinate_scope=protein_only means unavailable because no genomic mapping was supplied; an empty strategy- or input-specific field means not applicable. Readers such as pandas normally convert both cases to NaN, so retain coordinate_scope, generation_strategy, and source_input_type when interpreting nulls.

Exclusion TSV columns

Noncanonical exclusions are aggregated to one row per source protein, strategy, and reason. counts.exclusions remains the exact number of excluded candidate occurrences; counts.exclusion_groups is the number of written exclusion rows.

Only the canonical 20 amino acids are admitted to predictor-facing peptide FASTA. Candidates containing other symbols are written to the exclusions table with a reason.

Input normalization and identifiers

FAA sequences are converted to uppercase and trailing terminal * symbols are removed before length filtering, generation, peptide hashing, and output. FNA sequences and Pyrodigal protein translations are also converted to uppercase. The run configuration records normalization.case: upper. Internal tab, control, and DEL characters in FASTA headers are rejected because they are not safely interoperable with tabular command-line tools. FASTA identifiers must be unique within each input file; duplicates are hard errors.

Integrity, payload identity, and resume

checksums.sha256 covers the complete run inventory, including timestamped provenance files. Manifest schema version 2 additionally records role-specific SHA-256 values over the exact compressed bytes of the three deterministic scientific payload files—peptide FASTA, occurrence TSV, and exclusion TSV—and a combined payload_sha256. Scientific payloads use gzip compresslevel 6, an empty embedded filename, and mtime 0; these settings are recorded in run-config.yaml. The combined digest processes roles in lexical order using UTF-8 role + NUL + lowercase ASCII SHA-256 hex + LF framing; it excludes output-directory paths, timestamps, and command text, so equivalent payloads can be compared across output directories under the same supported compression environment.

Compressed-byte hashes can depend on the Python/zlib compression implementation. For cross-environment scientific comparison, compare decompressed payload content in addition to the recorded exact compressed-byte hashes.

These self-authored hashes detect accidental corruption and verify internal consistency; they are not authentication against deliberate modification. Anyone who can rewrite the complete output directory can replace payloads and recompute every manifest and checksum. Adversarial tamper resistance requires an external trust anchor such as a separately stored expected digest, signed manifest, or immutable repository record.

Outputs are written through an incomplete sibling directory and renamed only after successful completion. Failed runs close the temporary SQLite database and remove their incomplete directory. Concurrent GenoToPep generation attempts for the same destination are serialized by an exclusive sibling lock file; existing outputs are checked again before publication and are not intentionally overwritten. A hard process kill can bypass Python cleanup and leave a stale lock or .incomplete.<UUID> directory. After confirming that no matching generation process is running, remove only the stale artifacts belonging to the terminated invocation. The lock coordinates GenoToPep invocations; it cannot prevent an unrelated external program from racing to create the same path. --resume accepts only an identical, complete run after schema, code identity, input type, input/GFF SHA-256, parameters, normalization, compression settings, payload identity, inventory, and output checksum verification. Recorded input paths are provenance only and do not prevent resume when identical content is staged at a different path. Without --sample-id, the input basename determines output filenames; use a stable --sample-id when renamed input or GFF files must resume the same run. Resume does not continue a partially written stage.

GFF3 mapping scope

The current mapper resolves a protein FASTA identifier against GFF3 ID, protein_id, locus_tag, or Parent. Missing and ambiguous mappings are errors. Single-CDS bacterial proteins are supported; multipart spliced CDS models are not yet supported.

FAA runs without GFF3 provide protein-relative coordinates and record that genomic coordinates are unavailable. Their short-ORF recall remains limited by the upstream annotation or gene-calling workflow.

References

1. Larralde M. Pyrodigal: Python bindings and interface to Prodigal, an efficient method for gene prediction in prokaryotes. Journal of Open Source Software. 2022;7(72):4296. https://doi.org/10.21105/joss.04296
2. Hyatt D, Chen G-L, LoCascio PF, Land ML, Larimer FW, Hauser LJ. Prodigal: prokaryotic gene recognition and translation initiation site identification. BMC Bioinformatics. 2010;11:119. https://doi.org/10.1186/1471-2105-11-119
3. Keil B. Specificity of Proteolysis. Springer Berlin Heidelberg; 1992. https://doi.org/10.1007/978-3-642-48380-6
4. Gasteiger E, Hoogland C, Gattiker A, Duvaud S, Wilkins MR, Appel RD, Bairoch A. Protein Identification and Analysis Tools on the ExPASy Server. In: The Proteomics Protocols Handbook. Humana Press; 2005:571-607. https://doi.org/10.1385/1-59259-890-0:571
5. SIB Swiss Institute of Bioinformatics. PeptideCutter: cleavage specificities of selected enzymes and chemicals. https://web.expasy.org/peptide_cutter/peptidecutter_enzymes.html
6. Gross E, Witkop B. Selective cleavage of the methionyl peptide bonds in ribonuclease with cyanogen bromide. Journal of the American Chemical Society. 1961;83(6):1510-1511. https://doi.org/10.1021/ja01467a052

License

Copyright (C) 2026 Jun Won Lee.

GenoToPep is licensed under GPL-3.0-or-later. Pyrodigal is a GPL-3.0-or-later dependency. License compatibility should be reviewed before public redistribution; this statement is not legal advice.