Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
24 changes: 24 additions & 0 deletions DEVELOPMENT.md
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,30 @@ for user-facing installation/usage docs.
missing or headerless), even when the block index matches the MAF file
perfectly.

- **The index header carries aggregate counts, and both files of a pair must
get the identical set.** `mafutils index` accumulates blocks / scaffolds /
ref_bases / aln_cols / seq_lines / max_seqs in its existing block loop and
writes them as extra `key=value` tokens, so `mafutils info` can report them
from one line instead of rescanning. Two constraints:
`validate` compares the block and scaffold headers with **exact equality**
(`validate.py`), so block-only keys would make every future validate report
MISMATCH -- the same dict goes into both. And the format stays at **2**:
the change is purely additive, `readIndexHeader` already parses arbitrary
keys, and `tests/test_validate.py` asserts `format == "2"`.
`process_maf_block` returns the integers alongside its stringified row so the
loop never int()s its own output back again; measured cost of the whole
change on a 1.5GB slice, 3 runs each: **4.71s -> 4.77s (~1.3%)**, and a full
42.7GB rebuild produced aggregates matching an independent awk pass exactly
(`blocks=8094203 scaffolds=61 ref_bases=2728222451 aln_cols=3031569326
seq_lines=103950448 max_seqs=15`) with all 8.09M rows byte-identical.
- **`mafutils info` never scans the MAF for counts.** Header first; failing
that, stream the *index* (~52 MB/s) with a warning naming its size; failing
that, report only what needs no index. Species are the exception and are
always labelled a sample -- they are not in the index, and collecting them at
index time costs 2.06x (0.79s -> 1.62s per 1.5GB), which would undo most of
the v0.6.0 speedup. Verified end to end: 1.3s on the 42.7GB file, and the
fallback's scanned values equal the header's exactly.

## Compression / index internals

- **`mafutils/lib/bgzf.py` is vendored from Biopython, not a dependency.**
Expand Down
66 changes: 65 additions & 1 deletion README.md
Original file line number Diff line number Diff line change
Expand Up @@ -3,9 +3,10 @@
`mafutils` is a command-line toolkit for indexing, extracting, and summarizing
MAF (Multiple Alignment Format) files.

It currently provides five commands:
It currently provides six commands:

- `mafutils index`
- `mafutils info`
- `mafutils fetch`
- `mafutils stats`
- `mafutils gc`
Expand Down Expand Up @@ -51,6 +52,12 @@ Create block and scaffold indexes for a MAF (defaults to
mafutils index input.maf
```

Print a quick overview of an indexed MAF (sub-second on any size):

```bash
mafutils info input.maf
```

Fetch trimmed MAF regions from a BED file (index defaults to
`input.maf.block.idx`):

Expand Down Expand Up @@ -154,6 +161,63 @@ Arguments:
| `BLOCK_INDEX` | Output block index path (default: `<MAF_FILE>.block.idx`) |
| `SCAFFOLD_INDEX` | Output scaffold index path (default: `<MAF_FILE>.scaffold.idx`) |

### `mafutils info`

Print a quick overview of a MAF without reading it.

```bash
mafutils info [OPTIONS] MAF_FILE
```

```
MAF : cricetid-15spec.Mus_musculus.nodupes.maf
size : 42.70 GB (none)
index : matches MAF (size + mtime)

ref scaffolds : 61
blocks : 8,094,203
ref bases : 2,728,222,451
aln columns : 3,031,569,326
seq lines : 103,950,448 (all species, including the reference)
max seqs/block : 15

species : 15 (sampled from the first 1,000 block(s) -- possibly incomplete)
Arvicola_amphibius, Cricetulus_griseus, Cricetus_cricetus, ...

For exact per-species presence across the whole file, run `mafutils stats`.
```

Measured on that 42.70 GB file: **1.3 s**, or **0.7 s** with `--sample-blocks 0`.

| Option | Description |
|---|---|
| `--index`, `-i` | Block index (default: `<MAF_FILE>.block.idx`) |
| `--scaffold-index` | Scaffold index (default: `<MAF_FILE>.scaffold.idx`) |
| `--sample-blocks` | Blocks to sample for species names (default 1000); `0` disables the sample and reads no MAF data at all |
| `--all-species` | List every sampled species instead of the first 20 |

**Where the numbers come from.** `mafutils index` records the block counts in
the index header, so `info` normally reads nothing but that one line. Three
tiers, and the MAF is never scanned for these counts:

1. Index built by mafutils 0.7 or later — instant, header only.
2. Older index — falls back to streaming the *index* (~52 MB/s, so ~9 s for a
445 MB index but ~12 min for a 35 GB one), with a warning naming the index
size and recommending a rebuild.
3. No index — reports what needs none and points at `mafutils index`.

**`species` is a sample, not a count.** Species names are not recorded in the
index, and collecting them during indexing would mean parsing every
non-reference `s` line — measured at 2.06x slower indexing. Sampling the first
1,000 blocks costs milliseconds instead. On a whole-genome MAF those blocks all
sit at the start of the first scaffold, so a species absent from that region is
missed. Use `mafutils stats` for an exhaustive per-species breakdown.

**Column meanings.** `ref scaffolds` and `ref bases` are measured on the
*reference* sequence. `seq lines` counts every `s` line including the
reference. `aln columns` is the block width, gap columns included, which all
species share.

### `mafutils fetch`

Fetch regions or scaffolds from a MAF using an existing index.
Expand Down
2 changes: 2 additions & 0 deletions mafutils/cli.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@
from . import fetch as fetch_mod
from . import gc as gc_mod
from . import index as index_mod
from . import info as info_mod
from . import stats as stats_mod
from . import validate as validate_mod

Expand All @@ -19,6 +20,7 @@
app.command("fetch", help="Fetch alignment blocks from a MAF using a BED file and an existing index.")(fetch_mod.fetch_command)
app.command("gc", help="Calculate per-species GC content from a MAF file.")(gc_mod.gc_command)
app.command("index", help="Create block and scaffold indexes for a MAF file.")(index_mod.index_command)
app.command("info", help="Print a quick overview of a MAF file from its index.")(info_mod.info_command)
app.command("stats", help="Summarize an indexed MAF at overall, species, and block levels.")(stats_mod.stats_command)
app.command("validate", help="Check whether a MAF file's index is still trustworthy.")(validate_mod.validate_command)

Expand Down
50 changes: 44 additions & 6 deletions mafutils/index.py
Original file line number Diff line number Diff line change
Expand Up @@ -17,12 +17,22 @@


def process_maf_block(block):
"""
Returns (row, ref_len, aln_cols, num_seqs) for one block: the index row as
strings, plus the three integers the header aggregates need.

The integers are returned rather than re-parsed from `row` on purpose --
this loop runs once per block (233M times on a whole-genome MAF) and is the
hot path the v0.6.0 indexing speedup optimised, so int()-ing its own output
back again would be a needless per-block cost.
"""
ref_seq = block[1].split()
ref_scaff = ref_seq[1].split(".", 1)[1]
line_len = str(len(block[1]))
num_seqs = str(len(block) - 1)
seq_len = str(len(ref_seq[6]))
return [ref_scaff, ref_seq[2], ref_seq[3], seq_len, line_len, num_seqs]
num_seqs = len(block) - 1 # every s line, including the reference
seq_len = len(ref_seq[6]) # block width, gap columns included
row = [ref_scaff, ref_seq[2], ref_seq[3], str(seq_len), line_len, str(num_seqs)]
return row, int(ref_seq[3]), seq_len, num_seqs


def run_index(maf_file, block_index_path, scaffold_index_path):
Expand Down Expand Up @@ -52,10 +62,27 @@ def run_index(maf_file, block_index_path, scaffold_index_path):
# (see iterMafBlocks/openMafHashing).
use_binary = maf_compression in ("none", "gz")

# Header aggregates, accumulated as we go so `mafutils info` can
# report them later from the header alone rather than rescanning.
n_blocks = 0
total_ref_bases = 0
total_aln_cols = 0
total_seq_lines = 0
max_seqs = 0
distinct_scaffolds = set()

for block, block_start, block_end in COMMON.iterMafBlocks(maf_stream, binary=use_binary):
block_info = process_maf_block(block)
block_info, ref_len, aln_cols, num_seqs = process_maf_block(block)
ref_scaffold = block_info[0]

n_blocks += 1
total_ref_bases += ref_len
total_aln_cols += aln_cols
total_seq_lines += num_seqs
if num_seqs > max_seqs:
max_seqs = num_seqs
distinct_scaffolds.add(ref_scaffold)

mdx_line = block_info + [str(block_start), str(block_end)]
block_stream.write("\t".join(mdx_line) + "\n")

Expand All @@ -76,13 +103,24 @@ def run_index(maf_file, block_index_path, scaffold_index_path):

content_hash = f"md5:{hash_obj.hexdigest()}"

# Identical in both headers on purpose -- validate cross-checks the two
# for exact equality, and these describe the MAF, not the index file.
aggregates = {
"blocks": n_blocks,
"scaffolds": len(distinct_scaffolds),
"ref_bases": total_ref_bases,
"aln_cols": total_aln_cols,
"seq_lines": total_seq_lines,
"max_seqs": max_seqs,
}

with open(block_index_path, "w") as block_stream:
COMMON.writeIndexHeader(block_stream, maf_file, maf_compression, size, mtime, content_hash)
COMMON.writeIndexHeader(block_stream, maf_file, maf_compression, size, mtime, content_hash, aggregates)
with open(block_tmp_path, "r") as tmp_fp:
shutil.copyfileobj(tmp_fp, block_stream)

with open(scaffold_index_path, "w") as scaffold_stream:
COMMON.writeIndexHeader(scaffold_stream, maf_file, maf_compression, size, mtime, content_hash)
COMMON.writeIndexHeader(scaffold_stream, maf_file, maf_compression, size, mtime, content_hash, aggregates)
with open(scaffold_tmp_path, "r") as tmp_fp:
shutil.copyfileobj(tmp_fp, scaffold_stream)

Expand Down
Loading
Loading