953e29496b012cc36ee7933ce94cefea2f66ec38 max Fri Sep 11 09:21:22 2026 -0700 uniprot otto: align proteins with miniprot instead of BLAT or Augustus Where an assembly has no gene models worth mapping through, the pipeline used "blat -q=prot -t=dnax", which is slow enough that hs1 was excluded from the whole job over it. On GenArk assemblies the alternative was Augustus, which is an ab initio prediction, so mapping UniProt through it stacks its errors on top of ours. Aligning the proteins straight to the genome avoids both. miniprot 0.18 built from github into /hive/data/outside/otto/uniprot/bin/miniprot, linked against system libraries only and runnable by otto. miniprotProteins() replaces blatProteinsKeepBest(). miniprot reads fasta rather than 2bit, so the genome is unpacked to a temp file and removed again once the alignment is done. Its GFF3 needs some care before gff3ToGenePred will take it: the ##PAF meta lines and the capitalised Rank/Identity attributes have to go, and alignments are named MP000001 with the UniProt accession hidden in Target=, so the ids are rewritten to the accession. The uniquifying suffix that GFF3 requires is stripped again afterwards, leaving the bare accession as qName, which is what the rest of the pipeline keys on. Augustus is dropped from the GenArk gene source list entirely, so the order is now catGenes, ncbiRefSeq, ncbiGene, then miniprot. The classic "blat" fallback becomes miniprot too. Checked on real data: 300 bonobo UniProt proteins against one 227 Mb chromosome give 138 alignments over 102 distinct proteins, pslCheck reports 138 checked and 0 failed, the temp genome fasta is cleaned up, and the aligner version is recorded in the mapping stats. GRCz12ab and calJac240_pri, the two assemblies that had nothing but Augustus, now resolve to miniprot. refs #38300 diff --git src/hg/utils/otto/uniprot/doUniprot src/hg/utils/otto/uniprot/doUniprot index 290b2efc3ca..aa3ecf69557 100755 --- src/hg/utils/otto/uniprot/doUniprot +++ src/hg/utils/otto/uniprot/doUniprot @@ -999,32 +999,41 @@ # GenArk assemblies keep their gene models as bigBed files in the hub, not as MySQL tables, # so they need their own search. In descending order of preference: # CAT - Comparative Annotation Toolkit, shipped as a contrib collection # (track hprcCatGenes, contrib/hprc2annot/catGenes.bb). Best annotation we # have on the assemblies that carry it. # ncbiRefSeq - RefSeq, the GenArk RefSeq gene track. On most GCF assemblies. # ncbiGene - the annotation the submitter sent to GenBank with the assembly. # augustus - ab initio prediction, present almost everywhere, so it is the last resort. # Each entry is (name, glob relative to the hub directory). Globs, not constructed names, # because the files under bbi/ carry the full asmId with its assembly-name suffix while the # hub directory is named with the short accession. genArkGeneSources = [ ("catGenes", "contrib/*/catGenes.bb"), ("ncbiRefSeq", "bbi/*.ncbiRefSeq.bb"), ("ncbiGene", "bbi/*.ncbiGene.bb"), - ("augustus", "bbi/*.augustus.bb"), ] +# Augustus is deliberately not in that list. It is an ab initio prediction, so mapping +# UniProt through it stacks its errors on top of ours; aligning the proteins straight to +# the genome with miniprot is better and much faster than the old BLAT protein search. +miniprotBin = "/hive/data/outside/otto/uniprot/bin/miniprot" +miniprotThreads = 16 + +def miniprotVersion(): + " version string of the miniprot binary we are using " + proc = subprocess.Popen([miniprotBin, "--version"], stdout=PIPE, encoding="utf8") + return proc.communicate()[0].strip() 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 " @@ -1042,66 +1051,66 @@ 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)) + logging.warning("%s: no gene models found in %s, will align the proteins directly with miniprot" % (db, hubDir)) return None, None def findBestGeneTable(db): " find the best gene table for a given organism and return it " hubDir = genArkHubDir(db) if hubDir is not None: name, _ = findBestGeneBigBed(db, hubDir) - return name if name else "blat" + return name if name else "miniprot" if db=="hg19": tables = ["refGene"] # because in 2021, refSeq maps NM_001129826.3 still to chrX_jh159150_fix, confirmed as a bug by Terence # CSAG is an important gene elif db=="hg38": #tables = [ "wgEncodeGencodeCompV37lift37" ] tables = ["ncbiRefSeq"] else: tables = [ "ncbiRefSeq", "ensGene", "augustusGene"] # refGene has only a few hundred transcripts on most obscure organisms -> refGene only covers the manual NM_ sequences! for table in tables: query = "DESCRIBE %s" % (table) cmd = "hgsql %s -e 'DESCRIBE %s' > /dev/null 2>&1" % (db, table) ret = os.system(cmd) if ret==0: logging.debug("Best gene table for db %s is %s" % (db, table)) if db in ["hg19", "hg38", "mm10", "mm39"]: logging.info("DB is human or mouse. Tolerating multi mappers = transcript IDs that appear multiple times in the genome.") return table if hasNoMultiMappers(db, table): logging.info("%s: Best gene table is %s. No multi mappers." % (db, table)) return table - assert(db not in ["hg19", "hg38", "mm10", "danRer10", "mm39"]) # never use BLAT on our main dbs - return "blat" + assert(db not in ["hg19", "hg38", "mm10", "danRer10", "mm39"]) # never align proteins directly on our main dbs + return "miniprot" def gz_is_empty(fname): # from https://stackoverflow.com/questions/37874936/how-to-check-empty-gzip-file-in-python ''' Test if gzip file fname is empty Return True if the uncompressed data in fname has zero length or if fname itself has zero length Raises OSError if fname has non-zero length and is not a gzip file ''' with gzip.open(fname, 'rb') as f: data = f.read(1) return len(data) == 0 #def parseFaIds(fname): #"parse fasta, return dict seqId -> sequence " #seqIds = set() @@ -1466,45 +1475,94 @@ """ fnames = [ (join(tabDir,"swissprot.%d.tab" % taxId)) ] if doTrembl: fnames.append( (join(tabDir,"trembl.%d.tab" % taxId)) ) ret = {} for fname in fnames: for row in iterTsvRows(open(fname)): for acc in getAllIsoAccs(row): ret[acc] = row return ret -def blatProteinsKeepBest(fullFaFname, db, mapFname, stats): - " blat the proteins directly, only good for small genomes where we don't have any gene models " - myDir = dirname(__file__) - # consider using miniprot here one day? - cmd = "blat -q=prot -t=dnax /gbdb/{db}/{db}.2bit {inf} stdout -noHead | sort -k10 | pslReps stdin stdout /dev/null " \ - "-minAli=0.95 -nohead | {myDir}/pslProtCnv > {out}".format(db=db, inf=fullFaFname, out=mapFname, myDir=myDir) - stats["minAli"]=0.95 - #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) +def twoBitFname(db): + " the assembly sequence, wherever this assembly keeps it " + hubDir = genArkHubDir(db) + if hubDir is not None: + return genArkTwoBit(db, hubDir) + return "/gbdb/{db}/{db}.2bit".format(db=db) + +def miniprotProteins(fullFaFname, db, mapFname, stats): + """ align the UniProt proteins straight to the genome with miniprot and write a PSL. + Used when the assembly has no gene models worth mapping through. This replaced both + "blat -q=prot" and mapping through Augustus: Augustus is ab initio, so going through it + stacks its errors on top of ours, and the BLAT protein search is very slow on a big + genome. + """ + workDir = mapFname+".miniprot.tmp" + if not isdir(workDir): + os.makedirs(workDir) + + # miniprot reads fasta, not 2bit + genomeFa = join(workDir, "genome.fa") + run(["twoBitToFa", twoBitFname(db), genomeFa]) + + gffName = join(workDir, "miniprot.gff") + run("%s -t %d --gff %s %s > %s" % (miniprotBin, miniprotThreads, genomeFa, fullFaFname, gffName)) + os.remove(genomeFa) + + # Name each alignment after the UniProt accession it came from. miniprot calls them + # MP000001 and puts the accession in Target=, and its ##PAF meta lines and capitalised + # attributes make gff3ToGenePred unhappy, so both are cleaned up here. The ID has to + # stay unique for GFF3, hence the .N suffix, which comes off again after the + # conversion so that qName is the bare accession. + namedName = join(workDir, "named.gff") + cmd = """awk -F'\\t' -v OFS='\\t' ' + /^##PAF/ {next} + /^#/ {print; next} + $3=="mRNA" { + id=""; acc="" + if (match($9, /ID=[^;]+/)) id = substr($9, RSTART+3, RLENGTH-3) + if (match($9, /Target=[^ ;]+/)) acc = substr($9, RSTART+7, RLENGTH-7) + if (id!="" && acc!="") { map[id]=acc; $9 = "ID=" acc "." (++n[acc]) } + print; next + } + { if (match($9, /Parent=[^;]+/)) { + p = substr($9, RSTART+7, RLENGTH-7) + if (p in map) $9 = "Parent=" map[p] "." n[map[p]] + } + print + }' %s > %s""" % (gffName, namedName) run(cmd) + gpName = join(workDir, "miniprot.gp") + run(["gff3ToGenePred", "-warnAndContinue", "-maxParseErrors=-1", "-maxConvertErrors=-1", + namedName, gpName]) + run("""awk -F'\\t' -v OFS='\\t' '{sub(/\\.[0-9]+$/, "", $1); print}' %s > %s.acc""" % (gpName, gpName)) + + cdsName = join(workDir, "miniprot.cds") + chromSizes = join(workDir, "chrom.sizes") + run("twoBitInfo %s %s" % (twoBitFname(db), chromSizes)) + run(["genePredToFakePsl", "-chromSize=%s" % chromSizes, "noDb", gpName+".acc", mapFname, cdsName]) + + alnCount = len(open(mapFname).readlines()) + logging.info("%s: miniprot aligned the UniProt proteins to %d genomic locations" % (db, alnCount)) + stats["aligner"] = "miniprot " + miniprotVersion() + 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! # 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)) @@ -1804,31 +1862,31 @@ 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") 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": + elif geneTable=="miniprot": 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: json.dump(stats, mapDescFh, indent=4) def listMd5(arr): " return hex md5 given list of strings " @@ -1869,45 +1927,45 @@ 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() stats["taxId"] = taxId stats["db"] = db geneTable = findBestGeneTable(db) metaMd5 = manyTabMd5(tabFnames) - if geneTable!="blat": + if geneTable!="miniprot": transcriptFa, transcriptPsl = makeTranscriptFiles(db, geneTable, taxId, mapDir) transMd5 = fastaMd5(transcriptFa) pslMd5 = tabMd5(transcriptPsl) else: # this usually only happens on viral genomes or other weird cases that don't have a single gene model track # here we blat the proteins directly onto the genome transMd5, pslMd5 = "directBlat-noTranscripts", "directBlat-noPsl" transcriptFa, transcriptPsl = None, None geneToRefSeqs = None selectFname = None protToTrans = None - if geneTable in ["blat"]: + if geneTable in ["miniprot"]: protMapSource = {"default":"direct"} elif geneTable in ["augustusGene"]: protMapSource = {"default":"best"} else: if geneTable in ["refGene", "ncbiRefSeq"]: geneToRefSeqs = readGeneToRefSeq(taxId) selectFname, stats, protMapSource = buildSelectFile(tabFnames, mapDir, taxId, db, \ geneTable, transcriptFa, geneToRefSeqs, metaMd5, transMd5, stats) protToTrans = parseKeyVal(selectFname) allMd5s = [protMd5, transMd5, pslMd5] if selectFname: pairMd5 = tabMd5(selectFname) allMd5s.append(pairMd5) stats["pairMd5"] = pairMd5[:10] @@ -1920,33 +1978,33 @@ 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) # 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=="blat": + if geneTable=="miniprot": logging.error("Could not find any gene table for %s, using BLAT to map proteins" % db) - blatProteinsKeepBest(protFa, db, mapFname, stats) + miniprotProteins(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), 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)