#!/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.csv (per-position context) clean_variants.csv / ffpe_suspect_variants.csv 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 multiprocessing import os from collections import Counter, 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 get_default_threads(cap=8): try: n = len(os.sched_getaffinity(0)) except AttributeError: n = os.cpu_count() or 1 if n is None: n = 1 n = max(1, n - 1) return min(n, cap) def _process_shard(task): bam_path, ref_path, chrom, mapq, baseq, include_dup = task bam = pysam.AlignmentFile(bam_path, "rb") fasta = pysam.FastaFile(ref_path) try: chrom_seq = fasta.fetch(chrom) except Exception: chrom_seq = None sub_counts = Counter() ref_base_counts = Counter() profile = {} profile3 = {} strand_profile = {} strand_profile3 = {} genomic = {} total_reads = 0 used_reads = 0 usable_bases = 0 def _inc(store, k1, k2, k3, k4): a = store.setdefault(k1, {}) b = a.setdefault(k2, {}) c = b.setdefault(k3, {}) c[k4] = c.get(k4, 0) + 1 for read in bam.fetch(chrom): total_reads += 1 if read.is_secondary or read.is_supplementary: continue if read.mapping_quality < mapq: continue if read.is_duplicate and not include_dup: continue seq = read.query_sequence quals = read.query_qualities if seq is None or quals is None: continue used_reads += 1 read_label = "R1" if read.is_read1 else "R2" if read.is_read2 else "single" strand = "-" if read.is_reverse else "+" read_len = read.query_length chrom_name = bam.get_reference_name(read.reference_id) seq_chrom = chrom_seq if chrom_name == chrom else None if seq_chrom is None: try: seq_chrom = fasta.fetch(chrom_name) except Exception: seq_chrom = None for query_pos, ref_pos in read.get_aligned_pairs(matches_only=True): if quals[query_pos] < baseq: continue if seq_chrom is not None and 0 <= ref_pos < len(seq_chrom): ref_base = seq_chrom[ref_pos].upper() else: try: ref_base = fasta.fetch(chrom_name, ref_pos, ref_pos + 1).upper() except Exception: continue 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 key = (chrom_name, ref_pos) g = genomic.get(key) if g is None: g = {"ref": ref_base, "depth": 0, "depth_plus": 0, "depth_minus": 0, "alts": {}} genomic[key] = g else: g["ref"] = ref_base g["depth"] += 1 if strand == "+": g["depth_plus"] += 1 else: g["depth_minus"] += 1 _inc(profile, read_label, pos5, ref_base, "total") _inc(profile3, read_label, pos3, ref_base, "total") _inc(strand_profile, strand, pos5, ref_base, "total") _inc(strand_profile3, strand, pos3, ref_base, "total") if ref_base == alt_base: continue sub = f"{ref_base}>{alt_base}" sub_counts[sub] += 1 _inc(profile, read_label, pos5, ref_base, sub) _inc(profile3, read_label, pos3, ref_base, sub) _inc(strand_profile, strand, pos5, ref_base, sub) _inc(strand_profile3, strand, pos3, ref_base, sub) alts = g["alts"] alt = alts.get(alt_base) if alt is None: alt = {"count": 0, "bq_sum": 0, "pos5_sum": 0, "pos3_sum": 0, "plus": 0, "minus": 0, "r1": 0, "r2": 0} alts[alt_base] = alt 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 bam.close() fasta.close() return { "total_reads": total_reads, "used_reads": used_reads, "usable_bases": usable_bases, "sub_counts": dict(sub_counts), "ref_base_counts": dict(ref_base_counts), "profile": profile, "profile3": profile3, "strand_profile": strand_profile, "strand_profile3": strand_profile3, "genomic": genomic, } def _merge_nested(target, source): for k1, v1 in source.items(): t1 = target.setdefault(k1, {}) for k2, v2 in v1.items(): t2 = t1.setdefault(k2, {}) for k3, v3 in v2.items(): t3 = t2.setdefault(k3, {}) for k4, cnt in v3.items(): t3[k4] = t3.get(k4, 0) + cnt def _merge_results(partials): merged = { "total_reads": 0, "used_reads": 0, "usable_bases": 0, "sub_counts": Counter(), "ref_base_counts": Counter(), "profile": {}, "profile3": {}, "strand_profile": {}, "strand_profile3": {}, "genomic": {}, } for p in partials: merged["total_reads"] += p["total_reads"] merged["used_reads"] += p["used_reads"] merged["usable_bases"] += p["usable_bases"] merged["sub_counts"].update(p["sub_counts"]) merged["ref_base_counts"].update(p["ref_base_counts"]) _merge_nested(merged["profile"], p["profile"]) _merge_nested(merged["profile3"], p["profile3"]) _merge_nested(merged["strand_profile"], p["strand_profile"]) _merge_nested(merged["strand_profile3"], p["strand_profile3"]) for key, g in p["genomic"].items(): mg = merged["genomic"].get(key) if mg is None: merged["genomic"][key] = g else: mg["depth"] += g["depth"] mg["depth_plus"] += g["depth_plus"] mg["depth_minus"] += g["depth_minus"] for alt_base, alt in g["alts"].items(): malt = mg["alts"].get(alt_base) if malt is None: mg["alts"][alt_base] = alt else: malt["count"] += alt["count"] malt["bq_sum"] += alt["bq_sum"] malt["pos5_sum"] += alt["pos5_sum"] malt["pos3_sum"] += alt["pos3_sum"] malt["plus"] += alt["plus"] malt["minus"] += alt["minus"] malt["r1"] += alt["r1"] malt["r2"] += alt["r2"] return merged 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") ap.add_argument("--vcf", default=None, help="optional VCF to filter: FFPE-suspect variants are " "marked FILTER=FFPE or removed") ap.add_argument("--keep-damaged", action="store_true", help="when --vcf is given, keep damaged variants with " "FILTER=FFPE instead of removing them") ap.add_argument("--filter-end", type=int, default=10, help="mean read-position distance from either end that " "counts as end-associated (default 10)") ap.add_argument("--filter-strand-p", type=float, default=0.05, help="Fisher p-value threshold for strand bias (default 0.05)") ap.add_argument("--threads", "-t", type=int, default=None, help="threads for BAM processing (default auto = " "detected cores -1, capped at 8; 1 = single-thread)") 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 = _get_nested(store, label, pos, ref, "total") if total == 0: continue for alt in "ACGT": if alt == ref: continue sub = f"{ref}>{alt}" count = _get_nested(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 _get_nested(store, label, pos, ref, field): return store.get(label, {}).get(pos, {}).get(ref, {}).get(field, 0) def cg_rows(store, pos_key): out = [] for label in store: for pos in sorted(store[label]): c_total = _get_nested(store, label, pos, "C", "total") ct = _get_nested(store, label, pos, "C", "C>T") g_total = _get_nested(store, label, pos, "G", "total") ga = _get_nested(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 += _get_nested(store, label, pos, "C", "total") ct_num += _get_nested(store, label, pos, "C", "C>T") ga_den += _get_nested(store, label, pos, "G", "total") ga_num += _get_nested(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 _is_ffpe_suspect(row, end_thresh, strand_p_thresh): sub = f"{row['ref']}>{row['alt']}" if sub not in ("C>T", "G>A"): return False, "" is_end = (row["mean_read_position_5prime"] <= end_thresh or row["mean_read_position_3prime"] <= end_thresh) if not is_end: return False, "" reasons = [f"end<={end_thresh}"] frac_plus = row["alt_fraction_plus"] if frac_plus <= 0.1 or frac_plus >= 0.9: reasons.append(f"strand_frac={frac_plus:.2f}") if row["strand_bias_pvalue_fisher"] < strand_p_thresh: reasons.append(f"strand_p={row['strand_bias_pvalue_fisher']:.2g}") return True, ";".join(reasons) 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: df = pd.DataFrame(var_rows) flags = df.apply(lambda r: _is_ffpe_suspect(r, args.filter_end, args.filter_strand_p), axis=1) df["is_ffpe_suspect"] = [x[0] for x in flags] df["ffpe_reason"] = [x[1] for x in flags] df.to_csv(os.path.join(args.outdir, "candidate_variants.csv"), index=False) clean = df[~df["is_ffpe_suspect"]] damaged = df[df["is_ffpe_suspect"]] clean.to_csv(os.path.join(args.outdir, "clean_variants.csv"), index=False) damaged.to_csv(os.path.join(args.outdir, "ffpe_suspect_variants.csv"), index=False) for name, subset, fname in [ ("clean", clean, "candidate_CtoT_clean.bed"), ("clean", clean, "candidate_GtoA_clean.bed"), ]: pass clean_ct = [] clean_ga = [] for _, r in clean.iterrows(): line = (f"{r['chrom']}\t{r['position']}\t{r['position']+1}\t" f"{r['ref']}>{r['alt']};depth={int(r['depth'])};" f"alt={int(r['alt_count'])};VAF={r['VAF']:.3f}\n") if r["ref"] == "C" and r["alt"] == "T": clean_ct.append(line) if r["ref"] == "G" and r["alt"] == "A": clean_ga.append(line) with open(os.path.join(args.outdir, "candidate_CtoT_clean.bed"), "w") as f: f.writelines(clean_ct) with open(os.path.join(args.outdir, "candidate_GtoA_clean.bed"), "w") as f: f.writelines(clean_ga) return df else: open(os.path.join(args.outdir, "clean_variants.csv"), "w").close() open(os.path.join(args.outdir, "ffpe_suspect_variants.csv"), "w").close() return pd.DataFrame() 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 _filter_vcf(args, flagged_df): vcf_in = args.vcf outdir = args.outdir base = os.path.splitext(os.path.basename(vcf_in))[0] if base.endswith(".vcf"): base = base[:-4] clean_path = os.path.join(outdir, f"{base}.clean.vcf") flagged_path = os.path.join(outdir, f"{base}.ffpe_flagged.vcf") flagged_keys = set() for _, r in flagged_df.iterrows(): if r.get("is_ffpe_suspect"): flagged_keys.add((str(r["chrom"]), int(r["position"]) + 1, str(r["ref"]), str(r["alt"]))) n_total = n_kept = n_flagged = 0 with open(vcf_in) as fin, \ open(clean_path, "w") as fout_clean, \ open(flagged_path, "w") as fout_flagged: ffpe_header = '##FILTER=T/G>A, strand-biased)">' header_done = False for line in fin: if line.startswith("##FILTER= 6 else "." if is_flagged: n_flagged += 1 if args.keep_damaged: new_filt = "FFPE" if filt_col in (".", "PASS", "") else filt_col + ";FFPE" parts[6] = new_filt fout_clean.write("\t".join(parts) + "\n") fout_flagged.write(line) else: fout_flagged.write(line) else: n_kept += 1 fout_clean.write(line) print(f"[vcf] total {n_total}, kept {n_kept}, " f"flagged {n_flagged} -> {clean_path} / {flagged_path}") def main(): args = parse_args() os.makedirs(args.outdir, exist_ok=True) threads = args.threads if threads is None: threads = get_default_threads() if threads < 1: threads = 1 bam_probe = pysam.AlignmentFile(args.bam, "rb") has_index = bam_probe.has_index() bam_probe.close() if not has_index: try: print(f"[info] BAM index not found, creating {args.bam}.bai ...") pysam.index(args.bam) except Exception as e: raise RuntimeError( f"BAM index (.bai) not found and auto-creation failed: {e}. " f"Run: samtools index {args.bam}") from e try: fasta_probe = pysam.FastaFile(args.reference) fasta_probe.close() except Exception as e: if "fai" in str(e).lower() or "index" in str(e).lower(): try: print(f"[info] FASTA index not found, creating {args.reference}.fai ...") pysam.faidx(args.reference) except Exception as e2: raise RuntimeError( f"FASTA index (.fai) not found and auto-creation failed: {e2}. " f"Run: samtools faidx {args.reference}") from e2 else: raise if threads == 1: bam = pysam.AlignmentFile(args.bam, "rb") fasta = pysam.FastaFile(args.reference) d = process_bam(bam, fasta, args) bam.close() fasta.close() else: print(f"[info] using {threads} threads (auto, capped at 8, 1 reserved)") bam_tmp = pysam.AlignmentFile(args.bam, "rb") chroms = list(bam_tmp.references) bam_tmp.close() if not chroms: bam = pysam.AlignmentFile(args.bam, "rb") fasta = pysam.FastaFile(args.reference) d = process_bam(bam, fasta, args) bam.close() fasta.close() else: tasks = [(args.bam, args.reference, c, args.mapq, args.baseq, args.include_duplicates) for c in chroms] with multiprocessing.Pool(threads) as pool: partials = pool.map(_process_shard, tasks) d = _merge_results(partials) try: bam_cnt = pysam.AlignmentFile(args.bam, "rb") total_via_idx = sum(s.mapped for s in bam_cnt.get_index_statistics()) total_via_idx += bam_cnt.unmapped bam_cnt.close() if total_via_idx > d["total_reads"]: d["total_reads"] = total_via_idx except Exception: pass 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) flagged_df = write_beds_and_variants(d, args) write_summary(d, args) plot_all12(sub_df, args.outdir) plot_cg_profiles(cg5, cg3, args.outdir) if args.vcf: _filter_vcf(args, flagged_df) 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 not flagged_df.empty: n_flag = int(flagged_df["is_ffpe_suspect"].sum()) n_clean = len(flagged_df) - n_flag print(f"variants flagged FFPE : {n_flag} / {len(flagged_df)} " f"-> ffpe_suspect_variants.csv / clean_variants.csv") if __name__ == "__main__": main()