FFPE deamination QC pipeline

- ffpe_damage_v2.py: normalized C>T/G>A profiling (all 12 substitutions,
  R1/R2 and strand profiles, 5'/3' read-end distance, BED tracks, plots)
- ffpe_compare.py: low/moderate/high classification vs control samples
- make_test_data.py: synthetic BAM with known damage for validation

.gitignore

@@ -0,0 +1,5 @@
+.venv/
+__pycache__/
+*.pyc
+test_data/
+FFPE_QC/

ffpe_compare.py

@@ -0,0 +1,179 @@
+#!/usr/bin/env python3
+"""Classify a query sample's FFPE damage relative to validated controls.
+
+The classification is made ONLY by comparing measurable end-enrichment
+frequencies of the query against the distribution of control samples that
+were produced with the same library protocol.  There is no universal
+absolute threshold: "low / moderate / high" is defined relative to the
+controls.
+
+Metrics used (from ffpe_damage_v2.py --outdir/end_enrichment.csv, window W):
+  * C_to_T_R1_5prime_W : C>T frequency in the first W bases of R1 (forward)
+  * G_to_A_R2_3prime_W : G>A frequency in the last W bases of R2 (reverse)
+  * combined           : mean of the two damage-relevant metrics
+
+Classification per metric (controls mean M, SD S, z = (query - M)/S):
+  high      : z >  k      (query clearly above the control distribution)
+  moderate  : -k <= z <= k
+  low       : z <  -k
+
+Outputs (--out):
+  comparison_table.csv
+  classification.csv
+  control_vs_query.png
+"""
+
+import argparse
+import os
+
+import pandas as pd
+import matplotlib
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+
+
+def parse_args():
+    ap = argparse.ArgumentParser(description=__doc__)
+    ap.add_argument("--controls", nargs="+", required=True,
+                    help="dirs with ffpe_damage_v2.py outputs for controls, "
+                         "or a file with one dir per line")
+    ap.add_argument("--query", required=True,
+                    help="dir with ffpe_damage_v2.py outputs for the sample")
+    ap.add_argument("--out", default="compare_out")
+    ap.add_argument("--window", type=int, default=5,
+                    help="end window size (bases from read end)")
+    ap.add_argument("--k", type=float, default=2.0,
+                    help="number of SDs defining moderate band")
+    return ap.parse_args()
+
+
+def load_samples(dirs):
+    paths = []
+    for d in dirs:
+        if os.path.isfile(d):
+            with open(d) as fh:
+                paths.extend(line.strip() for line in fh if line.strip())
+        else:
+            paths.append(d)
+    out = {}
+    for p in paths:
+        csv_path = os.path.join(p, "end_enrichment.csv")
+        if not os.path.isfile(csv_path):
+            raise FileNotFoundError(f"missing {csv_path}")
+        e = pd.read_csv(csv_path)
+        parent = os.path.basename(os.path.dirname(os.path.normpath(p)))
+        name = parent or os.path.basename(os.path.normpath(p))
+        out[name] = e
+    return out
+
+
+def metric_value(e, read, end, col, window):
+    m = e[(e.read == read) & (e.end == end) & (e.window_size == window)]
+    if m.empty:
+        return float("nan")
+    return float(m[col].iloc[0])
+
+
+def build_metrics(e, window):
+    return {
+        f"C_to_T_R1_5prime_{window}":
+            metric_value(e, "R1", "5prime",
+                         "C_to_T_frequency_percent", window),
+        f"G_to_A_R2_3prime_{window}":
+            metric_value(e, "R2", "3prime",
+                         "G_to_A_frequency_percent", window),
+    }
+
+
+def classify(z, k):
+    if z > k:
+        return "high"
+    if z < -k:
+        return "low"
+    return "moderate"
+
+
+def main():
+    args = parse_args()
+    os.makedirs(args.out, exist_ok=True)
+
+    controls = load_samples(args.controls)
+    query_e = load_samples([args.query])
+    query_name = list(query_e)[0]
+    query_e = query_e[query_name]
+
+    rows = []
+    metric_names = [f"C_to_T_R1_5prime_{args.window}",
+                    f"G_to_A_R2_3prime_{args.window}"]
+    query_values = {}
+
+    for name, e in controls.items():
+        m = build_metrics(e, args.window)
+        rows.append({"sample": name, "type": "control", **m,
+                     "combined": sum(m.values()) / 2})
+    qm = build_metrics(query_e, args.window)
+    q_combined = sum(qm.values()) / 2
+    rows.append({"sample": query_name, "type": "query", **qm,
+                 "combined": q_combined})
+    for mn in metric_names:
+        query_values[mn] = qm[mn]
+
+    df = pd.DataFrame(rows)
+
+    result_rows = []
+    for mn in metric_names + ["combined"]:
+        vals = df.loc[df.type == "control", mn].dropna()
+        mean = vals.mean()
+        sd = vals.std(ddof=1) if len(vals) > 1 else 0.0
+        qv = query_values.get(mn, None) if mn in metric_names else q_combined
+        if sd == 0 or pd.isna(sd):
+            if pd.isna(qv) or qv == mean:
+                z = 0.0
+            else:
+                z = float("inf") if qv > mean else float("-inf")
+        else:
+            z = (qv - mean) / sd
+        cls = classify(z, args.k)
+        result_rows.append({
+            "metric": mn,
+            "control_mean": mean,
+            "control_sd": sd,
+            "query_value": qv,
+            "z_score": z,
+            "classification": cls,
+        })
+
+    result_df = pd.DataFrame(result_rows)
+    result_df.to_csv(os.path.join(args.out, "classification.csv"), index=False)
+    df.to_csv(os.path.join(args.out, "comparison_table.csv"), index=False)
+
+    n_metrics = len(metric_names)
+    fig, axes = plt.subplots(1, n_metrics, figsize=(6 * n_metrics, 5),
+                             squeeze=False)
+    for ax, mn in zip(axes[0], metric_names):
+        cvals = df.loc[df.type == "control", mn].dropna()
+        qv = query_values[mn]
+        mean = cvals.mean()
+        sd = cvals.std(ddof=1) if len(cvals) > 1 else 0.0
+        ax.scatter(range(len(cvals)), cvals, color="#1f77b4", s=60,
+                   label="controls")
+        ax.axhline(mean, color="#7f7f7f", linestyle="--", label="control mean")
+        if sd > 0:
+            ax.axhline(mean + args.k * sd, color="#ff7f0e", linestyle=":",
+                       label=f"mean ± k·SD (k={args.k:g})")
+            ax.axhline(mean - args.k * sd, color="#ff7f0e", linestyle=":")
+        ax.scatter([len(cvals)], [qv], color="#d62728", s=90, marker="x",
+                   label="query")
+        ax.set_xticks([])
+        ax.set_ylabel("frequency (%)")
+        ax.set_title(mn)
+        ax.legend(fontsize=8)
+    fig.tight_layout()
+    fig.savefig(os.path.join(args.out, "control_vs_query.png"), dpi=300)
+    plt.close(fig)
+
+    print(result_df.to_string(index=False))
+
+
+if __name__ == "__main__":
+    main()

