pythonforbio.
[WASM idle]
Applied Workflows06/13

Applied - Genomic File Formats

Starting sandbox…
Beginnerlesson

Applied - Genomic File Formats

In one line: what each format holds, its coordinate convention, and which Python library actually reads it.

The table

Format Holds Coords Read with
FASTA sequence — SeqIO, pyfastx, pysam.FastaFile
FASTQ sequence + quality — SeqIO, pysam.FastxFile, pyfastx
SAM/BAM/CRAM read alignments 1-based text / 0-based binary pysam
VCF/BCF variants 1-based pysam, cyvcf2, Reading and Writing Data
BED intervals 0-based half-open pandas, pyranges, bioframe
GFF3/GTF annotation 1-based closed pandas, gffutils, pyranges
bigWig/bigBed dense signal 0-based pyBigWig, deeptools
GenBank/EMBL sequence + features 1-based → 0-based in memory SeqIO
PDB/mmCIF structure residue numbering Bio.PDB
HDF5/.h5ad matrices + metadata — h5py, anndata
Parquet tables — Reading and Writing Data

The coordinate column is the one that will cost you. See Coordinate Systems.

Biopython does not read alignment files

Bio.SeqIO handles sequence formats. It does not read BAM, CRAM, VCF, BCF, tabix indexes or bigWig. For those, pysam is the answer — it wraps htslib, the same C library samtools and bcftools use.

import pysam

bam = pysam.AlignmentFile("sample.bam", "rb")
for read in bam.fetch("chr17", 7_668_000, 7_688_000):     # 0-based!
    read.query_name, read.reference_start, read.mapping_quality, read.cigarstring

vcf = pysam.VariantFile("cohort.vcf.gz")
for rec in vcf.fetch("chr17", 7_668_000, 7_688_000):
    rec.chrom, rec.pos, rec.ref, rec.alts       # rec.pos is 1-based
    rec.info["AF"], rec.samples["S1"]["GT"]

fa = pysam.FastaFile("hg38.fa")
fa.fetch("chr17", 7_676_153, 7_676_154)          # 0-based half-open

Note the mixed conventions inside pysam itself: fetch() takes 0-based coordinates, but rec.pos on a VCF record returns the 1-based POS. This is deliberate — the region query follows the API convention, the record field follows the file format. It is also a reliable source of off-by-one bugs.

FASTA

>NM_000546.6 Homo sapiens tumor protein p53 (TP53), mRNA
GATGGGATTGGGGTTTTCCCCTCCCATGTGCTCAAGACTGGCGCTAAAAGTTTTGAGCTT
  • Header line starts with >; the ID is up to the first whitespace.
  • Lowercase = soft-masked (repeats). .upper() before composition analysis.
  • N = unknown base.
  • .fai index (from samtools faidx) enables random access.

FASTQ

@read_001 1:N:0:ATCACG
GATTACAGATTACAGATTACA
+
IIIIIIIIIIIIIIIIIIIII

Four lines per record. Quality is ASCII with offset 33 (Phred+33) for anything modern.

Q = ord(char) - 33
P_error = 10 ** (-Q / 10)

Q20 = 1% error, Q30 = 0.1%, Q40 = 0.01%. Legacy Illumina (pre-1.8) used offset 64 — a file read with the wrong offset gives quality scores off by 31, which looks like uniformly terrible data.

VCF

##fileformat=VCFv4.3
##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency">
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT SAMPLE1
chr17 7676154 . G C 99 PASS AF=0.001;DP=44 GT:AD:DP 0/1:22,22:44
  • 1-based.
  • INFO packs many variables into one ;-delimited string.
  • FORMAT defines the :-delimited per-sample subfields.
  • Indels anchor on the preceding base: an insertion is REF=A ALT=AT.
  • Left-aligned by convention — the opposite of HGVS. See hgvs Pitfalls.

See Applied - Reading VCF with pandas.

GFF3 / GTF

chr17  RefSeq  gene  7668402  7687550  .  -  .  ID=gene-TP53;Name=TP53
chr17  RefSeq  CDS   7673535  7673837  .  -  0  Parent=rna-NM_000546.6

Nine tab-separated columns, 1-based closed, # comments. Column 9 packs arbitrary key–value attributes — GFF3 uses key=value;, GTF uses key "value";.

See Applied - GFF and Genomic Intervals.

BED

chr17  7668401  7687550  TP53  0  -

0-based half-open — note it is one less than the GFF start for the same feature. Minimum three columns; up to twelve. The chrom, start, end triple is the lingua franca of interval tools.

Compression and indexing

bgzip -c file.vcf > file.vcf.gz     # bgzip, NOT gzip — block-compressed
tabix -p vcf file.vcf.gz            # creates .tbi for random access
samtools faidx genome.fa            # .fai
samtools index sample.bam           # .bai

bgzip is not gzip. A bgzip file is a valid gzip file, but it is written in independently-decompressible blocks, which is what makes tabix random access possible. A plain gzip'd VCF cannot be tabix-indexed, and this catches people out constantly.

In Python, pd.read_csv and gzip.open read bgzip files fine (they are gzip-compatible); you just cannot query a region without going through pysam.

Chromosome naming

The single most common cause of a join returning nothing:

Source Style
UCSC chr1, chrM
Ensembl 1, MT
NCBI RefSeq NC_000001.11 
df["chrom"] = df["chrom"].str.removeprefix("chr")        # → Ensembl style
df["chrom"] = "chr" + df["chrom"].astype(str)            # → UCSC style
df["chrom"] = df["chrom"].replace({"MT": "chrM", "M": "chrM"})

Normalise at load time, in one place, and assert afterwards. See Merging and Joining.

A loader module worth writing

from pathlib import Path
import pandas as pd

GFF_COLS = ["seqid","source","type","start","end","score","strand","phase","attributes"]
BED_COLS = ["chrom","start","end","name","score","strand"]

def read_gff(path):
    df = pd.read_csv(path, sep="\t", comment="#", names=GFF_COLS,
                     dtype={"seqid":"category","type":"category",
                            "start":"int32","end":"int32"})
    df["start0"] = df["start"] - 1        # ← convert ONCE, name it explicitly
    df["end0"] = df["end"]
    return df

def read_bed(path, n_cols=6):
    return pd.read_csv(path, sep="\t", header=None, names=BED_COLS[:n_cols],
                       dtype={"chrom":"category","start":"int32","end":"int32"})

The start0/end0 naming convention is the discipline that prevents coordinate bugs. A bare start column in a frame that mixes BED and GFF data is an accident waiting to happen.

Reproducibility

Record for every reference file: the assembly (GRCh38, and which patch), the annotation release (RefSeq 110, GENCODE v45), the source URL, and a checksum. "We used hg38" is not enough — patch releases add and remove contigs.

See Applied - Reproducible Environment.

Common mistakes

  • Coordinate convention mixups. BED vs GFF vs VCF.
  • chr1 vs 1 producing silent empty joins.
  • gzip instead of bgzip, so tabix indexing fails.
  • Wrong FASTQ quality offset.
  • Ignoring soft-masking.
  • Expecting Biopython to read BAM/VCF.
  • Assuming one row per variant in a VCF. Multi-allelic sites pack several ALTs into one row.
  • Not recording the annotation release.

See also

Coordinate Systems · Applied - Reading VCF with pandas · Applied - GFF and Genomic Intervals · Applied - FASTA and FASTQ Workflows · SeqIO · Reading and Writing Data · Applied - Reproducible Environment

lesson example

No output yet — run the code to populate this drawer.

Applied - GFF and Genomic Intervals←Applied - Heatmaps and Clustermaps→