ProteinSearch

Script: proteinsearch.sh Package: prot Class: ProteinSearch.java

Blastp-style protein similarity search that reads query and database protein FASTA files and writes BLAST-tab (outfmt 6) TSV results. Seeds candidates with exact or amino8 reduced-alphabet k-mers, then scores with BLOSUM62 affine-gap local alignment (gap-open 11, gap-extend 1). This is a Phase-1 MVP implementation providing rigorous bitscores and approximate E-values.

Basic Usage

proteinsearch.sh query=<query.faa> db=<database.faa> out=<hits.tsv>

At minimum, you must provide a query protein FASTA file and a database protein FASTA file. The tool will search all query sequences against the database and write results as a TSV file in BLAST-tab outfmt 6 format (11 columns: query, target, pident, length, mismatch, gapopen, qstart, qend, tstart, tend, evalue, bitscore).

Output Files

  • Main output: The TSV file specified by out=, or stdout if not specified
  • Metadata sidecar: A .meta file written alongside the output recording parameters, scoring matrix details, counts, and a flag indicating that E-values are approximate (edge-length correction omitted)

Coordinate System

All coordinates in the output are 1-based inclusive, following BLAST convention. Query start/end and target start/end indicate the aligned region in each sequence.

Scoring Details

Bitscores are rigorously calculated using the gapped BLOSUM62 matrix with gap-open penalty 11 and gap-extend penalty 1. E-values are approximate and do not include Karlin-Altschul edge-length correction; this limitation is documented in the .meta sidecar file.

Parameters

Required Parameters

query=<file> (or in=, q=)
Query protein FASTA file (amino acids). Must contain valid protein sequences in FASTA format. The identifier is the first whitespace-delimited token of each header.
db=<file> (or ref=, database=, d=)
Database protein FASTA file (amino acids). Must contain valid protein sequences. Identifiers must be unique within the database.

Output Parameters

out=stdout
Output file for TSV results (outfmt 6). Default is stdout (write to console). When a file path is specified, a .meta metadata sidecar is automatically written alongside it. Use out=stdout to disable file output.
ow=t (or overwrite=)
Overwrite existing output files. Values: t (true, default) or f (false). If false and the output file already exists, an error is raised.

Seeding Parameters

k=5
Seed k-mer length for candidate selection. Valid range: 1-15. Longer k-mers are more specific (fewer false candidates) but may miss true homologs with high divergence. Default 5 provides a good balance for most protein searches.
reduced=f (or reducedseed=)
Use amino8 reduced-alphabet seeds for higher sensitivity. Values: t (true) or f (false, default). The amino8 alphabet collapses chemically similar residues into 8 groups, allowing seeds to match divergent homologs at reduced specificity. Useful for remote homology detection.
minseedhits=1
Minimum number of distinct shared k-mers required for a target to be considered a candidate for alignment. A target must share at least this many distinct k-mers with the query to proceed to the scoring phase. Default 1 means any k-mer match qualifies. Higher values reduce false candidates but may miss weak homologs.

Filtering Parameters

evalue=10
E-value significance cutoff. Hits with E-values greater than this threshold are not reported. Default 10 is permissive; use smaller values (e.g., 1e-5) for stringent filtering. Note: E-values are approximate and do not include edge-length correction.
minid=0 (or minpident=)
Minimum percent identity to report. Valid range: 0-100. Default 0 disables this filter. Specified as a percentage (e.g., minid=30 means 30% identity or better).
minscore=0
Minimum raw BLOSUM62 alignment score to report. Default 0 disables this filter. Raw score is the sum of residue match/mismatch scores plus gap penalties before normalization to bits.
mts=<N> (or maxtargetseqs=)
Maximum number of distinct target sequences to report per query, ranked by best HSP (highest bitscore first). Default is unlimited (Integer.MAX_VALUE). Use this to limit output volume when searching against large databases.

Java Runtime Parameters

-Xmx<memory>
Set maximum Java heap memory (e.g., -Xmx4g for 4 gigabytes). By default, the wrapper script auto-detects available system memory and sets reasonable defaults.
-Xms<memory>
Set initial Java heap memory. Usually not needed; auto-detected by the wrapper script.
-eoom
Exit immediately on out-of-memory exception, rather than attempting recovery. Useful for scripted pipelines that need clean failure modes.
-da
Disable Java assertions. Not recommended for normal use; assertions help catch internal logic errors.

Algorithm

Search Strategy

ProteinSearch uses a two-phase approach: seed-and-extend. First, k-mers extracted from the query are matched against a hash index of k-mers from all database sequences. A target sequence becomes a candidate if it shares at least minseedhits distinct k-mers with the query. For each candidate, a Smith-Waterman local alignment is computed using BLOSUM62 scores with affine gaps.