ffpe_damage_v2.py

@@ -0,0 +1,494 @@
+#!/usr/bin/env python3
+"""FFPE deamination profiling from an aligned BAM.
+
+Methodology
+-----------
+For every read position i the normalized frequency is:
+
+    C>T_i = #(C->T at position i) / #(C observations at position i)
+    G>A_i = #(G->A at position i) / #(G observations at position i)
+
+where the denominator counts ALL high-quality observations of that
+reference base at position i, including matches.
+
+Separate profiles are produced for R1 / R2 and for forward / reverse
+alignment strand, using both distance-from-5'-end and distance-from-3'-end.
+
+All 12 substitution types are counted globally so that C>T / G>A can be
+judged against the other mismatch classes.
+
+This tool reports measurable frequencies only; it does NOT assign a
+low/moderate/high classification.  Classification requires comparison
+against validated control samples (see ffpe_compare.py).
+
+Outputs (in --outdir):
+  sample_summary.csv
+  substitution_summary.csv
+  substitution_all12.png
+  normalized_damage_by_read_position.csv          (5' profile)
+  normalized_damage_by_read_position_3prime.csv   (3' profile)
+  CtoT_GtoA_normalized_profile.csv                (5' profile)
+  CtoT_GtoA_normalized_profile_3prime.csv         (3' profile)
+  strand_damage_profile.csv                       (5' profile)
+  strand_damage_profile_3prime.csv                (3' profile)
+  end_enrichment.csv
+  candidate_CtoT.bed / candidate_GtoA.bed         (IGV tracks)
+  candidate_variants.tsv                          (per-position context)
+  CtoT_R1.png CtoT_R2.png GtoA_R1.png GtoA_R2.png (5' plots)
+  CtoT_R1_3prime.png ... GtoA_R2_3prime.png       (3' plots)
+"""
+
+import argparse
+import os
+from collections import defaultdict
+
+import pysam
+import pandas as pd
+import matplotlib
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+from scipy.stats import fisher_exact
+
+SUBSTITUTIONS = [f"{r}>{a}" for r in "ACGT" for a in "ACGT" if r != a]
+END_SIZES = [1, 3, 5, 10, 20]
+
+
+def parse_args():
+    ap = argparse.ArgumentParser(description="FFPE damage analysis")
+    ap.add_argument("--bam", required=True)
+    ap.add_argument("--reference", required=True)
+    ap.add_argument("--outdir", default="FFPE_QC")
+    ap.add_argument("--mapq", type=int, default=20)
+    ap.add_argument("--baseq", type=int, default=20)
+    ap.add_argument("--min-depth", type=int, default=20)
+    ap.add_argument("--min-alt-count", type=int, default=5)
+    ap.add_argument("--include-duplicates", action="store_true")
+    return ap.parse_args()
+
+
+def process_bam(bam, fasta, args):
+    total_reads = 0
+    used_reads = 0
+    usable_bases = 0
+
+    sub_counts = defaultdict(int)
+    ref_base_counts = defaultdict(int)
+
+    profile = defaultdict(lambda: defaultdict(lambda: defaultdict(
+        lambda: defaultdict(int))))
+    profile3 = defaultdict(lambda: defaultdict(lambda: defaultdict(
+        lambda: defaultdict(int))))
+    strand_profile = defaultdict(lambda: defaultdict(lambda: defaultdict(
+        lambda: defaultdict(int))))
+    strand_profile3 = defaultdict(lambda: defaultdict(lambda: defaultdict(
+        lambda: defaultdict(int))))
+
+    genomic = defaultdict(lambda: {
+        "ref": None,
+        "depth": 0,
+        "depth_plus": 0,
+        "depth_minus": 0,
+        "alts": defaultdict(lambda: {
+            "count": 0,
+            "bq_sum": 0,
+            "pos5_sum": 0,
+            "pos3_sum": 0,
+            "plus": 0,
+            "minus": 0,
+            "r1": 0,
+            "r2": 0,
+        }),
+    })
+
+    for read in bam.fetch(until_eof=True):
+        total_reads += 1
+        if read.is_unmapped or read.is_secondary or read.is_supplementary:
+            continue
+        if read.mapping_quality < args.mapq:
+            continue
+        if read.is_duplicate and not args.include_duplicates:
+            continue
+        seq = read.query_sequence
+        quals = read.query_qualities
+        if seq is None or quals is None:
+            continue
+        used_reads += 1
+
+        if read.is_read1:
+            read_label = "R1"
+        elif read.is_read2:
+            read_label = "R2"
+        else:
+            read_label = "single"
+        strand = "-" if read.is_reverse else "+"
+        read_len = read.query_length
+        chrom = bam.get_reference_name(read.reference_id)
+
+        for query_pos, ref_pos in read.get_aligned_pairs(matches_only=True):
+            if quals[query_pos] < args.baseq:
+                continue
+            ref_base = fasta.fetch(chrom, ref_pos, ref_pos + 1).upper()
+            if ref_base not in "ACGT":
+                continue
+            alt_base = seq[query_pos].upper()
+            if alt_base not in "ACGT":
+                continue
+            usable_bases += 1
+
+            pos5 = query_pos + 1
+            pos3 = read_len - query_pos
+
+            ref_base_counts[ref_base] += 1
+
+            g = genomic[(chrom, ref_pos)]
+            g["ref"] = ref_base
+            g["depth"] += 1
+            if strand == "+":
+                g["depth_plus"] += 1
+            else:
+                g["depth_minus"] += 1
+
+            profile[read_label][pos5][ref_base]["total"] += 1
+            profile3[read_label][pos3][ref_base]["total"] += 1
+            strand_profile[strand][pos5][ref_base]["total"] += 1
+            strand_profile3[strand][pos3][ref_base]["total"] += 1
+
+            if ref_base == alt_base:
+                continue
+
+            sub = f"{ref_base}>{alt_base}"
+            sub_counts[sub] += 1
+
+            profile[read_label][pos5][ref_base][sub] += 1
+            profile3[read_label][pos3][ref_base][sub] += 1
+            strand_profile[strand][pos5][ref_base][sub] += 1
+            strand_profile3[strand][pos3][ref_base][sub] += 1
+
+            alt = g["alts"][alt_base]
+            alt["count"] += 1
+            alt["bq_sum"] += quals[query_pos]
+            alt["pos5_sum"] += pos5
+            alt["pos3_sum"] += pos3
+            if strand == "+":
+                alt["plus"] += 1
+            else:
+                alt["minus"] += 1
+            if read_label == "R1":
+                alt["r1"] += 1
+            elif read_label == "R2":
+                alt["r2"] += 1
+
+    return {
+        "total_reads": total_reads,
+        "used_reads": used_reads,
+        "usable_bases": usable_bases,
+        "sub_counts": sub_counts,
+        "ref_base_counts": ref_base_counts,
+        "profile": profile,
+        "profile3": profile3,
+        "strand_profile": strand_profile,
+        "strand_profile3": strand_profile3,
+        "genomic": genomic,
+    }
+
+
+def write_substitution_summary(d, args):
+    rows = []
+    for sub in SUBSTITUTIONS:
+        ref = sub[0]
+        count = d["sub_counts"][sub]
+        denom = d["ref_base_counts"][ref]
+        freq = count / denom if denom else 0.0
+        rows.append({
+            "substitution": sub,
+            "count": count,
+            "reference_base_observations": denom,
+            "frequency": freq,
+            "frequency_percent": freq * 100,
+        })
+    df = pd.DataFrame(rows)
+    df.to_csv(os.path.join(args.outdir, "substitution_summary.csv"),
+              index=False)
+    return df
+
+
+def write_read_position_profiles(d, args):
+    def rows_for(store, pos_key, side):
+        out = []
+        for label in store:
+            for pos in sorted(store[label]):
+                for ref in "ACGT":
+                    total = store[label][pos][ref]["total"]
+                    if total == 0:
+                        continue
+                    for alt in "ACGT":
+                        if alt == ref:
+                            continue
+                        sub = f"{ref}>{alt}"
+                        count = store[label][pos][ref][sub]
+                        out.append({
+                            "read": label,
+                            pos_key: pos,
+                            "end_side": side,
+                            "reference_base": ref,
+                            "substitution": sub,
+                            "count": count,
+                            "denominator": total,
+                            "frequency": count / total,
+                            "frequency_percent": count / total * 100,
+                        })
+        return out
+
+    df5 = pd.DataFrame(rows_for(d["profile"], "position_5prime", "5prime"))
+    df3 = pd.DataFrame(rows_for(d["profile3"], "position_3prime", "3prime"))
+    df5.to_csv(os.path.join(args.outdir,
+                            "normalized_damage_by_read_position.csv"),
+               index=False)
+    df3.to_csv(os.path.join(args.outdir,
+                            "normalized_damage_by_read_position_3prime.csv"),
+               index=False)
+    return df5, df3
+
+
+def cg_rows(store, pos_key):
+    out = []
+    for label in store:
+        for pos in sorted(store[label]):
+            c_total = store[label][pos]["C"]["total"]
+            ct = store[label][pos]["C"]["C>T"]
+            g_total = store[label][pos]["G"]["total"]
+            ga = store[label][pos]["G"]["G>A"]
+            out.append({
+                "read": label,
+                pos_key: pos,
+                "C_total": c_total,
+                "C_to_T": ct,
+                "C_to_T_frequency_percent": ct / c_total * 100 if c_total else 0,
+                "G_total": g_total,
+                "G_to_A": ga,
+                "G_to_A_frequency_percent": ga / g_total * 100 if g_total else 0,
+            })
+    return out
+
+
+def write_cg_profiles(d, args):
+    df5 = pd.DataFrame(cg_rows(d["profile"], "position_5prime"))
+    df3 = pd.DataFrame(cg_rows(d["profile3"], "position_3prime"))
+    df5.to_csv(os.path.join(args.outdir,
+                            "CtoT_GtoA_normalized_profile.csv"), index=False)
+    df3.to_csv(os.path.join(args.outdir,
+                            "CtoT_GtoA_normalized_profile_3prime.csv"),
+               index=False)
+    return df5, df3
+
+
+def write_strand_profiles(d, args):
+    df5 = pd.DataFrame(cg_rows(d["strand_profile"], "position_5prime"))
+    df3 = pd.DataFrame(cg_rows(d["strand_profile3"], "position_3prime"))
+    df5.to_csv(os.path.join(args.outdir, "strand_damage_profile.csv"),
+               index=False)
+    df3.to_csv(os.path.join(args.outdir,
+                            "strand_damage_profile_3prime.csv"), index=False)
+    return df5, df3
+
+
+def write_end_enrichment(d, args):
+    rows = []
+    for store, side in ((d["profile"], "5prime"), (d["profile3"], "3prime")):
+        for label in store:
+            for size in END_SIZES:
+                ct_num = ct_den = ga_num = ga_den = 0
+                for pos in range(1, size + 1):
+                    ct_den += store[label][pos]["C"]["total"]
+                    ct_num += store[label][pos]["C"]["C>T"]
+                    ga_den += store[label][pos]["G"]["total"]
+                    ga_num += store[label][pos]["G"]["G>A"]
+                rows.append({
+                    "read": label,
+                    "end": side,
+                    "window_size": size,
+                    "C_to_T_frequency_percent":
+                        ct_num / ct_den * 100 if ct_den else 0,
+                    "G_to_A_frequency_percent":
+                        ga_num / ga_den * 100 if ga_den else 0,
+                })
+    df = pd.DataFrame(rows)
+    df.to_csv(os.path.join(args.outdir, "end_enrichment.csv"), index=False)
+    return df
+
+
+def write_beds_and_variants(d, args):
+    ct_lines = []
+    ga_lines = []
+    var_rows = []
+
+    for (chrom, pos), data in d["genomic"].items():
+        depth = data["depth"]
+        if depth < args.min_depth:
+            continue
+        ref = data["ref"]
+        for alt_base, alt in data["alts"].items():
+            count = alt["count"]
+            if count < args.min_alt_count:
+                continue
+            vaf = count / depth
+            name = f"{ref}>{alt_base};depth={depth};alt={count};VAF={vaf:.3f}"
+            start, end = pos, pos + 1
+            line = f"{chrom}\t{start}\t{end}\t{name}\n"
+            if ref == "C" and alt_base == "T":
+                ct_lines.append(line)
+            if ref == "G" and alt_base == "A":
+                ga_lines.append(line)
+
+            plus = alt["plus"]
+            minus = alt["minus"]
+            depth_plus = data["depth_plus"]
+            depth_minus = data["depth_minus"]
+            ref_plus = depth_plus - plus
+            ref_minus = depth_minus - minus
+            table = [[plus, ref_plus], [minus, ref_minus]]
+            fisher_p = 1.0
+            try:
+                if plus + ref_plus > 0 and minus + ref_minus > 0:
+                    _, fisher_p = fisher_exact(table)
+            except Exception:
+                fisher_p = 1.0
+
+            alt_count = alt["count"]
+            mean_pos5 = alt["pos5_sum"] / alt_count
+            mean_pos3 = alt["pos3_sum"] / alt_count
+
+            var_rows.append({
+                "chrom": chrom,
+                "position": pos,
+                "ref": ref,
+                "alt": alt_base,
+                "depth": depth,
+                "alt_count": count,
+                "VAF": vaf,
+                "mean_base_quality": alt["bq_sum"] / alt_count,
+                "mean_read_position_5prime": mean_pos5,
+                "mean_read_position_3prime": mean_pos3,
+                "alt_fraction_plus": plus / alt_count,
+                "strand_bias_pvalue_fisher": fisher_p,
+                "alt_fraction_R1": alt["r1"] / alt_count,
+                "alt_fraction_R2": alt["r2"] / alt_count,
+            })
+
+    with open(os.path.join(args.outdir, "candidate_CtoT.bed"), "w") as f:
+        f.writelines(ct_lines)
+    with open(os.path.join(args.outdir, "candidate_GtoA.bed"), "w") as f:
+        f.writelines(ga_lines)
+    if var_rows:
+        pd.DataFrame(var_rows).to_csv(
+            os.path.join(args.outdir, "candidate_variants.tsv"),
+            index=False, sep="\t")
+
+
+def write_summary(d, args):
+    ct = d["sub_counts"]["C>T"]
+    ga = d["sub_counts"]["G>A"]
+    ct_den = d["ref_base_counts"]["C"]
+    ga_den = d["ref_base_counts"]["G"]
+    row = {
+        "total_reads": d["total_reads"],
+        "used_reads": d["used_reads"],
+        "usable_bases": d["usable_bases"],
+        "C_to_T_count": ct,
+        "C_observations": ct_den,
+        "C_to_T_frequency_percent": ct / ct_den * 100 if ct_den else 0,
+        "G_to_A_count": ga,
+        "G_observations": ga_den,
+        "G_to_A_frequency_percent": ga / ga_den * 100 if ga_den else 0,
+        "MAPQ_threshold": args.mapq,
+        "BQ_threshold": args.baseq,
+        "duplicates_included": args.include_duplicates,
+    }
+    df = pd.DataFrame([row])
+    df.to_csv(os.path.join(args.outdir, "sample_summary.csv"), index=False)
+    return df
+
+
+def plot_all12(sub_df, outdir):
+    x = sub_df["substitution"]
+    y = sub_df["frequency_percent"]
+    plt.figure(figsize=(10, 5))
+    bars = plt.bar(x, y)
+    for b, s in zip(bars, sub_df["substitution"]):
+        if s in ("C>T", "G>A"):
+            b.set_color("#d62728")
+        else:
+            b.set_color("#7f7f7f")
+    plt.ylabel("frequency (%)")
+    plt.title("All 12 substitution types, global frequency")
+    plt.xticks(rotation=45)
+    plt.tight_layout()
+    plt.savefig(os.path.join(outdir, "substitution_all12.png"), dpi=300)
+    plt.close()
+
+
+def plot_cg_profiles(cg5, cg3, outdir):
+    for label in sorted(set(cg5["read"]) | set(cg3["read"])):
+        for sub, col5, col3, prefix in (
+                ("C>T", "C_to_T_frequency_percent",
+                 "C_to_T_frequency_percent", "CtoT"),
+                ("G>A", "G_to_A_frequency_percent",
+                 "G_to_A_frequency_percent", "GtoA")):
+            for side, df in (("5prime", cg5), ("3prime", cg3)):
+                d = df[df["read"] == label]
+                if d.empty:
+                    continue
+                pos_col = ("position_5prime" if side == "5prime"
+                           else "position_3prime")
+                col = col5 if side == "5prime" else col3
+                plt.figure(figsize=(8, 5))
+                plt.plot(d[pos_col], d[col], marker=".", label=f"{label} {sub}")
+                plt.xlabel(f"position in read ({side})")
+                plt.ylabel(f"{sub} frequency (%)")
+                plt.title(f"FFPE damage profile — {label} — {sub} — {side}")
+                plt.legend()
+                plt.tight_layout()
+                suffix = "" if side == "5prime" else "_3prime"
+                plt.savefig(os.path.join(
+                    outdir, f"{prefix}_{label}{suffix}.png"), dpi=300)
+                plt.close()
+
+
+def main():
+    args = parse_args()
+    os.makedirs(args.outdir, exist_ok=True)
+
+    bam = pysam.AlignmentFile(args.bam, "rb")
+    if not bam.has_index():
+        bam.close()
+        raise RuntimeError("BAM index (.bai) not found.")
+    fasta = pysam.FastaFile(args.reference)
+
+    d = process_bam(bam, fasta, args)
+    bam.close()
+    fasta.close()
+
+    sub_df = write_substitution_summary(d, args)
+    write_read_position_profiles(d, args)
+    cg5, cg3 = write_cg_profiles(d, args)
+    write_strand_profiles(d, args)
+    write_end_enrichment(d, args)
+    write_beds_and_variants(d, args)
+    write_summary(d, args)
+    plot_all12(sub_df, args.outdir)
+    plot_cg_profiles(cg5, cg3, args.outdir)
+
+    print(f"total reads     : {d['total_reads']}")
+    print(f"used reads      : {d['used_reads']}")
+    print(f"usable bases    : {d['usable_bases']}")
+    print(f"C>T global      : "
+          f"{d['sub_counts']['C>T']} / {d['ref_base_counts']['C']} = "
+          f"{d['sub_counts']['C>T'] / d['ref_base_counts']['C'] * 100:.3f}%")
+    print(f"G>A global      : "
+          f"{d['sub_counts']['G>A']} / {d['ref_base_counts']['G']} = "
+          f"{d['sub_counts']['G>A'] / d['ref_base_counts']['G'] * 100:.3f}%")
+    print(f"outputs         : {args.outdir}/")
+
+
+if __name__ == "__main__":
+    main()

make_test_data.py

@@ -0,0 +1,214 @@
+#!/usr/bin/env python3
+"""Generate synthetic reference + paired-end BAM with a known damage profile.
+
+Two scenarios:
+  * damaged:  plus-strand deamination C>T at the 5' end of forward reads (R1)
+              and minus-strand deamination G>A at the 3' end of reverse
+              reads (R2, fragment far end) -- the two strands of the same
+              cytosine-deamination event, appearing at opposite read ends
+              + one "real" heterozygous SNP with no read-position / strand
+              bias, which the pipeline must NOT mistake for damage
+  * clean:    background sequencing error only, no end-associated damage
+              (used as a control sample)
+
+Ground truth for end-to-end validation of ffpe_damage_v2.py / ffpe_compare.py.
+"""
+
+import argparse
+import os
+import random
+
+import pysam
+
+
+COMP = str.maketrans("ACGT", "TGCA")
+
+
+def rc(seq):
+    return seq.translate(COMP)[::-1]
+
+
+def make_reference(outdir, n_chroms, length, seed):
+    rng = random.Random(seed)
+    fa_path = os.path.join(outdir, "synthetic_reference.fa")
+    with open(fa_path, "w") as fh:
+        for i in range(n_chroms):
+            name = f"chr{i + 1}"
+            seq = "".join(rng.choice("ACGT") for _ in range(length))
+            fh.write(f">{name}\n")
+            for j in range(0, length, 80):
+                fh.write(seq[j:j + 80] + "\n")
+    pysam.faidx(fa_path)
+    return fa_path
+
+
+def mutate_sequence(seq, positions, base, alt, rng, p):
+    """Set bases at given query positions to alt with probability p."""
+    s = list(seq)
+    for pos in positions:
+        if s[pos] == base and rng.random() < p:
+            s[pos] = alt
+    return "".join(s)
+
+
+def build_bam(fa_path, bam_path, rng, chroms, length, n_frags, read_len,
+              frag_len, damage_rate, snp_pos, snp_alt, background_error,
+              mapq, baseq):
+    fasta = pysam.FastaFile(fa_path)
+
+    header = {
+        "HD": {"VN": "1.6", "SO": "unsorted"},
+        "SQ": [{"SN": c, "LN": length} for c in chroms],
+    }
+
+    with pysam.AlignmentFile(bam_path, "wb", header=header) as out:
+        read_id = 0
+        for chrom in chroms:
+            for _ in range(n_frags):
+                start = rng.randint(0, length - frag_len - 1)
+                frag_start = start
+                frag_end = start + frag_len
+                ref_r1 = fasta.fetch(chrom, frag_start, frag_start + read_len)
+                ref_r2 = fasta.fetch(chrom, frag_end - read_len, frag_end)
+
+                q1 = ref_r1
+                q2 = ref_r2
+
+                # --- damage injection: fragment ends ---
+                # forward-strand read (R1): 5' end of stored read
+                # (plus-strand deamination) -> ref C shows as read T
+                q1 = mutate_sequence(q1, range(min(5, read_len)),
+                                     "C", "T", rng, damage_rate)
+                # reverse-strand read (R2): sequenced 5' end is the fragment's
+                # far end, which lands on the 3' end of the stored read.
+                # (minus-strand deamination) -> ref plus-strand G shows as read A
+                q2 = mutate_sequence(q2, range(read_len - 5, read_len),
+                                     "G", "A", rng, damage_rate)
+
+                # --- "real" SNP, heterozygous, no positional/strand bias ---
+                if frag_start <= snp_pos < frag_end:
+                    if rng.random() < 0.5:
+                        qp_r1 = snp_pos - frag_start
+                        if 0 <= qp_r1 < read_len:
+                            q1 = q1[:qp_r1] + snp_alt + q1[qp_r1 + 1:]
+                    qp_r2 = snp_pos - (frag_end - read_len)
+                    if 0 <= qp_r2 < read_len:
+                        if rng.random() < 0.5:
+                            q2 = q2[:qp_r2] + snp_alt + q2[qp_r2 + 1:]
+
+                # --- background error ---
+                q1 = _background_error(q1, rng, background_error)
+                q2 = _background_error(q2, rng, background_error)
+
+                quals1 = pysam.qualitystring_to_array(
+                    bytes([baseq + 33]) * read_len).tobytes()
+                quals2 = pysam.qualitystring_to_array(
+                    bytes([baseq + 33]) * read_len).tobytes()
+
+                r1 = pysam.AlignedSegment()
+                r1.query_name = f"{chrom}_frag{read_id}_R1"
+                r1.query_sequence = q1
+                r1.flag = 99          # paired, proper, mate reverse, read1
+                r1.reference_id = out.get_tid(chrom)
+                r1.reference_start = frag_start
+                r1.mapping_quality = mapq
+                r1.cigar = ((0, read_len),)
+                r1.query_qualities = quals1
+                r1.next_reference_id = out.get_tid(chrom)
+                r1.next_reference_start = frag_end - read_len
+                r1.template_length = frag_len
+
+                r2 = pysam.AlignedSegment()
+                r2.query_name = f"{chrom}_frag{read_id}_R2"
+                r2.query_sequence = q2
+                r2.flag = 147         # paired, proper, reverse, read2
+                r2.reference_id = out.get_tid(chrom)
+                r2.reference_start = frag_end - read_len
+                r2.mapping_quality = mapq
+                r2.cigar = ((0, read_len),)
+                r2.query_qualities = quals2
+                r2.next_reference_id = out.get_tid(chrom)
+                r2.next_reference_start = frag_start
+                r2.template_length = -frag_len
+
+                out.write(r1)
+                out.write(r2)
+                read_id += 1
+
+    fasta.close()
+
+
+def _background_error(seq, rng, p):
+    s = list(seq)
+    bases = "ACGT"
+    for i in range(len(s)):
+        if rng.random() < p:
+            alt = rng.choice(bases.replace(s[i], ""))
+            s[i] = alt
+    return "".join(s)
+
+
+def main():
+    ap = argparse.ArgumentParser(description=__doc__)
+    ap.add_argument("--outdir", default="test_data")
+    ap.add_argument("--bam-out", default=None,
+                    help="output BAM path (default <outdir>/sample.bam)")
+    ap.add_argument("--seed", type=int, default=42)
+    ap.add_argument("--n-chroms", type=int, default=3)
+    ap.add_argument("--chrom-length", type=int, default=5000)
+    ap.add_argument("--depth", type=int, default=30,
+                    help="approximate per-base read coverage")
+    ap.add_argument("--read-length", type=int, default=100)
+    ap.add_argument("--fragment-length", type=int, default=250)
+    ap.add_argument("--damage-rate", type=float, default=0.5,
+                    help="0.5 = damaged FFPE-like, 0.0 = clean control")
+    ap.add_argument("--background-error", type=float, default=0.001)
+    ap.add_argument("--mapq", type=int, default=60)
+    ap.add_argument("--baseq", type=int, default=35)
+    args = ap.parse_args()
+
+    os.makedirs(args.outdir, exist_ok=True)
+    rng = random.Random(args.seed)
+
+    n_chroms = args.n_chroms
+    length = args.chrom_length
+    chroms = [f"chr{i + 1}" for i in range(n_chroms)]
+
+    fa_path = make_reference(args.outdir, n_chroms, length, args.seed)
+
+    frags = args.depth * length // (2 * args.read_length)
+    snp_pos = None
+    snp_ref = None
+    snp_alt = "T"
+    with pysam.FastaFile(fa_path) as f:
+        for offset in range(length // 2, length):
+            base = f.fetch(chroms[0], offset, offset + 1).upper()
+            if base == "C":
+                snp_pos = offset
+                snp_ref = base
+                break
+        if snp_ref is None:
+            snp_pos = length // 2
+            snp_ref = f.fetch(chroms[0], snp_pos, snp_pos + 1).upper()
+            snp_alt = "A" if snp_ref != "A" else "G"
+
+    bam_path = args.bam_out or os.path.join(args.outdir, "sample.bam")
+    build_bam(fa_path, bam_path, rng, chroms, length, frags,
+              args.read_length, args.fragment_length, args.damage_rate,
+              snp_pos, snp_alt, args.background_error, args.mapq, args.baseq)
+
+    sorted_bam = bam_path.replace(".bam", ".sorted.bam")
+    pysam.sort("-o", sorted_bam, bam_path)
+    pysam.index(sorted_bam)
+
+    print(f"reference : {fa_path}")
+    print(f"bam       : {sorted_bam}")
+    print(f"fragments : {frags} per chromosome")
+    print(f"real SNP  : {chroms[0]}:{snp_pos} {snp_ref}>{snp_alt} (heterozygous)")
+    print(f"damage    : C>T at 5' of R1/forward (fragment 5' end), "
+          f"G>A at 3' of R2/reverse (fragment far end)"
+          f" (rate {args.damage_rate})")
+
+
+if __name__ == "__main__":
+    main()