Applied - Genomic File Formats

Core Principle: Standard bioinformatics file formats (VCF, BED, GFF/GTF, SAM/BAM) encode genomic intervals and annotations under strict coordinate indexing conventions and tabular specifications.

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

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

Exercise

Write the loader the lesson recommends: one function that recognises which of two formats it was handed, converts both to a single coordinate convention, and normalises chromosome names on the way.

Write convert(text), returning a list of output lines. The first input line is a feature type to keep, or * for all of them. Everything after it is either a GFF3 file or a BED file. Ignore blank lines and lines beginning with #. Decide the format from the first data line: nine tab-separated fields with integers in the fourth and fifth is gff3; three to six fields with integers in the second and third is bed; anything else, and a file with no data line at all, is unknown. GFF3 is 1-based closed and BED is 0-based half-open, so convert the GFF3 start by subtracting one and leave its end alone.

The type filter applies only to GFF3, whose third column carries the feature type; BED has no type column, so keep every BED record regardless. Take the GFF3 name from the ID attribute, falling back to Name, then to .; take the BED name from the fourth column and its strand from the sixth, using . when either is absent. Normalise every chromosome name: prefix chr when it is missing, and rewrite a name of chrMT as chrM. Print format=<name>, then filter=<type> -- none when the filter is *, and ignored for a BED input -- then records=<n>, then one line per record holding the chromosome, the 0-based start, the end, the name and the strand, tab separated, sorted by chromosome, then start, then end, then name. When no record survives, print no records instead of the record lines.