Reject soft clipping that comes from only one read strand

[jeffreybarrick]

Context

SC evidence has been unusable on real data. Across the 29 LTEE Ara-2 clones in claude-clone-curation, breseq accepted 1040 SC positions — 5 to 620 per sample — and nearly all of them are one artifact: Illumina dark-cycle poly-G read tails.

• 929/1040 (89%) had a homopolymer run ≥8 in the 12-base consensus tail; 802 were pure GGGGGGGGGGGG or CCCCCCCCCCCC.
• Direction and base are perfectly coupled (+1→poly-G, −1→poly-C), because BAM stores SEQ reference-forward.
• The walls pile up at reproducible positions because the aligner extends the poly-G into a reference G-run and clips at its last base (REL606 8154, 16852, 39932 all end …GGCGG on the read's own strand).
• Only 16/1040 sat within 20 bp of any JC.

The existing tail-consensus test is structurally blind to this: tails that are all G agree with each other perfectly.

The discriminator

Strand. The artifact is always the read's 3' end, so for a given clip direction it can only come from one strand — direction +1 from forward reads, direction −1 from reverse reads. Reads clipped at a real breakpoint come from both.

993 of 1040 accepted calls (95.5%) had every clipped read on a single strand. Verified read-by-read in the BAMs:

position	call	fw	rv
11832 3530893 +1	poly-G	10	0
11832 4291586 −1	poly-C	0	9
11830 3949767 −1	GGTGCAGGTGCG, f=1.00	31	34
965A 2563919 +1	TACTTTATGTAC	11	22
JEB1698 1270661 −1	GGGTGCAGTACA, f=1.00	6	4

Changes

Two new Tier-2 soft rejects in add_sc_evidence():

• FISHER_STRAND (--soft-clipping-fisher-strand-p-value-cutoff, default 0.05) — Fisher's exact test of the agreeing clipped reads' strand split against the strand split of the reads that read through the same position, mirroring the RA polymorphism test and reusing fisher_exact_test_2x2. Comparing against local coverage rather than 50/50 keeps a genuinely strand-skewed pileup from reading as a strand-skewed clip. Where a position has no read-through at all — the frequency = 1.000 case, which carries the highest clip counts — the contingency row is empty and Fisher returns 1 for anything, so it falls back to the genome-wide spanning strand ratio.
• LOW_COMPLEXITY_TAIL (--soft-clipping-maximum-tail-homopolymer-fraction 0.66, --soft-clipping-maximum-tail-base-fraction 0.75) — reaches the 5–7 read positions where the strand test has no power.

Neither changes p0, rho or the score. The artifact clips stay in the null, which keeps it conservative; removing them would instead make everything else more significant.

Supporting work:

• Tabulation splits both spanning coverage and tail agreement by read strand. The strand rides in bit 63 of the existing packed tail (21 bases × 3 bits fills only bits 0–62), so it costs no extra memory. Counts file bumped to #sc_format=3.
• SC evidence table gains agree +/− and strand p columns.
• New SC gates table in summary.html. Its load-bearing row is the fraction of clip events that saw only one read strand — 98.1% and 99.9% in the long tests here. That number is what explains a library predicting implausibly many SC items, and nothing in the report showed it before.

Verification

All 95 tests pass (make test-long, 23:17 total).

Long-test effect:

test	accepted SC	removed
long_ltee_ara_m1_40k_pe36	3 → 0	all three homopolymer tails
long_ltee_ara_p1_50k_pe101	6 → 2	4 poly-G/poly-C
long_ltee_ara_p3_30k_pe150	8 → 2	5 poly-G + CTGTCTCTTATA (Nextera adapter, 26 fw / 0 rv)

Replayed over the Ara-2 set: 1040 → 40 accepted, 0–9 per sample.

The tabulation change is provably inert: across all four SC tests every pre-existing SC field (read_count, total_count, agree_read_count, log10_e_value, frequency, clipped_sequence, consensus_fraction) is byte-identical, no SC row appeared or disappeared, both strand splits sum exactly to their totals, and all 1941/1103/2915 non-SC golden lines are unchanged.

Note

p3_30k_pe150 position 2005207 survives at 8 fw / 0 rv (p = 0.077, just above 0.05) because its tail isn't low-complexity. Tightening to p < 0.10 would catch it, at some cost in genuine low-count calls.

🤖 Generated with Claude Code