Bio core chromatogram analysis

Chromatogram Analysis: Sanger Sequencing

When to Use

• Verifying a plasmid, PCR product, or CRISPR edit from a Sanger .ab1 trace file
• Trimming low-quality primer/dye-blob regions before using a Sanger read downstream
• Visualizing the four-color chromatogram trace alongside base calls
• Screening a trace for heterozygous positions (double peaks) or contamination
• Building a consensus sequence from overlapping forward/reverse Sanger reads

Version Compatibility

• biopython ≥ 1.81 (Bio.SeqIO ABI parser), numpy ≥ 1.24, matplotlib ≥ 3.7, Python ≥ 3.10

Prerequisites

• pip install biopython numpy matplotlib
• Familiarity with Phred quality scores and basic sequence I/O (see bio-sequence-io-read-sequences)

Phred Quality Score Reference

Q Score	Error Probability	Accuracy	Interpretation
10	1/10	90%	Poor — avoid
20	1/100	99%	Minimum acceptable
30	1/1,000	99.9%	High confidence
40	1/10,000	99.99%	Excellent

Formula: Q = -10 × log₁₀(P_error). Practical thresholds: Q<20 = trim/re-sequence; Q20–29 = acceptable; Q≥30 = high confidence.

Sanger read quality profile: bases 1–25 are poor (primer artifact, dye blobs) — always trim; bases 25–700 form a high-quality plateau (Q>30 typical); bases 700+ decay sharply — trim the tail. Max practical read length is ~800 bp.

Goal: load an .ab1 trace and get sequence + per-base quality + raw fluorescence channels. Approach: SeqIO.read(..., "abi") parses the ABIF container into a SeqRecord; quality lives in letter_annotations, raw trace data lives in annotations["abif_raw"] under standard ABIF tags.

from Bio import SeqIO
import numpy as np


def read_ab1(filepath):
    """Read an .ab1 chromatogram and return (SeqRecord, quality array)."""
    record = SeqIO.read(filepath, "abi")
    quals = np.array(record.letter_annotations["phred_quality"])
    print(f"{filepath}: {len(record.seq)} bp, mean Q={quals.mean():.1f}, median Q={np.median(quals):.0f}")
    return record, quals


def extract_traces(record):
    """Pull processed four-color trace channels and peak locations from an .ab1 SeqRecord."""
    raw = record.annotations["abif_raw"]
    base_order = raw.get("FWO_1", b"GATC").decode()  # e.g. "GATC"
    channels = {base: np.array(raw[f"DATA{9 + i}"]) for i, base in enumerate(base_order)}  # DATA9-12 = processed
    peak_locs = list(raw["PLOC1"])  # scan position of each called base
    return channels, peak_locs

Key .ab1 Tags

Tag	Content
DATA9–DATA12	Processed trace channels (baseline-corrected, order from FWO_1)
DATA1–DATA4	Raw trace channels
PLOC1 / PLOC2	Peak locations (scan numbers)
PBAS1 / PBAS2	Called bases (raw / edited)
PCON1 / PCON2	Quality values
FWO_1	Base order for channels (e.g. b'GATC')

Note: .ab1 and .abi are the same format; .scf is an older MRC-LMB format. BioPython reads both via "abi".

Goal: visualize a chromatogram and trim low-quality ends. Approach: plot each color channel over a scan-position window, label called bases at their peak positions; trim by sliding a window of quality scores in from each end until the mean clears the threshold.

import matplotlib.pyplot as plt

TRACE_COLORS = {"A": "green", "C": "blue", "G": "black", "T": "red"}


def plot_chromatogram(channels, peak_locs, sequence, start=0, end=50, title="Chromatogram"):
    """Plot a four-color Sanger trace between base indices [start, end)."""
    end = min(end, len(peak_locs))
    scan_start = max(0, peak_locs[start] - 10)
    scan_end = min(len(channels["A"]), peak_locs[end - 1] + 10)
    x = np.arange(scan_start, scan_end)

    fig, ax = plt.subplots(figsize=(14, 3.5))
    for base, color in TRACE_COLORS.items():
        ax.plot(x, channels[base][scan_start:scan_end], color=color, lw=0.8, label=base)
    for i in range(start, end):
        base = sequence[i]
        ax.text(peak_locs[i], 0, base, ha="center", fontsize=7, fontweight="bold",
                color=TRACE_COLORS.get(base, "gray"))
    ax.legend(ncol=4, fontsize=8)
    ax.set_xlabel("Scan position"); ax.set_ylabel("Fluorescence"); ax.set_title(title)
    plt.tight_layout(); plt.show()


def trim_by_quality(sequence, quals, min_q=20, window=10):
    """Trim low-quality 5'/3' ends using a sliding-window mean quality (Sanger convention)."""
    quals = np.asarray(quals)
    n = len(quals)

    start = 0
    for i in range(n - window + 1):
        if quals[i:i + window].mean() >= min_q:
            start = i
            break

    end = n
    for i in range(n - window, -1, -1):
        if quals[i:i + window].mean() >= min_q:
            end = i + window
            break

    return sequence[start:end], quals[start:end]

Goal: flag heterozygous positions (double peaks) directly from trace intensities. Approach: at each called base's peak scan position, compare the fluorescence of the called (primary) channel to the other three; if a secondary channel's intensity exceeds a ratio threshold of the primary, report it and map the base pair to an IUPAC ambiguity code.

IUPAC_AMBIGUITY = {
    frozenset("AG"): "R", frozenset("CT"): "Y", frozenset("GC"): "S",
    frozenset("AT"): "W", frozenset("GT"): "K", frozenset("AC"): "M",
}


def detect_heterozygous(channels, peak_locs, sequence, ratio_threshold=0.35):
    """Flag positions with a significant secondary peak (candidate het SNP or contamination)."""
    results = []
    for i, (pos, base) in enumerate(zip(peak_locs, sequence.upper())):
        intensities = {b: channels[b][min(pos, len(channels[b]) - 1)] for b in "ATGC"}
        primary = intensities.get(base, 0)
        if primary == 0:
            continue
        for alt_base, alt_val in intensities.items():
            if alt_base == base:
                continue
            ratio = alt_val / primary
            if ratio >= ratio_threshold:
                iupac = IUPAC_AMBIGUITY.get(frozenset([base, alt_base]), "N")
                results.append({"position": i + 1, "called": base, "alt": alt_base,
                                 "ratio": ratio, "iupac": iupac})
    return results

Pitfalls

• Double peaks = heterozygosity OR contamination: two overlapping peaks of comparable height at one position — check surrounding context; genuine het SNPs form a consistent pattern across the read; contamination creates irregular double peaks throughout
• First 20–30 bases always low quality: this is normal (primer signal + dye artifacts) — never report variants from this region
• Why 4 dyes: each ddNTP (ddA, ddC, ddG, ddT) carries a different fluorescent dye; capillary electrophoresis separates by size, laser detects dye — the instrument reconstructs sequence from size+color
• .ab1 stores raw fluorescence, not just sequence: base calls are computed by KB Basecaller software; you can re-call from raw traces with third-party tools
• DATA9–DATA12 are processed; DATA1–DATA4 are raw: use processed channels for visualization (baseline correction and color deconvolution already applied)
• Sanger max read length ~800 bp: beyond this, band resolution collapses — design primers to place your region of interest in the Q>30 window (bases ~30–700)
• ratio_threshold in detect_heterozygous needs tuning per instrument: 0.35 works for typical ABI 3730 traces but noisy runs may need 0.4–0.5 to avoid false positives

See Also

• bio-sequence-io-read-sequences — general sequence file I/O with BioPython
• bio-primer-design-primer-validation — validating primers used to generate the Sanger read
• bio-variant-calling-vcf-basics — downstream representation once a het call is confirmed
• bio-sequence-manipulation-seq-objects — working with the resulting Seq/SeqRecord objects