af670dcecc6911452f3e1a5c38cfea1ef0067978 max Wed Sep 9 07:31:07 2026 -0700 UniProt otto: submit the parasol batch on hgwdev, drop the ssh to ku The protein-to-transcript BLAST batch was submitted with ssh ku "cd <workdir> && para make jobList" ku has been decommissioned for about two years. It still resolves in DNS, but "ssh ku" is an immediate "No route to host", so a run would have died at the mapping stage, which comes only after days of XML parsing. hgwdev is the parasol head node now, so "para make" runs here with no ssh hop. Verified with a one-job batch using the same bare-command jobList this script writes: it lands on a compute node with its cwd set to the batch directory, the same as the old ku jobs did. Also checked from a compute node that tclsh and the blast-2.2.16 blastall/formatdb that mapUniprot_doBlast needs are still there. The cluster name is gone rather than redirected. doUniprot no longer reads /cluster/bin/scripts/cluster.txt, a file written in 2017 that still says "ku" and that nothing else in the tree read, and makeUniProtPsl.sh no longer takes a head node argument, so its positional parameters shift down by one. refs #38300 diff --git src/hg/utils/otto/uniprot/doUniprot src/hg/utils/otto/uniprot/doUniprot index b62b82036e3..bcbdb3cc1f0 100755 --- src/hg/utils/otto/uniprot/doUniprot +++ src/hg/utils/otto/uniprot/doUniprot @@ -55,34 +55,30 @@ "wuhCor1", # no monkeypox yet "mpxvRivers", # do not run on hs1, has no good gene model and BLAT protein takes forever "hs1", # sea urchin doesn't have any normal gene model track, not even Augustus "strPur2" ] # the script will stop if the difference between the last uniprot version and the current version is > than <maxDiffVersions> # but this check is only done for some important assemblies. For fringe assemblies, the numbers can easily # change 10x, as trembl predictions are the majority of annotations maxDiffVersions = 5.0 checkDbs = ["hg19", "hg38", "mm10", "rn6", "danRer10", "danRer7", "ce11"] -# when we rename the parasol cluster, no need to change this file -# probably overkill in retrospect but seemed like a good idea back then -clusterFname = "/cluster/bin/scripts/cluster.txt" - TREMBLCOLOR="0,150,250" # light blue SWISSPCOLOR="12,12,120" # dark blue asDir = dirname(__file__) # some feature types should not go into the bed name field. # for these features, we use the 'comment' as the bed name # e.g. "region of interest" is not very interesting, the actual # description is usually much more interesting. useComment = set(["domain", "chain","region of interest","topological domain","short sequence motif"]) # mostly for the "other annotations" subtrack, to show the different types featTypeColors = { "modified residue" : "200,200,0", "glycosylation site" : "0,100,100", @@ -1756,31 +1752,31 @@ geneId, transList = line.rstrip("\n").split('\t') transList = transList.split(',') g2refseq[geneId] = transList return g2refseq def parseKeyVal(fname): " given a file with key<tab>val return a dict with key -> list of values " if fname is None: return None vals = defaultdict(list) for line in open(fname): key, val = line.rstrip("\n").split("\t") vals[key].append(val) return vals -def makeProteinGenomePsl(taxId, protFa, tabFnames, db, cluster, doForce, mapDir, mapDescFname): +def makeProteinGenomePsl(taxId, protFa, tabFnames, db, doForce, mapDir, mapDescFname): """ map proteins of a taxon stored in protFa to db via geneTable and create a PSL file for the mapping Will try to reuse the old PSL mapping file if possible, as it is expensive to build. Need to BLAST the entire file: BLAST has no option to subset the target (unlike blast or lastz), verified with an email to NCBI BLAST support team. A central point of this function is buildSelectFile(). The psl select pair files have a huge influence on the result, so it is rebuilt every time, as any change in UniProt will change it. It's also quick to create. Its MD5 is part of the PSL mapping file. Any change to the pslSelect pair file will trigger a realignment/rebuild of the pslMap PSL file. """ # use MD5s of the input files to determine if the protein -> genome map has to be rebuilt global MINALI protMd5 = fastaMd5(protFa) stats = OrderedDict() stats["user"] = os.getlogin() @@ -1839,31 +1835,31 @@ stats["metaMd5"] = metaMd5[:10] stats["transMd5"] = transMd5[:10] stats["fullMd5"] = fullMd5[:10] stats["mapFname"] = basename(mapFname) logging.debug("%s does not exist" % mapFname) if geneTable=="blat": logging.error("Could not find any gene table for %s, using BLAT to map proteins" % db) blatProteinsKeepBest(protFa, db, mapFname, stats) else: if db in ["ci3"]: # wow: ciona's genome is from a different organism than the refseq transcripts. MINALI = 0.85 stats["minAli"] = MINALI workDir = "clusterRun-map-%s-%s-%s.tmp" % (db, geneTable, fullMd5) - scriptArgs = [protFa, transcriptFa, transcriptPsl, str(MINALI), cluster, workDir, mapFname] + scriptArgs = [protFa, transcriptFa, transcriptPsl, str(MINALI), workDir, mapFname] if selectFname is not None: scriptArgs.append(selectFname) # This is where the BLAST alignment-meat happens cmd = ["time", "./makeUniProtPsl.sh"] cmd.extend(scriptArgs) logging.info("Running %s" % cmd) run(cmd) writeMapDesc(stats, db, geneTable, mapDescFname) assert(os.path.getsize(mapFname)!=0) # if you modify this line, also modify the other return statement above return mapFname, geneTable, fullMd5, protMapSource, protToTrans, True @@ -2029,34 +2025,30 @@ myDir = dirname(__file__) # directory of doUniprot # takes an hour or so cmd="%s/uniprotToTab %s %s %s" % (myDir, uprotDir, taxIdStr, tabDir) run(cmd) # takes 2-3 days on hgwdev - big data and XML don't mix... cmd="%s/uniprotToTab %s %s %s --trembl" % (myDir, uprotDir, taxIdStr, tabDir) run(cmd) writeReleaseString(uprotDir, tabDir) # if the uniprot update changed the sequences, update the corresponding pslMap files of that genome logging.info("checking/creating pslMap files") run("mkdir -p %s" % mapDir) - # parasol cluster name - assert(isfile(clusterFname)) # not running at UCSC? - cluster = open(clusterFname).read().strip() - # get the uniProt version for the trackVersion table that we will update later relFname = join(tabDir, "version.txt") versionString = open(relFname).read() for taxId, dbs in taxIdDbs.items(): if onlyDbs is not None and len(set(dbs).intersection(onlyDbs))==0: continue logging.info("Working on taxon ID %d" % taxId) faFnames = [ ("swissprot", join(tabDir,"swissprot.%d.fa.gz" % taxId)), ] if doTrembl: faFnames.append( ("trembl", join(tabDir,"trembl.%d.fa.gz" % taxId)) ) @@ -2075,31 +2067,31 @@ # find the best gene table for each database, create a mapping # protein -> genome and lift the uniprot annotations to bigBed files for db in dbs: if onlyDbs is not None and db not in onlyDbs: continue logging.debug("Annotating assembly %s" % db) chromSizesFname = "/hive/data/genomes/%s/chrom.sizes" % db accToMeta = parseRecordMeta(tabDir, taxId, doTrembl) dbMapDir = makeSubDir(mapDir, db) mapDescFname = join(dbMapDir, "liftInfo.json") mapFname, geneTable, mapMd5, accToMapSource, accToTrans, isNewMap = \ - makeProteinGenomePsl(taxId, fullFaFname, tabFnames, db, cluster, options.force, dbMapDir, mapDescFname) + makeProteinGenomePsl(taxId, fullFaFname, tabFnames, db, options.force, dbMapDir, mapDescFname) dbBigBedDir = makeSubDir(bigBedDir, db) convMapToBigPsl(fullFaFname, mapFname, dbBigBedDir, accToDb, accToMeta, accToMapSource, accToTrans, chromSizesFname, doTrembl) #shutil.copy(mapDescFname, mapDir) if options.onlyMap: continue annotFnames = ["tab/swissprot.%d.annots.tab" % taxId] if doTrembl: annotFnames.append( "tab/trembl.%d.annots.tab" % taxId ) uniprotLift(fullFaFname, annotFnames, chromSizesFname, mapFname, dbBigBedDir, accToDb, accToMeta, options)