rfmix_reader

Fast, memory-light access to local-ancestry output — RFMix (.msp.tsv, .fb.tsv), FLARE (.anc.vcf.gz) and haptools simulations — as one lazily-evaluated xarray.Dataset backed by a per-chromosome Zarr cache.

• One streaming pass over the source (constant memory), then instant reopen.
• Haplotype ancestry codes stored as int8 (2 bytes per sample-locus, a few MB per chromosome on disk); posteriors kept on request.
• The same Dataset for every format, with position queries, BED intervals, Parquet export, interpolation onto a variant grid, and gnomix-style phasing.

Installation

pip install rfmix-reader                 # core: numpy, pandas, dask, xarray, zarr, cyvcf2
pip install "rfmix-reader[viz,io]"       # + matplotlib/seaborn/cairosvg, pyarrow

Parsing is CPU-only by design (it is I/O bound). When CuPy is importable the imputation and plotting helpers use it as the array backend (rfmix_reader.backends.use_gpu()).

Quickstart

from rfmix_reader import open_rfmix, open_local_ancestry

# parse once (one streaming pass), cache as la_cache/<chrom>.zarr
ds = open_rfmix("two_pops/out/", cache_dir="la_cache/")

# later sessions: instant, lazy
ds = open_local_ancestry("la_cache/")
ds = open_local_ancestry("la_cache/", chrom="21")

ds.la.counts            # (variant, sample, ancestry) int8 diploid counts, lazy dask
ds.la.haplotypes        # (variant, sample, ploidy) int8 ancestry codes
ds.la.global_ancestry   # DataFrame: sample_id, <ancestries>, chrom
ds.la.samples, ds.la.ancestries, ds.la.chromosomes

rfmix-reader convert fb two_pops/out/ la_cache/ --keep-posteriors
rfmix-reader info la_cache/

Readers

path may be a directory, a single file, or a path prefix (/out/run_ matches /out/run_chr1.msp.tsv, ...). chrom="21" restricts to one chromosome. Without cache_dir the Dataset is built in memory, which is fine for .msp.tsv; for .fb.tsv with posteriors use a cache (a chr1-sized cohort of 500 samples is ~11 GB of float32 posteriors).

The Dataset

dims:   variant, sample, ploidy (=2), ancestry, contig
vars:   haplotype_ancestry (variant, sample, ploidy)            int8  code into `ancestry`, -1 missing
        posterior          (variant, sample, ploidy, ancestry)  float32  optional
        global_ancestry    (contig, sample, ancestry)           float32
coords: chromosome, variant_position, segment_end (variant); sample_id; ancestry; contig
attrs:  source_format, source_files, rfmix_reader_version

Conventions that hold for every reader:

• Ancestries are in the tool's own order (RFMix #reference_panel_population / #Subpopulation order/codes, FLARE ##ANCESTRY), and global_ancestry uses the same order.
• A sample/locus with no call (e.g. an all-zero RFMix posterior) is -1 in the codes and in every count of that row.
• Counts are derived lazily from the codes (ds.la.counts), so they are never stored twice.

ds.la is the accessor (registered on import of rfmix_reader.core, which every open_* call does):

Operations

All operations are lazy or streaming; sample is an index or a sample ID.

import pandas as pd

loci = pd.DataFrame({"chrom": ["chr21", "chr21"], "pos": [15_000_000, 30_000_000]})
ds.la.at_positions(loci)                       # AFR_haplotypes, AFR_fraction, ...
ds.la.at_positions(loci, samples=["NA19700"], aggregate=False)

# QTL-style access: one chromosome in memory, then index by genotype position
chrom = ds.la.sel_chrom("chr21")                        # contiguous slice, still lazy
counts = chrom.la.counts.values                          # (variants, samples, ancestries) int8
idx = chrom.la.locus_index("chr21", genotype_positions)  # segment per SNP, -1 = uncovered
H = counts[idx[idx >= 0]]

bed = ds.la.to_bed("NA19700", min_segment=3)
ds.la.to_parquet("out/", prefix="la", rows_per_file=100_000)

