78cdae7249c8609dcbc743e996ea7e5eec33d75a max Mon Aug 17 08:15:39 2026 -0700 lrSv: fix off-by-one anchor base in deletion coordinates across converters, refs #38099 VCF/pangenome deletions carry a non-deleted anchor (padding) base at POS. Several lrSv converters set chromStart = pos-1, which includes that anchor, so each deletion was 1 bp too wide on the left and svLen was 1 too big. Callsets handled this inconsistently, so the same deletion appeared at offset coordinates and failed to merge in lrSvAll. For deletions only (INS/INV/CPX unchanged), advance chromStart past the anchor so the interval covers exactly the deleted bases (svLen == |SVLEN|). Verified against the hg38 reference: the old left base is present in both REF and ALT (i.e. retained by the sample), so it should not be inside the deletion. Fixed 11 converters: lrSv1kLin1218VcfToBed, lrSv1kgOntVcfToBed, lrSvGustafsonVcfToBed, lrSvGa4kSvVcfToBed, lrSvDecodeVcfToBed, lrSvAou1kCsvToBed, lrSvColorsDbSvVcfToBed, lrSvCardBbToBed, lrSvAprVcfToBed, lrSvCpc1VcfToBed, lrSvVcfToBed (generic, used by han945). Left unchanged, verified already anchor-correct: hgsvc3 and hgsvc2 (0-based source), hprc2v21 (Ro converter prefix-trims), noyvert/tommoJp (POS is the first deleted base), chirmade101 (1-based-closed source). Rebuilt all affected bigBeds (hg38 + hs1 where present) and the lrSvAll merge: 3,111,026 -> 2,963,093 rows as ~148k duplicate deletions now merge. diff --git src/hg/makeDb/scripts/lrSv/lrSvVcfToBed.py src/hg/makeDb/scripts/lrSv/lrSvVcfToBed.py index 6c48315c4d2..65eb708f1af 100644 --- src/hg/makeDb/scripts/lrSv/lrSvVcfToBed.py +++ src/hg/makeDb/scripts/lrSv/lrSvVcfToBed.py @@ -1,135 +1,140 @@ #!/usr/bin/env python3 """Convert a SURVIVOR-merged SV VCF (site-only) to BED9+ for bigBed. Usage: lrSvVcfToBed.py input.vcf.gz output.bed """ import gzip import os import sys sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) from lrSvCommon import svName, normalizeSvType, insLenFor, svColor def parseInfo(infoStr): """Parse INFO field into a dict.""" d = {} for item in infoStr.split(";"): if "=" in item: k, v = item.split("=", 1) d[k] = v else: d[item] = True return d def suppVecToList(suppVec): """Convert binary support vector to comma-separated 1-based sample indices.""" indices = [] for i, c in enumerate(suppVec): if c == "1": indices.append(str(i + 1)) return ",".join(indices) if indices else "" def main(): if len(sys.argv) != 3: print(__doc__, file=sys.stderr) sys.exit(1) inFile, outFile = sys.argv[1], sys.argv[2] opener = gzip.open if inFile.endswith(".gz") else open with opener(inFile, "rt") as fIn, open(outFile, "w") as fOut: for line in fIn: if line.startswith("#"): continue fields = line.rstrip("\n").split("\t") chrom = fields[0] pos = int(fields[1]) qual = fields[5] info = parseInfo(fields[7]) svTypeRaw = info.get("SVTYPE", ".") svType = normalizeSvType(svTypeRaw) end = int(info.get("END", pos)) svLenRaw = int(float(info.get("SVLEN", "0"))) af = float(info.get("AF", "0")) supp = int(info.get("SUPP", "0")) ciPos = info.get("CIPOS", "0,0") ciEnd = info.get("CIEND", "0,0") chr2 = info.get("CHR2", ".") strands = info.get("STRANDS", "+-") suppVec = info.get("SUPP_VEC", "") # BED is 0-based half-open chromStart = pos - 1 # For INS, END == POS so the item has zero width; expand by 1 bp chromEnd = end + # For symbolic , VCF POS is the padding base before the event + # (not deleted); the deleted region is [POS+1, END]. Drop the padding + # base from the left so svLen == |SVLEN|. Only DEL moves. + if svType == "DEL": + chromStart += 1 if chromEnd <= chromStart: chromEnd = chromStart + 1 # Score: map QUAL to 0-1000 try: score = min(int(round(float(qual) * 2)), 1000) except ValueError: score = 0 # Strand from first character of STRANDS field strand = strands[0] if strands and strands[0] in "+-" else "." color = svColor(svType) # sampleList from SUPP_VEC sampleList = suppVecToList(suppVec) # end2 for TRA; empty for non-TRA so skipEmptyFields hides them end2 = str(end) if svType == "TRA" else "" chr2Out = chr2 if svType == "TRA" else "" # For TRA, chromEnd is the position on chr1 side, not chr2 if svType == "TRA": chromEnd = chromStart + 1 # svLen: length on reference svLen = chromEnd - chromStart # insLen: for INS use abs(SVLEN); else 0 (except TRA which is 0) if svType in ("INS", "MEI"): insLen = abs(svLenRaw) else: insLen = 0 # AC: SURVIVOR input doesn't have AC, use supp*2 as approximation # (SUPP is number of samples carrying; use 2*SUPP as proxy for diploid AC) ac = supp * 2 featLen = insLen if svType in ("INS", "MEI") else svLen name = svName(svType, featLen, ac) row = [ chrom, str(chromStart), str(chromEnd), name, str(score), strand, str(chromStart), # thickStart str(chromEnd), # thickEnd color, svType, str(svLen), str(insLen), str(ac), f"{af:.6f}", str(supp), ciPos, ciEnd, chr2Out, end2, sampleList, ] fOut.write("\t".join(row) + "\n") if __name__ == "__main__": main()