diff --git a/.gitignore b/.gitignore index e60fc67..39b15e1 100644 --- a/.gitignore +++ b/.gitignore @@ -3,3 +3,14 @@ __pycache__/ *.pyc test_data/ FFPE_QC/ +FFPE_sample/ +*.bam +*.bai +*.fai +hg38.fa +hg38.fa.gz +*.vcf.gz +*.tbi +*.zip +README.html +124/ diff --git a/ffpe_damage_v2.py b/ffpe_damage_v2.py index cc9ea89..fdaad61 100644 --- a/ffpe_damage_v2.py +++ b/ffpe_damage_v2.py @@ -63,6 +63,17 @@ def parse_args(): 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)") return ap.parse_args() @@ -317,6 +328,23 @@ def write_end_enrichment(d, args): 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 = [] @@ -380,9 +408,46 @@ def write_beds_and_variants(d, args): 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") + 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.tsv"), + index=False, sep="\t") + + clean = df[~df["is_ffpe_suspect"]] + damaged = df[df["is_ffpe_suspect"]] + clean.to_csv(os.path.join(args.outdir, "clean_variants.tsv"), + index=False, sep="\t") + damaged.to_csv(os.path.join(args.outdir, "ffpe_suspect_variants.tsv"), + index=False, sep="\t") + + 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.tsv"), "w").close() + open(os.path.join(args.outdir, "ffpe_suspect_variants.tsv"), "w").close() + return pd.DataFrame() def write_summary(d, args): @@ -454,15 +519,106 @@ def plot_cg_profiles(cg5, cg3, outdir): 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) bam = pysam.AlignmentFile(args.bam, "rb") if not bam.has_index(): + try: + print(f"[info] BAM index not found, creating {args.bam}.bai ...") + pysam.index(args.bam) + bam.close() + bam = pysam.AlignmentFile(args.bam, "rb") + except Exception as e: + bam.close() + raise RuntimeError( + f"BAM index (.bai) not found and auto-creation failed: {e}. " + f"Run: samtools index {args.bam}") from e + try: + fasta = pysam.FastaFile(args.reference) + except Exception as e: bam.close() - raise RuntimeError("BAM index (.bai) not found.") - fasta = pysam.FastaFile(args.reference) + 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) + fasta = pysam.FastaFile(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 d = process_bam(bam, fasta, args) bam.close() @@ -473,11 +629,14 @@ def main(): cg5, cg3 = write_cg_profiles(d, args) write_strand_profiles(d, args) write_end_enrichment(d, args) - write_beds_and_variants(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']}") @@ -488,6 +647,11 @@ def main(): 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.tsv / clean_variants.tsv") if __name__ == "__main__":