opensource/galaxy
3
0
Fork 0
mirror of https://github.com/galaxyproject/galaxy.git synced 2026-09-24 16:30:27 +08:00
Code Issues Packages Projects Releases Wiki Activity

Adding tool to fetch closest upstream and downstream features for a given set of intervals.

Browse Source 
Added: Tool xml and py files, test data.
Modified: tool_conf.xml.sample.
This commit is contained in:
Guruprasad Anada
2008-03-05 20:08:01 +00:00
parent f04eb986f3
commit b032b71bc0
3 changed files with 278 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files
tool_conf.xml.sample
+1
View File
@@ -103,6 +103,7 @@
<tool file="new_operations/cluster.xml" id="cluster" />
<tool file="new_operations/join.xml" />
<tool file="new_operations/get_flanks.xml" />
<tool file="new_operations/flanking_features.xml" />
</section>
<section name="Statistics" id="stats">
<tool file="stats/gsummary.xml" />
tools/new_operations/flanking_features.py
+208
View File
@@ -0,0 +1,208 @@
#! /usr/bin/python
#by: Guruprasad Ananda
"""
This tool finds the closest up- and/or down-stream feature in input2 for every interval in input1.
usage: %prog input1 input2 out_file direction
-1, --cols1=N,N,N,N: Columns for chrom, start, end, strand in file1
-2, --cols2=N,N,N,N: Columns for chrom, start, end, strand in file2
"""
import sys, os, tempfile, commands
import pkg_resources; pkg_resources.require( "bx-python" )
from bx.cookbook import doc_optparse
from galaxy.tools.util.galaxyops import *
def stop_err( msg ):
sys.stderr.write( msg )
sys.exit()
def main():
# Parsing Command Line here
options, args = doc_optparse.parse( __doc__ )
try:
chr_col_1, start_col_1, end_col_1, strand_col_1 = parse_cols_arg( options.cols1 )
chr_col_2, start_col_2, end_col_2, strand_col_2 = parse_cols_arg( options.cols2 )
infile1, infile2, out_file, direction = args
if strand_col_1 < 0:
strand = "+" #if strand is not defined, default it to +
if strand_col_2 < 0:
fstrand = ""
except:
stop_err( "Metadata issue, correct the metadata attributes by clicking on the pencil icon in the history item." )
try:
fi = open(infile1,'r') #fi: primary interval file
ff = open(infile2,'r') #ff: feature file
except:
stop_err( "Unable to open input file" )
try:
fo = open(out_file,'w')
except:
stop_err( "Unable to open output file" )
tmpfile1 = tempfile.NamedTemporaryFile()
tmpfile2 = tempfile.NamedTemporaryFile()
try:
#Sort the features file based on decreasing end positions
command_line1 = "sort -f -n -r -k " + str(end_col_2+1) + " -o " + tmpfile1.name + " " + infile2
#Sort the features file based on increasing start positions
command_line2 = "sort -f -n -k " + str(start_col_2+1) + " -o " + tmpfile2.name + " " + infile2
except Exception, exc:
stop_err( 'Initialization error -> %s' %str(exc) )
error_code1, stdout = commands.getstatusoutput(command_line1)
error_code2, stdout = commands.getstatusoutput(command_line2)
if error_code1 != 0:
stop_err( "Sorting input dataset resulted in error: %s: %s" %( error_code1, stdout ))
if error_code2 != 0:
stop_err( "Sorting input dataset resulted in error: %s: %s" %( error_code2, stdout ))
skipped_lines = 0
first_invalid_line = 0
invalid_line = None
elems = []
j=0
for i, line in enumerate( file(infile1) ):
line = line.strip('\r\n')
if line and (not line.startswith( '#' )) and line != '':
j+=1
try:
elems = line.split('\t')
#if the start and/or end columns are not numbers, skip that line.
chr = elems[chr_col_1]
start = int(elems[start_col_1])
end = int(elems[end_col_1])
if strand_col_1 != -1: #Strand column is defined
strand = elems[strand_col_1]
#if the stand value is not + or -, skip that line.
assert strand in ['+', '-']
if direction == 'Upstream' or direction == 'Both':
if strand == '+':
for fline in file( tmpfile1.name ):
fline = fline.strip('\r\n')
if fline and not fline.startswith( '#' ):
try:
felems = fline.split('\t')
fchr = felems[chr_col_2]
if fchr != chr:
continue
if strand_col_1 != -1: #Strand column is defined
try:
fstrand = felems[strand_col_2]
except:
fstrand = ""
else:
fstrand = ""
if fstrand != "":
if strand != fstrand:
continue
fstart = int(felems[start_col_2])
fend = int(felems[end_col_2])
if fend < start: #Highest feature end value encountered i.e. the closest upstream feature found
print >>fo, "%s\t%s" %(line, fline)
break
except:
continue
elif strand == '-':
for fline in file( tmpfile2.name ):
fline = fline.strip('\r\n')
if fline and not fline.startswith( '#' ):
try:
felems = fline.split('\t')
fchr = felems[chr_col_2]
if fchr != chr:
continue
if strand_col_1 != -1: #Strand column is defined
try:
fstrand = felems[strand_col_2]
except:
fstrand = ""
else:
fstrand = ""
if fstrand != "":
if strand != fstrand:
continue
fstart = int(felems[start_col_2])
fend = int(felems[end_col_2])
if fstart > end: #Lowest feature start value encountered i.e. the closest upstream feature found
print >>fo, "%s\t%s" %(line, fline)
break
except:
continue
if direction == 'Downstream' or direction == 'Both':
if strand == '-':
for fline in file( tmpfile1.name ):
fline = fline.strip('\r\n')
if fline and not fline.startswith( '#' ):
try:
felems = fline.split('\t')
fchr = felems[chr_col_2]
if fchr != chr:
continue
if strand_col_1 != -1: #Strand column is defined
try:
fstrand = felems[strand_col_2]
except:
fstrand = ""
else:
fstrand = ""
if fstrand != "":
if strand != fstrand:
continue
fstart = int(felems[start_col_2])
fend = int(felems[end_col_2])
if fend < start: #Highest feature end value encountered i.e. the closest DOWNstream feature found
print >>fo, "%s\t%s" %(line, fline)
break
except:
continue
elif strand == '+':
for fline in file( tmpfile2.name ):
fline = fline.strip('\r\n')
if fline and not fline.startswith( '#' ):
try:
felems = fline.split('\t')
fchr = felems[chr_col_2]
if fchr != chr:
continue
if strand_col_1 != -1: #Strand column is defined
try:
fstrand = felems[strand_col_2]
except:
fstrand = ""
else:
fstrand = ""
if fstrand != "":
if strand != fstrand:
continue
fstart = int(felems[start_col_2])
fend = int(felems[end_col_2])
if fstart > end: #Lowest feature start value encountered i.e. the closest DOWNstream feature found
print >>fo, "%s\t%s" %(line, fline)
break
except:
continue
except Exception, exo:
skipped_lines += 1
if not invalid_line:
first_invalid_line = i + 1
invalid_line = line
fo.close()
fi.close()
#If number of skipped lines = num of lines in the file, inform the user to check metadata attributes of the input file.
if skipped_lines == j:
print 'Data issue: Skipped all lines in your input. Check the metadata attributes of the chosen input by clicking on the pencil icon next to it.'
sys.exit()
elif skipped_lines > 0:
print '(Data issue: skipped %d invalid lines starting at line #%d which is "%s")' % ( skipped_lines, first_invalid_line, invalid_line )
print 'Location : %s' %(direction)
if __name__ == "__main__":
main()
tools/new_operations/flanking_features.xml
+69
View File
@@ -0,0 +1,69 @@
<tool id="flanking_features1" name="Fetch closest feature">
<description> for every interval</description>
<command interpreter="python2.4">flanking_features.py $input1 $input2 $out_file1 $direction -1 $input1_chromCol,$input1_startCol,$input1_endCol,$input1_strandCol -2 $input2_chromCol,$input2_startCol,$input2_endCol,$input2_strandCol</command>
<inputs>
<param format="interval" name="input1" type="data" label="Select data containing primary intervals"/>
<param format="interval" name="input2" type="data" label="Select data containing features"/>
<param name="direction" type="select" label="Location of the flanking features">
<option value="Both">Both Upstream and Downstream</option>
<option value="Upstream">Upstream</option>
<option value="Downstream">Downstream</option>
</param>
</inputs>
<outputs>
<data format="interval" name="out_file1" />
</outputs>
<tests>
<test>
<param name="input1" value="flanks_inp.bed"/>
<param name="input2" value="11.bed"/>
<param name="direction" value="Both"/>
<output name="out_file1" file="closest_features.interval"/>
</test>
</tests>
<help>
.. class:: infomark
**What it does**
For every interval in the primary input file (input 1), this tool fetches the **closest** upstream and/or downstream features from feature file (input 2).
-----
.. class:: warningmark
**Note:**
Every line should contain at least 3 columns: Chromosome number, Start and Stop co-ordinates. If any of these columns is missing or if start and stop co-ordinates are not numerical, the tool may encounter exceptions and such lines are skipped as invalid. The number of invalid skipped lines is documented in the resulting history item as a "Data issue".
-----
**Example**
If the **primary intervals** are as follows::
chr1 10 100 Query1.1
chr1 500 1000 Query1.2
chr1 1100 1250 Query1.3
and the **features** are as follows::
chr1 120 180 Query2.1
chr1 140 200 Query2.2
chr1 580 1050 Query2.3
chr1 2000 2204 Query2.4
chr1 2500 3000 Query2.5
Running this tool with direction as **Both Upstream and Downstream** will return the following::
chr1 10 100 Query1.1 chr1 120 180 Query2.1
chr1 500 1000 Query1.2 chr1 140 200 Query2.2
chr1 500 1000 Query1.2 chr1 2000 2204 Query2.4
chr1 1100 1250 Query1.3 chr1 580 1050 Query2.3
chr1 1100 1250 Query1.3 chr1 2000 2204 Query2.4
</help>
</tool>