variants = pd.read_parquet("genotypes/variants.parquet")   # chrom, pos
dense = ds.la.interpolate(variants, "imputed/", method="stepwise")

Phasing

ds.la.phase() corrects switch errors between the two haplotypes per sample the way gnomix does: heterozygous blocks are found from the haplotype codes, each window inside a block is scored by how well the two posterior tracks match the block-start orientation versus the swapped one, uninformative windows inherit the previous state, and the two haplotypes (codes and posteriors) are exchanged wherever the track says "switched". It needs no reference panel; open .fb.tsv output with keep_posteriors=True so the posteriors are available (without them the hard calls are used). The result adds a phase_swapped (variant, sample) mask.

from rfmix_reader import open_rfmix
from rfmix_reader.processing.phase import PhasingConfig, merge_phased_zarrs

ds = open_rfmix("two_pops/out/", source="fb", keep_posteriors=True,
                cache_dir="la_cache/", chrom="21")
phased = ds.la.phase(config=PhasingConfig(window_size=50, min_block_len=20, posterior_margin=0.2))
phased["phase_swapped"].sum("variant")         # exchanged loci per sample

Whole chromosomes phase in seconds (chr21, 175k loci by 500 samples: about 4 s). phase_rfmix_chromosome_to_zarr(prefix, None, None, "phased_chr21.zarr", chrom="21") does read, phase and write in one call, and merge_phased_zarrs concatenates per-chromosome stores; both outputs reopen with open_local_ancestry.

The previous reference-panel matcher remains as method="reference" (with ref_zarr_root / sample_annot_path) for comparison only: it compares ancestry labels against reference allele codes, which is not a sound test of phase. To build its inputs, convert bgzipped, indexed reference VCFs with prepare-reference refs/ 1kg_chr21.vcf.gz (needs the reference extra) and give a two-column sample_id<TAB>group annotation table.

Haptools simulations

haptools does not write the chromosome length into the ##contig header, but the tabix region pulls need it. Reheader each file with the contig entry from the contigs.txt haptools produces, for example:

CONTIG_LINE=$(grep -w "ID=${CHR}" contigs.txt)
bcftools view -h "$IN" | sed "s/^##contig=<ID=${CHR}>.*/${CONTIG_LINE}/" > header.${CHR}.tmp
bcftools reheader -h header.${CHR}.tmp -o "$OUT" "$IN"
tabix -p vcf "$OUT"

Visualization

from rfmix_reader import plot_global_ancestry, plot_ancestry_by_chromosome, plot_local_ancestry_tagore

g_anc = ds.la.global_ancestry
plot_global_ancestry(g_anc, save_path="global")            # writes global.png and global.pdf
plot_ancestry_by_chromosome(g_anc, save_path="by_chrom")
plot_local_ancestry_tagore(ds.la.to_tagore("NA19700"), prefix="NA19700", build="hg38", oformat="png")

Migrating from 0.x

See MIGRATION.md. Every legacy (0.5 and earlier) name maps to one call on the Dataset; ds.la.to_legacy() returns the old triple and rfmix_reader.from_legacy(loci_df, g_anc, admix) builds a Dataset from one.

Development

git clone https://github.com/heart-gen/rfmix_reader.git
cd rfmix_reader
poetry install --with test --extras "viz io"
poetry run pytest             # fast suite (seconds)
poetry run pytest --run-slow  # also the chr21 tests (needs git-LFS data)

Test fixtures for every format live under tests/data/ (regenerate with python tests/data/make_fixtures.py).

Citation

If you use this software, please cite:

Benjamin, K. J. M. (2024). RFMix-reader [Computer software]. https://github.com/heart-gen/rfmix_reader

Kynon JM Benjamin. "RFMix-reader: Accelerated reading and processing for local ancestry studies." bioRxiv (2024). DOI: 10.1101/2024.07.13.603370.

Funding

This work was supported by the National Institutes of Health, National Institute on Minority Health and Health Disparities (NIMHD) K99MD016964 / R00MD016964.