55a768d2e0ced94dc3ba7ab322b24daa67ea3570 max Fri Sep 11 10:09:16 2026 -0700 uniprot otto: resolve the 2bit and chrom.sizes for every kind of assembly The first GenArk run stopped on hs1 with "expected exactly one chrom.sizes file in /gbdb/hs1/hubs, found 0". hs1 is served as a hub but keeps its 2bit at /gbdb/hs1/hs1.2bit and its chrom.sizes at /hive/data/genomes/hs1/chrom.sizes, exactly where a classic assembly keeps them; only a real GenArk assembly keeps them in the hub directory. twoBitFname now keys on isGenArk rather than on being a hub, and a matching chromSizesFile does the same, so both are resolved the same way everywhere. The transcript building still uses the bigBed path for any hub assembly, because hs1 has no MySQL tables either, but it asks these two helpers for the supporting files. Checked that all four resolve to files that exist: hs1 and hg38 to the classic locations, mPanPan1_v2.0 and GRCz12ab to their hub directories. Also fixed the README, which documented the option to limit the run as --onlyDbs. It is --dbs; --onlyDbs is the internal dest name and the run fails with "no such option". refs #38300 diff --git src/hg/utils/otto/uniprot/doUniprot src/hg/utils/otto/uniprot/doUniprot index 8fe7dd65e6e..2735e6d5f7b 100755 --- src/hg/utils/otto/uniprot/doUniprot +++ src/hg/utils/otto/uniprot/doUniprot @@ -1502,36 +1502,44 @@ fnames = [ (join(tabDir,"swissprot.%d.tab" % taxId)) ] if doTrembl: fnames.append( (join(tabDir,"trembl.%d.tab" % taxId)) ) ret = {} for fname in fnames: for row in iterTsvRows(open(fname)): for acc in getAllIsoAccs(row): ret[acc] = row return ret def twoBitFname(db): - " the assembly sequence, wherever this assembly keeps it " - hubDir = genArkHubDir(db) - if hubDir is not None: - return genArkTwoBit(db, hubDir) + """ 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. """ 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 @@ -1872,33 +1880,33 @@ 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") cdsName = join(dbDir, "transcripts.cds") run(["bigGenePredToGenePred", geneBb, gpName]) - run(["genePredToFakePsl", "-chromSize=%s" % genArkChromSizes(db, hubDir), + run(["genePredToFakePsl", "-chromSize=%s" % chromSizesFile(db), "noDb", gpName, pslName, cdsName]) - run(["getRnaPred", "-genomeSeqs=%s" % genArkTwoBit(db, hubDir), + run(["getRnaPred", "-genomeSeqs=%s" % twoBitFname(db), "noDb", gpName, "all", faName]) logging.info("Created %s and %s for %s/%s from %s" % (faName, pslName, db, geneTable, geneBb)) return faName, pslName if geneTable=="ncbiRefSeq": origFa = "/gbdb/%s/ncbiRefSeq/seqNcbiRefSeq.rna.fa" % db shutil.copyfile(origFa, faName) pslTable = "ncbiRefSeqPsl" query = "SELECT "+pslFields+" from "+pslTable+" WHERE tName NOT LIKE '%_hap%' AND tName not like '%_alt%' AND tNAME NOT LIKE '%_fix%'" runQueryToFile(db, query, pslName) elif geneTable=="refGene": # ChrisL told me where these files can be found @@ -2277,35 +2285,31 @@ continue 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 - hubDir = genArkHubDir(db) - if hubDir is not None: - chromSizesFname = genArkChromSizes(db, hubDir) - else: - chromSizesFname = "/hive/data/genomes/%s/chrom.sizes" % db + chromSizesFname = chromSizesFile(db) accToMeta = parseRecordMeta(tabDir, taxId, doTrembl) dbMapDir = makeSubDir(mapDir, db) mapDescFname = join(dbMapDir, "liftInfo.json") mapFname, geneTable, mapMd5, accToMapSource, accToTrans, isNewMap = \ makeProteinGenomePsl(taxId, fullFaFname, tabFnames, db, options.force, dbMapDir, mapDescFname) dbBigBedDir = makeSubDir(bigBedDir, db) convMapToBigPsl(fullFaFname, mapFname, dbBigBedDir, accToDb, accToMeta, accToMapSource, accToTrans, chromSizesFname, doTrembl) #shutil.copy(mapDescFname, mapDir) if options.onlyMap: