Skip to content

Commit 1a6318a

Browse files
author
Jeremiah Wala
committed
fixed VCF 0-based to 1-based
1 parent 7fd1208 commit 1a6318a

2 files changed

Lines changed: 21 additions & 5 deletions

File tree

CLAUDE.md

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -459,6 +459,17 @@ Gotcha: the BAM is required only for the chromosome name/length table
459459
BreakPoint parser to turn chrom-name strings into chr IDs). No reads are
460460
actually read from it — any BAM that shares the reference is fine.
461461

462+
POS convention (fixed): breakend positions are stored 0-based internally
463+
(htslib), 1-based in bps.txt.gz (`BreakPoint::toFileString` adds +1) and
464+
1-based in the VCF. `VCFEntry::toFileString` / `getAltString` (`vcf.cpp`)
465+
previously emitted the raw 0-based `gr.pos1` for the POS column and the
466+
BND mate locus, so tovcf VCFs were off by one (one *less* than bps.txt) —
467+
even though the symbolic `END` INFO field already added +1, so symbolic
468+
records had a POS/END mismatch too. All three now add +1 and agree. The
469+
*internal* `gr.pos1` uses in `vcf.cpp` (dedup interval tree, id-hash
470+
strings, the entry sort comparator) stay 0-based on purpose — only the
471+
emitted POS / ALT-mate / END are 1-based.
472+
462473
Not-yet-done on this subcommand: bgzip-proper (current `--plain=false`
463474
output is plain gzip, which bcftools accepts but tabix doesn't index
464475
correctly). For now, pipe through `bcftools sort -Oz` + `tabix -p vcf`

src/svaba/vcf.cpp

Lines changed: 10 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -230,11 +230,13 @@ std::string VCFEntry::getAltString(const SeqLib::BamHeader& header) const {
230230
// SV breakends get the BND-style `N]chr:pos]` mate notation.
231231
const std::string ref = getRefString();
232232

233+
// +1: positions are 0-based internally (htslib); the mate locus in BND
234+
// notation is 1-based per the VCF spec, matching POS in toFileString.
233235
std::stringstream ptag;
234236
if (id_num == 1) {
235-
ptag << bp->b2.gr.ChrName(header) << ':' << bp->b2.gr.pos1;
237+
ptag << bp->b2.gr.ChrName(header) << ':' << (bp->b2.gr.pos1 + 1);
236238
} else {
237-
ptag << bp->b1.gr.ChrName(header) << ':' << bp->b1.gr.pos1;
239+
ptag << bp->b1.gr.ChrName(header) << ':' << (bp->b1.gr.pos1 + 1);
238240
}
239241

240242
std::stringstream alt;
@@ -435,17 +437,20 @@ std::string VCFEntry::toFileString(const SeqLib::BamHeader& header,
435437
// positions (the event's 5' anchor); the INFO/END field carries the
436438
// other boundary. For BND records, POS is just the breakend this
437439
// entry represents (legacy behavior).
440+
// +1: breakend positions are stored 0-based internally (htslib
441+
// convention); VCF POS is 1-based. This mirrors BreakPoint::toFileString
442+
// (bps.txt.gz) and the END INFO field below, both of which add +1.
438443
int pos;
439444
std::string chr_name;
440445
if (symbolic_rep) {
441-
const int p1 = bp->b1.gr.pos1;
442-
const int p2 = bp->b2.gr.pos1;
446+
const int p1 = bp->b1.gr.pos1 + 1;
447+
const int p2 = bp->b2.gr.pos1 + 1;
443448
pos = std::min(p1, p2);
444449
// Both breakends live on the same chrom when symbolic; use b1's.
445450
chr_name = bp->b1.gr.ChrName(header);
446451
} else {
447452
const BreakEnd* be = (id_num == 1) ? &bp->b1 : &bp->b2;
448-
pos = be->gr.pos1;
453+
pos = be->gr.pos1 + 1;
449454
chr_name = be->gr.ChrName(header);
450455
}
451456

0 commit comments

Comments
 (0)