17fb5d876fd59c39645d7589ed115e23208f8743
max
  Thu Sep 17 05:27:00 2026 -0700
uniprot otto: fewer, larger BLAST jobs

#Preview2 week - bugs introduced now will need a build patch to fix
The protein fasta was split into 2500 byte chunks regardless of its size, which for
zebrafish meant 18894 cluster jobs averaging 13 seconds each. At that length parasol
overhead and the creation of 18894 tiny result files cost more than the BLAST does.

The cluster was never the problem: it finished 57 CPU hours of zebrafish alignment
in 12 minutes of wall clock, a speedup of about 290. The cost is afterwards, in the
single-threaded "find aligns | xargs cat | pslReps" that has to open every one of
those files again over NFS, which is the slowest part of a per-assembly run.

Split for a job count instead, about a thousand, keeping 2500 bytes as a floor so a
small protein set still splits sensibly. Zebrafish goes from 18894 jobs to about
1000, each around five minutes, with the same total CPU and a twentieth of the files
to concatenate.

Also corrected the argument documentation at the top, which still described the
parasol cluster argument that went away with ku.

refs #38300

diff --git src/hg/utils/otto/uniprot/makeUniProtPsl.sh src/hg/utils/otto/uniprot/makeUniProtPsl.sh
index f28828c7131..f913e590f8d 100755
--- src/hg/utils/otto/uniprot/makeUniProtPsl.sh
+++ src/hg/utils/otto/uniprot/makeUniProtPsl.sh
@@ -2,34 +2,34 @@
 # 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
+# $5 = temporary workdir
+# $6 = OUTPUT: the nucleotide PSL output file to create
+# $7 = optional: tsv table with mapping from uniprot to transcript
+# (there used to be a parasol cluster argument here, removed when ku went away)
 
 # 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
@@ -95,31 +95,43 @@
 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 -p $WORKDIR/queries
         mkdir -p $WORKDIR/aligns
-        faSplit about $WORKDIR/uniProt.fa 2500 $WORKDIR/queries/
+        # Aim for a job count rather than a fixed chunk size. At a flat 2500 bytes this made
+        # 18894 jobs for zebrafish that averaged 13 seconds each, so parasol overhead and
+        # creating 18894 tiny result files cost more than the BLAST did. The cluster is not
+        # the bottleneck either way - it absorbed 57 CPU hours in 12 minutes - but every one
+        # of those files then has to be opened again by the single-threaded pslReps below,
+        # over NFS, which is the slowest part of the whole per-assembly run. refs #38300
+        targetJobs=1000
+        faBytes=`stat -c %s $WORKDIR/uniProt.fa`
+        chunkSize=`expr $faBytes / $targetJobs`
+        # keep the old size as a floor, so a small protein set still splits sensibly
+        if [ $chunkSize -lt 2500 ]; then chunkSize=2500; fi
+        echo "splitting `expr $faBytes / 1000000` MB of protein into chunks of $chunkSize bytes"
+        faSplit about $WORKDIR/uniProt.fa $chunkSize $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/
         # hgwdev is the parasol head node, so "para make" here talks to the hub directly
         ( cd $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