befc3bcf5e350b84ccce24c7aba1acbb0f073c4f max Sat Sep 19 09:07:54 2026 -0700 uniprot otto: per-taxon NCBI tables, and lift info next to a reused mapping Two failures on the GenArk run, both of which would have hit about a quarter of the 664 taxa. The per-taxon NCBI gene to RefSeq tables were only built inside the download step, which --skipDownload turns off. Those tables are keyed on the plan rather than on the UniProt release, so a plan with taxa that have never been run needs them even when nothing has to be downloaded: 104 existed, for the old dbDb plan of 109 taxa, against 664 here. The download and the split are now separate, and the split runs whenever a taxon in the plan has no table yet. Separating them turned up a third thing. The download unpacked to ncbi/gene2refseq.tsv, 16 GB, while the split read ncbi/gene2refseq.gz, so the two named different files and the tables were being built from whatever .gz happened to be on disk. That was a copy from last November. It now downloads to the file the split reads, through a .tmp so an interrupted download cannot be mistaken for a complete one, and stays gzipped, which also saves the 16 GB. Separately, makeProteinGenomePsl writes the lift info only when it builds the mapping, not when it reuses one. A run interrupted between building the PSL and writing that file leaves the PSL on its own, and the next run reuses the PSL and then fails copying a file that was never written. Twenty assemblies were in that state after the aborted start earlier today. It is now written on the reuse path as well, when it is missing. refs #38300 diff --git src/hg/utils/otto/uniprot/doUniprot src/hg/utils/otto/uniprot/doUniprot index a47eff5e64e..140fd7bb29b 100755 --- src/hg/utils/otto/uniprot/doUniprot +++ src/hg/utils/otto/uniprot/doUniprot @@ -2317,30 +2317,37 @@ stats["pairMd5"] = pairMd5[:10] stats["pairFname"] = selectFname # if either the protein sequences, the transcripts, their alignment, the # protein/transcript mapping or the alignment threshold changes -> create a new PSL # mapping file. minAli belongs in here: it decides which alignments survive, so a # mapping built at one threshold must not be silently reused at another. That is # exactly what happened after the MINALI global leaked across threads. refs #38300 fullMd5 = listMd5(allMd5s)[:10] mapFname = join(mapDir, "%(geneTable)s_%(fullMd5)s.psl" % locals()) if isfile(mapFname) and not doForce: logging.info("%s already exists, not rebuilding the protein -> genome mapping PSL" % mapFname) assert(os.path.getsize(mapFname)!=0) + # The lift info describes the mapping, so it has to exist wherever the mapping + # does. A run that was interrupted between building the PSL and writing this file + # leaves the PSL behind on its own, and the next run then reuses the PSL and falls + # over copying the file that was never written. Write it here too. refs #38300 + if not isfile(mapDescFname): + logging.info("%s is missing next to the reused mapping, writing it" % mapDescFname) + writeMapDesc(stats, db, geneTable, mapDescFname) # careful: if you modify this, also modify the other return statement below return mapFname, geneTable, fullMd5, protMapSource, protToTrans, False stats["protMd5"] = protMd5[:10] 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=="miniprot": logging.error("Could not find any gene table for %s, using miniprot to map proteins" % db) miniprotProteins(protFa, db, mapFname, stats) else: stats["minAli"] = minAli @@ -2489,30 +2496,36 @@ " This is the main function that runs a single uniprot update for a list of DBs " uprotDir = options.uniprotDir tabDir = options.tabDir mapDir = options.mapDir bigBedDir = options.bigBedDir faDir = options.faDir doTrembl = not options.skipTrembl if not options.skipParse: checkParserDeps() if not options.skipDownload and not options.skipParse and not options.onlyDbs: # download NCBI -> refseq tables downloadAndSplitNcbi(taxIdDbs) downloadUniprot(uprotDir) + else: + # the per-taxon NCBI tables are keyed on the plan, not on the UniProt release, so a + # plan with taxa that have never been run needs them even when the download is + # skipped. This is the GenArk case: 104 tables existed for the dbDb plan of 109 + # taxa, against 664 taxa here. + splitNcbiIfNeeded(taxIdDbs) downVersion = readRelStringUniprot(uprotDir) # No version file means nothing has been parsed into this directory yet, so parse. # Do not assume it is there: it is absent on a first run, and removing it is the way to # force a reparse when the release has not changed but the taxon list has. versionFname = join(tabDir, "version.txt") tabVersion = open(versionFname).read() if isfile(versionFname) else None doParse = True if downVersion==tabVersion: logging.info("Not converting to tab again, found same version is in tab directory") doParse = False if options.skipParse: logging.info("Not converting to tab, was switched off by option") @@ -3074,41 +3087,68 @@ for row in iterTsvRows(gzip.open(fname, "rt")): taxId = row.tax_id if not int(taxId) in uniprotTaxIds: continue geneId = row.GeneID transId = row.RNA_nucleotide_accession_version if transId!="-": geneToTrans[geneId].add(transId) if lastTax!=taxId and lastTax is not None: writeGeneTsv(geneToTrans, lastTax, outDir) geneToTrans = defaultdict(set) lastTax = taxId writeGeneTsv(geneToTrans, lastTax, outDir) -def downloadAndSplitNcbi(taxIdDbs): - " download and split the NCBI genes file with a mapping NCBI gene -> RefSeq transcripts " - logging.info("Downloading NCBI gene2refseq file gene2refseq.tsv") +ncbiGeneFname = join(NCBIDIR, "gene2refseq.gz") + +def downloadNcbi(): + " download the NCBI gene -> RefSeq transcript file " + logging.info("Downloading NCBI gene2refseq") # --no-progress-meter, not -s: it drops the progress bar, which otherwise writes a few # hundred lines of percentages into lastRun.log, but still shows errors. Plain -s would # hide the errors too, which is why they were dropped here in the first place. - # --fail so an HTTP error is reported as one, instead of piping NCBI's error page into - # zcat and failing with a confusing "not in gzip format" - cmd = "curl --no-progress-meter --show-error --fail https://ftp.ncbi.nlm.nih.gov/gene/DATA/gene2refseq.gz | zcat > ncbi/gene2refseq.tsv" - run(cmd) - splitGeneRefseq("ncbi/gene2refseq.gz", NCBIDIR, taxIdDbs) + # --fail so an HTTP error is reported as one, instead of writing NCBI's error page to + # the file and failing later with a confusing "not in gzip format". + # Keep it gzipped: this used to unpack to gene2refseq.tsv, 16 GB, while the split below + # went on reading gene2refseq.gz, so the download and the split named different files + # and the per-taxon tables were built from whichever .gz happened to be lying around. + tmpName = ncbiGeneFname + ".tmp" + run("curl --no-progress-meter --show-error --fail " + "https://ftp.ncbi.nlm.nih.gov/gene/DATA/gene2refseq.gz > %s" % tmpName) + os.rename(tmpName, ncbiGeneFname) + +def splitNcbiIfNeeded(taxIdDbs): + """ make sure there is a per-taxon NCBI gene -> RefSeq table for every taxon we will + map through RefSeq models. Splitting is cheap next to downloading, and the tables are + per taxon, so a plan with new taxa needs it even when the download is skipped. + """ + missing = [t for t in taxIdDbs if not isfile(join(NCBIDIR, "%s.tsv" % t))] + if not missing: + return + if not isfile(ncbiGeneFname): + errAbort("%s is missing and the download was skipped, so the per-taxon NCBI tables " + "for %d taxa cannot be built. Run without --skipDownload." % + (ncbiGeneFname, len(missing))) + logging.info("%d of %d taxa have no NCBI gene table yet, splitting %s" % \ + (len(missing), len(taxIdDbs), ncbiGeneFname)) + splitGeneRefseq(ncbiGeneFname, NCBIDIR, taxIdDbs) + +def downloadAndSplitNcbi(taxIdDbs): + " download and split the NCBI genes file with a mapping NCBI gene -> RefSeq transcripts " + downloadNcbi() + splitGeneRefseq(ncbiGeneFname, NCBIDIR, taxIdDbs) def delFlag(): global flagFname if isfile(flagFname): os.remove(flagFname) def main(): global flagFname args, options = parseArgs() onlyDbs = None if options.onlyDbs: onlyDbs = set(options.onlyDbs.split(","))