Read and write the core sequence file formats

Across the whole spine you learned what DNA, reads, and variants actually are. Now you cross into the track where you touch the data, and here is the first surprise: in day-to-day bioinformatics you almost never see a cell. You see files. A small handful of plain-text formats carries nearly all of genomics, and each one is a thin, honest encoding of a specific biological fact. Learn to read the file and you can see the biology straight through it. This lesson teaches the formats you will meet every day, and, more importantly, what each one physically means.

One pipeline, one analogy

The formats are not a random pile. They line up as a pipeline, and it maps cleanly onto something you already know: compile and diff.

• A sequencer produces raw, noisy output, one line per fragment it read, each with a confidence value. That is FASTQ. Think of it as an instrument's log stream, closer to sensor output than to source code.
• Somewhere there is a canonical string everyone agrees to measure against. That is FASTA, the plain sequence. Think of it as the committed source everyone builds against.
• You place your noisy reads back onto that canonical string, recording where each one landed. That is SAM (and its compressed binary twin BAM), like a debugger mapping raw addresses back to source lines through a sourcemap.
• You summarize how one genome differs from the reference. That is VCF, a diff against a specific base commit.
• You pin labels onto coordinates (this range is a gene, that range is an exon). Those are GFF/GTF and BED, like line-range annotations layered over source.

Use the analogy, then respect its failure edge, because an analogy without its limit is a latent bug. In software the coordinate system is authoritative and lossless: apply a diff to the exact base and you reconstruct the exact file, and a byte is a byte. None of that holds here. The reference is a consensus ruler, not ground truth (recall S6.3, no real person matches it everywhere). The reads are physical measurements with per-base uncertainty, not clean bytes. Placement is a probabilistic best guess, not a deterministic lookup, because repetitive genome regions let one read map equally well to several spots. And the same coordinate means different physical bases in different reference builds. Hold the pipeline picture, but never forget you are diffing noisy measurements against a shifting ruler.

The raw measurement: a read, and FASTQ

A read is one physical measurement. The machine shears your DNA sample into fragments, then reads the base sequence off each fragment, typically 100 to 300 bases for short-read Illumina. One fragment gives one read. It is a real observation of one molecule, and like any measurement it carries error.

FASTQ stores exactly that: the sequence plus a per-base confidence. Every record is four lines.

@SEQ_ID:run7:lane1:1234 read 1 of a pair
ACGTTAGCCGATTACAGGCA
+
IIIIFFFF????5555++##

Line one starts with @ and names the read. Line two is the 20 bases the machine called. Line three is a lone + (a separator, historically it could repeat the name). Line four is the quality string, one character per base, aligned position for position with line two.

Each quality character encodes a Phred quality score, and this is the beat that trips up newcomers. The score is defined as Q equals minus ten times log base ten of the error probability. Invert it and the error probability is ten to the power of minus Q over ten. The character is the score plus 33, mapped to ASCII (the Phred+33 offset). So decode by taking the character code, subtracting 33, then reading the log scale.

def phred_to_error(qchar):
    q = ord(qchar) - 33
    return q, 10 ** (-q / 10)

for c in "I?5+#":
    q, p = phred_to_error(c)
    print(c, "Q" + str(q), round(p, 4))

Read this panel, do not expect a Run button yet (in-browser execution is a later phase). Reason it out. The character I is code 73, so Q40, an error probability of 0.0001, one wrong base in ten thousand. ? is Q30, one in a thousand. 5 is Q20, one in a hundred. + is Q10, one in ten. # is Q2, an error near 0.63, so that base is more likely wrong than right. Our example read is trustworthy at the front and collapses at the tail, which is the ordinary Illumina pattern: the chemistry degrades as the read goes on, so the last bases are the least certain.

The canonical string: FASTA

FASTA is the plainest format in biology: a header line beginning with >, then the sequence. No quality, no coordinates, just letters. It holds a reference genome, a single gene, a set of protein sequences, anything sequence-shaped.

>chr17 Homo sapiens chromosome 17, GRCh38 reference
ACGTTAGCCGATTACAGGCATTGACCGGTTAACGGGCATA

One file can hold many records, each its own > header followed by its sequence. Real reference sequences wrap the bases at 60 or 80 characters per line, and real chromosome 17 is about 83 million bases, so this snippet is a keyhole view. The key contrast: FASTA carries the sequence with no confidence attached, while FASTQ carries the sequence plus the confidence. When you have a fixed truth to measure against, you want FASTA. When you have noisy observations, you want FASTQ.

Placing reads on the map: SAM/BAM and coverage

A pile of reads is useless until you know where each one belongs. Alignment takes every read and finds its best position on the reference, recording the result in SAM (Sequence Alignment/Map). BAM is the same content in a compressed binary form, which is what you actually store and index because SAM text is enormous.

QNAME  FLAG  RNAME  POS       MAPQ  CIGAR  SEQ                   QUAL
read1  99    chr17  43044300  60    20M    ACGTTAGCCGATTACAGGCA  IIIIFFFF????5555++##

