Applied - Quality Control Plots

Core Principle: Diagnostic quality control visualization identifies experimental batch effects, library complexity limitations, and outlier samples through PCA, read-distribution, and depth metrics.

Uses: Seaborn · Faceting · Window Functions · SeqIO · Axes and Reductions · groupby

QC plots are the highest-return figures you will ever make, and the ones most often skipped. Every one below has, at some point, saved someone from publishing a result that was an artefact.

Per-base sequence quality

import gzip, numpy as np
from Bio import SeqIO

def quality_matrix(path, n_reads=100_000, max_len=300):
    q = np.full((n_reads, max_len), np.nan, dtype=np.float32)
    with gzip.open(path, "rt") as fh:
        for i, rec in enumerate(SeqIO.parse(fh, "fastq")):
            if i >= n_reads:
                break
            s = rec.letter_annotations["phred_quality"][:max_len]
            q[i, :len(s)] = s
    return q[:i]

q = quality_matrix("reads.fastq.gz")

pos = np.arange(q.shape[1])
med = np.nanmedian(q, axis=0)                        # axis=0 → per POSITION
q25, q75 = np.nanpercentile(q, [25, 75], axis=0)
q10, q90 = np.nanpercentile(q, [10, 90], axis=0)

fig, ax = plt.subplots(figsize=(8, 3.5))
ax.fill_between(pos, q10, q90, alpha=0.2, color="steelblue", label="10–90%")
ax.fill_between(pos, q25, q75, alpha=0.4, color="steelblue", label="25–75%")
ax.plot(pos, med, c="navy", lw=1.2, label="median")
ax.axhspan(0, 20, color="crimson", alpha=0.08)
ax.axhspan(20, 28, color="orange", alpha=0.08)
ax.set(xlabel="Position in read (bp)", ylabel="Phred quality",
       ylim=(0, 42), title="Per-base quality")
ax.legend(frameon=False, fontsize=8)

axis=0 because reads are rows and positions are columns. This is the FastQC "per base sequence quality" plot, and the expected shape is a gentle decline toward the 3' end. A cliff, or a dip in the middle, indicates a run problem.

Per-base composition

def base_composition(path, n_reads=100_000, max_len=300):
    counts = {b: np.zeros(max_len) for b in "ACGTN"}
    total = np.zeros(max_len)
    with gzip.open(path, "rt") as fh:
        for i, rec in enumerate(SeqIO.parse(fh, "fastq")):
            if i >= n_reads:
                break
            for j, b in enumerate(str(rec.seq)[:max_len]):
                counts.get(b, counts["N"])[j] += 1
                total[j] += 1
    return pd.DataFrame({b: c / np.maximum(total, 1) for b, c in counts.items()})

comp = base_composition("reads.fastq.gz")
comp.plot(figsize=(8, 3), color={"A":"#3DA853","C":"#4285F4",
                                 "G":"#F9AB00","T":"#EA4335","N":"grey"})
plt.ylabel("Fraction"); plt.xlabel("Position in read")

The four lines should be roughly flat and roughly equal. Divergence in the first ~12 bases is normal for RNA-seq: random hexamer priming bias: and is a well-known FastQC "failure" that is not a failure. Divergence elsewhere means adapter contamination or an overrepresented sequence.

Library size and detection rate

qc = pd.DataFrame({
    "sample": samples,
    "lib_size": counts.sum(axis=0),
    "n_detected": (counts > 0).sum(axis=0),
}).merge(meta, on="sample")

fig, axes = plt.subplots(1, 3, figsize=(13, 3.5), layout="constrained")

sns.barplot(data=qc, x="sample", y="lib_size", hue="condition", ax=axes[0])
axes[0].tick_params(axis="x", rotation=90, labelsize=6)
axes[0].set(ylabel="Library size (reads)")

sns.scatterplot(data=qc, x="lib_size", y="n_detected", hue="batch",
                style="condition", s=60, ax=axes[1])
axes[1].set(xlabel="Library size", ylabel="Genes detected")

sns.boxplot(data=log_cpm_long, x="sample", y="log_cpm", hue="condition",
            showfliers=False, ax=axes[2])
axes[2].tick_params(axis="x", rotation=90, labelsize=6)

The middle plot is the informative one: library size and genes detected should rise together and saturate. A sample off that curve: high depth, few genes: is a low-complexity library, usually meaning too little input material or over-amplification.

PCA coloured by every covariate

The single most valuable QC plot.

X = log_cpm[hvg_idx].T
X = X - X.mean(axis=0, keepdims=True)
U, S, Vt = np.linalg.svd(X, full_matrices=False)
ve = S**2 / (S**2).sum()

pcs = pd.DataFrame(U[:, :4] * S[:4], columns=["PC1","PC2","PC3","PC4"])
pcs["sample"] = samples
pcs = pcs.merge(meta, on="sample")

fig, axes = plt.subplots(1, 4, figsize=(16, 3.8), layout="constrained")
for ax, colour in zip(axes, ["condition", "batch", "sex", "lib_size"]):
    sns.scatterplot(data=pcs, x="PC1", y="PC2", hue=colour, s=70, ax=ax)
    ax.set(xlabel=f"PC1 ({ve[0]:.1%})", ylabel=f"PC2 ({ve[1]:.1%})", title=colour)
    ax.set_aspect("equal")

Two non-negotiables:

• Variance explained on the axis labels. PC1 at 45% and PC1 at 4% look identical and mean completely different things.
• set_aspect("equal"). PCA distances are only interpretable when the axes share a scale.

Read it like this: if samples separate by condition, good. If they separate by batch, you have a batch effect. If they separate by lib_size, your normalisation is inadequate. Colour by every covariate you have, including the boring ones: RIN, extraction date, technician, plate position.

