29afbb7595289b958dfb4fe966dc51c288d260cb mspeir Sat Aug 15 13:16:06 2026 -0700 Count unrecognized G2P confidence values per record, not per output line, refs #38070 An HGNC ID can have several coordinate rows, so the old counter reported one record as many. Also fold case/whitespace so one new value is not logged as several, rename the remaining snake_case locals, and make makeSingleCellSignalsPeaksRa.py executable. Output is unchanged. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> diff --git src/hg/utils/otto/g2p/doG2p.py src/hg/utils/otto/g2p/doG2p.py index f23a5a82395..8e1145b23cc 100755 --- src/hg/utils/otto/g2p/doG2p.py +++ src/hg/utils/otto/g2p/doG2p.py @@ -89,133 +89,142 @@ # Confidence value -> itemRgb color. Unrecognized values fall back to DEFAULT_COLOR # (black) and are counted/logged by joinAndWrite so a source change is visible. CONFIDENCE_COLORS = { "definitive": "39,103,73", # dark green "strong": "56,161,105", # green "moderate": "104,211,145", # light green "limited": "252,129,129", # pink "disputed": "229,62,62", # red "refuted": "155,44,44", # dark red } DEFAULT_COLOR = "0,0,0" # black, for unrecognized confidence values +def normalizeConfidence(confidence): + """Fold a confidence string to its lookup form, so that case and stray + whitespace do not make one value look like several.""" + return confidence.lower().strip() + + def confidenceToColor(confidence): """Return the itemRgb color for a confidence string, or None if unrecognized.""" - return CONFIDENCE_COLORS.get(confidence.lower().strip()) + return CONFIDENCE_COLORS.get(normalizeConfidence(confidence)) -def loadG2p(file_path): +def loadG2p(filePath): """Load G2P CSV into a dict keyed by HGNC ID (each value is a list of rows).""" - g2p_map = {} + g2pMap = {} numOfRows = 0 - with open(file_path, newline="", encoding="utf-8") as csvfile: + with open(filePath, newline="", encoding="utf-8") as csvfile: reader = csv.DictReader(csvfile) for row in reader: numOfRows += 1 - hgnc_id = row["hgnc id"].strip() - g2p_map.setdefault(hgnc_id, []).append(row) + hgncId = row["hgnc id"].strip() + g2pMap.setdefault(hgncId, []).append(row) print("Number of rows in file: %s" % numOfRows) - return g2p_map + return g2pMap -def loadCoordinates(db, hgnc_ids): +def loadCoordinates(db, hgncIds): """Build a dict of gene coordinates for the given HGNC IDs from the HGNC bigBed. One bigBedToBed pass over the whole track (~49k rows) instead of one bigBedNamedItems subprocess per HGNC ID. The bigBed name field is "HGNC:<id>"; the G2P CSV stores the bare numeric id, so we key on that. """ - wanted = set(hgnc_ids) - coord_map = {} + wanted = set(hgncIds) + coordMap = {} hgncBB = "/gbdb/%s/hgnc/hgnc.bb" % db for line in bash("bigBedToBed %s stdout" % hgncBB).split("\n"): if not line.strip(): continue fields = line.split("\t")[:8] name = fields[3] # e.g. "HGNC:36036" - id = name.split("HGNC:")[-1] - if id in wanted: - coord_map.setdefault(id, []).append(fields) - return coord_map + hgncId = name.split("HGNC:")[-1] + if hgncId in wanted: + coordMap.setdefault(hgncId, []).append(fields) + return coordMap -def joinAndWrite(g2p_data, coords, output_file): - """Join G2P records and HGNC coordinates into BED 9+20 and write to output_file. +def joinAndWrite(g2pData, coords, outputFile): + """Join G2P records and HGNC coordinates into BED 9+20 and write to outputFile. - Returns a stats dict: + Returns a stats dict, both counts in G2P records so they are comparable: "unmatched" -> count of G2P records whose HGNC ID had no coordinate match in this assembly's HGNC track (they are skipped). - "unknownConfidence" -> {confidence value: count} for values not in - CONFIDENCE_COLORS (colored black). + "unknownConfidence" -> {normalized confidence value: count of records} for + values not in CONFIDENCE_COLORS (colored black). """ unmatched = 0 unknownConfidence = {} - with open(output_file, "w", newline="", encoding="utf-8") as out: + with open(outputFile, "w", newline="", encoding="utf-8") as out: writer = csv.writer(out, delimiter="\t") - for hgnc_id, rows in g2p_data.items(): - matches = coords.get(hgnc_id, []) + for hgncId, rows in g2pData.items(): + matches = coords.get(hgncId, []) if not matches: unmatched += len(rows) continue for row in rows: - for coord in matches: - # BED 9 fields - chrom = coord[0] - chromStart = coord[1] - chromEnd = coord[2] - name = row["gene symbol"] - score = coord[4] - strand = coord[5] - thickStart = coord[6] - thickEnd = coord[7] + # Counted once per G2P record, not once per output line: an HGNC ID + # can carry several coordinate rows, which would inflate the tally. rgb = confidenceToColor(row["confidence"]) if rgb is None: - unknownConfidence[row["confidence"]] = \ - unknownConfidence.get(row["confidence"], 0) + 1 + key = normalizeConfidence(row["confidence"]) + unknownConfidence[key] = unknownConfidence.get(key, 0) + 1 rgb = DEFAULT_COLOR # G2P 20 fields - g2p_id = row["g2p id"] - gene_mim = row["gene mim"] - hgnc_id_val = row["hgnc id"] - prev_symbols = row["previous gene symbols"].replace(";", ",") - disease_name = row["disease name"] - disease_mim = row["disease mim"] - disease_MONDO = row["disease MONDO"] - allelic_req = row["allelic requirement"] - cross_mod = row["cross cutting modifier"] + g2pId = row["g2p id"] + geneMim = row["gene mim"] + hgncIdVal = row["hgnc id"] + prevSymbols = row["previous gene symbols"].replace(";", ",") + diseaseName = row["disease name"] + diseaseMim = row["disease mim"] + diseaseMondo = row["disease MONDO"] + allelicReq = row["allelic requirement"] + crossMod = row["cross cutting modifier"] confidence = row["confidence"] - var_conseq = row["variant consequence"] - var_types = row["variant types"] - mol_mech = row["molecular mechanism"] - mol_mech_cat = row["molecular mechanism categorisation"] - mol_mech_ev = row["molecular mechanism evidence"] + varConseq = row["variant consequence"] + varTypes = row["variant types"] + molMech = row["molecular mechanism"] + molMechCat = row["molecular mechanism categorisation"] + molMechEv = row["molecular mechanism evidence"] phenotypes = row["phenotypes"].replace(";", ",") publications = row["publications"].replace(";", ",") panel = row["panel"] comments = row["comments"] - date_review = row["date of last review"] + dateReview = row["date of last review"] + + for coord in matches: + # BED 9 fields + chrom = coord[0] + chromStart = coord[1] + chromEnd = coord[2] + name = row["gene symbol"] + score = coord[4] + strand = coord[5] + thickStart = coord[6] + thickEnd = coord[7] writer.writerow([ chrom, chromStart, chromEnd, name, score, strand, thickStart, thickEnd, - rgb, g2p_id, gene_mim, hgnc_id_val, prev_symbols, disease_name, disease_mim, - disease_MONDO, allelic_req, cross_mod, confidence, var_conseq, var_types, - mol_mech, mol_mech_cat, mol_mech_ev, phenotypes, publications, panel, - comments, date_review, + rgb, g2pId, geneMim, hgncIdVal, prevSymbols, diseaseName, diseaseMim, + diseaseMondo, allelicReq, crossMod, confidence, varConseq, varTypes, + molMech, molMechCat, molMechEv, phenotypes, publications, panel, + comments, dateReview, ]) return {"unmatched": unmatched, "unknownConfidence": unknownConfidence} def itemCount(bb): line = bash('bigBedInfo %s | grep "itemCount"' % bb) return int(line.rstrip().split("itemCount:")[1].replace(",", "").strip()) def checkItemCount(db, newBb): """Abort if the item count moved more than COUNT_TOLERANCE vs the live track.""" liveBb = GBDB_BB % db if not Path(liveBb).exists(): print("%s: no live bigBed yet, skipping item-count check" % db) return @@ -240,44 +249,44 @@ print("Installed %s -> %s" % (liveBb, newBb)) def main(): if not updateNeeded(): # Silent no-op: nothing new from G2P this run. return validateColumns(NEW_CSV) date = str(datetime.now()).split(" ")[0] buildDir = "%s/%s" % (WORKDIR, date) bash("mkdir -p %s" % buildDir) bash("cp %s %s/AllG2P.csv" % (NEW_CSV, buildDir)) - g2p_data = loadG2p(NEW_CSV) - hgnc_ids = list(g2p_data.keys()) - print("Number of HGNC IDs found: %s" % len(hgnc_ids)) + g2pData = loadG2p(NEW_CSV) + hgncIds = list(g2pData.keys()) + print("Number of HGNC IDs found: %s" % len(hgncIds)) - coordsByDb = {db: loadCoordinates(db, hgnc_ids) for db in DBS} + coordsByDb = {db: loadCoordinates(db, hgncIds) for db in DBS} for db in DBS: print("Loaded %s %s HGNC IDs" % (len(coordsByDb[db]), db)) builtBb = {} for db in DBS: bedFile = "%s/%s_g2p_all.bed" % (buildDir, db) bbFile = "%s/%s_g2p.bb" % (buildDir, db) twoBit = "/gbdb/%s/%s.2bit" % (db, db) - stats = joinAndWrite(g2p_data, coordsByDb[db], bedFile) + stats = joinAndWrite(g2pData, coordsByDb[db], bedFile) print("Wrote %s" % bedFile) if stats["unmatched"]: print("%s: %d G2P record(s) had no HGNC coordinate match and were skipped" % (db, stats["unmatched"])) for conf, n in sorted(stats["unknownConfidence"].items()): print("%s: unrecognized confidence value %r on %d record(s); colored black" % (db, conf, n)) bash("bedToBigBed -type=bed9+20 -tab -sort " "-as=%s -sizesIs2Bit -extraIndex=name,g2p_id,gene_mim,hgnc_id %s %s %s" % (AS_FILE, bedFile, twoBit, bbFile)) print("Built %s" % bbFile) builtBb[db] = bbFile # Safety check before swapping anything live. for db in DBS: