ProteinSearch - BBTools

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 (or o=): 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. 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=t (or reducedseed=): Use amino8 reduced-alphabet seeds for higher sensitivity. Values: t (true, default) or f (false). The command-line default of t comes from the frozen v1 parameter set; the in-memory API default is false. 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=0.001 (or e=): E-value significance cutoff. Hits with E-values greater than this threshold are not reported. Default 0.001 (frozen v1 parameter set applied by the command-line front end); use larger values (e.g., 10) to be more permissive, or smaller values (e.g., 1e-5) for stringent filtering. Note: E-values are approximate and do not include edge-length correction.
minid=30 (or minpident=): Minimum percent identity to report. Valid range: 0-100. Default 30 (frozen v1 parameter set); set minid=0 to disable 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.
mincoverage=0.7 (or mincov=): Minimum bidirectional coverage to report. The alignment must cover at least this fraction of both the query and the target: alnLen/qLen >= mincoverage AND alnLen/tLen >= mincoverage, where alnLen counts gap columns. Default 0.7 (frozen v1 parameter set); 0 disables this filter. The in-memory API default is 0. This matches mmseqs --cov-mode 0.
mts=<N> (or maxtargetseqs=): Maximum number of distinct target sequences to report per query, ranked by best HSP (highest bitscore first). Default 25 (frozen v1 parameter set on the command line); the in-memory API default is unlimited (Integer.MAX_VALUE). Use this to limit output volume when searching against large databases.

Threading and Query Loading Parameters

threads=1 (or t=): Number of worker threads. Work is parallelized over queries, with exactly this many threads regardless of query count. Default 1 (single-threaded). The source documents that output is byte-identical for any thread count, so this option changes speed only.
queryloader=generic (or qloader=): How query IDs are assigned. Values: generic (default; first whitespace-delimited token of each header) or callgenes. Use callgenes when the query is raw CallGenes outa= output, where every gene on one contig shares an identical header; this mode appends a per-gene _gN suffix so IDs are unique.

Shortlist Parameters (Advanced)

shortlist=0: Number of families to align per query. Default 0 keeps the legacy behavior, in which every seed-passing candidate is aligned. When N is greater than 0, each query is scored against every family by the F4 fusion of composition, covering-set, and length terms, and only the top N families are aligned. The k-mer seed index is not built in this mode. Requires sidecar= and kmersets=.
sidecar=<file>: Precomputed per-family data file written by FamilyShortlistSidecarBuilder. Required when shortlist>0. Its hashes are verified against db= and kmersets= at load time; a stale or rebuilt sidecar is refused rather than silently scored.
kmersets=<file>: Covering-set file the sidecar was built from. Required when shortlist>0, for the hash cross-check.
aligner=aaaligner: Aligner used in shortlist mode. Values: aaaligner (default, legacy), d55 (canonical BLOSUM62 with linear gap penalty 4), glocal, or blosum (historical aliases for AAAligner.alignGlocal). Only used when shortlist>0.
kmernorm=none: Normalization of the k-mer term in the F4 score. Values: none (default; raw covering-set hits), sub (subtracts each query's expected coincidental hit count), or cs300 (hits divided by covering-set size plus 300). Only used when shortlist>0.
triage=f: Applies the v3 sidecar's per-family thresholds: a length-ratio window, a dimer-cosine floor, and a covering-k-mer floor as a survivor mask before top-N selection, plus identity, score-ratio, raw-score, and coverage thresholds inside the hit filter. Values: t or f (default). Requires a v3 sidecar (schema version 2 or higher); a v2 sidecar causes an error. Only used when shortlist>0.
covdef=column: Coverage definition used when triage=t. Values: column (default; alignment length including gap columns, divided by query and target length, matching mincoverage) or span (aligned span on each side, without gap-column credit).

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. If the query is shorter than k and yields no k-mers, every database sequence is aligned instead. When shortlist is greater than 0, the seed index is not built; candidates are chosen by the F4 family ranking described under Shortlist Parameters.

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. Filters are then applied in this order: minimum raw score, minimum percent identity, minimum coverage (when greater than 0), and E-value cutoff. The mts limit is applied last, keeping the highest-bitscore targets per query.

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

query: Query sequence identifier
target: Database sequence identifier
pident: Percent identity (0-100)
length: Alignment length (number of aligned positions)
mismatch: Number of mismatched residues
gapopen: Number of gap-opening events
qstart: Query sequence start coordinate (1-based)
qend: Query sequence end coordinate (1-based, inclusive)
tstart: Target sequence start coordinate (1-based)
tend: Target sequence end coordinate (1-based, inclusive)
evalue: Approximate E-value (not edge-corrected)
bitscore: Normalized bit score (rigorous, gapped BLOSUM62)

Metadata Sidecar (.meta)

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

Query and database filenames, SHA-256 digests of both files, and the query loader used
Scoring matrix (BLOSUM62) and gap penalties (11/1)
Statistical parameters (lambda, K from Blosum62)
Actual k-mer size, reduced-alphabet flag, and min-seed-hits threshold
Applied significance cutoffs (E-value, min identity, min coverage, min raw score, max-target-seqs)
A flag indicating E-values are approximate
Counts: queries, database sequences, and total hits
Shortlist settings (size, aligner, k-mer normalization, triage) and per-stage triage counters when shortlist is greater than 0

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. Set threads to parallelize over queries; the default is 1.

Limitations (Phase-1 MVP)

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

Indexed fast search (rebuilds index every run)
Edge-corrected E-values (current implementation is approximate)
Multiple HSPs per query-target pair (reports best HSP only)
BlastX or translated search (proteins only)

Troubleshooting

Empty Results

If no hits are reported, check:

Are your query and database sequences valid FASTA with protein characters?
Is the E-value cutoff too stringent? Try evalue=1000 temporarily to see if anything aligns.
Are the default filters too stringent for your sequences? The command-line defaults are minid=30 and mincoverage=0.7. Try minid=0 mincoverage=0 temporarily to see if anything aligns.
Are k-mers too long (k=5 may be too specific for distant sequences)? Try k=3 or reduced=t.

Out of Memory

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

Check the size of your database. Very large databases (e.g., GenBank nr) may exceed available RAM.
Increase Java heap with -Xmx8g (or more if available) before running.
Consider splitting your database into smaller chunks.

Slow Execution

If the tool runs slowly:

Many candidates per query slow down alignment scoring. Try stricter filtering: k=6 (more specific), reduced=f, or higher minseedhits.
Permissive E-value cutoffs produce large output. The default is evalue=0.001; keep it or lower it.

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:

Email: bbushnell@lbl.gov
Documentation: bbmap.org
GitHub: BBTools repository