Walk the columns that matter (a real SAM line has 11 mandatory columns, and this panel shows a representative subset, eliding the mate-pair fields RNEXT, PNEXT, and TLEN for clarity). RNAME and POS say the read landed on chr17 at 1-based position 43044300. CIGAR 20M says all 20 bases aligned in a row (M means aligned, later you will meet I and D for insertions and deletions, the indels from S6.3, and S for soft-clipped ends). MAPQ is 60, and it is Phred-scaled again: it is the confidence in the placement, not the bases, so MAPQ 60 means the aligner is very sure this is the right spot. FLAG 99 is a bitfield packing several yes/no facts, including whether the read aligned to the forward or reverse strand. That last point reaches back to S3.2. Because DNA has two strands and fragments come off both, roughly half your reads are the reverse complement of the reference, so the aligner tries both orientations and the FLAG records which one won.

Once many reads are placed, two physical quantities fall out. Depth at a position is simply how many reads overlap it. Coverage usually means the average depth across a region, so "30x coverage" means each base was read about 30 times on average. Each read is an independent sample of one molecule, so depth is literally how many independent measurements you have of that spot, and that is why it governs how much you can trust a call there.

Below, treat the top strand as a slice of the reference and the bottom as a single read, and align them. Every mismatch is a candidate SNP, every gap is a candidate indel. Predict where they land before you run it.

matrix[i][j] = max(
  matrix[i-1][j-1] + s(a[i-1], b[j-1]),   // diagonal: align the two residues
  matrix[i-1][j]   + gap,                 // up: a gap in sequence b
  matrix[i][j-1]   + gap                  // left: a gap in sequence a
)
s(x, y) = match when x == y, else mismatch

// A thin sketch of needlemanWunsch (global). Local alignment is the same
// loop with a 0 floor and a traceback that starts at the largest cell.
function fill(a: string, b: string, match: number, mismatch: number, gap: number) {
  const m = a.length;
  const n = b.length;
  const matrix: number[][] = [];
  for (let i = 0; i <= m; i++) matrix.push(new Array(n + 1).fill(0));

  // Seed the borders: a run of k gaps costs k * gap. (Local leaves these at 0.)
  for (let i = 1; i <= m; i++) matrix[i][0] = i * gap;
  for (let j = 1; j <= n; j++) matrix[0][j] = j * gap;

  // Fill. Row i indexes a, column j indexes b.
  for (let i = 1; i <= m; i++) {
    for (let j = 1; j <= n; j++) {
      const s = a[i - 1] === b[j - 1] ? match : mismatch;
      const diag = matrix[i - 1][j - 1] + s;   // align a[i-1] with b[j-1]
      const up = matrix[i - 1][j] + gap;        // a gap in b
      const left = matrix[i][j - 1] + gap;      // a gap in a
      matrix[i][j] = Math.max(diag, up, left);
      // Smith-Waterman (local) instead: Math.max(0, diag, up, left)
    }
  }

  // The global score is the bottom-right corner.
  return matrix[m][n];
}

The diff: VCF

You could store a person's genome as three billion letters. Nobody does. You store it as a VCF (Variant Call Format), a list of differences from the reference, exactly the diff-against-a-base-commit idea from S6.3 made into a real file.

##fileformat=VCFv4.2
##reference=GRCh38
#CHROM  POS       ID          REF  ALT  QUAL  FILTER  INFO   FORMAT    SAMPLE1
17      43044295  rs28897696  G    A    60    PASS    DP=54  GT:DP:GQ  0/1:54:99

The data row reads: at chromosome 17, position 43044295, the reference base is G and this sample carries an A. QUAL 60 is Phred-scaled confidence in the variant. The INFO column packs key-and-value tags (here just DP=54, the depth) and in real files it holds several tags joined by a separator (this course's own style forbids that separator character, so picture it between the tags). FORMAT names per-sample fields and SAMPLE1 fills them in: GT is 0/1, meaning one reference allele and one alt allele, a heterozygous variant. Now depth pays off. This site has DP 54, plenty of independent reads, so seeing about half of them carry the A is convincing. At depth 4 it would be a guess: by pure sampling you might see zero alt reads and miss a real heterozygous variant, or all four and wrongly call it homozygous. That is why whole-genome sequencing commonly targets around 30x.

The header carries the fact that makes or breaks the whole file. ##reference=GRCh38 says which build the positions are counted against.

Annotations on the coordinate system: GFF/GTF and BED

The last group does not store sequence at all. It pins labels onto coordinate ranges: this stretch is a gene, that stretch is an exon, this one is a regulatory region. GFF (and its close relative GTF) is nine tab-separated columns. BED is a leaner interval format, as few as three columns.

# GFF3 (1-based, inclusive)
chr17  RefSeq  gene  43044295  43125483  .  -  .  ID=BRCA1

# BED (0-based, half-open)
chr17  43044294  43125483  BRCA1  0  -

Both lines describe the same gene, BRCA1, which sits on chromosome 17's minus strand across roughly 43,044,295 to 43,125,483 in GRCh38. GFF columns are sequence name, source, feature type, start, end, score, strand, phase, and attributes (the ID=BRCA1 tag). BED columns are chrom, start, end, and then optional name, score, strand. Look hard at the two start numbers. GFF says 43044295 and BED says 43044294 for the same first base. That is not a typo.