faba

Feature statistics Accumulator for Base-pair-level Analysis

faba extracts per-cell genomic features directly from alignment (BAM) files: RNA modifications (DART-seq m6A, A-to-I editing), alternative polyadenylation (APA), gene counts, read depth, and SNP genotypes.

A quick look

Run in a terminal, run, qc, pileup, metagene, pwm and cluster-view open full-screen views; --batch-process skips them for scripts and clusters.

Each view runs as a detachable session: Ctrl+Alt+D leaves it running (as does a dropped connection), and starting the command again with nothing named, or faba attach, returns to it; Ctrl+Alt+S switches to another session. faba sessions lists them, faba kill stops one as its quit key does (a write under way finishes first) and faba rename labels one. --no-detach runs a view in the terminal only.

faba run: pick the BAMs (signal and control), steps and flags, preview the command, and run it.

faba qc: set site thresholds, see what each one keeps, and write a filtered copy.

faba pileup: browse a gene's modification sites, one track per batch.

faba pwm: the sequence logo around the called sites.

Installation

Prerequisites

These must be installed before running cargo install — Cargo only builds Rust crates and cannot pull in a C toolchain or system libraries on its own. If they are missing, the build fails midway in the vendored-htslib build script.

• Rust (stable, edition 2021+) — install via rustup.

• A C/C++ toolchain and the system libraries that rust-htslib builds against (it vendors and compiles htslib, which needs the compression dev headers, and runs bindgen, which needs libclang):

  # Debian / Ubuntu
  sudo apt-get install build-essential clang libclang-dev \
      zlib1g-dev libbz2-dev liblzma-dev pkg-config
  
  # macOS (Homebrew)
  brew install llvm xz bzip2 zlib

Install from crates.io

cargo install faba

Optional HDF5 (.h5 / .h5ad) support:

cargo install faba --features hdf5

Install from GitHub

cargo install --git https://github.com/causalpathlab/faba.git

Pin to a specific tag or branch if you need a reproducible build:

cargo install --git https://github.com/causalpathlab/faba.git --tag v0.14.0
cargo install --git https://github.com/causalpathlab/faba.git --branch main

Build from a local clone

git clone https://github.com/causalpathlab/faba.git
cd faba
cargo build --release
cargo install --path .

Optional features

Feature	Enables	Extra requirement
mcp (default)	faba mcp, the commands as tools for an agent	none
hdf5	HDF5 backend + .h5/.h5ad readers	libhdf5 on the build host

There is no cuda / metal feature: nothing on the BAM-to-matrix path touches a GPU.

Add them with --features:

cargo install --git https://github.com/causalpathlab/legume-rs.git faba --features hdf5

Verify

faba --help

Troubleshooting

Almost all install failures come from the hts-sys build step (it compiles a vendored copy of htslib and generates bindings with bindgen). The fixes:

• thread 'main' panicked ... Unable to find libclang — bindgen can't locate libclang. Install libclang-dev (Debian/Ubuntu) or brew install llvm (macOS), then point to it:

  # Linux: the unversioned symlink ships with libclang-dev
  export LIBCLANG_PATH=$(llvm-config --libdir 2>/dev/null || echo /usr/lib/llvm-18/lib)
  
  # macOS (Homebrew llvm is keg-only)
  export LIBCLANG_PATH="$(brew --prefix llvm)/lib"

• fatal error: zlib.h / bzlib.h / lzma.h: No such file or directory — the compression dev headers are missing (the runtime .so alone is not enough). Install zlib1g-dev libbz2-dev liblzma-dev (Debian/Ubuntu) or brew install xz bzip2 zlib (macOS).

• error: linker 'cc' not found — no C toolchain. Install build-essential (Debian/Ubuntu) or Xcode Command Line Tools (xcode-select --install) on macOS.

• failed to run custom build command for 'hdf5-sys' — you passed --features hdf5 without libhdf5 on the host. Install libhdf5-dev (Debian/Ubuntu) or brew install hdf5, or drop the feature.

If a build fails after fixing a prerequisite, re-run with a clean rebuild of the C bits: cargo install --git ... faba --force.

Usage

faba <COMMAND> [OPTIONS]

Command	Purpose
Feature profiling — BAM → per-cell features
dartseq (dart, m6a)	Call DART-seq m6A sites by a WT-vs-MUT control contrast on C-to-T conversions
atoi (a2i, editing)	Detect and quantify A-to-I RNA editing sites
apa (polya)	Quantify alternative polyadenylation sites per cell
count (genes)	Count reads per gene and call cells (single-cell or bulk RNA-seq)
depth (rd)	Compute read depth over genomic intervals
snp (genotype)	Discover and genotype SNP variants from BAM pileup
run	Set up the pipeline in a view (inputs, steps, flags), preview and save the command, and run it: SNP → count → ATOI → m6A → APA
QC — choose and apply thresholds after profiling
qc	Pick thresholds in a full-screen view, then filter a faba output directory into a new fileset
Clustering: cells by their own sites
cluster	Cluster cells by their site reads: cell and gene modules with a methylation shift per pair, fitted by variational inference. --modality count fits cell modules by gene counts instead; --within its run record then fits the sites within each of those cell modules
cluster-view (modules, blocks)	View a cluster fit from its run record {m}_blockmodel.run.json, or from its input directory (it never fits: run cluster first): its gene modules against its cell modules, each module's sites, the genes you pin side by side, and a gene's full pileup, one track per cell module; --html writes it as one self-contained page to share. A --within record ({m}_within.run.json) opens as the list of its fits, the count fit first; --group C3 opens the one within a cell module
Inspection & reference
pwm	Build a position weight matrix around genomic sites
pileup (inspect)	ASCII pileup, or a faceted Miami plot, for one gene
metagene (mg)	Metagene histogram of site positions across gene features
docs	Print the method write-ups compiled into this binary
mcp	Serve the commands to an agent over the Model Context Protocol
attach, sessions, kill, rename	Return to, list, stop and label the views' detached sessions

Run faba <COMMAND> --help for the detailed options of each subcommand.

For an agent (MCP)

faba mcp serves the commands as Model Context Protocol tools over stdio, so any MCP client can drive them, with any model behind it (Claude, GPT, Gemini, or a local one). faba sends nothing anywhere itself: the model sees what the tools return, and where that goes is up to the client and model you choose.

Connect a client

The server is one command, faba mcp. A client also needs to know where the project is: the server keeps its records (below) in its working directory, or in LEGUME_MCP_PROJECT when set. Set it whenever the client may start servers elsewhere, as desktop apps usually do. If the client cannot find faba, give its full path (which faba).

Most clients (Claude Desktop, Cursor's .cursor/mcp.json, Gemini CLI's settings.json, and others) take this shape:

{
  "mcpServers": {
    "faba": {
      "command": "faba",
      "args": ["mcp"],
      "env": { "LEGUME_MCP_PROJECT": "/path/to/project" }
    }
  }
}

Claude Code: claude mcp add faba -e LEGUME_MCP_PROJECT="$PWD" -- faba mcp. VS Code takes "servers" with "type": "stdio" in .vscode/mcp.json; Codex CLI takes [mcp_servers.faba] with command, args and env in ~/.codex/config.toml. Key names differ between clients; the server needs only the command, its argument and, optionally, the environment.

What the agent does

The server tells the agent how to work when it connects (the MCP instructions). Some clients do not pass those to the model; with one of them, put this in the agent's own instructions (an AGENTS.md, a system prompt):

faba turns coordinate-sorted BAM alignments into per-cell feature matrices: gene counts, read depth, A-to-I editing, DART-seq m6A, alternative polyadenylation and SNP genotypes. Row names are {gene_key}/{modality}/{detail}.

Producers write every candidate and apply no p-value cutoff. qc is where cells, features and sites are cut, and it writes a new directory.

1. docs lists the method write-ups. With `query` it prints the sections named (a topic such as profiling, a command such as qc, or a section number).
2. count, atoi, dartseq, apa, depth and snp each write one modality. run writes the pipeline (SNP, then count, A-to-I, m6A, APA) into one directory. dartseq needs control BAMs; without them, run skips m6A.
3. qc filters a directory into a new one, with the site thresholds given. It does not edit the input. With `dry_run` it writes nothing and returns what the cell QC keeps and what each threshold removes: try thresholds that way first.
4. pileup, metagene and pwm read a directory's sites. metagene and pwm return their table and its figure as an image. pileup returns a binned table for a gene or region (and, with `figure`, a Miami SVG, as an image too). cluster_view shows a cluster fit level by level, from its run record ({m}_blockmodel.run.json) or the directory holding its tables: the modules, then a gene module, a gene, the cell modules, the trace and verdict, the settings.
5. cluster fits cells by their own site reads (a blockmodel of cell and gene modules, by variational inference) and writes its tables to `out`. In two rounds: first `modality` count, cell modules by gene counts (it writes count_blockmodel.run.json); then, with `within` naming that record, m6a or atoi within each of its cell modules. The second round writes {m}_within.run.json, which cluster_view opens: the list of its fits, and with `group` the fit within one cell module.

The producers (count, depth, atoi, dartseq, apa, snp, run, qc, cluster) take `background`: on large BAMs set it, and the call returns a job id at once.

Every write goes to `out`, which must not exist. Leave `out` out and the output goes in today's results folder, named after the input and the tool. Paths resolve against the server's working directory. Common flags have their own arguments; further flags go in `options`, as `faba <command> --help` lists them; `options` cannot name the output or open a view. Results give the paths written and what the run recorded (`results`: a fit's ELBO verdict, its cluster sizes).

legume-mcp adds where the records are, the notebook tool (open with a plan: the question, the data, how the result will be checked) and the background jobs below.

Tool	Does
docs	the method write-ups
count, atoi, dartseq, apa, depth, snp	one modality each, into a new directory
run	the pipeline (SNP, then count, A-to-I, m6A, APA) into one directory
qc	filters a directory into a new one
cluster	fits cells by their own site reads, by gene counts (modality count), or by site reads within a count fit's cell modules (within); writes its tables to a new directory
pileup, metagene, pwm, cluster_view	read what was written
notebook, job_status, job_wait, job_cancel	from legume-mcp

Long runs in the background

On large BAMs a producer can outlast a client's wait on one call. Give it background: true and it returns at once with {job, written}; the agent follows the job with job_wait or job_status (or stops it with job_cancel), and uses written only once it is done. A job that succeeds is recorded like any step. Jobs end with their server: a run longer than a client session is for faba on the command line.

Records and watching

Every write goes to an out that must not exist yet: no tool replaces a directory or a file, and no tool opens a view. The server keeps the project's records the way legume-mcp lays them out, after Noble (2009), A quick guide to organizing computational biology projects, PLoS Comput Biol 5(7):e1000424: every step that writes goes into the day's results/<date>/runall.sh, which replays the session, and into results/INDEX.tsv; outputs the agent does not name go in results/<date>/. LEGUME_MCP_RECORDS=off turns the records off. While an agent works, the server also writes a live feed a person can follow and steer (stop a call or a job, leave the agent a note) with legume-mcp's watch; LEGUME_MCP_WATCH=off writes none.

Clustering on real data

faba cluster (cell and gene modules of the site reads) fits by variational inference, as the count fit does: each restart (--chains, 4) starts from k-means++ on the cells' residuals over their genes, runs --svi-epochs minibatch epochs, then coordinate ascent until an epoch moves the ELBO by less than --tol per read, and the best restart by ELBO is kept. A read's likelihood is bounded by Pólya-Gamma (Jaakkola and Jordan), one bound per read and block, so each link (a cell's reads at one gene) weighs every block at a few sums; the shifts and the cells' levels trade a constant through their priors after every global step. On the paper data (11,919 cells, 2.35M links) a restart takes about two minutes. The MCMC fits faba had before 0.19 (--model link and crp, the sampler's --sweeps and --burn-in, --holdout) are gone, and their flags are refused by name; their records still open in cluster-view.

# Cell and gene modules from a faba output directory, watching the restarts as they run
faba cluster out/
# Open the fit, as it is
faba cluster-view out/m6a_blockmodel.run.json

In two rounds. Counts and site reads are not fitted jointly: a cell has far more UMIs than site reads, so the counts would decide the modules alone. Instead the cells are first clustered by their gene counts, and the sites are then fitted within each of those cell modules:

# Round 1: cell modules by gene counts (writes out/count_blockmodel.run.json)
faba cluster out/ --modality count
# Round 2: m6A (or A-to-I: --modality atoi) within each count cell module
faba cluster --within out/count_blockmodel.run.json
# The count fit and the fits within it; --group C3 opens one of them
faba cluster-view out/m6a_within.run.json

The count fit is by variational inference (after Gopalan, Hofman and Blei's hierarchical Poisson factorization): each restart (--chains, 4) starts from k-means++ on the cells' normalized profiles, then the genes' profiles over those clusters, then runs online minibatch epochs (--svi-epochs 10 of --batch-size 1000 cells), then coordinate ascent until an epoch moves the ELBO by less than --tol (1e-5) per count, at most --epochs (200) in all. Only nonzero counts are visited, and nothing cells × genes is held dense. The tables are the best restart's by ELBO; the trace is each restart's ELBO by epoch, and the verdict says whether the best one settled.

The second round reads the site matrices beside the count fit's inputs, unless some are named, and matches cells by batch and barcode. Each cell goes to its main count module (C1, C2, …). A group with fewer than --within-min-cells (50) cells to fit is skipped. Each group's fit is an ordinary one under within/{m}/{group}/, with its own run record and verdict; {m}_within.run.json lists them, and {m}_within_cells.parquet gives each cell its group and its site module there. cluster-view --html fit.html writes the count fit's page and fit.{group}.html beside it for each fitted group, linked to each other.

Batch effects, after the fit. Batches fitted together can come out as modules of their own. cluster-view reads each cell module's spread over the batches off the cell table, each batch's cells weighing alike, with the batches as replicates of conditions: a batch's name less its replicate token (rep1_mut and mut_r2 are replicates of mut; a folder, runA/mut, is one too), or as --batch-groups rep1_mut=mut,... says. A module even over each condition's replicates (normalized entropy at least 0.8, each condition weighing alike however many replicates it has) is shared, or that condition's when 80% of it is in one condition; uneven, it is batch-specific, as a condition's replicates should agree. A module mostly where no condition has two replicates, or with under 10 cells' worth, is not judged. In the browser, b opens the modules against the batches and B has the panes show the modules that are not batch-specific, then those that are, then all; --batches prints the same.

faba cluster-view out/count_blockmodel.run.json --batches

GO terms of the gene modules. cluster-view --go tests each gene module's genes (those whose main gene module it is, at a weight of at least --go-min-weight, 0.5) against the fit's own genes that GO annotates: the hypergeometric tail of the overlap, Haldane's odds ratio, and BH over every (module, term) pair, for terms of 5 to 500 annotated genes, propagated up the GO graph (is_a, part_of) with NAS and ND evidence left out. A module's terms are listed by p, a parent left out when it holds the same genes. The annotations default to lupin's cached GO Consortium file for the genes' species (human or mouse); --gaf and --obo name others, --gmt any gene sets. In the browser, a opens the gene module in focus's terms (↑/↓ steps the module); MCP's cluster_view takes go.

faba cluster-view out/m6a_blockmodel.run.json --go

Watching the restarts. On a terminal faba cluster draws each restart's ELBO live while it fits (q cancels, Ctrl+C stops the restarts early and writes what they reached). When it is done it prints the same traces as one sparkline per restart on one scale, then a verdict and one line of advice naming the flags to run again with. Once the tables are open, the browser's header carries an ELBO ... ok or check badge, and t opens the trace view; --html pages carry the trace chart too. Each restart's trace is kept in {m}_trace.parquet.

Judging convergence. A fit has converged when the best restart's ELBO settled within --tol per unit (count or read) before --epochs ran out; the verdict also gives the spread of the restarts' ELBO per unit, how far apart their optima are. {m}_blockmodel.run.json keeps each restart's ELBO and epochs, the best, and whether it settled, under results.fit. When it has not settled, run again with more --epochs; restarts far apart point to more roles than the data support: try fewer --cell-roles or --gene-roles, or more --chains.

Figures for a paper

Every figure faba writes takes its size on the page in inches: --fig-width and --fig-height, either or both (the other is then the figure's own). It is laid out for that size, its text set in points as it prints, rather than scaled. This holds for cluster-view --figure (where 0.19.0's --width and --height still work), metagene --figure, pwm --figure, qc --plots, and pileup's Miami plot (8 inches wide unless given; its tracks share the height given). Over MCP, each of these tools takes width and height. A size too small for a figure's labels is refused, saying the least it needs.

# A column-wide modules grid and a page-wide logo
faba cluster-view out/m6a_blockmodel.run.json --figure grid --fig-width 3.5 --fig-height 3
faba pwm out/ -f genome.fa --batch-process --figure logo --fig-width 7

Methods

The method write-ups are compiled into the binary — include_str!, not files read at runtime — so they travel with a cargo installed faba or one copied to a cluster with no checkout beside it, and the build fails if one of them goes missing.

faba docs                # list what there is
faba docs profiling      # BAM -> per-cell features: m6A, A-to-I, APA, counts, SNPs

The same files live in docs/, which also carries the design notes for work that is planned but not implemented (kept separate on purpose — reading a plan as though it described the code is how people end up debugging things that were never built).

Examples

# Gene counts and cell calling from a single-cell BAM
faba count sample.bam -g genes.gff -o out/

# A-to-I editing sites (every putative site is written; `faba qc` decides)
faba atoi sample.bam -g genes.gff -f genome.fa -o out/

# DART-seq m6A: signal (WT APOBEC1-YTH) vs catalytically-dead control (YTHmut).
# A control is REQUIRED — m6A can't be told apart from genomic C/T variation
# without it (--mut / --control / --background are accepted aliases).
faba dartseq wt.bam --control-bam ctrl.bam -g genes.gff -f genome.fa -o out/

# Everything in one pass (the m6A step runs only when --control-bam is given,
# otherwise it is skipped; the other steps need no control).
# Set up and run everything in a view (plain `faba run`), or straight through:
faba run
faba run --batch-process sample.bam -g genes.gff -f genome.fa -o out/ --control-bam ctrl.bam

# The producers apply no p-value / effect-size / reproducibility cutoff.
# Pick the thresholds in a full-screen view, then cut into a NEW directory
# (plain `faba qc` asks for both directories in a pop-up):
faba qc out/ -o out_qc/
# ...or, in a script, cut with the thresholds as given:
faba qc out/ -o out_qc/ --batch-process --site-max-pv 0.05 --site-min-cells 10 --auto-cutoff

License

MIT — see the workspace license.