Reading files

This page covers how polars-bio loads bioinformatics files — the supported formats, the performance machinery shared across them (indexed reads, predicate/projection pushdown, parallel reads), per-format specifics, coordinate-system handling, and the metadata attached to every DataFrame. Files can be read from local disk or streamed directly from cloud storage.

On this page: File formats · Performance features · Format-specific notes · Schema inspection · Coordinate systems · File metadata

File formats support

For every bioinformatic format there are always three methods available: read_* (eager), scan_* (lazy) and register_* that can be used to either read the file into a Polars DataFrame/LazyFrame or register it as a DataFusion table for further processing using SQL or built-in interval methods. In either case, local and/or cloud storage files can be used as an input. Please refer to the cloud storage section for more details.

The matrix below summarizes which performance features each format supports. Format-specific options and behaviors are documented under format-specific notes.

Format	Single-threaded	Parallel (indexed)	Limit pushdown	Predicate pushdown	Projection pushdown
BED		❌		❌	❌
VCF		(TBI/CSI)
BCF		(CSI)
VCF Zarr		(region index)	❌
BAM		(BAI/CSI)
CRAM		(CRAI)
FASTQ		(GZI)		❌
FASTA		❌		❌	❌
GFF3		(TBI/CSI)
GTF		(TBI/CSI)
Pairs		(TBI/CSI)
BGEN		(BGI)
PGEN		(embedded/PGI)
BigWig		(built-in BBI index)	❌
BigBed		(built-in BBI index)	❌
Cooler (.cool/.mcool)		(built-in CSR index)	❌	(first axis)

Performance features

polars-bio applies the same performance machinery — indexing, pushdown, and parallel reads — across most formats. The capability matrix above shows which format supports what; this section explains each feature and how to use it.

Progress output

When polars-bio streams batches from DataFusion into Polars, it displays the number of processed rows and throughput on stderr. To suppress this progress output for the whole process, set TQDM_DISABLE=1 when starting Python:

TQDM_DISABLE=1 python analysis.py

You can also set the environment variable in Python. Set it at the start of the process, before importing polars_bio or any other module that imports tqdm, because tqdm reads environment overrides when it is imported:

import os

os.environ["TQDM_DISABLE"] = "1"

import polars_bio as pb

Indexed reads & random access

When an index file is present alongside the data file (BAI/CSI for BAM, CRAI for CRAM, TBI/CSI for VCF, CSI for BCF, and TBI/CSI for GFF and Pairs), polars-bio can push genomic region filters down to the DataFusion execution layer. This enables index-based random access — only the relevant genomic regions are read from disk, dramatically improving performance for selective queries on large files.

Index files are auto-discovered by convention. Predicate pushdown is enabled by default for BAM, CRAM, VCF, BCF, GFF, and Pairs formats — no extra configuration is needed.

Supported index formats

Unmapped reads and indexed scans

A whole-file scan returns the same records whether or not an index is present, including the unplaced, unmapped reads at the end of a file. Those reads carry no reference, so they surface with a null chrom:

import polars as pl
import polars_bio as pb

# sample.cram holds 300 mapped reads and 200 unmapped ones
pb.scan_cram("sample.cram").collect().height    # 500, with or without sample.cram.crai

# the unmapped ones
(
    pb.scan_cram("sample.cram")
    .filter(pl.col("chrom").is_null())
    .collect()
)                                               # 200

A region query asks for placed reads by definition, so it never returns the unmapped tail:

pb.scan_cram("sample.cram").filter(pl.col("chrom") == "chr1").collect().height   # 150

BAM behaves the same way, and always has.

Region queries with the DataFrame API

Simply use .filter() — predicate pushdown is enabled by default for BAM, CRAM, VCF, BCF, GFF, and Pairs:

import polars as pl
import polars_bio as pb

# Single chromosome filter — only chr1 data is read from disk
df = (
    pb.scan_bam("alignments.bam")
    .filter(pl.col("chrom") == "chr1")
    .collect()
)

# Multi-chromosome filter
df = (
    pb.scan_vcf("variants.vcf.gz")
    .filter(pl.col("chrom").is_in(["chr21", "chr22"]))
    .collect()
)

# BCF has a dedicated lazy API; a neighboring variants.bcf.csi is automatic
df = (
    pb.scan_bcf("variants.bcf")
    .filter(
        (pl.col("chrom") == "chr21")
        & (pl.col("start") >= 1_000_000)
        & (pl.col("start") <= 2_000_000)
    )
    .collect()
)

# Region query — combines chromosome and coordinate filters
df = (
    pb.scan_bam("alignments.bam")
    .filter(
        (pl.col("chrom") == "chr1")
        & (pl.col("start") >= 10000)
        & (pl.col("end") <= 50000)
    )
    .collect()
)

# CRAM with predicate pushdown
df = (
    pb.scan_cram("alignments.cram")
    .filter(pl.col("chrom") == "chr1")
    .collect()
)

Region queries with SQL

The SQL path works automatically — DataFusion parses the WHERE clause and uses the index without any extra flags:

import polars_bio as pb

pb.register_bam("alignments.bam", "reads")

# Single chromosome
result = pb.sql("SELECT * FROM reads WHERE chrom = 'chr1'").collect()

# Region query
result = pb.sql(
    "SELECT * FROM reads WHERE chrom = 'chr1' AND start >= 10000 AND \"end\" <= 50000"
).collect()

# Combined genomic and record filters
result = pb.sql(
    "SELECT * FROM reads WHERE chrom = 'chr1' AND mapping_quality >= 30"
).collect()

Predicate pushdown

All formats support record-level predicate evaluation — filters on columns like mapping_quality, flag, or strand are evaluated per-record during the scan, with or without an index file. When an index is present, genomic-coordinate filters additionally drive index-based random access.

See the Developers Guide for the translation pipeline internals and examples.

Projection pushdown

BAM, CRAM, VCF, BCF, and Pairs formats support parsing-level projection pushdown — unprojected fields are skipped entirely during record parsing. Enabled by default (projection_pushdown=True). See the Developers Guide for internals and execution plan inspection.

Parallel reads & partitioning

This section covers how reading a file is partitioned. The degree of parallelism is the same global datafusion.execution.target_partitions knob that governs every operation — see Parallel engine for the setting itself and its default.

When an index file is present, DataFusion distributes genomic regions across balanced partitions using index-derived size estimates, enabling parallel execution. Formats with known contig lengths (BAM, CRAM) can split large regions into sub-regions for full parallelism even on single-chromosome queries. For FASTQ files, a GZI index alongside a BGZF-compressed file enables parallel decoding of compressed blocks. This is controlled by the global target_partitions setting:

import polars_bio as pb

pb.set_option("datafusion.execution.target_partitions", "8")
df = pb.read_bam("large_file.bam")  # 8 partitions will be used for parallel execution
df = pb.read_fastq("reads.fastq.bgz")  # parallel BGZF decoding when .gzi index is present

Partitioning behavior (BAM, CRAM, VCF, BCF, GFF):

Partitioning behavior (FASTQ):

Generating index files

# BAM: sort and index
samtools sort input.bam -o sorted.bam
samtools index sorted.bam                # creates sorted.bam.bai

# CRAM: sort and index
samtools sort input.cram -o sorted.cram --reference ref.fa
samtools index sorted.cram               # creates sorted.cram.crai

# VCF: sort, compress, and index
bcftools sort input.vcf -Oz -o sorted.vcf.gz
bcftools index -t sorted.vcf.gz          # creates sorted.vcf.gz.tbi

# BCF: coordinate-sort and create a CSI index
bcftools sort input.bcf -Ob -o sorted.bcf
bcftools index --csi sorted.bcf           # creates sorted.bcf.csi

# GFF: sort, compress, and index
(grep "^#" input.gff; grep -v "^#" input.gff | sort -k1,1 -k4,4n) | bgzip > sorted.gff.gz
tabix -p gff sorted.gff.gz               # creates sorted.gff.gz.tbi

# Pairs: sort, compress, and index (col 2=chr1, col 3=pos1)
sort -k2,2 -k3,3n contacts.pairs | bgzip > contacts.pairs.gz
tabix -s 2 -b 3 -e 3 contacts.pairs.gz   # creates contacts.pairs.gz.tbi

# FASTQ: BGZF compress and create GZI index for parallel reads
bgzip reads.fastq                         # creates reads.fastq.bgz
bgzip -r reads.fastq.bgz                 # creates reads.fastq.bgz.gzi

Format-specific notes

Most formats work through the generic read_*/scan_*/register_* API with no extra options. The formats below expose additional capabilities or behaviors worth knowing about.

BGEN

BGEN 1.2 and 1.3 genotype files use read_bgen / scan_bgen, with register_bgen and describe_bgen for registration and schema inspection. One row is one BGEN variant. Encoded alleles stay ordered in alleles and are not assigned reference/alternate semantics, because BGEN does not define them.

