Bio Applied Bio Data Formats

Applied Bioinformatics Data Formats

When to Use

• Writing or debugging a custom parser for FASTA/FASTQ/SAM/VCF/BED instead of reaching for a full pipeline tool.
• Tracking down an off-by-one coordinate bug when converting between BED (0-based) and VCF/GFF/SAM (1-based).
• Decoding a SAM FLAG value or a CIGAR string by hand to understand what an alignment record means.
• Building random-access lookups into a FASTA (.fai) or indexed VCF/BAM without loading the whole file.
• Onboarding onto a new NGS dataset and needing a quick reference for what each format's columns mean.

Version Compatibility

• Python ≥3.10, BioPython ≥1.83, pysam ≥0.22 (htslib ≥1.19)
• Works with SAM spec v1.6, VCF spec v4.2/v4.3, GFF3/GTF2.2

Prerequisites

• pip install biopython pysam (pysam requires htslib; on conda use conda install -c bioconda pysam)
• Familiarity with basic Python file I/O and generators
• Related: bio-applied-ngs-fundamentals for sequencing background, bio-core-biopython-essentials for Bio.SeqIO

Format Quick Reference

Coordinate trap, side by side

BED:  chr1  0   100   → covers bases 1–100 (0-based, half-open)
GFF:  chr1  1   100   → covers bases 1–100 (1-based, inclusive)
VCF:  chr1  925952    → position 925952 (1-based)

Same biological interval, different numbers on disk — this is the single most common source of off-by-one bugs in genomics code.

Core Parsers

Goal: Parse FASTA/FASTQ without pulling in a dependency, and understand what BioPython/pysam do under the hood. Approach: Stream line-by-line with generators so multi-GB files never load fully into memory; gzip-transparent via extension sniffing.

import gzip
import io
from pathlib import Path


def parse_fasta(source):
    """Yield (header, sequence) tuples from a FASTA path or in-memory string."""
    if isinstance(source, (str, Path)) and Path(source).exists():
        opener = gzip.open if str(source).endswith(".gz") else open
        with opener(source, "rt") as fh:
            yield from _fasta_iter(fh)
    else:
        yield from _fasta_iter(io.StringIO(str(source)))


def _fasta_iter(fh):
    header, seq_chunks = None, []
    for line in fh:
        line = line.rstrip("\n")
        if line.startswith(">"):
            if header is not None:
                yield header, "".join(seq_chunks)
            header, seq_chunks = line[1:], []
        elif header is not None:
            seq_chunks.append(line)
    if header is not None:
        yield header, "".join(seq_chunks)


def parse_fastq(source, max_reads=None):
    """Yield (header, sequence, quality_string) from a FASTQ path or string."""
    opener = gzip.open if isinstance(source, (str, Path)) and str(source).endswith(".gz") else None
    fh = opener(source, "rt") if opener else (io.StringIO(source) if "\n" in str(source) else open(source, "rt"))
    try:
        count = 0
        it = iter(fh)
        for line in it:
            header = line.rstrip("\n")[1:]        # drop '@'
            seq = next(it).rstrip("\n")            # must be single line (no wrapping)
            next(it)                                # skip '+' separator
            qual = next(it).rstrip("\n")
            yield header, seq, qual
            count += 1
            if max_reads and count >= max_reads:
                break
    finally:
        fh.close()


def phred(char, offset=33):
    """Convert one FASTQ quality character to a Phred score (Illumina ≥1.8 uses offset 33)."""
    return ord(char) - offset

FASTA Random Access with .fai

Goal: Fetch a genomic region from a multi-GB FASTA without reading the whole file. Approach: Parse the 5-column samtools faidx index (NAME LENGTH OFFSET BASES_PER_LINE BYTES_PER_LINE), compute the byte offset, and seek directly.

import re


