947dcb525b76d579bd6ab425cf71a221d9bab347 max Fri Sep 11 12:50:09 2026 -0700 uniprot otto: give the unfiltered gene-track case a map source hs1 cleared pslMap after the duplicate-name fix, ran the full 91 minute alignment, and then died writing the bigPsl with KeyError: 'default' at mapSource = accToMapSource["default"] buildSelectFile returns an empty dict when no UniProt cross-reference matches the ids in the gene track, and pslToBigPsl looks up accToMapSource["default"] for every protein it writes. That branch was unreachable before: every gene table findBestGeneTable could return was handled somewhere else, with augustus and the direct protein alignment each setting their own default. catGenes and ncbiGene are the first to reach it. It now returns {"default":"best"}, the same answer the augustus case gives, which is the honest description: the alignment is the best match we found, we just cannot name the transcript evidence behind it. Hardened the lookup as well. protMapSource on the supported path only ever gets per-accession entries and never a "default" key, so any accession missing from it would have raised the same KeyError on hg38. It now falls back to "best" rather than throwing away hours of cluster work over one protein. Also reworded the log line, which said "No supported gene track found" while a perfectly good gene track was in use; the missing thing is a cross-reference that matches its ids. refs #38300 diff --git src/hg/utils/otto/uniprot/doUniprot src/hg/utils/otto/uniprot/doUniprot index 0ce0a7e180b..c017c44c890 100755 --- src/hg/utils/otto/uniprot/doUniprot +++ src/hg/utils/otto/uniprot/doUniprot @@ -1370,31 +1370,35 @@ acc = row[3] recAnnot = accToMeta.get(acc, None) if recAnnot is None: logging.error("Accession: %s - no record info?" % acc) assert(False) # add the transcripts that were used for mapping this protein to the genome chrom, start, end = row[0], row[1], row[2] if mapInfo: transListStr = ", ".join(sorted(list(mapInfo[(acc, chrom, start, end)]))) mapSource = accToMapSource.get(acc) if mapSource is None: - mapSource = accToMapSource["default"] + # "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" 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) @@ -1674,34 +1678,38 @@ for upSeqId in upSeqIds: ofh.write("%s\t%s\n" % (upSeqId, shortToLong[shortId])) foundTransCount+=1 return foundTransCount, verDiffCount, pairCount def buildSelectFile(tabFnames, mapDir, taxId, db, geneTable, transcriptFa, geneToRefSeqs, uniprotMd5, transMd5, stats): """ make a two-column file uniprotId <-> transcript ID for filtering the PSL file. The tool for this is called pslSelect. This helps force the alignment to the correct transcript. It turned out to be absolutely necessary. see makeUniProtPsl.sh """ # get the UniProt field that holds the list of annotated transcript IDs if "encode" in geneTable or "ensGene" in geneTable: transField = "ensemblTrans" elif "refGene" in geneTable or "ncbiRefSeq" in geneTable: transField = "refSeq" else: - logging.info("No supported gene track found. The protein/transcript alignments will not be filtered using pslSelect.") + logging.info("%s: no UniProt cross-reference matches the ids in %s, so the " + "protein/transcript alignments are not filtered with pslSelect" % (db, geneTable)) stats["noGeneTrack"] = True - # anoGam1 for example, has only augustus genes as a track - return None, stats, {} + # Every alignment that survives is then just the best match we found, so say so: + # the same map source the augustus case uses. Returning an empty dict here used to + # be safe only because nothing reached this branch; catGenes and ncbiGene do, and + # pslToBigPsl looks up accToMapSource["default"] for every protein it writes. + return None, stats, {"default":"best"} allTransIds = getTransIds(db, geneTable, transcriptFa) # create a mapping from no-version (short) to with-version accession (long), so we can fix them up later # to the accessions that we actually have in our fasta file shortToLong = {} for acc in allTransIds: shortToLong[acc.split(".")[0]] = acc logging.debug("Example transcript IDs: %s" % list(allTransIds)[:10]) uniprotMd5 = uniprotMd5[:10] transMd5 = transMd5[:10] selectFname = join(mapDir, "%(geneTable)s_%(uniprotMd5)s_%(transMd5)s.upToTrans.tab" % locals())