From 42c5339825a2267a9a9804012ca5f747cad9f460 Mon Sep 17 00:00:00 2001 From: abonney Date: Tue, 28 Jul 2026 15:40:57 -0500 Subject: [PATCH] offtarget: optional per-read XC tags on a target-window BAM for IGV review The ECS caller decides per read whether it supports an edit, aggregates to per-site counts and throws the per-read decision away, so a reviewer sees indel_fraction = 0.42 and has to re-derive by eye which reads made it up. Add --tagged-bam-out to bin/find_edited_reads.py: every read over a target window gets a string tag (default XC) naming its classification, so IGV can colour the pileup by the caller's own reasoning. Two constraints shaped the implementation: * The per-target read loop fetches with a +/-150 bp pad, so overlapping targets make it visit the same alignment record repeatedly. Writing in-loop would emit duplicate records, which IGV renders as inflated depth. Tags are therefore only accumulated in the loop, keyed on (query_name, flag, reference_start) -- query_name alone collides between mates and between primary/supplementary records -- and the BAM is written in a second pass, sorted and indexed. Conflicting calls for one record resolve by a documented precedence (edit > reference > unevaluable > padding). * Output is window-restricted and off by default. On one AAVS1 site14 ECS sample (1149 targets, ~11,000x) it is 800 MB / 21.8 M reads and lifts peak RSS from 1.4 GB to 6.3 GB, so ECS_INDELS gets 24 GB instead of 8 GB while the flag is set. The dominant classification is left implicit rather than stored, which keeps the tag map roughly 40% smaller. Verified on real data: with tagging off the TSV is byte-identical to the unmodified script; with it on, the full-panel BAM has 21,792,054 records, zero missing tags, zero duplicate keys, and is coordinate-sorted and indexable. Note the tag counts do not equal the TSV columns and are not meant to: tags are per alignment record, indel_reads is per fragment and post-filter. Documented. Wires params.offtarget_tagged_bam through ECS_INDELS and GET_INDELS, adds the schema entry, docs with the IGV recipe, and six tests on a new synthetic aligned-read fixture covering dedup, sort order, skip reasons and precedence. --- bin/find_edited_reads.py | 358 ++++++++++++++++++++++++++++++++--- conf/modules.config | 6 +- docs/OFFTARGET_WORKFLOW.md | 58 ++++++ modules/local/ecs_indels.nf | 9 +- modules/local/get_indels.nf | 11 +- nextflow.config | 1 + nextflow_schema.json | 5 + tests/conftest.py | 171 +++++++++++++++++ tests/test_offtarget_glue.py | 146 ++++++++++++++ 9 files changed, 738 insertions(+), 27 deletions(-) diff --git a/bin/find_edited_reads.py b/bin/find_edited_reads.py index a723b56..6abe6a4 100755 --- a/bin/find_edited_reads.py +++ b/bin/find_edited_reads.py @@ -2,7 +2,7 @@ from __future__ import division import edlib -import argparse, re, string, csv, sys +import argparse, re, string, csv, os, sys import scipy.stats as stats import joblib import pandas as pd @@ -1182,8 +1182,15 @@ def calculate_distance_to_closest_pam(position, targets_df, chrom): distances = np.abs(targets_df['Start'] - position) return distances.min() if len(distances) > 0 else -1 -def predict_reads_at_position(bam_file, chrom, start, end, pampos, model, fasta, is_on_target=0, control_bam=None, threshold=0.80): - +def predict_reads_at_position(bam_file, chrom, start, end, pampos, model, fasta, is_on_target=0, control_bam=None, threshold=0.80, read_probs=None): + """Score reads at one locus with the CRISPR model. + + Returns (reads at or above threshold, mean probability as a percentage). If + ``read_probs`` is a dict it is additionally filled with per-read + probabilities keyed by read_tag_key(), keeping the highest score seen for a + record across overlapping loci -- used to emit the optional XP BAM tag. + """ + reads = [] for read in bam_file.fetch(chrom, start, end): if read.is_unmapped or read.is_duplicate: @@ -1293,9 +1300,16 @@ def predict_reads_at_position(bam_file, chrom, start, end, pampos, model, fasta, # Get model predictions preds = model.predict_proba(features_df)[:, 1] + if read_probs is not None: + for read, prob in zip(reads, preds): + key = read_tag_key(read) + prob = float(prob) + if prob > read_probs.get(key, -1.0): + read_probs[key] = prob + # Calculate results avg_probability = float(preds.mean() * 100) - + return int((preds >= threshold).sum()), avg_probability def merge_dicts_to_tuples(data): @@ -1481,6 +1495,227 @@ def write_vcf_output(df, outfile_name, vcf_header=None, sample_name="EDITED"): vcf_out.close() +# ============================================================================ +# SECTION 2b: READ-LEVEL TAGS (optional, for IGV review) +# ============================================================================ +# +# The per-target read loop in main() already decides, for every read it sees, +# whether that read supports an edit -- and then throws the decision away, so a +# reviewer only ever sees the aggregated indel_fraction. These helpers keep the +# per-read verdict so it can be written back out as a string BAM tag that IGV +# can colour by (Color alignments by -> tag -> XC). +# +# Two properties of the read loop shape this code: +# +# 1. The loop fetches per target with a +/- target_window pad and +# multiple_iterators=True, so overlapping targets make it visit the *same* +# alignment record more than once. Writing from inside the loop would emit +# that record repeatedly, which IGV renders as inflated depth. So tags are +# only accumulated during the loop; the BAM is written in a second pass +# that visits each record exactly once (see write_tagged_bam). +# +# 2. query_name is not a unique key -- it collides between mates of a pair and +# between the primary and supplementary records of one fragment. The key is +# (query_name, flag, reference_start). + +# Tag precedence, most specific first. When the same record is classified +# differently at two overlapping targets, the lowest-index prefix wins: real edit +# evidence beats a plain reference call, which beats a read a classifier looked at +# but could not place, which beats a read that simply sat in the padding. Matching +# is by prefix and first-match, so the specific Skipped_ entries must precede the +# generic one. +TAG_PRECEDENCE = ( + 'Edited_BND', + 'Edited_Deletion', + 'Edited_Insertion', + 'Edited_Duplication', + 'Edited_Complex', + 'Edited_SoftClip', + 'Unedited_WT', + 'Skipped_Unevaluable', + 'Skipped_NoSpan', + 'Skipped_', +) + +# The tag assumed for any record not present in the tag map. This is the single +# biggest category by far -- at a 1 bp target the +/-150 bp fetch pad means most +# records neither span the target nor carry an event (69% of records at the AAVS1 +# on-target) -- so leaving it implicit rather than storing it cuts the tag map to +# roughly a third of its size on a real panel. It is safe to leave implicit only +# because it is *last* in TAG_PRECEDENCE among the tags a read can also receive +# elsewhere: anything else recorded for the same read outranks it anyway. +DEFAULT_TAG = 'Skipped_NoSpan' + + +def read_tag_rank(tag): + """Precedence rank of a tag string; lower wins. Unknown tags rank last.""" + for i, prefix in enumerate(TAG_PRECEDENCE): + if tag.startswith(prefix): + return i + return len(TAG_PRECEDENCE) + + +def read_tag_key(read): + """Unique key for one alignment record. + + query_name alone is not unique: read1/read2 of a pair share it, as do the + primary and supplementary records of a split read. Including the flag and + the start position makes the key identify exactly one record. + """ + return (read.query_name, read.flag, read.reference_start) + + +def classify_read_tag(vcf_dict, source): + """Map one per-read classification onto an IGV-colourable tag string. + + ``source`` is the branch of the read loop that produced ``vcf_dict``: + 'CIGAR', 'SA', 'SOFTCLIP' or 'REF'. + """ + if vcf_dict is None: + return 'Skipped_Unevaluable' + + alttype = vcf_dict.get('alttype') + + if alttype == 'REF': + return 'Unedited_WT' + + # A breakend is named by its partner contig -- that is what the reviewer is + # looking for (e.g. a junction into the transgene contig). + if alttype == 'BND': + return 'Edited_BND_{}'.format(vcf_dict.get('chrom2') or 'NA') + + # Soft-clip-derived calls are realignments rather than direct observations, + # so they are reported by mechanism and not by an implied exact size. + if source == 'SOFTCLIP': + return 'Edited_SoftClip' + + ref = vcf_dict.get('ref') or '' + alt = vcf_dict.get('alt') or '' + + # Symbolic alleles (, , ) carry no length; fall back to the + # reference span between the two breakpoints. + if alt.startswith('<') or ref.startswith('<'): + try: + size = abs(int(vcf_dict['pos2']) - int(vcf_dict['pos'])) + except (KeyError, TypeError, ValueError): + size = 0 + else: + size = abs(len(ref) - len(alt)) + + if alttype == 'DEL': + return 'Edited_Deletion_{}bp'.format(size) + if alttype == 'INS': + return 'Edited_Insertion_{}bp'.format(size) + if alttype == 'DUP': + return 'Edited_Duplication_{}bp'.format(size) + + return 'Edited_Complex' + + +def record_read_tag(read_tags, read, tag): + """Keep the highest-precedence tag seen for this record across all targets.""" + key = read_tag_key(read) + previous = read_tags.get(key) + if previous is None or read_tag_rank(tag) < read_tag_rank(previous): + read_tags[key] = tag + + +def merge_windows(windows): + """Collapse (chrom, start, end) windows into disjoint, sorted intervals. + + Padded target windows overlap each other; fetching them as-is would visit + some records twice. Merging first keeps the second pass close to a single + linear sweep and makes the output naturally near-coordinate-order. + """ + by_chrom = defaultdict(list) + for chrom, start, end in windows: + by_chrom[chrom].append((max(0, int(start)), int(end))) + + merged = [] + for chrom in sorted(by_chrom): + current_start, current_end = None, None + for start, end in sorted(by_chrom[chrom]): + if current_end is not None and start <= current_end: + current_end = max(current_end, end) + continue + if current_end is not None: + merged.append((chrom, current_start, current_end)) + current_start, current_end = start, end + if current_end is not None: + merged.append((chrom, current_start, current_end)) + + return merged + + +def write_tagged_bam(source_bamfile, out_path, windows, read_tags, tag_name, + probabilities=None, prob_tag='XP', verbose=False): + """Second pass: emit a sorted, indexed BAM of the target windows, tagged. + + Restricted to the target windows on purpose. Tagging whole CRAMs would be + enormous, and the windows are all IGV needs to show the reviewer why a site + was called. Returns the number of records written. + """ + merged = merge_windows(windows) + if not merged: + print("No target windows were processed; not writing a tagged BAM.", file=sys.stderr) + return 0 + + if out_path.endswith('.bam'): + final_path = out_path + else: + final_path = out_path + '.bam' + unsorted_path = final_path[:-4] + '.unsorted.bam' + + written = 0 + # Merged windows are disjoint and ascending, so the only record that can be + # fetched twice is one that runs past the end of a window into the next. + # Carrying just those keys (key -> reference_end) is enough to write each + # record exactly once, and costs far less than remembering every key written. + carry = {} + current_chrom = None + + out_bam = pysam.AlignmentFile(unsorted_path, 'wb', template=source_bamfile) + try: + for chrom, start, end in merged: + if chrom != current_chrom: + carry.clear() + current_chrom = chrom + elif carry: + # anything ending at or before this window cannot re-appear + carry = {k: e for k, e in carry.items() if e > start} + + for read in source_bamfile.fetch(chrom, start, end, multiple_iterators=True): + key = read_tag_key(read) + if key in carry: + # already emitted from the previous window -- writing it again + # would show up in IGV as doubled depth + continue + + read.set_tag(tag_name, read_tags.get(key, DEFAULT_TAG), value_type='Z') + if probabilities is not None and key in probabilities: + read.set_tag(prob_tag, float(probabilities[key]), value_type='f') + out_bam.write(read) + written += 1 + + read_end = read.reference_end + if read_end is not None and read_end > end: + carry[key] = read_end + finally: + out_bam.close() + + # Windows are visited in target order, which is not guaranteed to be + # coordinate order, and IGV needs coordinate-sorted + indexed input. + pysam.sort('-o', final_path, unsorted_path) + pysam.index(final_path) + os.remove(unsorted_path) + + if verbose: + print(f"Wrote {written} tagged reads over {len(merged)} merged windows to {final_path}", + file=sys.stderr) + + return written + + # ============================================================================ # SECTION 3: MAIN FUNCTION # ============================================================================ @@ -1515,6 +1750,16 @@ def main(): parser.add_argument('-o','--outfile',type=str,help='Output file (optional)') parser.add_argument('-u','--unevaluable-reads-logfile',type=str,help='File with information on reads that were not evaluable.') parser.add_argument('-V','--vcf-out',type=str,help="VCF output file with all passing events.") + parser.add_argument('--tagged-bam-out',type=str,default=None, + help=('Write a coordinate-sorted, indexed BAM in which every read carries a string tag ' + 'recording how this script classified it (Edited_Deletion_5bp, Unedited_WT, ...), ' + 'for review in IGV via Color alignments by -> tag. Off by default. Output is ' + 'restricted to the target windows only (target +/- --target-window), not the whole ' + 'genome: expect roughly (number of targets) x (2 x window + read length) x depth ' + 'bytes, i.e. a few hundred MB for a typical deep ECS panel, but it scales with your ' + 'target count -- check the size before enabling it across a cohort.')) + parser.add_argument('--tagged-bam-tag',type=str,default='XC', + help='Two-character tag name to write the per-read classification into (default: XC)') parser.add_argument('-v','--verbose',action='store_true',help='Print verbose output') # add vv option for debugging parser.add_argument('-vv','--debug',action='store_true',help='Print debug output') @@ -1665,10 +1910,17 @@ def main(): region_list = [parse_region_string(region) for region in args.regions.split(',')] if args.regions else None total_intervals = len(mergedBedDf) - + # This stores all indel records to print as a VCF at the end. all_indel_records = [] + # Optional read-level tagging (see SECTION 2b). Tags are only accumulated + # here; the BAM is written afterwards in a second pass, because this loop + # visits overlapping targets and would otherwise emit duplicate records. + read_tags = {} if args.tagged_bam_out else None + read_probs = {} if (args.tagged_bam_out and args.enable_crispr_prediction) else None + tagged_windows = [] + # Use enumerate(..., start=1) to keep track of the current loop index for i, (_, row) in enumerate(mergedBedDf.iterrows(), 1): @@ -1707,16 +1959,36 @@ def main(): if args.verbose: print(f"\tGetting reads that align within window of {row['Chromosome']}:{row['Start']}-{row['End']}", file=sys.stderr) + window_start = max(0, row['Start'] - args.target_window) + window_end = row['End'] + args.target_window + if read_tags is not None: + tagged_windows.append((row['Chromosome'], window_start, window_end)) + # get reads that align within a defined region containing the merged target interval - for read in edited_bamfile.fetch(row['Chromosome'], max(0, row['Start']-args.target_window), row['End']+args.target_window, multiple_iterators = True): - - # skip if not primary alignment or a duplicate or poor mapping quality - if (read.is_mapped is False or - read.is_duplicate is True or - read.is_secondary is True or - read.is_supplementary is True or - read.mapping_quality < args.min_mapqual or - count_mismatches_fast(read) > args.max_read_mismatches): + for read in edited_bamfile.fetch(row['Chromosome'], window_start, window_end, multiple_iterators = True): + + # skip if not primary alignment or a duplicate or poor mapping quality. + # Same predicates and same short-circuit order as a single boolean + # chain; split out so a skipped read can say which filter caught it. + skip_reason = None + if read.is_mapped is False: + skip_reason = 'Skipped_Unmapped' + elif read.is_duplicate is True: + skip_reason = 'Skipped_Duplicate' + elif read.is_secondary is True: + skip_reason = 'Skipped_Secondary' + elif read.is_supplementary is True: + skip_reason = 'Skipped_Supplementary' + elif read.mapping_quality < args.min_mapqual: + skip_reason = 'Skipped_LowMapQ' + elif count_mismatches_fast(read) > args.max_read_mismatches: + skip_reason = 'Skipped_Mismatches' + + if skip_reason is not None: + # tag it rather than dropping it, so a reviewer sees that the + # read was excluded on purpose instead of just missing + if read_tags is not None: + record_read_tag(read_tags, read, skip_reason) continue # Determine whether the read pair has the proper orientation, that is: ---> <--- @@ -1734,9 +2006,12 @@ def main(): # dict to store mutation info in VCF record format vcf_dict = None + # which classifier produced vcf_dict; only used for read-level tags + vcf_source = None # if cigar has any D/I operations if any(op in (1, 2) for op, _ in cigar) and proper_paired_read: + vcf_source = 'CIGAR' if args.verbose: print(f"\tAnalyzing cigars in {read.query_name} from {row['Chromosome']}:{row['Start']-args.target_window}-{row['End']+args.target_window}", file=sys.stderr) @@ -1744,16 +2019,20 @@ def main(): vcf_dict = get_cigar_indel_vcf(read, refFasta, row['Pos'], target_index) # if no indel is found or its too far away from the closest PAM position - if (vcf_dict is None or + if (vcf_dict is None or vcf_dict['distance'] > args.max_mutation_distance and vcf_dict['distance2'] > args.max_mutation_distance): # print abbreviated read info (name, cigar, mapping info, sequence) to unevaluable read log if unevaluable_read_log: print(f"{str(read)}\t{vcf_dict}", file=unevaluable_read_log) + if read_tags is not None: + record_read_tag(read_tags, read, 'Skipped_Unevaluable') + continue # read has supplementary alignments - elif read.has_tag('SA'): + elif read.has_tag('SA'): + vcf_source = 'SA' if args.verbose: print(f"\tAnalyzing SA in {read.query_name}, {read.get_tag('SA') if read.has_tag('SA') else 'None'} from {row['Chromosome']}:{row['Start']-args.target_window}-{row['End']+args.target_window}", file=sys.stderr) @@ -1769,8 +2048,11 @@ def main(): if unevaluable_read_log: print(f"{str(read)}\t{vcf_dict}", file=unevaluable_read_log) + if read_tags is not None: + record_read_tag(read_tags, read, 'Skipped_Unevaluable') + continue - + # if SA is an indel then it should be bounded by the read pair ends. If not then continue. if (vcf_dict['alttype'] in ['DEL','DUP','INS'] and proper_paired_read and (min([vcf_dict['pos'],vcf_dict['pos2']]) < min([read.reference_start,read.next_reference_start]) or @@ -1780,8 +2062,11 @@ def main(): if unevaluable_read_log: print(f"{str(read)}\t{vcf_dict}", file=unevaluable_read_log) + if read_tags is not None: + record_read_tag(read_tags, read, 'Skipped_Unevaluable') + continue - + # if SA is a BND, check to see if the other end is in the target list if (vcf_dict['alttype']=='BND' and (read.mapping_quality < args.min_bnd_mapqual or vcf_dict['info']['SAMAPQ'] < args.min_bnd_mapqual and @@ -1791,28 +2076,36 @@ def main(): if unevaluable_read_log: print(f"{str(read)}\t{vcf_dict}", file=unevaluable_read_log) + if read_tags is not None: + record_read_tag(read_tags, read, 'Skipped_Unevaluable') + continue # read has softclips elif cigar[0][0] >= 4 and cigar[0][1] >= args.min_softclip_length or \ cigar[-1][0] >= 4 and cigar[-1][1] >= args.min_softclip_length: + vcf_source = 'SOFTCLIP' if args.verbose: print(f"\tAnalyzing softclips in {read.query_name}, {read.cigarstring}, from {row['Chromosome']}:{row['Start']-args.target_window}-{row['End']+args.target_window}", file=sys.stderr) vcf_dict = get_softclip_indel_vcf(read, refFasta, row['Pos'], args.mutation_search_window, args.min_softclip_length) - if (vcf_dict is None or - (vcf_dict['distance'] > args.max_mutation_distance and vcf_dict['distance2'] > args.max_mutation_distance)): + if (vcf_dict is None or + (vcf_dict['distance'] > args.max_mutation_distance and vcf_dict['distance2'] > args.max_mutation_distance)): # print abbreviated read info (name, cigar, mapping info, sequence) to unevaluable read log if unevaluable_read_log: print(f"{str(read)}\t{vcf_dict}", file=unevaluable_read_log) + if read_tags is not None: + record_read_tag(read_tags, read, 'Skipped_Unevaluable') + continue - # read spans start and end and has no softclips + # read spans start and end and has no softclips elif read.reference_start < row['End'] and read.reference_end > row['Start'] and cigar[0][0] == 0 and cigar[-1][0] == 0: + vcf_source = 'REF' if args.verbose: print(f"\tThis read is reference {read.query_name}, {read.cigarstring}, from {row['Chromosome']}:{row['Start']-args.target_window}-{row['End']+args.target_window}", file=sys.stderr) @@ -1834,8 +2127,19 @@ def main(): # If read doesnt meet any of the criteria then skip it. else: + # Sits in the padded window but carries no event and does not span + # the target, so it is evidence for neither call. This is the bulk + # of the records in a window, and it is exactly DEFAULT_TAG -- so + # it is deliberately NOT recorded, which keeps the tag map small + # enough to survive a full panel. See DEFAULT_TAG. continue - + + # Record the per-read verdict *before* the alleles below are truncated + # for display, since the truncation overwrites alt with a length label + # and would make the indel size in the tag meaningless. + if read_tags is not None: + record_read_tag(read_tags, read, classify_read_tag(vcf_dict, vcf_source)) + # truncate ref or alt allele for readbility in Excel, etc. if len(vcf_dict['ref']) > 20: vcf_dict['alt'] = f"DEL{len(vcf_dict['ref'])-1}" @@ -1972,7 +2276,8 @@ def main(): crispr_predicted_reads, crispr_prediction_probability = predict_reads_at_position( edited_bamfile, row['Chromosome'], row['Start'], row['End'], positions, - crispr_model, refFasta, is_on_target=ontarget, control_bam=control_bamfile, threshold=args.crispr_threshold + crispr_model, refFasta, is_on_target=ontarget, control_bam=control_bamfile, threshold=args.crispr_threshold, + read_probs=read_probs ) crispr_prediction_fraction = round(crispr_predicted_reads/total_reads, 4) if total_reads > 0 else 0.0 @@ -2010,6 +2315,13 @@ def main(): print("\t".join([str(field) for field in output_fields]), file=fp, flush=True) + # Second pass: write the read-level tagged BAM, if requested. This has to be + # a separate sweep -- the loop above revisits the same records at overlapping + # targets, so writing there would emit duplicates and inflate IGV depth. + if read_tags is not None: + write_tagged_bam(edited_bamfile, args.tagged_bam_out, tagged_windows, read_tags, + args.tagged_bam_tag, probabilities=read_probs, verbose=True) + # Write VCF output if requested. if args.vcf_out: vcf_out_df = pd.concat(all_indel_records, axis=0, ignore_index=True) diff --git a/conf/modules.config b/conf/modules.config index ea319e2..82b5ec3 100644 --- a/conf/modules.config +++ b/conf/modules.config @@ -70,9 +70,13 @@ process { // Python scripts are single-threaded, so they get 1 CPU (the process_medium label gave 4). // ECS truth arm — long-running (~15 min–2.5 h, no checkpoint) but light on RAM/CPU. + // --offtarget_tagged_bam is the exception: it holds a per-read tag map for the whole + // target panel until the second (write) pass, which measured ~4 GB on top of the 1.4 GB + // baseline for the 1149-interval AAVS1 site14 panel at ~11,000x. Give it room rather + // than let an opt-in review artifact OOM the truth arm. withName: 'ECS_INDELS' { cpus = { 1 } - memory = { check_max( 8.GB * task.attempt, 'memory' ) } // was process_highmem = 64 GB + memory = { check_max( (params.offtarget_tagged_bam ? 24.GB : 8.GB) * task.attempt, 'memory' ) } time = { check_max( 8.h * task.attempt, 'time' ) } // longest observed run ~2.5 h } diff --git a/docs/OFFTARGET_WORKFLOW.md b/docs/OFFTARGET_WORKFLOW.md index 070704a..bf1c5f3 100644 --- a/docs/OFFTARGET_WORKFLOW.md +++ b/docs/OFFTARGET_WORKFLOW.md @@ -51,6 +51,63 @@ no IGV needed. Off by default (`offtarget_snapshots = false`) because it renders pile up at the cut site; right (unedited normal) — clean. Pass `--snapshots` to `run_offtarget.sh` to turn this on, so the confirmed off-targets (PLCB2, CNNM3) get the same tumor/normal packet.* +### Read-level tags for IGV (optional) + +`.offtarget_analysis.tsv` gives you `indel_fraction = 0.42` and leaves you to work out by eye +which reads made up the 0.42. Add `--offtarget_tagged_bam true` and the ECS caller also writes +`.tagged.bam` (+ `.bai`), in which **every read carries an `XC` string tag naming how the caller +classified it** — so the pileup shows you its reasoning directly. + +In IGV: load the BAM → right-click the track → **Color alignments by** → **tag** → type `XC`. + +| tag | meaning | +|---|---| +| `Edited_Deletion_bp` | deletion of N bp called from the CIGAR or a split alignment | +| `Edited_Insertion_bp` | insertion of N bp | +| `Edited_Duplication_bp` | tandem duplication from a split alignment | +| `Edited_BND_` | breakend — the read's other end maps to `` (e.g. the transgene contig) | +| `Edited_SoftClip` | event recovered by realigning a soft clip | +| `Unedited_WT` | spans the target cleanly, no event — this is the denominator | +| `Skipped_Duplicate` / `_LowMapQ` / `_Mismatches` / `_Secondary` / `_Supplementary` / `_Unmapped` | excluded by a read filter, shown so you can see *why* it isn't counted | +| `Skipped_Unevaluable` | a classifier ran but found nothing near enough to a PAM to call | +| `Skipped_NoSpan` | inside the padded window but doesn't span the target and carries no event | + +Each alignment appears **exactly once**, so the depth IGV shows is real. Where a read falls in two +overlapping target windows and gets classified differently, the most specific call wins — the table +above is in precedence order, top to bottom. + +**Don't expect the tag counts to equal the TSV columns** — they count different things, and the BAM +is the more literal of the two: + +- tags are per **alignment record**, while `indel_reads` is per **fragment** (the caller collapses + R1/R2 by read name) and is taken *after* the site-level filters. So `Edited_*` records run higher + than `indel_reads` — roughly 2× where both mates cover the cut site. On the AAVS1 site14 + on-target that is 17,396 `Edited_*` records against `indel_reads` = 10,821. +- the BAM covers the whole padded window, so it also holds reads that never entered `total_reads`. + Across the AAVS1 site14 panel that is 65% `Unedited_WT`, 35% `Skipped_*` and 0.7% `Edited_*`; at + a *1 bp* target specifically, the ±150 bp pad means most records are `Skipped_NoSpan` (69% at the + on-target). They are kept deliberately — a pileup cropped to reads spanning a single base is + unreadable in IGV. + +Use the tags to see *which* reads drove a call and why; use the TSV for the number. + +**Off by default, and worth keeping that way for routine runs.** Output is restricted to the target +windows (target ±150 bp) rather than the whole CRAM, but it still scales with target count and +depth. Measured on one AAVS1 site14 ECS sample (1,149 target intervals → 1,145 merged windows, +~11,000× at the on-target): + +| | without | with `--offtarget_tagged_bam` | +|---|---|---| +| output | — | **800 MB** `.bam` + 1.5 MB `.bai`, 21.8 M reads | +| peak RSS | 1.4 GB | **6.3 GB** (the tag map is held until the write pass) | +| wall clock | ~30 min | ~46 min | + +`ECS_INDELS` is given 24 GB instead of 8 GB when the flag is set, so enabling it does not OOM the +truth arm. Still check the size on one sample before turning it on across a cohort — and note the +debug log this pipeline learned that lesson from, `offtarget_ecs_unevaluable_log`, reached +0.1–1 TB/sample and filled the work directory. Restricting the run with `--regions` keeps both +numbers small when you only want to review a handful of loci. + ## What runs depends on the samplesheet The `datatype` column decides: @@ -111,6 +168,7 @@ You hand it one samplesheet with these columns: | `offtarget_germline_max_ctrl_if` | 0.05 | matched-normal indel frac above this = germline/artifact, not a somatic edit (`label` 0) | | `offtarget_hotspot_pad` | 25 | bp window to match a worklist hit to a predicted hotspot | | `offtarget_snapshots` | false | render IGV-style pileup PNGs for LIKELY EDITs | +| `offtarget_tagged_bam` | false | emit `.tagged.bam` with per-read `XC` tags for IGV review (above) | | `offtarget_rescue` | true | high-evidence rescue (below); `false` = call on model score alone | | `offtarget_rescue_min_ifrac` | 0.15 | rescue: minimum indel fraction in the edited sample | | `offtarget_rescue_min_conc` | 0.5 | rescue: minimum positional concordance (clonality) | diff --git a/modules/local/ecs_indels.nf b/modules/local/ecs_indels.nf index ebc27fc..e8aa716 100644 --- a/modules/local/ecs_indels.nf +++ b/modules/local/ecs_indels.nf @@ -13,6 +13,7 @@ process ECS_INDELS { output: tuple val(meta), path("${meta.id}.offtarget_analysis.tsv"), emit: indels_file tuple val(meta), path("${meta.id}.offtarget_edits.vcf"), emit: indels_vcf + tuple val(meta), path("${meta.id}.tagged.bam*"), optional: true, emit: tagged_bam path "versions.yml", emit: versions script: @@ -20,6 +21,10 @@ process ECS_INDELS { // output and is not consumed downstream, yet it grows to ~0.1-1 TB per sample and // was the sole cause of multi-TB work-dir bloat / ENOSPC. Off unless explicitly asked. def unevaluable = params.offtarget_ecs_unevaluable_log ? "-u ${meta.id}.unevaluable_reads.txt" : "" + // Review aid, not a pipeline input: a BAM of the target windows in which every read + // carries an XC tag naming the per-read call, for colouring the pileup in IGV. Window- + // restricted, but it still scales with target count -- off unless explicitly asked. + def tagged_bam = params.offtarget_tagged_bam ? "--tagged-bam-out ${meta.id}.tagged.bam" : "" """ python ${projectDir}/bin/find_edited_reads.py \\ --fasta ${reference} \\ @@ -27,6 +32,7 @@ process ECS_INDELS { --control-bam ${control_cram} \\ --target-file ${target_file} \\ ${unevaluable} \\ + ${tagged_bam} \\ --vcf-out ${meta.id}.offtarget_edits.vcf \\ -o ${meta.id}.offtarget_analysis.tsv @@ -37,5 +43,6 @@ process ECS_INDELS { """ stub: - "touch ${meta.id}.offtarget_analysis.tsv ${meta.id}.offtarget_edits.vcf versions.yml" + def tagged_bam = params.offtarget_tagged_bam ? "${meta.id}.tagged.bam ${meta.id}.tagged.bam.bai" : "" + "touch ${meta.id}.offtarget_analysis.tsv ${meta.id}.offtarget_edits.vcf ${tagged_bam} versions.yml" } diff --git a/modules/local/get_indels.nf b/modules/local/get_indels.nf index c6557f7..6c9db7c 100644 --- a/modules/local/get_indels.nf +++ b/modules/local/get_indels.nf @@ -12,6 +12,7 @@ process GET_INDELS { output: tuple val(meta), path("${meta.id}.offtarget_analysis.tsv"), emit: indels_file tuple val(meta), path("${meta.id}.offtarget_edits.vcf"), emit: indels_vcf + tuple val(meta), path("${meta.id}.tagged.bam*"), optional: true, emit: tagged_bam //tuple val(meta), path("${meta.id}.ml_results.txt"), emit: ml_results //tuple val(meta), path("${meta.id}.fp_filtered.txt"), emit: fp_log path "versions.yml", emit: versions @@ -25,8 +26,12 @@ process GET_INDELS { hotspot_file ? "--target-file ${hotspot_file}" : "" ].join(' ').trim() - """ - find_edited_reads.py ${inputs} -u ${meta.id}.unevaluable_reads.txt --vcf-out ${meta.id}.offtarget_edits.vcf -o ${meta.id}.offtarget_analysis.tsv + // Optional review aid: a window-restricted BAM whose reads carry an XC tag naming + // the per-read call, for colouring the pileup in IGV. Off unless explicitly asked. + def tagged_bam = params.offtarget_tagged_bam ? "--tagged-bam-out ${meta.id}.tagged.bam" : "" + + """ + find_edited_reads.py ${inputs} ${tagged_bam} -u ${meta.id}.unevaluable_reads.txt --vcf-out ${meta.id}.offtarget_edits.vcf -o ${meta.id}.offtarget_analysis.tsv cat <<-END_VERSIONS > versions.yml ${task.process}: @@ -35,9 +40,11 @@ process GET_INDELS { """ stub: + def tagged_bam = params.offtarget_tagged_bam ? "touch ${meta.id}.tagged.bam ${meta.id}.tagged.bam.bai" : "" """ touch ${meta.id}.offtarget_analysis.tsv touch ${meta.id}.offtarget_edits.vcf + ${tagged_bam} cat <<-END_VERSIONS > versions.yml ${task.process}: diff --git a/nextflow.config b/nextflow.config index fc05047..7a6a701 100644 --- a/nextflow.config +++ b/nextflow.config @@ -59,6 +59,7 @@ params { offtarget_top = 60 // top-N to print / snapshot offtarget_snapshots = false // render IGV-style pileup PNGs for LIKELY EDITs offtarget_ecs_unevaluable_log = false // write per-read unevaluable-reads debug log (very large: ~0.1-1 TB/sample); off by default + offtarget_tagged_bam = false // emit a target-window BAM whose reads carry an XC tag with the per-read call, for IGV review; off by default offtarget_ecs_edit_threshold = 0.0 // ECS indel_fraction strictly above this = an ECS edit (ecs_is_edit) offtarget_germline_max_ctrl_if = 0.05 // matched-normal (WGS) indel frac above this = germline/artifact, not a somatic edit (label 0) offtarget_hi_score = 0.60 // WGS shape score >= this = detected (recall curve) diff --git a/nextflow_schema.json b/nextflow_schema.json index 0ae44ef..7b8599f 100644 --- a/nextflow_schema.json +++ b/nextflow_schema.json @@ -166,6 +166,11 @@ "help_text": "Off by default. This log reaches roughly 0.1-1 TB per sample and was the sole cause of work-directory ENOSPC failures; enable only when debugging read evaluation.", "hidden": true }, + "offtarget_tagged_bam": { + "type": "boolean", + "description": "Emit a per-read tagged BAM from the ECS caller for review in IGV", + "help_text": "Off by default. Writes .tagged.bam(.bai) in which every read carries an XC string tag naming how the caller classified it (Edited_Deletion_5bp, Unedited_WT, Skipped_LowMapQ, ...), so a reviewer can load it in IGV and use Color alignments by -> tag -> XC. Output is restricted to the target windows (target +/- 150 bp), not the whole genome, but it still scales with target count and depth: one AAVS1 site14 ECS sample (1149 targets, ~11000x) produced a 800 MB BAM of 21.8 M reads and pushed ECS_INDELS from 1.4 GB to 6.3 GB peak RSS, so the process is given 24 GB instead of 8 GB while this is set. Check the size on one sample before enabling it across a cohort. See docs/OFFTARGET_WORKFLOW.md." + }, "offtarget_ecs_edit_threshold": { "type": "number", "default": 0.0, diff --git a/tests/conftest.py b/tests/conftest.py index 9416b1b..f5c7c92 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -80,3 +80,174 @@ def workspace(tmp_path): ], columns=["rank", "sample", "chrom", "start", "alt", "dragen_af", "score", "verdict_pon"] ).to_csv(d / "worklist_pon.csv", index=False) return d + + +# --------------------------------------------------------------------------- +# Synthetic aligned-read workspace for find_edited_reads.py +# +# The glue fixtures above are pure tables; the read-level tagging path needs +# real alignments. This builds the smallest workspace that reproduces the two +# things that make tagging non-trivial: +# +# * two targets 200 bp apart -- far enough that pyranges clusters them +# separately (slack = --target-window = 150), close enough that their +# +/-150 fetch windows overlap, so some reads are visited twice; +# * a read pair sitting in that overlap which one target calls Unedited_WT +# and the other cannot evaluate, i.e. the tag-conflict case. +# --------------------------------------------------------------------------- + +ECS_CHROM = "chr1" +ECS_CONTIG_LEN = 3000 +ECS_TARGET_A = 1001 # 1-based; the edited site +ECS_TARGET_B = 1201 # 1-based; a quiet site +ECS_DEL_LEN = 5 +ECS_N_EDIT = 12 # read pairs carrying the deletion +ECS_N_WT = 8 # read pairs spanning target A cleanly +ECS_N_DUP = 2 # duplicate-flagged pairs +ECS_N_LOWMAPQ = 2 # MAPQ below the default floor of 20 +ECS_N_MISMATCH = 2 # NM above the default ceiling of 4 +ECS_N_OVERLAP = 6 # pairs inside the window overlap, spanning neither target +ECS_N_SPAN_B = 3 # pairs in the overlap that DO span target B (conflict case) + + +def _write_ecs_fasta(path): + import random + random.seed(7) + seq = "".join(random.choice("ACGT") for _ in range(ECS_CONTIG_LEN)) + with open(path, "w") as fh: + fh.write(f">{ECS_CHROM}\n") + for i in range(0, ECS_CONTIG_LEN, 60): + fh.write(seq[i:i + 60] + "\n") + import pysam + pysam.faidx(str(path)) + return pysam.FastaFile(str(path)) + + +def _ecs_pair(fasta, name, r1_start, r1_cigar, r2_start, mapq=60, dup=False, + r1_nm=0, r2_nm=0): + """One properly-paired FR pair (R1 forward, R2 reverse) with real sequence.""" + import pysam + + def qseq(start, cigar): + out, ref = [], start + for op, ln in cigar: + if op == 0: # M + out.append(fasta.fetch(ECS_CHROM, ref, ref + ln)); ref += ln + elif op == 2: # D + ref += ln + elif op == 1: # I + out.append("A" * ln) + elif op == 4: # S + out.append("T" * ln) + return "".join(out) + + span = r2_start + 100 - r1_start + reads = [] + for is_read1, start, cigar, mate_start, nm in ( + (True, r1_start, r1_cigar, r2_start, r1_nm), + (False, r2_start, [(0, 100)], r1_start, r2_nm)): + a = pysam.AlignedSegment() + a.query_name = name + a.query_sequence = qseq(start, cigar) + a.reference_id = 0 + a.reference_start = start + a.mapping_quality = mapq + a.cigar = cigar + a.next_reference_id = 0 + a.next_reference_start = mate_start + a.template_length = span if is_read1 else -span + a.query_qualities = pysam.qualitystring_to_array("I" * len(a.query_sequence)) + a.is_paired = True + a.is_proper_pair = True + a.is_read1 = is_read1 + a.is_read2 = not is_read1 + a.is_reverse = not is_read1 + a.mate_is_reverse = is_read1 + a.is_duplicate = dup + a.set_tag("NM", nm, value_type="i") + a.set_tag("MC", "100M", value_type="Z") + reads.append(a) + return reads + + +def _build_ecs_cram(fasta, fasta_path, out_bam, reads): + import pysam + header = {"HD": {"VN": "1.6", "SO": "coordinate"}, + "SQ": [{"SN": ECS_CHROM, "LN": ECS_CONTIG_LEN}]} + reads.sort(key=lambda r: r.reference_start) + tmp = str(out_bam) + ".tmp.bam" + with pysam.AlignmentFile(tmp, "wb", header=header) as out: + for r in reads: + out.write(r) + pysam.sort("-o", str(out_bam), tmp) + Path(tmp).unlink() + pysam.index(str(out_bam)) + # production opens its inputs as CRAM, so hand the caller a CRAM + cram = str(out_bam)[:-4] + ".cram" + with pysam.AlignmentFile(str(out_bam)) as src, \ + pysam.AlignmentFile(cram, "wc", template=src, + reference_filename=str(fasta_path)) as out: + for r in src: + out.write(r) + pysam.index(cram) + return cram + + +@pytest.fixture +def ecs_reads_workspace(tmp_path): + """Synthetic reference + targets VCF + edited/control CRAMs for the ECS caller. + + Returns a dict of paths plus the read counts the tags should reconcile with. + """ + pysam = pytest.importorskip("pysam") + d = tmp_path + fasta_path = d / "ref.fa" + fasta = _write_ecs_fasta(fasta_path) + + header = pysam.VariantHeader() + header.contigs.add(ECS_CHROM, length=ECS_CONTIG_LEN) + targets = d / "targets.vcf" + with pysam.VariantFile(str(targets), "w", header=header) as vout: + for pos in (ECS_TARGET_A, ECS_TARGET_B): + rec = vout.new_record() + rec.chrom = ECS_CHROM + rec.pos = pos + rec.id = "." + rec.ref = fasta.fetch(ECS_CHROM, pos - 1, pos) + rec.alts = ("N",) + rec.filter.add("PASS") + vout.write(rec) + + del_cigar = [(0, 60), (2, ECS_DEL_LEN), (0, 40)] + edited = [] + for i in range(ECS_N_EDIT): + edited += _ecs_pair(fasta, f"edit{i}", 940, del_cigar, 1180, r1_nm=ECS_DEL_LEN) + for i in range(ECS_N_WT): + edited += _ecs_pair(fasta, f"wt{i}", 940, [(0, 100)], 1180) + for i in range(ECS_N_DUP): + edited += _ecs_pair(fasta, f"dup{i}", 940, [(0, 100)], 1180, dup=True) + for i in range(ECS_N_LOWMAPQ): + edited += _ecs_pair(fasta, f"lowq{i}", 940, [(0, 100)], 1180, mapq=3) + for i in range(ECS_N_MISMATCH): + edited += _ecs_pair(fasta, f"mm{i}", 940, [(0, 100)], 1180, r1_nm=6, r2_nm=6) + for i in range(ECS_N_OVERLAP): + edited += _ecs_pair(fasta, f"both{i}", 1080, [(0, 100)], 1220) + # R1 lands in both padded windows: target A cannot evaluate it, target B + # calls it reference. Unedited_WT must win, and it must be written once. + for i in range(ECS_N_SPAN_B): + edited += _ecs_pair(fasta, f"spanb{i}", 1150, [(0, 100)], 1220) + + control = [] + for i in range(ECS_N_EDIT + ECS_N_WT): + control += _ecs_pair(fasta, f"ctl{i}", 940, [(0, 100)], 1180) + + return { + "dir": d, + "fasta": fasta_path, + "targets": targets, + "edited": _build_ecs_cram(fasta, fasta_path, d / "edited.bam", edited), + "control": _build_ecs_cram(fasta, fasta_path, d / "control.bam", control), + "n_edit_reads": ECS_N_EDIT, + "del_len": ECS_DEL_LEN, + "total_records": len(edited), + } diff --git a/tests/test_offtarget_glue.py b/tests/test_offtarget_glue.py index ed34b5f..af47e78 100644 --- a/tests/test_offtarget_glue.py +++ b/tests/test_offtarget_glue.py @@ -8,6 +8,7 @@ import sys import pandas as pd +import pytest from conftest import BIN, run @@ -163,3 +164,148 @@ def test_verdict_callers_unpack_three_values(): offenders.append(f"{path.name}:{node.lineno} unpacks {len(tgt.elts)}, want 3") assert not offenders, "verdict() arity mismatch: " + "; ".join(offenders) + + +# -------------------------------------------------------------------------- +# find_edited_reads.py --tagged-bam-out (read-level tags for IGV review) +# +# The caller visits the same alignment record once per overlapping target, so +# the risk here is a BAM with duplicate records, which IGV silently renders as +# doubled depth. These run the real caller on the synthetic CRAM workspace. +# -------------------------------------------------------------------------- +import collections + +import conftest as C + + +def _load_find_edited_reads(): + """Import bin/find_edited_reads.py as a module (it is a script, not a package).""" + import importlib.util + spec = importlib.util.spec_from_file_location("find_edited_reads", + BIN / "find_edited_reads.py") + module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(module) + return module + + +def _run_caller(ws, *extra): + run("find_edited_reads.py", + "--fasta", ws["fasta"], "--edited-bam", ws["edited"], + "--control-bam", ws["control"], "--target-file", ws["targets"], + "-o", ws["dir"] / "out.tsv", *extra, cwd=ws["dir"]) + return pd.read_csv(ws["dir"] / "out.tsv", sep="\t") + + +def test_tagged_bam_off_by_default(ecs_reads_workspace): + ws = ecs_reads_workspace + tsv = _run_caller(ws) + assert tsv["indel_reads"].sum() == ws["n_edit_reads"] + # the feature is opt-in: nothing extra on disk unless asked for + assert not list(ws["dir"].glob("*.tagged.bam*")) + + +def test_tagged_bam_tags_every_read_exactly_once(ecs_reads_workspace): + pysam = pytest.importorskip("pysam") + ws = ecs_reads_workspace + out = ws["dir"] / "tagged.bam" + tsv = _run_caller(ws, "--tagged-bam-out", out) + + assert out.exists() and (ws["dir"] / "tagged.bam.bai").exists() + + keys, tags, positions = [], collections.Counter(), [] + spanb_read1 = [] + with pysam.AlignmentFile(str(out)) as bam: + assert bam.header.to_dict()["HD"]["SO"] == "coordinate" + for read in bam: + assert read.has_tag("XC"), f"{read.query_name} has no XC tag" + keys.append((read.query_name, read.flag, read.reference_start)) + tags[read.get_tag("XC")] += 1 + positions.append(read.reference_start) + if read.query_name.startswith("spanb") and read.is_read1: + spanb_read1.append(read.get_tag("XC")) + + # (a) one record per alignment: the whole point of the second pass + duplicated = [k for k, n in collections.Counter(keys).items() if n > 1] + assert not duplicated, f"duplicate records in tagged BAM: {duplicated[:5]}" + assert len(keys) == ws["total_records"] + + # (b) coordinate-sorted, so IGV will load it + assert positions == sorted(positions) + + # (c) tag counts reconcile with the TSV. Tags are per alignment record while + # indel_reads is per fragment (the caller collapses mates by read name), so + # these are only equal because just R1 carries the deletion in this fixture. + # On real data the record count runs ~2x the TSV number. + edit_tag = f"Edited_Deletion_{ws['del_len']}bp" + assert tags[edit_tag] == ws["n_edit_reads"] == tsv["indel_reads"].sum() + + # (d) skipped reads are visibly skipped, not silently absent + assert tags["Skipped_Duplicate"] == 2 * C.ECS_N_DUP + assert tags["Skipped_LowMapQ"] == 2 * C.ECS_N_LOWMAPQ + assert tags["Skipped_Mismatches"] == 2 * C.ECS_N_MISMATCH + # reads sitting in the +/-150 bp pad that span neither target: the implicit + # default, never stored in the tag map (see DEFAULT_TAG) + assert tags["Skipped_NoSpan"] == 2 * C.ECS_N_OVERLAP + C.ECS_N_SPAN_B + + # (e) conflict resolution: these reads fall in both padded windows, where one + # target cannot evaluate them and the other calls them reference + assert spanb_read1 == ["Unedited_WT"] * C.ECS_N_SPAN_B + + +def test_tagged_bam_honours_custom_tag_name(ecs_reads_workspace): + pysam = pytest.importorskip("pysam") + ws = ecs_reads_workspace + out = ws["dir"] / "tagged.bam" + _run_caller(ws, "--tagged-bam-out", out, "--tagged-bam-tag", "YC") + with pysam.AlignmentFile(str(out)) as bam: + read = next(iter(bam)) + assert read.has_tag("YC") and not read.has_tag("XC") + + +def test_read_tag_precedence_is_order_independent(): + """A read seen at two overlapping targets keeps the most specific call.""" + m = _load_find_edited_reads() + + ordered = ["Edited_BND_chr19", "Edited_Deletion_5bp", "Edited_SoftClip", + "Unedited_WT", "Skipped_NoSpan", "Skipped_LowMapQ"] + ranks = [m.read_tag_rank(t) for t in ordered] + assert ranks == sorted(ranks) and len(set(ranks)) == len(ranks) + + class FakeRead: + query_name, flag, reference_start = "r1", 99, 100 + + for first, second in ((("Unedited_WT"), "Edited_Deletion_5bp"), + ("Edited_Deletion_5bp", "Unedited_WT"), + ("Skipped_NoSpan", "Unedited_WT")): + tags = {} + m.record_read_tag(tags, FakeRead(), first) + m.record_read_tag(tags, FakeRead(), second) + winner = min([first, second], key=m.read_tag_rank) + assert tags[("r1", 99, 100)] == winner + + +def test_classify_read_tag_vocabulary(): + m = _load_find_edited_reads() + assert m.classify_read_tag(None, "CIGAR") == "Skipped_Unevaluable" + assert m.classify_read_tag({"alttype": "REF"}, "REF") == "Unedited_WT" + assert m.classify_read_tag( + {"alttype": "BND", "chrom2": "PLVM_CD19_CARv4_cd34"}, "SA" + ) == "Edited_BND_PLVM_CD19_CARv4_cd34" + assert m.classify_read_tag( + {"alttype": "DEL", "ref": "ATTTTT", "alt": "A"}, "CIGAR") == "Edited_Deletion_5bp" + assert m.classify_read_tag( + {"alttype": "INS", "ref": "A", "alt": "ACCC"}, "CIGAR") == "Edited_Insertion_3bp" + # symbolic alleles carry no length, so fall back to the breakpoint span + assert m.classify_read_tag( + {"alttype": "DUP", "ref": "A", "alt": "", "pos": 100, "pos2": 112}, + "SA") == "Edited_Duplication_12bp" + # soft-clip calls are realignments, reported by mechanism not implied size + assert m.classify_read_tag( + {"alttype": "DEL", "ref": "ATTTTT", "alt": "A"}, "SOFTCLIP") == "Edited_SoftClip" + + +def test_merge_windows_collapses_overlapping_targets(): + m = _load_find_edited_reads() + merged = m.merge_windows([("chr1", 850, 1151), ("chr1", 1050, 1351), + ("chr1", 5000, 5300), ("chr2", 10, 300)]) + assert merged == [("chr1", 850, 1351), ("chr1", 5000, 5300), ("chr2", 10, 300)]