9efe3bfe8af965d5e183e05615bc337035d8b283 hiram Fri Sep 18 13:45:24 2026 -0700 consider NCBI refSeq "reference" status in the priorities and fix the error prone manual maintained topPriorities list refs #38365 diff --git src/hg/hubApi/assemblyList.py src/hg/hubApi/assemblyList.py index 74fefc12705..0cc29c0e602 100755 --- src/hg/hubApi/assemblyList.py +++ src/hg/hubApi/assemblyList.py @@ -2,71 +2,83 @@ import subprocess import locale import gzip import sys import re import os import csv import requests from pathlib import Path from io import StringIO from datetime import datetime, UTC # priorities derived from the 'Monthly usage stats' report # from the qateam cron job running on the first day of the month -# special top priorities -topPriorities = { - 'hg38': 1, - 'mm39': 2, - 'hs1': 3, - 'hg19': 4, - 'GCF_028858775.2' : 5, # mPanTro3_v2.0 Chimp - 'GCF_049350105.2' : 6, # T2T_MMU8v2.0 Rhesus - 'GCF_029289425.2' : 7, # mPanPan1_v2.0 Bonobo - 'GCF_029281585.2' : 8, # mGorGor1_v2.1 Gorilla - 'GCF_028885655.2' : 9, # mPonAbe1_v2.0 Orangutan - 'GCF_037993035.2' : 10, # T2T_MFA8v1.1 Crab-eating macaque - 'GCA_049354715.1' : 11, # calJac240_pri Marmoset - 'GCF_011100555.1' : 12, # mCalJa1.2 Marmoset - 'GCF_040939455.1' : 13, # Inina_mat1.0 Mouse lemur - 'GCF_036323735.1' : 14, # rn8 rat - 'GCF_041296265.1' : 15, # TB_T2T horse - 'GCF_016772045.1' : 16, # ARS_UI_Ramb_v2.0 Sheep - 'GCF_002263795.3' : 17, # ARS_UCD2.0 Cow - 'GCF_018350175.1' : 18, # Fca126_mat1.0 Cat - 'GCF_016699485.2' : 19, # GRCg7b Chicken - 'GCF_003957565.2' : 20, # bTaeGut1.4 Zebra finch - 'GCF_049306965.1' : 21, # GRCz12tu Zebrafish - 'GCA_052040795.1' : 22, # GRCz12ab Zebrafish - 'mm10': 23, - 'dm6': 24, - 'danRer11': 25, - 'mm9': 26, - 'hetGla2': 27, - 'rn6': 28, - 'hg18': 29, - 'galGal6': 30, - 'bosTau9': 31, - 'ce11': 32, - 'canFam4': 33, -} +# special top priorities -- ORDER MATTERS: this is a plain list, not +# a dict, precisely so a new identifier can just be inserted wherever +# it belongs without renumbering anything by hand. initTopPriorities() +# turns list position into the actual priority number (1, 2, 3, ...). +topPriorityNames = [ + 'hg38', + 'mm39', + 'hs1', + 'hg19', + 'GCA_018852605.3', # human (NA24385 HG002 pat 2024) + 'GCA_018852615.3', # human (NA24385 HG002 mat 2024) + 'GCA_054883195.1', # human (H9 T2T hap1 2026) + 'GCA_054883165.1', # human (H9 T2T hap2 2026) + 'GCF_028858775.2', # mPanTro3_v2.0 Chimp + 'GCF_049350105.2', # T2T_MMU8v2.0 Rhesus + 'GCF_029289425.2', # mPanPan1_v2.0 Bonobo + 'GCF_029281585.2', # mGorGor1_v2.1 Gorilla + 'GCF_028885655.2', # mPonAbe1_v2.0 Orangutan + 'GCF_037993035.2', # T2T_MFA8v1.1 Crab-eating macaque + 'GCF_049354715.1', # calJac240_pri Marmoset + 'GCF_011100555.1', # mCalJa1.2 Marmoset + 'GCF_040939455.1', # Inina_mat1.0 Mouse lemur + 'GCF_036323735.1', # rn8 rat + 'GCF_041296265.1', # TB_T2T horse + 'GCF_016772045.1', # ARS_UI_Ramb_v2.0 Sheep + 'GCF_002263795.3', # ARS_UCD2.0 Cow + 'GCF_018350175.1', # Fca126_mat1.0 Cat + 'GCF_016699485.2', # GRCg7b Chicken + 'GCF_003957565.2', # bTaeGut1.4 Zebra finch + 'GCF_049306965.1', # GRCz12tu Zebrafish + 'GCA_052040795.1', # GRCz12ab Zebrafish + 'mm10', + 'dm6', + 'danRer11', + 'mm9', + 'hetGla2', + 'rn6', + 'hg18', + 'galGal6', + 'bosTau9', + 'ce11', + 'canFam4', +] + +### key will be dbDb/GCx name, value will be priority number. +### Populated by initTopPriorities() from topPriorityNames above. +topPriorities = {} ### key will be dbDb/GCx name, value will be priority number allPriorities = {} -priorityCounter = len(topPriorities) + 1 +### set for real by initTopPriorities(), once topPriorityNames is known +priorityCounter = 1 # key is clade, value is priority, will be initialized # by initCladePriority() function cladePrio = {} #################################################################### ### this is kinda like an environment setting, it gets everything ### into a UTF-8 reading mode #################################################################### def set_utf8_encoding(): """ Set UTF-8 encoding for stdin, stdout, and stderr in Python. """ if sys.stdout.encoding != 'utf-8': sys.stdout = open(sys.stdout.fileno(), mode='w', encoding='utf-8', buffering=1) @@ -222,30 +234,45 @@ hasn't changed. 2175 archaea 360585 bacteria 607 fungi 414 invertebrate 184 plant 96 protozoa 231 vertebrate_mammalian 405 vertebrate_other 14992 viral """ +#################################################################### +### turn topPriorityNames (list order = priority order) into the +### topPriorities dict, and set priorityCounter to continue right +### after it. This is the only place priority numbers 1..N get +### assigned for topPriorityNames -- add/reorder identifiers there, +### nothing here needs to change. +#################################################################### +def initTopPriorities(): + global topPriorities + global priorityCounter + + for i, name in enumerate(topPriorityNames, start=1): + topPriorities[name] = i + priorityCounter = len(topPriorityNames) + 1 + #################################################################### ### the various listings are ordered by these clade priorities to get ### primates first, mammals second, and so on #################################################################### def initCladePriority(): global cladePrio keys = [ "primates", "mammals", "vertebrate_mammalian", "birds", "fish", "vertebrate", "vertebrate_other", @@ -600,54 +627,117 @@ pat = r'\b' + re.escape(year) + r'\b' if not re.search(pat, item['commonName']): if not re.search(pat, item['taxId']): item['taxId'] += " " + year if gcAccession in status: stat = status[gcAccession] item['refSeqCategory'] = stat['refSeqCategory'].lower() item['versionStatus'] = stat['versionStatus'].lower() item['assemblyLevel'] = stat['assemblyLevel'].lower() ## pat = r'\b' + re.escape(stat) + r'\b' ## if not re.search(pat, item['taxId']): ## item['taxId'] += " " + stat return +#################################################################### +### GenArk accessions NCBI flags as the species' reference genome, +### read directly from the genark database's assemblySummary tables +### rather than re-parsing the assembly_summary_*.txt flat files a +### second time. Needed early, before establishPriorities(), so +### GenArk 'reference' assemblies can be ranked ahead of their +### non-reference siblings within each clade grouping. +#################################################################### +def readGenArkReferenceSet(): + referenceSet = set() + tables = [ + "assemblySummaryGenbank", + "assemblySummaryGenbankHistorical", + "assemblySummaryRefseq", + "assemblySummaryRefseqHistorical", + ] + for table in tables: + result = subprocess.run( + ["hgsql", "-N", "-e", + f"SELECT assemblyAccession FROM {table} WHERE refseqCategory='reference genome';", + "genark"], + stdout=subprocess.PIPE, stderr=subprocess.PIPE + ) + if result.returncode != 0: + print(f"Error executing MySQL command on {table}: {result.stderr.decode('utf-8')}") + exit(1) + for line in result.stdout.decode('utf-8').strip().split('\n'): + if line: + referenceSet.add(line) + print(f"# genArk reference-genome accessions: {len(referenceSet)}") + return referenceSet + +#################################################################### +### give key the next priority value if it doesn't have one yet. +### Returns 1 if assigned, 0 if key already had a priority. +#################################################################### +def assignPriority(key): + global allPriorities, priorityCounter + if key in allPriorities: + return 0 + allPriorities[key] = priorityCounter + priorityCounter += 1 + return 1 + +#################################################################### +### assign priorities to a list of GenArk accessions, giving every +### 'reference' genome in the list a priority ahead of the rest of +### the list, while preserving the incoming relative order within +### each of those two groups. Returns the number of accessions that +### were actually assigned a new priority. +#################################################################### +def assignGenArkBucket(gcAccessions, genArkRefCategory): + referenceFirst = [g for g in gcAccessions if genArkRefCategory.get(g, "") == "reference"] + theRest = [g for g in gcAccessions if genArkRefCategory.get(g, "") != "reference"] + itemCount = 0 + for gcAcc in referenceFirst + theRest: + itemCount += assignPriority(gcAcc) + return itemCount + #################################################################### ### for the genArk set, establish some ad-hoc priorities #################################################################### def establishPriorities(dbDb, genArk): global topPriorities global allPriorities global priorityCounter totalItemCount = 0 expectedTotal = len(dbDb) + len(genArk) print(f"### expected total: {expectedTotal:4} = {len(dbDb):4} dbDb genomes + {len(genArk):4} genArk genomes") # first priority are the specific manually selected top items itemCount = 0 for name, priority in topPriorities.items(): allPriorities[name] = priority itemCount += 1 totalItemCount += itemCount print(f"{totalItemCount:4} - total\ttopPriorities count: {itemCount:4}") primateList = extractClade('primates', genArk) mammalList = extractClade('mammals', genArk) + # lookup used by assignGenArkBucket() to rank 'reference' genArk + # assemblies ahead of their non-reference siblings within a bucket + genArkRefCategory = {item['gcAccession']: item.get('refSeqCategory', '') for item in genArk} + versionScan = {} # key is dbDb name without number version extension, # value is highest version number seen for this bare # name highestVersion = {} # key is dbDb name without number version extension, # value is the full dbDb name for the highest version # of this dbDb name allDbDbNames = {} # key is the full dbDb name, value is its version itemCount = 0 # scanning the dbDb entries, figure out the highest version number # of each name for item in dbDb: dbDbName = item['name'] splitMatch = re.match(r"([a-zA-Z]+)(\d+)", dbDbName) if splitMatch: @@ -674,166 +764,134 @@ for key in sortByValue: highVersion = highestVersion[key[0]] if highVersion not in allPriorities: # find the element in the dbDb list that matches this highVersion name highDict = next((d for d in dbDb if d.get('name') == highVersion), None) if highDict['clade'] == "primates": if highDict['scientificName'].lower() != "homo sapiens": allPriorities[highVersion] = priorityCounter priorityCounter += 1 itemCount += 1 totalItemCount += itemCount print(f"{totalItemCount:4} - total\tdbDb highest version primates count: {itemCount:4}") itemCount = 0 # and now the GenArk GCF/RefSeq homo sapiens should be lined up here next - for item in genArk: - gcAccession = item['gcAccession'] - if not gcAccession.startswith("GCF_"): - continue - if gcAccession not in allPriorities: - sciName = item['scientificName'] - if sciName.lower() == "homo sapiens": - allPriorities[gcAccession] = priorityCounter - priorityCounter += 1 - itemCount += 1 + candidates = [item['gcAccession'] for item in genArk + if item['gcAccession'].startswith("GCF_") + and item['scientificName'].lower() == "homo sapiens"] + itemCount = assignGenArkBucket(candidates, genArkRefCategory) totalItemCount += itemCount print(f"{totalItemCount:4} - total\tgenArk GCF homo sapiens count: {itemCount:4}") itemCount = 0 # and now the GenArk GCA/GenBank homo sapiens should be lined up here next # GCA/GenBank second - for item in genArk: - gcAccession = item['gcAccession'] - if not gcAccession.startswith("GCA_"): - continue - if gcAccession not in allPriorities: - sciName = item['scientificName'] - if sciName.lower() == "homo sapiens": - allPriorities[gcAccession] = priorityCounter - priorityCounter += 1 - itemCount += 1 + candidates = [item['gcAccession'] for item in genArk + if item['gcAccession'].startswith("GCA_") + and item['scientificName'].lower() == "homo sapiens"] + itemCount = assignGenArkBucket(candidates, genArkRefCategory) totalItemCount += itemCount print(f"{totalItemCount:4} - total\tgenArk GCA homo sapiens count: {itemCount:4}") itemCount = 0 # the primates, GCF/RefSeq first - for asmId, commonName in primateList.items(): + candidates = [] + for asmId in primateList: gcAcc = asmId.split('_')[0] + "_" + asmId.split('_')[1] - if not gcAcc.startswith("GCF_"): - continue - if gcAcc not in allPriorities: - allPriorities[gcAcc] = priorityCounter - priorityCounter += 1 - itemCount += 1 + if gcAcc.startswith("GCF_"): + candidates.append(gcAcc) + itemCount = assignGenArkBucket(candidates, genArkRefCategory) totalItemCount += itemCount print(f"{totalItemCount:4} - total\tgenArk GCF primates count: {itemCount:4}") itemCount = 0 # and the GCA/GenBank primates - for asmId, commonName in primateList.items(): + candidates = [] + for asmId in primateList: gcAcc = asmId.split('_')[0] + "_" + asmId.split('_')[1] - if not gcAcc.startswith("GCA_"): - continue - if gcAcc not in allPriorities: - allPriorities[gcAcc] = priorityCounter - priorityCounter += 1 - itemCount += 1 + if gcAcc.startswith("GCA_"): + candidates.append(gcAcc) + itemCount = assignGenArkBucket(candidates, genArkRefCategory) totalItemCount += itemCount print(f"{totalItemCount:4} - total\tgenArk GCA primates count: {itemCount:4}") # next are the highest versioned database mammals itemCount = 0 sortByValue = sorted(versionScan.items(), key=lambda x: x[1], reverse=True) for key in sortByValue: highVersion = highestVersion[key[0]] if highVersion not in allPriorities: # find the element in the dbDb list that matches this highVersion name highDict = next((d for d in dbDb if d.get('name') == highVersion), None) if highDict['clade'] == "mammals": allPriorities[highVersion] = priorityCounter priorityCounter += 1 itemCount += 1 totalItemCount += itemCount print(f"{totalItemCount:4} - total\tdbDb highest version mammals count: {itemCount:4}") itemCount = 0 # the mammals, GCF/RefSeq first - for asmId, commonName in mammalList.items(): + candidates = [] + for asmId in mammalList: gcAcc = asmId.split('_')[0] + "_" + asmId.split('_')[1] - if not gcAcc.startswith("GCF_"): - continue - if gcAcc not in allPriorities: - allPriorities[gcAcc] = priorityCounter - priorityCounter += 1 - itemCount += 1 + if gcAcc.startswith("GCF_"): + candidates.append(gcAcc) + itemCount = assignGenArkBucket(candidates, genArkRefCategory) totalItemCount += itemCount print(f"{totalItemCount:4} - total\tgenArk GCF mammals count: {itemCount:4}") itemCount = 0 # and the GCA/GenBank mammals - for asmId, commonName in mammalList.items(): + candidates = [] + for asmId in mammalList: gcAcc = asmId.split('_')[0] + "_" + asmId.split('_')[1] - if not gcAcc.startswith("GCA_"): - continue - if gcAcc not in allPriorities: - allPriorities[gcAcc] = priorityCounter - priorityCounter += 1 - itemCount += 1 + if gcAcc.startswith("GCA_"): + candidates.append(gcAcc) + itemCount = assignGenArkBucket(candidates, genArkRefCategory) totalItemCount += itemCount print(f"{totalItemCount:4} - total\tgenArk GCA mammals count: {itemCount:4}") itemCount = 0 # dbDb is in cladeOrder, process in that order, find highest versions for item in dbDb: dbDbName = item['name'] if dbDbName not in allPriorities: splitMatch = re.match(r"([a-zA-Z]+)(\d+)", dbDbName) if splitMatch: noVersion = splitMatch.group(1) version = allDbDbNames[dbDbName] highVersion = versionScan[noVersion] if highVersion == version: allPriorities[dbDbName] = priorityCounter priorityCounter += 1 itemCount += 1 totalItemCount += itemCount print(f"{totalItemCount:4} - total\tdbDb highest versions count: {itemCount:4}") itemCount = 0 # GCF RefSeq from GenArk next priority - for item in genArk: - gcAccession = item['gcAccession'] - if not gcAccession.startswith("GCF_"): - continue - if gcAccession not in allPriorities: - allPriorities[gcAccession] = priorityCounter - priorityCounter += 1 - itemCount += 1 + candidates = [item['gcAccession'] for item in genArk if item['gcAccession'].startswith("GCF_")] + itemCount = assignGenArkBucket(candidates, genArkRefCategory) totalItemCount += itemCount print(f"{totalItemCount:4} - total\tgenArk GCF count: {itemCount:4}") itemCount = 0 # GCA GenBank from GenArk next priority - for item in genArk: - gcAccession = item['gcAccession'] - if not gcAccession.startswith("GCA_"): - continue - if gcAccession not in allPriorities: - allPriorities[gcAccession] = priorityCounter - priorityCounter += 1 - itemCount += 1 + candidates = [item['gcAccession'] for item in genArk if item['gcAccession'].startswith("GCA_")] + itemCount = assignGenArkBucket(candidates, genArkRefCategory) totalItemCount += itemCount print(f"{totalItemCount:4} - total\tgenArk GCA count: {itemCount:4}") itemCount = 0 for entry in dbDb: dbName = entry['name'] if dbName not in allPriorities: allPriorities[dbName] = priorityCounter priorityCounter += 1 itemCount += 1 totalItemCount += itemCount print(f"{totalItemCount:4} - total\tthe rest of dbDb count: {itemCount:4}") #################################################################### @@ -846,46 +904,56 @@ print(" dbDb.hgcentral and UCSC_GI.assemblyHubList.txt file.\n") print("Usage: assemblyList.py dbDb.clade.year.acc.tsv\n") print("the dbDb.clade.year.tsv file is a manually curated file to relate") print(" UCSC database names to GenArk clades, in source tree hubApi/") print("This script is going to read the dbDb.hgcentral table, and the file") print(" UCSC_GI.assemblyHubList.txt from hgdownload.") print("Writing an output file assemblyList.tsv to be loaded into") print(" assemblyList.hgcentral. The output file needs to be sorted") print(" sort -k2,2n assemblyList.tsv before table load.") sys.exit(255) # Ensure stdout and stderr use UTF-8 encoding set_utf8_encoding() initCladePriority() + initTopPriorities() dbDbNameCladeFile = sys.argv[1] csvHaplotypes = readHaplotypes("/hive/data/outside/ncbi/genomes/reports/haploscan/csvHaplotypes.tsv.gz") # the correspondence of dbDb names to GenArk clade categories dbDbClades, dbDbYears, dbDbNcbi = dbDbCladeList(dbDbNameCladeFile) # Get the dbDb.hgcentral table data rawData = dbDbData() dbDbItems = processDbDbData(rawData, dbDbClades, dbDbYears, dbDbNcbi) dbDbItems = dropGenarkCuratedDuplicates(dbDbItems) aliasData = asmAliasData() # read the GenArk data from hgdownload into a list of dictionaries genArkUrl = "https://hgdownload.soe.ucsc.edu/hubs/UCSC_GI.assemblyHubList.txt" genArkItems = readGenArkData(genArkUrl) + # stamp refSeqCategory early so establishPriorities() can rank + # 'reference' genArk assemblies ahead of their non-reference + # siblings within each clade grouping; readAsmSummary()/ + # addYearsStatus() below re-stamp this later from the flat files + # for the final output value, same as before this feature existed + genArkReferenceSet = readGenArkReferenceSet() + for item in genArkItems: + item['refSeqCategory'] = "reference" if item['gcAccession'] in genArkReferenceSet else "" + establishPriorities(dbDbItems, genArkItems) asmIdClade = readAsmIdClade() commonNames = allCommonNames() print("# all common names: ", len(commonNames)) refSeqList, refSeqYears, refSeqStatus = readAsmSummary("refseq.txt", allPriorities, commonNames, asmIdClade) print("# refSeq assemblies: ", len(refSeqList)) refSeqListHist, refSeqYearsHist, refSeqStatusHist = readAsmSummary("refseq_historical.txt", allPriorities, commonNames, asmIdClade) print("# refSeq historical assemblies: ", len(refSeqListHist)) genBankList, genBankYears, genBankStatus = readAsmSummary("genbank.txt", allPriorities, commonNames, asmIdClade) print("# genBank assemblies: ", len(genBankList)) genBankListHist, genBankYearsHist, genBankStatusHist = readAsmSummary("genbank_historical.txt", allPriorities, commonNames, asmIdClade) print("# genBank historical assemblies: ", len(genBankListHist)) ### dictionary unpacking, combine both dictionaries (Python 3.5+)