9541aea4204c8079ea2e464403f2225f6aa33cfd max Mon Sep 14 07:59:00 2026 -0700 uniprot: show the splice variant track, and filter the CAT alignments Two things that were quietly missing. unipSplice has been built since the pipeline rewrite in #19351 and never had a trackDb stanza - "git log -S unipSplice" on uniprot.ra returns nothing, so it was never wired up rather than deliberately dropped. It holds UniProt's splice variant features, 28962 of them on hg38 and 29195 on hs1, and has been invisible on every assembly for years. Added to uniprot.ra and to the archive/contrib template, which is where the hub-served assemblies get their trackDb. The file exists on the same 131 assemblies as unipDomain.bb, so no stanza points at anything missing. The alignments on a CAT assembly were not being filtered with pslSelect. That filter maps a UniProt accession to the transcripts UniProt cross-references for it, and the README is blunt about how much it matters for protein families with nearly identical transcripts. It understood two kinds of id, Ensembl and RefSeq, and CAT names a transcript after its source gene, so hs1 matched neither and fell through unfiltered. Every row of the CAT bigBed carries the Ensembl transcript it was lifted from, so catSourceTransMap joins on that and UniProt's Ensembl cross-reference does the rest. The mapping had to become one-to-many for this: one Ensembl transcript can name several of ours, both because paralogs are lifted from the same source and because duplicate CAT names were given -dup suffixes earlier. It reproduces those suffixes by walking the bigBed in the same order rather than attaching every paralog to every source, which measured 234903 of 234903 hs1 transcripts mapped, no duplicates, against 330749 for the loose version. For every other gene track the lists hold one element and the behaviour is unchanged, which is checked: a single-valued entry still writes exactly one pair line and still counts a version difference, a multi-valued one writes a line per transcript, and an id we do not have is still skipped so pslSelect -qPass passes it through. refs #38300 diff --git src/hg/utils/otto/uniprot/doUniprot src/hg/utils/otto/uniprot/doUniprot index 34cfc56281f..c2858fca1a0 100755 --- src/hg/utils/otto/uniprot/doUniprot +++ src/hg/utils/otto/uniprot/doUniprot @@ -1627,85 +1627,148 @@ logging.debug("Reading transcript IDs from %s" % transcriptFa) allTrans = set(parseFasta(open(transcriptFa)).keys()) return allTrans logging.info("Getting transcript IDs for db=%s, table=%s, field=name" % (db, geneTable)) sql = "SELECT DISTINCT name from "+geneTable cmd = ["hgsql", db, "-NBe", sql] proc = subprocess.Popen(cmd, encoding="latin1", stdout=PIPE) allTrans = set() for line in proc.stdout: idStr = line.rstrip() allTrans.add(idStr) logging.info("Got %d distinct transcript IDs from MySQL table" % len(allTrans)) return allTrans +def catSourceTransMap(db, allTransIds): + """ map Ensembl transcript accession, without version, to the CAT transcript names we + have in our fasta. + + CAT names a transcript after the source gene, e.g. AL627309.3-1, which matches nothing + UniProt knows. Every row of the CAT bigBed carries a sourceTranscript field holding the + Ensembl transcript it was lifted from, so that is the join. One Ensembl transcript can + name several of ours, because paralogs are lifted from the same source, and because + duplicate CAT names were made unique with a -dup suffix earlier. + """ + hubDir = genArkHubDir(db) + geneBb = findBestGeneBigBed(db, hubDir)[1] if hubDir else None + if geneBb is None: + logging.warning("%s: no CAT bigBed found, cannot map to Ensembl transcripts" % db) + return {} + + # which column is sourceTranscript + asText = subprocess.Popen(["bigBedInfo", "-as", geneBb], stdout=PIPE, encoding="utf8").communicate()[0] + fields = re.findall(r"^\s+\S+\s+\**(\w+)\s*;", asText, flags=re.MULTILINE) + if "sourceTranscript" not in fields or "name" not in fields: + logging.warning("%s: %s has no sourceTranscript field, cannot map to Ensembl" % (db, geneBb)) + return {} + nameIdx, srcIdx = fields.index("name"), fields.index("sourceTranscript") + + # Reproduce the -dup suffixes exactly rather than attaching every paralog to every + # source transcript. makeTranscriptFiles walks the bigBed in the same order and renames + # the second and later rows that share a name, so the same running count recovers which + # of our transcripts each bigBed row became. + shortToLong = {} + seen = {} + proc = subprocess.Popen(["bigBedToBed", geneBb, "stdout"], stdout=PIPE, encoding="utf8") + for line in proc.stdout: + row = line.rstrip("\n").split("\t") + if len(row) <= max(nameIdx, srcIdx): + continue + name = row[nameIdx] + seen[name] = seen.get(name, 0) + 1 + ourId = name if seen[name]==1 else "%s-dup%d" % (name, seen[name]) + if ourId not in allTransIds: + continue + src = row[srcIdx].split(".")[0] + if not src or src in ("None", "NA"): + continue + shortToLong.setdefault(src, []).append(ourId) + proc.stdout.close() + + logging.info("%s: %d Ensembl transcripts map to %d of our %d CAT transcripts" % \ + (db, len(shortToLong), sum(len(v) for v in shortToLong.values()), len(allTransIds))) + return shortToLong + def writeTransToSelect(upSeqIds, transIds, shortToLong, notFound, outTransIds, commonTrans, \ ofh, verDiffCount, pairCount): """ given a list of transcript IDs, write them to select file, if their base (before dot) is found in shortToLong. Return count written """ foundTransCount = 0 for transId in transIds: # remove the version for the comparison, but keep the version in the pslSelect file shortId = transId.split(".")[0] if shortId not in shortToLong: # if we do not have this transcript, we don't write anything into the pair select file # This means that pslSelect falls back to the defaults, passes through everything and then # pslCdnaFilter picks the best ones notFound.append(transId) continue commonTrans.add(transId) - if transId!=shortToLong[shortId]: + # ourIds is a list: usually one transcript, but a CAT gene set can name several, + # since paralogs are lifted from the same source transcript + ourIds = shortToLong[shortId] + if transId not in ourIds: verDiffCount+=1 pairCount += 1 outTransIds.add(transId) # here we're fixing up the select file: we replace the acc.version that UniProt has # with the acc.version that we have in our transcript file. This saves a few hundred # transcripts, sometimes thousands, from getting removed for upSeqId in upSeqIds: - ofh.write("%s\t%s\n" % (upSeqId, shortToLong[shortId])) + for ourId in ourIds: + ofh.write("%s\t%s\n" % (upSeqId, ourId)) 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: + if "encode" in geneTable or "ensGene" in geneTable or geneTable=="catGenes": + # CAT is lifted from Ensembl/GENCODE, so its transcripts carry an Ensembl source + # transcript and UniProt's Ensembl cross-reference is the one that matches transField = "ensemblTrans" elif "refGene" in geneTable or "ncbiRefSeq" in geneTable: transField = "refSeq" else: 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 # 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 + if geneTable=="catGenes": + # CAT names its transcripts after the source gene, so the ids in our fasta match no + # UniProt cross-reference at all. The CAT bigBed does carry, per transcript, the + # Ensembl transcript it was lifted from, so map through that and use UniProt's + # Ensembl cross-reference. refs #38300 + shortToLong = catSourceTransMap(db, allTransIds) + else: shortToLong = {} for acc in allTransIds: - shortToLong[acc.split(".")[0]] = acc + shortToLong.setdefault(acc.split(".")[0], []).append(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()) logging.info("Building pair file %s" % selectFname) ofh = open(selectFname, "w") notFound = [] # the following loop is FULL of stat tracking. It took me a long time to get my head around # all the possible ways that the mapping can go astray and I wanted to keep them all # This makes the code harder to read but simplifies debugging a lot when ML questions come in.