dcf390fefe68c32f26bd32717bd858190e8875a0 max Fri Sep 11 06:06:58 2026 -0700 uniprot otto: build the transcript files for a GenArk assembly from the hub A GenArk assembly has no MySQL database, so everything makeTranscriptFiles did had to come from somewhere else. It all exists in the hub: bigGenePredToGenePred gene bigBed -> transcripts.gp genePredToFakePsl -chromSize=<hub> -> transcripts.psl getRnaPred -genomeSeqs=<2bit> -> transcripts.fa Both of those options exist precisely so the tools do not have to go through chromInfo, so no database is involved anywhere in the chain. Also made hub-aware, in the same call chain: - chrom.sizes for the bigPsl and bigBed steps, which used /hive/data/genomes/<db>/chrom.sizes and there is no such directory - getTransIds, which read the transcript IDs out of a MySQL table; the fasta we just built holds them all, so read them from there as the refGene case does - writeMapDesc, which read ncbiRefSeqVersion.txt; date the models by the gene bigBed we actually read instead Helpers genArkTwoBit and genArkChromSizes glob for the single matching file and abort if there is not exactly one, since the hub names those with the short accession while the files under bbi/ carry the full asmId. Verified on bonobo mPanPan1_v2.0: 95418 transcripts, identical ID sets across the fasta, genePred and PSL; getTransIds returns those 95418 without touching MySQL; writeMapDesc dates the models 2025-06-02, which matches the date embedded in the build's own ncbiRefSeq GTF filename. hg38 still reads its real ncbiRefSeqVersion.txt and is otherwise untouched. Still to do: the outputs. /gbdb/<db>/uniprot and a trackDb .ra stanza are not how a GenArk assembly is served. refs #38300 diff --git src/hg/utils/otto/uniprot/doUniprot src/hg/utils/otto/uniprot/doUniprot index 81963a28f3e..290b2efc3ca 100755 --- src/hg/utils/otto/uniprot/doUniprot +++ src/hg/utils/otto/uniprot/doUniprot @@ -1014,30 +1014,47 @@ ("augustus", "bbi/*.augustus.bb"), ] dbIsHubCache = {} def genArkHubDir(db): """ return the /gbdb/genark directory for a GenArk assembly, or None for a classic db. dbDb.nibPath is "hub:/gbdb/genark/GCF/029/289/425/GCF_029289425.2" for these. """ if db not in dbIsHubCache: rows = list(runQuery("hgcentral", "select nibPath from dbDb where name='%s'" % db, usePublic=True)) nibPath = rows[0][0] if rows else "" dbIsHubCache[db] = nibPath[len("hub:"):] if nibPath.startswith("hub:") else None return dbIsHubCache[db] +def genArkFile(db, hubDir, fileGlob, what): + " return the single file matching fileGlob inside a GenArk hub directory " + matches = glob.glob(join(hubDir, fileGlob)) + if len(matches)!=1: + errAbort("%s: expected exactly one %s in %s, found %d: %s" % \ + (db, what, hubDir, len(matches), matches)) + return matches[0] + +def genArkTwoBit(db, hubDir): + " the assembly sequence of a GenArk assembly " + # named with the short accession, unlike the files under bbi/ + return genArkFile(db, hubDir, "*.2bit", "2bit file") + +def genArkChromSizes(db, hubDir): + " the chrom.sizes of a GenArk assembly. Note the .txt suffix, which /hive does not use " + return genArkFile(db, hubDir, "*.chrom.sizes.txt", "chrom.sizes file") + def findBestGeneBigBed(db, hubDir): """ find the best gene model bigBed for a GenArk assembly. Returns (sourceName, fileName), or (None, None) if the hub has no gene models at all. """ for name, fileGlob in genArkGeneSources: matches = sorted(glob.glob(join(hubDir, fileGlob))) if matches: if len(matches) > 1: logging.warning("%s: %d files match %s, using %s" % (db, len(matches), fileGlob, matches[0])) logging.info("%s: best gene models are %s, from %s" % (db, name, matches[0])) return name, matches[0] logging.warning("%s: no gene models found in %s, will have to BLAT the proteins" % (db, hubDir)) return None, None def findBestGeneTable(db): @@ -1471,32 +1488,34 @@ #nuclFname = mapFname+".nucl" # originally used "| {myDir}/pslProtCnv " but tried to get rid of the markd-script dependency #pslToProtPsl(nuclFname, mapFname) # then later figured that pslToProtPsl() is not the same as pslProtCnv and went back the MarkD way # then in 2024 finally used pslProtToRnaCoords #os.remove(nuclFname) run(cmd) def getTransIds(db, geneTable, transcriptFa): """ return all possible transcript IDs given a gene table. The reason that we need is that Uniprot often refers to transcripts that don't exist and we are using a select file to map only to those. In these cases, we want to make sure that our select file contains only transcript that we actually have, so we can log how many transcript can possibly mapped (and which ones). This is important for debugging. """ # no idea how to get the version suffix out of our mysql tables - thanks Obama! - if geneTable=="refGene": - logging.debug("Gene table is refGene, reading IDs from %s" % transcriptFa) + # A GenArk assembly has no MySQL tables, but the fasta we built from its gene bigBed + # holds every transcript ID already, so read them from there as refGene does. + if geneTable=="refGene" or genArkHubDir(db) is not None: + logging.debug("Reading transcript IDs from %s" % transcriptFa) allTrans = set(parseFasta(open(transcriptFa)).keys()) return allTrans logging.info("Getting transcript IDs for db=%s, table=%s, field=name" % (db, geneTable)) sql = "SELECT DISTINCT name from "+geneTable cmd = ["hgsql", db, "-NBe", sql] proc = subprocess.Popen(cmd, encoding="latin1", stdout=PIPE) allTrans = set() for line in proc.stdout: idStr = line.rstrip() allTrans.add(idStr) logging.info("Got %d distinct transcript IDs from MySQL table" % len(allTrans)) return allTrans def writeTransToSelect(upSeqIds, transIds, shortToLong, notFound, outTransIds, commonTrans, \ @@ -1721,30 +1740,49 @@ def makeTranscriptFiles(db, geneTable, taxId, dbDir): """ given a db and a gene table, create two files, one for the transcript fasta sequences and one for the transcript -> genome PSL. Return tuple of the file names ALWAYS REMOVES _ALT/_FIX/_HAP results from the PSL! """ #dbDir = join(mapDir, db+"_"+geneTable) #if not isdir(dbDir): #os.makedirs(dbDir) faName = join(dbDir, "transcripts.fa") pslName = join(dbDir, "transcripts.psl") pslFields = "matches,misMatches,repMatches,nCount,qNumInsert,qBaseInsert,tNumInsert," \ "tBaseInsert,strand,qName,qSize,qStart,qEnd,tName,tSize,tStart,tEnd,blockCount,blockSizes,qStarts,tStarts" + hubDir = genArkHubDir(db) + if hubDir is not None: + # GenArk assembly: no MySQL tables at all, the gene models are a bigGenePred in the + # hub. bigGenePredToGenePred gets us back to a genePred, and both of the tools that + # follow can work off the hub's own chrom.sizes and 2bit rather than chromInfo, so + # nothing here needs a database. + geneBb = findBestGeneBigBed(db, hubDir)[1] + gpName = join(dbDir, "transcripts.gp") + cdsName = join(dbDir, "transcripts.cds") + + run(["bigGenePredToGenePred", geneBb, gpName]) + run(["genePredToFakePsl", "-chromSize=%s" % genArkChromSizes(db, hubDir), + "noDb", gpName, pslName, cdsName]) + run(["getRnaPred", "-genomeSeqs=%s" % genArkTwoBit(db, hubDir), + "noDb", gpName, "all", faName]) + + logging.info("Created %s and %s for %s/%s from %s" % (faName, pslName, db, geneTable, geneBb)) + return faName, pslName + if geneTable=="ncbiRefSeq": origFa = "/gbdb/%s/ncbiRefSeq/seqNcbiRefSeq.rna.fa" % db shutil.copyfile(origFa, faName) pslTable = "ncbiRefSeqPsl" query = "SELECT "+pslFields+" from "+pslTable+" WHERE tName NOT LIKE '%_hap%' AND tName not like '%_alt%' AND tNAME NOT LIKE '%_fix%'" runQueryToFile(db, query, pslName) elif geneTable=="refGene": # ChrisL told me where these files can be found #cmd = "curl https://hgdownload.cse.ucsc.edu/goldenpath/%s/bigZips/refMrna.fa.gz | zcat > %s" % (db, faName) cmd = "zcat /hive/data/outside/genbank/data/ftp/%s/bigZips/refMrna.fa.gz | tr ' ' '.' > %s" % (db, faName) run(cmd) pslFields = "matches,misMatches,repMatches,nCount,qNumInsert,qBaseInsert,tNumInsert," \ "tBaseInsert,strand,concat(qName, '.', gbSeq.version),qSize,qStart,qEnd,tName,tSize,tStart,tEnd,blockCount,blockSizes,qStarts,tStarts" @@ -1755,31 +1793,38 @@ query = "SELECT name, chrom, strand, txStart, txEnd, cdsStart, cdsEnd, exonCount, exonStarts, exonEnds from %s where cdsStart < cdsEnd" % (geneTable) gpName = join(dbDir, "transcripts.gp") cdsName = join(dbDir, "transcripts.cds") runQueryToFile(db, query, gpName) cmd = ["genePredToFakePsl", db, gpName, pslName, cdsName] run(cmd) cmd = ["getRnaPred", db, gpName, "all", faName] run(cmd) logging.info("Created %s and %s for %s/%s" % (faName, pslName, db, geneTable)) return faName, pslName def writeMapDesc(stats, db, geneTable, mapDescFname): " write a little json file that points to the map file used and some basic stats about it " logging.debug("making json file with version and other lift info") - if geneTable == "ncbiRefSeq": + hubDir = genArkHubDir(db) + if hubDir is not None: + # A GenArk assembly has no ncbiRefSeqVersion.txt. Date the gene models by the + # bigBed we actually read, which is what the user needs to tie the mapping to a + # particular annotation; the source name is carried separately in geneTable. + geneBb = findBestGeneBigBed(db, hubDir)[1] + version = datetime.datetime.fromtimestamp(os.path.getmtime(geneBb)).strftime('%Y-%m-%d') + elif geneTable == "ncbiRefSeq": version = open("/gbdb/"+db+"/ncbiRefSeq/ncbiRefSeqVersion.txt").read().strip() elif geneTable in ['refGene', 'augustusGene', 'knownGene']: version = datetime.datetime.today().strftime('%Y-%m-%d') elif geneTable=="blat": version = "direct" elif geneTable in ["ensGene"]: version = list(runQuery("hgFixed", "select version from trackVersion where db='%s' and name='ensGene' order by ix desc limit 1;" % db))[0][0] else: assert(False) stats["geneTable"] = geneTable stats["version"] = version stats["createdDate"] = datetime.datetime.today().strftime('%Y-%m-%d') with open(mapDescFname, "w") as mapDescFh: @@ -2109,30 +2154,36 @@ logging.warn("File %s is empty. Taxon ID %s does not have any UniProt annotations." \ " Skipping this organism." % (fullFaFname, taxId)) continue tabFnames = ["tab/swissprot.%d.tab" % taxId] if doTrembl: tabFnames.append( "tab/trembl.%d.tab" % taxId ) # 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) + # a GenArk assembly has no /hive/data/genomes directory, its chrom.sizes + # lives in the hub + hubDir = genArkHubDir(db) + if hubDir is not None: + chromSizesFname = genArkChromSizes(db, hubDir) + else: 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, options.force, dbMapDir, mapDescFname) dbBigBedDir = makeSubDir(bigBedDir, db) convMapToBigPsl(fullFaFname, mapFname, dbBigBedDir, accToDb, accToMeta, accToMapSource, accToTrans, chromSizesFname, doTrembl) #shutil.copy(mapDescFname, mapDir)