a2be3b4aee520bfa47602e50b4825b43e9854690 max Thu Dec 11 17:06:43 2014 -0800 small formatting changes in last commit diff --git src/hg/hgBeacon/hgBeacon src/hg/hgBeacon/hgBeacon index ab3e4bc..6e9b759 100755 --- src/hg/hgBeacon/hgBeacon +++ src/hg/hgBeacon/hgBeacon @@ -1,133 +1,133 @@ #!/usr/bin/env python # A beacon for tracks we can't distribute: LOVD and HGMD # see ga4gh.org/#/beacon and redmine 14393 # Please do read the GA4GH website on the background of this CGI, otherwise # the functionality provided does not make any sense, without the context import cgi, subprocess, sys, cgitb, os, datetime from os.path import join, isfile cgitb.enable() # deactivate this on the RR, but useful for debugging: prints a http header on errors # dir for bigBed files and insdcToUcsc.tab DATADIR = "/gbdb/hg19/beacon" 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. + A Genomics Alliance Beacon at UCSC. - Data provided by Leiden Open Variation (Ivo Fokkema), Biobase/HGMD. + In collaboration with the Leiden Open Variation Database and Biobase/HGMD. parameters: - genome: ignored, always hg19/GRCh37 - chrom: either a number or X/Y or UCSC format like "chr1" or INSDC accession like CM000663.1 - dataset: either lovd or hgmd. This parameter can also be called "track" (backwards compatible) - pos: 0-based position, like 10042537 - allele: A,C,T,G. Ignored for HGMD (HGMD gives out only positions, not the - alleles, they are on the website) + alleles, they are only on their website) If the allele was found, returns the string "Yes", otherwise "No". For HGMD, an overlap with a variant is enough to trigger a "Yes", we don't have the alleles. For LOVD, an overlap with a short deletion or duplication returns "Maybe". Examples: http://%(host)s/cgi-bin/hgBeacon/hgBeacon?dataset=lovd&chrom=chr1&pos=808921&allele=A - returns "Yes" http://%(host)s/cgi-bin/hgBeacon/hgBeacon?dataset=hgmd?chrom=1&pos=2160493&allele=T - returns "Maybe", because the allele is in a small deletion Number of variants: - LOVD: 180k, version from Nov 2014 - HGMD: 97k, "public" version from end of 2014 """ % locals() print "</body></pre></html>" sys.exit(0) # parse INSDC -> ucsc chrom name table insdcToUcsc = {} for line in open(join(DATADIR, "insdcToUcsc.tab")): name, chrom = line.rstrip("\n").split("\t") insdcToUcsc[name] = chrom print "Content-Type: text/html\n" form = cgi.FieldStorage() track = form.getfirst("dataset") if track is None: track = form.getfirst("track") chrom = form.getfirst("chrom") pos = form.getfirst("pos") allele = form.getfirst("allele") # check parameters if chrom==None: errAbort() if pos==None or not pos.isdigit(): errAbort("pos parameter is not an int") if allele==None: errAbort("need allele parameter") if chrom.isdigit() or chrom in "XY": chrom = "chr"+chrom if track is None: errAbort("missing dataset parameter") if not track.isalpha(): errAbort("dataset/track parameter must contain only alphanumberic characters") # 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) trackFname = join(DATADIR, track+".bb") if not isfile(trackFname): errAbort("illegal dataset/track name %s" % track) if allele==None or len(allele)!=1 or allele[0].lower() not in "actg": errAbort("invalid allele, can only be a single letter and ACTG") # change parameters slightly allele = allele.upper() pos = int(pos) # run tool args = ["./bigBedToBed", trackFname, "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 bedLen = int(fs[2])-int(fs[1]) # uniprot features are longer than one basepair if bedLen != 1 and track != "uniprot": continue # hgmd and uniprot do not have alleles if fs[3][-1]==allele or track in ["hgmd", "uniprot"]: print "Yes" sys.exit(0) if isMaybe: print "Maybe" sys.exit(0) # cannot check exit code, no overlap makes bigBedToBed return an error, too # so the following does not work #if p.returncode != 0: #print("Unable to run command") #sys.exit(0) print "No"