index

Structure search and annotation

Domainator can use protein structures as input, either annotating them with reference structures or selecting matching records from a structure database. Both produce ordinary Domainator genbank files, so the rest of the suite — domainator extract_domains, domainator plot_contigs, domainator enum_report, domainator summary_report, domainator color_genbank, domainator build_ssn, domainator seq_dist, domainator select_by_contig — works on the results unchanged. domainator structure_dist is the exception: it writes a data matrix rather than genbank, feeding domainator build_ssn, domainator build_tree and domainator build_projection the way domainator seq_dist does.

There are four roles, mirroring the sequence-based tools:

	sequence input	structure input
annotate (all records, Domainator features)	domainator domainate	domainator structure_domainate
search (hit records, Domain_Search features)	domainator domain_search	domainator structure_search
compare (pairwise score/distance matrix)	domainator seq_dist	domainator structure_dist
convert	—	domainator structure_to_genbank

Command roles and output shapes

Annotation and search both perform similarity searches against references. Choose between them by the required output, not merely the input size or reference count:

Command	Records retained and output shape	Annotation features
domainate	All input sequence records, including those without qualifying hits; no hit-centered extraction.	Domain annotations (Domainator) for retained hits.
structure_domainate	All input structure-derived protein records, including chains without qualifying hits; full observed sequences.	Domain annotations (Domainator) for retained hits.
domain_search	Matching protein records, nucleotide hit regions/neighborhoods, or whole contigs; protein-reference CDS hits can be translated with --translate.	Search best-hit annotations (Domain_Search); --add_annotations also adds domain annotations (Domainator).
structure_search	Matching structure-derived protein records, with full observed target sequences rather than alignment crops.	Search best-hit annotations (Domain_Search); --add_annotations also adds domain annotations (Domainator).

A best-hit feature marks an aligned region; its location is not necessarily the extent of the emitted sequence. Structure coordinates refer to observed residues, not original PDB residue numbering or residues absent from the structure. Additional domain annotations are subject to the tool's annotation filters and limits; --add_annotations does not mean that every raw backend hit is written.

In all four commands, -i identifies the records/database being annotated or searched, and -r supplies references being searched for. For the structural aligners, those references are backend queries and the input database contains backend targets; these alignment terms do not change the CLI convention. Search tools are optimized for large databases with relatively few references, but that performance guidance is not the definition of the annotation/search roles. Other commands use their own output contracts: hmmer_search, for example, selects HMM profiles rather than sequence neighborhoods.

Requirements

The structure tools shell out to a structural aligner. Two backends are supported, selected with --algorithm:

backend	install	notes
foldseek (default)	installed with the conda environment	full feature set
reseek	build from source, see below	no alignment types or structural metrics

Neither needs PyTorch, esmologs, or the ESM-2 3B checkpoint. That stack is only needed for annotating sequences that have no structure.

Check availability with foldseek version / reseek, or pass an explicit path with --foldseek_path / --reseek_path.

Building reseek

reseek is not on bioconda yet, and two things make the published v3.01 artifacts unusable on many machines, so it has to be built from source:

• The precompiled reseek-v3.01-linux-x86 release binary is compiled with -march=native on a machine with AVX-512. On any CPU without AVX-512 it dies immediately with Illegal instruction (SIGILL) before printing anything.
• The v3.01 source tag does not compile with GCC 13 or later: it copy-initializes std::atomic (atomic<uint> m_NextQueryIdx = 0;), whose copy constructor is deleted.

Both are fixed in commits after the v3.01 tag, so build from the default branch:

git clone https://github.com/rcedgar/reseek.git
cd reseek/src
bash gitver.bash          # generates git_hash.h, which the checked-in Makefile does not
make -j$(nproc) \
    CFLAGS="-flto -ffast-math -march=x86-64-v3 -O3 -DNDEBUG" \
    CPPFLAGS="-flto -ffast-math -march=x86-64-v3 -O3 -DNDEBUG" \
    LDFLAGS="-flto -ffast-math -march=x86-64-v3 -O3 -static"
cp ../bin/reseek /somewhere/on/your/PATH

