# Bio Core Chromatogram Analysis

> Parse Sanger .ab1/.abi chromatograms with BioPython, extract Phred quality/trace channels, plot traces, quality-trim, and flag het double-peaks. Use for .ab1/.scf files or detecting het SNPs/mixed peaks.

- Skill: `pavel-kravchenko/bio-core-chromatogram-analysis` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/bio-core-chromatogram-analysis`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/bio-core-chromatogram-analysis/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-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.

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

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

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