bb3efd1b30b2fe47ce45c3c336f8982b5c7a7463 max Thu Feb 18 06:32:41 2021 -0800 tiny change to uniprot pipeline, no redmine diff --git src/hg/utils/otto/uniprot/doUniprot src/hg/utils/otto/uniprot/doUniprot index 1c59145..ebdadc6 100755 --- src/hg/utils/otto/uniprot/doUniprot +++ src/hg/utils/otto/uniprot/doUniprot @@ -1,1532 +1,1536 @@ #!/usr/bin/env python3 import urllib.request, urllib.error, urllib.parse, time, os, atexit import datetime, optparse, sys, logging, subprocess, glob, shutil import sys, gzip, logging,re from urllib.parse import urlparse from os.path import * from os import system, makedirs, mkdir, remove, listdir from shutil import move from subprocess import PIPE from collections import defaultdict, namedtuple, Counter # main driver script for uniprot updates # This script has evolved over a long time. Parts of it do not have an optimal structure. # One reason is that the script used to parse only mutations. Even normal annotations # are still put into a structure that resembles mutations. # Globals --- #upUrl = "ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/complete/uniprot_sprot.xml.gz" upUrl = "https://ftp.expasy.org/databases/uniprot/current_release/knowledgebase/complete/uniprot_sprot.xml.gz" # manual mapping from uniProt taxon IDs to UCSC databases, the default mapping is taken from dbDb # these are added in addition, I obtained them in 2018 using some Apache log mining # the list could probably be a lot shorter now manualTaxIdDbs = { 9606 : ["hg19", "hg38"], # human 10090 : ["mm10"], # mouse 10116 : ["rn6"], # rat 8364 : ["xenTro9"], # xenopus 6239 : ["ce11", "ce10"], # c elegans 7227 : ["dm6"], # fruitfly 7955 : ["danRer7", "danRer10"], # zebrafish 4932 : ["sacCer3"], # yeast 9940 : ["oviAri3"], # sheep 9823 : ["susScr3"], # pig 2697049: None # ["wuhCor1"] # do not run on wuhCor1 - wuhCor1 has a special daily otto job in /hive/data/genomes/wuhCor1/bed/uniprot } # when we update the parasol cluster, no need to patch this script clusterFname = "/cluster/bin/scripts/cluster.txt" TREMBLCOLOR="0,150,250" # light blue SWISSPCOLOR="12,12,120" # dark blue asDir = dirname(__file__) #urls = { ##"var": "http://web.expasy.org/cgi-bin/variant_pages/get-sprot-variant.pl?", #"var": "http://www.uniprot.org/uniprot/", #"uniProt" : "http://www.uniprot.org/uniprot/", #"pubmed" : "http://www.ncbi.nlm.nih.gov/pubmed/" #} # 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"]) # the feature type for the 'faked' full sequence features for the unipFullSeq track fullSeqType = "UniProt sequence" # mostly for the "other annotations" subtrack, to show the different types featTypeColors = { "modified residue" : "200,200,0", "glycosylation site" : "0,100,100", "disulfide bond" : "100,100,100", "topological domain" : "100,0,0", "zinc finger region" : "100,100,0" } # some very common disease codes are not in the diseases.txt file disShortNames = { "hepatocellular carcinoma" : "HepatocC", "a hepatocellular carcinoma" : "HepatocC", "a hepatocellular carcinoma sample" : "HepatocC", "hepatocellular carcinomas" : "HepatocC", "ovarian cancer" : "OvC", "colon cancer" : "ColC", "colorectal cancer" : "ColrectC", "hepatoblastoma" : "HepatoBl", "a sporadic cancer": "sporadic", "sporadic cancers": "sporadic", "non-small cell lung cancer cell lines": "NSCLC-CL", "a colon cancer cell line":"ColC-CL" } def getTaxIdDbs(): " return a dict with taxonId -> list of most recent dbs (e.g. for human, it's hg19 and hg38) " query = "select taxId, name from dbDb where active=1 order by orderKey;" rows = runQuery("hgcentral", query, usePublic=True) taxIdToDbs = defaultdict(list) # take only the first database for every taxonId for taxId, dbCode in rows: taxId = int(taxId) if taxId not in taxIdToDbs: taxIdToDbs[taxId].append(dbCode) for taxId, manNames in manualTaxIdDbs.items(): # "None" means: do not do anything for this taxonId if manNames is None: del taxIdToDbs[taxId] else: for manName in manNames: if manName not in taxIdToDbs[taxId]: taxIdToDbs[taxId].append(manName) return taxIdToDbs.items() def parseArgs(): parser = optparse.OptionParser("""usage: %prog [options] [run] - run uniprot update: download UniProt, parse the annotations to .tab files, build liftOver .psls and lift the annotations to bigBed files, one set per genome db""") parser.add_option("-d", "--debug", dest="debug", action="store_true", \ help="show debug messages") parser.add_option("", "--dbs", dest="dbs", action="store", \ help="only run on certain databases, comma-sep list of UCSC database names, e.g. 'hg19,mm10' (debugging)") parser.add_option("-l", "--skipDownload", dest="skipDownload", action="store_true", \ help="Skip the download step (debugging)") parser.add_option("", "--skipTrembl", dest="skipTrembl", action="store_true", \ help="Skip processing Trembl") parser.add_option("-p", "--skipParse", dest="skipParse", action="store_true", \ help="Skip the download and parsing step (debugging)") parser.add_option("", "--onlyMap", dest="onlyMap", action="store_true", \ help="only create the pslMap files (debugging)") parser.add_option("-u", "--uniprotDir", dest="uniprotDir", action="store", \ help="UniProt download data directory, default %default", default="/hive/data/outside/uniProt/current/") parser.add_option("-t", "--tabDir", dest="tabDir", action="store", \ help="directory for Uniprot tab-sep and fasta files, default %default", default="tab") parser.add_option("-f", "--faDir", dest="faDir", action="store", \ help="directory for full fasta files, default %default", default="fasta") parser.add_option("-m", "--mapDir", dest="mapDir", action="store", \ help="directory for pslMap (~liftOver) psl files, default %default", default="map") parser.add_option("-b", "--bigBedDir", dest="bigBedDir", action="store", \ help="directory for bigBed files, one subdirectory per db will be created", default="bigBed") parser.add_option("", "--force", dest="force", action="store_true", \ help="skip the check of differences against the previous version") parser.add_option("", "--onlyLinks", dest="onlyLinks", action="store_true", \ help="only create the /gbdb/ symlinks") parser.add_option("", "--skipLinks", dest="skipLinks", action="store_true", \ help="do not create the /gbdb/ symlinks") parser.add_option("", "--mapQa", dest="mapQa", action="store_true", \ help="output some QA stats for the maps") parser.add_option("", "--db", dest="db", action="store_true", \ help="output the trackDb make command and uniprot <-> UCSC db assignments for debugging trackDb problems and showing which UCSC databases will be processed by the otto job and why") (options, args) = parser.parse_args() if options.db: taxIdDbs = getTaxIdDbs() print(("Current TaxId<->database assignment: %s" % str(taxIdDbs))) allDbs = [] for taxId, dbs in taxIdDbs: allDbs.extend(dbs) print(("To make the uniProt trackDbs for all DBs, run this in kent/src/hg/makeDb/trackDb: make alpha DBS='%s'" % " ".join(allDbs))) taxIdStr = ",".join([str(x) for (x,y) in taxIdDbs]) # just the taxIds themselves print(("to just convert the XML files (debugging?), run this: uniprotToTab %s %s %s" % (options.uniprotDir, taxIdStr, options.tabDir))) sys.exit(1) if len(args)==0 and not options.mapQa: print("To actually run the pipeline, you need to specify the argument 'run'.") parser.print_help() print("To actually run the pipeline, you need to specify the argument 'run'.") + sys.exit(0) if options.debug: consLevel = logging.DEBUG else: consLevel = logging.INFO # '' is the root logger logger = logging.getLogger('') logger.setLevel(logging.DEBUG) #By default, logs all messages # log to console ch = logging.StreamHandler() #StreamHandler logs to console ch.setLevel(consLevel) ch_format = logging.Formatter('%(asctime)s - %(message)s') ch.setFormatter(ch_format) # also always log everything to file, appends by default fname = abspath("doUniprot.log") fh = logging.FileHandler(fname) fh.setLevel(logging.DEBUG) fh_format = logging.Formatter('%(asctime)s - %(lineno)d - %(levelname)-8s - %(message)s') fh.setFormatter(fh_format) logger.addHandler(fh) logger.addHandler(ch) logging.info("Logging debug messages to %s" % fname) return args, options # ------------- def iterTsvRows(inFile,encoding="utf8", fieldSep="\t", isGzip=False): """ parses tab-sep file with headers as field names yields collection.namedtuples strips "#"-prefix from header line """ fh = inFile line1 = fh.readline() line1 = line1.strip("\n").strip("#") headers = line1.split(fieldSep) headers = [re.sub("[^a-zA-Z0-9_]","_", h) for h in headers] Record = namedtuple('tsvRec', headers) for line in fh: line = line.rstrip("\n") fields = line.split(fieldSep) if encoding!=None: fields = [f for f in fields] try: rec = Record(*fields) except Exception as msg: logging.error("Exception occured while parsing line, %s" % msg) logging.error("Filename %s" % fh.name) logging.error("Line was: %s" % line) logging.error("Does number of fields match headers?") logging.error("Headers are: %s" % headers) raise Exception("wrong field count in line %s" % line) # convert fields to correct data type yield rec threeToOne = \ {'Cys': 'C', 'Asp': 'D', 'Ser': 'S', 'Gln': 'Q', 'Lys': 'K', 'Ile': 'I', 'Pro': 'P', 'Thr': 'T', 'Phe': 'F', 'Asn': 'N', 'Gly': 'G', 'His': 'H', 'Leu': 'L', 'Arg': 'R', 'Trp': 'W', 'Ala': 'A', 'Val':'V', 'Glu': 'E', 'Tyr': 'Y', 'Met': 'M', 'Glx' : 'Z', 'Asx' : 'B', 'Pyl' : 'O', "Xaa" : "X", 'Sec': 'U' # very very rare amino acid (=stop codon) } oneToThree = dict([[v,k] for k,v in list(threeToOne.items())]) def aaToLong(seq): " convert amino acid to three letter code " res = [] for aa in seq: longAa = oneToThree.get(aa, aa) if longAa==aa: logging.error("unknown iupac "+ aa) res.append(longAa) return "-".join(res) def shortenDisCode(code): " some disease names are not shortened yet by UniProt. Do this now. " if code=="": return [] newCodes = [] if " and " in code: splitCodes = code.split(" and ") newCodes.append(disShortNames.get(splitCodes[0], splitCodes[0])) newCodes.append(disShortNames.get(splitCodes[1], splitCodes[1])) else: newCodes.append(disShortNames.get(code, code)) return newCodes def htmlLink(urlType, accs): # at the moment, not creating any html links, but this may be needed one day later, when # the links are not as well-formed as they are at the moment #strList = [] #for acc in accs: #strList.append('<a href="%s%s">%s</a>' % (urls[urlType], acc, acc)) #return ", ".join(strList) return ",".join(accs) def makeBedLine(dbName, firstAnnot, bed, isMut, color): """ convert a list of annotation objects and a bed with their positions to a single BED line with extra fields """ # try to create a more human-readable version of the disease codes disCodes = shortenDisCode(firstAnnot.disCode) disName = ",".join(disCodes) if len(disName)>30: disName = disName[:30]+"..." acc = firstAnnot.acc #dbSnpIds = [a.dbSnpId for a in annots] # set the bed name field to a disease, to the mutation or something else if isMut and firstAnnot.origAa=="": bed[3] = "%s-%sdel" % (firstAnnot.begin, firstAnnot.end) else: bed[3] = "%s%s%s" % (firstAnnot.origAa, firstAnnot.begin, firstAnnot.mutAa) if len(disCodes)>0 and not "strain" in disName: bed[3] += " in %s" % disName color = None if firstAnnot.featType=="sequence variant": annoType = "Naturally occurring sequence variant" color = "100,0,0" elif firstAnnot.featType=="mutagenesis site": annoType = "Experimental mutation of amino acids" color = "0,0,100" else: bed[3] = firstAnnot.shortFeatType if firstAnnot.featType in useComment: comment = firstAnnot.comment # trembl doesn't have comments for chains so make one up if comment=="": comment = "%s:%s-%s" % (acc, str(firstAnnot.begin), str(firstAnnot.end)) bed[3] = comment if firstAnnot.featType=="signal peptide": bed[3] = "Signal peptide" if firstAnnot.featType=="lipid moiety-binding region": bed[3] = "Lipidation" if firstAnnot.featType=="transmembrane region": bed[3] = "Transmembrane" if firstAnnot.featType==fullSeqType: bed[3] = firstAnnot.acc if firstAnnot.comment=="Nuclear localization signal": bed[3] = "Nuclear loc" if firstAnnot.comment=="Involved in receptor recognition and/or post-binding events": bed[3] = "Recept recog" if firstAnnot.comment=="Fibronectin type-III": bed[3] = "FibronectIII" # more general rules if firstAnnot.comment.startswith("Zinc finger protein "): bed[3] = firstAnnot.comment.replace("Zinc finger protein ", "ZF-") if firstAnnot.comment.startswith("Necessary for the"): bed[3] = firstAnnot.comment.replace("Necessary for the ", "") if firstAnnot.comment.startswith("Interaction with "): bed[3] = firstAnnot.comment.replace("Interaction with ", "Int:") if firstAnnot.comment.startswith("Involved in the"): bed[3] = firstAnnot.comment.replace("Involved in the ", "") if firstAnnot.comment.startswith("Required for "): bed[3] = firstAnnot.comment.replace("Required for ", "") if firstAnnot.comment.startswith("Cleavage; by host "): bed[3] = "Cleave:"+firstAnnot.comment.split()[-1] if firstAnnot.comment=="Receptor-binding motif; binding to human ACE2": bed[3] = "binds ACE2" # chain annotations if firstAnnot.longName!="": bed[3] = firstAnnot.longName if firstAnnot.shortName!="": bed[3] = firstAnnot.shortName if bed[3] == "Intrinsically disordered": bed[3] = "Disordered" if len(bed[3])>17: bed[3] = bed[3][:14]+"..." annoType = firstAnnot.featType if color is None: bed[8] = featTypeColors.get(annoType, SWISSPCOLOR) else: bed[8] = color if annoType=="strand": annoType="beta strand" # bigGenePred format makes no sense for mutations but a lot for the other annotations # add the bigGenePred fields: name2, cdsStartStat, cdsEndStat, exonFrames, type, geneName, geneName2, geneType if not isMut: bed.append("") # name2 bed.append("cmpl") # cdsStartStat bed.append("cmpl") # cdsEndStat blockCount = int(bed[9]) bed.append(",".join(blockCount*["0"])) #exonFrames bed.append("") # type bed.append("") # geneName bed.append("") # geneName2 bed.append("") # geneType # status field if dbName=="trembl": bed.append("Unreviewed (TrEMBL)") else: bed.append("Manually reviewed (Swiss-Prot)") bed.append(annoType) if isMut: bed.append(firstAnnot.disease) if isMut: # create description of mutation if int(firstAnnot.end)-int(firstAnnot.begin)==1: posStr = "position %s" % firstAnnot.begin else: posStr = "position %s-%s" % (firstAnnot.begin, firstAnnot.end) if firstAnnot.origAa=="": bed.append("%s, removal of amino acids" % (posStr)) else: bed.append("%s, %s changed to %s" % \ (posStr, aaToLong(firstAnnot.origAa), aaToLong(firstAnnot.mutAa))) else: # just position of feature if int(firstAnnot.end)-int(firstAnnot.begin)==1: posStr = "amino acid %s" % str(int(firstAnnot.begin)) # show a 1-based coordinate else: posStr = "amino acids %s-%s" % (firstAnnot.begin, str(int(firstAnnot.end)-1)) posStr += " on protein %s" % firstAnnot.acc bed.append(posStr) if not isMut: bed.append(firstAnnot.longName) bed.append(firstAnnot.syns.replace("|", "; ")) bed.append(firstAnnot.subCellLoc) # construct the comment field, which is also used for the mouse over commentField = firstAnnot.comment # add the disease to the mouse over if firstAnnot.disease!="": commentField = firstAnnot.disease + "; " + commentField if firstAnnot.featType == "sequence conflict": commentField = "%s, %s changed to %s, %s" % \ (posStr, aaToLong(firstAnnot.origAa), aaToLong(firstAnnot.mutAa), commentField) bed.append(commentField) # fields that are only used for mutations if isMut: varId = firstAnnot.varId if varId!="": varId = "%s#%s|%s" % (acc, varId, varId) # Andrew Nightingale requested that we link to UniProt bed.append(varId) bed.append(firstAnnot.dbSnpId) bed.append(htmlLink('uniProt', [acc])) pmids = firstAnnot.pmid.split(",") bed.append(htmlLink('pubmed', pmids)) # there is no filtering option yet for bigBed at the time of writing # allow to filter by database type based on the score if dbName=="trembl": bed[4] = "0" bed[8] = TREMBLCOLOR else: bed[4] = "1000" return bed def parseSizes(inFname): # keep protein lengths for later, we'll add one feature for every protein sequence seqSizes = {} for line in open(inFname): seq, size = line.rstrip("\n").split() seqSizes[seq] = int(int(size)/3) logging.debug("Parsed %s, found %d sequence sizes" % (inFname, len(seqSizes))) return seqSizes def makeSeqSizes(uniprotFa, outFname): " create a file with the lengths of the uniprot seqs " if not isfile(uniprotFa): raise Exception("Not found: %s" % uniprotFa) logging.debug("writing nucleotide lengths of %s to %s" % (uniprotFa, outFname)) cmd = "faSize %s -detailed | gawk '{$2=$2*3; print}'> %s" % (uniprotFa, outFname) #cmd = "faSize %s -detailed > %s" % (uniprotFa, outFname) run(cmd) def getMainIsoAcc(annot): acc = annot.mainIsoAcc return acc #if acc!="": #logging.debug("XX iso acc is %s" % acc) #return acc #if "-" in annot.acc: #return annot.acc #else: #logging.debug("XX faked iso acc is %s-1" % annot.acc) #return annot.acc+"-1" def writeProtAnnot(annotFnames, protBedName, seqSizes): " write bed file in prot coordinates, one bed per uniprot annotation " allFnames = ", ".join(annotFnames) logging.debug("converting UniProt tab %s to bed %s" % (allFnames, protBedName)) ofh = open(protBedName, "w") # index data by pos and write coords to BED doneAnnots = set() Record = None #bedNameToDb = {} for annotFname in annotFnames: logging.debug("Reading %s" % annotFname) for annot in iterTsvRows(open(annotFname)): acc = getMainIsoAcc(annot) if acc not in seqSizes: logging.warn("skipping %s, not in seqSizes" % annot.mainIsoAcc) continue annotName = annotToString(annot) annotPos = 3*(int(annot.begin)-1) annotPosEnd = 3*(int(annot.end)-1) if annotName not in doneAnnots: isoAcc = getMainIsoAcc(annot) newLine = "\t".join([isoAcc, str(int(annotPos)), str(int(annotPosEnd)), annotName]) ofh.write(newLine+"\n") doneAnnots.add(annotName) if Record is None: Record = namedtuple('fullSeqRec', annot._fields) # add one line for every accession, for accId, accSize in seqSizes.items(): annot = fakeAnnot(Record, accId, accSize) annotName = annotToString(annot) ofh.write("\t".join([accId, str(0), str(int(accSize*3)), annotName])) ofh.write("\n") ofh.close() def annotToString(annot): " create a unique key for an annotation " annotName = "|".join(annot).replace(" ", "_") assert("\t" not in annotName) return annotName def fakeAnnot(Record, seqId, seqSize): " create a faked uniProt annotation with just the full sequence " empty = {k:"" for k in Record._fields} #empty["comment"] = "Complete sequence %s mapped to the genome" % seqId empty["acc"] = seqId empty["begin"] = str(1) # UniProt is 1-based, like most protein resources empty["end"] = str(int(seqSize)) # and inclusive? empty["featType"] = fullSeqType fakeAnnot = Record(**empty) return fakeAnnot def parseAnnots(fnames, seqSizes): """ return uniprot annots indexed by unique string key adds pseudo features of featType fullSeq from first to last position. """ logging.debug("Parsing %s" % ",".join(fnames)) annotData = defaultdict(list) Record = None for fname in fnames: for annot in iterTsvRows(open(fname)): annotName = annotToString(annot) annotData[annotName].append(annot) # keep a copy of the fields for later if Record is None: fields = annot._fields if "isBlacklisted" in fields: fields = fields[:-1] Record = namedtuple('fullSeqRec', annot._fields) logging.debug("Found %d annotations" % len(annotData)) # now add entries for all protein sequences for seqId, seqSize in seqSizes.items(): annot = fakeAnnot(Record, seqId, seqSize) annotData[annotToString(annot)] = [annot] for key, annots in annotData.items(): if len(annots)!=1: # they do have complete dupes, maybe typos, e.g. on O44342 at position 50 logging.warn("Duplicated annotation, will be collapsed: %s" % repr( annots)) logging.debug("Added entries for full sequences, now at %d annotations" % len(annotData)) return annotData def convertToBigBed(bedFiles, genomeChromSizesFname, bigBedDir): for f in list(bedFiles.values()): f.close() myDir = dirname(__file__) bbFnames = [] for trackName, bedFile in bedFiles.items(): asFname = join(asDir, "bed12UniProtAnnotBgp.as") if trackName=="unipMut": asFname = join(asDir, "bed12UniProtMut.as") bedFname = bedFile.name assert(isfile(bedFname)) bbFname = join(bigBedDir, trackName+".new.bb") bbFnames.append(bbFname) cmd = "bedSort %(bedFname)s %(bedFname)s" % locals() run(cmd) addOpt = "" if "FullSeq" in bedFname: addOpt = "-extraIndex=uniProtId" cmd = "bedToBigBed -as=%(asFname)s %(bedFname)s %(genomeChromSizesFname)s %(bbFname)s -type=bed12+ -tab %(addOpt)s" % locals() run(cmd) assert(isfile(bbFname)) return bbFnames def moveBigBedToOutDir(bbFnames, outDir): " move all bigBeds over to the output directory, assures some leve of atomicity " for bbFname in bbFnames: targetFname = join(outDir, basename(bbFname)) logging.debug("moving %s to %s" % (bbFname, targetFname)) if isfile(targetFname): logging.debug("Removing old file %s" % targetFname) remove(targetFname) move(bbFname, targetFname) commentRe = re.compile(r' [(][^)]+[)]') def parseGroupAndAnnotate(genomeBedName, annotData, accToDb): """ parse bed file, pull out original UniProt fields, and group by genomePosition+category: group all beds with same coordinates and uniprot features together. this is important for two reasons: 1) the pslMap file may use multiple transcripts for each protein and as such can result in multiple mappings that all are exactly at the same place in the genome 2) trembl often annotates the same identical thing on two proteins that are otherwise almost identical we do not want to show it twice on the genome. """ logging.debug("parsing %s and removing overlaps (same uniProt feature 'category' and same chrom/start/end/blocks)" % genomeBedName) bedGroups = defaultdict(list) for line in open(genomeBedName): bed = line.rstrip("\n").split("\t") bedName = bed[3] annots = annotData[bedName] if len(annots)!=1: logging.debug("? duplicate annots: %s" % repr(annots)) annot = annots[0] # group identifier is: genome position and block structure + # feature type + comment + disease + mutation from -> to # we need the genome block structure, as a feature could be split over different exons # depending on the splicing keyFields = [bed[0], bed[1], bed[2], bed[10], bed[11]] # for mutations, we need the from -> to amino acids, otherwise two mutations at the same position # would get filtered out keyFields.append(annot.origAa) keyFields.append(annot.mutAa) keyFields.append(annot.featType) # we don't want to remove identical annotations with different comments # but UniProt maps things to other proteins and then adds (By similarity) # to avoid duplicates, remove these comments in comments comment = commentRe.sub("", annot.comment) # strips things like: ... (By similarity) (PubMed:32198291, PubMed:32272481) keyFields.append(comment) # I hope the comment does not include the accession # keyFields.append(annot.disease) keyFields.append(annot.disCode) #keyFields.append(annot.varId) # for chains this is the unique chain identifier, so cannot use it keyFields.append(annot.dbSnpId) bedKey = tuple(keyFields) dbName = accToDb.get(getMainIsoAcc(annot)) if dbName is None: dbName = accToDb[annot.acc] bedGroups[bedKey].append((dbName, bed, annot)) return bedGroups def uniprotLift(seqFname, cdnaSizesFname, annotFnames, genomeChromSizesFname, liftFname, outDir, accToDb, protSizes, options): " lift uniProt annotations from annotFname using liftFname and write bigBed files to outDir " annotFnameStr = ", ".join(annotFnames) logging.debug("lifting %(annotFnameStr)s using %(liftFname)s and writing bigBeds to %(outDir)s" % locals()) tempDir = "uniprotLift-"+basename(liftFname).split(".")[0]+".tmp" genomeBedName = join(tempDir, "genomeCoords.bed") protBedName = join(tempDir, "protCoords.bed") logging.debug("Cleaning directory %s, used for all temp files" % tempDir) cmd = "rm -rf %s; mkdir %s" % (tempDir, tempDir) run(cmd) writeProtAnnot(annotFnames, protBedName, protSizes) # lift the raw start-end information from the protein to the genome if (options.debug and isfile(genomeBedName)): logging.info("%s already exists, not remaking" % genomeBedName) else: logging.debug("lifting to genome coords") # added tawk '($16!=$17)' because in strPur2, with lots of tiny contigs, just a few amino acids got mapped # and after pslMap, nothing was left of a domain on the contig cmd = "bedToPsl %(cdnaSizesFname)s %(protBedName)s stdout | pslMap stdin %(liftFname)s stdout | tawk '($16!=$17)' | pslToBed stdin stdout | LC_COLLATE=C sort -k1,1 -k2,2n | bedFixBlockOverlaps /dev/stdin > %(genomeBedName)s" % locals() run(cmd) annotData = parseAnnots(annotFnames, protSizes) logging.debug("Adding extra fields to lifted coords") tracks = [ "unipFullSeq", # full protein sequence on genome "unipMut", # variants "unipStruct", # structural features "unipSplice", # splice variants "unipRepeat", # repeated regions "unipConflict", # sequence conflicts #"unipBlacklist", #blacklisted high-throughput-only features "unipLocSignal", "unipLocExtra", "unipLocTransMemb", "unipLocCytopl", # localization-related "unipDomain", "unipModif", "unipChain", "unipDisulfBond", # domains, etc "unipInterest", # regions of interest "unipOther", # all other features ] bedFiles = {} for t in tracks: bedFiles[t] = open(join(tempDir, t+".bed"), "w") bedGroups = parseGroupAndAnnotate(genomeBedName, annotData, accToDb) count = 0 # total number of annotations for bedName, beds in bedGroups.items(): beds.sort() # because "swissprot" comes before trembl alphabetically, all swissprot annotations are now first in this list # features are grouped by bedName, bedName includes chrom,start,end,type and various other fields # the aim is to not show something twice. UniProt projects features from one protein to others # and so often we risk ending up with duplicated annotation, when we project to the genome. # The grouping of features by (start, end, bedName, ...) removes the duplicates. # It is also dangerous, as we may drop some annotations accidentally. # The current system looks OK but may drop too much # if UniProt ever changes the way it uses the fields. dbName, bed, annot = beds[0] # of a group of identical annotations, we only show the first BED feature if len(beds)>1: logging.debug("Of the group '%s', only kept this feature: %s" % (bedName, repr(beds[0]))) logging.debug("Dropped these features because they have the same key: %s" % (repr(beds[1:]))) # pull out the original full Uniprot fields of this BED group bedName = bed[3] acc = annot.acc isMut = False color = None if annot.featType in ["mutagenesis site", "sequence variant"]: isMut = True ofh = bedFiles["unipMut"] elif annot.featType in ["splice variant"]: ofh = bedFiles["unipSplice"] elif annot.featType in ["strand","helix", "coiled-coil region", "turn"]: ofh = bedFiles["unipStruct"] # suggested by Regeneron: create four subtracks, for the subcell. localization indicators elif annot.featType in ["transmembrane region"]: ofh = bedFiles["unipLocTransMemb"] color = "0,150,0" # dark green, sent by Alejo Mujica in email from Oct 4 '17 elif annot.comment in ["Extracellular"]: ofh = bedFiles["unipLocExtra"] color = "0,110,180" # light-blue, alejo elif annot.comment in ["Cytoplasmic"]: ofh = bedFiles["unipLocCytopl"] color = "255,150,0" # light orange, alejo elif annot.featType in ["signal peptide"]: ofh = bedFiles["unipLocSignal"] color = "255,0,150" # light-red, alejo elif annot.featType in ["chain"]: ofh = bedFiles["unipChain"] elif annot.featType in ["repeat"]: ofh = bedFiles["unipRepeat"] elif annot.featType == "sequence conflict": ofh = bedFiles["unipConflict"] elif annot.featType in ["disulfide bond"]: ofh = bedFiles["unipDisulfBond"] elif annot.featType in ["domain", "zinc finger region", "topological domain"]: ofh = bedFiles["unipDomain"] elif annot.featType in ["glycosylation site", "modified residue", "lipid moiety-binding region"]: ofh = bedFiles["unipModif"] elif annot.featType in ["region of interest"]: ofh = bedFiles["unipInterest"] elif annot.featType==fullSeqType: ofh = bedFiles["unipFullSeq"] color = "0,0,170" # everything else else: ofh = bedFiles["unipOther"] bedFields = makeBedLine(dbName, annot, bed, isMut, color) # trembl features are always light-blue, no matter what, this overrides all other colors set before here if dbName=="trembl": bedFields[8] = TREMBLCOLOR bedLine = "\t".join(bedFields)+"\n" ofh.write(bedLine) count += 1 logging.info("%d features written to bedfiles in %s" % (count, tempDir)) bigBedDir = join(tempDir, "bigBed") if not isdir(bigBedDir): logging.debug("Making dir %s" % bigBedDir) mkdir(bigBedDir) bbFnames = convertToBigBed(bedFiles, genomeChromSizesFname, bigBedDir) moveBigBedToOutDir(bbFnames, outDir) cmd = "rm -rf %s" % tempDir run(cmd) def isNewer(url, localFname): " return true if url is more recent than localFname " # little helper class to download only the http HEAD class HeadRequest(urllib.request.Request): def get_method(self): return "HEAD" # get the mtime of the file on the server if url.startswith("http"): response = urllib.request.urlopen(HeadRequest(url)) timeStr = response.headers.get("Last-Modified") serverTime = datetime.datetime.fromtimestamp(time.mktime(time.strptime(timeStr, "%a, %d %b %Y %H:%M:%S GMT"))) elif url.startswith("ftp://"): # https://stackoverflow.com/questions/25724497/getting-servers-date-using-ftp u = urlparse(url) from ftputil import FTPHost # not found? run 'pip install ftputil' with FTPHost(u.netloc, "anonymous", "genome-www@soe.ucsc.edu") as host: mod_time = host.path.getmtime(u.path) logging.debug("mod time %s" % mod_time) serverTime = datetime.datetime.fromtimestamp(mod_time) # get the mtime of the flag file localTime = datetime.datetime.fromtimestamp(0) if os.path.isfile(localFname): t = os.path.getmtime(localFname) localTime = datetime.datetime.fromtimestamp(t) logging.debug("server time %s" % serverTime) logging.debug("local time %s" % localTime) if serverTime > localTime: logging.debug("URL %s is newer than %s" % (url, localFname)) return True else: logging.debug("URL %s is not newer than %s" % (url, localFname)) return False def run(cmd, ignoreErr=False): " run command, stop on error " logging.debug(cmd) cmd = "set -o pipefail; "+cmd ret = os.system(cmd) logging.info("Running: %s" % (cmd)) if ret!=0 and not ignoreErr: logging.error("command returned error: %s" % cmd) sys.exit(ret) def fastaMd5(fname): " sort gzipped fasta by id and return its md5 " logging.debug("Getting md5 of sorted seqs in %s" % fname) if fname.endswith(".gz"): catCmd = "zcat" else: catCmd = "cat" cmd = """%s %s | awk 'BEGIN{RS=">"}{print $1"\t"$2;}' | sort -k1,1 | md5sum""" % (catCmd, fname) proc = subprocess.Popen(cmd, stdout=PIPE, shell=True, encoding="latin1") md5 = proc.stdout.readline().split()[0] return md5[:15] def runQuery(db, query, usePublic=False): " run mysql query, yield tuples " logging.debug("Running query on db %s: %s" % (db, query)) cmd = ["hgsql", db, "-NB", "-e", query] if usePublic: cmd.extend(["-h", "genome-mysql.cse.ucsc.edu", "-u", "genome", "--password="]) proc = subprocess.Popen(cmd, encoding="latin1", stdout=subprocess.PIPE) try: for line in iter(proc.stdout.readline,''): yield line.rstrip().split("\t") except: while proc.poll() is None: # Process hasn't exited yet, let's wait some time.sleep(0.5) logging.error("could not run query: %s " %query ) assert(False) def hasNoMultiMappers(db, table): " return true if the gene table has at least one geneId that is mapped to two chroms " # this happens with Ensembl scaffold assemblies. It breaks our whole data model " logging.debug("Checking for multi-mapping genes") cmd = "hgsql %s -N -e 'select * from %s' | cut -f2- | sort -k1 | cut -f1 | uniq -c | sort -rn | head -n1" % \ (db, table) proc = subprocess.Popen(cmd, stdout=PIPE, shell=True, encoding="latin1") row1 = proc.stdout.readline().strip().split() count, gene = row1 count = int(count) if count!=1: logging.warn("Gene %s appears on %d chromosomes - this assembly has multi-mappers." % (gene, count)) return False return True def findBestGeneTable(db): " find the best gene table for a given organism and return it " for table in ["knownGene", "ensGene", "augustusGene"]: query = "DESCRIBE %s" % (table) cmd = "hgsql %s -e 'DESCRIBE %s'" % (db, table) ret = os.system(cmd) if ret==0: logging.debug("Best gene table for db %s is %s" % (db, table)) if hasNoMultiMappers(db, table): logging.info("%s: Gene table %s. No multi mappers." % (db, table)) return table return "blat" 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): "concat files and write to outFname. Also create a dict with seqId -> dbName " seqIds = set() logging.debug("Opening %s" % fname) ifh = gzip.open(fname, "rt") for line in ifh: if line.startswith(">"): seqId = line.lstrip(">").rstrip("\n").split()[0] seqIds.add(seqId) return seqIds def concatFiles(fnames, outFname): "concat files and write to outFname. Also create a dict with seqId -> dbName " logging.debug("infiles: %s, outfile: %s" % (fnames, outFname)) accToDb = {} ofh = open(outFname, "wt") for dbName, fname in fnames: if not isfile(fname): raise Exception("%s not found" % fname) ifh = gzip.open(fname, "rt") try: for line in ifh: ofh.write(line) if line.startswith(">"): seqId = line.lstrip(">").rstrip("\n").split()[0] accToDb[seqId] = dbName except: logging.error("Error reading %s" % fname) sys.exit(1) return accToDb def pslToProtPsl(inFname, outFname): " convert psl to prot coordinates " ofh = open(outFname, "w") for line in open(inFname): row = line.rstrip("\n").split("\t") matches, misMatches, repMatches, nCount, qNumInsert, qBaseInsert, tNumInsert, \ tBaseInsert, strand, qName, qSize, qStart, qEnd, tName, tSize, tStart, \ tEnd, blockCount, blockSizes, qStarts, tStarts = row qSize = int(qSize)/3 qStarts = ",".join([str(int(x)/3) for x in qStarts.rstrip(",").split(",")]) blockSizes = ",".join([str(int(x)/3) for x in blockSizes.rstrip(",").split(",")]) qStart = int(qStart)/3 qEnd = int(qEnd)/3 row = [matches, misMatches, repMatches, nCount, qNumInsert, qBaseInsert, tNumInsert, \ tBaseInsert, strand, qName, qSize, qStart, qEnd, \ tName, tSize, tStart, tEnd, \ blockCount, blockSizes, qStarts, tStarts] row = [str(x) for x in row] ofh.write("\t".join(row)) ofh.write("\n") ofh.close() def pslToBigPsl(mapFname, protFaFname, cdnaSizes, accToDb, accToAnnot, chromSizesFname, bigPslSwissFname, bigPslTremblFname): " rewrite the mapping psl to a bigPsl with extra fields, change colors and score " bpInputFname = mapFname.replace(".psl", ".pslInput") bpInputSwissFname = mapFname.replace(".psl", ".swissprot.pslInput") bpInputTremblFname = mapFname.replace(".psl", ".trembl.pslInput") # convert to bigPslInput cmd = "pslToBigPsl %(mapFname)s stdout | sort -k1,1 -k2,2n | bedFixBlockOverlaps /dev/stdin > %(bpInputFname)s" % locals() run(cmd) # rewrite the bigPslInput file: # - change colors and score # - add extra fields # - split into two subtracks: swissprot and trembl ofhSwiss = open(bpInputSwissFname, "w") ofhTrembl = None if bigPslTremblFname: ofhTrembl = open(bpInputTremblFname, "w") for line in open(bpInputFname): # chr1 11994 12144 B4E2Z4 1000 + 11994 12144 0 1 150, 0, 186 336 + 336 186, 249250621 150 0 0 0 0 row = line.rstrip("\n").split("\t") # change color and score acc = row[3] db = accToDb[acc] if db=="swissprot": color = SWISSPCOLOR score = "1000" ofh = ofhSwiss elif db=="trembl": color = TREMBLCOLOR score = "0" ofh = ofhTrembl else: assert(False) row[8] = color row[4] = score recAnnot = accToAnnot.get(acc, None) if recAnnot is None: # many trembl sequences have no annotations at all # logging.debug("No annotations for %s" % acc) continue row.append( recAnnot.acc ) # uniprot accessions, field 26 row.append( recAnnot.name ) # uniprot curated name, field 27 if db=="trembl": dbDesc = "Unreviewed (TrEMBL)" else: dbDesc = "Manually reviewed (Swiss-Prot)" row.append( dbDesc ) # db source, field 28 row.append( recAnnot.accList.replace("|", ", ")) # uniprot alternative IDs, #29 row.append( recAnnot.isoIds.replace("|", ", ")) # isoform IDs, #30 fullNames = recAnnot.protFullNames row.append( fullNames ) # uniprot full names, #31 # merge all these fields into one list of "gene synonyms" #14 geneSynonyms #15 isoNames #16 geneOrdLocus #17 geneOrf syns = recAnnot.geneSynonyms.split("|") syns.extend(recAnnot.isoNames.split("|")) syns.extend(recAnnot.geneOrdLocus.split("|")) syns.extend(recAnnot.geneOrf.split("|")) syns = [x.strip() for x in syns if x.strip()!=""] geneNames = recAnnot.geneName.replace("; ", ", ").replace("|", "; ") # relevant later but need it now shortNames = recAnnot.protShortNames.replace("|", "; ") if shortNames=="": shortNames = geneNames # in cov-2 pre-release, one protein only has a geneName, nothing else if shortNames=="" and len(syns)>0: shortNames = syns[0] if shortNames=="": shortNames = fullNames # in cov-2 pre-release, some proteins have no short name row.append( shortNames ) # uniprot short names #32 row.append( recAnnot.protAltFullNames.replace("; ", ", ").replace("|", "; ")) # uniprot alt names #33 row.append( recAnnot.protAltShortNames.replace("; ", ", ").replace("|", "; ")) # uniprot alt names #34 row.append( geneNames ) # #35 synStr = "; ".join(syns) row.append( synStr ) # #field 36: gene synonyms funcText = recAnnot.functionText.replace("|", "<li>").replace("-<li>-", "-|-") # address: [ILMVF]-Q-|-[SGACN] row.append( "<ul><li>"+funcText+"</ul>" ) # #37 row.append( recAnnot.hgncSym.replace("|", ", ")) row.append( recAnnot.hgncId.replace("|", ", ")) row.append( recAnnot.refSeq.replace("|", ", ")) row.append( recAnnot.refSeqProt.replace("|", ", ")) row.append( recAnnot.entrezGene.replace("|", ", ")) row.append( recAnnot.ensemblGene.replace("|", ", ")) row.append( recAnnot.ensemblProt.replace("|", ", ")) row.append( recAnnot.ensemblTrans.replace("|", ", ")) ofh.write("\t".join(row)) ofh.write("\n") ofhSwiss.close() if ofhTrembl: ofhTrembl.close() asFname = join(asDir, "bigPslUniprot.as") cmd = "bedToBigBed -extraIndex=acc -type=bed12+ -as=%s -tab %s %s %s" % \ (asFname, bpInputSwissFname, chromSizesFname, bigPslSwissFname) run(cmd) if ofhTrembl: cmd = "bedToBigBed -extraIndex=acc -type=bed12+ -as=%s -tab %s %s %s" % \ (asFname, bpInputTremblFname, chromSizesFname, bigPslTremblFname) run(cmd) def parseAnnot(tabDir, taxId, doTrembl): " parse the uniprot record-level annotations, like xrefs etc. return dict acc -> namedtuple " 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)): acc = getMainIsoAcc(row) ret[acc] = row return ret def blatProteinsKeepBest(fullFaFname, db, mapFname): " blat the proteins directly, only good for small genomes where we don't have any gene models " cmd = "blat -q=prot -t=dnax /gbdb/{db}/{db}.2bit {inf} stdout -noHead | sort -k10 | pslReps stdin stdout /dev/null " \ "-minAli=0.99 -nohead | pslProtCnv > {out} ".format(db=db, inf=fullFaFname, out=mapFname) run(cmd) def updateUniprot(args, onlyDbs, taxIdDbs, options): uprotDir = options.uniprotDir tmpDir = join(uprotDir, "download.tmp") tabDir = options.tabDir mapDir = options.mapDir bigBedDir = options.bigBedDir faDir = options.faDir doTrembl = not options.skipTrembl relFname = join(tabDir, "version.txt") if not options.skipDownload and not options.skipParse: + if not isdir(tmpDir): + os.makedirs(tmpDir) + localFname = join(uprotDir, "uniprot_sprot.xml.gz") if not isNewer(upUrl, localFname): logging.info("files at %s are not newer than file in %s, nothing to do. Specify -l to skip this check." % (upUrl, localFname)) sys.exit(0) # use lftp to update uniprot and download only changed files cmd = 'lftp ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/complete/ -e "lcd %s && mirror . -P 3 --use-pget-n=4 --exclude-glob *.dat.gz && exit"' % tmpDir run(cmd) # make sure that download interruptions don't mess up the files in the target dir logging.info("Moving files from %s to %s" % (tmpDir, uprotDir)) file_names = os.listdir(tmpDir) for file_name in file_names: shutil.move(os.path.join(tmpDir, file_name), uprotDir) # only use the first four words of version string # 'UniProt Knowledgebase Release 2017_09 consists of:' relString = open(join(uprotDir, "reldate.txt")).read().splitlines()[0] relString = " ".join(relString.split()[:4]) relFh = open(relFname, "w") relFh.write(relString) relFh.close() logging.debug("Wrote release version string '%s' to %s" % (relString, relFname)) logging.info("Converting uniprot from %s to %s, using maps in %s" % (uprotDir, tabDir, mapDir)) if not options.skipParse: # parsing the UniProt XML files logging.info("Converting uniprot XML from %s to tab-sep in %s, this can take 2-3 days!! Specify -p to skip this." % (uprotDir, tabDir)) dbToTaxIds = {} if onlyDbs: for taxId, dbs in taxIdDbs: for db in dbs: dbToTaxIds[db] = taxId taxIds = [dbToTaxIds[db] for db in onlyDbs] else: taxIds = [x for (x,y) in taxIdDbs] # just the taxIds themselves taxIdStr = ",".join([str(x) for x in taxIds]) 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) # 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() userName = os.getlogin() updateOfh = open(join(bigBedDir, "trackVersionUpdate.sh"), "w") for taxId, dbs in taxIdDbs: if onlyDbs is not None and len(set(dbs).intersection(onlyDbs))==0: continue logging.debug("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)) ) fullFaFname = join(faDir, str(taxId)+".fa") accToDb = concatFiles(faFnames, fullFaFname) md5 = fastaMd5(fullFaFname) logging.debug("md5 is %s" % md5) # 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 chromSizesFname = "/hive/data/genomes/%s/chrom.sizes" % db logging.debug("Working on %s" % db) geneTable = findBestGeneTable(db) mapFname = join(mapDir, "%(taxId)s_%(db)s_%(geneTable)s_%(md5)s.psl" % locals()) if isfile(mapFname): logging.debug("%s already exists" % mapFname) else: logging.debug("%s does not exist" % mapFname) if geneTable=="blat": logging.error("Could not find a gene table for %s, using BLAT to map proteins" % db) blatProteinsKeepBest(fullFaFname, db, mapFname) else: cmd = "time makeUniProtPsl.sh %(fullFaFname)s %(db)s %(geneTable)s %(cluster)s %(mapFname)s" % locals() run(cmd) # get the sizes of the cDNAs and write them to a file cdnaSizesFname = fullFaFname.replace(".fa", ".cdnaSize") makeSeqSizes(fullFaFname, cdnaSizesFname) protSizes = parseSizes(cdnaSizesFname) bigPslDir = join(bigBedDir, db) if not isdir(bigPslDir): os.makedirs(bigPslDir) bigPslSwissFname = join(bigPslDir, "%(geneTable)s_%(md5)s.swissprot.new.bb" % locals()) if doTrembl: bigPslTremblFname = join(bigBedDir, db, "%(geneTable)s_%(md5)s.trembl.new.bb" % locals()) else: bigPslTremblFname = None accToAnnot = parseAnnot(tabDir, taxId, doTrembl) pslToBigPsl(mapFname, fullFaFname, protSizes, accToDb, accToAnnot, chromSizesFname, bigPslSwissFname, bigPslTremblFname) if options.onlyMap: continue dbBigBedDir = join(bigBedDir, db) if not isdir(dbBigBedDir): mkdir(dbBigBedDir) tabFnames = ["tab/swissprot.%d.annots.tab" % taxId] if doTrembl: tabFnames.append( "tab/trembl.%d.annots.tab" % taxId ) uniprotLift(fullFaFname, cdnaSizesFname, tabFnames, chromSizesFname, mapFname, dbBigBedDir, accToDb, protSizes, options) # add a trackVersion entry assert("'" not in versionString) assert('"' not in versionString) sql = '''INSERT INTO trackVersion values (NULL, '%s', 'uniprot', '%s', '%s', now(), '%s', 'doUniprot Otto', '%s');''' % (db, userName, versionString, taxId, upUrl) cmd = """hgsql hgFixed -e "%s";""" % sql run(cmd) updateOfh.write("%s\n"%cmd) updateOfh.close() logging.debug("Created %s for cluster-admin's auto-pusher" % updateOfh.name) def makeSymlink(target, linkName): if isfile(linkName): #logging.debug("%s already exists" % linkName) os.remove(linkName) if target is None or not isfile(target): logging.error("Cannot symlink: %s does not exist" % str(target)) return targetPath = abspath(target) linkPath = abspath(linkName) cmd = "ln -sf %(targetPath)s %(linkPath)s " % locals() run(cmd) def findForMask(fileMask): fnames = glob.glob(fileMask) if len(fnames)==0: logging.error("NOT FOUND: %s" % fileMask) return None assert(len(fnames)==1) fname = fnames[0] return fname def makeLinks(bigBedDir, faDir, onlyDbs, taxIdDbs): " check the /gbdb symlinks " for taxId, dbs in taxIdDbs: fullFaFname = join(faDir, str(taxId)+".fa") md5 = None for db in dbs: if onlyDbs is not None and db not in onlyDbs: continue if md5 is None: md5 = fastaMd5(fullFaFname) # do only once per organism, seqs don't change dbBigBedDir = join(bigBedDir, db) # find the bigBed files bbTargetDir = join("/gbdb", db, "uniprot") if not isdir(bbTargetDir): logging.info("Making %s" % bbTargetDir) os.mkdir(bbTargetDir) bbFnames = glob.glob(join(dbBigBedDir, "*.bb")) logging.debug("Found %d bigBed files in %s" % (len(bbFnames), dbBigBedDir)) # and create links to them for bbFname in bbFnames: if "unipFullSeq" in bbFname: continue # this has now been replaced by the bigPsl file if "_" in bbFname: continue # the bigPsl files need special treatment later bbLink = join(bbTargetDir, basename(bbFname)) makeSymlink(bbFname, bbLink) # find the bigPsl files bigPslMaskSwiss = join(bigBedDir, db, "*_%(md5)s.swissprot.bb" % locals()) bigPslSwissprotFname = findForMask(bigPslMaskSwiss) bigPslMaskTrembl = join(bigBedDir, db, "*_%(md5)s.trembl.bb" % locals()) bigPslTremblFname = findForMask(bigPslMaskTrembl) # and create links to them bbLinkSwiss = join(bbTargetDir, "unipAliSwissprot.bb") makeSymlink(bigPslSwissprotFname, bbLinkSwiss) bbLinkTrembl = join(bbTargetDir, "unipAliTrembl.bb") makeSymlink(bigPslTremblFname, bbLinkTrembl) def checkPsl(pslName, faFname): " checking the qNames in a psl input (!) file " if pslName is None: logging.debug("Not found: %s" % pslName) return seqIds = parseFaIds(faFname) accs = set([x.split("-")[0] for x in seqIds]) qNames = Counter() pslCount = 0 for line in open(pslName): row = line.rstrip("\n").split("\t") qName = row[3] qNames[qName]+=1 pslCount += 1 accCount = len(accs) if accCount==0: logging.error("No single sequence found in %s" % faFname) return unmappedIds = seqIds - set(qNames) unmapCount = len(unmappedIds) unmappedAccs = set([x.split("-")[0] for x in seqIds]) - set([x.split("-")[0] for x in qNames]) unmappedAccCount = len(unmappedAccs) multiMapCount = 0 multiMapIds = {} for qName, qCount in qNames.items(): if qCount > 1: multiMapCount += 1 multiMapIds[qName] = qCount worstMultis = [] for qName, qCount in qNames.most_common()[:10]: if qCount > 1: worstMultis.append("%s (%d times)" % (qName, qCount)) qNameCount = len(qNames) seqCount = len(seqIds) unmappedAccPerc = 100.0 * float(unmappedAccCount) / accCount unmappedPerc = 100.0 * float(unmapCount) / seqCount multiPerc = 100.0 * float(len(multiMapIds)) / seqCount multiMapCount = len(multiMapIds) logging.info("%(faFname)s %(pslName)s: %(pslCount)d psls, %(accCount)s accessions, %(seqCount)d seqs, %(qNameCount)d mapped seqs" % locals()) #logging.info("%s: %d lines in PSL, %d different query names mapped" % (pslName, pslCount, len(qNames))) logging.info("Unmapped accessions: %d (%0.02f %%) Unmapped seqs: %d (%.02f %%) Multi-mapped seqs: %d (%0.02f %%)" % (unmappedAccCount, unmappedAccPerc, unmapCount, unmappedPerc, multiMapCount, multiPerc)) logging.info("10 examples of unmapped sequences: %s" % " ".join(list(unmappedIds)[:10])) logging.info("10 examples of unmapped accessions: %s" % " ".join(list(unmappedAccs)[:10])) logging.info("10 worst multi-mappers: %s" % (", ".join(worstMultis))) def bedInfo(fname): " return number of features and coverage in bigBed file " cmd = "bigBedInfo %s" % fname proc = subprocess.Popen(cmd, encoding="latin1", stdout=PIPE, shell=True) for line in proc.stdout: # itemCount: 2,662 # basesCovered: 341,994 if line.startswith("itemCount"): itemCount = int(line.split()[1].replace(",","")) if line.startswith("basesCovered"): basesCovered = int(line.split()[1].replace(",","")) return (itemCount, basesCovered) def flipNewVersionCheckMaxDiff(bigBedDir, maxDiff, force, onlyDbs, taxIdDbs): " rename all .new.bb files to .bb files, but only if they do not diff in more than 10% " for taxId, dbs in taxIdDbs: for db in dbs: if onlyDbs is not None and db not in onlyDbs: continue for newFname in glob.glob(join(bigBedDir, db, "*.new.bb")): if force: logging.debug("Not doing maxDiff check, --force specified") continue oldFname = newFname.replace(".new.bb", ".bb") if not isfile(oldFname): logging.warn("%s does not exist, cannot do max-diff check" % oldFname) else: newCount, newCov = bedInfo(newFname) oldCount, oldCov = bedInfo(oldFname) if oldCount > 1500 and (newCount-oldCount) / oldCount > maxDiff: logging.error("%s: old feature count %d, new feature count %d. Too big difference. Stopping now." % (oldFname, oldCount, newCount)) sys.exit(1) if oldCov > 10000 and (oldCov-newCov) / oldCov > maxDiff: logging.error("%s: old coverage %d, new coverage %d. Too big difference. Stopping now." % (oldFname, oldCov, newCov)) sys.exit(1) # only when everything is OK, do the flip for taxId, dbs in taxIdDbs: for db in dbs: if onlyDbs is not None and db not in onlyDbs: continue for fname in glob.glob(join(bigBedDir, db, "*.new.bb")): if isfile(fname): oldFname = fname.replace(".new.bb", ".bb") logging.debug("Renaming %s -> %s" % (fname, oldFname)) shutil.move(fname, oldFname) def mapQa(tabDir, faDir, mapDir, onlyDbs, taxIdDbs): " check the pslMap files " for taxId, dbs in taxIdDbs: fullFaFname = join(faDir, str(taxId)+".fa") md5 = None for db in dbs: if onlyDbs is not None and db not in onlyDbs: continue if md5 is None: md5 = fastaMd5(fullFaFname) # do only once per organism, seqs don't change # find the psl file pslMask = join(mapDir, "%(taxId)s_%(db)s_*_%(md5)s.swissprot.pslInput" % locals()) pslName = findForMask(pslMask) faFname = join(tabDir, "swissprot.%s.fa.gz" % taxId) checkPsl(pslName, faFname) pslMask = join(mapDir, "%(taxId)s_%(db)s_*_%(md5)s.trembl.pslInput" % locals()) pslName = findForMask(pslMask) faFname = join(tabDir, "trembl.%s.fa.gz" % taxId) checkPsl(pslName, faFname) def createVersionFiles(tabDir, bigBedDir, onlyDbs, taxIdDbs): " update the release version file. This is the indicator for the auto-archiver to create a new archive " relFname = join(tabDir, "version.txt") versionString = open(relFname).read() shortVersion = versionString.split()[-1] if (len(shortVersion)!=7): logging.warning("UniProt Version string is not seven characters long") for taxId, dbs in taxIdDbs: for db in dbs: if onlyDbs is not None and db not in onlyDbs: continue dbDir = join(bigBedDir, db) if not isdir(dbDir): continue versionOfh = open(join(dbDir, "version.txt"), "w") versionOfh.write(versionString) versionOfh.close() linkName = join("/gbdb", db, "uniprot", "version.txt") makeSymlink(versionOfh.name, linkName) logging.debug("Wrote release string to %s, symlink from %s" % (versionOfh.name, linkName)) def delFlag(): os.remove(flagFname) def main(): global flagFname args, options = parseArgs() onlyDbs = None if options.dbs: onlyDbs = set(options.dbs.split(",")) taxIdDbs = getTaxIdDbs() if options.mapQa: mapQa(options.tabDir, options.faDir, options.mapDir, onlyDbs, taxIdDbs) sys.exit(0) # create a lock file flagFname = join(options.uniprotDir, "doUniprot.lock") atexit.register(delFlag) if isfile(flagFname) and not options.debug: logging.error("%s exists. Is a doUniprot process already running ? " "If not, remove the flag file and restart this script." % flagFname) sys.exit(1) open(flagFname, "w") if not options.onlyLinks: updateUniprot(args, onlyDbs, taxIdDbs, options) flipNewVersionCheckMaxDiff(options.bigBedDir, 0.1, options.force, onlyDbs, taxIdDbs) if not options.skipLinks: makeLinks(options.bigBedDir, options.faDir, onlyDbs, taxIdDbs) createVersionFiles(options.tabDir, options.bigBedDir, onlyDbs, taxIdDbs) os.remove(flagFname) main()