Two notes on that command. bash gitver.bash is required because the checked-in Makefile includes git_hash.h but has no rule to generate it — only reseek's own vcxproj_make.py build path wires the two together, so a plain make fails with fatal error: git_hash.h: No such file or directory. And the -march=x86-64-v3 override replaces the Makefile's -march=native; x86-64-v3 is what upstream switched its default to after v3.01 "for wider CPU compatibility", and it requires only AVX2, so the resulting binary runs on anything from Haswell onward.

Note that the binary reports itself as reseek v3.0 even when built from the v3.01 tag or later — the internal version string was not bumped. The bracketed git hash in the banner is what actually identifies the build.

Accepted inputs

Both -i and -r accept, in any of the four tools:

• individual structure files — .pdb, .cif, .mmcif, .ent, .bcif, optionally gzipped,
• a directory, searched recursively for structure files,
• a glob such as structures/*.pdb.gz (expanded internally, so it also works from a --config file, where the shell would not expand it),
• the prefix of a prebuilt foldseek database.

One record is written per structure chain. Chain splitting, structure parsing, and entry naming are all done by foldseek, so records are named exactly as they appear in a hand-run foldseek search: <file basename>_<chain>, for example 1ubq_A.

Annotating structures

domainator structure_domainate -i my_structures/ -r reference_structures/ -o annotated.gb
domainator plot_contigs -i annotated.gb --html_out annotated.html

Every input chain is annotated and written, whether or not it had a hit. Use structure_search to keep only matching chains; add --add_annotations to add domain annotations (Domainator features) alongside the search best-hit annotation (Domain_Search feature) already present. structure_domainate no longer accepts --hits_only or --no_annotations, in either its CLI or its public Python API.

Because both databases are built from real structures they carry C-alpha coordinates, so the full range of foldseek alignment modes is available:

domainator structure_domainate -i my_structures/ -r reference_structures/ \
    --alignment_type 1 --metrics tmscore lddt rmsd -o annotated.gb

--alignment_type is 0 for 3Di only, 1 for TM-align, and 2 (the default) for 3Di+AA. --metrics adds tmscore, lddt, rmsd and prob qualifiers to the annotations. tmscore, lddt and rmsd need coordinates in both databases; asking for them against a database built from sequences rather than structures is an error rather than a silent degradation.

Searching a large structure database

domainator structure_search -i /path/to/afdb_foldseek_db -r query.pdb \
    --max_hits 1000 --max_seqs 1000000 --device cuda:0 --cpu 6 -o hits.gb
domainator enum_report -i hits.gb --columns domains -o hits.tsv

Only records that matched are written, each containing the full observed target sequence, not a crop to its alignment. A search best-hit annotation (Domain_Search feature) marks the region homologous to its best-matching reference. --add_annotations also writes domain annotations (Domainator features) for retained hits.

Hit records are built from the sequences foldseek reports for each hit, so the searched database is never read end to end and its size does not bound memory usage. Hits are grouped by database record using an external sort, so memory scales with the hits for one record rather than with the whole result set.

Two knobs matter at scale:

• --max_seqs bounds how many hits are kept per reference. foldseek defaults to 1000, which is low for a database of hundreds of millions of structures; reseek applies no cap unless asked. A warning is printed whenever a reference saturates the cap, so truncation is never silent.
• --max_hits bounds how many records are returned, keeping the best-scoring ones. It applies to the whole search, not per reference, and makes the output sorted best-first and written after the search completes.

Comparing structures to each other

domainator structure_dist is the structure-input counterpart of domainator seq_dist: an all-vs-all structural alignment whose output is a data matrix rather than annotations.

domainator structure_dist -i my_structures/ --mode tmscore --alignment_type 1 --hdf_out tm.hdf5
domainator build_tree -i tm.hdf5 --newick_out tree.nwk

Omitting -r compares the input against itself, which is the usual case; supplying it gives a rectangular matrix of inputs against references. Every input chain is a row whether or not it aligned to anything, and the row and column labels are exactly the record ids domainator structure_to_genbank writes for the same structures, so a matrix and a genbank of the same inputs line up.

Two things differ from the other structure tools:

• The alignment direction is reversed. Here the input is the query and the reference is the target, so rows are inputs and --evalue is computed against the reference set — the same meaning it has in domainator seq_dist, where diamond runs with the input as -q.
• The whole input database is read, since a structure with no hits still needs its row. domainator structure_search avoids that, and remains the right tool for finding a handful of matches in a very large database.

--mode chooses what the cells hold. The score-based modes behave exactly as they do in domainator seq_dist, and the structural modes — tmscore, lddt, fident — put a similarity already bounded in [0, 1] in the cell, so their _dist forms are simply 1 - value and need none of the min(row_max, col_max) normalization score_dist does. Because an absent pair in a distance matrix means "maximally distant" rather than zero, the _dist modes are dense-only.

foldseek scores each direction of a pair separately, so the raw matrix is slightly asymmetric. --symmetrize reconciles the two directions, and the choice matters when a pair was found in only one direction, because the missing direction is a stored zero rather than a missing value: max keeps such a pair (the convention domainator build_ssn, --mst, and --knn use), min deletes it, making it a mutual-hit filter, and mean halves it. --mst and --knn prune the matrix as it is built and are always max-symmetric.

Reusing a database

Building a structure database is the expensive part of a run, so it can be persisted and reused:

# build the reference database once
domainator structure_to_genbank -i reference_structures/ --db_out refdb -o refs.gb

# reuse it, and keep the input database too
domainator structure_domainate -i my_structures/ -r refdb --db_out mydb -o annotated.gb
domainator structure_search    -i mydb -r refdb -o hits.gb

A database written by --db_out carries C-alpha coordinates, so it supports every alignment mode and metric. Databases built from amino acid and 3Di sequences instead — including the ones described in the ESM-2 3B 3Di documentation — do not, and are limited to --alignment_type 0 and 2.

Using structures with the sequence-based tools

domainator structure_to_genbank converts structures without searching, which makes them ordinary Domainator inputs:

domainator structure_to_genbank -i my_structures/ -o chains.gb

domainator domainate -i chains.gb -r pfam.hmm -o annotated.gb        # annotate with HMM profiles
domainator domain_search -i chains.gb -r query.hmm -o hits.gb        # search them by sequence
domainator seq_dist -i chains.gb -r chains.gb --hdf_out distances.hdf5   # by sequence, not fold
domainator deduplicate_genbank -i chains.gb -o nonredundant.gb

The three genbank-producing structure tools write to stdout with -o -, so they can be piped into any tool that reads genbank from stdin with -i - (domainator structure_dist writes matrix files through --hdf_out and --dense_text_out; --dense_text_out - writes to stdout):

domainator structure_to_genbank -i my_structures/ -o - | domainator color_genbank --target domains -o colored.gb

Note that domainator domainate and domainator domain_search read from files rather than stdin, so write an intermediate file for those (as above).

The domainator seq_dist line above compares the amino acid sequences foldseek extracted from the structures, not the structures themselves, so it will miss remote homologs that share a fold but little sequence. Use domainator structure_dist when you want the fold compared.

--store_3di additionally records each chain's 3Di structural-alphabet string in a /threedi qualifier, which is useful for rebuilding a foldseek database from a genbank file. Note that it is one opaque string per record, so it does not survive tools that slice records.

How the reseek backend differs

reseek is an independent implementation, and the interface reflects what it can actually do rather than pretending the two backends are interchangeable. Requesting something reseek cannot do is an error, not a silent no-op.

	foldseek	reseek
--alignment_type	0/1/2, default 2	not supported (no such concept)
--metrics	tmscore, lddt, rmsd, prob	none
--mode (domainator structure_dist)	every mode	score-based modes plus fident/fident_dist; tmscore, lddt and efi_score rejected
--store_3di	supported	not supported
database	a set of files sharing a prefix	one .bcb file
minimum chain length	none	32 residues; shorter chains are dropped
provenance in record descriptions	yes	no source map available
sensitivity knobs	--alignment_type, --foldseek_args	--reseek_sensitivity, --reseek_stats
--max_seqs	bounds the search, default 1000	bounds the output only, no cap by default

Two reseek-specific options:

• --reseek_sensitivity {fast,sensitive} (default sensitive).
• --reseek_stats {family,superfamily,fold} (default superfamily) — the SCOP level reseek's p-value model is calibrated against. This is not cosmetic: it shifts p-values by orders of magnitude and therefore changes which hits pass --evalue. superfamily is the usual remote-homology setting.

--max_seqs bounds output, not work. foldseek's --max-seqs bounds its prefilter, so raising it costs time and improves coverage. reseek has no equivalent option and always does the full search, so Domainator applies the cap itself, keeping the most significant hits per reference in a single streaming pass over reseek's output. This matters because reseek's output is not ordered by significance — taking the first N lines would discard the best hits. The knob that actually limits reseek's work is --evalue. There is no cap by default.

Short chains. reseek drops any chain shorter than 32 residues when building its database, so those structures are simply absent from the output. A warning is printed whenever this happens, since otherwise the records would vanish with no explanation.

E-values. reseek reports a per-comparison p-value, not an E-value. Its own evalue output column is a rescaling against a fixed SCOP40 size (evalue = pvalue * 8290) and so ignores the database actually being searched, which would break the invariant described below. Domainator therefore derives the E-value itself from the p-value and the real database size, E = p * N, and converts the --evalue threshold the same way. That keeps --evalue meaning exactly what it means for foldseek. N is read from the .bcb header, so this costs nothing.

Scores. reseek reports no bitscore, so the /score qualifier holds -log10(p). It is monotonic in significance, which is what the overlap filtering needs, but it is not comparable to a foldseek bitscore — don't mix annotations from the two backends and then filter on score.

Choosing an alignment direction, and what --evalue means

References are always aligned as the query and the input as the target. This is what foldseek is optimized for, and it means --evalue has the same meaning it has everywhere else in Domainator: it is computed against the database being searched or annotated. Running it the other way would make --evalue depend on how many reference structures you happened to supply.

One consequence worth knowing: E-values shift when the input database changes size, so the same reference against the same structure will score differently in a 9-record and a 9-million-record database. This is the same behavior as hmmsearch, where -Z exists to normalize it; foldseek has no equivalent knob.

Annotation coordinates are positions on the input record. The matching positions on the reference are recorded in the /rstart, /rend and /rlen qualifiers. See File Formats for the full qualifier list.

--min_evalue excludes hits at or below a given E-value, which is the usual way to drop a reference's match to its own copy inside the database. foldseek has no equivalent option, so Domainator applies this bound itself.

Cross-checking against foldseek directly

--hits_tsv_out writes the aligner's own hit table alongside the genbank output, which makes the results directly comparable to a hand-run search. Note that with reseek this is the table as reseek produced it, before Domainator applies --max_seqs, so it can contain more rows than the genbank output has records.

foldseek easy-search refs/*.pdb $db ref.tsv tmp --format-mode 4 -e 0.001 \
    --format-output query,target,fident,alnlen,qstart,qend,tstart,tend,evalue,bits,qlen,tlen
domainator structure_search -i $db -r refs/ -e 0.001 --hits_tsv_out ours.tsv -o hits.gb

Passing options through to foldseek

--foldseek_args (or --reseek_args) appends options to the aligner's search command line, for options Domainator does not wrap. Give them as one quoted string, which is split like a shell command line, and use the = form so they are not mistaken for Domainator's own options:

domainator structure_search -i $db -r refs/ --foldseek_args="--prefilter-mode 2" -o hits.gb

Options that Domainator sets itself (for example -e or --max-seqs) are refused, so use Domainator's own option for them instead.

Limitations

• Residue numbering. Coordinates are indices into the observed residue sequence foldseek extracts from the structure, not original PDB residue numbers. A chain with a disordered loop has a sequence shorter than its residue-number span, so hit positions cannot be mapped back to the structure's own numbering.
• No coordinates in genbank. Structure coordinates are deliberately not stored in the genbank output — they would inflate records several-fold, exceed what standard genbank parsers accept on one line, and would be silently wrong after any tool that slices records. Persist the aligner database with --db_out instead.
• Backend-specific gaps. See the comparison table above: reseek has no alignment types, no structural metrics, no 3Di output, no provenance, and a 32-residue minimum chain length.
• LoLalign is not in any foldseek release yet. When it lands, --foldseek_args should make it reachable without a Domainator release.
• Scores are not comparable across backends. foldseek reports a bitscore and reseek reports -log10(p). Annotations from the two are individually meaningful but must not be pooled and then filtered or sorted on /score.