5783988b48ed0722ce2398b50bc1a422911c1fa7
max
  Fri Sep 11 06:02:22 2026 -0700
uniprot otto: teach the gene-model search about GenArk assemblies

findBestGeneTable only knew how to look for MySQL tables, so on a GenArk assembly
it found nothing and fell through to BLAT, which then died on a 2bit path that
does not exist for a hub assembly. GenArk keeps its gene models as bigBed files
in the hub instead.

Search order, best first:

catGenes    Comparative Annotation Toolkit, shipped as a contrib collection
(track hprcCatGenes, contrib/hprc2annot/catGenes.bb)
ncbiRefSeq  the GenArk RefSeq gene track, bbi/*.ncbiRefSeq.bb
ncbiGene    the annotation the submitter sent to GenBank with the assembly
augustus    ab initio, present nearly everywhere, so it is the last resort

Matched by glob rather than by constructed name: files under bbi/ carry the full
asmId including its assembly-name suffix, while the hub directory is named with
the short accession. A dbDb row is recognised as GenArk by nibPath starting with
"hub:", cached so the lookup happens once per db.

Checked against the assemblies that are actually blocking the update: the twelve
GCF ones resolve to ncbiRefSeq, GRCz12ab and calJac240_pri have no RefSeq or
submitted annotation and fall to augustus, and on an HPRC assembly that carries
both, catGenes wins over augustus as intended. Classic assemblies are untouched -
genArkHubDir returns None for them and hg38/hg19/mm39/panPan3/rn6 still resolve
exactly as before through the MySQL path.

This is the search only. Building from these still needs the transcript fasta and
PSL to come from the hub rather than from MySQL. refs #38300

diff --git src/hg/utils/otto/uniprot/doUniprot src/hg/utils/otto/uniprot/doUniprot
index bcbdb3cc1f0..81963a28f3e 100755
--- src/hg/utils/otto/uniprot/doUniprot
+++ src/hg/utils/otto/uniprot/doUniprot
@@ -984,32 +984,81 @@
     if table=="knownGene":
         cmd = "hgsql %s -N -e 'select * from %s' | sort -k1 | cut -f1 | uniq -c | sort -rn | head -n1" % \
                 (db, table)
     else:
         cmd = "hgsql %s -N -e 'select * from %s' | cut -f2- | sort -k1 | cut -f1 | uniq -c | sort -rn | head -n1" % \
                 (db, table)
     proc = subprocess.Popen(cmd, stdout=PIPE, shell=True, encoding="latin1")
     row1 = proc.stdout.readline().strip().split()
     count, gene = row1
     count = int(count)
     if count!=1:
         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.
+genArkGeneSources = [
+    ("catGenes",   "contrib/*/catGenes.bb"),
+    ("ncbiRefSeq", "bbi/*.ncbiRefSeq.bb"),
+    ("ncbiGene",   "bbi/*.ncbiGene.bb"),
+    ("augustus",   "bbi/*.augustus.bb"),
+]
+
+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.
+    """
+    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 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)))
+        if matches:
+            if len(matches) > 1:
+                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 have to BLAT the proteins" % (db, hubDir))
+    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 "blat"
+
     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" ]
         tables = ["ncbiRefSeq"]
     else:
         tables = [ "ncbiRefSeq", "ensGene", "augustusGene"]
         # refGene has only a few hundred transcripts on most obscure organisms -> refGene only covers the manual NM_ sequences!
 
     for table in tables:
         query = "DESCRIBE %s" % (table)
         cmd = "hgsql %s -e 'DESCRIBE %s' > /dev/null 2>&1" % (db, table)
         ret = os.system(cmd)