da869954bf88c67e25a8b68cfd49ad50c5cec7d0 max Mon Sep 14 05:58:57 2026 -0700 uniprot otto: optionally work on several taxa at once The assemblies were processed strictly one after another, and that leaves the cluster mostly idle. Each assembly submits its BLAST batch, waits for it to drain, and then spends the best part of an hour in the single-threaded pslReps that follows before the next assembly submits anything. Measured on the hs1 batch: 332 hours of CPU finished in 30 minutes of wall clock, a speedup of about 660, and then 58 minutes of one core concatenating and filtering 34335 PSL files. With 119 assemblies in the plan and a 1568 CPU cluster, that ordering costs about a day. --taxonThreads=N runs N taxa at a time. It defaults to 1, so nothing changes unless it is asked for; 4 to 6 is a reasonable range. Taxa, not assemblies. The assemblies of one taxon share fasta/<taxId>.fa, which is rebuilt at the start of each taxon, so two threads on one taxon would race on it. Below that everything is per-assembly - protToGenome/<db>, bigBed/<db>, the cluster batch directory - so separate taxa do not share files. os.makedirs calls are now exist_ok, since two taxa starting together can both find a directory missing. A failing worker is escalated, not swallowed. run() reports a failed command with sys.exit(), which raises SystemExit; in a worker thread that would kill only that thread and leave the run looking successful, which is the failure mode this pipeline has a long history of. runTaxa catches BaseException per taxon, names each one that failed, and aborts the run at the end. Checked: a worker calling sys.exit aborts the run, a worker raising an ordinary exception aborts the run, both after the other taxa have still been attempted, an all-good set returns normally, --taxonThreads=1 keeps the old sequential order, and --dbs still selects which taxa run. refs #38300 diff --git src/hg/utils/otto/uniprot/doUniprot src/hg/utils/otto/uniprot/doUniprot index 38a0f6d6631..2379ed0671c 100755 --- src/hg/utils/otto/uniprot/doUniprot +++ src/hg/utils/otto/uniprot/doUniprot @@ -1,19 +1,20 @@ #!/usr/bin/env python3 import urllib.request, urllib.error, urllib.parse, time, os, atexit import datetime, optparse, sys, logging, subprocess, glob, shutil import sys, gzip, logging,re, json +import concurrent.futures from urllib.parse import urlparse from os.path import * from os import system, makedirs, mkdir, remove, listdir from shutil import move from subprocess import PIPE from collections import defaultdict, namedtuple, Counter, OrderedDict # main driver script for uniprot updates # This script has evolved over a long time. Parts of it do not have an optimal structure. # One reason is that the script used to parse only mutations. Even normal annotations # are still put into a structure that resembles mutations. # the script first uses makeUniProtPsl.sh to create a mapping from UniProt Sequence to genome as a PSL file. @@ -188,30 +189,36 @@ parser.add_option("-f", "--faDir", dest="faDir", action="store", \ help="directory for full fasta files, default %default", default="fasta") parser.add_option("-m", "--mapDir", dest="mapDir", action="store", \ help="directory for pslMap (~liftOver) psl files, default %default", default="protToGenome") parser.add_option("-b", "--bigBedDir", dest="bigBedDir", action="store", \ help="directory for bigBed files, one subdirectory per db will be created", default="bigBed") parser.add_option("", "--force", dest="force", action="store_true", \ help="skip the check of differences against the previous version") parser.add_option("", "--onlyLinks", dest="onlyLinks", action="store_true", \ help="only create the /gbdb/ symlinks") parser.add_option("", "--skipLinks", dest="skipLinks", action="store_true", \ 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=1, + 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. 4 to 6 is a reasonable range. Assemblies " + "of the same taxon always run one after the other, they share a fasta file.") 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() if options.db: taxIdDbs = getTaxIdDbs(None) print("-- Current TaxId<->database assignment:") for key, val in taxIdDbs.items(): print (key, val) allDbs = [] @@ -1540,62 +1547,62 @@ sequence only, with headroom, and never below a small floor. """ 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) 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) jobDir = join(workDir, "cluster") if not isdir(jobDir): - os.makedirs(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)) 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) + os.makedirs(workDir, exist_ok=True) # miniprot reads fasta, not 2bit genomeFa = join(workDir, "genome.fa") run(["twoBitToFa", twoBitFname(db), genomeFa]) gffName = join(workDir, "miniprot.gff") runMiniprotOnCluster(db, genomeFa, fullFaFname, gffName, workDir) os.remove(genomeFa) # Name each alignment after the UniProt accession it came from. miniprot calls them # MP000001 and puts the accession in Target=, and its ##PAF meta lines and capitalised # attributes make gff3ToGenePred unhappy, so both are cleaned up here. The ID has to # stay unique for GFF3, hence the .N suffix, which comes off again after the # conversion so that qName is the bare accession. namedName = join(workDir, "named.gff") @@ -1874,31 +1881,31 @@ stats["notFound_examples"] = list(notFound[:10]) stats["notFoundClasses"] = dict(prefixCounts) stats["outPairCount"] = pairCount stats["outPairUniprot"] = len(outUpIds) stats["outPairTrans"] = len(outTransIds) return selectFname, stats, protMapSource def makeTranscriptFiles(db, geneTable, taxId, dbDir): """ given a db and a gene table, create two files, one for the transcript fasta sequences and one for the transcript -> genome PSL. Return tuple of the file names ALWAYS REMOVES _ALT/_FIX/_HAP results from the PSL! """ #dbDir = join(mapDir, db+"_"+geneTable) #if not isdir(dbDir): - #os.makedirs(dbDir) + #os.makedirs(dbDir, exist_ok=True) faName = join(dbDir, "transcripts.fa") pslName = join(dbDir, "transcripts.psl") pslFields = "matches,misMatches,repMatches,nCount,qNumInsert,qBaseInsert,tNumInsert," \ "tBaseInsert,strand,qName,qSize,qStart,qEnd,tName,tSize,tStart,tEnd,blockCount,blockSizes,qStarts,tStarts" hubDir = genArkHubDir(db) if hubDir is not None: # GenArk assembly: no MySQL tables at all, the gene models are a bigGenePred in the # hub. bigGenePredToGenePred gets us back to a genePred, and both of the tools that # follow can work off the hub's own chrom.sizes and 2bit rather than chromInfo, so # nothing here needs a database. geneBb = findBestGeneBigBed(db, hubDir)[1] gpName = join(dbDir, "transcripts.gp") @@ -2037,32 +2044,32 @@ 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) 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 blat the proteins directly onto the genome - transMd5, pslMd5 = "directBlat-noTranscripts", "directBlat-noPsl" + # 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"]: protMapSource = {"default":"best"} else: if geneTable in ["refGene", "ncbiRefSeq"]: geneToRefSeqs = readGeneToRefSeq(taxId) selectFname, stats, protMapSource = buildSelectFile(tabFnames, mapDir, taxId, db, \ geneTable, transcriptFa, geneToRefSeqs, metaMd5, transMd5, stats) @@ -2148,45 +2155,45 @@ chainFname = join(bigPslDir, "unipToGenome.over.chain") cmd = ["pslToChain",mapSwapFname, chainFname] run(cmd) cmd = ["gzip","-f",chainFname] run(cmd) logging.info("Created %s" % liftFname) chainFname += ".gz" os.remove(mapSwapFname) def makeSubDir(parent, child): " mkdir child under parent and return " newDir = join(parent, child) if not isdir(newDir): logging.info("Making dir: %s" % newDir) - os.makedirs(newDir) + os.makedirs(newDir, exist_ok=True) return newDir def downloadUniprot(uprotDir): " if necessary, download UniProt using LFTP in an atomic way, update release info file " localFname = join(uprotDir, "uniprot_sprot.xml.gz") if not isNewer(upUrl, localFname): logging.info("files at %s are not newer than file in %s, nothing to do. Specify -l to skip this check." % (upUrl, localFname)) sys.exit(0) # use lftp to update uniprot and download only changed files # in late 2020, pget started to trigger errors, so removed -P 3 and --use-pget-n=4 tmpDir = join(uprotDir, "download.tmp") if not isdir(tmpDir): - os.makedirs(tmpDir) + os.makedirs(tmpDir, exist_ok=True) # the uniprot FTP server is too slow and too unreliable. Switching to EBI's which has rsync !! count = 0 while True: #cmd = "rsync -av --partial rsync://ftp.ebi.ac.uk/pub/databases/uniprot/current_release/knowledgebase/complete/ %s/" % tmpDir # ftp.uniprot.org (USA,, PIR) # ftp.ebi.ac.uk (UK, EBI) # ftp.expasy.org (Switzerland, SIB) cmd = 'lftp ftp://ftp.expasy.org/databases/uniprot/current_release/knowledgebase/complete/ -e "lcd %s && mirror . -P 3 --use-pget-n=4 --exclude-glob *.dat.gz && exit"' % tmpDir ret = run(cmd, ignoreErr=True) if ret==0: break logging.info("Problem with rsync, waiting for 5 hours, then retrying") @@ -2285,49 +2292,51 @@ # takes 2-3 days on hgwdev - big data and XML don't mix... cmd="%s/uniprotToTab %s %s %s --trembl" % (myDir, uprotDir, taxIdStr, tabDir) run(cmd) writeReleaseString(uprotDir, tabDir) # if the uniprot update changed the sequences, update the corresponding pslMap files of that genome logging.info("checking/creating pslMap files") run("mkdir -p %s" % mapDir) # get the uniProt version for the trackVersion table that we will update later relFname = join(tabDir, "version.txt") versionString = open(relFname).read() - for taxId, dbs in taxIdDbs.items(): - if onlyDbs is not None and len(set(dbs).intersection(onlyDbs))==0: - continue + runTaxa(taxIdDbs, onlyDbs, options, tabDir, faDir, mapDir, bigBedDir, doTrembl) + +def oneTaxon(taxId, dbs, onlyDbs, options, tabDir, faDir, mapDir, bigBedDir, doTrembl): + " do all the assemblies of one taxon. Safe to run several of these at once, see runTaxa " + if True: logging.info("Working on taxon ID %d" % taxId) faFnames = [ ("swissprot", join(tabDir,"swissprot.%d.fa.gz" % taxId)), ] if doTrembl: faFnames.append( ("trembl", join(tabDir,"trembl.%d.fa.gz" % taxId)) ) fullFaFname = join(faDir, str(taxId)+".fa") accToDb = concatFiles(faFnames, fullFaFname) if os.path.getsize(fullFaFname)==0: logging.warn("File %s is empty. Taxon ID %s does not have any UniProt annotations." \ " Skipping this organism." % (fullFaFname, taxId)) - continue + return tabFnames = ["tab/swissprot.%d.tab" % taxId] if doTrembl: tabFnames.append( "tab/trembl.%d.tab" % taxId ) # find the best gene table for each database, create a mapping # protein -> genome and lift the uniprot annotations to bigBed files for db in dbs: if onlyDbs is not None and db not in onlyDbs: continue logging.debug("Annotating assembly %s" % db) # a GenArk assembly has no /hive/data/genomes directory, its chrom.sizes # lives in the hub chromSizesFname = chromSizesFile(db) @@ -2344,30 +2353,74 @@ convMapToBigPsl(fullFaFname, mapFname, dbBigBedDir, accToDb, accToMeta, accToMapSource, accToTrans, chromSizesFname, doTrembl) #shutil.copy(mapDescFname, mapDir) if options.onlyMap: continue annotFnames = ["tab/swissprot.%d.annots.tab" % taxId] if doTrembl: annotFnames.append( "tab/trembl.%d.annots.tab" % taxId ) uniprotLift(fullFaFname, annotFnames, chromSizesFname, mapFname, dbBigBedDir, accToDb, accToMeta, options) shutil.copyfile(mapDescFname, join(dbBigBedDir, "liftInfo.json")) +def runTaxa(taxIdDbs, onlyDbs, options, tabDir, faDir, mapDir, bigBedDir, doTrembl): + """ work through the taxa, optionally several at a time. + + One taxon at a time leaves the cluster mostly idle: each assembly submits its BLAST + batch, waits for it to drain, and then spends the best part of an hour in the + single-threaded pslReps that follows before the next assembly submits anything. + Running several taxa at once keeps jobs queued while others are in that serial phase. + + Taxa, not assemblies: the assemblies of one taxon share fasta/<taxId>.fa, which is + rebuilt at the start of each taxon, so two threads on the same taxon would race on it. + Everything below that point is per-assembly - protToGenome/<db>, bigBed/<db>, and the + cluster batch directory - so separate taxa do not touch the same files. + """ + todo = [(taxId, dbs) for taxId, dbs in taxIdDbs.items() + if onlyDbs is None or len(set(dbs).intersection(onlyDbs))>0] + + if options.taxonThreads <= 1 or len(todo) <= 1: + for taxId, dbs in todo: + oneTaxon(taxId, dbs, onlyDbs, options, tabDir, faDir, mapDir, bigBedDir, doTrembl) + return + + logging.info("Working on %d taxa, %d at a time" % (len(todo), options.taxonThreads)) + failed = [] + with concurrent.futures.ThreadPoolExecutor(max_workers=options.taxonThreads) as pool: + futures = {} + for taxId, dbs in todo: + fut = pool.submit(oneTaxon, taxId, dbs, onlyDbs, options, tabDir, faDir, + mapDir, bigBedDir, doTrembl) + futures[fut] = taxId + for fut in concurrent.futures.as_completed(futures): + taxId = futures[fut] + try: + fut.result() + except BaseException as ex: + # BaseException, not Exception: run() reports a failed command with + # sys.exit(), which raises SystemExit. In a worker thread that would + # otherwise kill just that thread and leave the run looking successful, + # which is exactly the kind of silent failure this pipeline has a history of. + logging.error("taxon %s failed: %s: %s" % (taxId, type(ex).__name__, ex)) + failed.append(taxId) + + if failed: + errAbort("%d of %d taxa failed: %s" % (len(failed), len(todo), ", ".join(str(t) for t in sorted(failed)))) + def makeSymlink(target, linkName): assert(".new" not in linkName) # make sure that this bug never happens again if target is None: logging.error("internal error, but not stopping: Symlink to a None target?") return targetPath = abspath(target) if not isfile(targetPath): logging.error("Cannot symlink: %s does not exist" % str(target)) return if isfile(linkName): currTarget = os.readlink(linkName) @@ -2638,31 +2691,31 @@ def installGenArkContrib(db, hubDir, dbBigBedDir, versionString, shortVersion): """ write this assembly's bigBeds and a trackDb into the GenArk contrib collection. A GenArk assembly is not served out of /gbdb/<db>/, so the /gbdb symlinks and the trackDb .ra stanza that classic assemblies get do not apply. This only fills in the collection under contrib/uniprot/. It deliberately does NOT symlink anything into the GenArk build directories and does not touch any hub.txt: installing the collection is a separate, deliberate step, done with "genark addContrib uniprot" when someone decides to, not something a data update should do behind your back. """ acc = basename(hubDir.rstrip("/")) collDir = join(genArkContribRoot, genArkContribName, acc) if not isdir(collDir): - os.makedirs(collDir) + os.makedirs(collDir, exist_ok=True) count = 0 for bbName in sorted(glob.glob(join(dbBigBedDir, "*.bb"))): # skip the historical per-geneTable files, only the current track files belong here if not basename(bbName).startswith("unip"): continue shutil.copyfile(bbName, join(collDir, basename(bbName))) count += 1 tdbFname = join(collDir, "trackDb.txt") with open(tdbFname, "w") as ofh: ofh.write(makeContribTrackDb(versionString, shortVersion)) logging.info("%s: wrote %d bigBeds and a trackDb to %s" % (db, count, collDir)) return collDir @@ -2685,31 +2738,31 @@ def copyToArchive(bigBedDir, archRoot, shortVersion, onlyDbs): " make copies of the track files under the archiveDir and adapt the 'current' symlink " for db in os.listdir(bigBedDir): if onlyDbs and db not in onlyDbs: continue if isGenArk(db): # goldenPath/archive/<db>/ is a classic-assembly download path; a GenArk # assembly is published under its accession in the hubs tree instead, so # archiving one here would just create a directory nothing can reach logging.debug("%s is a GenArk assembly, not archiving under goldenPath" % db) continue archDir = join(archRoot, db, "uniprot", shortVersion) if not isdir(archDir): - os.makedirs(archDir) + os.makedirs(archDir, exist_ok=True) inDir = join(bigBedDir, db) count = 0 for inFname in glob.glob(join(inDir, "*")): shutil.copyfile(inFname, join(archDir, basename(inFname))) count += 1 logging.info("Archive: Copied %d files from %s to %s" % (count, inDir, archDir)) makeTrackDb(archDir, shortVersion) currLink = join(archRoot, db, "uniprot", "current") if islink(currLink): os.remove(currLink) os.symlink(archDir, currLink)