CAT2 (Comparative Annotation Toolkit, v2) is an improved version of CAT, a pipeline that annotates one or more target genomes from a Cactus HAL alignment plus a reference annotation. It combines several complementary gene-finding strategies and reconciles them into a single consensus gene set per genome:
- transMap — projects the reference annotation across the HAL alignment (recovers orthologs of reference genes).
- AugustusTM/TMR — refines transMap projections with Augustus using RNA-seq and annotation hints.
- AugustusPB — Augustus with PacBio/IsoSeq hints for isoform discovery.
- augMP (miniprot) — aligns a protein database to each genome; See Protein set.
- StringTie — reference-guided transcript assembly from RNA-seq or Iso-seq.
This is the current version of CAT2 and is under active development.
- Linux x86-64,
conda/mamba(miniforge recommended) - A Cactus install (for
the HAL tools:
halStats,hal2fasta,halLiftover, …) - A cluster (SLURM or SGE) is recommended for real genomes; small inputs run locally.
git clone https://github.com/ph09/CAT2_smk.git
cd CAT2_smk
conda env create -f environment.yaml
conda activate cat2
pip install --no-deps -e . # installs the cat2 package itself
chmod +x install_standalones.sh
./install_standalones.sh # fetches UCSC Kent binaries into ./standalones/Point CACTUS_BIN at a Cactus binary install that provides HAL tools
(halStats, hal2fasta, …), then source the env helper. Conda stays first on
PATH so a Cactus-bundled Python cannot shadow the cat2 env:
export CACTUS_BIN=/path/to/cactus-bin-v3.2.1
source setup_env.sh
which python # should be .../envs/cat2/bin/python
which halStats # should be under $CACTUS_BINOn sites that ship Cactus as an environment module, load the module (or set
CACTUS_BIN to that prefix), then conda activate cat2 again so conda Python
wins. Set site scheduler knobs in your config (slurm.partition,
slurm.exclude_nodes, or cluster.sge.*); empty partition/exclude means “use
the cluster default”.
The repo ships small parser/test fixtures, but the test HAL and BAMs are hosted
separately. Download them into test_data/ first:
- HAL: https://public.gi.ucsc.edu/~pnhebbar/share/testData/vertebrates.hal
- BAMs: https://public.gi.ucsc.edu/~pnhebbar/share/testData/bams/
Then dry-run, then run:
# dry-run: parse Snakefile, validate config, build the DAG
snakemake --configfile input.yaml --cores 4 -n all
# run locally
snakemake --configfile input.yaml --cores 4 allinput.yaml is the template config you will copy and modify for your own data.
If you are looking to also include liftoff in the CAT2 consensus (recommended only when annotating same species or closely related species of the reference), use the optional-liftoff branch.
input.yaml is heavily commented; this section highlights the parts you will
most often change. Every key not shown here has a sensible built-in default.
| Key | Meaning |
|---|---|
work_dir |
Output directory for this run (created if missing). |
hal |
Path to the Cactus HAL alignment. |
annotation |
Reference annotation (GFF3) for ref_genome. |
ref_genome |
Name of the reference genome as it appears in the HAL. |
genomes |
List of target genome names (HAL names) to annotate. |
work_dir: "my_run"
hal: "/data/primates.hal"
annotation: "/data/GRCh38.gff3"
ref_genome: "hg38"
genomes: ["hg38", "rheMac10", "calJac4"]augMP aligns a protein database to every target genome. With a single-species (e.g. human-only) protein set it can only re-find genes that already exist in the reference. To discover lineage-specific / non-reference genes you must give it proteins from those lineages. There are two ways to provide proteins.
(A) A protein FASTA you already have:
protein_fasta: "/data/my_proteins.fa"(B) Auto-built multi-species DB — list species, NCBI taxon IDs, and/or clades,
and CAT2 downloads UniProt reference proteomes, merges and de-duplicates them, and
uses the result as the miniprot input. This runs once on the controller machine
(needs outbound internet) and is cached. If protein_db is set it takes
precedence over protein_fasta.
protein_db:
species: # names OR HAL genome names (auto-normalised)
- "Homo sapiens"
- "Macaca mulatta"
- "Callithrix jacchus"
- "Mus musculus" # an outgroup broadens coverage
# taxa: [9606, 9544] # alternative/supplement: NCBI taxon IDs
# clades: pull proteomes for members of a genus/family (or any higher rank).
# Use this to cover poorly annotated species via their better-annotated
# relatives. CAT2 finds clade members that have a proteome and keeps the best
# (highest quality first, up to max_per_clade).
clades:
- "Cercopithecidae" # Old World monkeys (family)
- "Callitrichidae" # marmosets/tamarins (family)
max_per_clade: 25 # cap proteomes per clade (best quality first)
clade_include_other: false # also accept non-reference proteomes when expanding clades
base_fasta: "" # optionally fold in your own proteins too
dedup: true # collapse identical sequences (keeps miniprot fast)
min_len: 20 # drop proteins shorter than this
out: "my_run/protein_db/protein_db.fa" # built DB path (this is the default location)Tip: pick a handful of well-spread representatives rather than every species —
most genes are shared, so dumping hundreds of proteomes is slower with little
recall gain. Genus-/family-level clades are the best way to get coverage for
species that lack their own good proteome.
Building the DB by hand (same machinery, standalone) and pointing
protein_fasta at the result:
python scripts/build_protein_db.py \
--out my_run/protein_db.fa \
--clades "Cercopithecidae,Callitrichidae" \
--species "Homo sapiens,Mus musculus" \
--base-fasta /data/GRCh38.prot.faUseful build_protein_db.py flags: --species, --taxa, --clades,
--max-per-clade (default 25), --clade-include-other, --base-fasta,
--min-len (default 20), --no-dedup, --cache-dir, --strict (fail on any
species that will not resolve), --summary (per-species TSV). Species names are
normalised leniently, so HAL genome names paste in directly
(PR00246~Eulemur_fulvus.pri → Eulemur fulvus).
miniprot sensitivity is tunable under a miniprot: block (all keys optional;
defaults are already tuned more permissively than miniprot's own to recover more
divergent/paralogous copies). Raising sensitivity finds more candidates at the
cost of runtime and false positives (which consensus filtering then prunes). See
the miniprot: block in input.yaml for every knob (splice_model,
max_intron, min_secondary_ratio, max_secondary, out_n/out_s/out_c, …).
Per-genome RNA-seq / long-read BAMs, plus which genomes use them:
transcriptomic_data:
hg38:
bam: # short-read RNA-seq
- "/data/hg38/rnaseq1.bam"
isoseq_bam: # long reads (PacBio/IsoSeq/ONT)
- "/data/hg38/isoseq.bam"
intronbam: # noisy short reads, used for intron hints only
- "/data/hg38/noisy.bam"
rnaseq_genomes: ["hg38"] # genomes with usable short-read RNA-seq
isoseq_genomes: ["hg38"] # genomes with usable long readsbam— short-read RNA-seq (splice + coverage hints)isoseq_bam— long reads (isoform evidence, drives AugustusPB / StringTie)intronbam— noisy short reads used only for intron hints
augustus: true # AugustusTM/TMR refinement
augustus_pb: true # AugustusPB (needs isoseq)
stringtie: true # StringTie assembly
stringtie_genomes: ["hg38"]
txTM: true # transMap at the transcript-level. Only helpful in genomes close to reference annotations set.
augustus_species: "human" # Augustus species parameter sethigh_recall is an opt-in master switch. When true, the recall-limiting gates
across every mode (transMap paralog/overlap filtering, coverage floors,
consensus length/CNV/fragment cutoffs, denovo support, postprocess low-support
drop) are loosened together so fewer genes are missed, at the cost of more false
positives. It overrides the individual knobs below.
high_recall: false # leave false and tune individual keys for fine controlFine-grained knobs (all documented inline in input.yaml) include transMap
filtering (global_near_best, tm_filter_overlapping, tm_min_cover, …),
consensus fragment reclassification (consensus_fragment_max_coverage/identity),
and postprocess drop tuning (postprocess_*).
Annotate reconstructed ancestral (internal HAL tree) genomes using alignment-only modes:
annotate_ancestors: true
ancestor_genomes: ["Anc0", "Anc1"] # optional; auto-detected if omitted
ancestor_modes: ["transMap", "transMap_pairwise", "txTM"]execution_mode: "auto" # auto | slurm | sge | localauto— detect the scheduler (sbatch→ slurm;qsub+$SGE_ROOT→ sge; otherwise local).slurm/sge— submit jobs to the cluster; per-rulemem/cpus/timecome from theslurm:→rules:block (shared by both schedulers). SGE-specific knobs (queue, parallel environment, memory flag) live undercluster.sge.local— run everything on the current machine; per-rule hints come from thelocal:block.
The slurm: block also has site-specific fields to review before running:
partition and exclude_nodes (the shipped value is UCSC-specific — change or
clear it).
source setup_env.sh # every new shell
# recommended: use the launcher (enforces snakemake >= 9, resumable)
./run_pipeline.sh --work-dir my_run --configfile input.yaml --cores 32
# long runs: keep alive under tmux/nohup
tmux new -s cat2 './run_pipeline.sh --work-dir my_run 2>&1 | tee my_run/run.log'
# or drive snakemake directly
snakemake --configfile input.yaml all --cores 32 --keep-going --rerun-incompletesnakemake must be ≥ 9. Snakemake 8.x corrupts multi-line f-strings inside rule
run:blocks under Python 3.12; the launcher refuses to start otherwise.
Final per-genome consensus annotations are written under work_dir (GFF3 / GenePred).
A fast, cluster-free smoke test (used by CI) checks the snakemake version, a
run:-block f-string regression guard, and that the cat2 package imports:
./scripts/smoke_test.shIt skips the full DAG dry-run automatically when halStats / the test HAL are
not present (e.g. in CI).