# Bio Applied Bio Data Formats

> Parse/write FASTA, FASTQ, SAM/BAM, VCF, BED, GFF/GTF with pysam and pure Python; decode SAM FLAG/CIGAR; reconcile 0-based vs 1-based coordinates. Use for custom format parsers or off-by-one coordinate bugs.

- Skill: `pavel-kravchenko/bio-applied-bio-data-formats` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/bio-applied-bio-data-formats`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/bio-applied-bio-data-formats/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: Coding & Dev Tools
- Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/pavel-kravchenko/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

| Format | Content | Coord system | Random access | Tools |
|--------|---------|-------------|--------------|-------|
| FASTA | Sequences | N/A | With `.fai` index | samtools faidx, BioPython |
| FASTQ | Reads + quality | N/A | No (stream) | fastp, cutadapt |
| SAM | Alignments (text) | 1-based | No | samtools, awk |
| BAM | Alignments (binary) | 1-based | With `.bai` | samtools, pysam |
| CRAM | Alignments (ref-compressed) | 1-based | With `.crai` | samtools, htslib |
| VCF | Variants | 1-based | With `.tbi` (tabix) | bcftools, cyvcf2 |
| BED | Intervals | **0-based half-open** | With `.bai`/tabix | bedtools, pybedtools |
| GFF3/GTF | Gene annotations | **1-based inclusive** | With tabix | gffutils |
| BigWig | Coverage/signal | 0-based | Yes (UCSC) | pyBigWig, deeptools |
| PDB/mmCIF | 3D structure | 1-based residues | No | Bio.PDB, MDAnalysis |

### Coordinate trap, side by side
```text
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.

```python
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.

```python
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.

```python
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.

```python
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)