badclip

Extract breakends and structural-variant signals from long-read alignments.

badclip reads read-to-reference alignments and reports breakends: the points where a read is soft-clipped away from an alignment (clips) and the junctions between chimeric alignments of the same read (joins). Joins reveal structural rearrangements (translocations, inversions, large indels); clips flag reads that end abruptly against nothing.

It has six subcommands, forming a pipeline:

Every subcommand reads its input from a file or from stdin via -, transparently decompresses gzip'd input, and prints its own --help (instead of waiting on stdin) when run with no input.

cargo build --release
# binary at target/release/badclip

Find breakends (clips and joins) in read-to-reference alignments.

badclip extract [OPTIONS] [INPUT]

• INPUT — an alignment file (SAM/BAM/CRAM, format auto-detected) by default, or PAF with --paf. Pass - to read from stdin (autodetection works on a pipe too). PAF input may be gzip'd (auto-detected). CRAM requires a reference via -r.

Options:

Examples:

badclip extract aln.bam                     # BAM (default), sorted or unsorted
badclip extract aln.cram -r ref.fa          # CRAM (needs faidx-indexed reference)
samtools view -b … | badclip extract -      # BAM/CRAM from stdin (autodetected)
badclip extract --paf aln.paf.gz            # PAF (gzip auto-detected)
minimap2 … | badclip extract --paf -        # PAF from stdin

At the end of a run, extract prints input statistics to stderr: the number of reads (primary records, including unmapped ones; secondary and supplementary records are ignored), the number of unmapped reads (for PAF, lines whose target name is *, as from minimap2 --paf-no-hit), the total bases in primary alignments (the primary hit's aligned query span per read; for PAF, which has no primary flag, the read's longest span), and the read N50 length. Reads are counted before the -a filter.

• Alignment file (SAM/BAM/CRAM): the container format is auto-detected. A read's chimeric hits are read from the primary alignment's SA:Z: tag; secondary/supplementary/unmapped records are ignored. Works on both coordinate-sorted and name-grouped files (no grouping required). Paired-end reads get a /1 or /2 suffix on the read name (first/last segment) so the mates stay distinct; unpaired reads keep the bare name. CRAM needs -r ref.fa (with ref.fa.fai) to decode read sequences; -r is ignored for BAM/SAM.
• PAF: alignments are grouped by read name (as produced directly by minimap2); primary/supplementary (tp:A:P) and inversion (tp:A:I) alignments are used, secondary (tp:A:S) are ignored.
• Within a read, hits are sorted by their start position along the read.

An alignment file and PAF produce the same output for the same alignments, with one caveat: aln_len for supplementary hits is a few bp smaller from BAM/CRAM because minimap2's SA CIGAR is collapsed (the primary hit's aln_len matches).

Tab-delimited, 9 columns:

ctg1   pos1   ori   ctg2   pos2   qname   mapq   strand   source=..;idx=..;n_aln=..;aln_len=..;qlen=..;mapq=..[;elen=..[;equal=..;eseq=..]]

ori is two characters, each > or < (or . for a missing mate). Positions are raw 0-based offsets that sit between bases. The INFO column is a list of key=value tags; their order is not significant (consumers parse by key). It is written as source, then idx, then three tags in output order (index 1 ↔ ctg1:pos1, index 2 ↔ ctg2:pos2):