def faidx_fetch(fasta_path, fai_path, chrom, start, end):
    """Fetch fasta[chrom][start:end], 0-based half-open, using a .fai index."""
    index = {}
    with open(fai_path) as fh:
        for line in fh:
            name, length, offset, bases_per_line, bytes_per_line = line.split()
            index[name] = (int(length), int(offset), int(bases_per_line), int(bytes_per_line))
    length, offset, bases_per_line, bytes_per_line = index[chrom]
    n_full_lines, remainder = divmod(start, bases_per_line)
    byte_start = offset + n_full_lines * bytes_per_line + remainder
    bases_needed = end - start
    bytes_needed = (bases_needed // bases_per_line + 2) * bytes_per_line
    with open(fasta_path, "rb") as fh:
        fh.seek(byte_start)
        raw = fh.read(bytes_needed).decode()
    return re.sub(r"\s", "", raw)[:bases_needed]

In practice, prefer pysam.FastaFile(path).fetch(chrom, start, end) — this hand-rolled version exists to show what the index actually encodes.

SAM FLAG, CIGAR, and pysam

Goal: Decode alignment metadata (FLAG bits, CIGAR ops) and pull reads from an indexed BAM. Approach: FLAG is a bitmask; CIGAR is a run-length-encoded list of (length, op) pairs where each op either consumes query bases, reference bases, both, or neither.

import pysam

FLAG_BITS = {
    1: "paired", 2: "proper_pair", 4: "unmapped", 8: "mate_unmapped",
    16: "reverse_strand", 32: "mate_reverse_strand", 64: "read1", 128: "read2",
    256: "secondary", 512: "qc_fail", 1024: "duplicate", 2048: "supplementary",
}

CIGAR_CONSUMES = {
    "M": "both", "I": "query", "D": "ref", "N": "ref",
    "S": "query", "H": "neither", "P": "neither", "=": "both", "X": "both",
}


def decode_flag(flag):
    """Return the list of set SAM FLAG bit names, e.g. decode_flag(99) -> ['paired', ...]."""
    return [name for bit, name in FLAG_BITS.items() if flag & bit]


def cigar_ref_span(cigar_str):
    """Reference bases consumed by a CIGAR string (the alignment's genomic span)."""
    import re as _re
    total = 0
    for length, op in _re.findall(r"(\d+)([MIDNSHP=X])", cigar_str):
        if CIGAR_CONSUMES[op] in ("both", "ref"):
            total += int(length)
    return total


def fetch_primary_reads(bam_path, chrom, start, end):
    """Yield primary, mapped, non-duplicate reads overlapping a region (0-based half-open)."""
    with pysam.AlignmentFile(bam_path, "rb") as bam:
        for read in bam.fetch(chrom, start, end):
            if read.is_unmapped or read.is_duplicate or read.is_secondary:
                continue
            yield read

VCF and BED Records

Goal: Parse VCF rows into typed records (SNP vs indel, genotypes) and BED intervals with overlap logic. Approach: Split on tabs, split INFO/FORMAT on ;/:, and keep REF/ALT/POS coupled since VCF indels encode an anchor base.

from dataclasses import dataclass, field


@dataclass
class VCFRecord:
    chrom: str
    pos: int          # 1-based
    id: str
    ref: str
    alt: list
    qual: float | None
    filter: list
    info: dict
    samples: list = field(default_factory=list)

    @property
    def is_snp(self):
        return all(len(a) == 1 and len(self.ref) == 1 for a in self.alt)

    @property
    def is_indel(self):
        return any(len(a) != len(self.ref) for a in self.alt)


def parse_vcf(text):
    """Yield VCFRecord objects from VCF text, skipping meta-information lines."""
    for line in text.splitlines():
        if line.startswith("##") or not line.strip():
            continue
        if line.startswith("#CHROM"):
            continue
        chrom, pos, vid, ref, alt_str, qual_str, filt_str, info_str, *rest = line.split("\t")
        alt = alt_str.split(",")
        qual = float(qual_str) if qual_str != "." else None
        filt = filt_str.split(";") if filt_str != "." else []
        info = dict(item.split("=", 1) if "=" in item else (item, True) for item in info_str.split(";"))
        samples = []
        if rest:
            fmt_keys = rest[0].split(":")
            samples = [dict(zip(fmt_keys, s.split(":"))) for s in rest[1:]]
        yield VCFRecord(chrom, int(pos), vid, ref, alt, qual, filt, info, samples)


@dataclass
class BEDRecord:
    chrom: str
    start: int   # 0-based
    end: int     # exclusive
    name: str | None = None

    def overlaps(self, other):
        """True if two 0-based half-open intervals on the same chrom share any base."""
        return self.chrom == other.chrom and self.start < other.end and other.start < self.end

Pitfalls

• Coordinate systems: BED is 0-based half-open; VCF/GFF/SAM are 1-based inclusive. Mixing them silently produces off-by-one errors — the most common bioinformatics bug.
• FASTQ must not wrap sequence/quality lines: unlike FASTA, each record's sequence and quality string must be exactly one line, or line-based parsers desync.
• BAM requires a sorted + indexed file for .fetch(): samtools sort then samtools index; calling fetch() on an unsorted BAM raises or silently returns nothing.
• CRAM needs the reference genome at decode time: set REF_PATH/REF_CACHE or ensure the @SQ UR: header points to a reachable FASTA.
• VCF indels include an anchor base: REF=GACT ALT=G is a 3-bp deletion, not 4-bp; POS points at the anchor, not the deleted bases.
• Index VCFs with bgzip, not gzip: bgzip file.vcf && tabix -p vcf file.vcf.gz — plain gzip blocks tabix random access.
• GTF/GFF attribute parsing: attributes are key "value" semicolon-delimited; never str.split(";") naively, since quoted values can themselves contain ;.

See Also

• bio-applied-ngs-fundamentals — sequencing background and read structure
• bio-applied-variant-calling-and-snp-analysis — calling and filtering variants into VCF
• bio-core-biopython-essentials — Bio.SeqIO/Bio.AlignIO for production-grade parsing
• bio-applied-advanced-ngs — pipeline-level use of these formats (alignment, coverage, duplicates)