Reduced-Alphabet Seeding

When reduced=t, k-mers are encoded using the amino8 reduced alphabet (groups of chemically similar amino acids map to the same code). This allows seeds to match more divergent sequences, improving sensitivity for remote homology at the cost of more candidates to score. Standard seeding uses 5 bits per residue (20 standard amino acids); reduced seeding uses 3 bits per residue (8 groups).

Alignment and Scoring

Each query-target pair proceeds to affine-gap local Smith-Waterman alignment. The BLOSUM62 substitution matrix is used with fixed penalties: gap-open 11, gap-extend 1. Only the best scoring HSP (local alignment) per query-target pair is retained. Scores are converted to bitscores using rigorous BLOSUM62 statistics.

E-value Calculation

E-values are calculated using Karlin-Altschul statistics without edge-length correction. This makes the E-values approximate and potentially overestimated at the margins of the search space. The sidecar .meta file explicitly flags that E-values are approximate. For publication-quality work, consider these E-values indicative only; more rigorous statistics are planned for future releases.

Output Ordering

Hits are sorted in a deterministic total order: query (ascending), E-value (ascending), bitscore (descending), target (ascending), target start (ascending), query start (ascending). This ensures reproducible results across runs.

Examples

Basic Search

Search a set of query sequences against a protein database, writing results to a file:

proteinsearch.sh query=query.faa db=uniref50.faa out=hits.tsv

This will produce hits.tsv (the search results) and hits.tsv.meta (metadata about the run).

Sensitive Remote Homology Search

Use reduced-alphabet seeding and a relaxed E-value cutoff to find more distant homologs:

proteinsearch.sh query=query.faa db=database.faa out=remote_hits.tsv reduced=t evalue=1 k=4

The combination of reduced=t (amino8 alphabet for seeds), smaller k=4 (shorter seeds for more matches), and higher evalue=1 makes the search more permissive, useful for detecting weakly conserved homologs.

Stringent High-Identity Search

Find only close homologs with high sequence identity:

proteinsearch.sh query=query.faa db=database.faa out=close_hits.tsv evalue=1e-10 minid=70 minscore=100

This restricts output to hits with E-value ≤ 1e-10, percent identity ≥ 70%, and raw alignment score ≥ 100, filtering for closely related sequences only.

Limited Output with max-target-seqs

Return only the top 5 matching targets per query to manage output volume on large databases:

proteinsearch.sh query=query.faa db=nr.faa out=top5.tsv mts=5

Each query will report at most 5 targets, ranked by best bitscore. Weaker matches are discarded.

To stdout (streaming)

Write results to standard output for piping to other tools:

proteinsearch.sh query=query.faa db=database.faa out=stdout | head -20

No .meta file is created when writing to stdout.

Output Format

Main Results (TSV)

Results are in BLAST-tab outfmt 6 format: 11 tab-separated columns per line.

query target pident length mismatch gapopen qstart qend tstart tend evalue bitscore

Metadata Sidecar (.meta)

A .meta file is written alongside the TSV output (when writing to a file, not stdout). It documents:

Performance Notes

Memory Usage

Memory is dominated by the target k-mer index, which stores one hash table entry per distinct k-mer in the database. For typical protein databases with reduced redundancy, expect approximately 1-2 MB per million residues. The script auto-detects available system memory and allocates up to 2 gigabytes by default; use -Xmx to override.

Runtime

Runtime depends on query size, database size, number of candidates per query, and sequence complexity. Seed-based filtering typically reduces candidates to 1-5% of the database, greatly accelerating alignment scoring. Short queries or highly selective k-mer parameters may produce fewer candidates and run faster.

Limitations (Phase-1 MVP)

This is a correctness-first MVP implementation. The following features are deferred:

Troubleshooting

Empty Results

If no hits are reported, check:

Out of Memory

If the tool fails with an out-of-memory error:

Slow Execution

If the tool runs slowly:

Workflow Integration

ProteinSearch is designed as a standalone tool but integrates naturally into larger pipelines:

Metagenomic Protein Binning

Search predicted proteins from metagenomic contigs against a curated marker gene set to assign taxonomy. ProteinSearcher (the in-memory API) is callable directly from Java code, allowing binning algorithms to assess in-memory proteins without disk round-trips.

Comparative Genomics

Search one organism's proteome against another to identify orthologs, paralogs, and lineage-specific proteins.

Functional Annotation

Search novel proteins against well-characterized databases (e.g., COG, Pfam) to infer function from homology.

Support

For questions and support: