af613a331e6839c6513c3e366abcb67af0fe8386 max Wed Sep 9 06:47:14 2026 -0700 UniProt otto: get the monthly update running again and make a stalled run visible The monthly UniProt job had produced nothing since January 2025. The tracks served release 2024_06 while the download sitting on disk was at 2026_02, on every assembly the job builds. Cause: uniprotToTab appended a personal conda site-packages directory to sys.path, and doUpdate.sh sourced a virtualenv, both built for python 3.6. A venv's python is only a symlink to the system one, so when hgwdev moved to python 3.9 the compiled lxml in there stopped loading and every run died at the parse step. Removed both. The system python3 has lxml from python3-lxml and the two are upgraded together, so there is nothing left here to go stale. Verified by parsing real 2026_02 records under python 3.9 with lxml 5.4. Why nobody noticed for nineteen months: - doUpdate.sh read $? after an intervening echo, so it captured the echo's exit code and mailed "Big Uniprot update OK" every month while the job was dying. It now reads the real exit code, says FAILED, prints the tail of the log and exits nonzero. A month with no new UniProt release stays silent, which is the normal otto behaviour, so silence again means "nothing to do". - The logs were overwritten on every run, so a failure left no trace on disk. doUpdate.sh now appends one line per run to runLog.txt, which is never truncated, and keeps a failing log as lastFail.log. - version.txt in each bigBed directory was rewritten on every run even when the release string was identical. That is the file the trackDb dataVersion setting shows, and its date is what people check to decide whether a pipeline is still alive, so a stalled track could look freshly updated. It is now written only when the release actually changes. Also, so this cannot come back: - doUniprot checks that uniprotToTab can start before the download, instead of finding out 35 minutes later. - pylint on hgwdev is itself pinned to pythons that no longer exist, so "make install" aborted on its first line and could not be used. Replaced with a syntax check that needs nothing but python3; pylint stays best-effort. - uniprotToTab, pslProtCnv, trackDb.template.txt and README.txt ran from /hive/data/outside/otto/uniprot without being in the makefile's copy list. The tree copy of uniprotToTab was still python 2 from 2021. All are now listed and in sync, and "make diff" reports drift. - Brought the two live-only fixes into the tree: mkdir -p in makeUniProtPsl.sh and the pslMap -inType/-mapType flags. refs #38300 diff --git src/utils/uniprotToTab src/utils/uniprotToTab index 47d06f742d5..88c20757797 100755 --- src/utils/uniprotToTab +++ src/utils/uniprotToTab @@ -1,45 +1,56 @@ -#!/usr/bin/env python +#!/usr/bin/env python3 # load default python packages import logging, optparse, sys, glob, gzip, collections, copy, gzip, os, doctest, re from os.path import * from collections import defaultdict +# lxml is not part of the python standard library. On hgwdev it comes from the +# system package python3-lxml, which is upgraded together with /usr/bin/python3. +# Do not put a private site-packages directory on sys.path here: this file used to +# append one from a personal conda environment built for python 3.6, and when the +# system python moved to 3.9 the compiled lxml in it stopped loading. Every monthly +# UniProt update then died at this line and the tracks sat unchanged for 19 months +# before anyone noticed (see redmine #38300). try: - from lxml import etree # if this fails, comment out this line and uncomment the next one. Or do the 'pip install' below. + from lxml import etree # if this fails, comment out this line and uncomment the next one #import xml.etree.cElementTree as etree # if using this line, search for cElementTree in this file and comment out the other part -except: - raise Exception("lxml library not found. install elementtree with 'sudo apt-get install libxml2-dev libxslt-dev python-dev; pip install lxml' or just 'apt-get install python-lxml'. Or read the source code to get rid of the dependency.") +except ImportError as ex: + raise Exception("Cannot import lxml.etree with %s: %s. On hgwdev this module comes " + "from the system package python3-lxml. Check with: python3 -c 'import lxml.etree'. " + "If the system python does not have it, install it with 'pip install --user lxml' or " + "build a private environment, see the README in the otto/uniprot directory." % + (sys.executable, ex)) debugMode=False # --- FASTA FILES --- class FastaReader: """ a class to parse a fasta file Example: fr = FastaReader(filename) for (id, seq) in fr.parse(): print id,seq """ def __init__(self, fname): if hasattr(fname, 'read'): self.f = fname elif fname=="stdin": self.f=sys.stdin elif fname.endswith(".gz"): - self.f=gzip.open(fname) + self.f=gzip.open(fname, "rt") else: self.f=open(fname) self.lastId=None def parse(self): """ Generator: returns sequences as tuple (id, sequence) """ lines = [] for line in self.f: if line.startswith("\n") or line.startswith("#"): continue elif not line.startswith(">"): lines.append(line.replace(" ","").strip()) continue else: @@ -60,59 +71,59 @@ yield faseq else: yield (None, None) def parseFastaAsDict(fname, inDict=None): if inDict==None: inDict = {} fname2 = fname.replace(".gz","") if isfile(fname2): logging.warn("Preferring unzipped file %s" % fname2) fname = fname2 fr = FastaReader(fname) for (id, seq) in fr.parse(): if id in inDict: - print inDict - print inDict[id] + print(inDict) + print(inDict[id]) raise Exception("%s already seen before" % id) inDict[id]=seq return inDict class ProgressMeter: """ prints a message "x%" every stepCount/taskCount calls of taskCompleted() """ def __init__(self, taskCount, stepCount=20, quiet=False): self.taskCount=taskCount self.stepCount=stepCount self.tasksPerMsg = taskCount/stepCount self.i=0 self.quiet = quiet #print "".join(9*["."]) def taskCompleted(self, count=1): if self.quiet and self.taskCount<=5: return #logging.debug("task completed called, i=%d, tasksPerMsg=%d" % (self.i, self.tasksPerMsg)) if self.tasksPerMsg!=0 and self.i % self.tasksPerMsg == 0: donePercent = (self.i*100) / self.taskCount #print "".join(5*[chr(8)]), sys.stderr.write("%.2d%% " % donePercent) sys.stderr.flush() self.i += count if self.i==self.taskCount: - print "" + print("") def setupLogging(progName, options, parser=None, logFileName=None, \ debug=False, fileLevel=logging.DEBUG, minimumLog=False, fileMode="w"): """ direct logging to a file and also to stdout, depending on options (debug, verbose, jobId, etc) """ assert(progName!=None) global debugMode stdoutLevel=logging.INFO if options==None: stdoutLevel=logging.DEBUG elif options.debug or debug: stdoutLevel=logging.DEBUG debugMode = True @@ -188,60 +199,60 @@ "non-terminal residue" : "nonTerm" } # main record info entryHeaders = ["dataset", "acc", "mainIsoAcc", "orgName", "orgCommon", "taxonId", "name", "accList", \ "protFullNames", "protShortNames", "protAltFullNames", "protAltShortNames", \ "geneName", "geneSynonyms", "isoNames", \ "geneOrdLocus", "geneOrf", \ "hgncSym", "hgncId", "refSeq", "refSeqProt", "entrezGene", "ensemblGene", "ensemblProt", "ensemblTrans", \ "kegg", "emblMrna", "emblMrnaProt", "emblDna", "emblDnaProt", \ "pdb", "ec", \ "uniGene", "omimGene", "omimPhenotype", "subCellLoc", "functionText", "isoIds"] EntryRec = collections.namedtuple("uprec", entryHeaders) # all annotations get parsed into this format -annotHeaders = ["acc", "mainIsoAcc", "varId", "featType", "shortFeatType", "begin", "end", "origAa", "mutAa", "dbSnpId", "disRelated", "disease", "disCode", "pmid", "comment"] -AnnotRec = collections.namedtuple("mutrec", annotHeaders) +annotHeaders = ["acc", "mainIsoAcc", "varId", "featType", "shortFeatType", "begin", "end", "origAa", "mutAa", "dbSnpId", "disRelated", "disease", "disCode", "pmid", "longName", "shortName", "syns", "subCellLoc","comment"] +AnnotRec = collections.namedtuple("annotRec", annotHeaders) # references from record refHeaders = ["name", "citType", "year", "journal", "vol", "page", \ "title", "authors", "doi", "pmid", "scopeList"] RefRec = collections.namedtuple("refRec", refHeaders) -emptyRef = dict(zip(refHeaders, len(refHeaders)*[""])) +emptyRef = dict(list(zip(refHeaders, len(refHeaders)*[""]))) def strip_namespace_inplace(etree, namespace=None,remove_from_attr=True): """ Takes a parsed ET structure and does an in-place removal of all namespaces, or removes a specific namespacem (by its URL). Can make node searches simpler in structures with unpredictable namespaces and in content given to be non-mixed. By default does so for node names as well as attribute names. (doesn't remove the namespace definitions, but apparently ElementTree serialization omits any that are unused) Note that for attributes that are unique only because of namespace, this may attributes to be overwritten. For example: <e p:at="bar" at="quu"> would become: <e at="bar"> I don't think I've seen any XML where this matters, though. """ if namespace==None: # all namespaces for elem in etree.getiterator(): tagname = elem.tag - if not isinstance(elem.tag, basestring): + if not isinstance(elem.tag, str): continue if tagname[0]=='{': elem.tag = tagname[ tagname.index('}',1)+1:] if remove_from_attr: to_delete=[] to_set={} for attr_name in elem.attrib: if attr_name[0]=='{': old_val = elem.attrib[attr_name] to_delete.append(attr_name) attr_name = attr_name[attr_name.index('}',1)+1:] to_set[attr_name] = old_val for key in to_delete: elem.attrib.pop(key) @@ -269,45 +280,57 @@ def parseDiseases(fname): " parse the file humanDiseases.txt from uniprot to resolve disease IDs to disease names " logging.info("Parsing %s" % fname) dis = {} for line in open(fname).read().splitlines(): if line.startswith("ID"): name = line[5:].strip(".") if line.startswith("AR"): code = line[5:].strip(".") dis[code]=name logging.info("read %d disease code -> disease name mappings" % len(dis)) return dis -def findSaveList(el, path, dataDict, key, attribKey=None, attribVal=None, useAttrib=None, subSubEl=None): +def findSaveList(el, path, dataDict, key, attribKey=None, attribVal=None, useAttrib=None, subSubEl=None, molAnnots=None): """ find all text of subelemets matching path with given optionally attrib and save into dataDict with key You can specify a subSubEl of the element to get the text from. + If you specify attribKey and attribVal, will only process elements with given key and given value. """ l = [] for se in el.findall(path): if attribKey!=None and se.attrib.get(attribKey, None)!=attribVal: continue if useAttrib: val = se.attrib[useAttrib] else: if subSubEl: val = se.find(subSubEl).text else: val = se.text + + # comments sometimes refer to a molecule, save these for later + molEl = se.find("molecule") + if molEl is not None: + molName = molEl.text + molAnnots[molName]["comment"] = val + val = "Molecule '"+molEl.text+"': "+val + + if attribVal=="ORF" and not val.lower().startswith("orf"): + val = "ORF"+val + l.append(val) s = "|".join(l) dataDict[key] = s def openOutTabFile(subDir, outName, headers): " create outdir and open outfile, write headers " #subDir = join(outDir, outSubDir) if not isdir(subDir): logging.info("Creating dir %s" % subDir) os.makedirs(subDir) outPath = join(subDir, outName) logging.debug("Writing output to %s" % outPath) ofh = open(outPath, "w") ofh.write("\t".join(headers)+"\n") return ofh @@ -516,31 +539,31 @@ val = propEl.attrib["value"] if val=="mRNA": # add now dbRefs["emblMrna"].append(emblId) dbRefs["emblMrnaProt"].append(emblProtId) else: dbRefs["emblDna"].append(emblId) dbRefs["emblDnaProt"].append(emblProtId) continue # don't add any id else: id = dbRefEl.attrib["id"] if id!=None: dbRefs[propDb].add(id) result = {} - for db, valList in dbRefs.iteritems(): + for db, valList in dbRefs.items(): result[db] = "|".join(valList) logging.debug("dbRefs: %s" % result) return result def splitAndResolve(disName, disCodes, splitWord): " split and split word, try to resolve via disCodes and rejoin again " subDises = disName.split(splitWord) newDises = [] for subDis in subDises: subDis = subDis.strip() if subDis in disCodes: newDises.append(disCodes[subDis]) else: newDises.append(subDis) @@ -590,31 +613,31 @@ disLongName = disCode # find snpId snpId = "" for m in re.finditer("dbSNP:(rs[0-9]+)", text): if m!=None: #assert(snpId=="") snpId = m.group(1) logging.debug("Disease: %s, snpId: %s" % (disLongName, snpId)) return disCode, disLongName, snpId, "; ".join(comments) ignoredTypes = collections.Counter() -def parseFeatures(entryEl, disRefs, defaultDisCodes, disToName, evidPmids, mainIsoAcc): +def parseFeatures(entryEl, disRefs, defaultDisCodes, disToName, evidPmids, mainIsoAcc, compAnnots): " go over features and yield annotation records " acc = entryEl.find("accession").text mutations = [] for featEl in entryEl.findall("feature"): featType = featEl.attrib["type"] if featType not in featTypes: ignoredTypes[featType] += 1 continue if featType in ["sequence variant"]: isVariant = True else: isVariant = False @@ -648,30 +671,46 @@ beginEl = featEl.find("location/begin") begin = beginEl.attrib.get("position", None) if begin==None: logging.debug("Unknown start, skipping a feature") continue endEl = featEl.find("location/end") end = endEl.attrib.get("position", None) if end==None: logging.debug("Unknown end, skipping a feature") continue end = str(int(end)+1) # UniProt is 1-based, open-end desc = featEl.attrib.get("description", None) if desc==None: desc = "" + + # for polypeptide chains. Their annotations were in the main protein entry, spread over + # the comments, cell loc and "components" sections + longName = "" + shortName = "" + cellLoc = "" + syns = "" + if desc in compAnnots: + compAnnot = compAnnots[desc] + longName = desc + desc = compAnnot.get("comment", "") + + shortName = compAnnot.get("shortName", "") + cellLoc = compAnnot.get("subCellLocs", "") + syns = compAnnot.get("syns", "") + if "sulfinic" in desc: shortFeatType = "sulfo" descWords = desc.split() if len(descWords)>0: desc1 = descWords[0].lower() if "phos" in desc1: shortFeatType = "phos" elif "acetyl" in desc1: shortFeatType = "acetyl" elif "methyl" in desc1: shortFeatType = "methyl" elif "lipo" in desc1: shortFeatType = "lipo" elif "hydroxy" in desc1: @@ -702,90 +741,184 @@ diseaseRelated = "noEvidence" for evidId in evidList: if evidId in disRefs: diseaseRelated="disRelated" else: diseaseRelated="notDisRelated" logging.debug("evidence is not a disease evidence or blacklisted, check description") pmids = evidPmids.get(evidId, []) assert(len(pmids)<=1) if len(pmids)>0: pmid = list(pmids)[0] annotPmids.append(pmid) - annot = AnnotRec(acc, mainIsoAcc, varId, featType, shortFeatType, begin, end, orig, variant, snpId, diseaseRelated, disName, disCode, ",".join(annotPmids), comments) + annot = AnnotRec(acc, mainIsoAcc, varId, featType, shortFeatType, begin, end, orig, variant, snpId, diseaseRelated, disName, disCode, ",".join(annotPmids), longName, shortName, syns, cellLoc, comments) logging.debug("Accepted annotation: %s" % str(annot)) yield annot def parseEvidence(entryEl): " return a dict with evidCode -> PMID " result = {} for evidEl in entryEl.findall("evidence"): evidCode = evidEl.attrib["key"] for dbRefEl in evidEl.findall("source/dbReference"): dbType = dbRefEl.attrib["type"] if dbType=="PubMed": pmid = dbRefEl.attrib["id"] result.setdefault(evidCode, []) result[evidCode].append(pmid) return result -def parseAnnotations(entryEl, mainIsoAcc, disToName): - " return MutRecs with disease associated variants " +def parseAnnotations(entryEl, mainIsoAcc, disToName, compAnnots): + " return features with features located on protein sequence " # parse the general record comment about diseases disRefs, allDiseaseCodes = parseDiseaseComment(entryEl, disToName) acc = entryEl.find("accession").text logging.debug("Diseases in %s" % acc) evidPmids = parseEvidence(entryEl) - annotRecs = list(parseFeatures(entryEl, disRefs, allDiseaseCodes, disToName, evidPmids, mainIsoAcc)) + annotRecs = list(parseFeatures(entryEl, disRefs, allDiseaseCodes, disToName, evidPmids, mainIsoAcc, compAnnots)) return annotRecs def parseRecInfo(entryEl, entry, isoSeqs): """parse uniprot general record info into entry dict use isoform sequences from isoSeqs only process certain taxonIds """ dataset = entryEl.attrib["dataset"] entry["dataset"] = dataset findSaveList(entryEl, "name", entry, "name") findSaveList(entryEl, "accession", entry, "accList") acc = entry["accList"].split("|")[0] entry["acc"] = acc logging.debug("Parsing rec info for acc %s" % acc) + compAnnots = defaultdict(dict) + + # human-readable gene names and syns findSaveList(entryEl, "protein/recommendedName/fullName", entry, "protFullNames") findSaveList(entryEl, "protein/recommendedName/shortName", entry, "protShortNames") findSaveList(entryEl, "protein/alternativeName/fullName", entry, "protAltFullNames") findSaveList(entryEl, "protein/alternativeName/shortName", entry, "protAltShortNames") findSaveList(entryEl, "gene/name", entry, "geneName", attribKey="type", attribVal="primary") findSaveList(entryEl, "gene/name", entry, "geneSynonyms", attribKey="type", attribVal="synonym") findSaveList(entryEl, "gene/name", entry, "geneOrdLocus", attribKey="type", attribVal="ordered locus") findSaveList(entryEl, "gene/name", entry, "geneOrf", attribKey="type", attribVal="ORF") + + # the rest findSaveList(entryEl, "organism/name", entry, "orgName", attribKey="type", attribVal="scientific") findSaveList(entryEl, "organism/name", entry, "orgCommon", attribKey="type", attribVal="common") findSaveList(entryEl, "organism/dbReference", entry, "taxonId", useAttrib="id") findSaveList(entryEl, "comment/isoform/id", entry, "isoIds") findSaveList(entryEl, "comment/isoform/name", entry, "isoNames") findSaveList(entryEl, "comment/subcellularLocation/location", entry, "subCellLoc") - findSaveList(entryEl, "comment", entry, "functionText", attribKey="type", attribVal="function", subSubEl="text") + findSaveList(entryEl, "comment", entry, "functionText", attribKey="type", attribVal="function", subSubEl="text", \ + molAnnots=compAnnots) + + # when we write out the chains later, some may only be copies of the main protein. So we're keeping a mapping + # of all synonyms now later for these features + synFields= ["protShortNames", "protAltShortNames", "geneOrf", "geneSynonyms", "protAltFullNames", "isoNames", "geneOrdLocus"] + syns = [] + for synField in synFields: + fieldSyns = entry[synField].split("|") + fieldSyns.sort(key=len) # prefer short synonyms + syns.extend(fieldSyns) + # remove empty ones + syns = [x.strip() for x in syns] + syns = [x for x in syns if x!=""] + + # the full protein is sometimes referenced in the record like a peptide + # later so create a faked protein component annotation for the full protein + protFullName = entry["protFullNames"].split("|")[0] + compAnnots[protFullName]["syns"] = "; ".join(syns) + if len(syns)>0: + compAnnots[protFullName]["shortName"] = syns[0] + + # protein components have their short names not stored in the features but only in the protein record + # so we make a translation table here and use it later when writing the features + for compEl in entryEl.findall("protein/component"): + recNameEl = compEl.find("recommendedName") + if recNameEl==None: + assert(False) # component without a name? + + # <fullName>2'-O-methyltransferase</fullName> + # <ecNumber>2.1.1.-</ecNumber> + # </recommendedName>k + #<alternativeName> + #<fullName>nsp16</fullName> + #</alternativeName> </component> + #</protein> + #<gene> + #<name type="primary">rep</name> + #<name type="ORF">1a-1b</name> + #</gene> + + shortName = "" + fullName = "" + fullNameEl = recNameEl.find("fullName") + shortNameEl = recNameEl.find("shortName") + if fullNameEl!=None: + fullName = fullNameEl.text + compAnnots[fullName]["fullName"] = fullName + if shortNameEl!=None: + shortName = shortNameEl.text + compAnnots[fullName]["shortName"] = shortName + + compSyns = [] + for altEl in compEl.findall("alternativeName"): + fullNameEl = altEl.find("fullName") + if fullNameEl is not None: + compSyns.append(fullNameEl.text) + shortNameEl = altEl.find("shortName") + if shortNameEl is not None: + compSyns.append(shortNameEl.text) + + # if we have no shortName, try to take the gene name + # -> not a good idea, as otherwise the components all have the same name + #if shortName=="" and entry["geneName"]!="": + #compAnnots[fullName]["shortName"] = entry["geneName"] + + # if we have no component shortName, take the shortest alternative name + if shortName=="" and len(compSyns)>0: + shortName = list(sorted(compSyns, key=len))[0] + compAnnots[fullName]["shortName"] = shortName + + if len(compSyns)!=0: + if ("compSyns" in compAnnots[fullName]) and fullName!=protFullName: + print("Error: multiple alternative names with synonyms?") + print(compAnnots) + print(compSyns) + print(entry) + assert(False) + compAnnots[fullName]["compSyns"] = "|".join(compSyns) + + # same for subcell localization + for commEl in entryEl.findall("comment"): + if commEl.attrib.get("type")=="subcellular location": + if commEl.find("molecule") is not None: + locs = [] + molName = commEl.find("molecule").text + for locEl in commEl.findall("subcellularLocation"): + locs.append(locEl.text) + compAnnots[molName]["cellLocs"] = "|".join(locs) + mainSeq = entryEl.find("sequence").text #entry["mainSeq"] = mainSeq dbRefs = parseDbRefs(entryEl) isoIds, isoNames, mainIsoId = parseIsoforms(entryEl, acc) entry["mainIsoAcc"] = mainIsoId seqs = [] seqs.append( (mainIsoId+" isRefOf "+acc, mainSeq) ) for isoId in isoIds: if isoId not in isoSeqs: @@ -805,31 +938,31 @@ entry["uniGene"] = dbRefs.get("UniGene", "") entry["omimGene"] = dbRefs.get("omimGene", "") entry["omimPhenotype"] = dbRefs.get("omimPhenotype", "") entry["emblMrna"] = dbRefs.get("emblMrna", "") # mrnas entry["emblMrnaProt"] = dbRefs.get("emblMrnaProt", "") # the protein accessions for mrnas entry["emblDna"] = dbRefs.get("EmblDna", "") # anything not an mrna entry["emblDnaProt"] = dbRefs.get("EmblDnaProt", "") # protein accessions for non-mrnas entry["pdb"] = dbRefs.get("PDB", "") entry["ec"] = dbRefs.get("EC", "") entry["isoIds"]="|".join(isoIds) #entry["isoSeqs"]="|".join(seqs) entry["isoNames"]="|".join(isoNames) entryRow = EntryRec(**entry) - return entryRow, seqs + return entryRow, seqs, compAnnots def parseRefInfo(entryEl, recName): for refEl in entryEl.findall("reference"): ref = copy.copy(emptyRef) ref["name"] = recName citEl = refEl.find("citation") ref["citType"] = citEl.attrib["type"] year = citEl.attrib.get("date", "") ref["year"] = year.split("-")[0] ref["journal"] = citEl.attrib.get("name", "") if ref["journal"]=="": ref["journal"] = citEl.attrib.get("db", "") # for submissions ref["vol"] = citEl.attrib.get("volume", "") ref["page"] = citEl.attrib.get("first", "") for titleEl in citEl.findall("title"): @@ -852,60 +985,60 @@ refRow = RefRec(**ref) yield refRow def readIsoforms(inDir, db): " return all isoform sequences as dict isoName (eg. P48347-2) -> sequence " if db=="swissprot": isoFname = join(inDir, "uniprot_sprot_varsplic.fasta.gz") elif db=="trembl": isoFname = join(inDir, "uniprot_trembl.fasta.gz") else: assert(False) logging.info("reading isoform sequences from %s (or non-gz version)" % isoFname) isoSeqs = parseFastaAsDict(isoFname) result = {} - for id, seq in isoSeqs.iteritems(): + for id, seq in isoSeqs.items(): idParts = id.split("|") isoName = idParts[1] result[isoName] = seq logging.info("Found %d isoform sequences" % len(result)) return result def writeFaSeqs(faFiles, taxonId, seqs): """ write main sequence to faFile with the right taxonId base sequence always has accession as ID """ #seqIds = entry.isoIds.split("|") #if allVariants: if "all" in faFiles: ofh = faFiles["all"] else: ofh = faFiles[taxonId] for seqId, seq in seqs: ofh.write(">%s\n%s\n" % (seqId.strip(), seq.strip())) def openFaFiles(taxonIds, outDir, outPrefix): faFiles = {} if taxonIds == None: taxonIds = ["all"] for taxonId in taxonIds: taxonId = str(taxonId) faFname = join(outDir, outPrefix+"."+taxonId+".fa.gz") - faFiles[int(taxonId)] = gzip.open(faFname, "w") + faFiles[int(taxonId)] = gzip.open(faFname, "wt") logging.debug("Writing fasta seqs for taxon %s to %s" % (taxonId, faFname)) return faFiles def stupidXmlFilter(xmlFile, taxonIds): " return only the uniprot XML lines that refer to one of the taxon Ids " lines = [] taxonOk = False for line in xmlFile: if line.startswith("<entry"): lines = [] lines.append(line) taxonOk = True # in case of doubt, it's True # parse line like <dbReference id="654924" type="NCBI Taxonomy"/> elif line.startswith(' <dbReference id="') and line.endswith('type="NCBI Taxonomy"/>\n'): recTax = int(line.split('"')[1]) @@ -917,147 +1050,140 @@ if taxonOk: yield lines else: yield None lines = None else: if taxonOk and lines is not None: lines.append(line) def parseUniprot(db, inDir, outDir, taxonIds): " parse uniprot, write records and refs to outdir " if options.parse: fname = options.parse logging.info("Debug parse of %s" % fname) - xmlFile = open(fname) + xmlFile = open(inDir) isoSeqs, recCount = {}, 1 - outDir = "." - outPrefix = "temp" + outDir = outDir + outPrefix = "swissprot" disToName = {} else: isoSeqs = readIsoforms(inDir, db) if db=="swissprot": xmlBase = "uniprot_sprot.xml.gz" outPrefix = "swissprot" recCount = 600000 elif db=="trembl": xmlBase = "uniprot_trembl.xml.gz" outPrefix = "trembl" - recCount = 600000*100 + recCount = 600000*700 else: raise Exception("unknown db") xmlBase2 = xmlBase.replace(".gz", "") if isfile(xmlBase2): logging.debug("Using non-gzipped file %s" % xmlBase2) xmlFile = open(join(inDir, xmlBase2)) else: - xmlFile = gzip.open(join(inDir, xmlBase)) + xmlFile = gzip.open(join(inDir, xmlBase), "rt") logging.info("Parsing main XML file %s" % xmlFile.name) disToName = parseDiseases(join(inDir, "docs", "humdisease.txt")) faFiles = openFaFiles(taxonIds, outDir, outPrefix) logging.debug("Only extracting taxon IDs %s" % str(taxonIds)) # create a dict taxonId -> output file handles for record info, pmid reference info and annotation info outFhs = {} for taxId in taxonIds: entryOf = openOutTabFile(outDir, "%s.%s.tab" % (outPrefix, taxId), entryHeaders) refOf = openOutTabFile(outDir, "%s.%s.refs.tab" % (outPrefix, taxId), refHeaders) annotOf = openOutTabFile(outDir, "%s.%s.annots.tab" % (outPrefix, taxId), annotHeaders) outFhs[taxId] = (entryOf, refOf, annotOf) - emptyEntry = dict(zip(entryHeaders, len(entryHeaders)*[""])) + emptyEntry = dict(list(zip(entryHeaders, len(entryHeaders)*[""]))) pm = ProgressMeter(recCount) for recLines in stupidXmlFilter(xmlFile, taxonIds): # the original solution below did a partial parse, but required 300GB of # RAM #for _, entryEl in etree.iterparse(xmlFile): #if entryEl.tag!="{http://uniprot.org/uniprot}entry": #continue - #strip_namespace_inplace(entryEl) # die, die stupid namespaces!! pm.taskCompleted() if recLines is None: continue entryEl = etree.fromstring("".join(recLines)) + strip_namespace_inplace(entryEl) # die, die stupid namespaces!! - entryTax = int(entryEl.find("organism/dbReference").attrib["id"]) + dbRefEl = entryEl.find("organism") + dbRefEl = entryEl.find("organism/dbReference") + entryTax = int(dbRefEl.attrib["id"]) if taxonIds==['all']: taxId = "all" else: if entryTax not in taxonIds: logging.debug("taxon ID %d not a target taxon" % entryTax) continue entryOf, refOf, annotOf = outFhs[entryTax] entry = copy.copy(emptyEntry) - entryRow, seqs = parseRecInfo(entryEl, entry, isoSeqs) + entryRow, seqs, compAnnots = parseRecInfo(entryEl, entry, isoSeqs) writeFaSeqs(faFiles, entryTax, seqs) entryOf.write("\t".join(entryRow)+"\n") recName = entryRow.name refRows = list(parseRefInfo(entryEl, recName)) for refRow in refRows: refOf.write("\t".join(refRow)+"\n") - annotRecs = parseAnnotations(entryEl, entryRow.mainIsoAcc, disToName) + annotRecs = parseAnnotations(entryEl, entryRow.mainIsoAcc, disToName, compAnnots) for annotRow in annotRecs: logging.debug("writing row %s" % str(annotRow)) annotOf.write("\t".join(annotRow)+"\n") - # clear some RAM (not all) - # not needed anymore, since I'm using stupidXmlFilter now - # https://stackoverflow.com/questions/12160418/why-is-lxml-etree-iterparse-eating-up-all-my-memory - # remove this if using cElementTree - #entryEl.clear() - #for ancestor in entryEl.xpath('ancestor-or-self::*'): - #while ancestor.getprevious() is not None: - #del ancestor.getparent()[0] - logging.info("Skipped annotation types: %s" % ignoredTypes.most_common()) def main(args, options): if options.test: import doctest doctest.testmod() sys.exit(0) setupLogging("pubParseDb", options) db = args[0] db = "swissprot" if options.trembl: db = "trembl" dbDir = args[0] - taxonIds = args[1] + refDir = args[2] + if taxonIds=="all": taxonIds = ['all'] else: taxonIds=[int(x) for x in taxonIds.split(",")] - refDir = args[2] if not isdir(refDir): logging.info("Making directory %s" % refDir) os.makedirs(refDir) if len(args)!=3: raise Exception("Invalid command line. Show help with -h") parseUniprot(db, dbDir, refDir, taxonIds) # === COMMAND LINE INTERFACE, OPTIONS AND HELP === parser = optparse.OptionParser("""usage: %prog [options] uniprotFtpDir taxonIds outDir - Convert UniProt to tab-sep files taxonIds can be "all" To download uniProt, this command is a good idea: @@ -1072,24 +1198,24 @@ - only gets a limited list of xrefs, but others are easy to add - can parse Trembl (Who came up with the idea of creating a 500GB XML file?) Example: %prog /hive/data/outside/uniProt/current 9606 tab/ If you get no results from this script, your species may be only in Trembl. Use the '--trembl' option to parse UniProt/Trembl instead of UniProt/SwissProt. Most organisms have entries in both databases and you have to run the script twice to get all entries. """) parser.add_option("-d", "--debug", dest="debug", action="store_true", help="show debug messages") parser.add_option("", "--test", dest="test", action="store_true", help="run tests") -parser.add_option("-p", "--parse", dest="parse", action="store", help="parse a single uniprot xml file (debugging)") +parser.add_option("-p", "--parse", dest="parse", action="store_true", help="parse a single uniprot xml file") parser.add_option("", "--trembl", dest="trembl", action="store_true", help="parse trembl. Default is to parse only the swissprot files.") (options, args) = parser.parse_args() if args==[] and not options.test: parser.print_help() exit(1) main(args, options)