Files
galaxy/tools/rgenetics/rgRegion.py
T

81 lines
2.8 KiB
Python

"""
released under the terms of the LGPL
copyright ross lazarus August 2007
for the rgenetics project
Special galaxy tool for the camp2007 data
Allows grabbing arbitrary columns from an arbitrary region
Needs a mongo results file in the location hardwired below or could be passed in as
a library parameter - but this file must have a very specific structure
rs chrom offset float1...floatn
called as
<command interpreter="python">
rsRegion.py $infile '$cols' $r $tag $out_file1
</command>
cols is a delimited list of chosen column names for the subset
r is a ucsc location region pasted into the tool
"""
import sys,string
trantab = string.maketrans(string.punctuation,'_'*len(string.punctuation))
print >> sys.stdout, '##rgRegion.py started'
if len(sys.argv) <> 6:
print >> sys.stdout, '##!expected params in sys.argv, got %d - %s' % (len(sys.argv),sys.argv)
sys.exit(1)
print '##got %d - %s' % (len(sys.argv),sys.argv)
# quick and dirty for galaxy - we always get something for each parameter
fname = sys.argv[1]
wewant = sys.argv[2].split(',')
region = sys.argv[3].lower()
tag = sys.argv[4].translate(trantab)
ofname = sys.argv[5]
myname = 'rgRegion'
if len(wewant) == 0: # no columns selected?
print >> sys.stdout, '##!%s: no columns selected - cannot run' % myname
sys.exit(1)
try:
f = open(fname,'r')
except: # bad input file name?
print >> sys.stdout, '##!%s unable to open file %s' % (myname, fname)
sys.exit(1)
try: # TODO make a regexp?
c,rest = region.split(':')
c = c.replace('chr','') # leave although will break strict genome graphs
rest = rest.replace(',','') # remove commas
spos,epos = rest.split('-')
spos = int(spos)
epos = int(epos)
except:
print >> sys.stdout, '##!%s unable to parse region %s - MUST look like "chr8:10,000-100,000' % (myname,region)
sys.exit(1)
print >> sys.stdout, '##%s parsing chrom %s from %d to %d' % (myname, c,spos,epos)
res = []
cnames = f.next().strip().split() # column titles for output
linelen = len(cnames)
wewant = [int(x) - 1 for x in wewant] # need col numbers base 0
for n,l in enumerate(f):
ll = l.strip().split()
thisc = ll[1]
thispos = int(ll[2])
if (thisc == c) and (thispos >= spos) and (thispos <= epos):
if len(ll) == linelen:
res.append([ll[x] for x in wewant]) # subset of columns!
else:
print >> sys.stdout, '##! looking for %d fields - found %d in ll=%s' % (linelen,len(ll),str(ll))
o = file(ofname,'w')
res = ['%s\n' % '\t'.join(x) for x in res] # turn into tab delim string
print >> sys.stdout, '##%s selected and returning %d data rows' % (myname,len(res))
head = [cnames[x] for x in wewant] # ah, list comprehensions - list of needed column names
o.write('%s\n' % '\t'.join(head)) # header row for output
o.write(''.join(res))
o.close()
f.close()