45eb30c0414d20710274af10680660b9d00f04b8 max Mon Oct 5 01:23:09 2026 -0700 uniprot otto: raise miniprot's memory when a job runs out, and skip an assembly with no transcripts Two of the eight taxa that failed the GenArk run died with "[morecore] insufficient memory" inside their reservation, so the sizing was wrong a second time. Measuring says it cannot be sized from the inputs at all: 0.46 Gb of genome with 78457 proteins fits in 16 GB, while 0.86 Gb with 50305 proteins peaks at 43.3 GB. Half the protein and twice the genome needs at least 2.7 times the memory, so no sum of a genome term and a protein term fits both. What drives it is how much alignment work the sequence generates, which the file sizes do not show. So stop predicting. The formula stays as a first guess and the batch is run again with four times the reservation when a job ran out of memory, up to three tries. A batch that failed for any other reason is not retried, so the 2.5 Gb genome that segfaults is still only run once. Checked the detector against both real batches: the one that ran out of memory reads as such, the one that segfaulted does not. para records the raw wait status, so 134 means the abort miniprot makes when malloc fails. The reservation lives in the batch, so freeing it and deleting the state files is what lets para make use a new one; checked that freeBatch does not prompt and that a cleared batch runs again. Three more of the eight had a gene track that produced no transcripts at all. BLAST was pointed at an empty database and all 929 of its jobs crashed, taking the taxon down after a long detour through the cluster. There is nothing to align against, so that assembly is now skipped the way one with no alignments already is. refs #38300 diff --git src/hg/utils/otto/uniprot/doUniprot src/hg/utils/otto/uniprot/doUniprot index b0a0b4b6b8b..e65937b79a4 100755 --- src/hg/utils/otto/uniprot/doUniprot +++ src/hg/utils/otto/uniprot/doUniprot @@ -1732,80 +1732,134 @@ """ 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, 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 + """ a first guess at how much RAM to reserve for a miniprot run, in whole GB. + + Only a guess. Two attempts to predict this from the input sizes were both wrong, and the + measurements say the quantity is not predictable from them at all: + + genome proteins peak RSS + 0.46 Gb 78457 (66 MB) at most 16 GB, the job completed in that + 0.86 Gb 50305 (37 MB) 43.3 GB, measured + + The second has half the protein and twice the genome, and needs at least 2.7 times the + memory, so no sum of a genome term and a protein term fits both. What actually drives it + is how much alignment work the sequence generates - repeats and paralogs give many + candidate alignments to hold at once - and that is not visible in the file sizes. + + So this stays a cheap starting estimate and runMiniprotOnCluster() raises it when a job + runs out of memory, rather than a third guess at a coefficient. refs #38300 """ 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 paraBatchOutOfMemory(jobDir): + """ did a job in this finished parasol batch die for lack of memory? + + para.results holds the raw wait status in its first column, so status >> 8 is the exit + code. miniprot prints "[morecore] insufficient memory" and aborts, which arrives as 134 + (128 + SIGABRT). Reading the file is local and cheap; para problems would have to be + asked the hub. + """ + resultsFname = join(jobDir, "para.results") + if not isfile(resultsFname): + return False + for line in open(resultsFname): + fields = line.split() + if not fields: + continue + try: + status = int(fields[0]) + except ValueError: + continue + if status >> 8 == 134: + return True + return False + +def clearParaBatch(jobDir): + """ forget a finished batch so it can be run again with a different reservation. + + The reservation is recorded in the batch, so para make on top of an existing one does not + change it. freeBatch tells the hub to drop it and the local state files have to go too. + """ + run("cd %s && para freeBatch" % jobDir, ignoreErr=True) + for fname in ("batch", "batch.bak", "para.results", "para.bookmark"): + path = join(jobDir, fname) + if isfile(path): + os.remove(path) + 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, 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))) os.chmod(jobSh, 0o755) jobList = join(jobDir, "jobList") with open(jobList, "w") as ofh: ofh.write("%s {check out exists %s}\n" % (abspath(jobSh), abspath(gffName))) - logging.info("%s: miniprot on the cluster, -cpu=%d -ram=%dg" % (db, miniprotThreads, ram)) - run("cd %s && para make -cpu=%d -ram=%dg jobList" % (jobDir, miniprotThreads, ram)) + # The first reservation is only an estimate (see miniprotRamGb), so when a job dies for + # lack of memory, raise it and run the batch again rather than losing the assembly. Three + # tries reach eight times the estimate, which covers the worst case measured so far by a + # wide margin. Anything that fails for another reason stops here, so a genome that + # segfaults does not get run four times over. refs #38300 + for attempt in range(3): + tryRam = ram * (4 ** attempt) + logging.info("%s: miniprot on the cluster, -cpu=%d -ram=%dg%s" % \ + (db, miniprotThreads, tryRam, "" if attempt == 0 else " (attempt %d)" % (attempt+1))) + ret = run("cd %s && para make -cpu=%d -ram=%dg jobList" % \ + (jobDir, miniprotThreads, tryRam), ignoreErr=True) + if ret == 0: + return + if not paraBatchOutOfMemory(jobDir): + break + logging.warning("%s: miniprot ran out of memory in %dg, trying again with more" % (db, tryRam)) + clearParaBatch(jobDir) + + errAbort("%s: miniprot on the cluster failed, see %s" % (db, jobDir)) def miniprotProteins(fullFaFname, db, mapFname, stats): """ align the UniProt proteins straight to the genome with miniprot and write a PSL. Used when the assembly has no gene models worth mapping through. This replaced both "blat -q=prot" and mapping through Augustus: Augustus is ab initio, so going through it stacks its errors on top of ours, and the BLAT protein search is very slow on a big genome. """ workDir = mapFname+".miniprot.tmp" if not isdir(workDir): os.makedirs(workDir, exist_ok=True) # miniprot reads fasta, not 2bit genomeFa = join(workDir, "genome.fa") run(["twoBitToFa", twoBitFname(db), genomeFa]) @@ -2304,30 +2358,40 @@ """ # use MD5s of the input files to determine if the protein -> genome map has to be rebuilt protMd5 = fastaMd5(protFa) stats = OrderedDict() stats["user"] = os.getlogin() stats["taxId"] = taxId stats["db"] = db geneTable = findBestGeneTable(db) metaMd5 = manyTabMd5(tabFnames) if geneTable!="miniprot": transcriptFa, transcriptPsl = makeTranscriptFiles(db, geneTable, taxId, mapDir) + # An assembly can carry a gene track that yields no transcripts at all, and three of + # them did. BLAST was then pointed at an empty database and every one of its 929 jobs + # crashed, which took the taxon down after a long detour through the cluster. There is + # nothing to align against, so say so and skip the assembly the way one with no + # alignments is skipped. refs #38300 + if os.path.getsize(transcriptFa) == 0: + logging.warning("%s: the %s track produced no transcript sequences, so there is " + "nothing for the proteins to be aligned to, skipping this assembly" % \ + (db, geneTable)) + return None, geneTable, None, {"default":"best"}, None, True transMd5 = fastaMd5(transcriptFa) pslMd5 = tabMd5(transcriptPsl) else: # this usually only happens on viral genomes or other weird cases that don't have a single gene model track # here we align the proteins directly onto the genome with miniprot transMd5, pslMd5 = "directMiniprot-noTranscripts", "directMiniprot-noPsl" transcriptFa, transcriptPsl = None, None geneToRefSeqs = None selectFname = None protToTrans = None if geneTable in ["miniprot"]: protMapSource = {"default":"direct"} elif geneTable in ["augustusGene"]: