af613a331e6839c6513c3e366abcb67af0fe8386
max
  Wed Sep 9 06:47:14 2026 -0700
UniProt otto: get the monthly update running again and make a stalled run visible

The monthly UniProt job had produced nothing since January 2025. The tracks
served release 2024_06 while the download sitting on disk was at 2026_02, on
every assembly the job builds.

Cause: uniprotToTab appended a personal conda site-packages directory to
sys.path, and doUpdate.sh sourced a virtualenv, both built for python 3.6. A
venv's python is only a symlink to the system one, so when hgwdev moved to
python 3.9 the compiled lxml in there stopped loading and every run died at the
parse step. Removed both. The system python3 has lxml from python3-lxml and the
two are upgraded together, so there is nothing left here to go stale. Verified
by parsing real 2026_02 records under python 3.9 with lxml 5.4.

Why nobody noticed for nineteen months:

- doUpdate.sh read $? after an intervening echo, so it captured the echo's exit
code and mailed "Big Uniprot update OK" every month while the job was dying.
It now reads the real exit code, says FAILED, prints the tail of the log and
exits nonzero. A month with no new UniProt release stays silent, which is the
normal otto behaviour, so silence again means "nothing to do".
- The logs were overwritten on every run, so a failure left no trace on disk.
doUpdate.sh now appends one line per run to runLog.txt, which is never
truncated, and keeps a failing log as lastFail.log.
- version.txt in each bigBed directory was rewritten on every run even when the
release string was identical. That is the file the trackDb dataVersion setting
shows, and its date is what people check to decide whether a pipeline is still
alive, so a stalled track could look freshly updated. It is now written only
when the release actually changes.

Also, so this cannot come back:

- doUniprot checks that uniprotToTab can start before the download, instead of
finding out 35 minutes later.
- pylint on hgwdev is itself pinned to pythons that no longer exist, so
"make install" aborted on its first line and could not be used. Replaced with
a syntax check that needs nothing but python3; pylint stays best-effort.
- uniprotToTab, pslProtCnv, trackDb.template.txt and README.txt ran from
/hive/data/outside/otto/uniprot without being in the makefile's copy list.
The tree copy of uniprotToTab was still python 2 from 2021. All are now
listed and in sync, and "make diff" reports drift.
- Brought the two live-only fixes into the tree: mkdir -p in makeUniProtPsl.sh
and the pslMap -inType/-mapType flags.

refs #38300

diff --git src/hg/utils/otto/uniprot/makeUniProtPsl.sh src/hg/utils/otto/uniprot/makeUniProtPsl.sh
index c2166844c38..61090c7db74 100755
--- src/hg/utils/otto/uniprot/makeUniProtPsl.sh
+++ src/hg/utils/otto/uniprot/makeUniProtPsl.sh
@@ -1,152 +1,154 @@
 #!/bin/bash
 # create mapping from UniProt protein sequences to genome
 # input is protein sequence file and transcript table, and optional a mapping proteinId <-> transcriptId
 
 # originally inspired/copied from Markd's script LS-SNP pipeline
 # original: snpProtein/build/Makefile
 
 # for human, this is now using RefSeq, see #27836
 # The reason is that for all proteins that are broken in the reference, the positions will be wrong 
 # for Gencode, but not with RefSeq (for which we have a proper PSL alignment).
 
 # parameters:
 # $1 = the fasta file with uniprot sequences
 # $2 = the transcript fasta file
 # $3 = the transcript->genome PSL file
 # $4 = MINALI, the minimum percent ID, e.g. 0.93
 # $5 = the parasol cluster
 # $6 = temporary workdir
 # $7 = OUTPUT: the nucleotide PSL output file to create
 # $8 = optional: tsv table with mapping from uniprot to transcript
 
 # Will always rm -rf the work directory, before and after a run
 
 if [ "$1" == "" ]; then
         echo Please specify a uniprotToTab fasta input file and a db, a gene table, a cluster name and the output file name
         exit 0
 fi
 
 # ideally, we have a table of uniprotID <-> gene model ID, so can decide where to map protein sequences that match identically twice.
 # but often we don't have known genes tables.
 # to get an idea of the impact, I compared hg38 with and without known gene tables
 #
 # 2385 out of 38931 PSLs are multi-mapping
 
 # most of them are mapping twice, but some 35 times
 #   2 ************************************************************ 1477
 #   3 ********** 238
 #   4 ***** 126
 #   5 **** 109
 #   6 *** 84
 #   7 ****** 142
 #   8 *** 82
 #   9  10
 #  10 *** 66
 #  11  1
 #  12  3
 #  13  6
 #  14  6
 #  15  6
 #  16  6
 #  17  7
 #  18  3
 #  19  3
 #  20  3
 # <minVal or >= 21  7
 
 # worst ones:
 #A6NER0  20      0.042%
 #Q99706-6        20      0.042%
 #P43629-2        20      0.042%
 #P43630-2        28      0.059%
 #P43630-1        28      0.059%
 #Q99706-3        30      0.063%
 #Q99706-4        30      0.063%
 #Q99706-2        30      0.063%
 #Q99706-1        30      0.063%
 #Q8N743  32      0.067%
 
 # the input fasta file has to be created by the uniprotToTab parser (originally from the publications track code)
 UNIPROTFAGZ=$1
 TRANSCRIPTFA=$2
 TRANSCRIPTPSL=$3
 MINALI=$4
 CLUSTER=$5
 WORKDIR=$6
 OUTFNAME=$7
 PAIRNAME=$8
 
 #if [[ "$DB" == "ci3" ]]; then
    #MINALI=0.85
 #fi
 
 # if you change BLASTDIR, also must change the cluster script mapUniprot_doBlast
 BLASTDIR=/cluster/bin/blast/x86_64/blast-2.2.16/bin
 
 # stop on errors
 set -e
 # show commands
 set -x 
 
 
 #WORKDIR=makeUniProtPsl-$2-$3-$4.tmp
 
 #rm -rf $WORKDIR
 
 mkdir -p $WORKDIR
 
 cp ${UNIPROTFAGZ} $WORKDIR/uniProt.fa
 cp $TRANSCRIPTFA $WORKDIR/transcripts.fa
 cp $TRANSCRIPTPSL $WORKDIR/transcripts.psl
 if [ "$PAIRNAME" != "" ] ; then
    cp $PAIRNAME $WORKDIR/upToTrans.pairs
 fi
 
 # setup blast uniprot -> transcript cDNAs
 if [ -f $WORKDIR/bestAln.psl ] ; then
         echo WARNING: re-using existing protein-transcript alignments to save time! see $WORKDIR/bestAln.psl
 else
-        mkdir $WORKDIR/queries 
-        mkdir $WORKDIR/aligns 
+        mkdir -p $WORKDIR/queries
+        mkdir -p $WORKDIR/aligns
         faSplit about $WORKDIR/uniProt.fa 2500 $WORKDIR/queries/
         ${BLASTDIR}/formatdb -i $WORKDIR/transcripts.fa -p F
 
         # create joblist and run
         set +x # quiet for now
         >$WORKDIR/jobList
         for i in $WORKDIR/queries/*.fa; do
                 echo "mapUniprot_doBlast transcripts.fa queries/`basename $i` {check out exists aligns/`basename $i .fa`.psl}" >> $WORKDIR/jobList
         done; 
         set -x
         cp mapUniprot_doBlast $WORKDIR/
         ssh $CLUSTER "cd `pwd`/$WORKDIR && para make jobList"
         echo Concatenating and filtering protein/transcript alignments
         # sort and pick the best alignments for each protein
         find $WORKDIR/aligns -name '*.psl' | xargs cat | pslReps -noIntrons -nohead -nearTop=0.01 -minAli=$MINALI stdin stdout /dev/null > $WORKDIR/bestAln.psl
 fi
 
 # when we have a list of pairs to filter on, then use it
 if [ -f $WORKDIR/upToTrans.pairs ]; then
     # these are very special pslSelect options:
     # - pass through all PSLs for queries that do not appear in our tables (=all Trembl and a few swissprot sequences missed by our tables)
     # - only compare the part before the dot
     # - use pslSwap as a workaround, as pslSelect doesn't have -tDelim yet and the query IDs don't have dots in them anyways
     #   too lazy to patch pslSelect now
     cat $WORKDIR/upToTrans.pairs | gawk '{OFS="\t"; print $2, $1; }' > $WORKDIR/transToUp.pairs
     pslSwap $WORKDIR/bestAln.psl stdout | pslSelect -qPass -qDelim=. -qtPairs=$WORKDIR/transToUp.pairs stdin stdout | pslSwap stdin $WORKDIR/uniProtVsTranscripts.psl
     # when we have no known gene table, blat the proteins directly and
     # take the top 1% alignments with 93% query coverage. Hopefully that gives similar results...
     # inspired by kent/src/hg/makeDb/doc/ucscGenes/hg19.ucscGenes13.csh
     # as requested by Alejo: lower to 93%, which are the pslReps defaults
 else
     cp $WORKDIR/bestAln.psl $WORKDIR/uniProtVsTranscripts.psl
 fi
 
 # now combine the two alignments with pslMap
-pslMap $WORKDIR/uniProtVsTranscripts.psl $WORKDIR/transcripts.psl $WORKDIR/uniProtVsGenome.psl -mapInfo=$WORKDIR/mapInfo.tab
+# the query is protein and the target is nucleotide, so pslMap has to be told the types,
+# otherwise it guesses and the block sizes come out in the wrong units
+pslMap $WORKDIR/uniProtVsTranscripts.psl $WORKDIR/transcripts.psl $WORKDIR/uniProtVsGenome.psl -mapInfo=$WORKDIR/mapInfo.tab -inType=prot_na -mapType=na_na
 # 2016: lowering to 95% identity due to hg38 alt loci sucking up our main (and more important) alignments from the
 # 2021: using MINALI is more consistent
 # normal chromosomes
 sort -k10,10 $WORKDIR/uniProtVsGenome.psl | pslCDnaFilter stdin -globalNearBest=$MINALI  -bestOverlap -filterWeirdOverlapped stdout | sort | uniq > $OUTFNAME
 cp $WORKDIR/mapInfo.tab ${OUTFNAME}.mapInfo
 #rm -rf $WORKDIR