e461209cf1fd3758d63641915cd91ca9c8ab3020 lrnassar Thu Sep 24 17:00:10 2026 -0700 Remove tool output accidentally left in the mei description page, per CR. refs #37524 getTrackReferences writes its diagnostics to stdout rather than stderr, so six "Failed to fetch complete links from NCBI" lines ended up in the References section of mei.html and rendered as visible text on the track description page. NCBI is still not answering, so rather than rerun the tool the references are now assembled from the citation blocks already present on the six subtrack pages. That also restores the publisher links for every paper, which the failed lookups had degraded to bare PubMed URLs. Also make the INFO SEQ guard in meiHgsvc3CsvToBed.py require a usable string, so an empty SEQ= value would fall back to the ALT-derived sequence instead of silently producing an empty one. No record in either callset carries an empty SEQ today and the rebuilt output is byte-identical. diff --git src/hg/makeDb/scripts/mei/meiHgsvc3CsvToBed.py src/hg/makeDb/scripts/mei/meiHgsvc3CsvToBed.py index 310200c8502..9de1c944c8b 100755 --- src/hg/makeDb/scripts/mei/meiHgsvc3CsvToBed.py +++ src/hg/makeDb/scripts/mei/meiHgsvc3CsvToBed.py @@ -1,205 +1,206 @@ #!/usr/bin/env python3 """Convert HGSVC3 MEI_Callset_GRCh38 CSV to a bed9+ file. The HGSVC3 MEI_Callset_GRCh38 callset (Logsdon et al. 2025, Nature) lists novel Mobile Element Insertions identified in 65 long-read assembled samples relative to the GRCh38 reference. Each row of the input CSV is VCF-like (CHROM, POS, REF, ALT, QUAL, FILTER, INFO, FORMAT, <65 sample columns>, Caller_Count, TE_Designation, L1ME-AID, PALMER, L1ME-AID_INFO, PALMER_INFO, PAVMergedCalls). For BED, insertions are drawn as a 1bp anchor block at the position where the insertion would attach in the reference: chromStart = POS - 1 (1-based VCF anchor -> 0-based) chromEnd = chromStart + 1 Genotypes are encoded "h1|h2|h3" where each haplotype is 0 (reference, no insertion), 1 (insertion present), or "." (missing). """ import csv import gzip import sys COLOR_BY_CLASS = { "Alu": "0,114,178", # blue (Okabe-Ito) "L1": "213,94,0", # vermillion (Okabe-Ito) "SVA": "0,158,115", # bluish green (Okabe-Ito) "HERVK": "204,121,167", # reddish purple (Okabe-Ito) "snRNA": "0,0,0", # black (Okabe-Ito) "Other": "0,0,0", } def shortClass(teDesignation): """Map 'SINE/Alu' -> 'Alu', 'LINE/L1' -> 'L1', 'Retroposon/SVA' -> 'SVA', 'HERVK' -> 'HERVK', 'snRNA' -> 'snRNA'.""" if "/" in teDesignation: return teDesignation.split("/", 1)[1] return teDesignation def parseInfo(infoField): fields = {} for tok in infoField.split(";"): if "=" in tok: k, v = tok.split("=", 1) fields[k] = v else: fields[tok] = True return fields def main(): if len(sys.argv) != 4: sys.exit("usage: meiHgsvc3CsvToBed.py input.csv.gz chrom.sizes output.bed") inFn, sizesFn, outFn = sys.argv[1:] chromSizes = {} with open(sizesFn) as f: for line in f: chrom, size = line.rstrip("\n").split("\t") chromSizes[chrom] = int(size) nIn = 0 nOut = 0 nSkippedChrom = 0 nSkippedBoundary = 0 opener = gzip.open if inFn.endswith(".gz") else open out = open(outFn, "w") with opener(inFn, "rt") as f: reader = csv.reader(f) header = next(reader) # Locate column indices by name (more robust than positional). idx = {name: i for i, name in enumerate(header)} # Sample columns are everything between FORMAT and Caller_Count. firstSampleCol = idx["FORMAT"] + 1 lastSampleCol = idx["Caller_Count"] sampleNames = header[firstSampleCol:lastSampleCol] for row in reader: nIn += 1 chrom = row[idx["CHROM"]] pos = int(row[idx["POS"]]) # 1-based VCF position of anchor ref = row[idx["REF"]] alt = row[idx["ALT"]] info = parseInfo(row[idx["INFO"]]) sampleGts = row[firstSampleCol:lastSampleCol] teDesignation = row[idx["TE_Designation"]] # Insertion length: prefer SVLEN from INFO, else len(ALT)-len(REF). if "SVLEN" in info: svLen = abs(int(info["SVLEN"])) else: svLen = len(alt) - len(ref) if svLen <= 0: sys.stderr.write(f"WARN: non-positive svLen at {chrom}:{pos}\n") continue # Anchor-base BED interval (1bp wide). chromStart = pos - 1 chromEnd = chromStart + 1 if chrom not in chromSizes: nSkippedChrom += 1 sys.stderr.write(f"WARN: skipping record at {chrom}:{pos}, chrom not in chrom.sizes\n") continue if chromEnd > chromSizes[chrom]: nSkippedBoundary += 1 sys.stderr.write( f"WARN: skipping record at {chrom}:{pos}, end {chromEnd} > chromSize {chromSizes[chrom]}\n") continue # Per-record genotype tallying. altAC = 0 AN = 0 carrierCount = 0 sampleCount = 0 carrierSamples = [] for sampleName, gt in zip(sampleNames, sampleGts): hapCalls = gt.split("|") hapAlt = 0 hapCalled = 0 for h in hapCalls: if h == ".": continue AN += 1 hapCalled += 1 if h == "1": altAC += 1 hapAlt += 1 elif h != "0": sys.stderr.write(f"WARN: unexpected GT '{h}' at {chrom}:{pos} sample {sampleName}\n") if hapCalled > 0: sampleCount += 1 if hapAlt > 0: carrierCount += 1 carrierSamples.append(sampleName) altAF = (altAC / AN) if AN > 0 else 0.0 score = max(0, min(1000, int(round(altAF * 1000)))) cls = shortClass(teDesignation) color = COLOR_BY_CLASS.get(cls, COLOR_BY_CLASS["Other"]) refSd = float(info.get("REF_SD", 0)) if info.get("REF_SD", False) else 0.0 refTrf = "True" if info.get("REF_TRF", False) else "False" sourceSample = info.get("SAMPLE", "") callerCount = int(row[idx["Caller_Count"]]) if row[idx["Caller_Count"]] else 0 l1meAid = "Yes" if row[idx["L1ME-AID"]] == "1" else "No" palmer = "Yes" if row[idx["PALMER"]] == "1" else "No" # Inserted DNA. Most records follow the VCF convention of carrying the # anchor base at ALT[0] (= REF[0]), so the element is ALT minus that # base. The PALMER-only records do not: ALT there is the element # itself (len(ALT) == SVLEN, and ALT[0] usually differs from REF[0]), # and some of them carry a truncated ALT. Those records supply the # full element in INFO SEQ, which always matches SVLEN, so prefer it. - if "SEQ" in info and info["SEQ"] is not True: - insertSeq = info["SEQ"] + seq = info.get("SEQ") + if isinstance(seq, str) and seq: + insertSeq = seq else: insertSeq = alt[1:] if len(alt) > 1 else "" # Name format: <class>-<svLen>:<carrierCount> (e.g. Alu-281:33). name = f"{cls}-{svLen}:{carrierCount}" out.write("\t".join([ chrom, str(chromStart), str(chromEnd), name, str(score), ".", str(chromStart), str(chromEnd), color, cls, teDesignation, str(svLen), str(altAC), str(AN), f"{altAF:.4f}", str(carrierCount), str(sampleCount), sourceSample, str(callerCount), l1meAid, palmer, f"{refSd:.2f}", refTrf, ",".join(carrierSamples), insertSeq, ]) + "\n") nOut += 1 out.close() sys.stderr.write( f"Read {nIn} records, wrote {nOut}, skipped {nSkippedChrom} (chrom missing) + " f"{nSkippedBoundary} (off chrom end)\n") if __name__ == "__main__": main()