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())