99a26061f6cf8d3ee709f3468dda88af74e4049d
max
  Mon Sep 14 06:17:21 2026 -0700
uniprot otto: convert miniprot alignments to PSL properly

GRCz12ab cleared miniprot and then died in the annotation lift:

Error: inPsl Q98TT6 tSize (336) != mapPsl Q98TT6 qSize (339)

The lift maps annotations that are given in protein coordinates, so the mapping PSL
has to have the protein as its query, with qSize three times the protein length.
Routing miniprot's output through gff3ToGenePred and genePredToFakePsl does not give
that: it makes the query the transcript implied by the alignment, so qSize comes out
as the aligned CDS length. Measured over 93518 zebrafish alignments, that was wrong
for 95% of them - 79% out by exactly one codon, the trailing stop, and 16% out by
other amounts where miniprot aligned only part of the protein.

pafToPsl cannot do it either: it rejects miniprot's CIGAR, which is splice aware and
uses operators it does not know.

But every CDS line of a miniprot GFF carries its own "Target=<acc> <start> <end>"
giving the protein range that block covers, so the alignment can be reconstructed
exactly. New miniprotToPsl does that.

Two things it has to get right. On the minus strand miniprot lists blocks in protein
order, which is descending genomic order, while a PSL lists them ascending and puts
the block starts on the reverse complemented query, with qStart and qEnd still in
forward coordinates. And a block whose genomic span is shorter than its protein range
implies, where miniprot placed a frameshift, is clamped to the genome, with the
overall ranges then derived from the blocks rather than from the protein ranges.

Verified against real bonobo alignments: pslCheck reports 65 checked, 0 failed; qSize
is three times the protein length for every row; qName, strand and tStart match
miniprot exactly for all 65; and tEnd matches for 43, is short by exactly 3 on 12,
and those 12 are the ones carrying a stop_codon feature, which is correct since the
protein has no stop codon. 26 of the rows are on the minus strand.

refs #38300

diff --git src/hg/utils/otto/uniprot/doUniprot src/hg/utils/otto/uniprot/doUniprot
index 2379ed0671c..8676593dab5 100755
--- src/hg/utils/otto/uniprot/doUniprot
+++ src/hg/utils/otto/uniprot/doUniprot
@@ -1588,63 +1588,39 @@
     stacks its errors on top of ours, and the BLAT protein search is very slow on a big
     genome.
     """
     workDir = mapFname+".miniprot.tmp"
     if not isdir(workDir):
         os.makedirs(workDir, exist_ok=True)
 
     # miniprot reads fasta, not 2bit
     genomeFa = join(workDir, "genome.fa")
     run(["twoBitToFa", twoBitFname(db), genomeFa])
 
     gffName = join(workDir, "miniprot.gff")
     runMiniprotOnCluster(db, genomeFa, fullFaFname, gffName, workDir)
     os.remove(genomeFa)
 
-    # Name each alignment after the UniProt accession it came from. miniprot calls them
-    # MP000001 and puts the accession in Target=, and its ##PAF meta lines and capitalised
-    # attributes make gff3ToGenePred unhappy, so both are cleaned up here. The ID has to
-    # stay unique for GFF3, hence the .N suffix, which comes off again after the
-    # conversion so that qName is the bare accession.
-    namedName = join(workDir, "named.gff")
-    cmd = """awk -F'\\t' -v OFS='\\t' '
-      /^##PAF/ {next}
-      /^#/     {print; next}
-      $3=="mRNA" {
-          id=""; acc=""
-          if (match($9, /ID=[^;]+/))      id  = substr($9, RSTART+3, RLENGTH-3)
-          if (match($9, /Target=[^ ;]+/)) acc = substr($9, RSTART+7, RLENGTH-7)
-          if (id!="" && acc!="") { map[id]=acc; $9 = "ID=" acc "." (++n[acc]) }
-          print; next
-      }
-      { if (match($9, /Parent=[^;]+/)) {
-            p = substr($9, RSTART+7, RLENGTH-7)
-            if (p in map) $9 = "Parent=" map[p] "." n[map[p]]
-        }
-        print
-      }' %s > %s""" % (gffName, namedName)
-    run(cmd)
-
-    gpName = join(workDir, "miniprot.gp")
-    run(["gff3ToGenePred", "-warnAndContinue", "-maxParseErrors=-1", "-maxConvertErrors=-1",
-         namedName, gpName])
-    run("""awk -F'\\t' -v OFS='\\t' '{sub(/\\.[0-9]+$/, "", $1); print}' %s > %s.acc""" % (gpName, gpName))
-
-    cdsName = join(workDir, "miniprot.cds")
+    # miniprot's GFF carries, on every CDS block, the protein range that block covers, so
+    # the alignment can be turned into a PSL exactly. Going via gff3ToGenePred and
+    # genePredToFakePsl instead makes the query the transcript implied by the alignment,
+    # which gives a qSize of the aligned CDS length rather than of the protein, and the
+    # annotation lift needs the protein. refs #38300
+    myDir = dirname(__file__)
     chromSizes = join(workDir, "chrom.sizes")
     run("twoBitInfo %s %s" % (twoBitFname(db), chromSizes))
-    run(["genePredToFakePsl", "-chromSize=%s" % chromSizes, "noDb", gpName+".acc", mapFname, cdsName])
+    run([join(myDir, "miniprotToPsl"), gffName, fullFaFname, chromSizes, mapFname])
 
     alnCount = len(open(mapFname).readlines())
     logging.info("%s: miniprot aligned the UniProt proteins to %d genomic locations" % (db, alnCount))
     stats["aligner"] = "miniprot " + miniprotVersion()
 
 def getTransIds(db, geneTable, transcriptFa):
     """ return all possible transcript IDs given a gene table. The reason that we need is that Uniprot often
       refers to transcripts that don't exist and we are using a select file to map only to those. In these
       cases, we want to make sure that our select file contains only transcript that we actually have,
       so we can log how many transcript can possibly mapped (and which ones). This is important for debugging.
     """
     # no idea how to get the version suffix out of our mysql tables - thanks Obama!
     # A GenArk assembly has no MySQL tables, but the fasta we built from its gene bigBed
     # holds every transcript ID already, so read them from there as refGene does.
     if geneTable=="refGene" or genArkHubDir(db) is not None: