4decf5fbe82051b2d5aff9cabce3feae9dfafae1
max
  Sat Sep 26 18:48:55 2026 -0700
uniprot: stop declaring the alignments as amino acid coordinates, and show them properly

The bigPsl seqType field describes the coordinates, not the letters stored beside
them, and the UniProt query side is in bases: these proteins reach the genome
through transcripts, so a query runs three bases to a residue. Declaring amino
acids made pslFromBigPsl divide the block sizes by three and leave the query
coordinates alone, so every reader got an alignment measured in two units at once.
That fed a heap overflow in the alignment page, drew blocks short in hgTracks, and
loaded sub-codon blocks as size 0, which aborted pslTransMap and took down the
lifted SwissProt track (#38249).

pslProtFromNaLike() converts such a psl to one counted in residues. Blocks are
trimmed to whole codons, and a residue whose codon straddles an exon junction sits
in two places in the genome at once, so it gets no column and the page says how
many are missing rather than dropping them silently: 0.9% of residues, though 87%
of alignments have at least one. A bigPsl also keeps the reference strand where a
psl reads the query strand, so minus-strand items arrived claiming the protein was
reversed and were rendered as reverse complemented nucleotide ambiguity codes;
pslRc moves them to the convention blat uses, query forward and the strand on the
target.

Checked on hg38 against the translated genome: 210,416 residues over both strands,
99.86% identical, the remainder real protein-vs-reference variation. All 29 items
of a test region render the stored protein at the right residues. Files already
published still say amino acid and must keep working until they are rebuilt, so the
conversion also requires the blocks to measure the target the way the target is
measured; verified that separates the two shapes on 3000 records each way, and that
all 29 render without crashing in the old format, where they now say plainly that
the coordinates and the sequence do not match.

refs #38300

diff --git src/hg/utils/otto/uniprot/doUniprot src/hg/utils/otto/uniprot/doUniprot
index 422749f44a0..e0b65d36ad3 100755
--- src/hg/utils/otto/uniprot/doUniprot
+++ src/hg/utils/otto/uniprot/doUniprot
@@ -1558,31 +1558,38 @@
     # - split into swissprot and trembl files
     bpInputSwissFname = mapFname.replace(".psl", ".swissprot.pslInput")
     bpInputTremblFname = mapFname.replace(".psl", ".trembl.pslInput")
 
     ofhSwiss = open(bpInputSwissFname, "w")
 
     ofhTrembl = None
     if bigPslTremblFname:
         ofhTrembl = open(bpInputTremblFname, "w")
 
     protToLocs = defaultdict(list)
     for line in open(bpInputFname):
         # chr1    11994   12144   B4E2Z4  1000    +       11994   12144   0       1       150,    0,      186     336     +       336     186,          249250621        150     0       0       0       0
         row = line.rstrip("\n").split("\t")
 
-        row[-1] = "2" # seqType field: 0=empty, 1=nucleotide, 2=amino_acid
+        # The seqType field (0=empty, 1=nucleotide, 2=amino_acid) used to be stamped "2" here.
+        # It describes the coordinates, not the letters in oSequence, and ours are nucleotides:
+        # these proteins reach the genome through transcripts, so a query is three bases per
+        # residue.  Declaring amino acids made pslFromBigPsl() divide the block sizes by three
+        # and leave the query coordinates alone, which left every reader with an alignment whose
+        # two sides were in different units - a heap overflow in the alignment page, blocks drawn
+        # short in hgTracks, and blocks shorter than a codon loading as size 0, which aborted
+        # pslTransMap and took the lifted SwissProt track down (#38249).  refs #38300
 
         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: