Converts a multiple alignment format (MAF) file into the tabix-indexable BED
that JBrowse 2's MafTabixAdapter reads. It streams: MAF on stdin, BED on
stdout, constant memory, so a whole-genome alignment costs no more RAM than a
small one.
MAF support is built into JBrowse 2 core now (a MafTrack with a
LinearMafDisplay) — you no longer need to install a plugin. See the
MAF track config guide
for the track JSON, and the
MAF track user guide for
what the display can do with it. The older standalone
jbrowse-plugin-mafviewer
still exists, and also reads this format.
Install
cargo install maf2bed
maf2bed --help
Usage
The first argument is the assembly name to use as the BED's reference — the
part of the MAF src token before the contig. Rows whose src starts with
<asm>. supply the BED's first three columns, so maf2bed hg38 turns a
hg38.chr1 row into a chr1 line.
export LC_ALL=C # sort and tabix have to agree on collation; also faster maf2bed hg38 < file.maf | sort -k1,1 -k2,2n | bgzip > file.bed.gz tabix -p bed file.bed.gz
It only reads stdin, which is the point — pipe a compressed MAF straight in and never write the plaintext:
pigz -dc file.maf.gz | maf2bed hg38 | sort -k1,1 -k2,2n | bgzip -@8 > file.bed.gz tabix -p bed file.bed.gz
Give the assembly name exactly as it appears in the MAF, dots included:
maf2bed Species1.1 for Species1.1.chr3 rows.
The output is emitted in input order, which is usually but not always sorted, so
the sort is what makes it safe to tabix. It can be dropped if you already
know the MAF is coordinate-sorted on the reference.
Output format
One line per MAF alignment block:
| col | contents |
|---|---|
| 1 | reference contig, with the <asm>. prefix stripped |
| 2 | start, 0-based (MAF is already 0-based, so this is copied through) |
| 3 | end, exclusive |
| 4 | <asm>_<n>, a block counter over the input |
| 5 | the block's score= from its a line, or 0 if it has none |
| 6 | every row of the block, comma separated, as src:start:size:strand:srcSize:seq |
Column 6 includes the reference row itself, and reports each row's fields
verbatim from the MAF — including a - row's start, which MAF measures on the
reverse complement of the source. Columns 2 and 3 are the exception: a
reverse-strand reference row is flipped to forward coordinates there, because a
BED interval has to be a forward interval.
The summary tier (--summary)
A tabix MAF carries every species' bases on a single BED line, so a zoomed-out
query downloads the whole alignment — JBrowse blocks it with a "too much data"
prompt, and the track has no zoom-out view at all. --summary writes a second,
tiny file that gives it one: per-species presence bars shaded by identity, with
no sequence in it.
It costs nothing extra to produce — it comes out of the same streaming pass:
export LC_ALL=C maf2bed hg38 --summary summary.bed < file.maf \ | sort -k1,1 -k2,2n | bgzip > file.bed.gz tabix -p bed file.bed.gz sort -k1,1 -k2,2n summary.bed | bgzip > summary.bed.gz tabix -p bed summary.bed.gz
Keep the # header line where it is — sort leaves it first because # sorts
below the contig names, tabix -p bed skips it, and JBrowse reads the column
names off it so the track needs no extra configuration. Then point the adapter's
summaryAdapter slot at it:
{
"type": "MafTabixAdapter",
"bedGzLocation": { "uri": "file.bed.gz" },
"index": { "location": { "uri": "file.bed.gz.tbi" } },
"summaryAdapter": {
"type": "BedTabixAdapter",
"bedGzLocation": { "uri": "summary.bed.gz" },
"index": { "location": { "uri": "summary.bed.gz.tbi" } }
}
}This is the same slot UCSC's bigMafSummary.bb goes in, so a track that already
has one does not need this. All four of JBrowse's MAF adapters take the slot and
read the same file, so the summary you write here also serves a BigMafAdapter,
BgzipMafAdapter or BgzipTaffyAdapter track over the same alignment.
What is in it
| col | contents |
|---|---|
| 1 | reference contig |
| 2 | run start, 0-based |
| 3 | run end, exclusive |
| 4 | src, the species id — the row name JBrowse matches on |
| 5 | score, percent identity to the reference over the run, 0..1 |
Consecutive blocks are merged into runs (--merge-gap, default 500
reference bases, matching UCSC hgLoadMafSummary -mergeGap). That merge is what
makes this a tier rather than a second copy of the alignment index: one row per
block per species would multiply the row count by the species count as fast as
dropping the sequence shrinks it, leaving a file only a few times smaller.
Merged runs make it orders smaller, which is what lets a whole chromosome
resolve.
Two details worth knowing:
- Column 4 is the species id, not the full
srctoken —hg38, nothg38.chr1. JBrowse matches it against the track's row names by equality and drops what it cannot find silently, so this has to agree with the browser's own rule, haplotype suffixes (Species1.1.chr3→Species1.1) included. It is also what UCSC'sbigMafSummarycarries. - Column 5 is percent identity, length-weighted across the merged blocks. A
UCSC
bigMafSummaryputs a normalized alignment score here instead. Both are 0..1 and both shade the bar the same way; neither is displayed as a number.
The reference row gets its own runs, scoring 1.0 against itself — the same visual anchor the browser's identity heatmap draws for it.
What it expects of the input
- Every block must be rooted on the assembly you name. Blocks with no row of
that assembly are skipped silently. A MAF from
hal2maf --refGenome <name>,cactus-hal2maf --refGenome <name>, or UCSC multiz is already rooted this way. A MAF whose blocks are each rooted on a different genome, aspggb -Mproduces, needs re-rooting first — see the pggb tutorial. - One BED line per MAF block. If a block carries several reference rows (a
duplication), the block is placed at the first of them; all the copies still
ride along in column 6. If you need one line per reference copy instead, use
maf_to_bed.py, which takes row 0 as the reference and so preserves a split. - Only
aandslines are read.i,e, andqlines are dropped, so the empty-region and quality annotations some MAFs carry do not survive the conversion. - Block numbering in column 4 counts blocks in the input, so skipped blocks leave gaps in the sequence. It is an identifier, not a row count.
- Malformed rows (truncated, or with no alignment text) are dropped rather than
emitted; an unparseable
score=becomes0. Nothing aborts the run.
Loading it in JBrowse 2
{
"type": "MafTrack",
"trackId": "my_maf",
"name": "My multiple alignment",
"assemblyNames": ["hg38"],
"adapter": {
"type": "MafTabixAdapter",
"samples": ["hg38", "panTro6", "mm39"],
"bedGzLocation": { "uri": "file.bed.gz" },
"index": { "indexType": "TBI", "location": { "uri": "file.bed.gz.tbi" } }
}
}samples lists the species in track order; an nhLocation Newick tree can be
given instead, which both supplies the species and orders the rows as a
dendrogram. The adapter keys rows by the species part of src, so at most one
row per species per block is displayed — where a block holds two copies of a
species, the last one wins.
Getting a MAF in the first place
There are many routes; one is to project a pangenome graph or a Cactus HAL onto a reference. The Cactus docs cover this in Using the HAL output. In short:
cactus-hal2maf ./js aln.hal aln.maf.gz --refGenome hg38 --noAncestors --chunkSize 1000000 pigz -dc aln.maf.gz | maf2bed hg38 | sort -k1,1 -k2,2n | bgzip > aln.bed.gz tabix -p bed aln.bed.gz
Cactus recommends cactus-hal2maf over calling hal2maf directly (it
normalizes with taffy and parallelizes). Note that the --dupeMode option the
older docs describe has been replaced by --outType, which takes one or more
of raw, norm (the default), single, single-ref, and consensus.
--outType single greedily filters each species down to one row per block. That
is what most MAF-reading tools want, and JBrowse only draws one row per species
anyway — but it discards paralogous alignments you would otherwise be able to
inspect in column 6. The default norm keeps them, at the cost of a larger file
and duplicate rows that overdraw each other in the display. --outType norm single writes both in one invocation, so you can compare.
Worked end-to-end examples of both graph routes:
- Minigraph-Cactus pangenome —
hal2maf --refGenome, no re-rooting needed - pggb pangenome —
pggb -M, needs re-rooting first
JBrowse also reads bgzipped TAF directly (BgzipTaffyAdapter), and
cactus-hal2maf writes TAF when the output path ends in .taf.gz. TAF is much
smaller than MAF and needs no conversion step, so if you are producing the
alignment yourself and only want it in JBrowse, that route skips maf2bed
entirely. maf2bed is for the case where a MAF is what you have.
A deep TAF still wants a summary, though. A .tai index bounds a read to the
span on screen rather than to the blocks that span lands in, which is the read
cost TAF fixes — but a read costs span times depth, and nothing bounds the
depth. On HPRC's 464-haplotype release-2 alignment a span-bounded TAF read
settles at about 2 compressed bytes per reference base, so a whole chromosome is
still hundreds of megabytes. taffy view back to MAF and pipe it through
--summary to get the tier; the alignment itself stays TAF.
Motivation
I wanted to try the bigMaf (bigBed based) format ecosystem with large MAF files,
but bedToBigBed doesn't support streaming or reading compressed files(?), so
that requires reading big files on disk and in memory. In contrast, the MAF
tabix approach implemented here can be compressed and streamed, which allows
much lower memory usage and disk space.
Footnote
Converted to Rust from Perl as a coding exercise mostly, gaining a modest
speedup on the way https://twitter.com/cmdcolin/status/1719608993310486883. The
original
maf2bed.pl
is still around and its output is interchangeable with this one, if you would
rather not install a Rust toolchain.