diff --git a/tool_conf.xml.sample b/tool_conf.xml.sample index 29f7fdbf81b..7ebfbf8516a 100644 --- a/tool_conf.xml.sample +++ b/tool_conf.xml.sample @@ -103,6 +103,7 @@ +
diff --git a/tools/new_operations/flanking_features.py b/tools/new_operations/flanking_features.py new file mode 100644 index 00000000000..a94facf6a11 --- /dev/null +++ b/tools/new_operations/flanking_features.py @@ -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() \ No newline at end of file diff --git a/tools/new_operations/flanking_features.xml b/tools/new_operations/flanking_features.xml new file mode 100644 index 00000000000..40e6d54e37b --- /dev/null +++ b/tools/new_operations/flanking_features.xml @@ -0,0 +1,69 @@ + + for every interval + 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 + + + + + + + + + + + + + + + + + + + + + + +.. 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 + + + + + \ No newline at end of file