• source=NAME — the -s dataset label (default foo), for telling reads from different datasets apart (tumor/normal, trio) in downstream tools. Written first only for readability; the same on every record.
• idx=N — 0-based index of this clip/breakend within the read, in emission order (left clip, joins, right clip); resets per read.
• n_aln=N — number of alignments (primary + supplementary) for the read that passed the -q/-a filters. Same on every breakend of the read (e.g. 1 for a single-alignment read, 3 for a 3-way chimera).
• aln_len=len1,len2 — alignment length (PAF column 11) of each hit; a clip has no second hit, so len2 = 0.
• qlen=left,middle,right — query lengths that sum to the read length. For a clip, right is the clipped length and middle is 0; for a join, middle is the (possibly negative) query gap between the two hits.
• mapq=mapq1,mapq2 — mapq of each hit; for a clip mapq2 = 0. (The mapq column above is the larger of a join's two mapqs; this tag keeps both.)

The last two tags come from alignment-file input only (PAF has no read sequence; for CRAM the sequence is reconstructed against -r), in read-forward orientation (not flipped with the output):

• elen=leftFlank,qdist,rightFlank — sizes of the extracted window: up to -f bases either side of the junction (qdist is the query gap, 0 for a clip), so |eseq| = leftFlank + |qdist| + rightFlank.

• equal=<phred> — the quality of the eseq window: each base quality is turned into an error rate, the expected errors are averaged over the window and scaled back to phred (-10·log10((Σ 10^(-Q/10))/|eseq|)). Sits between elen and eseq, emitted only when eseq is (and when the read has base qualities). For a constant quality Q it is Q.

• eseq=<bases> — the read sequence over that window, on the original read strand. Omitted (while elen stays) when the window exceeds -e.

• Clip — one side is a real breakend, the other is nothing (.):

  chr1   123569841   >.   .   .   read/ccs   1   -   source=foo;idx=0;aln_len=14732,0;qlen=14925,0,141;mapq=1,0
  chr1   123555134   <.   .   .   read/ccs   1   -   source=foo;idx=1;aln_len=14732,0;qlen=14774,0,292;mapq=1,0
  

• Join — two mapped ends meet. The smaller (contig, position) is written first; the two arrows show how the sides are oriented:

  chr1    57375269   >>   chr21   32069271   read/ccs   60   -   source=foo;idx=0;aln_len=16879,26284;qlen=16808,1,26163;mapq=60,60
  chr13   51911798   <>   chr2    76026062   read/ccs   60   +   source=foo;idx=0;aln_len=29505,5668;qlen=29436,1,5661;mapq=60,60
  

  chr1:57375269 >> chr21:32069271 means the right side of chr1:57375269 joins the left side of chr21:32069271. chr13:51911798 <> chr2:76026062 means the left side of chr13:51911798 joins the left side of chr2:76026062.

Convert extract output into a FASTA of the extracted breakend sequences.

badclip geteseq [OPTIONS] [INPUT]

• INPUT — extract output (gzip ok; - or omit for stdin).

Options:

For each record that has an eseq tag, it writes one FASTA entry named readName_idx_leftFlank_rightFlank (from the idx and elen tags) with the eseq bases; records without eseq are skipped.

badclip extract aln.bam | badclip geteseq - > eseq.fa   # "-" = stdin
badclip geteseq extract.txt                             # or a file

Example output:

>m84039_.../234884533/ccs_0_250_250
ACTTTGGGAGGCCAAGGCAGGCGGATCACCTG...

Keep only breakends whose junction sequence is not already explained by a pangenome. Align the geteseq FASTA against a pangenome index (e.g. ropebwt3 sw) to get a PAF, then:

badclip flteseq [OPTIONS] <extractOut> <ropebwt3.paf>

• <extractOut> — extract output (gzip ok; - for stdin).
• <ropebwt3.paf> — the ropebwt3 sw PAF of the geteseq FASTA (gzip ok).

Options:

A record is dropped if it has no eseq tag, or if some ropebwt3 alignment's query interval [qs,qe] on the eseq contains the junction interval [max(0, elen0-l), min(elen0+|elen1|+l, |eseq|)] (i.e. the alignment spans the breakend, with -l margin on each side — the breakend is "protected"). Surviving extract lines are printed verbatim.

With -s STR, the output instead contains every input line, with the survivors' source= tag rewritten to STR (dropped lines keep their original source). This marks novel breakends distinctly so a downstream merge can tell them from pangenome-explained ones.

The PAF's query names must be the geteseq names, and in the same order as the eseq records (as produced by the pipeline), so the two files are streamed without loading into memory.

badclip extract aln.bam > clip.txt
badclip geteseq clip.txt | ropebwt3 sw -d pangenome.fmd - > rb3.paf
badclip flteseq clip.txt rb3.paf > novel.txt

Cluster per-read breakends from extract into consensus SV calls, each with a supporting-read count. Unlike the streaming subcommands above, merge loads its input into memory and sorts it by (ctg1, pos1), so the input need not be pre-sorted.

badclip merge [OPTIONS] [INPUT]...

• INPUT — extract output (gzip ok; - for stdin). Several files are merged as if concatenated, e.g. badclip merge a.clip.gz <(badclip flteseq ...).

Options:

Breakends within -w bp on both endpoints (with ></<> treated as the same inversion) are grouped; a group is reported only if it has at least -c reads with at least -s on each strand. Output keeps the 9-column, TAB-delimited extract shape, with the supporting-read count in the mapq slot and a merge-derived INFO column:

ctg1   pos1   ori   ctg2   pos2   .   count   strand   avg_mapq=q1,q2;count=..[;count_fr=..;count_rf=..[;foldback]];reads=..

• avg_mapq=q1,q2 — mean per-hit mapq of the cluster's reads, split by output endpoint: q1 on the ctg1:pos1 side, q2 on the ctg2:pos2 side (q2 = 0 for a clip cluster).
• count=src:fwd,rev|… — for each source= dataset present in the cluster (from the reads' source= tag), its forward,reverse read counts; |-joined and sorted alphabetically by source. Only sources actually present are listed (never x:0,0), e.g. count=foo:2,3|retain:4,5, or count=foo:2,2 for a single source. The column-7 count is the grand total across all sources and strands.
• count_fr, count_rf — reads oriented >< and <>; present only when the cluster has at least one such read. foldback is added when only one of the two orientations is present (count_fr or count_rf is 0) and both endpoints are on the same contig.
• reads=src:name,…|… — the supporting read names, stratified by source= the same way as count=: |-joined and alphabetical by source (names in member order within each source), e.g. reads=foo:r1,r2|retain:r3,r4. Lists the same sources, in the same order, as count=.

badclip extract aln.bam | badclip merge - > sv.txt

With -m, merge re-runs the same clustering on the (filtered) combination of several samples' merge outputs — a second-level, population merge. Only the INFO parsing and emission change; the clustering is identical. mapq is read from avg_mapq=q1,q2 and per-strand read counts from the count= tag (summed over sources); the -c/-s cluster filters and -q/-p stay read-based (-Q is a no-op — merge output has no equal tag). The per-line -C/-S filters drop under-supported input lines before clustering.

Output stays in the 9-column shape but with column 7 = the number of distinct sources (samples) in the cluster and a compact INFO: avg_mapq=q1,q2;count=F,R (aggregate forward,reverse reads, no per-source breakdown) plus count_fr/count_rf/foldback for inversions — no reads= tag. The output is itself -m-parseable, so merges chain.

# combine per-sample merge outputs, then re-merge across samples
cat sample*.merge.txt | badclip merge -m - > population.txt

Drop extract/merge breakends that fall in BED regions (e.g. blacklists, centromeres, low-complexity regions).

badclip fltreg [-l INT] [INPUT] [BED]

• INPUT — extract or merge output (gzip ok; - or omit for stdin).
• BED — BED file of regions to filter against (gzip ok).
• -l INT — pad each BED interval to [start-l, end+l) before testing (default 0, i.e. the raw intervals).

Both output formats share the 9-column layout, so one filter serves both. A line is dropped if either breakend endpoint — (ctg1, pos1) (columns 1,2) or, when present, (ctg2, pos2) (columns 4,5) — falls in a region; a clip (ctg2 = ".") tests only its one endpoint. Survivors are printed verbatim; a line whose endpoints can't be parsed is kept.

Positions are raw 0-based offsets and BED intervals are half-open, so an offset p is inside [start, end) iff start <= p < end (or start-l <= p < end+l with -l). The BED is loaded into per-contig intervals, padded by -l, sorted by start and merged into a disjoint list, so each lookup is a single binary search (no new dependency; overlapping BED input is handled).

badclip merge - < sv.txt | badclip fltreg - blacklist.bed.gz > sv.filtered.txt

Keep merge calls with enough read support from selected sources — by default the flteseq-filtered (*.flt) ones.

badclip fltcnt [OPTIONS] [INPUT]

• INPUT — merge output (gzip ok; - or omit for stdin).
• -c INT — minimum reads from the --src sources (default 3).
• -s INT — minimum --src reads on each strand (default 0).
• -r INT — maximum reads from the --rest sources (default 0, i.e. every read must come from --src; -1 = no limit).
• -l NUM — minimum pos2 - pos1 (columns 5 and 2). Any value >= 0 also drops clips (column 4 = .), so the default 0 drops clips only; -1 disables the filter. Inter-contig joins have no distance and pass. Accepts a k/m/g suffix, e.g. -l 100k, -l 1m.
• --src LIST — comma-separated sources tested by -c/-s (default: every source ending with .flt).
• --rest LIST — comma-separated sources tested by -r (default: every source not in --src).

Read counts come from the count= INFO tag (count=S.flt:f,r|S.raw:f,r); each group's counts are summed over its sources. The typical input merges the raw breakends (e.g. S.raw) with the flteseq -s S.flt survivors, so by default -c/-s test the novel S.flt reads and -r 0 drops any call that also has reads from another source (S.raw, or other samples in a multi-sample merge). A line with no source in either group (e.g. merge -m output, whose count= carries no source names) aborts the run, and a warning is printed if no line has a --src source at all (the *.flt default doesn't fit — name the sources with --src). Survivors are printed verbatim.

badclip extract -s S.raw in.bam > S.clip
badclip flteseq -s S.flt S.clip S.rb3.paf | badclip merge - | badclip fltcnt - > S.sv

All six subcommands are implemented: extract (SAM/BAM/CRAM auto-detected, or PAF via --paf) finds breakends; geteseq turns its output into FASTA; flteseq filters breakends against a pangenome; merge clusters them into consensus SV calls (with -m for a second-level, cross-sample merge); fltreg drops breakends in BED regions; and fltcnt keeps calls with enough raw/filtered support. See CLAUDE.md for internals and the interval/offset convention.