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 @@ -1,154 +1,166 @@ #!/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 +# $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 # 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 # = 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 WORKDIR=$5 OUTFNAME=$6 PAIRNAME=$7 #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 -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 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 # 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