90abb8f3321fc8cc0feb99ec196cd692dd55a5b0 max Fri Sep 25 14:42:10 2026 -0700 uniprot otto: size miniprot's memory from the proteins too, and survive a missing NCBI table Twelve taxa died with "[morecore] insufficient memory" inside their reservation. The formula sized the job from the genome alone, on a curve I measured by varying the genome (227, 424 and 811 Mb) while holding the protein set at 125 sequences. That characterised one input and generalised as though it had characterised the function: the real runs align whole proteomes, 78457 sequences for zebrafish, and the alignment work scales with that as well. A 0.46 Gb genome that the genome term put at 8 GB ran out of memory. There is now a term for the protein input and the floor is 16 GB rather than 8. Checked by rerunning the exact job that failed, at the new reservation: it completes and writes a 74 MB GFF. Reserving too much costs queue slots on a cluster that absorbs these jobs easily; reserving too little throws the assembly away after minutes of work. Six more taxa died on a missing ncbi/<taxId>.tsv. Not every organism is in NCBI's gene2refseq, so for some of them that file cannot exist. It is only used to resolve UniProt's Entrez cross-references, and the caller already checks whether it got anything, so a missing table now returns nothing and costs a little filtering accuracy rather than the assembly. Those taxa would also have sent every future run back to re-split 2.4 GB looking for a table that will never appear, so the split now leaves an empty file behind to record that it looked. refs #38300 diff --git src/hg/utils/otto/uniprot/doUniprot src/hg/utils/otto/uniprot/doUniprot index 140fd7bb29b..9e178188f3b 100755 --- src/hg/utils/otto/uniprot/doUniprot +++ src/hg/utils/otto/uniprot/doUniprot @@ -1717,48 +1717,58 @@ def twoBitFname(db): """ the assembly sequence, wherever this assembly keeps it. Only a real GenArk assembly keeps it in the hub directory. hs1 is served as a hub but keeps its 2bit at the classic /gbdb/<db>/<db>.2bit, same as a classic assembly. """ if isGenArk(db): return genArkTwoBit(db, genArkHubDir(db)) return "/gbdb/{db}/{db}.2bit".format(db=db) def chromSizesFile(db): " the chrom.sizes of this assembly, wherever it keeps it " if isGenArk(db): return genArkChromSizes(db, genArkHubDir(db)) return "/hive/data/genomes/%s/chrom.sizes" % db -def miniprotRamGb(genomeFa): - """ how much RAM to reserve for a miniprot run on this genome, in whole GB. - Measured on real assemblies: peak RSS is linear in genome size at about 11 GB per Gb - of sequence and is flat in the thread count, because the index is built once and shared - (2.49 GB at -t 1 versus 2.60 GB at -t 16 on the same 227 Mb). So this is sized off the - sequence only, with headroom, and never below a small floor. +def miniprotRamGb(genomeFa, protFa): + """ how much RAM to reserve for a miniprot run, in whole GB. + + Two terms, and the second one was learned the hard way. The genome index is linear in + sequence and flat in the thread count, measured at about 11 GB per Gb (2.49 GB at -t 1 + against 2.60 GB at -t 16 on the same 227 Mb, then 2.60 / 4.69 / 8.93 GB for 227 / 424 / + 811 Mb). But that was measured with 125 proteins held constant, so it missed that the + alignment work scales with the protein set as well: the real runs align whole proteomes, + 78457 sequences for zebrafish against 125 in the measurement, and a 0.46 Gb genome that + the genome term alone put at 5 GB died with "[morecore] insufficient memory" inside an + 8 GB reservation. + + So add a term for the protein input and keep the generous floor. Reserving too much + only costs queue slots on a cluster that absorbs these jobs easily; reserving too little + throws away the whole assembly after it has already run for minutes. refs #38300 """ - bases = os.path.getsize(genomeFa) # near enough, fasta is one byte per base plus headers - gb = int(bases / 1e9 * 14) + 1 - return max(gb, 8) + genomeBytes = os.path.getsize(genomeFa) # fasta is about one byte per base plus headers + protBytes = os.path.getsize(protFa) + gb = int(genomeBytes / 1e9 * 14) + int(protBytes / 1e9 * 100) + 1 + return max(gb, 16) def runMiniprotOnCluster(db, genomeFa, protFa, gffName, workDir): """ run one miniprot job on the parasol cluster. It has to be told both numbers: para's default RAM is the node's RAM divided by its CPU count, which for a 16 CPU job is far less than miniprot needs, and without -cpu parasol would pack more of these onto a node than it has cores for. """ - ram = miniprotRamGb(genomeFa) + ram = miniprotRamGb(genomeFa, protFa) jobDir = join(workDir, "cluster") if not isdir(jobDir): os.makedirs(jobDir, exist_ok=True) # a wrapper, so the jobList line stays free of the redirection and quoting that # parasol's job parser does not accept # Every path here has to be absolute. A parasol job runs with its working directory # set to the batch directory, not to the directory the pipeline runs in, so a path # like "fasta/7955.fa" simply does not exist from the job's point of view and miniprot # exits without writing anything. jobSh = join(jobDir, "runMiniprot.sh") with open(jobSh, "w") as ofh: ofh.write("#!/bin/sh\nset -e\n") ofh.write("%s -t %d --gff %s %s > $1\n" % \ (miniprotBin, miniprotThreads, abspath(genomeFa), abspath(protFa))) @@ -2227,33 +2237,43 @@ stats["version"] = version stats["createdDate"] = datetime.datetime.today().strftime('%Y-%m-%d') with open(mapDescFname, "w") as mapDescFh: json.dump(stats, mapDescFh, indent=4) def listMd5(arr): " return hex md5 given list of strings " import hashlib m = hashlib.md5() for a in arr: m.update(a.encode("utf8")) return m.hexdigest() def readGeneToRefSeq(taxId): - " return entrezId -> list of Refseq from tsv file " + """ return entrezId -> list of Refseq from tsv file, empty if we have none for this taxon. + + Not every organism is in NCBI's gene2refseq, so for some taxa this file cannot exist. + The caller only uses it to resolve UniProt's Entrez cross-references to transcripts and + already checks whether it got anything, so an empty answer costs a little filtering + accuracy rather than the assembly. It used to raise and take the taxon down. refs #38300 + """ g2refseq = {} fname = join(NCBIDIR, str(taxId)+".tsv") + if not isfile(fname): + logging.warning("taxon %s is not in the NCBI gene2refseq tables, so UniProt's Entrez " + "cross-references cannot be resolved for it" % taxId) + return g2refseq for line in open(fname): geneId, transList = line.rstrip("\n").split('\t') transList = transList.split(',') g2refseq[geneId] = transList return g2refseq def parseKeyVal(fname): " given a file with key<tab>val return a dict with key -> list of values " if fname is None: return None vals = defaultdict(list) for line in open(fname): key, val = line.rstrip("\n").split("\t") vals[key].append(val) return vals @@ -3121,30 +3141,40 @@ """ make sure there is a per-taxon NCBI gene -> RefSeq table for every taxon we will map through RefSeq models. Splitting is cheap next to downloading, and the tables are per taxon, so a plan with new taxa needs it even when the download is skipped. """ missing = [t for t in taxIdDbs if not isfile(join(NCBIDIR, "%s.tsv" % t))] if not missing: return if not isfile(ncbiGeneFname): errAbort("%s is missing and the download was skipped, so the per-taxon NCBI tables " "for %d taxa cannot be built. Run without --skipDownload." % (ncbiGeneFname, len(missing))) logging.info("%d of %d taxa have no NCBI gene table yet, splitting %s" % \ (len(missing), len(taxIdDbs), ncbiGeneFname)) splitGeneRefseq(ncbiGeneFname, NCBIDIR, taxIdDbs) + # Some organisms are not in gene2refseq at all, so the split cannot produce a table for + # them and they would send every future run back to re-split 2.4 GB. Leave an empty file + # to record that we looked; readGeneToRefSeq treats it as "nothing known", which is true. + stillMissing = [t for t in missing if not isfile(join(NCBIDIR, "%s.tsv" % t))] + for taxId in stillMissing: + open(join(NCBIDIR, "%s.tsv" % taxId), "w").close() + if stillMissing: + logging.info("%d taxa are not in gene2refseq at all, wrote empty tables for them" % \ + len(stillMissing)) + def downloadAndSplitNcbi(taxIdDbs): " download and split the NCBI genes file with a mapping NCBI gene -> RefSeq transcripts " downloadNcbi() splitGeneRefseq(ncbiGeneFname, NCBIDIR, taxIdDbs) def delFlag(): global flagFname if isfile(flagFname): os.remove(flagFname) def main(): global flagFname args, options = parseArgs()