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/lrSv1kLin1218VcfToBed.py src/hg/makeDb/scripts/lrSv/lrSv1kLin1218VcfToBed.py index 631841cb88b..d58656e4dcb 100644 --- src/hg/makeDb/scripts/lrSv/lrSv1kLin1218VcfToBed.py +++ src/hg/makeDb/scripts/lrSv/lrSv1kLin1218VcfToBed.py @@ -1,131 +1,138 @@ #!/usr/bin/env python3 """Convert 1000 Genomes linear long-read SV VCF (1,218 samples) to BED9+ for bigBed. Usage: lrSv1kLin1218VcfToBed.py input.vcf.gz output.bed Input VCF was merged with Truvari v5.2.0 and annotated with bcftools fill-tags for population-level AC/AN/AF per 1000 Genomes super-populations (EUR, AMR, EAS, AFR, SAS). Only DEL and INS variant types are present. """ import gzip import os import sys sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) from lrSvCommon import svName, normalizeSvType SV_COLORS = { "DEL": "200,0,0", "INS": "0,0,200", } def parseInfo(infoStr): d = {} for item in infoStr.split(";"): if "=" in item: k, v = item.split("=", 1) d[k] = v else: d[item] = True return d def toInt(s, default=0): if not s or s in (".", "NA"): return default try: return int(float(s)) except ValueError: return default def toFloat(s, default=0.0): if not s or s in (".", "NA"): return default try: return float(s) except ValueError: return default 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 nWritten = 0 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]) info = parseInfo(fields[7]) svType = normalizeSvType(info.get("SVTYPE", ".")) svLenRaw = toInt(info.get("SVLEN", "0")) end = toInt(info.get("END", str(pos))) chromStart = pos - 1 chromEnd = end + # VCF POS is the anchor base (present in both REF and ALT); it is not + # part of a deletion, so the deleted region starts one base to the + # right. Drop the anchor from the left of DEL intervals so svLen == + # |SVLEN| and coordinates match anchor-excluded callsets (e.g. + # HGSVC3). INS keeps the anchor-based position. + if svType == "DEL": + chromStart += 1 if chromEnd <= chromStart: chromEnd = chromStart + 1 svLen = chromEnd - chromStart insLen = abs(svLenRaw) if svType in ("INS", "MEI") else 0 ac = toInt(info.get("AC", "0")) an = toInt(info.get("AN", "0")) af = toFloat(info.get("AF", "0")) afAfr = toFloat(info.get("AF_AFR", "0")) afAmr = toFloat(info.get("AF_AMR", "0")) afEas = toFloat(info.get("AF_EAS", "0")) afEur = toFloat(info.get("AF_EUR", "0")) afSas = toFloat(info.get("AF_SAS", "0")) ns = toInt(info.get("NS", "0")) numConsolidated = toInt(info.get("NumConsolidated", "0")) color = SV_COLORS.get(svType, "100,100,100") featLen = insLen if svType in ("INS", "MEI") else svLen name = svName(svType, featLen, ac) row = [ chrom, str(chromStart), str(chromEnd), name, "0", ".", str(chromStart), str(chromEnd), color, svType, str(svLen), str(insLen), str(ac), str(an), f"{af:.6f}", f"{afAfr:.6f}", f"{afAmr:.6f}", f"{afEas:.6f}", f"{afEur:.6f}", f"{afSas:.6f}", str(ns), str(numConsolidated), ] fOut.write("\t".join(row) + "\n") nWritten += 1 print(f"Wrote {nWritten} records to {outFile}", file=sys.stderr) if __name__ == "__main__": main()