# quantify it
for cov in ["condition", "batch", "sex"]:
    for pc in ["PC1", "PC2", "PC3"]:
        groups = [g[pc].values for _, g in pcs.groupby(cov)]
        if len(groups) > 1:
            f, p = stats.f_oneway(*groups)
            if p < 0.05:
                print(f"{cov} associates with {pc}: p={p:.2g}")

Sample correlation heatmap

corr = pd.DataFrame(np.corrcoef(log_cpm, rowvar=False),
                    index=samples, columns=samples)

cg = sns.clustermap(corr, cmap="rocket", vmin=0.7, vmax=1.0,
                    row_colors=col_colors, col_colors=col_colors,
                    figsize=(8, 8), cbar_kws=dict(label="Pearson r"))

Replicates should correlate ~0.95+. A sample correlating better with the other group than its own is a swap. This plot catches sample swaps, and sample swaps are more common than anyone likes to admit.

Coverage uniformity

cov = pd.read_csv("coverage.bed", sep="\t",
                  names=["chrom", "pos", "depth"], dtype={"chrom": "category"})
cov = cov.sort_values(["chrom", "pos"])

# smooth WITHIN chromosome: never across boundaries
cov["smoothed"] = (cov.groupby("chrom", observed=True)["depth"]
                      .transform(lambda s: s.rolling(1000, center=True,
                                                     min_periods=1).mean()))

fig, ax = plt.subplots(figsize=(11, 3))
sub = cov[cov["chrom"] == "chr17"]
ax.fill_between(sub["pos"], sub["smoothed"], alpha=0.6, step="mid")
ax.axhline(sub["smoothed"].median(), ls="--", c="crimson", lw=0.8)
ax.xaxis.set_major_formatter(plt.FuncFormatter(lambda v, _: f"{v/1e6:.1f}"))
ax.set(xlabel="Position on chr17 (Mb)", ylabel="Depth")

# dropout regions via run-length encoding
low = cov["smoothed"] < 0.2 * cov["smoothed"].median()
runs = (low != low.shift()).cumsum()
dropouts = (cov[low].groupby(runs)
            .agg(chrom=("chrom","first"), start=("pos","min"),
                 end=("pos","max"), n=("pos","size"))
            .query("n >= 50"))

The (flag != flag.shift()).cumsum() idiom is the run-length encoding pattern from Window Functions: it is how you turn a boolean mask into contiguous segments.

groupby("chrom").transform(...) is essential. A plain rolling mean averages the end of chr1 with the start of chr2.

Insert size (paired-end)

import pysam
bam = pysam.AlignmentFile("sample.bam", "rb")
sizes = [abs(r.template_length) for r in bam.head(500_000)
         if r.is_proper_pair and 0 < abs(r.template_length) < 1000]

fig, ax = plt.subplots(figsize=(6, 3))
sns.histplot(sizes, bins=100, ax=ax)
ax.set(xlabel="Insert size (bp)", ylabel="Read pairs")
print(f"median {np.median(sizes):.0f}, IQR {np.percentile(sizes,[25,75])}")

A single clean peak is expected. Bimodality, or a peak below the read length, indicates adapter dimers or library prep problems.

A QC dashboard

Generate this automatically for every run and look at it before doing anything else. See Subplots and Layout.

What each plot catches

Plot	Catches
Per-base quality	run failure, cycle problems
Base composition	adapter contamination, priming bias
Library size / detection	low-input, over-amplified libraries
PCA by covariate	batch effects, confounding, outliers
Sample correlation	sample swaps, failed samples
Coverage uniformity	capture dropout, GC bias, CNVs
Insert size	adapter dimers, prep failures
p-value histogram	model misspecification (see Applied - Differential Expression Volcano Plot)

The two bolded rows are the ones that most often change what you do next.

Common mistakes

• Skipping QC and discovering the batch effect after analysis.
• Only colouring PCA by the variable of interest. Colour by everything.
• No variance explained on PCA axes.
• PCA without equal aspect.
• Rolling means across chromosome boundaries.
• Ignoring an outlier because removing it weakens the result. Decide the exclusion rule before looking.
• Not saving QC plots. They belong in the supplement.
• Treating rendered plots as the record. Save the underlying metrics table too.

See also

Applied - Expression Matrices with NumPy · Applied - Heatmaps and Clustermaps · Applied - FASTA and FASTQ Workflows · Window Functions · Faceting · Subplots and Layout · Statistical Estimation in seaborn

Exercise

Build the per-base quality figure the lesson opens with, using matplotlib on the Agg backend and never calling show. seaborn is not installed here, and Biopython is not either, so parse the FASTQ four lines at a time yourself. Assert on what the Axes ended up holding, not on pixels.

Write qc(text), returning a list of output lines. The first input line is the mean-quality cutoff. The records follow, four lines each, with quality in Phred+33. Skip a record whose sequence and quality strings differ in length and count it as malformed. Reads may have different lengths, so the mean at each position must be taken over only the reads that actually reach that position -- averaging over all reads silently drags the tail down.

Plot the per-position means as one line, add a horizontal line at the cutoff, and set the y limits to 0 and 45. Print reads=<n> malformed=<n> max_len=<n>, then first_fail=<p> giving the 1-based position of the first position whose mean is below the cutoff, or none when every position passes, then means=<v,v,...> with each mean rounded to one decimal, then lines=<n> ylim=<low>,<high> where the line count comes from len(ax.get_lines()) and the limits are read back from the Axes and rounded to one decimal. If no usable record is found, print the first line with max_len=0 and then no data, and nothing else.