---
title: "08 · VCF Mutation Analysis"
subtitle: "Mutation frequency, Shannon entropy heatmaps, and tissue-type comparisons"
description: |
Per-sample VCF files from a variant caller. Coverage filtering, Shannon entropy per
position, mutation frequency bar charts and entropy heatmaps per tissue type,
Kruskal-Wallis comparison across tissue sources.
categories: [genomics, virology, VCF, mutation, entropy]
jupyter: python3
---
## Overview
| Item | Details |
|------|---------|
| **Input** | Directory of `.vcf` files: `aa_change_{ID}-{Source}_{Segment}_{Threshold}.vcf` |
| **Key packages** | `vcfpy`, `pandas`, `numpy`, `scipy`, `seaborn`, `matplotlib` |
| **Statistics** | Shannon entropy (natural log) · Kruskal–Wallis with Benjamini–Hochberg FDR |
| **Output** | Mutation bar charts · entropy heatmaps (SVG) · frequency tables · KW results CSV |
| **Download** | [template.ipynb](template.ipynb) |
Edit `SEGMENT_CONFIG` to match your pathogen's genome (total length, CDS boundaries, AA count). `TISSUE_TYPES` controls which sources are included in statistical comparisons.
`template.ipynb` is generated from this file with `quarto convert template.qmd` — edit the `.qmd` and re-convert rather than editing the notebook, so the two cannot drift apart.
[← Gallery](../../index.html){.btn .btn-outline-secondary}
[Download .ipynb](template.ipynb){.btn .btn-primary}
---
## User Configuration
```{python}
#| label: user-config
# ── USER CONFIGURATION ────────────────────────────────────────────────────────
# Directory containing VCF files (relative to this file)
VCF_DIR = "data/aa_change_vcf"
# VCF filename format (regex). Named groups: id, source, segment, threshold
# Default: aa_change_{id}-{source}_{segment}_{threshold}.vcf
VCF_FILENAME_PATTERN = r"aa_change_(?P<id>[^-]+)-(?P<source>[^_]+)_(?P<segment>[SM])_(?P<threshold>\d+pct)\.vcf"
# Map threshold labels in filenames to display strings
THRESHOLD_LABEL_MAP = {"5pct": "5%", "2pct": "2%", "1pct": "1%"}
# Minimum read depth to retain a variant (coverage QC)
DEPTH_THRESHOLD = 500
# Tissue / source types to INCLUDE in the final analysis
# (all others present in VCF filenames will be loaded but filtered out)
TISSUE_TYPES = ["Lung", "Saliva", "Urine"]
# Segment geometry — edit to match your pathogen's genome.
# The CDS must be codon-aligned: (cds_end - cds_start + 1) == aa_length * 3.
# That is checked below, because a CDS that is not a whole number of codons
# means one of the three numbers is wrong.
SEGMENT_CONFIG = {
"S": {
"length": 1902, # total nucleotide length
"cds_start": 250, # first CDS position (1-based)
"cds_end": 1329, # last CDS position (1-based, inclusive)
"aa_length": 360, # total amino acids in CDS
},
"M": {
"length": 3675,
"cds_start": 52,
"cds_end": 3456,
"aa_length": 1135,
},
}
# VAF thresholds to display in bar charts (ordered low → high)
VAF_THRESHOLDS = ["1%", "2%", "5%"]
# Colors and transparency for each VAF threshold in bar charts
VAF_COLORS = {
"5%": "#3594cc", # Blue
"2%": "#ea801c", # Orange
"1%": "#a00000", # Dark red
}
VAF_ALPHA = {"5%": 0.3, "2%": 0.6, "1%": 0.9}
# Family-wise error rate for the Benjamini–Hochberg correction on the
# Kruskal–Wallis tests
FDR_ALPHA = 0.05
# Output directory for SVG plots and derived tables (relative to this file)
OUTPUT_DIR = "Plots"
# ─────────────────────────────────────────────────────────────────────────────
```
## Setup
```{python}
#| label: setup
import os
import re
import warnings
import numpy as np
import pandas as pd
import vcfpy
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.stats import entropy, kruskal, false_discovery_control
# Narrow suppression only: seaborn's violinplot still passes `palette` without
# `hue` internally on some versions. Everything else stays visible on purpose.
warnings.filterwarnings("ignore", category=FutureWarning, module="seaborn")
# ── CDS geometry sanity check ────────────────────────────────────────────────
# Mutations_per_100aa and the coding-region filter both assume the CDS is a whole
# number of codons; a mismatch here means the per-100aa rate is silently wrong.
for _seg, _cfg in SEGMENT_CONFIG.items():
_cds_len = _cfg["cds_end"] - _cfg["cds_start"] + 1
assert _cds_len == _cfg["aa_length"] * 3, (
f"Segment {_seg}: CDS spans {_cds_len} nt, but aa_length="
f"{_cfg['aa_length']} implies {_cfg['aa_length'] * 3} nt. "
f"Fix cds_start / cds_end / aa_length in SEGMENT_CONFIG."
)
assert _cfg["cds_end"] <= _cfg["length"], (
f"Segment {_seg}: cds_end ({_cfg['cds_end']}) exceeds the segment "
f"length ({_cfg['length']})."
)
# Cached so the mutation-rate denominators cannot drift from the config.
_cfg["cds_length"] = _cds_len
# ── Global plot style ─────────────────────────────────────────────────────────
# Reproducible SVG output: matplotlib stamps a wall-clock date and randomly
# generated element IDs into every SVG, so re-rendering an unchanged figure
# produced a large but meaningless git diff. A fixed salt and a fixed
# SOURCE_DATE_EPOCH make repeat renders byte-identical.
os.environ.setdefault("SOURCE_DATE_EPOCH", "1735689600") # 2025-01-01 UTC
plt.rcParams["svg.hashsalt"] = "lab-bioinfo-templates"
plt.rcParams.update({
'axes.titlesize': 12,
'axes.labelsize': 10,
'xtick.labelsize': 10,
'ytick.labelsize': 10,
'legend.fontsize': 10,
})
pd.set_option('display.width', 1000)
os.makedirs(OUTPUT_DIR, exist_ok=True)
print(f"Output will be saved to: {os.path.abspath(OUTPUT_DIR)}")
```
## Load & Parse VCF Files
```{python}
#| label: load-vcf
pattern = re.compile(VCF_FILENAME_PATTERN)
records = []
vcf_files = [f for f in os.listdir(VCF_DIR) if f.endswith(".vcf")]
print(f"Found {len(vcf_files)} VCF files in {VCF_DIR}")
for fname in vcf_files:
m = pattern.match(fname)
if not m:
print(f" Skipping (name mismatch): {fname}")
continue
sample_id = m.group("id")
source = m.group("source")
segment = m.group("segment")
thresh_raw = m.group("threshold")
threshold = THRESHOLD_LABEL_MAP.get(thresh_raw, thresh_raw)
fpath = os.path.join(VCF_DIR, fname)
try:
reader = vcfpy.Reader.from_path(fpath)
for rec in reader:
if not rec.ALT:
continue
alt_allele = rec.ALT[0].value if rec.ALT else "."
aa_change = rec.INFO.get("AA_CHANGE", "N/A")
if isinstance(aa_change, list):
aa_change = aa_change[0]
sample_call = (
rec.call_for_sample["SAMPLE"]
if "SAMPLE" in rec.call_for_sample
else rec.calls[0]
)
ad = sample_call.data.get("AD", [0, 0])
dp = sample_call.data.get("DP", 0)
records.append({
"Position": rec.POS,
"Reference_Allele": rec.REF,
"Alternate_Allele": alt_allele,
"Quality": rec.QUAL,
"Amino_Acid_Change": aa_change,
"Allelic_Depth_Ref": ad[0] if len(ad) > 0 else 0,
"Allelic_Depth_Alt": ad[1] if len(ad) > 1 else 0,
"Total_Depth": dp,
"Source": source,
"Segment": segment,
"Variant_Threshold": threshold,
"ID": sample_id,
})
except Exception as e:
print(f" Error reading {fname}: {type(e).__name__}: {e}")
combined_df = pd.DataFrame(records)
print(f"Total records parsed: {len(combined_df)}")
if combined_df.empty:
raise RuntimeError(
f"No VCF records were parsed from '{os.path.abspath(VCF_DIR)}'. "
f"Check that simulate_data.py ran successfully and that files match "
f"VCF_FILENAME_PATTERN. Files found: {vcf_files[:5]}"
)
REQUIRED_COLUMNS = [
"Position", "Amino_Acid_Change", "Allelic_Depth_Ref", "Allelic_Depth_Alt",
"Total_Depth", "Source", "Segment", "Variant_Threshold", "ID",
]
_missing = set(REQUIRED_COLUMNS) - set(combined_df.columns)
if _missing:
raise ValueError(f"Parsed variant table is missing column(s): {sorted(_missing)}")
unknown_segments = sorted(set(combined_df["Segment"]) - set(SEGMENT_CONFIG))
if unknown_segments:
raise ValueError(
f"VCF filenames reference segment(s) absent from SEGMENT_CONFIG: "
f"{unknown_segments}. Add them, or narrow VCF_FILENAME_PATTERN."
)
combined_df.head()
```
## Quality Control: Coverage Filtering
```{python}
#| label: qc-filter
pre_filter = len(combined_df)
combined_df = combined_df[combined_df["Total_Depth"] >= DEPTH_THRESHOLD].copy()
print(f"Retained {len(combined_df)} / {pre_filter} records (depth ≥ {DEPTH_THRESHOLD}x)")
# Fix Allelic_Depth_Ref when reported as 0 (some callers omit it)
mask = combined_df["Allelic_Depth_Ref"] == 0
combined_df.loc[mask, "Allelic_Depth_Ref"] = (
combined_df.loc[mask, "Total_Depth"] - combined_df.loc[mask, "Allelic_Depth_Alt"]
)
combined_df["VAF"] = combined_df["Allelic_Depth_Alt"] / combined_df["Total_Depth"]
# Save the combined dataset. This is a derived table, so it belongs beside the
# other outputs — data/ holds only the inputs the template reads.
mutation_table_path = os.path.join(OUTPUT_DIR, "mutation_df.csv")
combined_df.to_csv(mutation_table_path, index=False)
print(f"Saved: {mutation_table_path}")
combined_df.describe()
```
## Mutation Distribution Bar Charts
Each bar spans one variant call position; transparency encodes VAF threshold.
```{python}
#| label: mutation-bars
#| fig-width: 14
#| fig-height: 3
# Bar widths per segment (in nucleotides)
bar_widths = {"S": 9.0, "M": 14.0}
for source in TISSUE_TYPES:
src_df = combined_df[combined_df["Source"] == source]
if src_df.empty:
continue
fig, axes = plt.subplots(1, 2, figsize=(14, 3), constrained_layout=True)
for ax, seg in zip(axes, ["S", "M"]):
cfg = SEGMENT_CONFIG[seg]
seg_df = src_df[src_df["Segment"] == seg]
for thresh in VAF_THRESHOLDS:
t_df = seg_df[seg_df["Variant_Threshold"] == thresh]
if t_df.empty:
continue
ax.bar(
t_df["Position"],
t_df["VAF"],
width=bar_widths[seg],
color=VAF_COLORS[thresh],
alpha=VAF_ALPHA[thresh],
label=thresh,
)
ax.set_xlim(1, cfg["length"])
# No fixed y ceiling: a hardcoded limit clips any variant above it off
# the chart without saying so. Let the data set the top.
ax.set_ylim(0, max(0.05, src_df["VAF"].max() * 1.1))
ax.set_xlabel("Genomic position (nt)")
ax.set_ylabel("Variant allele frequency (VAF)")
ax.set_title(f"{seg} segment — {source}")
ax.axvspan(cfg["cds_start"], cfg["cds_end"], alpha=0.06, color="grey", label="CDS")
ax.legend(loc="upper right", fontsize=8)
fig.savefig(os.path.join(OUTPUT_DIR, f"mutation_map_{source}.svg"),
format="svg", bbox_inches="tight")
plt.show()
print(f"Saved: mutation_map_{source}.svg")
```
## Shannon Entropy
Entropy (natural log base) calculated from reference and alternate allele frequencies.
```{python}
#| label: entropy-calc
# Allele frequencies
combined_df["Freq_Ref"] = combined_df["Allelic_Depth_Ref"] / combined_df["Total_Depth"]
combined_df["Freq_Alt"] = combined_df["Allelic_Depth_Alt"] / combined_df["Total_Depth"]
def calculate_entropy(row):
freqs = np.nan_to_num([row["Freq_Ref"], row["Freq_Alt"]], nan=0.0)
return entropy(freqs, base=np.e) # Shannon entropy, natural log
combined_df["Entropy"] = combined_df.apply(calculate_entropy, axis=1)
print(f"Entropy range: {combined_df['Entropy'].min():.4f} – {combined_df['Entropy'].max():.4f}")
combined_df[["Position", "Segment", "Source", "Variant_Threshold", "Entropy"]].head(10)
```
## Entropy Heatmaps
Rows = tissue source, columns = genomic position. Colour = mean Shannon entropy.
```{python}
#| label: entropy-heatmaps
custom_cmap = sns.color_palette("Blues", as_cmap=True)
for seg in SEGMENT_CONFIG:
seg_df = combined_df[
(combined_df["Segment"] == seg) &
(combined_df["Source"].isin(TISSUE_TYPES))
]
# Pivot: mean entropy per source × position
pivot = (
seg_df
.groupby(["Source", "Position"])["Entropy"]
.mean()
.unstack(fill_value=0)
)
if pivot.empty:
print(f"No data for {seg} segment — skipping heatmap")
continue
g = sns.clustermap(
pivot,
cmap=custom_cmap,
row_cluster=True,
col_cluster=False,
figsize=(12, max(2, len(pivot) * 0.8 + 1)),
linewidths=0,
cbar_pos=(1.05, 0.2, 0.03, 0.6),
xticklabels=False,
)
g.ax_heatmap.set_xlabel("Genomic position (nt)")
g.ax_heatmap.set_ylabel("Tissue source")
g.fig.suptitle(f"{seg} segment — Shannon Entropy by Source", y=1.02, fontsize=12)
out_svg = os.path.join(OUTPUT_DIR, f"entropy_map_{seg}.svg")
g.fig.savefig(out_svg, format="svg", bbox_inches="tight")
plt.show()
print(f"Saved: {out_svg}")
```
## Mutation Frequency per Sample
Normalised to mutations per 1,000 nt of **CDS** (nucleotide level) and per 100 amino
acids (non-synonymous only). Both numerators are restricted to coding positions, so
both denominators are coding lengths.
```{python}
#| label: mutation-freq
def filter_coding(df):
"""Retain only CDS positions for each segment."""
masks = []
for seg, cfg in SEGMENT_CONFIG.items():
m = (
(df["Segment"] == seg) &
(df["Position"].between(cfg["cds_start"], cfg["cds_end"]))
)
masks.append(m)
combined_mask = masks[0]
for mask in masks[1:]:
combined_mask = combined_mask | mask
return df[combined_mask]
coding_df = filter_coding(combined_df)
analysis_df = coding_df[coding_df["Source"].isin(TISSUE_TYPES)].copy()
freq_rows = []
for (sample_id, source, seg, thresh), grp in analysis_df.groupby(
["ID", "Source", "Segment", "Variant_Threshold"]
):
cfg = SEGMENT_CONFIG[seg]
n_total = len(grp)
n_nonsyn = (grp["Amino_Acid_Change"] != "N/A").sum()
freq_rows.append({
"ID": sample_id,
"Source": source,
"Segment": seg,
"Threshold": thresh,
# Numerator counts CDS variants only (filter_coding above), so the
# denominator must be the CDS length — dividing by the full segment
# length understates every rate, and by a different factor per segment.
"Mutations_per_1000bp": (n_total / cfg["cds_length"]) * 1000,
"Mutations_per_100aa": (n_nonsyn / cfg["aa_length"]) * 100,
})
freq_df = pd.DataFrame(freq_rows)
print(freq_df.groupby(["Segment", "Threshold"])[
["Mutations_per_1000bp", "Mutations_per_100aa"]
].mean().round(2))
freq_df.head()
```
## Mutation Frequency Violin Plots
```{python}
#| label: mutation-violin
#| fig-width: 10
#| fig-height: 4
for seg in SEGMENT_CONFIG:
seg_freq = freq_df[freq_df["Segment"] == seg]
if seg_freq.empty:
continue
fig, axes = plt.subplots(1, 2, figsize=(10, 4), constrained_layout=True)
fig.suptitle(f"{seg} segment — Mutation Frequency by Source", fontsize=12)
metrics = [
("Mutations_per_1000bp", "Mutations per 1,000 bp"),
("Mutations_per_100aa", "Non-syn. mutations per 100 aa"),
]
for ax, (col, label) in zip(axes, metrics):
sns.violinplot(
data=seg_freq[seg_freq["Source"].isin(TISSUE_TYPES)],
x="Source", y=col, hue="Threshold",
palette=VAF_COLORS, order=TISSUE_TYPES,
inner="quartile", ax=ax,
)
ax.set_xlabel("Tissue type")
ax.set_ylabel(label)
ax.set_title(label)
fig.savefig(os.path.join(OUTPUT_DIR, f"mutation_freq_{seg}.svg"),
format="svg", bbox_inches="tight")
plt.show()
```
## Kruskal–Wallis Tests
Non-parametric comparison of mutation frequencies across tissue types. Up to twelve
tests are run (three thresholds × two segments × two metrics), so the raw p-values
are corrected together with Benjamini–Hochberg
(`scipy.stats.false_discovery_control`). Significance is judged on the corrected
FDR, not on the raw p-value.
```{python}
#| label: kruskal-wallis
kruskal_rows = []
metrics = [
("Mutations_per_1000bp", "Nucleotide (per 1,000 bp)"),
("Mutations_per_100aa", "Amino acid (per 100 aa)"),
]
for thresh in VAF_THRESHOLDS:
t_df = freq_df[freq_df["Threshold"] == thresh]
for seg in SEGMENT_CONFIG:
s_df = t_df[t_df["Segment"] == seg]
for col, metric_label in metrics:
groups = [s_df[s_df["Source"] == src][col].dropna()
for src in TISSUE_TYPES]
groups = [g for g in groups if len(g) >= 2]
if len(groups) < 2:
continue
stat, pval = kruskal(*groups)
kruskal_rows.append({
"Threshold": thresh,
"Segment": seg,
"Metric": metric_label,
"Statistic": round(stat, 4),
"p-value": pval,
})
kruskal_df = pd.DataFrame(kruskal_rows)
# ── Benjamini–Hochberg across the whole family of tests run above ────────────
if len(kruskal_df):
kruskal_df["BH_FDR"] = false_discovery_control(
kruskal_df["p-value"].to_numpy(), method="bh"
)
kruskal_df["Significant"] = kruskal_df["BH_FDR"] < FDR_ALPHA
kruskal_df["p-value"] = kruskal_df["p-value"].round(4)
kruskal_df["BH_FDR"] = kruskal_df["BH_FDR"].round(4)
print(f"{len(kruskal_df)} Kruskal–Wallis tests, "
f"{int(kruskal_df['Significant'].sum())} significant at BH FDR < {FDR_ALPHA} "
f"({int((kruskal_df['p-value'] < 0.05).sum())} would pass an uncorrected "
f"p < 0.05).")
else:
print("No Kruskal–Wallis tests could be run (fewer than two groups with n ≥ 2).")
kruskal_df.to_csv(os.path.join(OUTPUT_DIR, "kruskal_results.csv"), index=False)
kruskal_df
```
## Summary Table
```{python}
#| label: summary-table
summary = (
freq_df
.groupby(["Segment", "Threshold", "Source"])[
["Mutations_per_1000bp", "Mutations_per_100aa"]
]
.agg(["mean", "std"])
.round(3)
)
summary.columns = ["_".join(c) for c in summary.columns]
summary.reset_index().to_csv(os.path.join(OUTPUT_DIR, "mutation_summary.csv"), index=False)
print("Saved: mutation_summary.csv")
summary
```
