mgnifam

Iterative HMM-based protein family generation over very large sequence databases.

Given a chunk of MMseqs2 clusters and a protein FASTA, mgnifam generate_families builds an HMM from each cluster, recruits new members from the whole database, re-aligns, and either converges on a family or discards the cluster. It is the core algorithm of the mgnifams Nextflow pipeline, extracted into a standalone, tested package.

Install

pip install mgnifam        # or: uv tool install mgnifam

Or from bioconda, into its own environment:

conda create -n mgnifam -c conda-forge -c bioconda mgnifam
conda activate mgnifam

Channel order matters — put conda-forge before bioconda, per the bioconda setup. The package is noarch, and conda pulls in a Python 3.13 interpreter itself, so it does not have to be the one already on your PATH. mamba/micromamba work the same way with the same flags.

Requires Python >= 3.13. Verify with mgnifam --version.

To work on the package itself, or to reproduce published results byte-for-byte, install from the repository against the committed lockfile instead — see Reproducibility, which is scoped to that resolved dependency set:

git clone https://github.com/vagkaratzas/mgnifam && cd mgnifam
uv sync --frozen

Usage

Commands below are written uv run mgnifam ... for the cloned checkout. On a pip or conda install, drop the uv run prefix.

uv run mgnifam generate_families \
    --clusters_chunk clusters.tsv \
    --fasta_file mgnifams_input.fa

Only those two are required. Every other flag defaults to the value below, so the run above is equivalent to spelling all of them out:

uv run mgnifam generate_families \
    --clusters_chunk clusters.tsv \
    --fasta_file mgnifams_input.fa \
    --output_dir output \
    --cpus 8 \
    --chunk_id 1 \
    --discard_min_rep_length 75 \
    --discard_max_rep_length 2000 \
    --discard_min_starting_membership 0.9 \
    --max_seq_identity 0.8 \
    --max_seed_seqs 2000 \
    --max_gap_occupancy 0.5 \
    --recruit_evalue_cutoff 0.001 \
    --recruit_hit_length_percentage 0.9

--clusters_chunk is a headerless TSV of representative<TAB>member. --fasta_file must be an uncompressed FASTA — Easel cannot seek within a gzip stream — and its sequence names must be unique.

Sequence names

A record that is a slice of a larger protein may say so in either of two spellings, and both are read identically:

The second is the form this tool emits, so <chunk>_reps.fasta from one release can be used directly as the database for the next without its coordinates being lost. The base keeps any slashes it carries: 3387826881/v1/356-472 is region 356–472 of the protein 3387826881/v1.

Bounds are read as coordinates only if they span the record exactly. scaffold_12_34 holding 15 residues is a whole protein named scaffold_12_34, not residues 12–34 of scaffold. Anything else is identity and is kept whole — 3387826881/356_472, 3387826881/356, 3387826881/v1 and 3387826881/356-472-243 are four distinct protein names, none of them carrying a region.

Any other character is allowed in a name, including further slashes. Names are never split on their first slash, so two records sharing a prefix stay distinct. X and a literal X/1_10 remain distinct even when the same family recruits residues 1–10 of X alongside the complete X/1_10 record. Literal percent sequences such as %2F are preserved too, in both update modes and in all emitted identities.

No name is reserved, but the slice spelling is not inert either. Whether a record is independent of 3387826881 depends on which spelling it uses and on its own length:

The middle row is the round-trip working as intended: a record that says it is a region of 3387826881 is reported at those parent coordinates, exactly as the corresponding residues of 3387826881 itself would be. If a database contains both, the same residues are the same protein region and get the same name — they are not two things. Include the parent and its own slices in one database only if that is what you mean.

Optional flags

Pass every threshold explicitly on a production run. The defaults exist for ad-hoc use; relying on them means a forgotten flag produces a plausible-looking family set instead of an error.

Streaming re-reads and re-parses the database once per query. --prefetch_targets parses it once and keeps it in RAM; the results are byte-identical either way, so the flag is purely a memory-vs-time dial. Leave it off unless the database fits comfortably in RAM.

On a production run, build the index once and share it. Every chunk task would otherwise re-index the whole database:

# once, upstream -- either of these
uv run python -c "from mgnifam.generate_families import build_ssi_index; \
                  build_ssi_index('db.fa', 'db.fa.ssi')"
