[semenko]

Summary

Adds a per-BAM check that catches wrongly strand-tagged reads, the situation found in #97: a non-directional library aligned with bwa-meth, where most complementary-strand pairs get a YD tag that contradicts their bases and every CpG in them is then called methylated. Until now that was silent.

How it works

A read reports on one strand, so only one base of each CpG can show a conversion. The other base never does: the G for a top-strand read, the C for a bottom-strand read. Wherever a read covers that other base, the extractor now looks at it and counts the CpGs where it shows the other strand's conversion product (A at the G, T at the C).

• Warning. If more than 2% of at least 1,000 checked CpGs contradict their read's tag, a RuntimeWarning is raised. The CLI prints it in that BAM's report. The output file is still written.
• Verbose. With --verbose the counts are printed for every BAM, warning or not.
• Scope. Applies to reads decoded from their bases: Biscuit, bwameth, gem3. Bismark reads are not checked, because Bismark resolves the strand itself.

Does it separate the cases?

Data	CpGs checked	Contradictions
Real bwa-meth output, directional pairs (bwameth.bam)	1,617	0.00%
Real bwa-meth output, complementary-strand pairs (bwameth_nondirectional.bam, new)	1,542	15.9%
Upstream nf-core file as a whole	24,325	6.6%
Real directional EM-seq cfDNA (tests/test.bam), sequencing errors only	18,406	0.05%

The 2% limit sits well above the error rate and well below the mis-tagged case.

Limits

• It warns, it does not repair. Recovering the right calls would need the reference sequence in memory, to infer each read's strand from its non-CpG bases.
• The SNP contribution is not measured. CpG SNPs also produce contradictions, and the real test BAM has no matching reference to measure them against. I expect well under 1%. Running once with --verbose on a real BAM with its reference will show the actual figure.
• The signal scales with unmethylated CpGs. A very heavily methylated sample would sit closer to the limit.
• Cost. About 6% more time on the base-decoding path, measured on 240k simulated reads. The output is unchanged.

Tests

378 pass. nox: all 12 sessions pass locally.

• Hand-built reads: exactly what is and is not counted, for both strands and both mates, and each threshold condition.
• Simulated reads with 40% of fragments deliberately mis-tagged: the warning's counts are checked against an independent count, in serial and parallel runs. Correctly tagged simulated reads, which include a realistic trickle of contradicting sequencing errors, raise nothing.
• Real aligner output: the directional fixtures pass, and the new fixture of 634 real mis-tagged pairs is flagged. That fixture is another subset of the same MIT-licensed nf-core data, about 86 KB, produced by the same documented script.

[github-actions]

Coverage report

Click to see where and how coverage changed

File	Statements	Missing	Coverage	Coverage (new stmts)	Lines missing
src/bam2tensor
functions.py					988
Project Total

This report was generated by python-coverage-comment-action