895 lines
33 KiB
Python
895 lines
33 KiB
Python
#!/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=<ID=FFPE,Description="FFPE deamination suspect (end-associated C>T/G>A, strand-biased)">'
|
|
header_done = False
|
|
for line in fin:
|
|
if line.startswith("##FILTER=<ID=FFPE"):
|
|
fout_clean.write(line)
|
|
fout_flagged.write(line)
|
|
header_done = True
|
|
continue
|
|
if line.startswith("#CHROM"):
|
|
if not header_done:
|
|
fout_clean.write(ffpe_header + "\n")
|
|
fout_flagged.write(ffpe_header + "\n")
|
|
fout_clean.write(line)
|
|
fout_flagged.write(line)
|
|
continue
|
|
if line.startswith("#"):
|
|
fout_clean.write(line)
|
|
fout_flagged.write(line)
|
|
continue
|
|
n_total += 1
|
|
parts = line.rstrip("\n").split("\t")
|
|
if len(parts) < 8:
|
|
fout_clean.write(line)
|
|
continue
|
|
chrom, pos_s, _, ref, alt_s = parts[0], parts[1], parts[2], parts[3], parts[4]
|
|
try:
|
|
pos = int(pos_s)
|
|
except ValueError:
|
|
fout_clean.write(line)
|
|
continue
|
|
alts = alt_s.split(",")
|
|
is_flagged = any((chrom, pos, ref, a) in flagged_keys for a in alts)
|
|
filt_col = parts[6] if len(parts) > 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()
|