b8320ab5d2dd4d520742faedc2def3144595c450 max Mon Sep 14 05:49:34 2026 -0700 uniprot otto: writeMapDesc crashed on a hub assembly with no gene models GRCz12ab got through miniprot on the cluster and then died writing its lift info: os.path.getmtime(geneBb) -> stat on None The hub branch dates the gene models by the bigBed that was read, but an assembly that has no gene models is precisely the one that reaches miniprot, and there is no bigBed to date. The miniprot case was already handled further down the chain, so the fix is only to test it first. Checked every combination that can occur: GRCz12ab and calJac240_pri, hub assemblies with no models, report "direct"; hs1 dates its CAT models 2022-03-15; mPanPan1_v2.0 dates its GenArk RefSeq 2025-06-02; hg38 still reads its real ncbiRefSeqVersion.txt and rn6 still queries trackVersion. refs #38300 diff --git src/hg/utils/otto/uniprot/doUniprot src/hg/utils/otto/uniprot/doUniprot index 97b2abd17d8..38a0f6d6631 100755 --- src/hg/utils/otto/uniprot/doUniprot +++ src/hg/utils/otto/uniprot/doUniprot @@ -1381,41 +1381,41 @@ if mapInfo: transListStr = ", ".join(sorted(list(mapInfo[(acc, chrom, start, end)]))) mapSource = accToMapSource.get(acc) if mapSource is None: # "best" is the honest answer when we have nothing better to say: the # alignment survived, but we cannot name the transcript evidence behind it. # Do not turn this into a KeyError - it would throw away hours of cluster # work over one protein. mapSource = accToMapSource.get("default", "best") protToLocs[acc].append( ("%s:%s-%s" % (chrom, str(start), str(end)), mapSource, transListStr.replace(" ", "")) ) if mapSource == "best": mapDesc = "Best match(es) when aligning protein against all transcripts" elif mapSource == "direct": - mapDesc = "Best BLAT match when aligning protein against entire genome" + mapDesc = "Best miniprot match when aligning protein against entire genome" elif mapSource == "uniprot": mapDesc = "Alignment of protein to transcript(s) annotated by UniProt" elif mapSource == "entrez": mapDesc = "Best match(es) when aligning protein against RefSeq transcripts of NCBI gene annotated by UniProt" else: assert(False) # bug mapInfoStr = "%s: %s" % (mapDesc, transListStr) else: - mapInfoStr = "no transcript: direct BLAT to genome" + mapInfoStr = "no transcript: direct miniprot alignment to genome" row.append(mapInfoStr) protDb = accToDb[acc] # there are two types of sequences in UniProt: the "main" ones, that # are annotated, and the others ("isoforms"), which are not annotated. I now color # the alternative isoforms like trembl sequences isMain = (acc==recAnnot.mainIsoAcc) isMainStr = "alternative isoform" if isMain: isMainStr = "primary sequence" if protDb=="swissprot": if isMain: color = "2,12,120" @@ -1953,42 +1953,45 @@ gpName = join(dbDir, "transcripts.gp") cdsName = join(dbDir, "transcripts.cds") runQueryToFile(db, query, gpName) cmd = ["genePredToFakePsl", db, gpName, pslName, cdsName] run(cmd) cmd = ["getRnaPred", db, gpName, "all", faName] run(cmd) logging.info("Created %s and %s for %s/%s" % (faName, pslName, db, geneTable)) return faName, pslName def writeMapDesc(stats, db, geneTable, mapDescFname): " write a little json file that points to the map file used and some basic stats about it " logging.debug("making json file with version and other lift info") hubDir = genArkHubDir(db) - if hubDir is not None: - # A GenArk assembly has no ncbiRefSeqVersion.txt. Date the gene models by the + if geneTable=="miniprot": + # no gene models were involved at all, the proteins went straight to the genome. + # This has to be tested before the hub case below: a hub assembly with no gene + # models is exactly the one that ends up here, and it has no bigBed to date. + version = "direct" + elif hubDir is not None: + # A hub assembly has no ncbiRefSeqVersion.txt. Date the gene models by the # bigBed we actually read, which is what the user needs to tie the mapping to a # particular annotation; the source name is carried separately in geneTable. geneBb = findBestGeneBigBed(db, hubDir)[1] version = datetime.datetime.fromtimestamp(os.path.getmtime(geneBb)).strftime('%Y-%m-%d') elif geneTable == "ncbiRefSeq": version = open("/gbdb/"+db+"/ncbiRefSeq/ncbiRefSeqVersion.txt").read().strip() elif geneTable in ['refGene', 'augustusGene', 'knownGene']: version = datetime.datetime.today().strftime('%Y-%m-%d') - elif geneTable=="miniprot": - version = "direct" elif geneTable in ["ensGene"]: version = list(runQuery("hgFixed", "select version from trackVersion where db='%s' and name='ensGene' order by ix desc limit 1;" % db))[0][0] else: assert(False) stats["geneTable"] = geneTable stats["version"] = version stats["createdDate"] = datetime.datetime.today().strftime('%Y-%m-%d') with open(mapDescFname, "w") as mapDescFh: json.dump(stats, mapDescFh, indent=4) def listMd5(arr): " return hex md5 given list of strings " import hashlib @@ -2080,31 +2083,31 @@ if isfile(mapFname) and not doForce: logging.info("%s already exists, not rebuilding the protein -> genome mapping PSL" % mapFname) assert(os.path.getsize(mapFname)!=0) # careful: if you modify this, also modify the other return statement below return mapFname, geneTable, fullMd5, protMapSource, protToTrans, False stats["protMd5"] = protMd5[:10] stats["metaMd5"] = metaMd5[:10] stats["transMd5"] = transMd5[:10] stats["fullMd5"] = fullMd5[:10] stats["mapFname"] = basename(mapFname) logging.debug("%s does not exist" % mapFname) if geneTable=="miniprot": - logging.error("Could not find any gene table for %s, using BLAT to map proteins" % db) + logging.error("Could not find any gene table for %s, using miniprot to map proteins" % db) miniprotProteins(protFa, db, mapFname, stats) else: if db in ["ci3"]: # wow: ciona's genome is from a different organism than the refseq transcripts. MINALI = 0.85 stats["minAli"] = MINALI workDir = "clusterRun-map-%s-%s-%s.tmp" % (db, geneTable, fullMd5) scriptArgs = [protFa, transcriptFa, transcriptPsl, str(MINALI), workDir, mapFname] if selectFname is not None: scriptArgs.append(selectFname) # This is where the BLAST alignment-meat happens cmd = ["time", "./makeUniProtPsl.sh"] cmd.extend(scriptArgs) logging.info("Running %s" % cmd)