Applied - Quality Control Plots

In one line: the figures you make before the analysis — because they are what catch the sample swap, the batch effect and the failed library.

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

fig, axd = plt.subplot_mosaic(
    [["qual",  "qual",  "comp"],
     ["libsize", "pca", "pca"],
     ["corr",  "corr",  "corr"]],
    figsize=(13, 12), layout="constrained",
)
# ... populate each panel ...
for k, ax in axd.items():
    ax.text(-0.06, 1.06, k[0].upper(), transform=ax.transAxes,
            fontsize=11, fontweight="bold")
fig.savefig("qc_report.pdf", bbox_inches="tight", dpi=200)

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.

Applied - Quality Control Plots

Applied - Quality Control Plots

Per-base sequence quality

Per-base composition

Library size and detection rate

PCA coloured by every covariate

Sample correlation heatmap

Coverage uniformity

Insert size (paired-end)

A QC dashboard

What each plot catches

Common mistakes

See also