esl-sfetch --index db.fa                     # HMMER/Easel, e.g. the nf-core module
# then, per chunk
uv run mgnifam generate_families --fasta_index db.fa.ssi ...

A supplied index is used as given and never rebuilt, so parallel chunk tasks can share one read-only index safely — including one staged as a symlink by a workflow manager. It is an error for it to be missing rather than a request to build one there, and one that does not match the FASTA fails at the first fetch instead of being silently replaced. The FASTA's filename need not match the one it was indexed under.

An index from esl-sfetch --index is interchangeable with one from build_ssi_index for whole-record fetches, which is all generate_families performs. The two are not byte-identical: esl-sfetch also records each record's data_offset and record_length, which enables esl-sfetch -c <from>..<to> subsequence fetches against its own index but not against ours, and it sizes the index's filename field from the path you typed, so its output is not reproducible across directories. Ours is.

Updating existing families

mgnifam update_families refreshes families that already exist as HMMs against a new database. It is the answer to "a new release came out" — you do not re-derive the families from their original clusters, you search the models you already have.

uv run mgnifam update_families \
    --hmm_input previous_output/hmm \
    --fasta_file new_release.fa

--hmm_input is either a directory of .hmm/.hmm.gz files or a single multi-model library (hmm.lib.gz works). Both forms produce identical output for the same models. --fasta_file is uncompressed, for the same Easel reason as above.

Every threshold flag from generate_families carries over with the same name and default.

What identity means here

A family keeps the name its model carries in its NAME field — 1_7 stays 1_7 across releases, in the filenames and in every identity-bearing field inside the outputs. Two consequences:

• <chunk>_updated_metadata.csv's family_id column holds 1_7, not a bare integer. That differs from generate_families, whose ids are a rank.
• Chunks sharing one output root must own disjoint family names. Nothing enforces it, because the names come from the input models rather than from --chunk_id.

A NAME must match [A-Za-z0-9._-]+ and be neither . nor ... It is interpolated into artifact paths and into CSV fields, and it arrives from a file this tool did not write.

Outputs

Per-family artifacts land in the same hmm/, full_msa/, seed_msa/ and rf/ directories, named by family. Aggregates are <chunk>_updated_*: families.tsv, metadata.csv, discarded.csv, successful.txt, converged.txt, reps.fasta.gz, delta.csv, and <chunk>_updated.log.

<chunk>_updated_delta.csv is what an update run is for — one row per family, whether it survived or not:

family_id,model_length_before,model_length_after,round1_recruits,full_msa_size,retention,rounds_run,converged,outcome

Every field but family_id, model_length_before and outcome may be empty, because a family discarded early never reached the stage that would produce one. model_length_after is the length of the model that recruited the final membership. retention is the fraction of round 1's own recruits still present at the end — under --skip_refine that is 1.0 by construction, since there are no later rounds to drift.

outcome is successful or the discard reason. no hits in the new database means the model found nothing at all in the new release; low complexity model - confounding cluster means it found hits and none cleared the envelope-length filter. For an update run that distinction is the point.

Give each run its own --output_dir. Re-running the same models into the same directory is allowed, so a failed chunk can be retried in place. Running a smaller set of models over a directory that still holds a larger one is refused rather than silently cleaned up: generate_families can clear its own past output because it derives names as <chunk>_<rank>, but an updated family keeps its model's name and --chunk_id never appears in a per-family filename, so nothing on disk says which run wrote hmm/1_7.hmm.gz.

Input files must not overlap output paths, including through symlinks or hard links. This is checked before writing, for both model directories and single-file libraries. On an accepted retry, previous artifacts for the input family names are removed before processing. Discarded families therefore leave no old models, and --skip_refine leaves no seed/RF files from a previous refine run. A cleanup failure aborts the run.

Cost

hmmsearch is O(n_families x database) and this command does not change that. --skip_refine is one database pass per family; refining is up to three. Against a billion-sequence release that term, not the alignment, is what to budget.

mgnifam --help lists the subcommands, and python -m mgnifam is equivalent to the console script.

Outputs

Written under --output_dir (default: output), keyed by --chunk_id:

One file per family, so one directory each:

One file per chunk, so flat in the output root:

Family ids are a 1-based rank among successful families, in cluster-file order.