• Genotype output — genotype_output="probability" (default) emits genotypes.GP, preserving every format-defined probability state, and genotypes.PLOIDY, the declared ploidy of each selected sample. genotype_output="dosage" emits genotypes.DS instead, the expected copy count of alleles[1], and rejects multiallelic variants.
• Probability layout — for an unphased biallelic record GP holds one value per allele-count state; for a phased record it holds one vector per haplotype, haplotype-major. phased and bits are columns, so a mixed file stays interpretable row by row.
• Probability storage — probability_layout="nested" (default) stores each sample's states as a variable-length list, which every BGEN file can use. probability_layout="fixed" stores them as a fixed-width list instead, dropping the per-sample offsets that are about a quarter of the emitted probability bytes for a diploid biallelic cohort. It requires every variant to store the same number of states and rejects a file that mixes them, so reach for it on whole-cohort imputed data and stay on the default otherwise. The option has no effect when genotype_output="dosage".
• Indexes — a neighbouring cohort.bgen.bgi is auto-discovered and used to push chrom, id, rsid, start, and end predicates into the scan. Pass bgi_path for an index stored elsewhere. Without an index, the provider builds a transient in-memory catalog by scanning variant metadata.
• Samples — identifiers come from the embedded sample block, an explicit sample_path, or generated sample_1…sample_N names. samples=[...] selects and reorders the emitted samples, and the emitted order is available via meta["header"]["sample_names"].
• Genotype fields — genotype_fields selects children of the genotypes struct by name, from the output mode's value child — "DS" for dosage, "GP" for probability — and "PLOIDY", emitted in the requested order. All of them are emitted by default. "PLOIDY" is a byte per genotype, 2.53 GB on a whole 1000 Genomes chromosome 22, and a NumPy view of the result keeps the whole struct alive, so pass ["DS"] when only the dosages are wanted.
• Projection — a scan that selects only metadata columns never reads or decompresses probability blocks.

import polars_bio as pb

dosage = pb.scan_bgen("cohort.bgen", genotype_output="dosage")
schema = pb.describe_bgen("cohort.bgen")

# Uniform-width probabilities without the per-sample offsets.
probabilities = pb.scan_bgen("cohort.bgen", probability_layout="fixed")

# Only the dosages, without the per-genotype ploidy byte.
dosage_only = pb.scan_bgen(
    "cohort.bgen", genotype_output="dosage", genotype_fields=["DS"]
)

For a whole-cohort dosage matrix, use read_bgen_matrix instead of consolidating Arrow batches yourself. It decodes each variant straight to its final address in a dense NumPy array, so the dosages are never copied, and rows stay in file order at every thread count:

matrix = pb.read_bgen_matrix("chr22.bgen")
matrix.values.shape         # (variants, samples), float32
matrix.values.mean(axis=1)  # per-variant mean dosage
matrix.positions            # one per row
matrix.sample_names         # one per column

Dosage only: BGEN probabilities are variable width and have no single dense shape, so read those with scan_bgen(genotype_output="probability").

BGEN is an input format; polars-bio does not write it.

PGEN

PLINK 2 filesets use read_pgen / scan_pgen, with register_pgen and describe_pgen for registration and schema inspection. One row is one PVAR variant, with ref and a list-typed alt carrying the declared alleles.

For a whole-cohort genotype matrix — what association testing, PCA and relatedness pipelines consume — use read_pgen_matrix instead. It returns a dense NumPy array rather than a DataFrame, and decodes straight into it, so the values are never copied:

matrix = pb.read_pgen_matrix("cohort.pgen", field="ALT_COUNT")
matrix.values.shape         # (variants, samples), int8
matrix.values.mean(axis=1)  # per-variant ALT frequency x 2
matrix.positions            # one per row
matrix.sample_names         # one per column

Rows are in PVAR order at every partition count, which a read_pgen scan does not guarantee. field takes "ALT_COUNT" or "DS"; fields with more than one value per sample have no dense form.

• Genotype fields — genotype_fields selects children of the genotypes struct by name, from "GT", "ALT_COUNT", "PHASED", "DS", "DS_STORED", and "HDS", emitted in the requested order. It defaults to ("GT",). This narrows the provider default, which emits all of them: reading several representations of the same genotypes is rarely what you want, so ask for the others explicitly. "ALT_COUNT" is the hardcall ALT allele count as int8, one byte per genotype rather than the four "DS" uses; prefer it when the fileset stores only hardcalls.
• Companions — the .pvar (then .pvar.zst) and .psam are discovered from the .pgen basename. Pass pvar_path, psam_path, or pgi_path for companions stored elsewhere.
• Samples — selectable names are built from PSAM identifiers under psam_id_mode: "iid" (default) uses IID alone and rejects duplicates, "fid_iid" uses FID:IID, and "fid_iid_sid" uses FID:IID:SID. A PSAM without FID or SID columns defaults those parts to 0. samples=[...] selects and reorders the emitted samples; missing_sample_policy="ignore" drops requested names the PSAM does not have, instead of raising.
• Read coalescing — max_range_gap bounds the run of unselected bytes bridged when merging reads. It defaults to 0, so a subset scan issues one read per contiguous run of selected variants. Raising it trades wasted bytes for fewer requests, which matters most on object storage. max_range_bytes and batch_soft_byte_limit bound the largest coalesced read and the soft genotype-byte target per batch. Leave any of them unset to keep the provider default.
• Projection — a scan that selects only metadata columns never reads genotype records.

import polars_bio as pb

dosages = pb.scan_pgen("cohort.pgen", genotype_fields=["DS"])
schema = pb.describe_pgen("cohort.pgen")

# Bridge gaps up to 64 KiB to cut request count on object storage.
subset = pb.scan_pgen(
    "cohort.pgen", samples=["NA12878", "NA12879"], max_range_gap=65536
)

PGEN is an input format; polars-bio does not write it.

VCF, BCF, and VCF Zarr

Text VCF (plain or compressed) uses read_vcf / scan_vcf, while binary BCF uses the dedicated read_bcf / scan_bcf methods. Registration and schema inspection follow the same split: use register_vcf / describe_vcf for text VCF and register_bcf / describe_bcf for BCF. Local VCF Zarr stores use the corresponding *_vcf_zarr functions. Default BCF output has the same rows, columns, data types, INFO handling, and FORMAT layout as an equivalent VCF. Key behaviors:

• Lazy BCF scans — use scan_bcf("cohort.bcf") to retain streaming execution, projection and predicate pushdown, and CSI-backed parallel partition processing. read_bcf is the eager convenience wrapper.
• BCF indexes — cohort.bcf.csi is auto-discovered. Genomic filters use CSI byte-range reads; full scans use CSI regions to fill datafusion.execution.target_partitions. Without CSI, BCF falls back to one sequential input partition.
• INFO fields — by default (info_fields=None) all header INFO fields are available in the schema. Pass an explicit list to select a subset, or info_fields=[] to exclude INFO columns entirely.
• Single-sample FORMAT — FORMAT fields are exposed as top-level columns (GT, DP, GQ, ...).
• Multisample FORMAT — exposed as a nested genotypes column (struct<GT: list, DP: list, ...>), where each FORMAT field is a list of values ordered by sample. Sample names are available via meta["header"]["sample_names"].
• Sample subset selection — pass samples=[...] to the corresponding VCF or BCF read/scan method to keep only selected samples in the nested genotypes output. Missing sample names are skipped with a warning.
• FORMAT metadata fidelity — meta["header"]["format_fields"] preserves each FORMAT field's number / type / description.
• Range operations — .bcf paths are accepted anywhere the range APIs accept .vcf paths and use the same logical VCF schema.

import polars as pl
import polars_bio as pb

# INFO selection: all fields (default) vs none
df_all_info = pb.read_vcf("variants.vcf")                  # all INFO fields
df_no_info  = pb.read_vcf("variants.vcf", info_fields=[])  # no INFO columns

# Multisample FORMAT is exposed as a nested `genotypes` column
df = pb.read_vcf("multisample.vcf", format_fields=["GT", "DP"])
df.select(["chrom", "start", "genotypes"])

# Restrict the nested genotypes output to selected samples
df_subset = pb.read_vcf(
    "multisample.vcf",
    format_fields=["GT"],
    samples=["NA12880", "NA12878"],
)

# BCF remains lazy until collect()
rare = (
    pb.scan_bcf("cohort.bcf", info_fields=["AF"], format_fields=["GT"])
    .filter((pl.col("chrom") == "chr21") & (pl.col("AF").list.first() < 0.01))
    .select(["chrom", "start", "ref", "alt", "AF", "genotypes"])
    .collect()
)

# Optional typed biallelic dosage avoids materializing GT strings. Multisample
# input returns nullable Int8 ALT-allele counts in genotypes.GT; single-sample
# input preserves the FORMAT layout and returns a top-level nullable Int8 GT.
dosage = pb.scan_bcf(
    "cohort.bcf",
    format_fields=["GT"],
    genotype_output="dosage",
)

On read_bcf and scan_bcf, genotype_output="string" is the default and preserves VCF-compatible GT values. genotype_output="dosage" requires GT to be the only selected FORMAT field and returns the number of ALT alleles per sample as nullable Int8—normally 0, 1, or 2 for diploid calls. When format_fields is omitted, all header-defined FORMAT fields are selected, so pass format_fields=["GT"] when the header declares additional fields. A fully or partially missing GT becomes null. Multiallelic records are rejected instead of silently collapsing non-reference alleles. Text VCF methods do not expose this BCF-only option. Neither the VCF nor BCF methods expose genotype_encoding_raw; that option remains specific to read_vcf_zarr and scan_vcf_zarr.

BCF is currently an input format. write_vcf / sink_vcf still write text VCF (optionally gzip/BGZF compressed). It shares the internal logical VCF schema, while BCF scan metadata reports source_format="bcf"; source_path retains the .bcf path.

BAM, SAM and CRAM

polars-bio supports reading BAM, SAM, and CRAM optional alignment tags as individual columns. Tags are only parsed when explicitly requested, ensuring zero overhead for standard reads.

Reading optional tags

import polars_bio as pb

# Read BAM with specific tags
df = pb.read_bam(
    "alignments.bam",
    tag_fields=["NM", "AS", "MD"]  # Edit distance, alignment score, mismatch string
)

# Tags appear as regular columns
print(df.select(["name", "chrom", "NM", "AS"]))

# Lazy scan with tag filtering
lf = pb.scan_bam("alignments.bam", tag_fields=["NM", "AS"])
high_quality = lf.filter((pl.col("NM") <= 2) & (pl.col("AS") >= 100)).collect()

# SQL queries (tags must be quoted)
pb.register_bam("alignments.bam", "reads", tag_fields=["NM", "RG"])
result = pb.sql('SELECT name, "NM" FROM reads WHERE "NM" <= 2').collect()

# Exact type hints for custom or array tags
typed = pb.read_bam(
    "alignments.bam",
    tag_fields=["tp", "ML", "FZ"],
    infer_tag_types=False,
    tag_type_hints=["tp:A", "ML:B:C", "FZ:B:S"],
)

tag_type_hints accepts scalar forms such as NM:i, de:f, tp:A, XH:H, plus array forms TAG:B and TAG:B:SUBTYPE such as ML:B:C or FZ:B:S. Bare TAG:B is treated as the default integer-array hint and normalized to TAG:B:i internally, so it reads back as list[i32].

Common tags

• NM (Int32): Edit distance to reference
• MD (Utf8): Mismatch positions string
• AS (Int32): Alignment score
• XS (Int32): Secondary alignment score
• RG (Utf8): Read group identifier
• CB (Utf8): Cell barcode (single-cell)
• UB (Utf8): UMI barcode (single-cell)

Full registry includes ~40 common SAM tags.

Tag reading performance

• Zero overhead when tag_fields=None (default)
• Projection pushdown: only selected tags are parsed
• Tags parsed once per batch, not per record

BigWig and BigBed

BigWig (continuous signal) and BigBed (feature intervals) are supported through the same eager/lazy/register access patterns. Predicate pushdown on the genomic coordinate columns and projection pushdown are enabled by default.

Parallel scans use each file's built-in cir-tree (R-tree) index, so BigWig and BigBed do not need a sidecar index. The reader balances compressed data blocks across up to datafusion.execution.target_partitions independent partitions; small files can expose fewer partitions than requested when they contain fewer independent indexed work units.

import polars as pl
import polars_bio as pb

# Lazy scan with a genomic range filter (predicate pushdown)
signal = (
    pb.scan_bigwig("signal.bw")
    .filter(pl.col("chrom") == "chr1")
    .collect()
)

# Eager read
features = pb.read_bigbed("features.bb")

# Register as a DataFusion table for SQL
pb.register_bigwig("signal.bw", "signal")
pb.sql("SELECT chrom, start, `end`, value FROM signal WHERE chrom = 'chr1'").collect()

Cooler (.cool/.mcool)

Cooler files store Hi-C contact matrices as HDF5 containers: .cool holds a single resolution, .mcool nests one data collection per resolution. polars-bio reads them natively (statically linked HDF5, no cooler/h5py dependency) through the same eager/lazy/register access patterns.

Each row is one stored (upper-triangle) pixel, joined with bin coordinates by default: chrom1, start1, end1, chrom2, start2, end2, count. The stored count dtype is preserved as Int32, Int64, UInt32, UInt64, or Float64, while joined start/end coordinates use UInt64 so values beyond the 32-bit range remain lossless. join_bins=False returns the raw COO triple (bin1_id, bin2_id, count), and include_weights=True adds the weight1/weight2 balancing weights of a balanced cooler (NaN marks bins filtered out by balancing). Coordinates are natively 0-based half-open and follow the standard use_zero_based handling.

Select an .mcool resolution with the resolution argument or the cooler URI syntax file.mcool::/resolutions/10000; a multi-resolution file with no selection raises an error listing the stored resolutions. describe_cool lists every data collection (resolution, nbins, nnz, assembly) without touching pixel data.

Predicate pushdown prunes pixel row ranges through the cooler chrom_offset/bin1_offset CSR indexes for filters on the first axis (chrom1, start1, end1); second-axis and count filters are applied client-side. Projection pushdown reads only the HDF5 datasets the requested columns need, and count(*) is served from the index without reading pixels. Parallel scans split the (pruned) pixel row space along bin1 boundaries across datafusion.execution.target_partitions; note that libhdf5 serializes raw reads behind a global lock, so parallel speedups flatten beyond a few partitions. Only local filesystem paths are supported in this version, and the single-cell .scool layout is not supported.

import polars as pl
import polars_bio as pb

# Available resolutions, without scanning pixels
pb.describe_cool("contacts.mcool")

# Lazy scan of one resolution with a first-axis region filter (pushdown)
cis = (
    pb.scan_cool("contacts.mcool", resolution=10000)
    .filter((pl.col("chrom1") == "chr2") & (pl.col("start1") >= 20_000_000))
    .collect()
)

# Balanced contact counts, computed client-side from the exposed weights
balanced = (
    pb.scan_cool("contacts.mcool", resolution=10000, include_weights=True)
    .with_columns((pl.col("count") * pl.col("weight1") * pl.col("weight2")).alias("balanced"))
    .collect()
)

# Hi-C end to end: raw pairs and the binned matrix side by side
pairs = pb.scan_pairs("sample.pairs.gz")
n_cis_pairs = (
    pairs.filter(pl.col("chr1") == pl.col("chr2")).select(pl.len()).collect().item()
)
n_cis_pixels = (
    pb.scan_cool("contacts.mcool", resolution=10000)
    .filter(pl.col("chrom1") == pl.col("chrom2"))
    .select(pl.len())
    .collect()
    .item()
)

# Register as a DataFusion table for SQL
pb.register_cool("contacts.mcool", "hic", resolution=10000)
pb.sql(
    "SELECT chrom1, start1, count FROM hic WHERE chrom1 = 'chr2' ORDER BY count DESC LIMIT 10"
).collect()

Schema inspection

Quickly inspect BAM/CRAM file schemas without reading the entire file:

import polars_bio as pb
import polars as pl

# Get schema information for BAM file
schema = pb.describe_bam("file.bam")
print(schema)
# shape: (11, 2)
# ┌─────────────────┬──────────┐
# │ column          ┆ datatype │
# │ ---             ┆ ---      │
# │ str             ┆ str      │
# ╞═════════════════╪══════════╡
# │ name            ┆ String   │
# │ chrom           ┆ String   │
# │ start           ┆ UInt32   │
# ...

# Include tag columns in schema
schema = pb.describe_bam("file.bam", tag_fields=["NM", "AS", "MD"])
print(schema)  # Shows 14 columns including tags

# CRAM schema
schema = pb.describe_cram("file.cram")

# VCF, BCF, and local VCF Zarr describe output includes INFO and FORMAT rows.
# Nested FORMAT data is reported by its selectable column name, `genotypes`.
vcf_schema = pb.describe_vcf("variants.vcf")
bcf_schema = pb.describe_bcf("variants.bcf")
vcz_schema = pb.describe_vcf_zarr("cohort.vcz")
format_fields = vcf_schema.filter(pl.col("field_type") == "FORMAT")

# VCF/BCF describe columns: name, field_type, data_type, description.

Coordinate systems support

polars-bio supports both 0-based half-open and 1-based closed coordinate systems for genomic ranges operations. By default, it uses 1-based closed coordinates, which is the native format for VCF, GFF, and SAM/BAM files.

How it works

The coordinate system is managed through DataFrame metadata that is set at I/O time and read by range operations. This ensures consistency throughout your analysis pipeline.

flowchart TB
    subgraph IO["I/O Layer"]
        scan["scan_vcf/bcf/gff/bam/cram/bed()"]
        read["read_vcf/bcf/gff/bam/cram/bed()"]
    end

    subgraph Config["Session Configuration"]
        zero_based["datafusion.bio.coordinate_system_zero_based<br/>(default: false = 1-based)"]
        check["datafusion.bio.coordinate_system_check<br/>(default: false = lenient)"]
    end

    subgraph DF["DataFrame with Metadata"]
        polars_meta["Polars DataFrame/LazyFrame<br/>coordinate_system_zero_based"]
        pandas_meta["Pandas DataFrame<br/>df.attrs"]
    end

    subgraph RangeOps["Range Operations"]
        overlap["overlap()"]
        nearest["nearest()"]
        count["count_overlaps()"]
        coverage["coverage()"]
        merge["merge()"]
        cluster["cluster()"]
        complement["complement()"]
        subtract["subtract()"]
    end

    subgraph Validation["Metadata Validation"]
        validate["validate_coordinate_systems()"]
        error1["MissingCoordinateSystemError"]
        error2["CoordinateSystemMismatchError"]
        fallback["Fallback to global config<br/>+ emit warning"]
    end

    scan --> |"sets metadata"| polars_meta
    read --> |"sets metadata"| polars_meta
    zero_based --> |"use_zero_based param<br/>or default"| scan
    zero_based --> |"use_zero_based param<br/>or default"| read

    polars_meta --> overlap
    polars_meta --> nearest
    polars_meta --> count
    polars_meta --> coverage
    polars_meta --> merge
    polars_meta --> cluster
    polars_meta --> complement
    polars_meta --> subtract
    pandas_meta --> overlap

    overlap --> validate
    nearest --> validate
    count --> validate
    coverage --> validate
    merge --> validate
    cluster --> validate
    complement --> validate
    subtract --> validate

    validate --> |"metadata missing"| check
    validate --> |"metadata mismatch"| error2
    check --> |"true (strict)"| error1
    check --> |"false (lenient)"| fallback
    fallback --> zero_based

Session parameters

polars-bio provides two session parameters to control coordinate system behavior:

import polars_bio as pb

# Check current settings
print(pb.get_option(pb.POLARS_BIO_COORDINATE_SYSTEM_ZERO_BASED))  # "false"
print(pb.get_option(pb.POLARS_BIO_COORDINATE_SYSTEM_CHECK))       # "false"

# Change to 0-based coordinates globally
pb.set_option(pb.POLARS_BIO_COORDINATE_SYSTEM_ZERO_BASED, True)

Reading files with coordinate system metadata

When you read genomic files using polars-bio I/O functions, the coordinate system metadata is automatically set on the returned DataFrame:

import polars_bio as pb

# Default: 1-based coordinates (use_zero_based=False)
df = pb.scan_vcf("variants.vcf")
# Metadata is automatically set: coordinate_system_zero_based=False

# Explicit 0-based coordinates
df_zero = pb.scan_bed("regions.bed", use_zero_based=True)
# Metadata is automatically set: coordinate_system_zero_based=True

# Range operations read coordinate system from metadata
result = pb.overlap(df, df_zero, ...)  # Raises CoordinateSystemMismatchError!

Setting metadata on DataFrames

For DataFrames not created via polars-bio I/O functions, you must set the coordinate system metadata manually:

import polars as pl

# Create a DataFrame
df = pl.DataFrame({
    "chrom": ["chr1", "chr1"],
    "start": [100, 200],
    "end": [150, 250]
}).lazy()

# Set coordinate system metadata (requires polars-config-meta)
df = df.config_meta.set(coordinate_system_zero_based=False)  # 1-based

# Now it can be used with range operations
result = pb.overlap(df, other_df, ...)

import pandas as pd

# Create a DataFrame
pdf = pd.DataFrame({
    "chrom": ["chr1", "chr1"],
    "start": [100, 200],
    "end": [150, 250]
})

# Set coordinate system metadata via df.attrs
pdf.attrs["coordinate_system_zero_based"] = False  # 1-based

# Now it can be used with range operations
result = pb.overlap(pdf, other_df, output_type="pandas.DataFrame", ...)

Error handling

polars-bio raises specific errors to prevent coordinate system mismatches:

MissingCoordinateSystemError

Raised when a DataFrame lacks coordinate system metadata:

import polars as pl
import polars_bio as pb

# DataFrame without metadata
df = pl.DataFrame({"chrom": ["chr1"], "start": [100], "end": [200]}).lazy()

# This raises MissingCoordinateSystemError
pb.overlap(df, other_df, ...)

How to fix: Set metadata on your DataFrame before passing it to range operations (see examples above).

CoordinateSystemMismatchError

Raised when two DataFrames have different coordinate systems:

import polars_bio as pb

# One DataFrame is 1-based, another is 0-based
df1 = pb.scan_vcf("file.vcf")                    # 1-based (default)
df2 = pb.scan_bed("file.bed", use_zero_based=True)  # 0-based

# This raises CoordinateSystemMismatchError
pb.overlap(df1, df2, ...)

How to fix: Ensure both DataFrames use the same coordinate system.

Default behavior (lenient validation)

By default, polars-bio uses lenient validation (coordinate_system_check=false). When a DataFrame lacks coordinate system metadata, it falls back to the global configuration and emits a warning:

import polars as pl
import polars_bio as pb

# DataFrames without metadata will use the global config with a warning
df = pl.DataFrame({"chrom": ["chr1"], "start": [100], "end": [200]}).lazy()
result = pb.overlap(df, other_df, ...)  # Uses global coordinate system setting
# Warning: Coordinate system metadata is missing. Using global config...

Strict mode

For production pipelines where coordinate system consistency is critical, you can enable strict validation:

import polars_bio as pb

# Enable strict coordinate system check
pb.set_option(pb.POLARS_BIO_COORDINATE_SYSTEM_CHECK, True)

# Now DataFrames without metadata will raise MissingCoordinateSystemError

Migration from previous versions

If you're upgrading from a previous version of polars-bio:

1. Range operations no longer accept use_zero_based parameter - coordinate system is read from DataFrame metadata
2. I/O functions use use_zero_based parameter (renamed from one_based with inverted logic)
3. Pandas DataFrames require explicit metadata - set df.attrs["coordinate_system_zero_based"] before range operations

# Before (old API)
result = pb.overlap(df1, df2, use_zero_based=True, ...)

# After (new API) - set metadata at I/O time or on DataFrames
df1 = pb.scan_vcf("file.vcf", use_zero_based=True)
df2 = pb.scan_bed("file.bed", use_zero_based=True)
result = pb.overlap(df1, df2, ...)  # Reads from metadata

File metadata

polars-bio automatically attaches comprehensive metadata to DataFrames when reading genomic files. This metadata includes format information, coordinate systems, and format-specific details like VCF header fields.

Metadata structure

The metadata is stored in a clean, user-friendly structure:

import polars_bio as pb

lf = pb.scan_vcf("variants.vcf")
meta = pb.get_metadata(lf)

# Returns:
{
  "format": "vcf",                           # File format
  "path": "variants.vcf",                    # Source file path
  "coordinate_system_zero_based": False,     # Coordinate system (VCF is 1-based)
  "header": {
    "version": "VCFv4.2",                    # VCF version
    "sample_names": ["Sample1", "Sample2"],  # Sample names
    "info_fields": {                         # INFO field definitions
      "AF": {
        "number": "A",
        "type": "Float",
        "description": "Allele Frequency",
        "id": "AF"
      }
    },
    "format_fields": {                       # FORMAT field definitions
      "GT": {
        "number": "1",
        "type": "String",
        "description": "Genotype"
      }
    },
    "contigs": [...],                        # Contig definitions
    "filters": [...],                        # Filter definitions
    "_datafusion_table_name": "variants"     # Internal table name (for debugging)
  }
}

Accessing metadata

polars-bio provides three main functions for working with metadata:

1. Get all metadata as a dictionary

import polars_bio as pb

lf = pb.scan_vcf("file.vcf")
meta = pb.get_metadata(lf)

# Access different parts
print(meta["format"])                       # "vcf"
print(meta["path"])                         # "file.vcf"
print(meta["coordinate_system_zero_based"]) # False (1-based)

# Access VCF-specific fields
print(meta["header"]["version"])            # "VCFv4.2"
print(meta["header"]["sample_names"])       # ["Sample1", "Sample2"]

# Access INFO field definitions
af_field = meta["header"]["info_fields"]["AF"]
print(af_field["type"])                     # "Float"
print(af_field["description"])              # "Allele Frequency"

# Access FORMAT field definitions
gt_field = meta["header"]["format_fields"]["GT"]
print(gt_field["type"])                     # "String"

2. Print metadata as formatted JSON

import polars_bio as pb

lf = pb.scan_vcf("file.vcf")

# Print as pretty JSON
pb.print_metadata_json(lf)

# Customize indentation
pb.print_metadata_json(lf, indent=4)

3. Print human-readable summary

import polars_bio as pb

lf = pb.scan_vcf("file.vcf")
pb.print_metadata_summary(lf)

Output:

======================================================================
Metadata Summary
======================================================================

Format: vcf
Path: file.vcf
Coordinate System: 1-based

Format-specific metadata:
----------------------------------------------------------------------
  VCF Version: VCFv4.2
  Samples (3): Sample1, Sample2, Sample3
  INFO fields: 5
    - AF: Float (Allele Frequency)
    - DP: Integer (Total Depth)
    - AC: Integer (Allele Count)
  FORMAT fields: 3
    - GT: String (Genotype)
    - DP: Integer (Read Depth)
    - GQ: Integer (Genotype Quality)

======================================================================

Format-specific metadata

Different file formats include different metadata:

lf = pb.scan_vcf("variants.vcf")
meta = pb.get_metadata(lf)

# VCF header metadata
meta["header"]["version"]          # VCF version
meta["header"]["sample_names"]     # Sample names
meta["header"]["info_fields"]      # INFO field definitions
meta["header"]["format_fields"]    # FORMAT field definitions
meta["header"]["contigs"]          # Contig definitions
meta["header"]["filters"]          # Filter definitions

lf = pb.scan_fastq("reads.fastq.gz")
meta = pb.get_metadata(lf)

# FASTQ-specific metadata
meta["format"]                     # "fastq"
meta["path"]                       # "reads.fastq.gz"
meta["coordinate_system_zero_based"] # None (N/A for FASTQ)

lf = pb.scan_bed("regions.bed")
meta = pb.get_metadata(lf)

# Basic metadata
meta["format"]                     # "bed"
meta["coordinate_system_zero_based"] # True (0-based)

Setting custom metadata

You can set metadata on DataFrames created from other sources:

import polars as pl
import polars_bio as pb

# Create a DataFrame
df = pl.DataFrame({
    "chrom": ["chr1", "chr1"],
    "start": [100, 200],
    "end": [150, 250]
}).lazy()

# Set metadata
pb.set_source_metadata(
    df,
    format="bed",
    path="custom.bed",
    header={"description": "Custom intervals"}
)

# Now metadata is available
meta = pb.get_metadata(df)
print(meta["format"])  # "bed"
print(meta["header"]["description"])  # "Custom intervals"

Metadata preservation

Metadata is preserved through Polars operations:

lf = pb.scan_vcf("variants.vcf")

# Metadata persists after operations
filtered = lf.filter(pl.col("qual") > 30)
selected = lf.select(["chrom", "start", "end"])
limited = lf.head(100)

# All have the same metadata
meta1 = pb.get_metadata(lf)
meta2 = pb.get_metadata(filtered)
meta3 = pb.get_metadata(selected)

assert meta1["format"] == meta2["format"] == meta3["format"]  # All "vcf"

Using metadata for debugging

The _datafusion_table_name field is useful for debugging DataFusion SQL queries:

lf = pb.scan_vcf("variants.vcf")
meta = pb.get_metadata(lf)

# Get internal table name
table_name = meta["header"]["_datafusion_table_name"]
print(f"Table name: {table_name}")  # "variants"

# Use it in SQL queries for debugging
result = pb.sql(f"SELECT COUNT(*) FROM {table_name}")

API reference