08221b327a1496fd1e4d3f88dc17e839b15b64f7 max Tue Sep 15 05:24:18 2026 -0700 uniprot otto: minAli must not be global, and must be part of the cache key max spotted that ci3's lower alignment threshold was assigned to the global MINALI. With --taxonThreads that leaks: ci3 sets 0.85 and every assembly whose thread reads the global afterwards uses it too. In the 2026_03 run it did exactly that - 89 of 95 alignments ran at 0.85 instead of 0.93, and which assemblies escaped depended purely on thread timing, with ce11 and cb3 getting 0.93 by luck. His fix makes it a local; this commit keeps that and closes the second half of the problem. The mapping files are cached under an md5 of the protein sequences, the transcripts, their alignment and the protein/transcript pair file. The threshold was not in it, so a mapping built at 0.85 would be silently reused on a run that asked for 0.93, and the bad alignments would have survived the rebuild that is meant to fix them. minAli is now part of the key, so changing it invalidates the mapping by itself. Checked that the key differs between 0.93 and 0.85. Also gave lastRun_noAnnotatedTranscriptIDs.txt a per-assembly name. It is written into the shared mapDir under a fixed name by every assembly, which several threads would write at once. Only a debugging artifact, nothing reads it, but it is the last shared write in the per-taxon path. Audited the rest of the module for state shared across threads: the only remaining global is flagFname, set once in main before the pool starts; dbIsHubCache can be written by two threads but only ever with the same value; and every temp directory and working file is keyed by assembly, md5 or taxon. refs #38300 diff --git src/hg/utils/otto/uniprot/doUniprot src/hg/utils/otto/uniprot/doUniprot index 96f1d9d05ad..7391929e714 100755 --- src/hg/utils/otto/uniprot/doUniprot +++ src/hg/utils/otto/uniprot/doUniprot @@ -1866,31 +1866,33 @@ protMapSource[seqId] = mapSource sourceCounts[mapSource] += 1 if mapSource=="best": notFoundAtAll.add(upAcc) sourceCounts = list(dict(sourceCounts).items()) stats["notFoundAtAll"] = len(notFoundAtAll) stats["notFoundAtAll_Examples"] = list(notFoundAtAll)[:10] stats["verDiffCount"] = verDiffCount stats["verDiffCount2"] = verDiffCount2 stats["sourceCounts"] = sourceCounts # for debugging if len(notFoundAtAll)!=0: - notFoundFname = join(mapDir, "lastRun_noAnnotatedTranscriptIDs.txt") + # per assembly: mapDir is shared, and with --taxonThreads several assemblies would + # otherwise write this same file at the same time + notFoundFname = join(mapDir, "lastRun_noAnnotatedTranscriptIDs.%s.txt" % db) notFoundFh = open(notFoundFname, "w") #notFoundFh.write(str(notFoundAtAll)+"\n") for acc in notFoundAtAll: notFoundFh.write(acc+"\n") notFoundFh.close() logging.info("Wrote UniProt IDs where no transcript could be found to %s" % notFoundFname) if pairCount==0: logging.warn("No single uniprot/transcript pair is left in select file. Not doing pair file filtering at all.") stats["outPairCount"]=0 return None, stats, protMapSource # make a little table with NM_ -> count, XM_ -> count, etc prefixCounts = defaultdict(int) for nfId in notFound: @@ -2065,31 +2067,30 @@ key, val = line.rstrip("\n").split("\t") vals[key].append(val) return vals def makeProteinGenomePsl(taxId, protFa, tabFnames, db, doForce, mapDir, mapDescFname): """ map proteins of a taxon stored in protFa to db via geneTable and create a PSL file for the mapping Will try to reuse the old PSL mapping file if possible, as it is expensive to build. Need to BLAST the entire file: BLAST has no option to subset the target (unlike blast or lastz), verified with an email to NCBI BLAST support team. A central point of this function is buildSelectFile(). The psl select pair files have a huge influence on the result, so it is rebuilt every time, as any change in UniProt will change it. It's also quick to create. Its MD5 is part of the PSL mapping file. Any change to the pslSelect pair file will trigger a realignment/rebuild of the pslMap PSL file. """ # use MD5s of the input files to determine if the protein -> genome map has to be rebuilt - global MINALI protMd5 = fastaMd5(protFa) 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) @@ -2102,66 +2103,72 @@ 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) protToTrans = parseKeyVal(selectFname) - allMd5s = [protMd5, transMd5, pslMd5] + # ciona's genome is from a different organism than the refseq transcripts, so it needs a + # lower threshold. Keep it a local and never assign to the global MINALI: with + # --taxonThreads several assemblies run at once and would see each other's value. + minAli = 0.85 if db in ["ci3"] else MINALI + + allMd5s = [protMd5, transMd5, pslMd5, str(minAli)] if selectFname: pairMd5 = tabMd5(selectFname) allMd5s.append(pairMd5) stats["pairMd5"] = pairMd5[:10] stats["pairFname"] = selectFname - # if either the protein sequences, the transcripts, their alignment or the - # protein/transcript mapping changes -> create a new PSL mapping file + # if either the protein sequences, the transcripts, their alignment, the + # protein/transcript mapping or the alignment threshold changes -> create a new PSL + # mapping file. minAli belongs in here: it decides which alignments survive, so a + # mapping built at one threshold must not be silently reused at another. That is + # exactly what happened after the MINALI global leaked across threads. refs #38300 fullMd5 = listMd5(allMd5s)[:10] mapFname = join(mapDir, "%(geneTable)s_%(fullMd5)s.psl" % locals()) 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 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 + stats["minAli"] = minAli workDir = "clusterRun-map-%s-%s-%s.tmp" % (db, geneTable, fullMd5) - scriptArgs = [protFa, transcriptFa, transcriptPsl, str(MINALI), workDir, mapFname] + 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) run(cmd) writeMapDesc(stats, db, geneTable, mapDescFname) if os.path.getsize(mapFname)==0: # Nothing aligned. On an organism with a handful of UniProt proteins this is a fact # about the data, not a failure of the run: taxon 216574, the golden eagle, offers # two proteins and neither of them aligns. Treat it the way an organism with no