[semenko]

Summary

Each row of the output matrix is now one sequenced DNA fragment instead of one read, and every read is decoded on the base that actually carries its call. This fixes two defects that each lost about half of the data, and makes the two aligner families behave the same way.

This changes the output. Files written by this code are not comparable row-for-row with earlier ones. They are marked with "format_version": 2 in metadata.json.

The two defects

1. Biscuit / bwameth / gem3: read 2 of every pair was discarded. A read was skipped when its alignment orientation disagreed with its bisulfite strand tag, on the premise that it was an uninformative "daughter strand". Both mates of a pair carry the same strand tag, and BAM stores both in top-strand orientation, so the skipped mate shows the same conversion pattern as the kept one.

Checked on tests/test.bam, which is real paired-end bwa-meth EM-seq data (16,686 usable pairs):

Check	Result
Which mate was kept	Always read 1. Read 2 was always dropped.
C content, top-strand fragments, kept vs dropped mate	0.8% vs 0.7%
G content, bottom-strand fragments, kept vs dropped mate	0.8% vs 0.7%
Same base in both mates where they overlap	99.21% of 1.3 M positions
Methylated CpG calls confirmed by the dropped mate	98.7%
Bases covered only by the dropped mate	38% more than the kept mates cover

2. Bismark: every bottom-strand call was recorded as -1. The XM call string was read at the C of each CpG. Bismark writes the call for a bottom-strand (XG:Z:GA) alignment on the G. Confirmed on real Bismark v0.23.0 output: all 7,230 top-strand CpG calls sit on the C, all 7,409 bottom-strand calls sit on the G.

Before and after, on real aligner output

The same paired-end FASTQs aligned by Bismark and by bwa-meth (from nf-core/test-datasets, full files):

	Before	After
Bismark BAM, rows	6,831 reads	4,277 fragments
Bismark, calls stored as -1	50.6%	4.4%
Bismark, rows holding nothing but -1	50.3%	under 3%
bwa-meth BAM, rows	5,864 (read 1 only)	7,353 fragments
bwa-meth, CpG calls	12,208	21,660
Fragments both files have in common that decode identically	not comparable	4,277 of 4,277

The last row is the strongest check. One file is decoded from Bismark's call string and XG tag, the other from the read bases and bwa-meth's YD tag. They now agree on every one of 12,818 calls.

The remaining -1 calls in this dataset are mates disagreeing in their overlap. These reads were simulated with each mate converted independently, so they disagree far more often than real mates do (0.3% in the real cfDNA data above).

How reads become rows

• Mates are merged when both pass the read filters and map to the same chromosome within 100 kb. Otherwise each is its own row. Single-end reads are one row each.
• A CpG covered by both mates is stored once: the shared call if they agree, the available call if only one mate has one, -1 if one says methylated and the other unmethylated.
• Strand, not orientation, decides which base is read: the C of the CpG for a top-strand fragment, the G for a bottom-strand fragment. A CpG counts as covered when that base lies within the alignment. If that base is deleted from the read, the CpG is -1.
• Row order: by chromosome, then by where each fragment's first CpG-covering read appears in the BAM. Serial and --threads output stay identical.
• --filter-non-converted still evaluates each read against the threshold, as the NEB tool does. A failing read now removes its whole fragment, since incomplete conversion is a property of the molecule.
• --filter-em-overconversion is evaluated on the merged fragment.
• tlen is the fragment length (absolute TLEN). The BAM's sign only says which mate is leftmost.

Output format flag

metadata.json gains format_version (now 2) and row_unit ("fragment"). Files without format_version are format 1. bam2tensor-inspect prints the format and labels rows Fragments or Reads accordingly. The README documents both versions and how to check.

Tests

• tests/data/sarscov2/: 557 read pairs aligned by Bismark and by bwa-meth, with their reference. A subset of nf-core/test-datasets, redistributed under its MIT licence (licence, provenance and the subsetting script are in that directory, about 180 KB in total). tests/test_real_aligners.py requires the two files to extract to identical matrices and to match an independent decode.
• tests/test_fragments.py: one hand-built case per rule above.
• tests/test_simulated.py: fragment-level truth for both dialects and both strands, serial and parallel, plus a check that the same simulated molecules give the same matrix in either dialect. The Bismark bottom-strand xfail from Robustness fixes: clear errors, atomic writes, no silent filter no-ops; simulated-truth tests #93 is now an ordinary passing test.
• The old tests asserting that "daughter strand" reads are skipped are replaced by tests asserting both orientations decode correctly.

344 tests pass. nox: all 12 sessions pass locally.

Also in this PR

Two analysis snippets in the README densified the matrix with .toarray(), which makes an uncovered CpG indistinguishable from an unmethylated one. On the fixture above the per-CpG example reported 1.7% methylation for a site that is 72.7% methylated. Both snippets now use the stored entries, and the over-conversion replay snippet was checked to reproduce --filter-em-overconversion exactly.

Performance

Same harness as #95. Where both versions decode every read (Bismark-style BAM), extraction takes about 12% longer and uses less memory. For Biscuit-style BAMs the old code skipped half the reads, so the new code does more work for the extra data: 56% more time for 75% more calls.

Not verified, and not done

• gem3 XB follows the same logic as YD, per its documented meaning, but I have no real gem3 output to check it against.
• Package version is unchanged. A bump publishes to PyPI. Given the output change, the next release probably deserves a major version.
• Mates whose read names differ (some old tools append /1 and /2) cannot be paired and stay as separate rows.

---

Note added after merge. The upstream nf-core reads turned out to be a non-directional library. The "bwa-meth BAM" rows in the before-and-after table count every fragment bwa-meth aligned, including about 3,700 complementary-strand pairs that bwa-meth tags with a YD value contradicting their bases. Calls for those pairs are unreliable in both the old and the new code, because bwa-meth assumes a directional library. The cross-aligner row (4,277 of 4,277) covers exactly the directional pairs, for which the tags are right, and the bundled fixture contains only those. See #97 for details and the documentation of this limit.

[github-actions]

Coverage report

Click to see where and how coverage changed

File	Statements	Missing	Coverage	Coverage (new stmts)	Lines missing
src/bam2tensor
__main__.py
functions.py					732, 799, 807, 811, 835, 894, 906
inspect.py
metadata.py
Project Total

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