Both CSVs carry a header row, so they load with pandas.read_csv as they are:

protein is quoted, with embedded quotes doubled; a protein or representative containing a comma or a quote is escaped, so both files parse with a standard CSV reader. Literal slashes stay in protein, including punctuation after a slash: protein/v1,variant is one protein field. Only a trailing coordinate range spanning the emitted sequence is separated into region. region is <start>-<end> on the parent protein, or - when the representative spans a whole unsliced record. Those two columns together are the <base>/<start>-<end> spelling above, which is also how <chunk>_reps.fasta names its records. The representative is the highest-scoring reported domain of HMMER's top-ranked hit. The header is written before the run starts, so a chunk that produces no families still yields a parseable file.

Exit status

The status describes whether the output is safe to consume, not only whether the process stopped:

Why this is fast now

The previous implementation took roughly eight months to process the full database. Three defects accounted for most of it:

1. run_initial_msa was O(database × members), per family. A map() iterator was rebuilt inside a comprehension's condition, turning a membership test into a full linear scan of the cluster for every one of the billions of database sequences. It is now a constant-time SSI lookup per member.
2. Cluster selection was O(N²) — the cluster table was boolean-masked and re-filtered once per family. It is now a single grouping pass.
3. The exit-branch hmmsearch re-ran a search that had just been performed with the identical HMM, differing only in a post-filter. Its hits are now cached and re-filtered, saving a full database pass per family.

On top of that, the entire FASTA was held in RAM twice — once as a DigitalSequenceBlock and once as a Python dict of DigitalSequence objects. In the default (streaming) mode both are gone: targets stream from disk, and random access goes through an Easel SSI index. Passing --prefetch_targets deliberately restores the first copy, trading that memory back for speed. Families are searched in batched waves, so hmmsearch uses up to --cpus workers whenever enough families remain in the wave.

Reproducibility

For the dependency set resolved in the committed uv.lock (install with uv sync --frozen), scientific outputs are byte-identical across repeated runs, across PYTHONHASHSEED values, across --batch_size, across --prefetch_targets, and — unlike the previous implementation — across --cpus. The contract is scoped to that lockfile: pyhmmer, pyfamsa and pytrimal decide hit retention, alignment and serialised bytes. (<chunk>.log carries timestamps and is excluded from that contract. HMM files omit the DATE and COM lines, which are otherwise a wall-clock and an argv dump.)

The old pipeline's recruitment depended on how many CPUs it was given. pyhmmer selects parallel="targets" whenever the query count is below the CPU count, which was every call in the old family-at-a-time loop. Each worker runs its own Pipeline over a slice of the database, and the merge concatenates each slice's stored hits while re-thresholding only the reporting flags. Z and domZ come out identical, but the stored list grows — and the old code iterated that raw list rather than .reported. Measured on the 50 000-sequence fixture under the old pinned pyhmmer==0.11.1: len(TopHits) goes 26/19/55 at --cpus 1 to 27/19/56 at --cpus 4, while .reported stays 26/19/54 throughout.

Forcing parallel="queries" fixes this: the answer is the same at any core count. Those extra stored hits were exactly the ones failing the reporting threshold, so reading .reported closes both halves of the problem at once.

Two bugs fixed, and what they change

Recruitment ignored --recruit_evalue_cutoff. The old code iterated the raw TopHits, which retains hits pyhmmer stored but did not report. Extraction now reads top_hits.reported. On the small fixture, family 4497037939_1_144 used to recruit sequence 6320430079, which is stored but below the reporting threshold. Families are correspondingly smaller: on that fixture, 32/19/65 members become 31/19/61. Same families, same representatives, fewer spurious members.

Recruitment depended on the CPU count, as described above. Both halves are fixed, so --recruit_evalue_cutoff now means what it says, on any machine.

Outputs are therefore not byte-compatible with the legacy script. Every difference is enumerated in CHANGELOG.md.

Indexing a very large database

Easel buffers up to 2 GB of keys in RAM before spilling to an external sort, which then needs scratch space in TMPDIR plus room for the final index (roughly n_sequences x (name_length + 16) bytes). Size TMPDIR accordingly before indexing a billion-record FASTA.

Development

uv lock --check && uv sync --frozen
uv run pre-commit install
uv run pre-commit run --all-files
uv run pytest