880ad5f5e5b759ac892553985d6ab3721452e870 max Thu Dec 11 16:58:34 2014 -0800 largely revamped beacon, refs #14393 diff --git src/beacon/query src/beacon/query deleted file mode 100755 index e951874..0000000 --- src/beacon/query +++ /dev/null @@ -1,156 +0,0 @@ -#!/usr/bin/env python -# a beacon for LOVD, ClinVar, UniProt and HGMD, see ga4gh.org/#/beacon -import cgi, subprocess, sys, cgitb, os, datetime, gdbm -from os.path import join -cgitb.enable() - -# dir for bigBed files -DATADIR = "/gbdb/hg19/beacon" - -trackFnames = { -"clinvar" : "clinvarMain.bb", -"uniprot" : "spMut.bb", -"hgmd" : "hgmd.bb", -"lovd" : "lovd.bb" -} - -MAXREQPERDAY=1000 - -def checkReqCount(): - " check if IP address has exceeded max number of requests today " - accCounts = gdbm.open("accCounts.gdbm", "c") - - ipAddr = os.environ["REMOTE_ADDR"] - dateStr = datetime.date.today().isoformat() - accKey = ipAddr+"|"+dateStr - try: - oldCount = int(accCounts[accKey]) - except KeyError: - oldCount = 0 - - newCount = oldCount+1 - if newCount > MAXREQPERDAY: - print("No more requests left today") - sys.exit(0) - - accCounts[accKey]=str(newCount) - accCounts.close() - -def getTrackFromScript(): - " check which track name occurs in script name, allows symlinks to cgi that include the track name " - name = sys.argv[0].lower() - for trackName in trackFnames: - if trackName in name: - return trackName - -def errAbort(errMsg=None): - print "<html><body><pre>" - if errMsg is not None: - print "error: "+errMsg - - host = os.environ.get("HTTP_HOST", "") - - print """ - A Genomics Alliance Beacon at UCSC. In collaboration with the Leiden Open Variation Database. - - Data provided by Leiden Open Variation (Ivo Fokkema), NCBI ClinVar, UniProt, Biobase/HGMD. - - parameters: - - chrom: either a number or X/Y or UCSC format like "chr1" or INSDC accession like CM000663.1 - - track: clinvar, uniprot or lovd - - pos: 0-based position, like 10042537 - - allele: A,C,T,G. Ignored for UniProt (UniProt stores only amino acid changes) - - If the allele was found in the track, returns the string "Yes", otherwise "No". - - For UniProt and HGMD, an overlap with an amino acid change is enough to trigger "Yes". - Therefore, the alleles are ignored for UniProt and HGMD. - For LOVD, an overlap with a deletion or duplication returns "Maybe". - - Examples: - http://%(host)s/cgi-bin/beacon/query?track=lovd&chrom=chr1&pos=808921&allele=T - http://%(host)s/cgi-bin/beacon/query?track=clinvar&chrom=chr1&pos=10042537&allele=T - http://%(host)s/cgi-bin/beacon/query?track=uniprot&chrom=CM000663.1&pos=977027&allele=T - - Alternative URLs are queryLovd, queryClinvar and queryUniprot, e.g. - http://%(host)s/cgi-bin/beacon/queryLovd?chrom=1&pos=2160493&allele=T - - returns "Maybe" - - Number of variants: - - LOVD: 180k - - ClinVar: 30k - - HGMD: 90k - - UniProt: 47k - """ % locals() - print "</body></pre></html>" - sys.exit(0) - -# parse INSDC -> ucsc chrom name table -insdcToUcsc = {} -for line in open("insdcToUcsc.tab"): - name, chrom = line.rstrip("\n").split("\t") - insdcToUcsc[name] = chrom - -print "Content-Type: text/html\n" -form = cgi.FieldStorage() - -track = form.getfirst("track") -if track is None: - track = getTrackFromScript() - -chrom = form.getfirst("chrom") -pos = form.getfirst("pos") -nucl = form.getfirst("allele") - -# check parameters -if chrom==None: - errAbort() -if not pos.isdigit(): - errAbort("pos is not an int") -if chrom.isdigit() or chrom in "XY": - chrom = "chr"+chrom - -# convert chrom from genbank acc to ucsc string like "chr1" -if chrom not in insdcToUcsc.values(): - chrom = insdcToUcsc.get(chrom, None) -if chrom==None: - errAbort( "ERROR: illegal chrom identifier %s" % chrom) - -if track not in trackFnames: - errAbort("illegal track name") - -if nucl==None or len(nucl)!=1 or nucl[0].lower() not in "actg": - errAbort("invalid allele, can only be ACTG") - -checkReqCount() - -# change parameters slightly -nucl = nucl.upper() -fname = join(DATADIR, trackFnames[track]) -pos = int(pos) - -# run tool -args = ["./bigBedToBed", fname, "stdout", "-chrom=%s"%chrom, "-start=%d"%pos, "-end=%d"%(pos+1)] -p = subprocess.Popen(args, stdout=subprocess.PIPE) - -# parse output -isMaybe = False -for line in p.stdout: - fs = line.rstrip("\n").split("\t") - if fs[-1]=="maybe": - isMaybe=True - if int(fs[2])-int(fs[1])!=1 and track != "uniprot": - continue - if fs[3][-1]==nucl or track in ["hgmd", "uniprot"]: - print "Yes" - sys.exit(0) -if isMaybe: - print "Maybe" - sys.exit(0) - -# cannot check exit code, no overlap triggers error -#if p.returncode != 0: - #print("Unable to run command") - #sys.exit(0) - -print "No"