GitHub

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 src tokenhg38, not hg38.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.chr3Species1.1) included. It is also what UCSC's bigMafSummary carries.
  • Column 5 is percent identity, length-weighted across the merged blocks. A UCSC bigMafSummary puts 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, as pggb -M produces, 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 a and s lines are read. i, e, and q lines 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= becomes 0. 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:

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.

Read the original on github.com ↗