62f11347129c7aa7067c03879d2b6f3d9b05426a max Thu Sep 17 10:53:15 2026 -0700 uniprot otto: a skip list of assemblies to leave alone #Preview2 week - bugs introduced now will need a build patch to fix --skipList names a file of assemblies not to build, one per line, with # comments so each entry can say why. Entries may carry the full asmId with its assembly-name suffix, which is how the GenArk orderLists write them, and only the accession part is matched, so an orderList can be handed over unchanged: --skipList=kent/src/hg/makeDb/doc/hprcAsmHub/hprc.orderList.tsv That is the case it was written for. GenArk lists 572 human assemblies, 464 of them HPRC haplotypes, and the pool works taxon by taxon, so all 572 would have run one after another in a single slot while the other 663 taxa finished in hours. Dropping the HPRC ones takes human to 108 assemblies and the whole 20000-protein plan from 2224 to 1760. Applies to the dbDb plan as well, alongside notAutoDbs, so a troublesome assembly can be kept out of the monthly run without editing the script. refs #38300 diff --git src/hg/utils/otto/uniprot/doUniprot src/hg/utils/otto/uniprot/doUniprot index 217673d3061..a31f64630bf 100755 --- src/hg/utils/otto/uniprot/doUniprot +++ src/hg/utils/otto/uniprot/doUniprot @@ -157,91 +157,124 @@ Tab separated: accession, asmName, scientificName, commonName, taxId, clade. """ accToTax = {} # latin1: the organism and common names carry bytes that are not valid UTF-8, and # python would stop on the first one. Only the accession and the taxon id are used. for line in open(fname, encoding="latin1"): if line.startswith("#"): continue f = line.rstrip("\n").split("\t") if len(f) < 5 or not f[4].isdigit(): continue accToTax[f[0]] = int(f[4]) logging.info("Read %d GenArk assemblies from %s" % (len(accToTax), fname)) return accToTax +def readSkipList(fname): + """ read a file of assemblies to leave alone, one name or accession per line. + + Blank lines and anything after a # are ignored, so the file can say why each entry is + there. Accessions may carry the full asmId with its assembly-name suffix, as the GenArk + orderLists do, and only the accession part is used, so + kent/src/hg/makeDb/doc/hprcAsmHub/hprc.orderList.tsv can be handed over unchanged. + """ + if not fname: + return set() + skip = set() + lineCount = 0 + for line in open(fname, encoding="latin1"): + line = line.split("#")[0].strip() + if not line: + continue + name = line.split()[0].split("\t")[0] + lineCount += 1 + skip.add(name) + # GCA_018466835.2_HG02257_mat_hprc_f2 -> GCA_018466835.2 + parts = name.split("_") + if len(parts) > 2 and accRe.match("%s_%s" % (parts[0], parts[1])): + skip.add("%s_%s" % (parts[0], parts[1])) + # len(skip) counts both name forms of an accession, so report the lines read + logging.info("Read %d assemblies to skip from %s" % (lineCount, fname)) + return skip + def getGenArkTaxIdDbs(onlyDbs, options): """ build the plan from the GenArk assembly list instead of from dbDb. dbDb holds the few thousand assemblies the browser serves directly; GenArk holds about fifty thousand. Most of those are not worth running: UniProt annotates most species barely at all, and an assembly whose organism has a handful of proteins cannot produce a useful track. --minProteins is the filter, and it wants to be high to start with. A bacterial proteome is only a few thousand proteins, so a threshold of a few thousand pulls in the whole of the bacteria clade; 20000 keeps it to the large eukaryotic proteomes. refs #38300 """ accToTax = readGenArkList(options.genArkList) counts = taxonProteinCounts(options.uniprotDir) minProt = options.minProteins + skipAsm = readSkipList(options.skipList) taxIdToDbs = defaultdict(list) - skippedTaxa, skippedAsm, noHub = set(), 0, 0 + skippedTaxa, skippedAsm, noHub, blacklisted = set(), 0, 0, 0 for acc, taxId in sorted(accToTax.items()): if onlyDbs is not None and acc not in onlyDbs: continue + if acc in skipAsm: + blacklisted += 1 + continue if counts.get(taxId, 0) < minProt: skippedTaxa.add(taxId) skippedAsm += 1 continue if genArkAccDir(acc) is None or not isdir(genArkAccDir(acc)): # listed but not served from /gbdb/genark here, nothing to build against noHub += 1 continue taxIdToDbs[taxId].append(acc) logging.info("GenArk plan: %d assemblies over %d taxa with at least %d UniProt proteins" % \ (sum(len(v) for v in taxIdToDbs.values()), len(taxIdToDbs), minProt)) logging.info("GenArk plan: skipped %d assemblies over %d taxa below the protein threshold, " - "and %d not present under %s" % (skippedAsm, len(skippedTaxa), noHub, genArkRoot)) + "%d on the skip list, and %d not present under %s" % \ + (skippedAsm, len(skippedTaxa), blacklisted, noHub, genArkRoot)) return taxIdToDbs def getTaxIdDbs(onlyDbs, options): """ return a dict with taxonId -> list of most recent dbs (e.g. for human, it's hg19 and hg38) if onlyDbs is a set, keep only these dbs. If onlyDbs is None, remove all nonAutoDbs. """ if options.genArkList: return getGenArkTaxIdDbs(onlyDbs, options) + skipAsm = readSkipList(options.skipList) query = "select taxId, name, nibPath from dbDb where active=1 order by orderKey;" rows = runQuery("hgcentral", query, usePublic=True) taxIdToDbs = defaultdict(list) # Take the first assembly of each kind for every taxon: the first GenArk one and the # first classic one. Taking only the very first, as this did, stopped updating the # classic assembly the moment a GenArk assembly for the same organism appeared and # sorted ahead of it. Cow is the example: ARS_UCD2.0 arrived, so bosTau9 silently # stayed on the UniProt release it had, while the new data went into a GenArk contrib # collection. Users are still on bosTau9. refs #38300 seen = defaultdict(set) for taxId, dbCode, nibPath in rows: if onlyDbs: if dbCode not in onlyDbs: continue else: - if dbCode in notAutoDbs: + if dbCode in notAutoDbs or dbCode in skipAsm: continue taxId = int(taxId) kind = "genark" if (nibPath or "").startswith("hub:") else "classic" if kind not in seen[taxId]: seen[taxId].add(kind) taxIdToDbs[taxId].append(dbCode) for taxId, manNames in manualTaxIdDbs.items(): if manNames is None: del taxIdToDbs[taxId] else: for dbCode in manNames: if onlyDbs: if dbCode not in onlyDbs: @@ -303,30 +336,35 @@ help="do not create the /gbdb/ symlinks") parser.add_option("", "--archiveDir", dest="archiveDir", action="store", \ default="/usr/local/apache/htdocs-hgdownload/goldenPath/archive/", help="Location of archive directory, default %default") parser.add_option("", "--taxonThreads", dest="taxonThreads", action="store", type="int", default=20, help="how many taxa to process at the same time, default %default. Raising this " "keeps the cluster busy: one taxon at a time leaves it idle during the long " "single-threaded steps between batches. Around 20 is a reasonable working value. Assemblies " "of the same taxon always run one after the other, they share a fasta file.") parser.add_option("", "--genArkList", dest="genArkList", action="store", help="build the plan from this GenArk assembly list instead of from dbDb, e.g. " "/hive/data/genomes/asmHubs/UCSC_GI.assemblyHubList.txt or a fresh copy of " "https://hgdownload.soe.ucsc.edu/hubs/UCSC_GI.assemblyHubList.txt . Use with " "--minProteins, which decides how much of GenArk is worth running.") + parser.add_option("", "--skipList", dest="skipList", action="store", + help="file of assemblies to leave alone, one name or accession per line, with " + "# comments. Entries may carry the full asmId, so a GenArk orderList works as " + "is: --skipList=kent/src/hg/makeDb/doc/hprcAsmHub/hprc.orderList.tsv drops the " + "464 HPRC haplotype assemblies, which take human from 572 assemblies to 108.") parser.add_option("", "--minProteins", dest="minProteins", action="store", type="int", default=1, help="skip a taxon with fewer than this many UniProt proteins, default %default, " "i.e. skip only the empty ones. UniProt annotates most species barely at all: of " "the 3974 taxa that have a GenArk assembly and a SwissProt entry, the median has " "four proteins and only 558 have more than a hundred. Raise this when running " "across many assemblies, to skip the ones that cannot produce a useful track.") parser.add_option("", "--mapQa", dest="mapQa", action="store_true", \ help="output some QA stats for the maps") parser.add_option("", "--db", dest="db", action="store_true", \ help="output the trackDb make command and uniprot <-> UCSC db assignments for debugging trackDb problems and showing which UCSC databases will be processed by the otto job and why") parser.add_option("", "--onlyFlip", dest="onlyFlip", action="store_true", \ help="After a run was aborted because of too many changes, now flip the files and ignore the size of the changes. Do not check for size increases anymore.") (options, args) = parser.parse_args()