89e1e4c6ebefe29b7219fc98417b73b6a1755f05 max Fri Sep 11 10:07:02 2026 -0700 uniprot otto: handle hub assemblies that are not GenArk, and bring hs1 back in Not every assembly with a "hub:" nibPath is a GenArk assembly. hs1 is served as a hub from /gbdb/hs1/hubs but keeps its track files in /gbdb/hs1/ exactly like a classic assembly, so treating every hub as GenArk put its gene model search in the wrong place and would have sent its output to a contrib collection it does not belong in. Split the two questions. genArkHubDir still answers "is this served as a hub, and from where", and a new isGenArk asks the narrower question the output routing actually cares about: is this under /gbdb/genark with the sharded layout. The /gbdb symlinks, the version.txt symlink, the goldenPath archive and the contrib collection now all key on isGenArk, so hs1 gets the classic treatment. The gene model search now takes several globs per source and looks through both the hub directory and /gbdb//, which covers both layouts: a GenArk assembly keeps its models under bbi/ with CAT in a contrib collection, hs1 keeps them in /gbdb/hs1//. That makes hs1 work, so it comes out of notAutoDbs. It was excluded for having "no good gene model" and because the BLAT protein search took forever; it actually has CAT genes, 234903 of them in catLiftOffGenesV1, which is the best source we look for, and BLAT is gone anyway. Resolutions checked: hs1 -> catGenes from /gbdb/hs1/catLiftOffGenesV1, not GenArk; mPanPan1_v2.0 -> ncbiRefSeq from its hub bbi/, GenArk; GRCz12ab -> miniprot, GenArk; hg38 and panPan3 unchanged on the MySQL path. refs #38300 diff --git src/hg/utils/otto/uniprot/doUniprot src/hg/utils/otto/uniprot/doUniprot index 66d6f243c5b..8fe7dd65e6e 100755 --- src/hg/utils/otto/uniprot/doUniprot +++ src/hg/utils/otto/uniprot/doUniprot @@ -43,32 +43,30 @@ 7955 : ["danRer7", "danRer10"], # zebrafish 559292 : ["sacCer3"], # yeast 9940 : ["oviAri3"], # sheep 9823 : ["susScr3"], # pig 2697049: ["wuhCor1"]} # see also covidCheck.sh notAutoDbs = [ # got errors due to short contigs in galGal6, just skipping "galGal6", # do not run on wuhCor1 - wuhCor1 has a special daily otto job in /hive/data/genomes/wuhCor1/bed/uniprot # because UniProt releases SARS-CoV-2 much more often "wuhCor1", # no monkeypox yet "mpxvRivers", -# do not run on hs1, has no good gene model and BLAT protein takes forever -"hs1", # sea urchin doesn't have any normal gene model track, not even Augustus "strPur2" ] # the script will stop if the difference between the last uniprot version and the current version is > than # but this check is only done for some important assemblies. For fringe assemblies, the numbers can easily # change 10x, as trembl predictions are the majority of annotations maxDiffVersions = 5.0 checkDbs = ["hg19", "hg38", "mm10", "rn6", "danRer10", "danRer7", "ce11"] TREMBLCOLOR="0,150,250" # light blue SWISSPCOLOR="12,12,120" # dark blue asDir = dirname(__file__) @@ -995,87 +993,115 @@ logging.warning("%s: Gene %s appears on %d chromosomes - this assembly has multi-mappers." % (db, gene, count)) return False return True # GenArk assemblies keep their gene models as bigBed files in the hub, not as MySQL tables, # so they need their own search. In descending order of preference: # CAT - Comparative Annotation Toolkit, shipped as a contrib collection # (track hprcCatGenes, contrib/hprc2annot/catGenes.bb). Best annotation we # have on the assemblies that carry it. # ncbiRefSeq - RefSeq, the GenArk RefSeq gene track. On most GCF assemblies. # ncbiGene - the annotation the submitter sent to GenBank with the assembly. # augustus - ab initio prediction, present almost everywhere, so it is the last resort. # Each entry is (name, glob relative to the hub directory). Globs, not constructed names, # because the files under bbi/ carry the full asmId with its assembly-name suffix while the # hub directory is named with the short accession. +# Several globs per source, because the layouts differ: a GenArk assembly keeps its gene +# models under bbi/ and its CAT genes in a contrib collection, while hs1 is a hub assembly +# that keeps them in /gbdb/hs1// like a classic assembly does. genArkGeneSources = [ - ("catGenes", "contrib/*/catGenes.bb"), - ("ncbiRefSeq", "bbi/*.ncbiRefSeq.bb"), - ("ncbiGene", "bbi/*.ncbiGene.bb"), + ("catGenes", ["contrib/*/catGenes.bb", "catLiftOffGenesV1/*.bb"]), + ("ncbiRefSeq", ["bbi/*.ncbiRefSeq.bb", "ncbiRefSeq/ncbiRefSeq.bb"]), + ("ncbiGene", ["bbi/*.ncbiGene.bb"]), ] # Augustus is deliberately not in that list. It is an ab initio prediction, so mapping # UniProt through it stacks its errors on top of ours; aligning the proteins straight to # the genome with miniprot is better and much faster than the old BLAT protein search. miniprotBin = "/hive/data/outside/otto/uniprot/bin/miniprot" miniprotThreads = 16 def miniprotVersion(): " version string of the miniprot binary we are using " proc = subprocess.Popen([miniprotBin, "--version"], stdout=PIPE, encoding="utf8") return proc.communicate()[0].strip() dbIsHubCache = {} def genArkHubDir(db): - """ return the /gbdb/genark directory for a GenArk assembly, or None for a classic db. - dbDb.nibPath is "hub:/gbdb/genark/GCF/029/289/425/GCF_029289425.2" for these. + """ return the hub directory of a hub assembly, or None for a classic db. + dbDb.nibPath is "hub:" for these, e.g. + hub:/gbdb/genark/GCF/029/289/425/GCF_029289425.2 for a GenArk assembly. """ if db not in dbIsHubCache: rows = list(runQuery("hgcentral", "select nibPath from dbDb where name='%s'" % db, usePublic=True)) nibPath = rows[0][0] if rows else "" dbIsHubCache[db] = nibPath[len("hub:"):] if nibPath.startswith("hub:") else None return dbIsHubCache[db] +def isGenArk(db): + """ true only for the GenArk assemblies, the ones under /gbdb/genark with the + accession-sharded layout. Not every hub assembly is one: hs1 is served as a hub from + /gbdb/hs1/hubs but keeps its track files in /gbdb/hs1/ like a classic assembly, and + its outputs belong there rather than in a contrib collection. + """ + hubDir = genArkHubDir(db) + return hubDir is not None and "/gbdb/genark/" in hubDir + +def geneSearchRoots(db): + " the directories that may hold gene model bigBeds for a hub assembly, best first " + roots = [] + hubDir = genArkHubDir(db) + if hubDir is not None: + roots.append(hubDir) + gbdbDir = join("/gbdb", db) + if isdir(gbdbDir): + roots.append(gbdbDir) + return roots + def genArkFile(db, hubDir, fileGlob, what): " return the single file matching fileGlob inside a GenArk hub directory " matches = glob.glob(join(hubDir, fileGlob)) if len(matches)!=1: errAbort("%s: expected exactly one %s in %s, found %d: %s" % \ (db, what, hubDir, len(matches), matches)) return matches[0] def genArkTwoBit(db, hubDir): " the assembly sequence of a GenArk assembly " # named with the short accession, unlike the files under bbi/ return genArkFile(db, hubDir, "*.2bit", "2bit file") def genArkChromSizes(db, hubDir): " the chrom.sizes of a GenArk assembly. Note the .txt suffix, which /hive does not use " return genArkFile(db, hubDir, "*.chrom.sizes.txt", "chrom.sizes file") def findBestGeneBigBed(db, hubDir): """ find the best gene model bigBed for a GenArk assembly. Returns (sourceName, fileName), or (None, None) if the hub has no gene models at all. """ - for name, fileGlob in genArkGeneSources: - matches = sorted(glob.glob(join(hubDir, fileGlob))) + for name, fileGlobs in genArkGeneSources: + for root in geneSearchRoots(db): + for fileGlob in fileGlobs: + matches = sorted(glob.glob(join(root, fileGlob))) if matches: if len(matches) > 1: - logging.warning("%s: %d files match %s, using %s" % (db, len(matches), fileGlob, matches[0])) + logging.warning("%s: %d files match %s, using %s" % \ + (db, len(matches), fileGlob, matches[0])) logging.info("%s: best gene models are %s, from %s" % (db, name, matches[0])) return name, matches[0] - logging.warning("%s: no gene models found in %s, will align the proteins directly with miniprot" % (db, hubDir)) + logging.warning("%s: no gene models found in %s, will align the proteins directly " + "with miniprot" % (db, geneSearchRoots(db))) return None, None def findBestGeneTable(db): " find the best gene table for a given organism and return it " hubDir = genArkHubDir(db) if hubDir is not None: name, _ = findBestGeneBigBed(db, hubDir) return name if name else "miniprot" if db=="hg19": tables = ["refGene"] # because in 2021, refSeq maps NM_001129826.3 still to chrX_jh159150_fix, confirmed as a bug by Terence # CSAG is an important gene elif db=="hg38": #tables = [ "wgEncodeGencodeCompV37lift37" ] @@ -2321,31 +2347,31 @@ if len(fnames)==0: logging.error("NOT FOUND: %s" % fileMask) return None assert(len(fnames)==1) fname = fnames[0] return fname def makeLinks(bigBedDir, onlyDbs, taxIdDbs): " check the /gbdb symlinks " for taxId, dbs in taxIdDbs.items(): for db in dbs: if onlyDbs is not None and db not in onlyDbs: continue - if genArkHubDir(db) is not None: + if isGenArk(db): # a GenArk assembly is served from its hub, not from /gbdb//, and # /gbdb/genark belongs to the GenArk build. Its files go into the contrib # collection instead, see installGenArkContrib. logging.debug("%s is a GenArk assembly, no /gbdb symlinks" % db) continue dbBigBedDir = join(bigBedDir, db) # find the bigBed files bbTargetDir = makeSubDir(join("/gbdb", db), "uniprot") bbFnames = glob.glob(join(dbBigBedDir, "*.bb")) logging.debug("Found %d bigBed files in %s" % (len(bbFnames), dbBigBedDir)) # and create links to them for bbFname in bbFnames: @@ -2520,31 +2546,31 @@ # trackDb dataVersion setting and its date is what people check to see whether # the pipeline is alive, so rewriting an identical file would make a stalled # track look freshly updated. refs #38300 oldVersion = None if isfile(versionFname): oldVersion = open(versionFname, encoding="utf8").read() if oldVersion == fullVersion: logging.debug("%s already holds the current release string, not rewriting" % versionFname) else: with open(versionFname, "w", encoding="utf8") as versionOfh: versionOfh.write(fullVersion) logging.info("Wrote release string to %s" % versionFname) hubDir = genArkHubDir(db) - if hubDir is not None: + if isGenArk(db): # a GenArk assembly has no /gbdb//uniprot to link into; its files and # its trackDb go into the contrib collection, which carries the version # string inline because a contrib trackDb cannot read a /gbdb path installGenArkContrib(db, hubDir, dbDir, fullVersion, shortVersion) continue linkName = join("/gbdb", db, "uniprot", "version.txt") makeSymlink(versionFname, linkName) logging.debug("Release string in %s, symlink from %s" % (versionFname, linkName)) return versionString, shortVersion genArkContribRoot = "/hive/data/genomes/asmHubs/contrib" genArkContribName = "uniprot" @@ -2613,31 +2639,31 @@ tdbLines = tdb.splitlines() tdbLines = [l for l in tdbLines if not l.startswith("#")] tdbFname = join(archDir, "trackDb.txt") with open(tdbFname, "w") as ofh: ofh.write("\n".join(tdbLines)) logging.info("Created %s" % tdbFname) 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 genArkHubDir(db) is not None: + if isGenArk(db): # goldenPath/archive// 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) inDir = join(bigBedDir, db) count = 0 for inFname in glob.glob(join(inDir, "*")): shutil.copyfile(inFname, join(archDir, basename(inFname)))