ClusterProteins

Script: clusterproteins.sh Package: prot Class: ClusterProteins.java

Greedy identity-threshold protein clustering following the CD-HIT and linclust approach. Reads a protein FASTA file and groups sequences into clusters based on identity and coverage thresholds using BLOSUM62 affine-gap local alignment. Sequences are processed longest-first; each new sequence is seeded against existing representatives using shared k-mers, aligned with BLOSUM62, and assigned to the best matching representative if it meets both identity and coverage requirements, or starts a new cluster.

Basic Usage

clusterproteins.sh in=<proteins.faa> out=<clusters.tsv>

ClusterProteins reads amino acid sequences from a FASTA file and produces a representative-to-member TSV output with one row per member. A metadata sidecar file (with .meta extension) is automatically written alongside the output file, recording run parameters and cluster statistics.

Input Format

Input must be a protein FASTA file (amino acids only). Standard FASTA format: header line starting with '>' followed by sequence identifier, then amino acid sequences on subsequent lines. Sequences with ambiguous residues (X) are handled gracefully — ambiguous positions interrupt k-mer extraction but do not fail processing.

Output Format

Output is a tab-separated values (TSV) file with the following columns:

cluster_id: Integer identifier for the cluster (assigned in order of representative creation)
representative: Sequence identifier of the cluster representative
member: Sequence identifier of this member
identity: Percent identity (0-100) between member and representative alignment
coverage: Fraction (0-1) of residues spanned by the alignment for both member and representative
is_representative: 1 if this sequence is the cluster representative, 0 otherwise

Important Notes

Every sequence is represented once as a member (representatives appear as members of their own cluster with identity=100.000, coverage=1.000)
Cluster IDs are stable within a single run only; cross-run stable identifiers are not implemented
The .meta sidecar file contains run parameters, gap penalties, counts, and cluster statistics

Parameters

ClusterProteins supports both required and optional parameters. Optional parameters control clustering stringency, seeding strategy, and output behavior. Identity thresholds accept either percent (e.g., 90) or fraction (e.g., 0.9) format.

Required Parameters

in=<file>: (input, i) Input protein FASTA file (amino acids). This parameter is required.

Output Parameters

out=stdout: (o) Output TSV file path. Default: stdout (write to standard output). A .meta sidecar manifest is written beside the output file (when out is not stdout). Default: stdout
ow=t: (overwrite) Overwrite existing output files (t/f). Set to 't' to overwrite, 'f' to fail if output already exists. Default: t

Clustering Thresholds

minid=90: (minidentity, id) Minimum percent identity threshold to join a cluster [0-100]. Also accepts fractional format (e.g., 0.9 is equivalent to 90). Default: 90
mincov=0.8: (mincoverage, cov) Minimum aligned fraction for both member and representative [0-1]. Both the member's coverage (aligned length / member length) and the representative's coverage must meet this threshold. Default: 0.8

K-mer Seeding Parameters

k=5: Seed k-mer length [1-15]. Smaller k values (3-5) are more sensitive but slower; larger values (7-10) reduce candidates at the cost of sensitivity. Default: 5
minseedhits=1: Minimum number of distinct shared k-mers required to align a candidate representative against a query sequence. This parameter gates alignment attempts to avoid unnecessary all-vs-all comparisons. Default: 1
reduced=f: (reducedseed) Use amino8 reduced-alphabet seeds for k-mer generation (t/f). When enabled, uses a more sensitive amino acid grouping (8 groups) instead of the standard 20-letter alphabet. This increases seed sharing and sensitivity at the cost of specificity. Default: f

Common Use Cases

Standard Protein Clustering at 90% Identity

Cluster a set of protein sequences at the default 90% identity threshold with default coverage and seeding parameters. This is the most common workflow for grouping proteins by similarity.

clusterproteins.sh in=proteins.faa out=clusters.tsv

This command clusters proteins with 90% minimum identity, 80% minimum coverage, k-mer length 5, and default seeding strategy. A clusters.tsv.meta file is written with run parameters.

Stringent Clustering at 95% Identity

For applications requiring very high sequence similarity (e.g., identifying near-duplicates or exact variants), increase the identity threshold.

clusterproteins.sh in=proteins.faa out=clusters_strict.tsv minid=95 mincov=0.9

This command uses 95% identity and 90% coverage thresholds, producing fewer clusters with more similar members.

Sensitive Clustering with Reduced Alphabet

For divergent protein families where standard exact k-mer matching misses related sequences, use the amino8 reduced alphabet to improve seed sharing between distant homologs.

clusterproteins.sh in=proteins.faa out=clusters_sensitive.tsv minid=70 k=4 reduced=t

This command uses more sensitive seeding (amino8 reduced alphabet with k=4) and relaxed identity (70%), allowing clustering of more distant homologous sequences. Smaller k values and reduced alphabet increase overlap likelihood.

High-Coverage Clustering to Avoid Fragment Matches

When sequences vary in length and you want to avoid spurious clustering of fragments to full-length proteins, require high coverage from both sides.

clusterproteins.sh in=proteins.faa out=clusters_full.tsv minid=85 mincov=0.95

This command requires that both member and representative have at least 95% coverage of their respective lengths, effectively filtering out fragment-to-full-length matches.

Algorithm Details

Clustering Strategy

ClusterProteins uses a greedy, deterministic clustering approach:

Ordering: Input sequences are sorted longest-first (ties broken by sequence identifier), following the CD-HIT convention. Longer sequences become representatives, which reduces fragmentation and improves biological coherence in many cases.
Seeding: A growing hash map index tracks which representatives contain each k-mer. For each new sequence, only representatives sharing at least minseedhits distinct k-mers are considered as candidates, avoiding expensive all-vs-all alignment.
Alignment: Each candidate representative is aligned to the query sequence using BLOSUM62 with affine gap penalties (gap-open=11, gap-extend=1).
Assignment: The sequence joins the best-scoring candidate (highest identity) that meets both identity and coverage thresholds; otherwise it becomes a new representative.

Identity and Coverage

Identity is the percent similarity in aligned columns (matches/aligned_length). Both member and representative must have at least mincov fraction of their respective lengths covered by the alignment to be considered a valid cluster member. This two-sided coverage requirement prevents spurious matches between full-length and fragment sequences.

K-mer Seeding and Reduced Alphabets

Standard k-mer seeding uses the full 20-letter amino acid alphabet; reduced alphabet seeds (amino8) group similar residues (e.g., I, L, M, V → 1 group), increasing seed overlap for distant homologs while reducing specificity. K-mer positions containing X (ambiguous residues) are skipped without failing the sequence.

Output Stability

Cluster IDs are assigned in order of representative creation and are stable within a run given identical input and parameters. Cross-run stability (resuming a run with new sequences) is not implemented; re-clustering full input will produce different IDs.

Alignment Parameters

All alignments use BLOSUM62 matrix with the following fixed parameters (not currently user-configurable):

Matrix: BLOSUM62
Alignment mode: Local (Smith-Waterman)
Gap open penalty: 11
Gap extend penalty: 1

Output Files

Main Output (TSV)

The primary output file contains tab-separated columns: cluster_id, representative, member, identity, coverage, is_representative. One row per input sequence.

Metadata Sidecar (.meta)

When output is written to a file (not stdout), a .meta file is automatically created (e.g., clusters.tsv.meta). This manifest records:

Input FASTA path
Alignment matrix and gap penalties
All clustering parameters (k, reduced_seed, min_seed_hits, minid, mincov)
Cross-run stability status (currently: false)
Total sequences and clusters formed
Singleton cluster count
Largest cluster size

Memory and Performance

ClusterProteins stores all sequences and cluster data in memory. Memory usage scales with input size and the diversity of sequences (more diverse input means more clusters and larger seed index). For typical proteomes (thousands to tens of thousands of sequences), memory usage is typically under 1-2 GB. The seed index is a hash map of (k-mer → list of representative cluster indices), so smaller k or higher redundancy increases index size proportionally.

Runtime scales approximately O(n log n) for sorting plus O(n × candidates × alignment_cost), where candidates depends on seed overlap (smaller k = more candidates = slower). Setting larger k or higher minseedhits reduces candidate count at the cost of sensitivity.

Performance Tuning

Use larger k (7-8) for faster clustering on large datasets, at the cost of sensitivity
Increase minseedhits (2-3) to reduce alignment attempts, trading sensitivity for speed
For highly redundant input, larger k and higher minseedhits have minimal sensitivity cost
For divergent families, reduce k to 3-4 or enable reduced=t to improve seed sharing

Java Parameters

ClusterProteins inherits standard Java command-line options that can be passed through the shell script:

-Xmx<size>: Set maximum heap memory (e.g., -Xmx4g for 4GB). By default, memory is auto-detected.
-eoom: Exit cleanly on out-of-memory errors instead of crashing.
-da: Disable Java assertions (not recommended for debugging).

Support

For questions, bug reports, or feature requests: