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/makefile src/hg/utils/otto/uniprot/makefile
index 48c58aefac1..0c53981cb08 100644
--- src/hg/utils/otto/uniprot/makefile
+++ src/hg/utils/otto/uniprot/makefile
@@ -1,21 +1,21 @@
 PREFIX=/hive/data/outside/otto/uniprot/
 PARSER=../../../../utils/uniprotToTab
 
 # Everything doUniprot calls at run time must be listed here, or the copy that cron runs
 # is not under version control. refs #38300
-FILES=*.sh *.as doUniprot mapUniprot_doBlast pslProtCnv trackDb.template.txt README.txt ${PARSER}
+FILES=*.sh *.as doUniprot miniprotToPsl mapUniprot_doBlast pslProtCnv trackDb.template.txt README.txt ${PARSER}
 
 # Not "pylint -E": every pylint on hgwdev is itself pinned to a python that no longer
 # exists, so "make install" used to abort on its first line.
 SYNTAXCHECK=python3 -c 'import ast,sys; [ast.parse(open(f).read(), f) for f in sys.argv[1:]]'
 
 install:
 	${SYNTAXCHECK} doUniprot ${PARSER}
 	-pylint -E doUniprot ${PARSER} # best effort, do not stop the install if pylint is broken
 	rsync -avpu ${FILES} ${PREFIX} # never overwrite newer files
 
 # rsync -u above refuses to overwrite a file that is newer in ${PREFIX} than in the tree,
 # so a hand edit of the live copy wins silently and forever. Run this to see what has
 # drifted apart.
 diff:
 	@for f in ${FILES}; do \