155514c1985d0fdc25f4a33f4e1786f35b47b09d
max
  Fri Sep 11 11:11:49 2026 -0700
uniprot otto: make CAT transcript names unique before deriving the fasta and PSL

The hs1 run got through a 52 minute BLAST batch with no crashes and then stopped in
pslMap with

Error: inPsl RBMY1F-1 tSize (1887) != mapPsl RBMY1F-1 qSize (1718)

CAT names its transcripts after the source gene, so paralogs share a name. The hs1
CAT set has 777 duplicated names among 234903 transcripts, and RBMY1F-1 is three Y
chromosome paralogs of 1810, 1718 and 1887 bases. The protein-to-transcript
alignment and the transcript-to-genome alignment then disagree about which
transcript the name refers to and pslMap refuses to map. This is a property of CAT,
not of hub assemblies: the GenArk ncbiRefSeq sets use accessions and have no
duplicates at all.

Uniquify the names on the genePred before the fasta and the PSL are derived from
it, so the two can never disagree.

Verified against the real hs1 CAT data: 234903 rows and 234903 distinct names, no
duplicates in the genePred, fasta or PSL, the three RBMY1F-1 paralogs now separate,
and every one of the 234903 PSL rows has a qSize equal to the length of its own
fasta sequence, which is the invariant pslMap enforces. On bonobo ncbiRefSeq it is
a no-op: 95418 rows, no name changed.

refs #38300

diff --git src/hg/utils/otto/uniprot/doUniprot src/hg/utils/otto/uniprot/doUniprot
index 2735e6d5f7b..91de301bdcf 100755
--- src/hg/utils/otto/uniprot/doUniprot
+++ src/hg/utils/otto/uniprot/doUniprot
@@ -1879,31 +1879,46 @@
     pslName = join(dbDir, "transcripts.psl")
 
     pslFields = "matches,misMatches,repMatches,nCount,qNumInsert,qBaseInsert,tNumInsert," \
         "tBaseInsert,strand,qName,qSize,qStart,qEnd,tName,tSize,tStart,tEnd,blockCount,blockSizes,qStarts,tStarts"
 
     hubDir = genArkHubDir(db)
     if hubDir is not None:
         # GenArk assembly: no MySQL tables at all, the gene models are a bigGenePred in the
         # hub. bigGenePredToGenePred gets us back to a genePred, and both of the tools that
         # follow can work off the hub's own chrom.sizes and 2bit rather than chromInfo, so
         # nothing here needs a database.
         geneBb = findBestGeneBigBed(db, hubDir)[1]
         gpName = join(dbDir, "transcripts.gp")
         cdsName = join(dbDir, "transcripts.cds")
 
-        run(["bigGenePredToGenePred", geneBb, gpName])
+        run(["bigGenePredToGenePred", geneBb, gpName+".raw"])
+
+        # CAT names its transcripts after the source gene, so paralogs share a name: the
+        # hs1 CAT set has 777 duplicated names in 234903 transcripts, RBMY1F-1 three times
+        # with three different lengths. pslMap then cannot tell which transcript a protein
+        # alignment meant and stops with "inPsl RBMY1F-1 tSize (1887) != mapPsl RBMY1F-1
+        # qSize (1718)". Make the names unique before the fasta and the PSL are derived
+        # from this file, so the two always agree. RefSeq accessions are unique already,
+        # so this is a no-op for them.
+        run("""awk -F'\t' -v OFS='\t' '{ seen[$1]++; if (seen[$1]>1) $1 = $1 "-dup" seen[$1]; print }' %s.raw > %s""" \
+                % (gpName, gpName))
+        dupCount = len(open(gpName+".raw").readlines()) - len(set(l.split("\t")[0] for l in open(gpName)))
+        if dupCount:
+            logging.info("%s: made %d duplicated transcript names unique" % (db, dupCount))
+        os.remove(gpName+".raw")
+
         run(["genePredToFakePsl", "-chromSize=%s" % chromSizesFile(db),
              "noDb", gpName, pslName, cdsName])
         run(["getRnaPred", "-genomeSeqs=%s" % twoBitFname(db),
              "noDb", gpName, "all", faName])
 
         logging.info("Created %s and %s for %s/%s from %s" % (faName, pslName, db, geneTable, geneBb))
         return faName, pslName
 
     if geneTable=="ncbiRefSeq":
         origFa = "/gbdb/%s/ncbiRefSeq/seqNcbiRefSeq.rna.fa" % db
         shutil.copyfile(origFa, faName)
         pslTable = "ncbiRefSeqPsl"
 
         query = "SELECT "+pslFields+" from "+pslTable+" WHERE tName NOT LIKE '%_hap%' AND tName not like '%_alt%' AND tNAME NOT LIKE '%_fix%'"
         runQueryToFile(db, query, pslName)