].
- Intervals are GenomicInterval objects.
- """
- primary = readers[0]
- features = readers[1]
- either = False
- if region == 'Upstream':
- up, down = True, False
- elif region == 'Downstream':
- up, down = False, True
- else:
- up, down = True, True
- if region == 'Either':
- either = True
-
- # Read features into memory:
- rightTree = quicksect.IntervalTree()
- for item in features:
- if type( item ) is GenomicInterval:
- rightTree.insert( item, features.linenum, item )
-
- for interval in primary:
- if type( interval ) is Header:
- yield interval
- if type( interval ) is Comment and comments:
- yield interval
- elif type( interval ) == GenomicInterval:
- chrom = interval.chrom
- start = int(interval.start)
- end = int(interval.end)
- strand = interval.strand
- if chrom not in rightTree.chroms:
- continue
- else:
- root = rightTree.chroms[chrom] #root node for the chrom tree
- result_up = []
- result_down = []
- if (strand == '+' and up) or (strand == '-' and down):
- #upstream +ve strand and downstream -ve strand cases
- get_closest_feature (root, 1, start, None, lambda node: result_up.append( node ), None)
-
- if (strand == '+' and down) or (strand == '-' and up):
- #downstream +ve strand and upstream -ve strand case
- get_closest_feature (root, 0, None, end-1, None, lambda node: result_down.append( node ))
-
- if result_up:
- if len(result_up) > 1: #The results_up list has a list of intervals upstream to the given interval.
- ends = []
- for n in result_up:
- ends.append(n.end)
- res_ind = ends.index(max(ends)) #fetch the index of the closest interval i.e. the interval with the max end from the results_up list
- else:
- res_ind = 0
- if not(either):
- yield [ interval, result_up[res_ind].other ]
-
- if result_down:
- if not(either):
- #The last element of result_down will be the closest element to the given interval
- yield [ interval, result_down[-1].other ]
-
- if either and (result_up or result_down):
- iter_val = []
- if result_up and result_down:
- if abs(start - int(result_up[res_ind].end)) <= abs(end - int(result_down[-1].start)):
- iter_val = [ interval, result_up[res_ind].other ]
- else:
- #The last element of result_down will be the closest element to the given interval
- iter_val = [ interval, result_down[-1].other ]
- elif result_up:
- iter_val = [ interval, result_up[res_ind].other ]
- elif result_down:
- #The last element of result_down will be the closest element to the given interval
- iter_val = [ interval, result_down[-1].other ]
- yield iter_val
-
-def main():
- 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 )
- in1_gff_format = bool( options.gff1 )
- in2_gff_format = bool( options.gff2 )
- in_fname, in2_fname, out_fname, direction = args
- except:
- doc_optparse.exception()
-
- # Set readers to handle either GFF or default format.
- if in1_gff_format:
- in1_reader_wrapper = GFFIntervalToBEDReaderWrapper
- else:
- in1_reader_wrapper = NiceReaderWrapper
- if in2_gff_format:
- in2_reader_wrapper = GFFIntervalToBEDReaderWrapper
- else:
- in2_reader_wrapper = NiceReaderWrapper
-
- g1 = in1_reader_wrapper( fileinput.FileInput( in_fname ),
- chrom_col=chr_col_1,
- start_col=start_col_1,
- end_col=end_col_1,
- strand_col=strand_col_1,
- fix_strand=True )
- g2 = in2_reader_wrapper( fileinput.FileInput( in2_fname ),
- chrom_col=chr_col_2,
- start_col=start_col_2,
- end_col=end_col_2,
- strand_col=strand_col_2,
- fix_strand=True )
-
- # Find flanking features.
- out_file = open( out_fname, "w" )
- try:
- for result in proximal_region_finder([g1,g2], direction):
- if type( result ) is list:
- line, closest_feature = result
- # Need to join outputs differently depending on file types.
- if in1_gff_format:
- # Output is GFF with added attribute 'closest feature.'
-
- # Invervals are in BED coordinates; need to convert to GFF.
- line = convert_bed_coords_to_gff( line )
- closest_feature = convert_bed_coords_to_gff( closest_feature )
-
- # Replace double quotes with single quotes in closest feature's attributes.
- out_file.write( "%s closest_feature \"%s\" \n" %
- ( "\t".join( line.fields ), \
- "\t".join( closest_feature.fields ).replace( "\"", "\\\"" )
- ) )
- else:
- # Output is BED + closest feature fields.
- output_line_fields = []
- output_line_fields.extend( line.fields )
- output_line_fields.extend( closest_feature.fields )
- out_file.write( "%s\n" % ( "\t".join( output_line_fields ) ) )
- else:
- out_file.write( "%s\n" % result )
- except ParseError, exc:
- fail( "Invalid file format: %s" % str( exc ) )
-
- print "Direction: %s" %(direction)
- if g1.skipped > 0:
- print skipped( g1, filedesc=" of 1st dataset" )
- if g2.skipped > 0:
- print skipped( g2, filedesc=" of 2nd dataset" )
-
-if __name__ == "__main__":
- main()
diff --git a/tools/new_operations/flanking_features.xml b/tools/new_operations/flanking_features.xml
deleted file mode 100644
index e8018b83677..00000000000
--- a/tools/new_operations/flanking_features.xml
+++ /dev/null
@@ -1,127 +0,0 @@
-
- for every interval
-
- flanking_features.py $input1 $input2 $out_file1 $direction
-
- #if isinstance( $input1.datatype, $__app__.datatypes_registry.get_datatype_by_extension('gff').__class__):
- -1 1,4,5,7 --gff1
- #else:
- -1 ${input1.metadata.chromCol},${input1.metadata.startCol},${input1.metadata.endCol},${input1.metadata.strandCol}
- #end if
-
- #if isinstance( $input2.datatype, $__app__.datatypes_registry.get_datatype_by_extension('gff').__class__):
- -2 1,4,5,7 --gff2
- #else:
- -2 ${input2.metadata.chromCol},${input2.metadata.startCol},${input2.metadata.endCol},${input2.metadata.strandCol}
- #end if
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-.. class:: infomark
-
-**What it does**
-
-For every interval in the **interval** dataset, this tool fetches the **closest non-overlapping** upstream and / or downstream features from the **features** dataset.
-
------
-
-.. class:: warningmark
-
-**Note:**
-
-Every line should contain at least 3 columns: chromosome number, start and stop coordinates. If any of these columns is missing or if start and stop coordinates are not numerical, the lines will be treated as invalid and skipped. The number of skipped lines is documented in the resulting history item as a "data issue".
-
-If the strand column is missing from your input interval dataset, the intervals will be considered to be on positive strand. You can add a strand column to your input dataset by using the *Text Manipulation->Add column* tool.
-
-For GFF files, features are added as a GTF-style attribute at the end of the line.
-
------
-
-**Example**
-
-If the **intervals** are::
-
- chr1 10 100 Query1.1
- chr1 500 1000 Query1.2
- chr1 1100 1250 Query1.3
-
-and the **features** are::
-
- 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 for **Both Upstream and Downstream** will return::
-
- 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
diff --git a/tools/new_operations/get_flanks.py b/tools/new_operations/get_flanks.py
deleted file mode 100644
index 05d6ef42605..00000000000
--- a/tools/new_operations/get_flanks.py
+++ /dev/null
@@ -1,191 +0,0 @@
-#!/usr/bin/env python
-#Done by: Guru
-
-"""
-Get Flanking regions.
-
-usage: %prog input out_file size direction region
- -l, --cols=N,N,N,N: Columns for chrom, start, end, strand in file
- -o, --off=N: Offset
-"""
-
-import sys, re, os
-from galaxy import eggs
-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():
- try:
- if int( sys.argv[3] ) < 0:
- raise Exception
- except:
- stop_err( "Length of flanking region(s) must be a non-negative integer." )
-
- # 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.cols )
- inp_file, out_file, size, direction, region = args
- if strand_col_1 <= 0:
- strand = "+" #if strand is not defined, default it to +
- except:
- stop_err( "Metadata issue, correct the metadata attributes by clicking on the pencil icon in the history item." )
- try:
- offset = int(options.off)
- size = int(size)
- except:
- stop_err( "Invalid offset or length entered. Try again by entering valid integer values." )
-
- fo = open(out_file,'w')
-
- skipped_lines = 0
- first_invalid_line = 0
- invalid_line = None
- elems = []
- j=0
- for i, line in enumerate( file( inp_file ) ):
- line = line.strip()
- 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.
- assert int(elems[start_col_1])
- assert int(elems[end_col_1])
- if strand_col_1 != -1:
- strand = elems[strand_col_1]
- #if the stand value is not + or -, skip that line.
- assert strand in ['+', '-']
- if direction == 'Upstream':
- if strand == '+':
- if region == 'end':
- elems[end_col_1] = str(int(elems[end_col_1]) + offset)
- elems[start_col_1] = str( int(elems[end_col_1]) - size )
- else:
- elems[end_col_1] = str(int(elems[start_col_1]) + offset)
- elems[start_col_1] = str( int(elems[end_col_1]) - size )
- elif strand == '-':
- if region == 'end':
- elems[start_col_1] = str(int(elems[start_col_1]) - offset)
- elems[end_col_1] = str(int(elems[start_col_1]) + size)
- else:
- elems[start_col_1] = str(int(elems[end_col_1]) - offset)
- elems[end_col_1] = str(int(elems[start_col_1]) + size)
- assert int(elems[start_col_1]) > 0 and int(elems[end_col_1]) > 0
- fo.write( "%s\n" % '\t'.join( elems ) )
-
- elif direction == 'Downstream':
- if strand == '-':
- if region == 'start':
- elems[end_col_1] = str(int(elems[end_col_1]) - offset)
- elems[start_col_1] = str( int(elems[end_col_1]) - size )
- else:
- elems[end_col_1] = str(int(elems[start_col_1]) - offset)
- elems[start_col_1] = str( int(elems[end_col_1]) - size )
- elif strand == '+':
- if region == 'start':
- elems[start_col_1] = str(int(elems[start_col_1]) + offset)
- elems[end_col_1] = str(int(elems[start_col_1]) + size)
- else:
- elems[start_col_1] = str(int(elems[end_col_1]) + offset)
- elems[end_col_1] = str(int(elems[start_col_1]) + size)
- assert int(elems[start_col_1]) > 0 and int(elems[end_col_1]) > 0
- fo.write( "%s\n" % '\t'.join( elems ) )
-
- elif direction == 'Both':
- if strand == '-':
- if region == 'start':
- start = str(int(elems[end_col_1]) - offset)
- end1 = str(int(start) + size)
- end2 = str(int(start) - size)
- elems[start_col_1]=start
- elems[end_col_1]=end1
- assert int(elems[start_col_1]) > 0 and int(elems[end_col_1]) > 0
- fo.write( "%s\n" % '\t'.join( elems ) )
- elems[start_col_1]=end2
- elems[end_col_1]=start
- assert int(elems[start_col_1]) > 0 and int(elems[end_col_1]) > 0
- fo.write( "%s\n" % '\t'.join( elems ) )
- elif region == 'end':
- start = str(int(elems[start_col_1]) - offset)
- end1 = str(int(start) + size)
- end2 = str(int(start) - size)
- elems[start_col_1]=start
- elems[end_col_1]=end1
- assert int(elems[start_col_1]) > 0 and int(elems[end_col_1]) > 0
- fo.write( "%s\n" % '\t'.join( elems ) )
- elems[start_col_1]=end2
- elems[end_col_1]=start
- assert int(elems[start_col_1]) > 0 and int(elems[end_col_1]) > 0
- fo.write( "%s\n" % '\t'.join( elems ) )
- else:
- start1 = str(int(elems[end_col_1]) - offset)
- end1 = str(int(start1) + size)
- start2 = str(int(elems[start_col_1]) - offset)
- end2 = str(int(start2) - size)
- elems[start_col_1]=start1
- elems[end_col_1]=end1
- assert int(elems[start_col_1]) > 0 and int(elems[end_col_1]) > 0
- fo.write( "%s\n" % '\t'.join( elems ) )
- elems[start_col_1]=end2
- elems[end_col_1]=start2
- assert int(elems[start_col_1]) > 0 and int(elems[end_col_1]) > 0
- fo.write( "%s\n" % '\t'.join( elems ) )
- elif strand == '+':
- if region == 'start':
- start = str(int(elems[start_col_1]) + offset)
- end1 = str(int(start) - size)
- end2 = str(int(start) + size)
- elems[start_col_1]=end1
- elems[end_col_1]=start
- assert int(elems[start_col_1]) > 0 and int(elems[end_col_1]) > 0
- fo.write( "%s\n" % '\t'.join( elems ) )
- elems[start_col_1]=start
- elems[end_col_1]=end2
- assert int(elems[start_col_1]) > 0 and int(elems[end_col_1]) > 0
- fo.write( "%s\n" % '\t'.join( elems ) )
- elif region == 'end':
- start = str(int(elems[end_col_1]) + offset)
- end1 = str(int(start) - size)
- end2 = str(int(start) + size)
- elems[start_col_1]=end1
- elems[end_col_1]=start
- assert int(elems[start_col_1]) > 0 and int(elems[end_col_1]) > 0
- fo.write( "%s\n" % '\t'.join( elems ) )
- elems[start_col_1]=start
- elems[end_col_1]=end2
- assert int(elems[start_col_1]) > 0 and int(elems[end_col_1]) > 0
- fo.write( "%s\n" % '\t'.join( elems ) )
- else:
- start1 = str(int(elems[start_col_1]) + offset)
- end1 = str(int(start1) - size)
- start2 = str(int(elems[end_col_1]) + offset)
- end2 = str(int(start2) + size)
- elems[start_col_1]=end1
- elems[end_col_1]=start1
- assert int(elems[start_col_1]) > 0 and int(elems[end_col_1]) > 0
- fo.write( "%s\n" % '\t'.join( elems ) )
- elems[start_col_1]=start2
- elems[end_col_1]=end2
- assert int(elems[start_col_1]) > 0 and int(elems[end_col_1]) > 0
- fo.write( "%s\n" % '\t'.join( elems ) )
- except:
- skipped_lines += 1
- if not invalid_line:
- first_invalid_line = i + 1
- invalid_line = line
- fo.close()
-
- if skipped_lines == j:
- stop_err( "Data issue: click the pencil icon in the history item to correct the metadata attributes." )
- if skipped_lines > 0:
- print 'Skipped %d invalid lines starting with #%dL "%s"' % ( skipped_lines, first_invalid_line, invalid_line )
- print 'Location: %s, Region: %s, Flank-length: %d, Offset: %d ' %( direction, region, size, offset )
-
-if __name__ == "__main__":
- main()
diff --git a/tools/new_operations/get_flanks.xml b/tools/new_operations/get_flanks.xml
deleted file mode 100644
index 2c4a0e5b852..00000000000
--- a/tools/new_operations/get_flanks.xml
+++ /dev/null
@@ -1,78 +0,0 @@
-
- returns flanking region/s for every gene
- get_flanks.py $input $out_file1 $size $direction $region -o $offset -l ${input.metadata.chromCol},${input.metadata.startCol},${input.metadata.endCol},${input.metadata.strandCol}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-This tool finds the upstream and/or downstream flanking region(s) of all the selected regions in the input file.
-
-**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 1**
-
-- For the following dataset::
-
- chr22 1000 7000 NM_174568 0 +
-
-- running get flanks with Region: Around start, Offset: -200, Flank-length: 300 and Location: Upstream will return **(Red: Dataset positive strand; Blue: Flanks output)**::
-
- chr22 500 800 NM_174568 0 +
-
-.. image:: ${static_path}/operation_icons/flanks_ex1.gif
-
-**Example 2**
-
-- For the following dataset::
-
- chr22 1000 7000 NM_028946 0 -
-
-- running get flanks with Region: Whole, Offset: 200, Flank-length: 300 and Location: Downstream will return **(Orange: Dataset negative strand; Magenta: Flanks output)**::
-
- chr22 500 800 NM_028946 0 -
-
-.. image:: ${static_path}/operation_icons/flanks_ex2.gif
-
-
-
-
-
diff --git a/tools/new_operations/gops_basecoverage.py b/tools/new_operations/gops_basecoverage.py
deleted file mode 100644
index da78701afa9..00000000000
--- a/tools/new_operations/gops_basecoverage.py
+++ /dev/null
@@ -1,50 +0,0 @@
-#!/usr/bin/env python
-"""
-Count total base coverage.
-
-usage: %prog in_file out_file
- -1, --cols1=N,N,N,N: Columns for start, end, strand in first file
-"""
-from galaxy import eggs
-import pkg_resources
-pkg_resources.require( "bx-python" )
-import sys, traceback, fileinput
-from warnings import warn
-from bx.intervals import *
-from bx.intervals.io import *
-from bx.intervals.operations.base_coverage import *
-from bx.cookbook import doc_optparse
-from galaxy.tools.util.galaxyops import *
-
-assert sys.version_info[:2] >= ( 2, 4 )
-
-def main():
- upstream_pad = 0
- downstream_pad = 0
-
- options, args = doc_optparse.parse( __doc__ )
- try:
- chr_col_1, start_col_1, end_col_1, strand_col_1 = parse_cols_arg( options.cols1 )
- in_fname, out_fname = args
- except:
- doc_optparse.exception()
-
- g1 = NiceReaderWrapper( fileinput.FileInput( in_fname ),
- chrom_col=chr_col_1,
- start_col=start_col_1,
- end_col=end_col_1,
- strand_col = strand_col_1,
- fix_strand=True )
-
- try:
- bases = base_coverage(g1)
- except ParseError, exc:
- fail( "Invalid file format: %s" % str( exc ) )
- out_file = open( out_fname, "w" )
- out_file.write( "%s\n" % str( bases ) )
- out_file.close()
- if g1.skipped > 0:
- print skipped( g1, filedesc="" )
-
-if __name__ == "__main__":
- main()
diff --git a/tools/new_operations/gops_cluster.py b/tools/new_operations/gops_cluster.py
deleted file mode 100644
index ac2f4f75d0a..00000000000
--- a/tools/new_operations/gops_cluster.py
+++ /dev/null
@@ -1,132 +0,0 @@
-#!/usr/bin/env python
-"""
-Cluster regions of intervals.
-
-usage: %prog in_file out_file
- -1, --cols1=N,N,N,N: Columns for start, end, strand in file
- -d, --distance=N: Maximum distance between clustered intervals
- -v, --overlap=N: Minimum overlap require (negative distance)
- -m, --minregions=N: Minimum regions per cluster
- -o, --output=N: 1)merged 2)filtered 3)clustered 4) minimum 5) maximum
-"""
-from galaxy import eggs
-import pkg_resources
-pkg_resources.require( "bx-python" )
-import sys, traceback, fileinput
-from warnings import warn
-from bx.intervals import *
-from bx.intervals.io import *
-from bx.intervals.operations.find_clusters import *
-from bx.cookbook import doc_optparse
-from galaxy.tools.util.galaxyops import *
-
-assert sys.version_info[:2] >= ( 2, 4 )
-
-def main():
- distance = 0
- minregions = 2
- output = 1
- upstream_pad = 0
- downstream_pad = 0
-
- options, args = doc_optparse.parse( __doc__ )
- try:
- chr_col_1, start_col_1, end_col_1, strand_col_1 = parse_cols_arg( options.cols1 )
- if options.distance: distance = int( options.distance )
- if options.overlap: distance = -1 * int( options.overlap )
- if options.output: output = int( options.output )
- if options.minregions: minregions = int( options.minregions )
- in_fname, out_fname = args
- except:
- doc_optparse.exception()
-
- g1 = NiceReaderWrapper( fileinput.FileInput( in_fname ),
- chrom_col=chr_col_1,
- start_col=start_col_1,
- end_col=end_col_1,
- strand_col=strand_col_1,
- fix_strand=True )
-
- # Get the cluster tree
- try:
- clusters, extra = find_clusters( g1, mincols=distance, minregions=minregions)
- except ParseError, exc:
- fail( "Invalid file format: %s" % str( exc ) )
-
- f1 = open( in_fname, "r" )
- out_file = open( out_fname, "w" )
-
- # If "merge"
- if output == 1:
- fields = ["." for x in range(max(g1.chrom_col, g1.start_col, g1.end_col)+1)]
- for chrom, tree in clusters.items():
- for start, end, lines in tree.getregions():
- fields[g1.chrom_col] = chrom
- fields[g1.start_col] = str(start)
- fields[g1.end_col] = str(end)
- out_file.write( "%s\n" % "\t".join( fields ) )
-
- # If "filtered" we preserve order of file and comments, etc.
- if output == 2:
- linenums = dict()
- for chrom, tree in clusters.items():
- for linenum in tree.getlines():
- linenums[linenum] = 0
- linenum = -1
- f1.seek(0)
- for line in f1.readlines():
- linenum += 1
- if linenum in linenums or linenum in extra:
- out_file.write( "%s\n" % line.rstrip( "\n\r" ) )
-
- # If "clustered" we output original intervals, but near each other (i.e. clustered)
- if output == 3:
- linenums = list()
- f1.seek(0)
- fileLines = f1.readlines()
- for chrom, tree in clusters.items():
- for linenum in tree.getlines():
- out_file.write( "%s\n" % fileLines[linenum].rstrip( "\n\r" ) )
-
- # If "minimum" we output the smallest interval in each cluster
- if output == 4 or output == 5:
- linenums = list()
- f1.seek(0)
- fileLines = f1.readlines()
- for chrom, tree in clusters.items():
- regions = tree.getregions()
- for start, end, lines in tree.getregions():
- outsize = -1
- outinterval = None
- for line in lines:
- # three nested for loops?
- # should only execute this code once per line
- fileline = fileLines[line].rstrip("\n\r")
- try:
- cluster_interval = GenomicInterval( g1, fileline.split("\t"),
- g1.chrom_col,
- g1.start_col,
- g1.end_col,
- g1.strand_col,
- g1.default_strand,
- g1.fix_strand )
- except Exception, exc:
- print >> sys.stderr, str( exc )
- f1.close()
- sys.exit()
- interval_size = cluster_interval.end - cluster_interval.start
- if outsize == -1 or \
- ( outsize > interval_size and output == 4 ) or \
- ( outsize < interval_size and output == 5 ) :
- outinterval = cluster_interval
- outsize = interval_size
- out_file.write( "%s\n" % outinterval )
-
- f1.close()
- out_file.close()
-
- if g1.skipped > 0:
- print skipped( g1, filedesc="" )
-
-if __name__ == "__main__":
- main()
diff --git a/tools/new_operations/gops_complement.py b/tools/new_operations/gops_complement.py
deleted file mode 100644
index 7615f83f4e2..00000000000
--- a/tools/new_operations/gops_complement.py
+++ /dev/null
@@ -1,98 +0,0 @@
-#!/usr/bin/env python
-"""
-Complement regions.
-
-usage: %prog in_file out_file
- -1, --cols1=N,N,N,N: Columns for chrom, start, end, strand in file
- -l, --lengths=N: Filename of .len file for species (chromosome lengths)
- -a, --all: Complement all chromosomes (Genome-wide complement)
-"""
-from galaxy import eggs
-import pkg_resources
-pkg_resources.require( "bx-python" )
-import sys, traceback, fileinput
-from warnings import warn
-from bx.intervals import *
-from bx.intervals.io import *
-from bx.intervals.operations.complement import complement
-from bx.intervals.operations.subtract import subtract
-from bx.cookbook import doc_optparse
-from galaxy.tools.util.galaxyops import *
-
-assert sys.version_info[:2] >= ( 2, 4 )
-
-def main():
- allchroms = False
- upstream_pad = 0
- downstream_pad = 0
-
- options, args = doc_optparse.parse( __doc__ )
- try:
- chr_col_1, start_col_1, end_col_1, strand_col_1 = parse_cols_arg( options.cols1 )
- lengths = options.lengths
- if options.all: allchroms = True
- in_fname, out_fname = args
- except:
- doc_optparse.exception()
-
- g1 = NiceReaderWrapper( fileinput.FileInput( in_fname ),
- chrom_col=chr_col_1,
- start_col=start_col_1,
- end_col=end_col_1,
- strand_col=strand_col_1,
- fix_strand=True )
-
- lens = dict()
- chroms = list()
- # dbfile is used to determine the length of each chromosome. The lengths
- # are added to the lens dict and passed copmlement operation code in bx.
- dbfile = fileinput.FileInput( lengths )
-
- if dbfile:
- if not allchroms:
- try:
- for line in dbfile:
- fields = line.split("\t")
- lens[fields[0]] = int(fields[1])
- except:
- # assume LEN doesn't exist or is corrupt somehow
- pass
- elif allchroms:
- try:
- for line in dbfile:
- fields = line.split("\t")
- end = int(fields[1])
- chroms.append("\t".join([fields[0],"0",str(end)]))
- except:
- pass
-
- # Safety...if the dbfile didn't exist and we're on allchroms, then
- # default to generic complement
- if allchroms and len(chroms) == 0:
- allchroms = False
-
- if allchroms:
- chromReader = GenomicIntervalReader(chroms)
- generator = subtract([chromReader, g1])
- else:
- generator = complement(g1, lens)
-
- out_file = open( out_fname, "w" )
-
- try:
- for interval in generator:
- if type( interval ) is GenomicInterval:
- out_file.write( "%s\n" % "\t".join( interval ) )
- else:
- out_file.write( "%s\n" % interval )
- except ParseError, exc:
- out_file.close()
- fail( "Invalid file format: %s" % str( exc ) )
-
- out_file.close()
-
- if g1.skipped > 0:
- print skipped( g1, filedesc="" )
-
-if __name__ == "__main__":
- main()
diff --git a/tools/new_operations/gops_concat.py b/tools/new_operations/gops_concat.py
deleted file mode 100644
index 9490c41401d..00000000000
--- a/tools/new_operations/gops_concat.py
+++ /dev/null
@@ -1,76 +0,0 @@
-#!/usr/bin/env python
-"""
-Concatenate two bed files. The concatenated files are returned in the
-same format as the first. If --sameformat is specified, then all
-columns will be treated as the same, and all fields will be saved,
-although the output will be trimmed to match the primary input. In
-addition, if --sameformat is specified, missing fields will be padded
-with a period(.).
-
-usage: %prog in_file_1 in_file_2 out_file
- -1, --cols1=N,N,N,N: Columns for chrom, start, end, strand in first file
- -2, --cols2=N,N,N,N: Columns for chrom, start, end, strand in second file
- -s, --sameformat: All files are precisely the same format.
-"""
-from galaxy import eggs
-import pkg_resources
-pkg_resources.require( "bx-python" )
-import sys, traceback, fileinput
-from warnings import warn
-from bx.intervals import *
-from bx.intervals.io import *
-from bx.intervals.operations.concat import *
-from bx.cookbook import doc_optparse
-from galaxy.tools.util.galaxyops import *
-
-assert sys.version_info[:2] >= ( 2, 4 )
-
-def main():
- sameformat=False
- upstream_pad = 0
- downstream_pad = 0
-
- 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 )
- if options.sameformat: sameformat = True
- in_file_1, in_file_2, out_fname = args
- except:
- doc_optparse.exception()
-
- g1 = NiceReaderWrapper( fileinput.FileInput( in_file_1 ),
- chrom_col=chr_col_1,
- start_col=start_col_1,
- end_col=end_col_1,
- strand_col=strand_col_1,
- fix_strand=True )
-
- g2 = NiceReaderWrapper( fileinput.FileInput( in_file_2 ),
- chrom_col=chr_col_2,
- start_col=start_col_2,
- end_col=end_col_2,
- strand_col=strand_col_2,
- fix_strand=True )
-
- out_file = open( out_fname, "w" )
-
- try:
- for line in concat( [g1, g2], sameformat=sameformat ):
- if type( line ) is GenomicInterval:
- out_file.write( "%s\n" % "\t".join( line.fields ) )
- else:
- out_file.write( "%s\n" % line )
- except ParseError, exc:
- out_file.close()
- fail( "Invalid file format: %s" % str( exc ) )
-
- out_file.close()
-
- if g1.skipped > 0:
- print skipped( g1, filedesc=" of 1st dataset" )
- if g2.skipped > 0:
- print skipped( g2, filedesc=" of 2nd dataset" )
-
-if __name__ == "__main__":
- main()
diff --git a/tools/new_operations/gops_coverage.py b/tools/new_operations/gops_coverage.py
deleted file mode 100644
index 91b3ad19803..00000000000
--- a/tools/new_operations/gops_coverage.py
+++ /dev/null
@@ -1,68 +0,0 @@
-#!/usr/bin/env python
-"""
-Calculate coverage of one query on another, and append the coverage to
-the last two columns as bases covered and percent coverage.
-
-usage: %prog bed_file_1 bed_file_2 out_file
- -1, --cols1=N,N,N,N: Columns for start, end, strand in first file
- -2, --cols2=N,N,N,N: Columns for start, end, strand in second file
-"""
-from galaxy import eggs
-import pkg_resources
-pkg_resources.require( "bx-python" )
-import sys, traceback, fileinput
-from warnings import warn
-from bx.intervals import *
-from bx.intervals.io import *
-from bx.intervals.operations.coverage import *
-from bx.cookbook import doc_optparse
-from galaxy.tools.util.galaxyops import *
-
-assert sys.version_info[:2] >= ( 2, 4 )
-
-def main():
- upstream_pad = 0
- downstream_pad = 0
-
- 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 )
- in_fname, in2_fname, out_fname = args
- except:
- doc_optparse.exception()
-
- g1 = NiceReaderWrapper( fileinput.FileInput( in_fname ),
- chrom_col=chr_col_1,
- start_col=start_col_1,
- end_col=end_col_1,
- strand_col=strand_col_1,
- fix_strand=True )
- g2 = NiceReaderWrapper( fileinput.FileInput( in2_fname ),
- chrom_col=chr_col_2,
- start_col=start_col_2,
- end_col=end_col_2,
- strand_col=strand_col_2,
- fix_strand=True )
-
- out_file = open( out_fname, "w" )
-
- try:
- for line in coverage( [g1,g2] ):
- if type( line ) is GenomicInterval:
- out_file.write( "%s\n" % "\t".join( line.fields ) )
- else:
- out_file.write( "%s\n" % line )
- except ParseError, exc:
- out_file.close()
- fail( "Invalid file format: %s" % str( exc ) )
-
- out_file.close()
-
- if g1.skipped > 0:
- print skipped( g1, filedesc=" of 1st dataset" )
- if g2.skipped > 0:
- print skipped( g2, filedesc=" of 2nd dataset" )
-
-if __name__ == "__main__":
- main()
diff --git a/tools/new_operations/gops_intersect.py b/tools/new_operations/gops_intersect.py
deleted file mode 100755
index 11a0c80f902..00000000000
--- a/tools/new_operations/gops_intersect.py
+++ /dev/null
@@ -1,98 +0,0 @@
-#!/usr/bin/env python
-"""
-Find regions of first interval file that overlap regions in a second interval file.
-Interval files can either be BED or GFF format.
-
-usage: %prog interval_file_1 interval_file_2 out_file
- -1, --cols1=N,N,N,N: Columns for start, end, strand in first file
- -2, --cols2=N,N,N,N: Columns for start, end, strand in second file
- -m, --mincols=N: Require this much overlap (default 1bp)
- -p, --pieces: just print pieces of second set (after padding)
- -G, --gff1: input 1 is GFF format, meaning start and end coordinates are 1-based, closed interval
- -H, --gff2: input 2 is GFF format, meaning start and end coordinates are 1-based, closed interval
-"""
-from galaxy import eggs
-import pkg_resources
-pkg_resources.require( "bx-python" )
-import sys, traceback, fileinput
-from warnings import warn
-from bx.intervals import *
-from bx.intervals.io import *
-from bx.intervals.operations.intersect import *
-from bx.cookbook import doc_optparse
-from galaxy.tools.util.galaxyops import *
-from galaxy.datatypes.util.gff_util import GFFFeature, GFFReaderWrapper, convert_bed_coords_to_gff
-
-assert sys.version_info[:2] >= ( 2, 4 )
-
-def main():
- mincols = 1
- upstream_pad = 0
- downstream_pad = 0
-
- 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 )
- if options.mincols: mincols = int( options.mincols )
- pieces = bool( options.pieces )
- in1_gff_format = bool( options.gff1 )
- in2_gff_format = bool( options.gff2 )
- in_fname, in2_fname, out_fname = args
- except:
- doc_optparse.exception()
-
- # Set readers to handle either GFF or default format.
- if in1_gff_format:
- in1_reader_wrapper = GFFReaderWrapper
- else:
- in1_reader_wrapper = NiceReaderWrapper
- if in2_gff_format:
- in2_reader_wrapper = GFFReaderWrapper
- else:
- in2_reader_wrapper = NiceReaderWrapper
-
- g1 = in1_reader_wrapper( fileinput.FileInput( in_fname ),
- chrom_col=chr_col_1,
- start_col=start_col_1,
- end_col=end_col_1,
- strand_col=strand_col_1,
- fix_strand=True )
- if in1_gff_format:
- # Intersect requires coordinates in BED format.
- g1.convert_to_bed_coord=True
- g2 = in2_reader_wrapper( fileinput.FileInput( in2_fname ),
- chrom_col=chr_col_2,
- start_col=start_col_2,
- end_col=end_col_2,
- strand_col=strand_col_2,
- fix_strand=True )
- if in2_gff_format:
- # Intersect requires coordinates in BED format.
- g2.convert_to_bed_coord=True
-
- out_file = open( out_fname, "w" )
- try:
- for feature in intersect( [g1,g2], pieces=pieces, mincols=mincols ):
- if isinstance( feature, GFFFeature ):
- # Convert back to GFF coordinates since reader converted automatically.
- convert_bed_coords_to_gff( feature )
- for interval in feature.intervals:
- out_file.write( "%s\n" % "\t".join( interval.fields ) )
- elif isinstance( feature, GenomicInterval ):
- out_file.write( "%s\n" % "\t".join( feature.fields ) )
- else:
- out_file.write( "%s\n" % feature )
- except ParseError, e:
- out_file.close()
- fail( "Invalid file format: %s" % str( e ) )
-
- out_file.close()
-
- if g1.skipped > 0:
- print skipped( g1, filedesc=" of 1st dataset" )
- if g2.skipped > 0:
- print skipped( g2, filedesc=" of 2nd dataset" )
-
-if __name__ == "__main__":
- main()
diff --git a/tools/new_operations/gops_join.py b/tools/new_operations/gops_join.py
deleted file mode 100644
index daab6cb18e0..00000000000
--- a/tools/new_operations/gops_join.py
+++ /dev/null
@@ -1,82 +0,0 @@
-#!/usr/bin/env python
-"""
-Join two sets of intervals using their overlap as the key.
-
-usage: %prog bed_file_1 bed_file_2 out_file
- -1, --cols1=N,N,N,N: Columns for start, end, strand in first file
- -2, --cols2=N,N,N,N: Columns for start, end, strand in second file
- -m, --mincols=N: Require this much overlap (default 1bp)
- -f, --fill=N: none, right, left, both
-"""
-from galaxy import eggs
-import pkg_resources
-pkg_resources.require( "bx-python" )
-import sys, traceback, fileinput
-from warnings import warn
-from bx.intervals import *
-from bx.intervals.io import *
-from bx.intervals.operations.join import *
-from bx.cookbook import doc_optparse
-from galaxy.tools.util.galaxyops import *
-
-assert sys.version_info[:2] >= ( 2, 4 )
-
-def main():
- mincols = 1
- upstream_pad = 0
- downstream_pad = 0
- leftfill = False
- rightfill = False
-
- 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 )
- if options.mincols: mincols = int( options.mincols )
- if options.fill:
- if options.fill == "both":
- rightfill = leftfill = True
- else:
- rightfill = options.fill == "right"
- leftfill = options.fill == "left"
- in_fname, in2_fname, out_fname = args
- except:
- doc_optparse.exception()
-
- g1 = NiceReaderWrapper( fileinput.FileInput( in_fname ),
- chrom_col=chr_col_1,
- start_col=start_col_1,
- end_col=end_col_1,
- strand_col=strand_col_1,
- fix_strand=True )
- g2 = NiceReaderWrapper( fileinput.FileInput( in2_fname ),
- chrom_col=chr_col_2,
- start_col=start_col_2,
- end_col=end_col_2,
- strand_col=strand_col_2,
- fix_strand=True )
-
- out_file = open( out_fname, "w" )
-
- try:
- for outfields in join(g1, g2, mincols=mincols, rightfill=rightfill, leftfill=leftfill):
- if type( outfields ) is list:
- out_file.write( "%s\n" % "\t".join( outfields ) )
- else:
- out_file.write( "%s\n" % outfields )
- except ParseError, exc:
- out_file.close()
- fail( "Invalid file format: %s" % str( exc ) )
- except MemoryError:
- out_file.close()
- fail( "Input datasets were too large to complete the join operation." )
-
- out_file.close()
-
- if g1.skipped > 0:
- print skipped( g1, filedesc=" of 1st dataset" )
- if g2.skipped > 0:
- print skipped( g2, filedesc=" of 2nd dataset" )
-
-if __name__ == "__main__":
- main()
diff --git a/tools/new_operations/gops_merge.py b/tools/new_operations/gops_merge.py
deleted file mode 100644
index 85d215f83f4..00000000000
--- a/tools/new_operations/gops_merge.py
+++ /dev/null
@@ -1,71 +0,0 @@
-#!/usr/bin/env python
-"""
-Merge overlaping regions.
-
-usage: %prog in_file out_file
- -1, --cols1=N,N,N,N: Columns for start, end, strand in first file
- -m, --mincols=N: Require this much overlap (default 1bp)
- -3, --threecol: Output 3 column bed
-"""
-from galaxy import eggs
-import pkg_resources
-pkg_resources.require( "bx-python" )
-import sys, traceback, fileinput
-from warnings import warn
-from bx.intervals import *
-from bx.intervals.io import *
-from bx.intervals.operations.merge import *
-from bx.cookbook import doc_optparse
-from galaxy.tools.util.galaxyops import *
-
-assert sys.version_info[:2] >= ( 2, 4 )
-
-def main():
- mincols = 1
- upstream_pad = 0
- downstream_pad = 0
-
- options, args = doc_optparse.parse( __doc__ )
- try:
- chr_col_1, start_col_1, end_col_1, strand_col_1 = parse_cols_arg( options.cols1 )
- if options.mincols: mincols = int( options.mincols )
- in_fname, out_fname = args
- except:
- doc_optparse.exception()
-
- g1 = NiceReaderWrapper( fileinput.FileInput( in_fname ),
- chrom_col=chr_col_1,
- start_col=start_col_1,
- end_col=end_col_1,
- strand_col = strand_col_1,
- fix_strand=True )
-
- out_file = open( out_fname, "w" )
-
- try:
- for line in merge(g1,mincols=mincols):
- if options.threecol:
- if type( line ) is GenomicInterval:
- out_file.write( "%s\t%s\t%s\n" % ( line.chrom, str( line.startCol ), str( line.endCol ) ) )
- elif type( line ) is list:
- out_file.write( "%s\t%s\t%s\n" % ( line[chr_col_1], str( line[start_col_1] ), str( line[end_col_1] ) ) )
- else:
- out_file.write( "%s\n" % line )
- else:
- if type( line ) is GenomicInterval:
- out_file.write( "%s\n" % "\t".join( line.fields ) )
- elif type( line ) is list:
- out_file.write( "%s\n" % "\t".join( line ) )
- else:
- out_file.write( "%s\n" % line )
- except ParseError, exc:
- out_file.close()
- fail( "Invalid file format: %s" % str( exc ) )
-
- out_file.close()
-
- if g1.skipped > 0:
- print skipped( g1, filedesc=" of 1st dataset" )
-
-if __name__ == "__main__":
- main()
diff --git a/tools/new_operations/gops_subtract.py b/tools/new_operations/gops_subtract.py
deleted file mode 100644
index 9a8aa0d66ad..00000000000
--- a/tools/new_operations/gops_subtract.py
+++ /dev/null
@@ -1,99 +0,0 @@
-#!/usr/bin/env python
-"""
-Find regions of first interval file that do not overlap regions in a second
-interval file. Interval files can either be BED or GFF format.
-
-usage: %prog interval_file_1 interval_file_2 out_file
- -1, --cols1=N,N,N,N: Columns for start, end, strand in first file
- -2, --cols2=N,N,N,N: Columns for start, end, strand in second file
- -m, --mincols=N: Require this much overlap (default 1bp)
- -p, --pieces: just print pieces of second set (after padding)
- -G, --gff1: input 1 is GFF format, meaning start and end coordinates are 1-based, closed interval
- -H, --gff2: input 2 is GFF format, meaning start and end coordinates are 1-based, closed interval
-"""
-from galaxy import eggs
-import pkg_resources
-pkg_resources.require( "bx-python" )
-import sys, traceback, fileinput
-from warnings import warn
-from bx.intervals import *
-from bx.intervals.io import *
-from bx.intervals.operations.subtract import *
-from bx.cookbook import doc_optparse
-from galaxy.tools.util.galaxyops import *
-from galaxy.datatypes.util.gff_util import GFFFeature, GFFReaderWrapper, convert_bed_coords_to_gff
-
-assert sys.version_info[:2] >= ( 2, 4 )
-
-def main():
- mincols = 1
- upstream_pad = 0
- downstream_pad = 0
-
- 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 )
- if options.mincols: mincols = int( options.mincols )
- pieces = bool( options.pieces )
- in1_gff_format = bool( options.gff1 )
- in2_gff_format = bool( options.gff2 )
- in_fname, in2_fname, out_fname = args
- except:
- doc_optparse.exception()
-
- # Set readers to handle either GFF or default format.
- if in1_gff_format:
- in1_reader_wrapper = GFFReaderWrapper
- else:
- in1_reader_wrapper = NiceReaderWrapper
- if in2_gff_format:
- in2_reader_wrapper = GFFReaderWrapper
- else:
- in2_reader_wrapper = NiceReaderWrapper
-
- g1 = in1_reader_wrapper( fileinput.FileInput( in_fname ),
- chrom_col=chr_col_1,
- start_col=start_col_1,
- end_col=end_col_1,
- strand_col=strand_col_1,
- fix_strand=True )
- if in1_gff_format:
- # Subtract requires coordinates in BED format.
- g1.convert_to_bed_coord=True
-
- g2 = in2_reader_wrapper( fileinput.FileInput( in2_fname ),
- chrom_col=chr_col_2,
- start_col=start_col_2,
- end_col=end_col_2,
- strand_col=strand_col_2,
- fix_strand=True )
- if in2_gff_format:
- # Subtract requires coordinates in BED format.
- g2.convert_to_bed_coord=True
-
- out_file = open( out_fname, "w" )
- try:
- for feature in subtract( [g1,g2], pieces=pieces, mincols=mincols ):
- if isinstance( feature, GFFFeature ):
- # Convert back to GFF coordinates since reader converted automatically.
- convert_bed_coords_to_gff( feature )
- for interval in feature.intervals:
- out_file.write( "%s\n" % "\t".join( interval.fields ) )
- elif isinstance( feature, GenomicInterval ):
- out_file.write( "%s\n" % "\t".join( feature.fields ) )
- else:
- out_file.write( "%s\n" % feature )
- except ParseError, exc:
- out_file.close()
- fail( "Invalid file format: %s" % str( exc ) )
-
- out_file.close()
-
- if g1.skipped > 0:
- print skipped( g1, filedesc=" of 2nd dataset" )
- if g2.skipped > 0:
- print skipped( g2, filedesc=" of 1st dataset" )
-
-if __name__ == "__main__":
- main()
diff --git a/tools/new_operations/intersect.xml b/tools/new_operations/intersect.xml
deleted file mode 100644
index 642a94d34cd..00000000000
--- a/tools/new_operations/intersect.xml
+++ /dev/null
@@ -1,143 +0,0 @@
-
- the intervals of two datasets
- gops_intersect.py
- $input1 $input2 $output
-
- #if isinstance( $input1.datatype, $__app__.datatypes_registry.get_datatype_by_extension('gff').__class__):
- -1 1,4,5,7 --gff1
- #else:
- -1 ${input1.metadata.chromCol},${input1.metadata.startCol},${input1.metadata.endCol},${input1.metadata.strandCol}
- #end if
-
- #if isinstance( $input2.datatype, $__app__.datatypes_registry.get_datatype_by_extension('gff').__class__):
- -2 1,4,5,7 --gff2
- #else:
- -2 ${input2.metadata.chromCol},${input2.metadata.startCol},${input2.metadata.endCol},${input2.metadata.strandCol}
- #end if
-
- -m $min $returntype
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-.. class:: infomark
-
-**TIP:** If your dataset does not appear in the pulldown menu, it means that it is not in interval format. Use "edit attributes" to set chromosome, start, end, and strand columns.
-
------
-
-**Screencasts!**
-
-See Galaxy Interval Operation Screencasts_ (right click to open this link in another window).
-
-.. _Screencasts: http://wiki.g2.bx.psu.edu/Learn/Interval%20Operations
-
------
-
-**Syntax**
-
-- **Where overlap is at least** sets the minimum length (in base pairs) of overlap between elements of the two datasets
-- **Overlapping Intervals** returns entire intervals from the first dataset that overlap the second dataset. The returned intervals are completely unchanged, and this option only filters out intervals that do not overlap with the second dataset.
-- **Overlapping pieces of Intervals** returns intervals that indicate the exact base pair overlap between the first dataset and the second dataset. The intervals returned are from the first dataset, and all fields besides start and end are guaranteed to remain unchanged.
-
------
-
-**Examples**
-
-Overlapping Intervals:
-
-.. image:: ${static_path}/operation_icons/gops_intersectOverlappingIntervals.gif
-
-Overlapping Pieces of Intervals:
-
-.. image:: ${static_path}/operation_icons/gops_intersectOverlappingPieces.gif
-
-
-
diff --git a/tools/new_operations/join.xml b/tools/new_operations/join.xml
deleted file mode 100644
index 4d5331ec995..00000000000
--- a/tools/new_operations/join.xml
+++ /dev/null
@@ -1,117 +0,0 @@
-
- the intervals of two datasets side-by-side
- gops_join.py $input1 $input2 $output -1 ${input1.metadata.chromCol},${input1.metadata.startCol},${input1.metadata.endCol},${input1.metadata.strandCol} -2 ${input2.metadata.chromCol},${input2.metadata.startCol},${input2.metadata.endCol},${input2.metadata.strandCol} -m $min -f $fill
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-.. class:: infomark
-
-**TIP:** If your dataset does not appear in the pulldown menu, it means that it is not in interval format. Use "edit attributes" to set chromosome, start, end, and strand columns.
-
------
-
-**Screencasts!**
-
-See Galaxy Interval Operation Screencasts_ (right click to open this link in another window).
-
-.. _Screencasts: http://wiki.g2.bx.psu.edu/Learn/Interval%20Operations
-
------
-
-**Syntax**
-
-- **Where overlap** specifies the minimum overlap between intervals that allows them to be joined.
-- **Return only records that are joined** returns only the records of the first dataset that join to a record in the second dataset. This is analogous to an INNER JOIN.
-- **Return all records of first dataset (fill null with ".")** returns all intervals of the first dataset, and any intervals that do not join an interval from the second dataset are filled in with a period(.). This is analogous to a LEFT JOIN.
-- **Return all records of second dataset (fill null with ".")** returns all intervals of the second dataset, and any intervals that do not join an interval from the first dataset are filled in with a period(.). **Note that this may produce an invalid interval file, since a period(.) is not a valid chrom, start, end or strand.**
-- **Return all records of both datasets (fill nulls with ".")** returns all records from both datasets, and fills on either the right or left with periods. **Note that this may produce an invalid interval file, since a period(.) is not a valid chrom, start, end or strand.**
-
------
-
-**Examples**
-
-.. image:: ${static_path}/operation_icons/gops_joinRecordsList.gif
-
-Only records that are joined (inner join):
-
-.. image:: ${static_path}/operation_icons/gops_joinInner.gif
-
-All records of first dataset:
-
-.. image:: ${static_path}/operation_icons/gops_joinLeftOuter.gif
-
-All records of second dataset:
-
-.. image:: ${static_path}/operation_icons/gops_joinRightOuter.gif
-
-All records of both datasets:
-
-.. image:: ${static_path}/operation_icons/gops_joinFullOuter.gif
-
-
-
-
diff --git a/tools/new_operations/merge.xml b/tools/new_operations/merge.xml
deleted file mode 100644
index 33c2abc096f..00000000000
--- a/tools/new_operations/merge.xml
+++ /dev/null
@@ -1,58 +0,0 @@
-
- the overlapping intervals of a dataset
- gops_merge.py $input1 $output -1 ${input1.metadata.chromCol},${input1.metadata.startCol},${input1.metadata.endCol},${input1.metadata.strandCol} $returntype
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-.. class:: infomark
-
-**TIP:** If your dataset does not appear in the pulldown menu, it means that it is not in interval format. Use "edit attributes" to set chromosome, start, end, and strand columns.
-
------
-
-**Screencasts!**
-
-See Galaxy Interval Operation Screencasts_ (right click to open this link in another window).
-
-.. _Screencasts: http://wiki.g2.bx.psu.edu/Learn/Interval%20Operations
-
------
-
-This operation merges all overlapping intervals into single intervals.
-
-**Example**
-
-.. image:: ${static_path}/operation_icons/gops_merge.gif
-
-
-
\ No newline at end of file
diff --git a/tools/new_operations/operation_filter.py b/tools/new_operations/operation_filter.py
deleted file mode 100644
index a44a043cbd6..00000000000
--- a/tools/new_operations/operation_filter.py
+++ /dev/null
@@ -1,99 +0,0 @@
-# runs after the job (and after the default post-filter)
-import os
-from galaxy import eggs
-from galaxy import jobs
-from galaxy.tools.parameters import DataToolParameter
-
-from galaxy.jobs.handler import JOB_ERROR
-
-# Older py compatibility
-try:
- set()
-except:
- from sets import Set as set
-
-#def exec_before_process(app, inp_data, out_data, param_dict, tool=None):
-# """Sets the name of the data"""
-# dbkeys = sets.Set( [data.dbkey for data in inp_data.values() ] )
-# if len(dbkeys) != 1:
-# raise Exception, 'Both Queries must be from the same genome build
'
-
-def validate_input( trans, error_map, param_values, page_param_map ):
- dbkeys = set()
- data_param_names = set()
- data_params = 0
- for name, param in page_param_map.iteritems():
- if isinstance( param, DataToolParameter ):
- # for each dataset parameter
- if param_values.get(name, None) != None:
- dbkeys.add( param_values[name].dbkey )
- data_params += 1
- # check meta data
- try:
- param = param_values[name]
- if isinstance( param.datatype, trans.app.datatypes_registry.get_datatype_by_extension( 'gff' ).__class__ ):
- # TODO: currently cannot validate GFF inputs b/c they are not derived from interval.
- pass
- else: # Validate interval datatype.
- startCol = int( param.metadata.startCol )
- endCol = int( param.metadata.endCol )
- chromCol = int( param.metadata.chromCol )
- if param.metadata.strandCol is not None:
- strandCol = int ( param.metadata.strandCol )
- else:
- strandCol = 0
- except:
- error_msg = "The attributes of this dataset are not properly set. " + \
- "Click the pencil icon in the history item to set the chrom, start, end and strand columns."
- error_map[name] = error_msg
- data_param_names.add( name )
- if len( dbkeys ) > 1:
- for name in data_param_names:
- error_map[name] = "All datasets must belong to same genomic build, " \
- "this dataset is linked to build '%s'" % param_values[name].dbkey
- if data_params != len(data_param_names):
- for name in data_param_names:
- error_map[name] = "A dataset of the appropriate type is required"
-
-# Commented out by INS, 5/30/2007. What is the PURPOSE of this?
-def exec_after_process(app, inp_data, out_data, param_dict, tool=None, stdout=None, stderr=None):
- """Verify the output data after each run"""
- items = out_data.items()
-
- for name, data in items:
- try:
- if stderr and len( stderr ) > 0:
- raise Exception( stderr )
-
- except Exception, exc:
- data.blurb = JOB_ERROR
- data.state = JOB_ERROR
-
-## def exec_after_process(app, inp_data, out_data, param_dict, tool=None, stdout=None, stderr=None):
-## pass
-
-
-def exec_after_merge(app, inp_data, out_data, param_dict, tool=None, stdout=None, stderr=None):
- exec_after_process(
- app, inp_data, out_data, param_dict, tool=tool, stdout=stdout, stderr=stderr)
-
- # strip strand column if clusters were merged
- items = out_data.items()
- for name, data in items:
- if param_dict['returntype'] == True:
- data.metadata.chromCol = 1
- data.metadata.startCol = 2
- data.metadata.endCol = 3
- # merge always clobbers strand
- data.metadata.strandCol = None
-
-
-def exec_after_cluster(app, inp_data, out_data, param_dict, tool=None, stdout=None, stderr=None):
- exec_after_process(
- app, inp_data, out_data, param_dict, tool=tool, stdout=stdout, stderr=stderr)
-
- # strip strand column if clusters were merged
- if param_dict["returntype"] == '1':
- items = out_data.items()
- for name, data in items:
- data.metadata.strandCol = None
diff --git a/tools/new_operations/subtract.xml b/tools/new_operations/subtract.xml
deleted file mode 100644
index 9ab69748252..00000000000
--- a/tools/new_operations/subtract.xml
+++ /dev/null
@@ -1,124 +0,0 @@
-
- the intervals of two datasets
- gops_subtract.py
- $input1 $input2 $output
-
- #if isinstance( $input1.datatype, $__app__.datatypes_registry.get_datatype_by_extension('gff').__class__):
- -1 1,4,5,7 --gff1
- #else:
- -1 ${input1.metadata.chromCol},${input1.metadata.startCol},${input1.metadata.endCol},${input1.metadata.strandCol}
- #end if
-
- #if isinstance( $input2.datatype, $__app__.datatypes_registry.get_datatype_by_extension('gff').__class__):
- -2 1,4,5,7 --gff2
- #else:
- -2 ${input2.metadata.chromCol},${input2.metadata.startCol},${input2.metadata.endCol},${input2.metadata.strandCol}
- #end if
-
- -m $min $returntype
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-.. class:: infomark
-
-**TIP:** If your dataset does not appear in the pulldown menu, it means that it is not in interval format. Use "edit attributes" to set chromosome, start, end, and strand columns.
-
------
-
-**Screencasts!**
-
-See Galaxy Interval Operation Screencasts_ (right click to open this link in another window).
-
-.. _Screencasts: http://wiki.g2.bx.psu.edu/Learn/Interval%20Operations
-
------
-
-**Syntax**
-
-- **Where overlap is at least** sets the minimum length (in base pairs) of overlap between elements of the two datasets.
-- **Intervals with no overlap** returns entire intervals from the first dataset that do not overlap the second dataset. The returned intervals are completely unchanged, and this option only filters out intervals that overlap with the second dataset.
-- **Non-overlapping pieces of intervals** returns intervals from the first dataset that have the intervals from the second dataset removed. Any overlapping base pairs are removed from the range of the interval. All fields besides start and end are guaranteed to remain unchanged.
-
------
-
-**Example**
-
-Intervals with no overlap:
-
-.. image:: ${static_path}/operation_icons/gops_subtractOverlappingIntervals.gif
-
-Non-overlapping pieces of intervals:
-
-.. image:: ${static_path}/operation_icons/gops_subtractOverlappingPieces.gif
-
-
-
diff --git a/tools/new_operations/subtract_query.py b/tools/new_operations/subtract_query.py
deleted file mode 100644
index b06440dd0ea..00000000000
--- a/tools/new_operations/subtract_query.py
+++ /dev/null
@@ -1,113 +0,0 @@
-#!/usr/bin/env python
-# Greg Von Kuster
-
-"""
-Subtract an entire query from another query
-usage: %prog in_file_1 in_file_2 begin_col end_col output
- --ignore-empty-end-cols: ignore empty end columns when subtracting
-"""
-import sys, re
-from galaxy import eggs
-import pkg_resources; pkg_resources.require( "bx-python" )
-from bx.cookbook import doc_optparse
-
-# Older py compatibility
-try:
- set()
-except:
- from sets import Set as set
-
-assert sys.version_info[:2] >= ( 2, 4 )
-
-def get_lines(fname, begin_col='', end_col='', ignore_empty_end_cols=False):
- lines = set([])
- i = 0
- for i, line in enumerate(file(fname)):
- line = line.rstrip('\r\n')
- if line and not line.startswith('#'):
- if begin_col and end_col:
- """Both begin_col and end_col must be integers at this point."""
- try:
- line = line.split('\t')
- line = '\t'.join([line[j] for j in range(begin_col-1, end_col)])
- if ignore_empty_end_cols:
- # removing empty fields, we do not compare empty fields at the end of a line.
- line = line.rstrip()
- lines.add( line )
- except: pass
- else:
- if ignore_empty_end_cols:
- # removing empty fields, we do not compare empty fields at the end of a line.
- line = line.rstrip()
- lines.add( line )
- if i: return (i+1, lines)
- else: return (i, lines)
-
-def main():
-
- # Parsing Command Line here
- options, args = doc_optparse.parse( __doc__ )
-
- try:
- inp1_file, inp2_file, begin_col, end_col, out_file = args
- except:
- doc_optparse.exception()
-
- begin_col = begin_col.strip()
- end_col = end_col.strip()
-
- if begin_col != 'None' or end_col != 'None':
- """
- The user selected columns for restriction. We'll allow default
- values for both begin_col and end_col as long as the user selected
- at least one of them for restriction.
- """
- if begin_col == 'None':
- begin_col = end_col
- elif end_col == 'None':
- end_col = begin_col
- begin_col = int(begin_col)
- end_col = int(end_col)
- """Make sure that begin_col <= end_col (switch if not)"""
- if begin_col > end_col:
- tmp_col = end_col
- end_col = begin_col
- begin_col = tmp_col
- else:
- begin_col = end_col = ''
-
- try:
- fo = open(out_file,'w')
- except:
- print >> sys.stderr, "Unable to open output file"
- sys.exit()
-
- """
- len1 is the number of lines in inp1_file
- lines1 is the set of unique lines in inp1_file
- diff1 is the number of duplicate lines removed from inp1_file
- """
- len1, lines1 = get_lines(inp1_file, begin_col, end_col, options.ignore_empty_end_cols)
- diff1 = len1 - len(lines1)
- len2, lines2 = get_lines(inp2_file, begin_col, end_col, options.ignore_empty_end_cols)
-
- lines1.difference_update(lines2)
- """lines1 is now the set of unique lines in inp1_file - the set of unique lines in inp2_file"""
-
- for line in lines1:
- print >> fo, line
-
- fo.close()
-
- info_msg = 'Subtracted %d lines. ' %((len1 - diff1) - len(lines1))
-
- if begin_col and end_col:
- info_msg += 'Restricted to columns c' + str(begin_col) + ' thru c' + str(end_col) + '. '
-
- if diff1 > 0:
- info_msg += 'Eliminated %d duplicate/blank/comment/invalid lines from first query.' %diff1
-
- print info_msg
-
-if __name__ == "__main__":
- main()
diff --git a/tools/new_operations/subtract_query.xml b/tools/new_operations/subtract_query.xml
deleted file mode 100644
index a603ce144be..00000000000
--- a/tools/new_operations/subtract_query.xml
+++ /dev/null
@@ -1,126 +0,0 @@
-
- from another dataset
-
- subtract_query.py $input1 $input2 $begin_col $end_col $output
- #if str($ignore_empty_end_cols) == 'true':
- --ignore-empty-end-cols
- #end if
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-.. class:: infomark
-
-**TIP:** This tool complements the tool in the **Operate on Genomic Intervals** tool set which subtracts the intervals of two datasets.
-
-
------
-
-**Syntax**
-
-This tool subtracts an entire dataset from another dataset.
-
-- Any text format is valid.
-- If both dataset formats are tabular, you may restrict the subtraction to specific columns **contained in both datasets** and the resulting dataset will include only the columns specified.
-- The begin column must be less than or equal to the end column. If it is not, begin column is switched with end column.
-- If begin column is specified but end column is not, end column will default to begin_column (and vice versa).
-- All blank and comment lines are skipped and not included in the resulting dataset (comment lines are lines beginning with a # character).
-- Duplicate lines are eliminated from both dataset prior to subtraction. If any duplicate lines were eliminated from the first dataset, the number is displayed in the resulting history item.
-
------
-
-**Example**
-
-If this is the **First dataset**::
-
- chr1 4225 19670
- chr10 6 8
- chr1 24417 24420
- chr6_hla_hap2 0 150
- chr2 1 5
- chr10 2 10
- chr1 30 55
- chrY 1 20
- chr1 1225979 42287290
- chr10 7 8
-
-and this is the **Second dataset**::
-
- chr1 4225 19670
- chr10 6 8
- chr1 24417 24420
- chr6_hla_hap2 0 150
- chr2 1 5
- chr1 30 55
- chrY 1 20
- chr1 1225979 42287290
-
-Subtracting the **Second dataset** from the **First dataset** (including all columns) will yield::
-
- chr10 7 8
- chr10 2 10
-
-Conversely, subtracting the **First dataset** from the **Second dataset** (including all columns) will result in an empty dataset.
-
-Subtracting the **Second dataset** from the **First dataset** (restricting to columns c1 and c2) will yield::
-
- chr10 7
- chr10 2
-
-
-
\ No newline at end of file
diff --git a/tools/new_operations/tables_arithmetic_operations.pl b/tools/new_operations/tables_arithmetic_operations.pl
deleted file mode 100644
index e5b0ce2e3f8..00000000000
--- a/tools/new_operations/tables_arithmetic_operations.pl
+++ /dev/null
@@ -1,117 +0,0 @@
-# A program to implement arithmetic operations on tabular files data. The program takes three inputs:
-# The first input is a TABULAR format file containing numbers only.
-# The second input is a TABULAR format file containing numbers only.
-# The two files must have the same number of columns and the same number of rows
-# The third input is an arithmetic operation: +, -, *, or / for addition, subtraction, multiplication, or division, respectively
-# The output file is a TABULAR format file containing the result of implementing the arithmetic operation on both input files.
-# The output file has the same number of columns and the same number of rows as each of the two input files.
-# Note: in case of division, none of the values in the second input file could be 0.
-
-use strict;
-use warnings;
-
-#variables to handle information of the first input tabular file
-my $lineData1 = "";
-my @lineDataArray1 = ();
-my $lineArraySize = 0;
-my $lineCounter1 = 0;
-
-#variables to handle information of the second input tabular file
-my $lineData2= "";
-my @lineDataArray2 = ();
-my $lineCounter2 = 0;
-
-my $result = 0;
-
-# check to make sure having the correct number of arguments
-my $usage = "usage: tables_arithmetic_operations.pl [TABULAR.in] [TABULAR.in] [ArithmeticOperation] [TABULAR.out] \n";
-die $usage unless @ARGV == 4;
-
-#variables to store the names of input and output files
-my $inputTabularFile1 = $ARGV[0];
-my $inputTabularFile2 = $ARGV[1];
-my $arithmeticOperation = $ARGV[2];
-my $outputTabularFile = $ARGV[3];
-
-#open the input and output files
-open (INPUT1, "<", $inputTabularFile1) || die("Could not open file $inputTabularFile1 \n");
-open (INPUT2, "<", $inputTabularFile2) || die("Could not open file $inputTabularFile2 \n");
-open (OUTPUT, ">", $outputTabularFile) || die("Could not open file $outputTabularFile \n");
-
-#store the first input file in the array @motifsFrequencyData1
-my @tabularData1 = ;
-
-#store the second input file in the array @motifsFrequencyData2
-my @tabularData2 = ;
-
-#reset the $lineCounter1 to 0
-$lineCounter1 = 0;
-
-#iterated through the lines of the first input file
-INDEL1:
-foreach $lineData1 (@tabularData1){
- chomp ($lineData1);
- $lineCounter1++;
-
- #reset the $lineCounter2 to 0
- $lineCounter2 = 0;
-
- #iterated through the lines of the second input file
- foreach $lineData2 (@tabularData2){
- chomp ($lineData2);
- $lineCounter2++;
-
- #check if the two motifs are the same in the two input files
- if ($lineCounter1 == $lineCounter2){
-
- @lineDataArray1 = split(/\t/, $lineData1);
- @lineDataArray2 = split(/\t/, $lineData2);
-
- $lineArraySize = @lineDataArray1;
-
- for (my $index = 0; $index < $lineArraySize; $index++){
-
- if ($arithmeticOperation eq "Addition"){
- #compute the additin of both values
- $result = $lineDataArray1[$index] + $lineDataArray2[$index];
- }
-
- if ($arithmeticOperation eq "Subtraction"){
- #compute the subtraction of both values
- $result = $lineDataArray1[$index] - $lineDataArray2[$index];
- }
-
- if ($arithmeticOperation eq "Multiplication"){
- #compute the multiplication of both values
- $result = $lineDataArray1[$index] * $lineDataArray2[$index];
- }
-
- if ($arithmeticOperation eq "Division"){
-
- #check if the denominator is 0
- if ($lineDataArray2[$index] != 0){
- #compute the division of both values
- $result = $lineDataArray1[$index] / $lineDataArray2[$index];
- }
- else{
- die("A denominator could not be zero \n");
- }
- }
-
- #store the result in the output file
- if ($index < $lineArraySize - 1){
- print OUTPUT $result . "\t";
- }
- else{
- print OUTPUT $result . "\n";
- }
- }
- next INDEL1;
- }
- }
-}
-
-#close the input and output files
-close(OUTPUT);
-close(INPUT2);
-close(INPUT1);
\ No newline at end of file
diff --git a/tools/new_operations/tables_arithmetic_operations.xml b/tools/new_operations/tables_arithmetic_operations.xml
deleted file mode 100644
index 0e2f891de95..00000000000
--- a/tools/new_operations/tables_arithmetic_operations.xml
+++ /dev/null
@@ -1,105 +0,0 @@
-
- on tables
-
-
- tables_arithmetic_operations.pl $inputFile1 $inputFile2 $inputArithmeticOperation3 $outputFile1
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-.. class:: infomark
-
-**What it does**
-
-This program implements arithmetic operations on tabular files data. The program takes three inputs:
-
-- The first input is a TABULAR format file containing numbers only.
-- The second input is a TABULAR format file containing numbers only.
-- The third input is an arithmetic operation: +, -, x, or / for addition, subtraction, multiplication, or division, respectively.
-- The output file is a TABULAR format file containing the result of implementing the arithmetic operation on both input files.
-
-
-Notes:
-
-- The two files must have the same number of columns and the same number of rows.
-- The output file has the same number of columns and the same number of rows as each of the two input files.
-- In case of division, none of the values in the second input file could be 0, otherwise the program will stop and report an error.
-
-**Example**
-
-Let us have the first input file as follows::
-
- 5 4 0
- 10 11 12
- 1 3 1
- 1 2 1
- 2 0 4
-
-And the second input file as follows::
-
- 5 4 4
- 2 5 8
- 1 2 1
- 3 2 5
- 2 4 4
-
-Running the program and choosing "Addition" as an arithmetic operation will give the following output::
-
- 10 8 4
- 12 16 20
- 2 5 2
- 4 4 6
- 4 4 8
-
-
-
-
-
diff --git a/tools/regVariation/WeightedAverage.py b/tools/regVariation/WeightedAverage.py
deleted file mode 100755
index 8c4c934ebc5..00000000000
--- a/tools/regVariation/WeightedAverage.py
+++ /dev/null
@@ -1,94 +0,0 @@
-#!/usr/bin/env python
-"""
-usage: %prog bed_file_1 bed_file_2 out_file
- -1, --cols1=N,N,N,N: Columns for chr, start, end, strand in first file
- -2, --cols2=N,N,N,N,N: Columns for chr, start, end, strand, name/value in second file
-"""
-
-import collections
-import sys
-#import numpy
-from galaxy import eggs
-import pkg_resources
-pkg_resources.require( "bx-python" )
-from galaxy.tools.util.galaxyops import *
-from bx.cookbook import doc_optparse
-
-
-#export PYTHONPATH=~/galaxy/lib/
-#running command python WeightedAverage.py interval_interpolate.bed value_interpolate.bed interpolate_result.bed
-
-def stop_err(msg):
- sys.stderr.write(msg)
- sys.exit()
-
-
-def FindRate(chromosome, start_stop, dictType):
- OverlapList = []
- for tempO in dictType[chromosome]:
- DatabaseInterval = [tempO[0], tempO[1]]
- Overlap = GetOverlap( start_stop, DatabaseInterval )
- if Overlap > 0:
- OverlapList.append([Overlap, tempO[2]])
-
- if len(OverlapList) > 0:
- SumRecomb = 0
- SumOverlap = 0
- for member in OverlapList:
- SumRecomb += member[0]*member[1]
- SumOverlap += member[0]
- averageRate = SumRecomb/SumOverlap
- return averageRate
- else:
- return 'NA'
-
-
-def GetOverlap(a, b):
- return min(a[1], b[1])-max(a[0], b[0])
-
-
-options, args = doc_optparse.parse( __doc__ )
-
-try:
- chr_col_1, start_col_1, end_col_1, strand_col1 = parse_cols_arg( options.cols1 )
- chr_col_2, start_col_2, end_col_2, strand_col2, name_col_2 = parse_cols_arg( options.cols2 )
- input1, input2, input3 = args
-except Exception, eee:
- print eee
- stop_err( "Data issue: click the pencil icon in the history item to correct the metadata attributes." )
-
-fd2 = open(input2)
-lines2 = fd2.readlines()
-RecombChrDict = collections.defaultdict(list)
-
-skipped = 0
-for line in lines2:
- temp = line.strip().split()
- try:
- assert float(temp[int(name_col_2)])
- except:
- skipped += 1
- continue
- tempIndex = [int(temp[int(start_col_2)]), int(temp[int(end_col_2)]), float(temp[int(name_col_2)])]
- RecombChrDict[temp[int(chr_col_2)]].append(tempIndex)
-
-print "Skipped %d features with invalid values" % (skipped)
-
-fd1 = open(input1)
-lines = fd1.readlines()
-finalProduct = ''
-for line in lines:
- temp = line.strip().split('\t')
- chromosome = temp[int(chr_col_1)]
- start = int(temp[int(start_col_1)])
- stop = int(temp[int(end_col_1)])
- start_stop = [start, stop]
- RecombRate = FindRate( chromosome, start_stop, RecombChrDict )
- try:
- RecombRate = "%.4f" % (float(RecombRate))
- except:
- RecombRate = RecombRate
- finalProduct += line.strip()+'\t'+str(RecombRate)+'\n'
-fdd = open(input3, 'w')
-fdd.writelines(finalProduct)
-fdd.close()
diff --git a/tools/regVariation/WeightedAverage.xml b/tools/regVariation/WeightedAverage.xml
deleted file mode 100755
index 9e9406985d2..00000000000
--- a/tools/regVariation/WeightedAverage.xml
+++ /dev/null
@@ -1,71 +0,0 @@
-
- of the values of features overlapping an interval
- WeightedAverage.py $genomic_interval $genomic_feature $out_file1 -1 ${genomic_interval.metadata.chromCol},${genomic_interval.metadata.startCol},${genomic_interval.metadata.endCol},${genomic_interval.metadata.strandCol} -2 ${genomic_feature.metadata.chromCol},${genomic_feature.metadata.startCol},${genomic_feature.metadata.endCol},${genomic_feature.metadata.strandCol},${genomic_feature.metadata.nameCol}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-.. class:: infomark
-
-**What it does**
-
-For each interval in your first dataset, this tool calculates the weighted average value of the overlapping features in your second dataset.
-
-- When a genomic interval partially or totally overlaps a single genomic feature, the value of that genomic feature is assigned to the genomic interval.
-- When a genomic interval partially or totally overlaps with more than one genomic features, the average of the values of the overlapping genomic features weighted by the corresponding number of overlapping bases is assigned to the genomic interval.
-- When a genomic interval does not overlap with any genomic feature, 'NA' will be assigned as it's value.
-
------
-
-.. class:: warningmark
-
-**Note**
-
-The input datasets should be in **bed** or **interval** format. Please use "edit attributes"/pencil icon to specify the column containing the values for the features in the second dataset as **name/identifier** column.
-
-The output will contain all the columns in the first input plus a new column containing the assigned value for each interval.
-
------
-
-**Example**
-
-- Suppose our first dataset contains the following **genomic intervals**::
-
- chr start stop
- chr1 1000 2000
- chr1 3000 5000
- chr1 8000 9000
-
-- and our second dataset contains the following **genomic features** each having an associated value (in fourth column) ::
-
- chr start stop name
- chr1 900 1200 0.5
- chr1 2900 3100 0.2
- chr1 4800 5100 0.8
-
-- For each **genomic interval** in our first dataset, this tool calculates the weighted average value of the overlapping **genomic features** in our second dataset ::
-
- chr1 1000 2000 0.5
- chr1 3000 5000 0.6
- chr1 8000 9000 NA
-
-
-
-
-
\ No newline at end of file
diff --git a/tools/regVariation/best_regression_subsets.py b/tools/regVariation/best_regression_subsets.py
deleted file mode 100644
index a4883e33c3d..00000000000
--- a/tools/regVariation/best_regression_subsets.py
+++ /dev/null
@@ -1,91 +0,0 @@
-#!/usr/bin/env python
-
-from galaxy import eggs
-
-import sys
-from rpy import *
-import numpy
-
-def stop_err(msg):
- sys.stderr.write(msg)
- sys.exit()
-
-
-infile = sys.argv[1]
-y_col = int(sys.argv[2])-1
-x_cols = sys.argv[3].split(',')
-outfile = sys.argv[4]
-outfile2 = sys.argv[5]
-print "Predictor columns: %s; Response column: %d" % ( x_cols, y_col+1 )
-fout = open(outfile,'w')
-
-for i, line in enumerate( file ( infile )):
- line = line.rstrip('\r\n')
- if len( line )>0 and not line.startswith( '#' ):
- elems = line.split( '\t' )
- break
- if i == 30:
- break # Hopefully we'll never get here...
-
-if len( elems )<1:
- stop_err( "The data in your input dataset is either missing or not formatted properly." )
-
-y_vals = []
-x_vals = []
-
-for k, col in enumerate(x_cols):
- x_cols[k] = int(col)-1
- x_vals.append([])
-
-NA = 'NA'
-for ind, line in enumerate( file( infile ) ):
- if line and not line.startswith( '#' ):
- try:
- fields = line.split("\t")
- try:
- yval = float(fields[y_col])
- except Exception, ey:
- yval = r('NA')
- y_vals.append(yval)
- for k, col in enumerate(x_cols):
- try:
- xval = float(fields[col])
- except Exception, ex:
- xval = r('NA')
- x_vals[k].append(xval)
- except:
- pass
-
-response_term = ""
-
-x_vals1 = numpy.asarray(x_vals).transpose()
-
-dat = r.list(x=array(x_vals1), y=y_vals)
-
-r.library("leaps")
-
-set_default_mode(NO_CONVERSION)
-try:
- leaps = r.regsubsets(r("y ~ x"), data= r.na_exclude(dat))
-except RException, rex:
- stop_err("Error performing linear regression on the input data.\nEither the response column or one of the predictor columns contain no numeric values.")
-set_default_mode(BASIC_CONVERSION)
-
-summary = r.summary(leaps)
-tot = len(x_vals)
-pattern = "["
-for i in range(tot):
- pattern = pattern + 'c' + str(int(x_cols[int(i)]) + 1) + ' '
-pattern = pattern.strip() + ']'
-print >> fout, "#Vars\t%s\tR-sq\tAdj. R-sq\tC-p\tbic" % (pattern)
-for ind, item in enumerate(summary['outmat']):
- print >> fout, "%s\t%s\t%s\t%s\t%s\t%s" % (str(item).count('*'), item, summary['rsq'][ind], summary['adjr2'][ind], summary['cp'][ind], summary['bic'][ind])
-
-
-r.pdf( outfile2, 8, 8 )
-r.plot(leaps, scale="Cp", main="Best subsets using Cp Criterion")
-r.plot(leaps, scale="r2", main="Best subsets using R-sq Criterion")
-r.plot(leaps, scale="adjr2", main="Best subsets using Adjusted R-sq Criterion")
-r.plot(leaps, scale="bic", main="Best subsets using bic Criterion")
-
-r.dev_off()
diff --git a/tools/regVariation/best_regression_subsets.xml b/tools/regVariation/best_regression_subsets.xml
deleted file mode 100644
index b93ae03ac15..00000000000
--- a/tools/regVariation/best_regression_subsets.xml
+++ /dev/null
@@ -1,66 +0,0 @@
-
-
-
- best_regression_subsets.py
- $input1
- $response_col
- $predictor_cols
- $out_file1
- $out_file2
- 1>/dev/null
- 2>/dev/null
-
-
-
-
-
-
-
-
-
-
-
-
-
- rpy
-
-
-
-
-
-
-.. class:: infomark
-
-**TIP:** If your data is not TAB delimited, use *Edit Datasets->Convert characters*
-
------
-
-.. class:: infomark
-
-**What it does**
-
-This tool uses the 'regsubsets' function from R statistical package for regression subset selection. It outputs two files, one containing a table with the best subsets and the corresponding summary statistics, and the other containing the graphical representation of the results.
-
------
-
-.. class:: warningmark
-
-**Note**
-
-- This tool currently treats all predictor and response variables as continuous variables.
-
-- Rows containing non-numeric (or missing) data in any of the chosen columns will be skipped from the analysis.
-
-- The 6 columns in the output are described below:
-
- - Column 1 (Vars): denotes the number of variables in the model
- - Column 2 ([c2 c3 c4...]): represents a list of the user-selected predictor variables (full model). An asterix denotes the presence of the corresponding predictor variable in the selected model.
- - Column 3 (R-sq): the fraction of variance explained by the model
- - Column 4 (Adj. R-sq): the above R-squared statistic adjusted, penalizing for higher number of predictors (p)
- - Column 5 (Cp): Mallow's Cp statistics
- - Column 6 (bic): Bayesian Information Criterion.
-
-
-
-
diff --git a/tools/regVariation/compute_q_values.pl b/tools/regVariation/compute_q_values.pl
deleted file mode 100644
index 6ff3e8f4881..00000000000
--- a/tools/regVariation/compute_q_values.pl
+++ /dev/null
@@ -1,95 +0,0 @@
-# A program to compute the q-values based on the p-values of multiple simultaneous tests.
-# The q-valules are computed using a specific R package created by John Storey called "qvalue".
-# The input is a TABULAR format file consisting of one column only that represents the p-values
-# of multiple simultaneous tests, one line for every p-value.
-# The first output is a TABULAR format file consisting of one column only that represents the q-values
-# corresponding to p-values, one line for every q-value.
-# the second output is a TABULAR format file consisting of three pages: the first page represents
-# the p-values histogram, the second page represents the q-values histogram, and the third page represents
-# the four Q-plots as introduced in the "qvalue" package manual.
-
-use strict;
-use warnings;
-use IO::Handle;
-use File::Temp qw/ tempfile tempdir /;
-my $tdir = tempdir( CLEANUP => 0 );
-
-# check to make sure having correct input and output files
-my $usage = "usage: compute_q_values.pl [TABULAR.in] [lambda] [pi0_method] [fdr_level] [robust] [TABULAR.out] [PDF.out] \n";
-die $usage unless @ARGV == 7;
-
-#get the input arguments
-my $p_valuesInputFile = $ARGV[0];
-my $lambdaValue = $ARGV[1];
-my $pi0_method = $ARGV[2];
-my $fdr_level = $ARGV[3];
-my $robustValue = $ARGV[4];
-my $q_valuesOutputFile = $ARGV[5];
-my $p_q_values_histograms_QPlotsFile = $ARGV[6];
-
-if($lambdaValue =~ /sequence/){
- $lambdaValue = "seq(0, 0.95, 0.05)";
-}
-
-#open the input files
-open (INPUT, "<", $p_valuesInputFile) || die("Could not open file $p_valuesInputFile \n");
-open (OUTPUT1, ">", $q_valuesOutputFile) || die("Could not open file $q_valuesOutputFile \n");
-open (OUTPUT2, ">", $p_q_values_histograms_QPlotsFile) || die("Could not open file $p_q_values_histograms_QPlotsFile \n");
-#open (ERROR, ">", "error.txt") or die ("Could not open file error.txt \n");
-
-#save all error messages into the error file $errorFile using the error file handle ERROR
-#STDERR -> fdopen( \*ERROR, "w" ) or die ("Could not direct errors to the error file error.txt \n");
-
-#warn "Hello Error File \n";
-
-#variable to store the name of the R script file
-my $r_script;
-
-# R script to implement the calcualtion of q-values based on multiple simultaneous tests p-values
-# construct an R script file and save it in a temp directory
-chdir $tdir;
-$r_script = "q_values_computation.r";
-
-open(Rcmd,">", $r_script) or die "Cannot open $r_script \n\n";
-print Rcmd "
- #options(show.error.messages = FALSE);
-
- #load necessary packages
- suppressPackageStartupMessages(library(tcltk));
- library(qvalue);
-
- #read the p-values of the multiple simultaneous tests from the input file $p_valuesInputFile
- p <- scan(\"$p_valuesInputFile\", quiet = TRUE);
-
- #compute the q-values that correspond to the p-values of the multiple simultaneous tests
- qobj <- qvalue(p, pi0.meth = \"$pi0_method\", lambda = $lambdaValue, fdr.level = $fdr_level, robust = $robustValue);
- #qobj <- qvalue(p, pi0.meth = \"smoother\", lambda = seq(0, 0.95, 0.05), fdr.level = 0.05);
- #qobj <- qvalue(p, pi0.meth = \"bootstrap\", fdr.level = 0.05);
-
- #draw the p-values histogram, the q-values histogram, and the four Q-plots
- # and save them on multiple pages of the output file $p_q_values_histograms_QPlotsFile
- pdf(file = \"$p_q_values_histograms_QPlotsFile\", width = 6.25, height = 6, family = \"Times\", pointsize = 12, onefile = TRUE)
- hist(qobj\$pvalues);
- #dev.off();
-
- hist(qobj\$qvalues);
- #dev.off();
-
- qplot(qobj);
- dev.off();
-
- #save the q-values in the output file $q_valuesOutputFile
- qobj\$pi0 <- signif(qobj\$pi0,digits=6)
- qwrite(qobj, filename=\"$q_valuesOutputFile\");
-
- #options(show.error.messages = TRUE);
- #eof\n";
-close Rcmd;
-
-system("R --no-restore --no-save --no-readline < $r_script > $r_script.out");
-
-#close the input and output and error files
-#close(ERROR);
-close(OUTPUT2);
-close(OUTPUT1);
-close(INPUT);
diff --git a/tools/regVariation/compute_q_values.xml b/tools/regVariation/compute_q_values.xml
deleted file mode 100644
index 075b1f090ad..00000000000
--- a/tools/regVariation/compute_q_values.xml
+++ /dev/null
@@ -1,155 +0,0 @@
-
- based on multiple simultaneous tests p-values
-
-
- compute_q_values.pl $inputFile1 $inputLambda2 $inputPI0_method3 $inputFDR_level4 $inputRobust5 $outputFile1 $outputFile2
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-.. class:: infomark
-
-**What it does**
-
-This program computes the q-values based on the p-values of multiple simultaneous tests. The q-values are computed using a specific R package, created by John Storey and Alan Dabney, called "qvalue". The program takes five inputs:
-
-- The first input is a TABULAR format file consisting of one column only that represents the p-values of multiple simultaneous tests, one line for every p-value.
-- The second input is the lambda parameter. The user can choose either the default: seq(0, 0.95, 0.05) or a decimal number between 0.0 and 1.0.
-- The third input is PI method which is either "smoother" or "bootstrap".
-- The fourth input is the FDR (false discovery rate) level which is a decimal number between 0.0 and 1.0.
-- The fifth input is either TRUE or FALSE for the estimate robustness.
-
-The program gives two outputs:
-
-- The first output is a TABULAR format file consisting of three columns:
-
- - the left column represents the p-values of multiple simultaneous tests, one line for every p-value
- - the middle column represents the q-values corresponding to the p-values
- - the third column represent the significance values, either 1 for significant or 0 for non-significant
-
-- The second output is a PDF format file consisting of three pages:
-
- - the first page represents the p-values histogram
- - the second page represents the q-values histogram
- - the third page represents the four Q-plots as introduced in the "qvalue" package manual.
-
-
-**Example**
-
-Let us have the first input file of p-values as follows::
-
- 0.140627492
- 0.432249886
- 0.122120877
- 0.142010182
- 0.012909858
- 0.000142807
- 0.039841941
- 0.035173303
- 0.011340057
- 1.01E-05
- 0.212738282
- 0.091256284
- 0.547375415
- 0.189589833
- 6.18E-12
- 0.001235875
- 1.10E-05
- 9.75E-07
- 2.13E-18
- 2.54E-16
- 1.20E-19
- 9.76E-14
- 0.359181534
- 0.03661672
- 0.400459987
- 0.387436466
- 0.342075061
- 0.904129283
- 0.031152635
-
-Running the program will give the following output::
-
- pi0: 0.140311054
-
- FDR level: 0.05
-
- p-value q-value significant
- 0.1406275 0.02889212 1
- 0.4322499 0.06514199 0
- 0.1221209 0.02760624 1
- 0.1420102 0.02889212 1
- 0.01290986 0.00437754 1
- 0.000142807 6.46E-05 1
- 0.03984194 0.01013235 1
- 0.0351733 0.009932946 1
- 0.01134006 0.004194811 1
- 1.01E-05 5.59E-06 1
- 0.2127383 0.03934711 1
- 0.09125628 0.02184257 1
- 0.5473754 0.07954578 0
- 0.1895898 0.03673547 1
- 6.18E-12 5.03E-12 1
- 0.001235875 0.00050288 1
- 1.10E-05 5.59E-06 1
- 9.75E-07 6.61E-07 1
- 2.13E-18 4.33E-18 1
- 2.54E-16 3.45E-16 1
- 1.20E-19 4.88E-19 1
- 9.76E-14 9.93E-14 1
- 0.3591815 0.06089654 0
- 0.03661672 0.009932946 1
- 0.40046 0.0626723 0
- 0.3874365 0.0626723 0
- 0.3420751 0.06051785 0
- 0.9041293 0.1268593 0
- 0.03115264 0.009750824 1
-
-
-.. image:: ${static_path}/operation_icons/p_hist.png
-
-
-.. image:: ${static_path}/operation_icons/q_hist.png
-
-
-.. image:: ${static_path}/operation_icons/Q_plots.png
-
-
-
-
-
diff --git a/tools/regVariation/featureCounter.py b/tools/regVariation/featureCounter.py
deleted file mode 100644
index 5cbfd428596..00000000000
--- a/tools/regVariation/featureCounter.py
+++ /dev/null
@@ -1,149 +0,0 @@
-#!/usr/bin/env python
-#Guruprasad Ananda
-"""
-Calculate count and coverage of one query on another, and append the Coverage and counts to
-the last four columns as bases covered, percent coverage, number of completely present features, number of partially present/overlapping features.
-
-usage: %prog bed_file_1 bed_file_2 out_file
- -1, --cols1=N,N,N,N: Columns for chr, start, end, strand in first file
- -2, --cols2=N,N,N,N: Columns for chr, start, end, strand in second file
-"""
-from galaxy import eggs
-import pkg_resources
-pkg_resources.require( "bx-python" )
-import sys, fileinput
-from bx.intervals.io import *
-from bx.cookbook import doc_optparse
-from bx.intervals.operations import quicksect
-from galaxy.tools.util.galaxyops import *
-
-assert sys.version_info[:2] >= ( 2, 4 )
-
-def stop_err(msg):
- sys.stderr.write(msg)
- sys.exit()
-
-def counter(node, start, end):
- global full, partial
- if node.start <= start and node.maxend > start:
- if node.end >= end or (node.start == start and end > node.end > start):
- full += 1
- elif end > node.end > start:
- partial += 1
- if node.left and node.left.maxend > start:
- counter(node.left, start, end)
- if node.right:
- counter(node.right, start, end)
- elif start < node.start < end:
- if node.end <= end:
- full += 1
- else:
- partial += 1
- if node.left and node.left.maxend > start:
- counter(node.left, start, end)
- if node.right:
- counter(node.right, start, end)
- else:
- if node.left:
- counter(node.left, start, end)
-
-def count_coverage( readers, comments=True ):
- primary = readers[0]
- secondary = readers[1]
- secondary_copy = readers[2]
-
- rightTree = quicksect.IntervalTree()
- for item in secondary:
- if type( item ) is GenomicInterval:
- rightTree.insert( item, secondary.linenum, item.fields )
-
- bitsets = secondary_copy.binned_bitsets()
-
- global full, partial
-
- for interval in primary:
- if type( interval ) is Header:
- yield interval
- if type( interval ) is Comment and comments:
- yield interval
- elif type( interval ) == GenomicInterval:
- chrom = interval.chrom
- start = int(interval.start)
- end = int(interval.end)
- full = 0
- partial = 0
- if chrom not in bitsets:
- bases_covered = 0
- percent = 0.0
- full = 0
- partial = 0
- else:
- bases_covered = bitsets[ chrom ].count_range( start, end-start )
- if (end - start) == 0:
- percent = 0
- else:
- percent = float(bases_covered) / float(end - start)
- if bases_covered:
- root = rightTree.chroms[chrom] #root node for the chrom tree
- counter(root, start, end)
- interval.fields.append(str(bases_covered))
- interval.fields.append(str(percent))
- interval.fields.append(str(full))
- interval.fields.append(str(partial))
- yield interval
-
-
-def main():
- 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 )
- in1_fname, in2_fname, out_fname = args
- except:
- stop_err( "Data issue: click the pencil icon in the history item to correct the metadata attributes." )
-
- g1 = NiceReaderWrapper( fileinput.FileInput( in1_fname ),
- chrom_col=chr_col_1,
- start_col=start_col_1,
- end_col=end_col_1,
- strand_col=strand_col_1,
- fix_strand=True )
- g2 = NiceReaderWrapper( fileinput.FileInput( in2_fname ),
- chrom_col=chr_col_2,
- start_col=start_col_2,
- end_col=end_col_2,
- strand_col=strand_col_2,
- fix_strand=True )
- g2_copy = NiceReaderWrapper( fileinput.FileInput( in2_fname ),
- chrom_col=chr_col_2,
- start_col=start_col_2,
- end_col=end_col_2,
- strand_col=strand_col_2,
- fix_strand=True )
-
-
- out_file = open( out_fname, "w" )
-
- try:
- for line in count_coverage([g1, g2, g2_copy]):
- if type( line ) is GenomicInterval:
- out_file.write( "%s\n" % "\t".join( line.fields ) )
- else:
- out_file.write( "%s\n" % line )
- except ParseError, exc:
- out_file.close()
- fail( str( exc ) )
-
- out_file.close()
-
- if g1.skipped > 0:
- print skipped( g1, filedesc=" of 1st dataset" )
- if g2.skipped > 0:
- print skipped( g2, filedesc=" of 2nd dataset" )
- elif g2_copy.skipped > 0:
- print skipped( g2_copy, filedesc=" of 2nd dataset" )
-
-
-if __name__ == "__main__":
- main()
diff --git a/tools/regVariation/featureCounter.xml b/tools/regVariation/featureCounter.xml
deleted file mode 100644
index b85c107815b..00000000000
--- a/tools/regVariation/featureCounter.xml
+++ /dev/null
@@ -1,75 +0,0 @@
-
-
- featureCounter.py $input1 $input2 $output -1 ${input1.metadata.chromCol},${input1.metadata.startCol},${input1.metadata.endCol},${input1.metadata.strandCol} -2 ${input2.metadata.chromCol},${input2.metadata.startCol},${input2.metadata.endCol},${input2.metadata.strandCol}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-.. class:: infomark
-
-**What it does**
-
-This tool finds the coverage of intervals in the first dataset on intervals in the second dataset. The coverage and count are appended as 4 new columns in the resulting dataset.
-
------
-
-**Example**
-
-- If **First dataset** consists of the following windows::
-
- chrX 1 10001 seg 0 -
- chrX 10001 20001 seg 0 -
- chrX 20001 30001 seg 0 -
- chrX 30001 40001 seg 0 -
-
-- and **Second dataset** consists of the following exons::
-
- chrX 5000 6000 seg2 0 -
- chrX 5500 7000 seg2 0 -
- chrX 9000 22000 seg2 0 -
- chrX 24000 34000 seg2 0 -
- chrX 36000 38000 seg2 0 -
-
-- the **Result** is the coverage of exons of the second dataset in each of the windows contained in first dataset::
-
- chrX 1 10001 seg 0 - 3001 0.3001 2 1
- chrX 10001 20001 seg 0 - 10000 1.0 1 0
- chrX 20001 30001 seg 0 - 8000 0.8 0 2
- chrX 30001 40001 seg 0 - 5999 0.5999 1 1
-
-- To clarify, the following line of output ( added columns are indexed by a, b and c )::
-
- a b c d
- chrX 1 10001 seg 0 - 3001 0.3001 2 1
-
- implies that 2 exons (c) fall fully in this window (chrX:1-10001), 1 exon (d) partially overlaps this window, and these 3 exons cover 30.01% (c) of the window size, spanning 3001 nucleotides (a).
-
- * a: number of nucleotides in this window covered by the features in (c) and (d) - features overlapping with each other will be merged to calculate (a)
- * b: fraction of window size covered by features in (c) and (d) - features overlapping with each other will be merged to calculate (b)
- * c: number of features in the 2nd dataset that fall **completely** within this window
- * d: number of features in the 2nd dataset that **partially** overlap this window
-
-
-
diff --git a/tools/regVariation/getIndelRates_3way.py b/tools/regVariation/getIndelRates_3way.py
deleted file mode 100755
index c209e2e1858..00000000000
--- a/tools/regVariation/getIndelRates_3way.py
+++ /dev/null
@@ -1,248 +0,0 @@
-#!/usr/bin/env python
-#Guruprasad Ananda
-
-from galaxy import eggs
-import pkg_resources
-pkg_resources.require( "bx-python" )
-
-import sys, os, tempfile
-import fileinput
-from warnings import warn
-
-from galaxy.tools.util.galaxyops import *
-from bx.intervals.io import *
-
-from bx.intervals.operations import quicksect
-
-def stop_err(msg):
- sys.stderr.write(msg)
- sys.exit()
-
-
-def counter(node, start, end, sort_col):
- global full, blk_len, blk_list
- if node.start < start:
- if node.right:
- counter(node.right, start, end, sort_col)
- elif start <= node.start <= end and start <= node.end <= end:
- full += 1
- if node.other[0] not in blk_list:
- blk_list.append(node.other[0])
- blk_len += int(node.other[sort_col+2])
- if node.left and node.left.maxend > start:
- counter(node.left, start, end, sort_col)
- if node.right:
- counter(node.right, start, end, sort_col)
- elif node.start > end:
- if node.left:
- counter(node.left, start, end, sort_col)
-
-
-infile = sys.argv[1]
-fout = open(sys.argv[2],'w')
-int_file = sys.argv[3]
-if int_file != "None": #User has specified an interval file
- try:
- fint = open(int_file, 'r')
- dbkey_i = sys.argv[4]
- chr_col_i, start_col_i, end_col_i, strand_col_i = parse_cols_arg( sys.argv[5] )
- except:
- stop_err("Unable to open input Interval file")
-
-
-def main():
- for i, line in enumerate( file ( infile )):
- line = line.rstrip('\r\n')
- if len( line )>0 and not line.startswith( '#' ):
- elems = line.split( '\t' )
- break
- if i == 30:
- break # Hopefully we'll never get here...
-
- if len( elems ) != 18:
- stop_err( "This tool only works on tabular data output by 'Fetch Indels from 3-way alignments' tool. The data in your input dataset is either missing or not formatted properly." )
-
- for i, line in enumerate( file ( infile )):
- line = line.rstrip('\r\n')
- elems = line.split('\t')
- try:
- assert int(elems[0])
- assert len(elems) == 18
- if int_file != "None":
- if dbkey_i not in elems[3] and dbkey_i not in elems[8] and dbkey_i not in elems[13]:
- stop_err("The species build corresponding to your interval file is not present in the Indel file.")
- if dbkey_i in elems[3]:
- sort_col = 4
- elif dbkey_i in elems[8]:
- sort_col = 9
- elif dbkey_i in elems[13]:
- sort_col = 14
- else:
- species = []
- species.append( elems[3].split('.')[0] )
- species.append( elems[8].split('.')[0] )
- species.append( elems[13].split('.')[0] )
- sort_col = 0 #Based on block numbers
- break
- except:
- continue
-
- fin = open(infile, 'r')
- skipped = 0
-
- if int_file == "None":
- sorted_infile = tempfile.NamedTemporaryFile()
- cmdline = "sort -n -k"+str(1)+" -o "+sorted_infile.name+" "+infile
- try:
- os.system(cmdline)
- except:
- stop_err("Encountered error while sorting the input file.")
- print >> fout, "#Block\t%s_InsRate\t%s_InsRate\t%s_InsRate\t%s_DelRate\t%s_DelRate\t%s_DelRate" % ( species[0], species[1], species[2], species[0], species[1], species[2] )
- prev_bnum = -1
- sorted_infile.seek(0)
- for line in sorted_infile.readlines():
- line = line.rstrip('\r\n')
- elems = line.split('\t')
- try:
- assert int(elems[0])
- assert len(elems) == 18
- new_bnum = int(elems[0])
- if new_bnum != prev_bnum:
- if prev_bnum != -1:
- irate = []
- drate = []
- for i, elem in enumerate(inserts):
- try:
- irate.append(str("%.2e" % (inserts[i]/blen[i])))
- except:
- irate.append('0')
- try:
- drate.append(str("%.2e" % (deletes[i]/blen[i])))
- except:
- drate.append('0')
- print >> fout, "%s\t%s\t%s" % ( prev_bnum, '\t'.join(irate) , '\t'.join(drate) )
- inserts = [0.0, 0.0, 0.0]
- deletes = [0.0, 0.0, 0.0]
- blen = []
- blen.append( int(elems[6]) )
- blen.append( int(elems[11]) )
- blen.append( int(elems[16]) )
- line_sp = elems[1].split('.')[0]
- sp_ind = species.index(line_sp)
- if elems[1].endswith('insert'):
- inserts[sp_ind] += 1
- elif elems[1].endswith('delete'):
- deletes[sp_ind] += 1
- prev_bnum = new_bnum
- except Exception, ei:
- #print >>sys.stderr, ei
- continue
- irate = []
- drate = []
- for i, elem in enumerate(inserts):
- try:
- irate.append(str("%.2e" % (inserts[i]/blen[i])))
- except:
- irate.append('0')
- try:
- drate.append(str("%.2e" % (deletes[i]/blen[i])))
- except:
- drate.append('0')
- print >> fout, "%s\t%s\t%s" % ( prev_bnum, '\t'.join(irate) , '\t'.join(drate) )
- sys.exit()
-
- inf = open(infile, 'r')
- start_met = False
- end_met = False
- sp_file = tempfile.NamedTemporaryFile()
- for n, line in enumerate(inf):
- line = line.rstrip('\r\n')
- elems = line.split('\t')
- try:
- assert int(elems[0])
- assert len(elems) == 18
- if dbkey_i not in elems[1]:
- if not(start_met):
- continue
- else:
- sp_end = n
- break
- else:
- print >> sp_file, line
- if not(start_met):
- start_met = True
- sp_start = n
- except:
- continue
-
- try:
- assert sp_end
- except:
- sp_end = n+1
-
- sp_file.seek(0)
- win = NiceReaderWrapper( fileinput.FileInput( int_file ),
- chrom_col=chr_col_i,
- start_col=start_col_i,
- end_col=end_col_i,
- strand_col=strand_col_i,
- fix_strand=True)
-
- indel = NiceReaderWrapper( fileinput.FileInput( sp_file.name ),
- chrom_col=1,
- start_col=sort_col,
- end_col=sort_col+1,
- strand_col=-1,
- fix_strand=True)
-
- indelTree = quicksect.IntervalTree()
- for item in indel:
- if type( item ) is GenomicInterval:
- indelTree.insert( item, indel.linenum, item.fields )
- result = []
-
- global full, blk_len, blk_list
- for interval in win:
- if type( interval ) is Header:
- pass
- if type( interval ) is Comment:
- pass
- elif type( interval ) == GenomicInterval:
- chrom = interval.chrom
- start = int(interval.start)
- end = int(interval.end)
- if start > end:
- warn( "Interval start after end!" )
- ins_chr = "%s.%s_insert" % ( dbkey_i, chrom )
- del_chr = "%s.%s_delete" % ( dbkey_i, chrom )
- irate = 0
- drate = 0
- if ins_chr not in indelTree.chroms and del_chr not in indelTree.chroms:
- pass
- else:
- if ins_chr in indelTree.chroms:
- full = 0.0
- blk_len = 0
- blk_list = []
- root = indelTree.chroms[ins_chr] #root node for the chrom insertion tree
- counter(root, start, end, sort_col)
- if blk_len:
- irate = full/blk_len
-
- if del_chr in indelTree.chroms:
- full = 0.0
- blk_len = 0
- blk_list = []
- root = indelTree.chroms[del_chr] #root node for the chrom insertion tree
- counter(root, start, end, sort_col)
- if blk_len:
- drate = full/blk_len
-
- interval.fields.append(str("%.2e" %irate))
- interval.fields.append(str("%.2e" %drate))
- print >> fout, "\t".join(interval.fields)
- fout.flush()
-
-
-if __name__ == "__main__":
- main()
diff --git a/tools/regVariation/getIndelRates_3way.xml b/tools/regVariation/getIndelRates_3way.xml
deleted file mode 100644
index 8f4fe3126ec..00000000000
--- a/tools/regVariation/getIndelRates_3way.xml
+++ /dev/null
@@ -1,61 +0,0 @@
-
- for 3-way alignments
-
- getIndelRates_3way.py $input1 $out_file1
- #if $region.type == "align"
- "None"
- #else
- $region.input2 $input2.dbkey $input2.metadata.chromCol,$input2.metadata.startCol,$input2.metadata.endCol,$input2.metadata.strandCol
- #end if
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-.. class:: infomark
-
-**What it does**
-
-This tool estimates the insertion and deletion rates for alignments in a window of specified size. Rates are computed over the total adjusted lengths (adjusted by disregarding masked bases) of all the alignments blocks from the indel file that fall within that window.
-
------
-
-.. class:: warningmark
-
-**Note**
-
-This tool only works on the output of the 'Estimate Indel Rates for 3-way alignments' tool.
-
-
-
-
-
diff --git a/tools/regVariation/getIndels.py b/tools/regVariation/getIndels.py
deleted file mode 100644
index e9ea3f28230..00000000000
--- a/tools/regVariation/getIndels.py
+++ /dev/null
@@ -1,123 +0,0 @@
-#!/usr/bin/env python
-
-"""
-Estimate INDELs for pair-wise alignments.
-
-usage: %prog maf_input out_file1 out_file2
-"""
-
-from __future__ import division
-from galaxy import eggs
-import pkg_resources
-pkg_resources.require( "bx-python" )
-try:
- pkg_resources.require("numpy")
-except:
- pass
-import sys
-from bx.cookbook import doc_optparse
-from galaxy.tools.exception_handling import *
-import bx.align.maf
-
-assert sys.version_info[:2] >= ( 2, 4 )
-
-def main():
- # Parsing Command Line here
- options, args = doc_optparse.parse( __doc__ )
-
- try:
- inp_file, out_file1 = args
- except:
- print >> sys.stderr, "Tool initialization error."
- sys.exit()
-
- try:
- open(inp_file, 'r')
- except:
- print >> sys.stderr, "Unable to open input file"
- sys.exit()
- try:
- fout1 = open(out_file1, 'w')
- #fout2 = open(out_file2, 'w')
- except:
- print >> sys.stderr, "Unable to open output file"
- sys.exit()
-
- try:
- maf_reader = bx.align.maf.Reader( open(inp_file, 'r') )
- except:
- print >> sys.stderr, "Your MAF file appears to be malformed."
- sys.exit()
-
- print >> fout1, "#Block\tSource\tSeq1_Start\tSeq1_End\tSeq2_Start\tSeq2_End\tIndel_length"
- for block_ind, block in enumerate(maf_reader):
- if len(block.components) < 2:
- continue
- seq1 = block.components[0].text
- src1 = block.components[0].src
- start1 = block.components[0].start
- if len(block.components) == 2:
- seq2 = block.components[1].text
- src2 = block.components[1].src
- start2 = block.components[1].start
- #for pos in range(len(seq1)):
- nt_pos1 = start1-1 #position of the nucleotide (without counting gaps)
- nt_pos2 = start2-1
- pos = 0 #character column position
- gaplen1 = 0
- gaplen2 = 0
- prev_pos_gap1 = 0
- prev_pos_gap2 = 0
- while pos < len(seq1):
- if prev_pos_gap1 == 0:
- gaplen1 = 0
- if prev_pos_gap2 == 0:
- gaplen2 = 0
-
- if seq1[pos] == '-':
- if seq2[pos] != '-':
- nt_pos2 += 1
- gaplen1 += 1
- prev_pos_gap1 = 1
- #write 2
- if prev_pos_gap2 == 1:
- prev_pos_gap2 = 0
- print >> fout1, "%d\t%s\t%s\t%s\t%s\t%s\t%s" % ( block_ind+1, src2, nt_pos1, nt_pos1+1, nt_pos2-1, nt_pos2-1+gaplen2, gaplen2 )
- if pos == len(seq1)-1:
- print >> fout1, "%d\t%s\t%s\t%s\t%s\t%s\t%s" % ( block_ind+1, src1, nt_pos1, nt_pos1+1, nt_pos2+1-gaplen1, nt_pos2+1, gaplen1 )
- else:
- prev_pos_gap1 = 0
- prev_pos_gap2 = 0
- """
- if prev_pos_gap1 == 1:
- prev_pos_gap1 = 0
- print >> fout1, "%d\t%s\t%s\t%s\t%s" % ( block_ind+1, src1, nt_pos1-1, nt_pos1, gaplen1 )
- elif prev_pos_gap2 == 1:
- prev_pos_gap2 = 0
- print >> fout1, "%d\t%s\t%s\t%s\t%s" % ( block_ind+1, src2, nt_pos2-1, nt_pos2, gaplen2 )
- """
- else:
- nt_pos1 += 1
- if seq2[pos] != '-':
- nt_pos2 += 1
- #write both
- if prev_pos_gap1 == 1:
- prev_pos_gap1 = 0
- print >> fout1, "%d\t%s\t%s\t%s\t%s\t%s\t%s" % ( block_ind+1, src1, nt_pos1-1, nt_pos1, nt_pos2-gaplen1, nt_pos2, gaplen1 )
- elif prev_pos_gap2 == 1:
- prev_pos_gap2 = 0
- print >> fout1, "%d\t%s\t%s\t%s\t%s\t%s\t%s" % ( block_ind+1, src2, nt_pos1-gaplen2, nt_pos1, nt_pos2-1, nt_pos2, gaplen2 )
- else:
- gaplen2 += 1
- prev_pos_gap2 = 1
- #write 1
- if prev_pos_gap1 == 1:
- prev_pos_gap1 = 0
- print >> fout1, "%d\t%s\t%s\t%s\t%s\t%s\t%s" % ( block_ind+1, src1, nt_pos1-1, nt_pos1, nt_pos2, nt_pos2+gaplen1, gaplen1 )
- if pos == len(seq1)-1:
- print >> fout1, "%d\t%s\t%s\t%s\t%s\t%s\t%s" % ( block_ind+1, src2, nt_pos1+1-gaplen2, nt_pos1+1, nt_pos2, nt_pos2+1, gaplen2 )
- pos += 1
-
-
-if __name__ == "__main__":
- main()
diff --git a/tools/regVariation/getIndels_2way.xml b/tools/regVariation/getIndels_2way.xml
deleted file mode 100644
index 1d4780ce49a..00000000000
--- a/tools/regVariation/getIndels_2way.xml
+++ /dev/null
@@ -1,59 +0,0 @@
-
- from pairwise alignments
-
- getIndels.py $input1 $out_file1
-
-
-
-
-
-
-
-
-
-
- numpy
-
-
-
-
-
-
-
-
-
-.. class:: infomark
-
-**What it does**
-
-This tool estimates the number of indels for every alignment block of the MAF file.
-
------
-
-.. class:: warningmark
-
-**Note**
-
-Any block/s not containing exactly 2 species will be omitted.
-
------
-
-**Example**
-
-- For the following alignment block::
-
- a score=7233.0
- s hg18.chr1 100 35 + 247249719 AT--GACTGAGGACTTAGTTTAAGATGTTCCTACT
- s rheMac2.chr11 200 31 + 134511895 ATAAG-CGGACGACTTAGTTTAAGATGTTCC----
-
-- running this tool will return::
-
- #Block Source Seq1_Start Seq1_End Seq2_Start Seq2_End Indel_length
- 1 hg18.chr1 101 102 202 204 2
- 1 rheMac2.chr11 103 104 204 205 1
- 1 rheMac2.chr11 129 133 229 230 4
-
-
-
-
-
diff --git a/tools/regVariation/linear_regression.py b/tools/regVariation/linear_regression.py
deleted file mode 100644
index cf7afdc5ced..00000000000
--- a/tools/regVariation/linear_regression.py
+++ /dev/null
@@ -1,147 +0,0 @@
-#!/usr/bin/env python
-
-from galaxy import eggs
-import sys
-from rpy import *
-import numpy
-
-def stop_err(msg):
- sys.stderr.write(msg)
- sys.exit()
-
-infile = sys.argv[1]
-y_col = int(sys.argv[2])-1
-x_cols = sys.argv[3].split(',')
-outfile = sys.argv[4]
-outfile2 = sys.argv[5]
-
-print "Predictor columns: %s; Response column: %d" % ( x_cols, y_col+1 )
-fout = open(outfile,'w')
-elems = []
-for i, line in enumerate( file ( infile )):
- line = line.rstrip('\r\n')
- if len( line )>0 and not line.startswith( '#' ):
- elems = line.split( '\t' )
- break
- if i == 30:
- break # Hopefully we'll never get here...
-
-if len( elems )<1:
- stop_err( "The data in your input dataset is either missing or not formatted properly." )
-
-y_vals = []
-x_vals = []
-
-for k, col in enumerate(x_cols):
- x_cols[k] = int(col)-1
- x_vals.append([])
-
-NA = 'NA'
-for ind, line in enumerate( file( infile )):
- if line and not line.startswith( '#' ):
- try:
- fields = line.split("\t")
- try:
- yval = float(fields[y_col])
- except:
- yval = r('NA')
- y_vals.append(yval)
- for k, col in enumerate(x_cols):
- try:
- xval = float(fields[col])
- except:
- xval = r('NA')
- x_vals[k].append(xval)
- except:
- pass
-
-x_vals1 = numpy.asarray(x_vals).transpose()
-
-dat = r.list(x=array(x_vals1), y=y_vals)
-
-set_default_mode(NO_CONVERSION)
-try:
- linear_model = r.lm(r("y ~ x"), data = r.na_exclude(dat))
-except RException, rex:
- stop_err("Error performing linear regression on the input data.\nEither the response column or one of the predictor columns contain only non-numeric or invalid values.")
-set_default_mode(BASIC_CONVERSION)
-
-coeffs = linear_model.as_py()['coefficients']
-yintercept = coeffs['(Intercept)']
-summary = r.summary(linear_model)
-
-co = summary.get('coefficients', 'NA')
-"""
-if len(co) != len(x_vals)+1:
- stop_err("Stopped performing linear regression on the input data, since one of the predictor columns contains only non-numeric or invalid values.")
-"""
-
-try:
- yintercept = r.round(float(yintercept), digits=10)
- pvaly = r.round(float(co[0][3]), digits=10)
-except:
- pass
-
-print >> fout, "Y-intercept\t%s" % (yintercept)
-print >> fout, "p-value (Y-intercept)\t%s" % (pvaly)
-
-if len(x_vals) == 1: #Simple linear regression case with 1 predictor variable
- try:
- slope = r.round(float(coeffs['x']), digits=10)
- except:
- slope = 'NA'
- try:
- pval = r.round(float(co[1][3]), digits=10)
- except:
- pval = 'NA'
- print >> fout, "Slope (c%d)\t%s" % ( x_cols[0]+1, slope )
- print >> fout, "p-value (c%d)\t%s" % ( x_cols[0]+1, pval )
-else: #Multiple regression case with >1 predictors
- ind = 1
- while ind < len(coeffs.keys()):
- try:
- slope = r.round(float(coeffs['x'+str(ind)]), digits=10)
- except:
- slope = 'NA'
- print >> fout, "Slope (c%d)\t%s" % ( x_cols[ind-1]+1, slope )
- try:
- pval = r.round(float(co[ind][3]), digits=10)
- except:
- pval = 'NA'
- print >> fout, "p-value (c%d)\t%s" % ( x_cols[ind-1]+1, pval )
- ind += 1
-
-rsq = summary.get('r.squared','NA')
-adjrsq = summary.get('adj.r.squared','NA')
-fstat = summary.get('fstatistic','NA')
-sigma = summary.get('sigma','NA')
-
-try:
- rsq = r.round(float(rsq), digits=5)
- adjrsq = r.round(float(adjrsq), digits=5)
- fval = r.round(fstat['value'], digits=5)
- fstat['value'] = str(fval)
- sigma = r.round(float(sigma), digits=10)
-except:
- pass
-
-print >> fout, "R-squared\t%s" % (rsq)
-print >> fout, "Adjusted R-squared\t%s" % (adjrsq)
-print >> fout, "F-statistic\t%s" % (fstat)
-print >> fout, "Sigma\t%s" % (sigma)
-
-r.pdf( outfile2, 8, 8 )
-if len(x_vals) == 1: #Simple linear regression case with 1 predictor variable
- sub_title = "Slope = %s; Y-int = %s" % ( slope, yintercept )
- try:
- r.plot(x=x_vals[0], y=y_vals, xlab="X", ylab="Y", sub=sub_title, main="Scatterplot with regression")
- r.abline(a=yintercept, b=slope, col="red")
- except:
- pass
-else:
- r.pairs(dat, main="Scatterplot Matrix", col="blue")
-try:
- r.plot(linear_model)
-except:
- pass
-r.dev_off()
diff --git a/tools/regVariation/linear_regression.xml b/tools/regVariation/linear_regression.xml
deleted file mode 100644
index d84639e0971..00000000000
--- a/tools/regVariation/linear_regression.xml
+++ /dev/null
@@ -1,71 +0,0 @@
-
-
-
- linear_regression.py
- $input1
- $response_col
- $predictor_cols
- $out_file1
- $out_file2
- 1>/dev/null
-
-
-
-
-
-
-
-
-
-
-
-
-
- rpy
-
-
-
-
-
-
-
-
-
-
-
-
-
-.. class:: infomark
-
-**TIP:** If your data is not TAB delimited, use *Edit Datasets->Convert characters*
-
------
-
-.. class:: infomark
-
-**What it does**
-
-This tool uses the 'lm' function from R statistical package to perform linear regression on the input data. It outputs two files, one containing the summary statistics of the performed regression, and the other containing diagnostic plots to check whether model assumptions are satisfied.
-
-*R Development Core Team (2009). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. ISBN 3-900051-07-0, URL http://www.R-project.org.*
-
------
-
-.. class:: warningmark
-
-**Note**
-
-- This tool currently treats all predictor and response variables as continuous numeric variables. Running the tool on categorical variables might result in incorrect results.
-
-- Rows containing non-numeric (or missing) data in any of the chosen columns will be skipped from the analysis.
-
-- The summary statistics in the output are described below:
-
- - sigma: the square root of the estimated variance of the random error (standard error of the residiuals)
- - R-squared: the fraction of variance explained by the model
- - Adjusted R-squared: the above R-squared statistic adjusted, penalizing for the number of the predictors (p)
- - p-value: p-value for the t-test of the null hypothesis that the corresponding slope is equal to zero against the two-sided alternative.
-
-
-
-
diff --git a/tools/regVariation/logistic_regression_vif.py b/tools/regVariation/logistic_regression_vif.py
deleted file mode 100755
index 68524893aca..00000000000
--- a/tools/regVariation/logistic_regression_vif.py
+++ /dev/null
@@ -1,168 +0,0 @@
-#!/usr/bin/env python
-
-from galaxy import eggs
-import sys
-from rpy import *
-import numpy
-
-def stop_err(msg):
- sys.stderr.write(msg)
- sys.exit()
-
-infile = sys.argv[1]
-y_col = int(sys.argv[2])-1
-x_cols = sys.argv[3].split(',')
-outfile = sys.argv[4]
-
-
-print "Predictor columns: %s; Response column: %d" % ( x_cols, y_col+1 )
-fout = open(outfile,'w')
-elems = []
-for i, line in enumerate( file( infile ) ):
- line = line.rstrip('\r\n')
- if len( line )>0 and not line.startswith( '#' ):
- elems = line.split( '\t' )
- break
- if i == 30:
- break # Hopefully we'll never get here...
-
-if len( elems )<1:
- stop_err( "The data in your input dataset is either missing or not formatted properly." )
-
-y_vals = []
-x_vals = []
-
-for k, col in enumerate(x_cols):
- x_cols[k] = int(col)-1
- x_vals.append([])
-
-NA = 'NA'
-for ind, line in enumerate( file( infile )):
- if line and not line.startswith( '#' ):
- try:
- fields = line.split("\t")
- try:
- yval = float(fields[y_col])
- except:
- yval = r('NA')
- y_vals.append(yval)
- for k, col in enumerate(x_cols):
- try:
- xval = float(fields[col])
- except:
- xval = r('NA')
- x_vals[k].append(xval)
- except:
- pass
-
-x_vals1 = numpy.asarray(x_vals).transpose()
-
-check1 = 0
-check0 = 0
-for i in y_vals:
- if i == 1:
- check1 = 1
- if i == 0:
- check0 = 1
-if check1 == 0 or check0 == 0:
- sys.exit("Warning: logistic regression must have at least two classes")
-
-for i in y_vals:
- if i not in [1, 0, r('NA')]:
- print >> fout, str(i)
- sys.exit("Warning: the current version of this tool can run only with two classes and need to be labeled as 0 and 1.")
-
-dat = r.list(x=array(x_vals1), y=y_vals)
-novif = 0
-set_default_mode(NO_CONVERSION)
-try:
- linear_model = r.glm(r("y ~ x"), data=r.na_exclude(dat), family="binomial")
-except RException, rex:
- stop_err("Error performing logistic regression on the input data.\nEither the response column or one of the predictor columns contain only non-numeric or invalid values.")
-if len(x_cols)>1:
- try:
- r('suppressPackageStartupMessages(library(car))')
- r.assign('dat', dat)
- r.assign('ncols', len(x_cols))
- vif = r.vif(r('glm(dat$y ~ ., data = na.exclude(data.frame(as.matrix(dat$x,ncol=ncols))->datx), family="binomial")'))
- except RException, rex:
- print rex
-else:
- novif = 1
-
-set_default_mode(BASIC_CONVERSION)
-
-coeffs = linear_model.as_py()['coefficients']
-null_deviance = linear_model.as_py()['null.deviance']
-residual_deviance = linear_model.as_py()['deviance']
-yintercept = coeffs['(Intercept)']
-summary = r.summary(linear_model)
-co = summary.get('coefficients', 'NA')
-"""
-if len(co) != len(x_vals)+1:
- stop_err("Stopped performing logistic regression on the input data, since one of the predictor columns contains only non-numeric or invalid values.")
-"""
-
-try:
- yintercept = r.round(float(yintercept), digits=10)
- pvaly = r.round(float(co[0][3]), digits=10)
-except:
- pass
-print >> fout, "response column\tc%d" % (y_col+1)
-tempP = []
-for i in x_cols:
- tempP.append('c'+str(i+1))
-tempP = ','.join(tempP)
-print >> fout, "predictor column(s)\t%s" % (tempP)
-print >> fout, "Y-intercept\t%s" % (yintercept)
-print >> fout, "p-value (Y-intercept)\t%s" % (pvaly)
-
-if len(x_vals) == 1: #Simple linear regression case with 1 predictor variable
- try:
- slope = r.round(float(coeffs['x']), digits=10)
- except:
- slope = 'NA'
- try:
- pval = r.round(float(co[1][3]), digits=10)
- except:
- pval = 'NA'
- print >> fout, "Slope (c%d)\t%s" % ( x_cols[0]+1, slope )
- print >> fout, "p-value (c%d)\t%s" % ( x_cols[0]+1, pval )
-else: #Multiple regression case with >1 predictors
- ind = 1
- while ind < len(coeffs.keys()):
- try:
- slope = r.round(float(coeffs['x'+str(ind)]), digits=10)
- except:
- slope = 'NA'
- print >> fout, "Slope (c%d)\t%s" % ( x_cols[ind-1]+1, slope )
- try:
- pval = r.round(float(co[ind][3]), digits=10)
- except:
- pval = 'NA'
- print >> fout, "p-value (c%d)\t%s" % ( x_cols[ind-1]+1, pval )
- ind += 1
-
-rsq = summary.get('r.squared','NA')
-
-try:
- rsq = r.round(float((null_deviance-residual_deviance)/null_deviance), digits=5)
- null_deviance = r.round(float(null_deviance), digits=5)
- residual_deviance = r.round(float(residual_deviance), digits=5)
-except:
- pass
-
-print >> fout, "Null deviance\t%s" % (null_deviance)
-print >> fout, "Residual deviance\t%s" % (residual_deviance)
-print >> fout, "pseudo R-squared\t%s" % (rsq)
-print >> fout, "\n"
-print >> fout, 'vif'
-
-if novif == 0:
- py_vif = vif.as_py()
- count = 0
- for i in sorted(py_vif.keys()):
- print >> fout, 'c'+str(x_cols[count]+1), str(py_vif[i])
- count += 1
-elif novif == 1:
- print >> fout, "vif can calculate only when model have more than 1 predictor"
diff --git a/tools/regVariation/logistic_regression_vif.xml b/tools/regVariation/logistic_regression_vif.xml
deleted file mode 100755
index 2b491bd1aac..00000000000
--- a/tools/regVariation/logistic_regression_vif.xml
+++ /dev/null
@@ -1,74 +0,0 @@
-
-
-
- logistic_regression_vif.py
- $input1
- $response_col
- $predictor_cols
- $out_file1
- 1>/dev/null
-
-
-
-
-
-
-
-
-
-
-
-
-
- rpy
-
-
-
-
-
-
-
-
-
-
-
-
-
-.. class:: infomark
-
-**TIP:** If your data is not TAB delimited, use *Edit Datasets->Convert characters*
-
------
-
-.. class:: infomark
-
-**What it does**
-
-This tool uses the **'glm'** function from R statistical package to perform logistic regression on the input data. It outputs one file containing the summary statistics of the performed regression. Also, it calculates VIF(Variance Inflation Factor) with **'vif'** function from library (car) in R.
-
-
-*R Development Core Team (2010). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. ISBN 3-900051-07-0, URL http://www.R-project.org.*
-
------
-
-.. class:: warningmark
-
-**Note**
-
-- This tool currently treats all predictor variables as continuous numeric variables and response variable as categorical variable. Currently, the response variable can have only two classes, namely 0 and 1. The program will take 0 as base class.
-
-- Rows containing non-numeric (or missing) data in any of the chosen columns will be skipped from the analysis.
-
-- The summary statistics in the output are described below:
-
-- Pseudo R-squared: the proportion of model improvement from null model
-- p-value: p-value for the z-test of the null hypothesis that the corresponding slope is equal to zero against the two-sided alternative.
-- Coefficient indicates log ratio of (probability to be class 1 / probability to be class 0)
-
-- This tool also provides **Variance Inflation Factor or VIF** which quantifies the level of multicollinearity. The tool will automatic generate VIF if the model has more than one predictor. The higher the VIF, the higher is the multicollinearity. Multicollinearity will inflate standard error and reduce level of significance of the predictor. In the worst case, it can reverse direction of slope for highly correlated predictors if one of them is significant. A general thumb-rule is to use those predictors having VIF lower than 10 or 5.
-- **vif** is calculated by
- - First, regressing each predictor over all other predictors, and recording R-squared for each regression.
- - Second, computing vif as 1/(1- R_squared)
-
-
-
diff --git a/tools/regVariation/maf_cpg_filter.py b/tools/regVariation/maf_cpg_filter.py
deleted file mode 100644
index 0e3e4cb0b5d..00000000000
--- a/tools/regVariation/maf_cpg_filter.py
+++ /dev/null
@@ -1,60 +0,0 @@
-#!/usr/bin/env python
-#Guruprasad Ananda
-#Adapted from bx/scripts/maf_mask_cpg.py
-"""
-Mask out potential CpG sites from a maf. Restricted or inclusive definition
-of CpG sites can be used. The total fraction masked is printed to stderr.
-
-usage: %prog < input > output restricted
- -m, --mask=N: Character to use as mask ('?' is default)
-"""
-
-from galaxy import eggs
-import pkg_resources
-pkg_resources.require( "bx-python" )
-try:
- pkg_resources.require( "numpy" )
-except:
- pass
-import bx.align
-import bx.align.maf
-from bx.cookbook import doc_optparse
-import sys
-import bx.align.sitemask.cpg
-
-assert sys.version_info[:2] >= ( 2, 4 )
-
-def main():
- options, args = doc_optparse.parse( __doc__ )
- try:
- inp_file, out_file, sitetype, definition = args
- if options.mask:
- mask = int(options.mask)
- else:
- mask = 0
- except:
- print >> sys.stderr, "Tool initialization error."
- sys.exit()
-
- reader = bx.align.maf.Reader( open(inp_file, 'r') )
- writer = bx.align.maf.Writer( open(out_file,'w') )
-
- mask_chr_dict = {0:'#', 1:'$', 2:'^', 3:'*', 4:'?', 5:'N'}
- mask = mask_chr_dict[mask]
-
- if sitetype == "CpG":
- if int(definition) == 1:
- cpgfilter = bx.align.sitemask.cpg.Restricted( mask=mask )
- defn = "CpG-Restricted"
- else:
- cpgfilter = bx.align.sitemask.cpg.Inclusive( mask=mask )
- defn = "CpG-Inclusive"
- else:
- cpgfilter = bx.align.sitemask.cpg.nonCpG( mask=mask )
- defn = "non-CpG"
- cpgfilter.run( reader, writer.write )
-
- print "%2.2f percent bases masked; Mask character = %s, Definition = %s" % ( float(cpgfilter.masked)/float(cpgfilter.total) * 100, mask, defn )
-
-if __name__ == "__main__":
- main()
diff --git a/tools/regVariation/maf_cpg_filter.xml b/tools/regVariation/maf_cpg_filter.xml
deleted file mode 100644
index 7d0d51ceda5..00000000000
--- a/tools/regVariation/maf_cpg_filter.xml
+++ /dev/null
@@ -1,87 +0,0 @@
-
- from MAF file
-
- maf_cpg_filter.py
- $input
- $out_file1
- $masksite.type
- #if $masksite.type == "CpG":
- $masksite.definition
- #else:
- "NA"
- #end if
- -m $mask_char
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
- numpy
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-.. class:: infomark
-
-**What it does**
-
-This tool takes a MAF file as input and masks CpG sites in every alignment block of the MAF file.
-
------
-
-.. class:: warningmark
-
-**Note**
-
-*Inclusive definition* defines CpG sites as those sites that are CG in at least one of the species.
-
-*Restricted definition* considers sites to be CpG if they are CG in at least one of the species, however, sites that are part of overlapping CpGs are excluded.
-
-For more information on CpG site definitions, please refer this article_.
-
-.. _article: http://mbe.oxfordjournals.org/cgi/content/full/23/3/565
-
-
-
-
-
diff --git a/tools/regVariation/microsats_alignment_level.py b/tools/regVariation/microsats_alignment_level.py
deleted file mode 100644
index 617261e7809..00000000000
--- a/tools/regVariation/microsats_alignment_level.py
+++ /dev/null
@@ -1,318 +0,0 @@
- #!/usr/bin/env python
-#Guruprasad Ananda
-"""
-Uses SPUTNIK to fetch microsatellites and extracts orthologous repeats from the sputnik output.
-"""
-from galaxy import eggs
-import os
-import re
-import string
-import sys
-import tempfile
-
-def reverse_complement(text):
- DNA_COMP = string.maketrans( "ACGTacgt", "TGCAtgca" )
- comp = [ch for ch in text.translate(DNA_COMP)]
- comp.reverse()
- return "".join(comp)
-
-
-def main():
- if len(sys.argv) != 8:
- print >> sys.stderr, "Insufficient number of arguments."
- sys.exit()
-
- infile = open(sys.argv[1],'r')
- separation = int(sys.argv[2])
- outfile = sys.argv[3]
- mono_threshold = int(sys.argv[5])
- non_mono_threshold = int(sys.argv[6])
- allow_different_units = int(sys.argv[7])
-
- print "Min distance = %d bp; Min threshold for mono repeats = %d; Min threshold for non-mono repeats = %d; Allow different motifs = %s" % ( separation, mono_threshold, non_mono_threshold, allow_different_units==1 )
- try:
- fout = open(outfile, "w")
- print >> fout, "#Block\tSeq1_Name\tSeq1_Start\tSeq1_End\tSeq1_Type\tSeq1_Length\tSeq1_RepeatNumber\tSeq1_Unit\tSeq2_Name\tSeq2_Start\tSeq2_End\tSeq2_Type\tSeq2_Length\tSeq2_RepeatNumber\tSeq2_Unit"
- #sputnik_cmd = os.path.join(os.path.split(sys.argv[0])[0], "sputnik")
- sputnik_cmd = "sputnik"
- input = infile.read()
- block_num = 0
- input = input.replace('\r','\n')
- for block in input.split('\n\n'):
- block_num += 1
- tmpin = tempfile.NamedTemporaryFile()
- tmpout = tempfile.NamedTemporaryFile()
- tmpin.write(block.strip())
- cmdline = sputnik_cmd + " " + tmpin.name + " > /dev/null 2>&1 >> " + tmpout.name
- try:
- os.system(cmdline)
- except Exception:
- continue
- sputnik_out = tmpout.read()
- tmpin.close()
- tmpout.close()
- if sputnik_out != "":
- if len(block.split('>')[1:]) != 2: #len(sputnik_out.split('>')):
- continue
- align_block = block.strip().split('>')
-
- lendict = {'mononucleotide':1, 'dinucleotide':2, 'trinucleotide':3, 'tetranucleotide':4, 'pentanucleotide':5, 'hexanucleotide':6}
- blockdict = {}
- r = 0
- namelist = []
- for k, sput_block in enumerate(sputnik_out.split('>')[1:]):
- whole_seq = ''.join(align_block[k+1].split('\n')[1:]).replace('\n','').strip()
- p = re.compile('\n(\S*nucleotide)')
- repeats = p.split(sput_block.strip())
- repeats_count = len(repeats)
- j = 1
- name = repeats[0].strip()
- try:
- coords = re.search('\d+[-_:]\d+', name).group()
- coords = coords.replace('_', '-').replace(':', '-')
- except Exception:
- coords = '0-0'
- r += 1
- blockdict[r] = {}
- try:
- sp_name = name[:name.index('.')]
- chr_name = name[name.index('.'):name.index('(')]
- namelist.append(sp_name + chr_name)
- except:
- namelist.append(name[:20])
- while j < repeats_count:
- try:
- if repeats[j].strip() not in lendict:
- j += 2
- continue
-
- if blockdict[r].has_key('types'):
- blockdict[r]['types'].append(repeats[j].strip()) #type of microsat
- else:
- blockdict[r]['types'] = [repeats[j].strip()] #type of microsat
-
- start = int(repeats[j+1].split('--')[0].split(':')[0].strip())
- #check to see if there are gaps before the start of the repeat, and change the start accordingly
- sgaps = 0
- ch_pos = start - 1
- while ch_pos >= 0:
- if whole_seq[ch_pos] == '-':
- sgaps += 1
- else:
- break #break at the 1st non-gap character
- ch_pos -= 1
- if blockdict[r].has_key('starts'):
- blockdict[r]['starts'].append(start+sgaps) #start co-ords adjusted with alignment co-ords to include GAPS
- else:
- blockdict[r]['starts'] = [start+sgaps]
-
- end = int(repeats[j+1].split('--')[0].split(':')[1].strip())
- #check to see if there are gaps after the end of the repeat, and change the end accordingly
- egaps = 0
- for ch in whole_seq[end:]:
- if ch == '-':
- egaps += 1
- else:
- break #break at the 1st non-gap character
- if blockdict[r].has_key('ends'):
- blockdict[r]['ends'].append(end+egaps) #end co-ords adjusted with alignment co-ords to include GAPS
- else:
- blockdict[r]['ends'] = [end+egaps]
-
- repeat_seq = ''.join(repeats[j+1].replace('\r','\n').split('\n')[1:]).strip() #Repeat Sequence
- repeat_len = repeats[j+1].split('--')[1].split()[1].strip()
- gap_count = repeat_seq.count('-')
- #print repeats[j+1].split('--')[1], len(repeat_seq), repeat_len, gap_count
- repeat_len = str(int(repeat_len) - gap_count)
-
- rel_start = blockdict[r]['starts'][-1]
- gaps_before_start = whole_seq[:rel_start].count('-')
-
- if blockdict[r].has_key('gaps_before_start'):
- blockdict[r]['gaps_before_start'].append(gaps_before_start) #lengths
- else:
- blockdict[r]['gaps_before_start'] = [gaps_before_start] #lengths
-
- whole_seq_start = int(coords.split('-')[0])
- if blockdict[r].has_key('whole_seq_start'):
- blockdict[r]['whole_seq_start'].append(whole_seq_start) #lengths
- else:
- blockdict[r]['whole_seq_start'] = [whole_seq_start] #lengths
-
- if blockdict[r].has_key('lengths'):
- blockdict[r]['lengths'].append(repeat_len) #lengths
- else:
- blockdict[r]['lengths'] = [repeat_len] #lengths
-
- if blockdict[r].has_key('counts'):
- blockdict[r]['counts'].append(str(int(repeat_len)/lendict[repeats[j].strip()])) #Repeat Unit
- else:
- blockdict[r]['counts'] = [str(int(repeat_len)/lendict[repeats[j].strip()])] #Repeat Unit
-
- if blockdict[r].has_key('units'):
- blockdict[r]['units'].append(repeat_seq[:lendict[repeats[j].strip()]]) #Repeat Unit
- else:
- blockdict[r]['units'] = [repeat_seq[:lendict[repeats[j].strip()]]] #Repeat Unit
-
- except Exception:
- pass
- j += 2
- #check the co-ords of all repeats corresponding to a sequence and remove adjacent repeats separated by less than the user-specified 'separation'.
- delete_index_list = []
- for ind, item in enumerate(blockdict[r]['ends']):
- try:
- if blockdict[r]['starts'][ind+1]-item < separation:
- if ind not in delete_index_list:
- delete_index_list.append(ind)
- if ind+1 not in delete_index_list:
- delete_index_list.append(ind+1)
- except Exception:
- pass
- for index in delete_index_list: #mark them for deletion
- try:
- blockdict[r]['starts'][index] = 'marked'
- blockdict[r]['ends'][index] = 'marked'
- blockdict[r]['types'][index] = 'marked'
- blockdict[r]['gaps_before_start'][index] = 'marked'
- blockdict[r]['whole_seq_start'][index] = 'marked'
- blockdict[r]['lengths'][index] = 'marked'
- blockdict[r]['counts'][index] = 'marked'
- blockdict[r]['units'][index] = 'marked'
- except Exception:
- pass
- #remove 'marked' elements from all the lists
- """
- for key in blockdict[r].keys():
- for elem in blockdict[r][key]:
- if elem == 'marked':
- blockdict[r][key].remove(elem)
- """
- #print blockdict
-
- #make sure that the blockdict has keys for both the species
- if (1 not in blockdict) or (2 not in blockdict):
- continue
-
- visited_2 = [0 for x in range(len(blockdict[2]['starts']))]
- for ind1, coord_s1 in enumerate(blockdict[1]['starts']):
- if coord_s1 == 'marked':
- continue
- coord_e1 = blockdict[1]['ends'][ind1]
- out = []
- for ind2, coord_s2 in enumerate(blockdict[2]['starts']):
- if coord_s2 == 'marked':
- visited_2[ind2] = 1
- continue
- coord_e2 = blockdict[2]['ends'][ind2]
- #skip if the 2 repeats are not of the same type or don't have the same repeating unit.
- if allow_different_units == 0:
- if (blockdict[1]['types'][ind1] != blockdict[2]['types'][ind2]):
- continue
- else:
- if (blockdict[1]['units'][ind1] not in blockdict[2]['units'][ind2]*2) and (reverse_complement(blockdict[1]['units'][ind1]) not in blockdict[2]['units'][ind2]*2):
- continue
- #print >> sys.stderr, (reverse_complement(blockdict[1]['units'][ind1]) not in blockdict[2]['units'][ind2]*2)
- #skip if the repeat number thresholds are not met
- if blockdict[1]['types'][ind1] == 'mononucleotide':
- if (int(blockdict[1]['counts'][ind1]) < mono_threshold):
- continue
- else:
- if (int(blockdict[1]['counts'][ind1]) < non_mono_threshold):
- continue
-
- if blockdict[2]['types'][ind2] == 'mononucleotide':
- if (int(blockdict[2]['counts'][ind2]) < mono_threshold):
- continue
- else:
- if (int(blockdict[2]['counts'][ind2]) < non_mono_threshold):
- continue
- #print "s1,e1=%s,%s; s2,e2=%s,%s" % ( coord_s1, coord_e1, coord_s2, coord_e2 )
- if (coord_s1 in range(coord_s2, coord_e2)) or (coord_e1 in range(coord_s2, coord_e2)):
- out.append(str(block_num))
- out.append(namelist[0])
- rel_start = blockdict[1]['whole_seq_start'][ind1] + coord_s1 - blockdict[1]['gaps_before_start'][ind1]
- rel_end = rel_start + int(blockdict[1]['lengths'][ind1])
- out.append(str(rel_start))
- out.append(str(rel_end))
- out.append(blockdict[1]['types'][ind1])
- out.append(blockdict[1]['lengths'][ind1])
- out.append(blockdict[1]['counts'][ind1])
- out.append(blockdict[1]['units'][ind1])
- out.append(namelist[1])
- rel_start = blockdict[2]['whole_seq_start'][ind2] + coord_s2 - blockdict[2]['gaps_before_start'][ind2]
- rel_end = rel_start + int(blockdict[2]['lengths'][ind2])
- out.append(str(rel_start))
- out.append(str(rel_end))
- out.append(blockdict[2]['types'][ind2])
- out.append(blockdict[2]['lengths'][ind2])
- out.append(blockdict[2]['counts'][ind2])
- out.append(blockdict[2]['units'][ind2])
- print >> fout, '\t'.join(out)
- visited_2[ind2] = 1
- out = []
-
- if 0 in visited_2: #there are still some elements in 2nd set which haven't found orthologs yet.
- for ind2, coord_s2 in enumerate(blockdict[2]['starts']):
- if coord_s2 == 'marked':
- continue
- if visited_2[ind] != 0:
- continue
- coord_e2 = blockdict[2]['ends'][ind2]
- out = []
- for ind1, coord_s1 in enumerate(blockdict[1]['starts']):
- if coord_s1 == 'marked':
- continue
- coord_e1 = blockdict[1]['ends'][ind1]
- #skip if the 2 repeats are not of the same type or don't have the same repeating unit.
- if allow_different_units == 0:
- if (blockdict[1]['types'][ind1] != blockdict[2]['types'][ind2]):
- continue
- else:
- if (blockdict[1]['units'][ind1] not in blockdict[2]['units'][ind2]*2):# and reverse_complement(blockdict[1]['units'][ind1]) not in blockdict[2]['units'][ind2]*2:
- continue
- #skip if the repeat number thresholds are not met
- if blockdict[1]['types'][ind1] == 'mononucleotide':
- if (int(blockdict[1]['counts'][ind1]) < mono_threshold):
- continue
- else:
- if (int(blockdict[1]['counts'][ind1]) < non_mono_threshold):
- continue
-
- if blockdict[2]['types'][ind2] == 'mononucleotide':
- if (int(blockdict[2]['counts'][ind2]) < mono_threshold):
- continue
- else:
- if (int(blockdict[2]['counts'][ind2]) < non_mono_threshold):
- continue
-
- if (coord_s2 in range(coord_s1, coord_e1)) or (coord_e2 in range(coord_s1, coord_e1)):
- out.append(str(block_num))
- out.append(namelist[0])
- rel_start = blockdict[1]['whole_seq_start'][ind1] + coord_s1 - blockdict[1]['gaps_before_start'][ind1]
- rel_end = rel_start + int(blockdict[1]['lengths'][ind1])
- out.append(str(rel_start))
- out.append(str(rel_end))
- out.append(blockdict[1]['types'][ind1])
- out.append(blockdict[1]['lengths'][ind1])
- out.append(blockdict[1]['counts'][ind1])
- out.append(blockdict[1]['units'][ind1])
- out.append(namelist[1])
- rel_start = blockdict[2]['whole_seq_start'][ind2] + coord_s2 - blockdict[2]['gaps_before_start'][ind2]
- rel_end = rel_start + int(blockdict[2]['lengths'][ind2])
- out.append(str(rel_start))
- out.append(str(rel_end))
- out.append(blockdict[2]['types'][ind2])
- out.append(blockdict[2]['lengths'][ind2])
- out.append(blockdict[2]['counts'][ind2])
- out.append(blockdict[2]['units'][ind2])
- print >> fout, '\t'.join(out)
- visited_2[ind2] = 1
- out = []
-
- #print >> fout, blockdict
- except Exception, exc:
- print >> sys.stderr, "type(exc),args,exc: %s, %s, %s" % ( type(exc), exc.args, exc )
-
-
-if __name__ == "__main__":
- main()
diff --git a/tools/regVariation/microsats_alignment_level.xml b/tools/regVariation/microsats_alignment_level.xml
deleted file mode 100644
index 949f79ca440..00000000000
--- a/tools/regVariation/microsats_alignment_level.xml
+++ /dev/null
@@ -1,61 +0,0 @@
-
- from pair-wise alignments
-
- microsats_alignment_level.py $input1 $separation $out_file1 "2way" $mono_threshold $non_mono_threshold $allow_different_units
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
- sputnik
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-.. class:: infomark
-
-**What it does**
-
-This tool uses a modified version of SPUTNIK to fetch microsatellite repeats from the input fasta sequences and extracts orthologous repeats from the sputnik output. The modified version allows detection of mononucleotide microsatellites. More information on SPUTNIK can be found on this website_. The modified version is available here_.
-
------
-
-.. class:: warningmark
-
-**Note**
-
-- Any block/s not containing exactly 2 species will be omitted.
-
-- This tool will filter out microsatellites based on the user input values for minimum distance and repeat number thresholds. Further, this tool will also filter out microsatellites that have no orthologous microsatellites in one of the species.
-
-.. _website: http://espressosoftware.com/pages/sputnik.jsp
-.. _here: http://www.bx.psu.edu/svn/universe/dependencies/sputnik/
-
-
-
-
diff --git a/tools/regVariation/microsats_mutability.py b/tools/regVariation/microsats_mutability.py
deleted file mode 100644
index e99101c896c..00000000000
--- a/tools/regVariation/microsats_mutability.py
+++ /dev/null
@@ -1,495 +0,0 @@
-#!/usr/bin/env python
-#Guruprasad Ananda
-"""
-This tool computes microsatellite mutability for the orthologous microsatellites fetched from 'Extract Orthologous Microsatellites from pair-wise alignments' tool.
-"""
-from galaxy import eggs
-import fileinput
-import string
-import sys
-import tempfile
-from galaxy.tools.util.galaxyops import *
-from bx.intervals.io import *
-from bx.intervals.operations import quicksect
-
-fout = open(sys.argv[2],'w')
-p_group = int(sys.argv[3]) #primary "group-by" feature
-p_bin_size = int(sys.argv[4])
-s_group = int(sys.argv[5]) #sub-group by feature
-s_bin_size = int(sys.argv[6])
-mono_threshold = 9
-non_mono_threshold = 4
-p_group_cols = [p_group, p_group+7]
-s_group_cols = [s_group, s_group+7]
-num_generations = int(sys.argv[7])
-region = sys.argv[8]
-int_file = sys.argv[9]
-if int_file != "None": #User has specified an interval file
- try:
- fint = open(int_file, 'r')
- dbkey_i = sys.argv[10]
- chr_col_i, start_col_i, end_col_i, strand_col_i = parse_cols_arg( sys.argv[11] )
- except:
- stop_err("Unable to open input Interval file")
-
-
-def stop_err(msg):
- sys.stderr.write(msg)
- sys.exit()
-
-
-def reverse_complement(text):
- DNA_COMP = string.maketrans( "ACGTacgt", "TGCAtgca" )
- comp = [ch for ch in text.translate(DNA_COMP)]
- comp.reverse()
- return "".join(comp)
-
-
-def get_unique_elems(elems):
- seen = set()
- return[x for x in elems if x not in seen and not seen.add(x)]
-
-
-def get_binned_lists(uniqlist, binsize):
- binnedlist = []
- uniqlist.sort()
- start = int(uniqlist[0])
- bin_ind = 0
- l_ind = 0
- binnedlist.append([])
- while l_ind < len(uniqlist):
- elem = int(uniqlist[l_ind])
- if elem in range(start, start+binsize):
- binnedlist[bin_ind].append(elem)
- else:
- start += binsize
- bin_ind += 1
- binnedlist.append([])
- binnedlist[bin_ind].append(elem)
- l_ind += 1
- return binnedlist
-
-
-def fetch_weight(H, C, t):
- if (H-(C-H)) < t:
- return 2.0
- else:
- return 1.0
-
-
-def mutabilityEstimator(repeats1, repeats2, thresholds):
- mut_num = 0.0 #Mutability Numerator
- mut_den = 0.0 #Mutability denominator
- for ind, H in enumerate(repeats1):
- C = repeats2[ind]
- t = thresholds[ind]
- w = fetch_weight(H, C, t)
- mut_num += ((H-C)*(H-C)*w)
- mut_den += w
- return [mut_num, mut_den]
-
-
-def output_writer(blk, blk_lines):
- global winspecies, speciesind
- all_elems_1 = []
- all_elems_2 = []
- all_s_elems_1 = []
- all_s_elems_2 = []
- for bline in blk_lines:
- if not(bline):
- continue
- items = bline.split('\t')
- seq1 = items[1]
- seq2 = items[8]
- if p_group_cols[0] == 6:
- items[p_group_cols[0]] = int(items[p_group_cols[0]])
- items[p_group_cols[1]] = int(items[p_group_cols[1]])
- if s_group_cols[0] == 6:
- items[s_group_cols[0]] = int(items[s_group_cols[0]])
- items[s_group_cols[1]] = int(items[s_group_cols[1]])
- all_elems_1.append(items[p_group_cols[0]]) #primary col elements for species 1
- all_elems_2.append(items[p_group_cols[1]]) #primary col elements for species 2
- if s_group_cols[0] != -1: #sub-group is not None
- all_s_elems_1.append(items[s_group_cols[0]]) #secondary col elements for species 1
- all_s_elems_2.append(items[s_group_cols[1]]) #secondary col elements for species 2
- uniq_elems_1 = get_unique_elems(all_elems_1)
- uniq_elems_2 = get_unique_elems(all_elems_2)
- if s_group_cols[0] != -1:
- uniq_s_elems_1 = get_unique_elems(all_s_elems_1)
- uniq_s_elems_2 = get_unique_elems(all_s_elems_2)
- mut1 = {}
- mut2 = {}
- count1 = {}
- count2 = {}
- """
- if p_group_cols[0] == 7: #i.e. the option chosen is group-by unit(AG, GTC, etc)
- uniq_elems_1 = get_unique_units(j.sort(lambda x, y: len(x)-len(y)))
- """
- if p_group_cols[0] == 6: #i.e. the option chosen is group-by repeat number.
- uniq_elems_1 = get_binned_lists( uniq_elems_1, p_bin_size )
- uniq_elems_2 = get_binned_lists( uniq_elems_2, p_bin_size )
-
- if s_group_cols[0] == 6: #i.e. the option chosen is subgroup-by repeat number.
- uniq_s_elems_1 = get_binned_lists( uniq_s_elems_1, s_bin_size )
- uniq_s_elems_2 = get_binned_lists( uniq_s_elems_2, s_bin_size )
-
- for pitem1 in uniq_elems_1:
- #repeats1 = []
- #repeats2 = []
- thresholds = []
- if s_group_cols[0] != -1: #Sub-group by feature is not None
- for sitem1 in uniq_s_elems_1:
- repeats1 = []
- repeats2 = []
- if type(sitem1) == type(''):
- sitem1 = sitem1.strip()
- for bline in blk_lines:
- belems = bline.split('\t')
- if type(pitem1) == list:
- if p_group_cols[0] == 6:
- belems[p_group_cols[0]] = int(belems[p_group_cols[0]])
- if belems[p_group_cols[0]] in pitem1:
- if belems[s_group_cols[0]] == sitem1:
- repeats1.append(int(belems[6]))
- repeats2.append(int(belems[13]))
- if belems[4] == 'mononucleotide':
- thresholds.append(mono_threshold)
- else:
- thresholds.append(non_mono_threshold)
- mut1[str(pitem1)+'\t'+str(sitem1)] = mutabilityEstimator( repeats1, repeats2, thresholds )
- if region == 'align':
- count1[str(pitem1)+'\t'+str(sitem1)] = min( sum(repeats1), sum(repeats2) )
- else:
- if winspecies == 1:
- count1["%s\t%s" % ( pitem1, sitem1 )] = sum(repeats1)
- elif winspecies == 2:
- count1["%s\t%s" % ( pitem1, sitem1 )] = sum(repeats2)
- else:
- if type(sitem1) == list:
- if s_group_cols[0] == 6:
- belems[s_group_cols[0]] = int(belems[s_group_cols[0]])
- if belems[p_group_cols[0]] == pitem1 and belems[s_group_cols[0]] in sitem1:
- repeats1.append(int(belems[6]))
- repeats2.append(int(belems[13]))
- if belems[4] == 'mononucleotide':
- thresholds.append(mono_threshold)
- else:
- thresholds.append(non_mono_threshold)
- mut1["%s\t%s" % ( pitem1, sitem1 )] = mutabilityEstimator( repeats1, repeats2, thresholds )
- if region == 'align':
- count1[str(pitem1)+'\t'+str(sitem1)] = min( sum(repeats1), sum(repeats2) )
- else:
- if winspecies == 1:
- count1[str(pitem1)+'\t'+str(sitem1)] = sum(repeats1)
- elif winspecies == 2:
- count1[str(pitem1)+'\t'+str(sitem1)] = sum(repeats2)
- else:
- if belems[p_group_cols[0]] == pitem1 and belems[s_group_cols[0]] == sitem1:
- repeats1.append(int(belems[6]))
- repeats2.append(int(belems[13]))
- if belems[4] == 'mononucleotide':
- thresholds.append(mono_threshold)
- else:
- thresholds.append(non_mono_threshold)
- mut1["%s\t%s" % ( pitem1, sitem1 )] = mutabilityEstimator( repeats1, repeats2, thresholds )
- if region == 'align':
- count1[str(pitem1)+'\t'+str(sitem1)] = min( sum(repeats1), sum(repeats2) )
- else:
- if winspecies == 1:
- count1["%s\t%s" % ( pitem1, sitem1 )] = sum(repeats1)
- elif winspecies == 2:
- count1["%s\t%s" % ( pitem1, sitem1 )] = sum(repeats2)
- else: #Sub-group by feature is None
- for bline in blk_lines:
- belems = bline.split('\t')
- if type(pitem1) == list:
- #print >> sys.stderr, "item: " + str(item1)
- if p_group_cols[0] == 6:
- belems[p_group_cols[0]] = int(belems[p_group_cols[0]])
- if belems[p_group_cols[0]] in pitem1:
- repeats1.append(int(belems[6]))
- repeats2.append(int(belems[13]))
- if belems[4] == 'mononucleotide':
- thresholds.append(mono_threshold)
- else:
- thresholds.append(non_mono_threshold)
- else:
- if belems[p_group_cols[0]] == pitem1:
- repeats1.append(int(belems[6]))
- repeats2.append(int(belems[13]))
- if belems[4] == 'mononucleotide':
- thresholds.append(mono_threshold)
- else:
- thresholds.append(non_mono_threshold)
- mut1["%s" % (pitem1)] = mutabilityEstimator( repeats1, repeats2, thresholds )
- if region == 'align':
- count1["%s" % (pitem1)] = min( sum(repeats1), sum(repeats2) )
- else:
- if winspecies == 1:
- count1[str(pitem1)] = sum(repeats1)
- elif winspecies == 2:
- count1[str(pitem1)] = sum(repeats2)
-
- for pitem2 in uniq_elems_2:
- #repeats1 = []
- #repeats2 = []
- thresholds = []
- if s_group_cols[0] != -1: #Sub-group by feature is not None
- for sitem2 in uniq_s_elems_2:
- repeats1 = []
- repeats2 = []
- if type(sitem2)==type(''):
- sitem2 = sitem2.strip()
- for bline in blk_lines:
- belems = bline.split('\t')
- if type(pitem2) == list:
- if p_group_cols[0] == 6:
- belems[p_group_cols[1]] = int(belems[p_group_cols[1]])
- if belems[p_group_cols[1]] in pitem2 and belems[s_group_cols[1]] == sitem2:
- repeats2.append(int(belems[13]))
- repeats1.append(int(belems[6]))
- if belems[4] == 'mononucleotide':
- thresholds.append(mono_threshold)
- else:
- thresholds.append(non_mono_threshold)
- mut2["%s\t%s" % ( pitem2, sitem2 )] = mutabilityEstimator( repeats2, repeats1, thresholds )
- #count2[str(pitem2)+'\t'+str(sitem2)]=len(repeats2)
- if region == 'align':
- count2["%s\t%s" % ( pitem2, sitem2 )] = min( sum(repeats1), sum(repeats2) )
- else:
- if winspecies == 1:
- count2["%s\t%s" % ( pitem2, sitem2 )] = len(repeats2)
- elif winspecies == 2:
- count2["%s\t%s" % ( pitem2, sitem2 )] = len(repeats1)
- else:
- if type(sitem2) == list:
- if s_group_cols[0] == 6:
- belems[s_group_cols[1]] = int(belems[s_group_cols[1]])
- if belems[p_group_cols[1]] == pitem2 and belems[s_group_cols[1]] in sitem2:
- repeats2.append(int(belems[13]))
- repeats1.append(int(belems[6]))
- if belems[4] == 'mononucleotide':
- thresholds.append(mono_threshold)
- else:
- thresholds.append(non_mono_threshold)
- mut2["%s\t%s" % ( pitem2, sitem2 )] = mutabilityEstimator( repeats2, repeats1, thresholds )
- if region == 'align':
- count2["%s\t%s" % ( pitem2, sitem2 )] = min( sum(repeats1), sum(repeats2) )
- else:
- if winspecies == 1:
- count2["%s\t%s" % ( pitem2, sitem2 )] = len(repeats2)
- elif winspecies == 2:
- count2["%s\t%s" % ( pitem2, sitem2 )] = len(repeats1)
- else:
- if belems[p_group_cols[1]] == pitem2 and belems[s_group_cols[1]] == sitem2:
- repeats1.append(int(belems[13]))
- repeats2.append(int(belems[6]))
- if belems[4] == 'mononucleotide':
- thresholds.append(mono_threshold)
- else:
- thresholds.append(non_mono_threshold)
- mut2["%s\t%s" % ( pitem2, sitem2 )] = mutabilityEstimator( repeats2, repeats1, thresholds )
- if region == 'align':
- count2["%s\t%s" % ( pitem2, sitem2 )] = min( sum(repeats1), sum(repeats2) )
- else:
- if winspecies == 1:
- count2["%s\t%s" % ( pitem2, sitem2 )] = len(repeats2)
- elif winspecies == 2:
- count2["%s\t%s" % ( pitem2, sitem2 )] = len(repeats1)
- else: #Sub-group by feature is None
- for bline in blk_lines:
- belems = bline.split('\t')
- if type(pitem2) == list:
- if p_group_cols[0] == 6:
- belems[p_group_cols[1]] = int(belems[p_group_cols[1]])
- if belems[p_group_cols[1]] in pitem2:
- repeats2.append(int(belems[13]))
- repeats1.append(int(belems[6]))
- if belems[4] == 'mononucleotide':
- thresholds.append(mono_threshold)
- else:
- thresholds.append(non_mono_threshold)
- else:
- if belems[p_group_cols[1]] == pitem2:
- repeats2.append(int(belems[13]))
- repeats1.append(int(belems[6]))
- if belems[4] == 'mononucleotide':
- thresholds.append(mono_threshold)
- else:
- thresholds.append(non_mono_threshold)
- mut2["%s" % (pitem2)] = mutabilityEstimator( repeats2, repeats1, thresholds )
- if region == 'align':
- count2["%s" % (pitem2)] = min( sum(repeats1), sum(repeats2) )
- else:
- if winspecies == 1:
- count2["%s" % (pitem2)] = sum(repeats2)
- elif winspecies == 2:
- count2["%s" % (pitem2)] = sum(repeats1)
- for key in mut1.keys():
- if key in mut2.keys():
- mut = (mut1[key][0]+mut2[key][0])/(mut1[key][1]+mut2[key][1])
- count = count1[key]
- del mut2[key]
- else:
- unit_found = False
- if p_group_cols[0] == 7 or s_group_cols[0] == 7: #if it is Repeat Unit (AG, GCT etc.) check for reverse-complements too
- if p_group_cols[0] == 7:
- this, other = 0, 1
- else:
- this, other = 1, 0
- groups1 = key.split('\t')
- mutn = mut1[key][0]
- mutd = mut1[key][1]
- count = 0
- for key2 in mut2.keys():
- groups2 = key2.split('\t')
- if groups1[other] == groups2[other]:
- if groups1[this] in groups2[this]*2 or reverse_complement(groups1[this]) in groups2[this]*2:
- #mut = (mut1[key][0]+mut2[key2][0])/(mut1[key][1]+mut2[key2][1])
- mutn += mut2[key2][0]
- mutd += mut2[key2][1]
- count += int(count2[key2])
- unit_found = True
- del mut2[key2]
- #break
- if unit_found:
- mut = mutn/mutd
- else:
- mut = mut1[key][0]/mut1[key][1]
- count = count1[key]
- mut = "%.2e" % (mut/num_generations)
- if region == 'align':
- print >> fout, str(blk) + '\t'+seq1 + '\t' + seq2 + '\t' +key.strip()+ '\t'+str(mut) + '\t'+ str(count)
- elif region == 'win':
- fout.write("%s\t%s\t%s\t%s\n" % ( blk, key.strip(), mut, count ))
- fout.flush()
-
- #catch any remaining repeats, for instance if the orthologous position contained different repeat units
- for remaining_key in mut2.keys():
- mut = mut2[remaining_key][0]/mut2[remaining_key][1]
- mut = "%.2e" % (mut/num_generations)
- count = count2[remaining_key]
- if region == 'align':
- print >> fout, str(blk) + '\t'+seq1 + '\t'+seq2 + '\t'+remaining_key.strip()+ '\t'+str(mut)+ '\t'+ str(count)
- elif region == 'win':
- fout.write("%s\t%s\t%s\t%s\n" % ( blk, remaining_key.strip(), mut, count ))
- fout.flush()
- #print >> fout, blk + '\t'+remaining_key.strip()+ '\t'+str(mut)+ '\t'+ str(count)
-
-
-def counter(node, start, end, report_func):
- if start <= node.start < end and start < node.end <= end:
- report_func(node)
- if node.right:
- counter(node.right, start, end, report_func)
- if node.left:
- counter(node.left, start, end, report_func)
- elif node.start < start and node.right:
- counter(node.right, start, end, report_func)
- elif node.start >= end and node.left and node.left.maxend > start:
- counter(node.left, start, end, report_func)
-
-
-def main():
- infile = sys.argv[1]
-
- for i, line in enumerate( file ( infile )):
- line = line.rstrip('\r\n')
- if len( line )>0 and not line.startswith( '#' ):
- elems = line.split( '\t' )
- break
- if i == 30:
- break # Hopefully we'll never get here...
-
- if len( elems ) != 15:
- stop_err( "This tool only works on tabular data output by 'Extract Orthologous Microsatellites from pair-wise alignments' tool. The data in your input dataset is either missing or not formatted properly." )
- global winspecies, speciesind
- if region == 'win':
- if dbkey_i in elems[1]:
- winspecies = 1
- speciesind = 1
- elif dbkey_i in elems[8]:
- winspecies = 2
- speciesind = 8
- else:
- stop_err("The species build corresponding to your interval file is not present in the Microsatellite file.")
-
- fin = open(infile, 'r')
- skipped = 0
- linestr = ""
-
- if region == 'win':
- msats = NiceReaderWrapper( fileinput.FileInput( infile ),
- chrom_col = speciesind,
- start_col = speciesind+1,
- end_col = speciesind+2,
- strand_col = -1,
- fix_strand = True)
- msatTree = quicksect.IntervalTree()
- for item in msats:
- if type( item ) is GenomicInterval:
- msatTree.insert( item, msats.linenum, item.fields )
-
- for iline in fint:
- try:
- iline = iline.rstrip('\r\n')
- if not(iline) or iline == "":
- continue
- ielems = iline.strip("\r\n").split('\t')
- ichr = ielems[chr_col_i]
- istart = int(ielems[start_col_i])
- iend = int(ielems[end_col_i])
- isrc = "%s.%s" % ( dbkey_i, ichr )
- if isrc not in msatTree.chroms:
- continue
- result = []
- root = msatTree.chroms[isrc] #root node for the chrom
- counter(root, istart, iend, lambda node: result.append( node ))
- if not(result):
- continue
- tmpfile1 = tempfile.NamedTemporaryFile('wb+')
- for node in result:
- tmpfile1.write("%s\n" % "\t".join( node.other ))
-
- tmpfile1.seek(0)
- output_writer(iline, tmpfile1.readlines())
- except:
- skipped += 1
- if skipped:
- print "Skipped %d intervals as invalid." % (skipped)
- elif region == 'align':
- if s_group_cols[0] != -1:
- print >> fout, "#Window\tSpecies_1\tSpecies_2\tGroupby_Feature\tSubGroupby_Feature\tMutability\tCount"
- else:
- print >> fout, "#Window\tSpecies_1\tWindow_Start\tWindow_End\tSpecies_2\tGroupby_Feature\tMutability\tCount"
- prev_bnum = -1
- try:
- for line in fin:
- line = line.strip("\r\n")
- if not(line) or line == "":
- continue
- elems = line.split('\t')
- try:
- assert int(elems[0])
- assert len(elems) == 15
- except:
- continue
- new_bnum = int(elems[0])
- if new_bnum != prev_bnum:
- if prev_bnum != -1:
- output_writer(prev_bnum, linestr.strip().replace('\r','\n').split('\n'))
- linestr = line + "\n"
- else:
- linestr += line
- linestr += "\n"
- prev_bnum = new_bnum
- output_writer(prev_bnum, linestr.strip().replace('\r','\n').split('\n'))
- except Exception, ea:
- print >> sys.stderr, ea
- skipped += 1
- if skipped:
- print "Skipped %d lines as invalid." % (skipped)
-
-
-if __name__ == "__main__":
- main()
diff --git a/tools/regVariation/microsats_mutability.xml b/tools/regVariation/microsats_mutability.xml
deleted file mode 100644
index ef65084e794..00000000000
--- a/tools/regVariation/microsats_mutability.xml
+++ /dev/null
@@ -1,121 +0,0 @@
-
- by specified attributes
-
- microsats_mutability.py
- $input1
- $out_file1
- ${pri_condition.primary_group}
- #if $pri_condition.primary_group == "6":
- ${pri_condition.binsize} ${pri_condition.subgroup} -1
- #else:
- 0 ${pri_condition.sub_condition.subgroup}
- #if $pri_condition.sub_condition.subgroup == "6":
- ${pri_condition.sub_condition.s_binsize}
- #else:
- -1
- #end if
- #end if
- $gens
- ${region.type}
- #if $region.type == "win":
- ${region.input2} $input2.dbkey $input2.metadata.chromCol,$input2.metadata.startCol,$input2.metadata.endCol,$input2.metadata.strandCol
- #else:
- "None"
- #end if
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-.. class:: infomark
-
-**What it does**
-
-This tool computes microsatellite mutability for the orthologous microsatellites fetched from 'Extract Orthologous Microsatellites from pair-wise alignments' tool.
-
-Mutability is computed according to the method described in the following paper:
-
-*Webster et al., Microsatellite evolution inferred from human-chimpanzee genomic sequence alignments, Proc Natl Acad Sci 2002 June 25; 99(13): 8748-8753*
-
------
-
-.. class:: warningmark
-
-**Note**
-
-The user selected group and subgroup by features, the computed mutability and the count of the number of repeats used to compute that mutability are added as columns to the output.
-
-
diff --git a/tools/regVariation/partialR_square.py b/tools/regVariation/partialR_square.py
deleted file mode 100755
index b4fbcd6dff1..00000000000
--- a/tools/regVariation/partialR_square.py
+++ /dev/null
@@ -1,147 +0,0 @@
-#!/usr/bin/env python
-
-from galaxy import eggs
-
-import sys
-from rpy import *
-import numpy
-
-#export PYTHONPATH=~/galaxy/lib/
-#running command python partialR_square.py reg_inp.tab 4 1,2,3 partialR_result.tabular
-
-def stop_err(msg):
- sys.stderr.write(msg)
- sys.exit()
-
-
-def sscombs(s):
- if len(s) == 1:
- return [s]
- else:
- ssc = sscombs(s[1:])
- return [s[0]] + [s[0]+comb for comb in ssc] + ssc
-
-
-infile = sys.argv[1]
-y_col = int(sys.argv[2])-1
-x_cols = sys.argv[3].split(',')
-outfile = sys.argv[4]
-
-print "Predictor columns: %s; Response column: %d" % ( x_cols, y_col+1 )
-fout = open(outfile,'w')
-
-for i, line in enumerate( file ( infile )):
- line = line.rstrip('\r\n')
- if len( line )>0 and not line.startswith( '#' ):
- elems = line.split( '\t' )
- break
- if i == 30:
- break # Hopefully we'll never get here...
-
-if len( elems )<1:
- stop_err( "The data in your input dataset is either missing or not formatted properly." )
-
-y_vals = []
-x_vals = []
-
-for k, col in enumerate(x_cols):
- x_cols[k] = int(col)-1
- x_vals.append([])
- """
- try:
- float( elems[x_cols[k]] )
- except:
- try:
- msg = "This operation cannot be performed on non-numeric column %d containing value '%s'." % ( col, elems[x_cols[k]] )
- except:
- msg = "This operation cannot be performed on non-numeric data."
- stop_err( msg )
- """
-NA = 'NA'
-for ind, line in enumerate( file( infile )):
- if line and not line.startswith( '#' ):
- try:
- fields = line.split("\t")
- try:
- yval = float(fields[y_col])
- except Exception, ey:
- yval = r('NA')
- #print >> sys.stderr, "ey = %s" %ey
- y_vals.append(yval)
- for k, col in enumerate(x_cols):
- try:
- xval = float(fields[col])
- except Exception, ex:
- xval = r('NA')
- #print >> sys.stderr, "ex = %s" %ex
- x_vals[k].append(xval)
- except:
- pass
-
-x_vals1 = numpy.asarray(x_vals).transpose()
-dat = r.list(x=array(x_vals1), y=y_vals)
-
-set_default_mode(NO_CONVERSION)
-try:
- full = r.lm(r("y ~ x"), data= r.na_exclude(dat)) #full model includes all the predictor variables specified by the user
-except RException, rex:
- stop_err("Error performing linear regression on the input data.\nEither the response column or one of the predictor columns contain no numeric values.")
-set_default_mode(BASIC_CONVERSION)
-
-summary = r.summary(full)
-fullr2 = summary.get('r.squared','NA')
-
-if fullr2 == 'NA':
- stop_err("Error in linear regression")
-
-if len(x_vals) < 10:
- s = ""
- for ch in range(len(x_vals)):
- s += str(ch)
-else:
- stop_err("This tool only works with less than 10 predictors.")
-
-print >> fout, "#Model\tR-sq\tpartial_R_Terms\tpartial_R_Value"
-all_combos = sorted(sscombs(s), key=len)
-all_combos.reverse()
-for j, cols in enumerate(all_combos):
- #if len(cols) == len(s): #Same as the full model above
- # continue
- if len(cols) == 1:
- x_vals1 = x_vals[int(cols)]
- else:
- x_v = []
- for col in cols:
- x_v.append(x_vals[int(col)])
- x_vals1 = numpy.asarray(x_v).transpose()
- dat = r.list(x=array(x_vals1), y=y_vals)
- set_default_mode(NO_CONVERSION)
- red = r.lm(r("y ~ x"), data= dat) #Reduced model
- set_default_mode(BASIC_CONVERSION)
- summary = r.summary(red)
- redr2 = summary.get('r.squared','NA')
- try:
- partial_R = (float(fullr2)-float(redr2))/(1-float(redr2))
- except:
- partial_R = 'NA'
- col_str = ""
- for col in cols:
- col_str = col_str + str(int(x_cols[int(col)]) + 1) + " "
- col_str.strip()
- partial_R_col_str = ""
- for col in s:
- if col not in cols:
- partial_R_col_str = partial_R_col_str + str(int(x_cols[int(col)]) + 1) + " "
- partial_R_col_str.strip()
- if len(cols) == len(s): #full model
- partial_R_col_str = "-"
- partial_R = "-"
- try:
- redr2 = "%.4f" % (float(redr2))
- except:
- pass
- try:
- partial_R = "%.4f" % (float(partial_R))
- except:
- pass
- print >> fout, "%s\t%s\t%s\t%s" % ( col_str, redr2, partial_R_col_str, partial_R )
diff --git a/tools/regVariation/partialR_square.xml b/tools/regVariation/partialR_square.xml
deleted file mode 100755
index 4068a07a374..00000000000
--- a/tools/regVariation/partialR_square.xml
+++ /dev/null
@@ -1,68 +0,0 @@
-
-
-
- partialR_square.py
- $input1
- $response_col
- $predictor_cols
- $out_file1
- 1>/dev/null
-
-
-
-
-
-
-
-
-
-
-
-
- rpy
-
-
-
-
-
-
-
-
-
-
-
-
-
-.. class:: infomark
-
-**TIP:** If your data is not TAB delimited, use *Edit Datasets->Convert characters*
-
------
-
-.. class:: infomark
-
-**What it does**
-
-This tool computes the Partial R squared for all possible variable subsets using the following formula:
-
-**Partial R squared = [SSE(without i: 1,2,...,p-1) - SSE (full: 1,2,..,i..,p-1) / SSE(without i: 1,2,...,p-1)]**, which denotes the case where the 'i'th predictor is dropped.
-
-
-
-In general, **Partial R squared = [SSE(without i: 1,2,...,p-1) - SSE (full: 1,2,..,i..,p-1) / SSE(without i: 1,2,...,p-1)]**, where,
-
-- SSE (full: 1,2,..,i..,p-1) = Sum of Squares left out by the full set of predictors SSE(X1, X2 … Xp)
-- SSE (full: 1,2,..,i..,p-1) = Sum of Squares left out by the set of predictors excluding; for example, if we omit the first predictor, it will be SSE(X2 … Xp).
-
-
-The 4 columns in the output are described below:
-
-- Column 1 (Model): denotes the variables present in the model
-- Column 2 (R-sq): denotes the R-squared value corresponding to the model in Column 1
-- Column 3 (Partial R squared_Terms): denotes the variable/s for which Partial R squared is computed. These are the variables that are absent in the reduced model in Column 1. A '-' in this column indicates that the model in Column 1 is the Full model.
-- Column 4 (Partial R squared): denotes the Partial R squared value corresponding to the variable/s in Column 3. A '-' in this column indicates that the model in Column 1 is the Full model.
-
-*R Development Core Team (2010). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. ISBN 3-900051-07-0, URL http://www.R-project.org.*
-
-
-
diff --git a/tools/regVariation/quality_filter.py b/tools/regVariation/quality_filter.py
deleted file mode 100644
index 12750060319..00000000000
--- a/tools/regVariation/quality_filter.py
+++ /dev/null
@@ -1,242 +0,0 @@
-#!/usr/bin/env python
-#Guruprasad Ananda
-"""
-Filter based on nucleotide quality (PHRED score).
-
-usage: %prog input out_file primary_species mask_species score mask_char mask_region mask_region_length
-"""
-
-
-from __future__ import division
-from galaxy import eggs
-import pkg_resources
-pkg_resources.require( "bx-python" )
-pkg_resources.require( "lrucache" )
-try:
- pkg_resources.require("numpy")
-except:
- pass
-
-import sys
-import os, os.path
-from UserDict import DictMixin
-from bx.binned_array import FileBinnedArray
-from bx.bitset import *
-from bx.bitset_builders import *
-from bx.cookbook import doc_optparse
-from galaxy.tools.exception_handling import *
-import bx.align.maf
-
-class FileBinnedArrayDir( DictMixin ):
- """
- Adapter that makes a directory of FileBinnedArray files look like
- a regular dict of BinnedArray objects.
- """
- def __init__( self, dir ):
- self.dir = dir
- self.cache = dict()
- def __getitem__( self, key ):
- value = None
- if key in self.cache:
- value = self.cache[key]
- else:
- fname = os.path.join( self.dir, "%s.qa.bqv" % key )
- if os.path.exists( fname ):
- value = FileBinnedArray( open( fname ) )
- self.cache[key] = value
- if value is None:
- raise KeyError( "File does not exist: " + fname )
- return value
-
-def stop_err(msg):
- sys.stderr.write(msg)
- sys.exit()
-
-def load_scores_ba_dir( dir ):
- """
- Return a dict-like object (keyed by chromosome) that returns
- FileBinnedArray objects created from "key.ba" files in `dir`
- """
- return FileBinnedArrayDir( dir )
-
-def bitwise_and ( string1, string2, maskch ):
- result = []
- for i, ch in enumerate(string1):
- try:
- ch = int(ch)
- except:
- pass
- if string2[i] == '-':
- ch = 1
- if ch and string2[i]:
- result.append(string2[i])
- else:
- result.append(maskch)
- return ''.join(result)
-
-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.cols )
- inp_file, out_file, pri_species, mask_species, qual_cutoff, mask_chr, mask_region, mask_length, loc_file = args
- qual_cutoff = int(qual_cutoff)
- mask_chr = int(mask_chr)
- mask_region = int(mask_region)
- if mask_region != 3:
- mask_length = int(mask_length)
- else:
- mask_length_r = int(mask_length.split(',')[0])
- mask_length_l = int(mask_length.split(',')[1])
- except:
- stop_err( "Data issue, click the pencil icon in the history item to correct the metadata attributes of the input dataset." )
-
- if pri_species == 'None':
- stop_err( "No primary species selected, try again by selecting at least one primary species." )
- if mask_species == 'None':
- stop_err( "No mask species selected, try again by selecting at least one species to mask." )
-
- mask_chr_count = 0
- mask_chr_dict = {0:'#', 1:'$', 2:'^', 3:'*', 4:'?', 5:'N'}
- mask_reg_dict = {0:'Current pos', 1:'Current+Downstream', 2:'Current+Upstream', 3:'Current+Both sides'}
-
- #ensure dbkey is present in the twobit loc file
- try:
- pspecies_all = pri_species.split(',')
- pspecies_all2 = pri_species.split(',')
- pspecies = []
- filepaths = []
- for line in open(loc_file):
- if pspecies_all2 == []:
- break
- if line[0:1] == "#":
- continue
- fields = line.split('\t')
- try:
- build = fields[0]
- for i, dbkey in enumerate(pspecies_all2):
- if dbkey == build:
- pspecies.append(build)
- filepaths.append(fields[1])
- del pspecies_all2[i]
- else:
- continue
- except:
- pass
- except Exception, exc:
- stop_err( 'Initialization errorL %s' % str( exc ) )
-
- if len(pspecies) == 0:
- stop_err( "Quality scores are not available for the following genome builds: %s" % ( pspecies_all2 ) )
- if len(pspecies) < len(pspecies_all):
- print "Quality scores are not available for the following genome builds: %s" % (pspecies_all2)
-
- scores_by_chrom = []
- #Get scores for all the primary species
- for file in filepaths:
- scores_by_chrom.append(load_scores_ba_dir( file.strip() ))
-
- try:
- maf_reader = bx.align.maf.Reader( open(inp_file, 'r') )
- maf_writer = bx.align.maf.Writer( open(out_file,'w') )
- except Exception, e:
- stop_err( "Your MAF file appears to be malformed: %s" % str( e ) )
-
- maf_count = 0
- for block in maf_reader:
- status_strings = []
- for seq in range (len(block.components)):
- src = block.components[seq].src
- dbkey = src.split('.')[0]
- chr = src.split('.')[1]
- if not (dbkey in pspecies):
- continue
- else: #enter if the species is a primary species
- index = pspecies.index(dbkey)
- sequence = block.components[seq].text
- s_start = block.components[seq].start
- size = len(sequence) #this includes the gaps too
- status_str = '1'*size
- status_list = list(status_str)
- if status_strings == []:
- status_strings.append(status_str)
- ind = 0
- s_end = block.components[seq].end
- #Get scores for the entire sequence
- try:
- scores = scores_by_chrom[index][chr][s_start:s_end]
- except:
- continue
- pos = 0
- while pos < (s_end-s_start):
- if sequence[ind] == '-': #No score for GAPS
- ind += 1
- continue
- score = scores[pos]
- if score < qual_cutoff:
- score = 0
-
- if not(score):
- if mask_region == 0: #Mask Corresponding position only
- status_list[ind] = '0'
- ind += 1
- pos += 1
- elif mask_region == 1: #Mask Corresponding position + downstream neighbors
- for n in range(mask_length+1):
- try:
- status_list[ind+n] = '0'
- except:
- pass
- ind = ind + mask_length + 1
- pos = pos + mask_length + 1
- elif mask_region == 2: #Mask Corresponding position + upstream neighbors
- for n in range(mask_length+1):
- try:
- status_list[ind-n] = '0'
- except:
- pass
- ind += 1
- pos += 1
- elif mask_region == 3: #Mask Corresponding position + neighbors on both sides
- for n in range(-mask_length_l, mask_length_r+1):
- try:
- status_list[ind+n] = '0'
- except:
- pass
- ind = ind + mask_length_r + 1
- pos = pos + mask_length_r + 1
- else:
- pos += 1
- ind += 1
-
- status_strings.append(''.join(status_list))
-
- if status_strings == []: #this block has no primary species
- continue
- output_status_str = status_strings[0]
- for stat in status_strings[1:]:
- try:
- output_status_str = bitwise_and (status_strings[0], stat, '0')
- except Exception, e:
- break
-
- for seq in range (len(block.components)):
- src = block.components[seq].src
- dbkey = src.split('.')[0]
- if dbkey not in mask_species.split(','):
- continue
- sequence = block.components[seq].text
- sequence = bitwise_and (output_status_str, sequence, mask_chr_dict[mask_chr])
- block.components[seq].text = sequence
- mask_chr_count += output_status_str.count('0')
- maf_writer.write(block)
- maf_count += 1
-
- maf_reader.close()
- maf_writer.close()
- print "No. of blocks = %d; No. of masked nucleotides = %s; Mask character = %s; Mask region = %s; Cutoff used = %d" % (maf_count, mask_chr_count, mask_chr_dict[mask_chr], mask_reg_dict[mask_region], qual_cutoff)
-
-
-if __name__ == "__main__":
- main()
diff --git a/tools/regVariation/quality_filter.xml b/tools/regVariation/quality_filter.xml
deleted file mode 100644
index 0ebedc448b4..00000000000
--- a/tools/regVariation/quality_filter.xml
+++ /dev/null
@@ -1,115 +0,0 @@
-
- based on quality scores
-
- quality_filter.py
- $input
- $out_file1
- $primary_species
- $mask_species
- $score
- $mask_char
- ${mask_region.region}
- #if $mask_region.region == "3"
- ${mask_region.lengthr},${mask_region.lengthl}
- #elif $mask_region.region == "0"
- 1
- #else
- ${mask_region.length}
- #end if
- ${GALAXY_DATA_INDEX_DIR}/quality_scores.loc
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
- numpy
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-.. class:: infomark
-
-**What it does**
-
-This tool takes a MAF file as input and filters nucleotides in every alignment block of the MAF file based on their quality/PHRED scores.
-
------
-
-.. class:: warningmark
-
-**Note**
-
-Any block/s not containing the primary species (species whose quality scores is to be used), will be omitted.
-Also, any primary species whose quality scores are not available in Galaxy will be considered as a non-primary species. This info will appear as a message in the job history panel.
-
------
-
-**Example**
-
-- For the following alignment block::
-
- a score=4050.0
- s hg18.chrX 3719221 48 - 154913754 tattttacatttaaaataaatatgtaaatatatattttatatttaaaa
- s panTro2.chrX 3560945 48 - 155361357 tattttatatttaaaataaagatgtaaatatatattttatatttaaaa
-
-- running this tool with **Primary species as panTro2**, **Mask species as hg18, panTro2**, **Quality cutoff as 20**, **Mask character as #** and **Mask region as only the corresponding position** will return::
-
- a score=4050.0
- s hg18.chrX 3719221 48 - 154913754 ###tttac#####a###a#atatgtaaat###tattt#####ttaaaa
- s panTro2.chrX 3560945 48 - 155361357 ###tttat#####a###a#agatgtaaat###tattt#####ttaaaa
-
- where, the positions containing # represent panTro2 nucleotides having quality scores less than 20.
-
-
diff --git a/tools/regVariation/qv_to_bqv.py b/tools/regVariation/qv_to_bqv.py
deleted file mode 100755
index 172dce0e4cb..00000000000
--- a/tools/regVariation/qv_to_bqv.py
+++ /dev/null
@@ -1,89 +0,0 @@
-#!/usr/bin/env python
-
-"""
-Adapted from bx/scripts/qv_to_bqv.py
-
-Convert a qual (qv) file to several BinnedArray files for fast seek.
-This script takes approximately 4 seconds per 1 million base pairs.
-
-The input format is fasta style quality -- fasta headers followed by
-whitespace separated integers.
-
-usage: %prog qual_file output_file
-"""
-
-import pkg_resources
-pkg_resources.require( "bx-python" )
-pkg_resources.require( "numpy" )
-import os
-import sys
-import tempfile
-from bx.binned_array import BinnedArrayWriter
-from bx.cookbook import *
-import fileinput
-
-def load_scores_ba_dir( dir ):
- """
- Return a dict-like object (keyed by chromosome) that returns
- FileBinnedArray objects created from "key.ba" files in `dir`
- """
- return FileBinnedArrayDir( dir )
-
-def main():
- args = sys.argv[1:]
- try:
- qual_file_dir = args[0]
- mydir = "/home/gua110/Desktop/rhesus_quality_scores/rheMac2.qual.qv"
- qual_file_dir = mydir.replace(mydir.split("/")[-1], "")
- output_file = args[1]
- fo = open(output_file, "w")
- except:
- print "usage: qual_file output_file"
- sys.exit()
-
- tmpfile = tempfile.NamedTemporaryFile()
- cmdline = "ls " + qual_file_dir + "*.qa | cat >> " + tmpfile.name
- os.system (cmdline)
- for qual_file in tmpfile.readlines():
- qual = fileinput.FileInput( qual_file.strip() )
- outfile = None
- outbin = None
- base_count = 0
- mega_count = 0
-
- for line in qual:
- line = line.rstrip("\r\n")
- if line.startswith(">"):
- # close old
- if outbin and outfile:
- print "\nFinished region " + region + " at " + str(base_count) + " base pairs."
- outbin.finish()
- outfile.close()
- # start new file
- region = line.lstrip(">")
- #outfname = output_file + "." + region + ".bqv" #CHANGED
- outfname = qual_file.strip() + ".bqv"
- print >> fo, "Writing region " + region + " to file " + outfname
- outfile = open( outfname , "wb")
- outbin = BinnedArrayWriter(outfile, typecode='b', default=0)
- base_count = 0
- mega_count = 0
- else:
- if outfile and outbin:
- nums = line.split()
- for val in nums:
- outval = int(val)
- assert outval <= 255 and outval >= 0
- outbin.write(outval)
- base_count += 1
- if (mega_count * 1000000) <= base_count:
- sys.stdout.write(str(mega_count)+" ")
- sys.stdout.flush()
- mega_count = base_count // 1000000 + 1
- if outbin and outfile:
- print "\nFinished region " + region + " at " + str(base_count) + " base pairs."
- outbin.finish()
- outfile.close()
-
-if __name__ == "__main__":
- main()
diff --git a/tools/regVariation/qv_to_bqv.xml b/tools/regVariation/qv_to_bqv.xml
deleted file mode 100644
index d899a8d5fc2..00000000000
--- a/tools/regVariation/qv_to_bqv.xml
+++ /dev/null
@@ -1,17 +0,0 @@
-
-
- qv_to_bqv.py "$input1" $output
-
-
-
-
-
-
-
-
-
-
-
-
-
-
\ No newline at end of file
diff --git a/tools/regVariation/rcve.py b/tools/regVariation/rcve.py
deleted file mode 100644
index 2e7113165a4..00000000000
--- a/tools/regVariation/rcve.py
+++ /dev/null
@@ -1,144 +0,0 @@
-#!/usr/bin/env python
-
-from galaxy import eggs
-
-import sys
-from rpy import *
-import numpy
-
-def stop_err(msg):
- sys.stderr.write(msg)
- sys.exit()
-
-
-def sscombs(s):
- if len(s) == 1:
- return [s]
- else:
- ssc = sscombs(s[1:])
- return [s[0]] + [s[0]+comb for comb in ssc] + ssc
-
-
-infile = sys.argv[1]
-y_col = int(sys.argv[2])-1
-x_cols = sys.argv[3].split(',')
-outfile = sys.argv[4]
-
-print "Predictor columns: %s; Response column: %d" % ( x_cols, y_col+1 )
-fout = open(outfile,'w')
-
-for i, line in enumerate( file ( infile )):
- line = line.rstrip('\r\n')
- if len( line )>0 and not line.startswith( '#' ):
- elems = line.split( '\t' )
- break
- if i == 30:
- break # Hopefully we'll never get here...
-
-if len( elems )<1:
- stop_err( "The data in your input dataset is either missing or not formatted properly." )
-
-y_vals = []
-x_vals = []
-
-for k, col in enumerate(x_cols):
- x_cols[k] = int(col)-1
- x_vals.append([])
- """
- try:
- float( elems[x_cols[k]] )
- except:
- try:
- msg = "This operation cannot be performed on non-numeric column %d containing value '%s'." % ( col, elems[x_cols[k]] )
- except:
- msg = "This operation cannot be performed on non-numeric data."
- stop_err( msg )
- """
-NA = 'NA'
-for ind, line in enumerate( file( infile )):
- if line and not line.startswith( '#' ):
- try:
- fields = line.split("\t")
- try:
- yval = float(fields[y_col])
- except Exception, ey:
- yval = r('NA')
- #print >>sys.stderr, "ey = %s" %ey
- y_vals.append(yval)
- for k, col in enumerate(x_cols):
- try:
- xval = float(fields[col])
- except Exception, ex:
- xval = r('NA')
- #print >>sys.stderr, "ex = %s" %ex
- x_vals[k].append(xval)
- except:
- pass
-
-x_vals1 = numpy.asarray(x_vals).transpose()
-dat = r.list( x=array(x_vals1), y=y_vals )
-
-set_default_mode(NO_CONVERSION)
-try:
- full = r.lm( r("y ~ x"), data=r.na_exclude(dat) ) #full model includes all the predictor variables specified by the user
-except RException, rex:
- stop_err("Error performing linear regression on the input data.\nEither the response column or one of the predictor columns contain no numeric values.")
-set_default_mode(BASIC_CONVERSION)
-
-summary = r.summary(full)
-fullr2 = summary.get('r.squared','NA')
-
-if fullr2 == 'NA':
- stop_err("Error in linear regression")
-
-if len(x_vals) < 10:
- s = ""
- for ch in range(len(x_vals)):
- s += str(ch)
-else:
- stop_err("This tool only works with less than 10 predictors.")
-
-print >> fout, "#Model\tR-sq\tRCVE_Terms\tRCVE_Value"
-all_combos = sorted(sscombs(s), key=len)
-all_combos.reverse()
-for j, cols in enumerate(all_combos):
- #if len(cols) == len(s): #Same as the full model above
- # continue
- if len(cols) == 1:
- x_vals1 = x_vals[int(cols)]
- else:
- x_v = []
- for col in cols:
- x_v.append(x_vals[int(col)])
- x_vals1 = numpy.asarray(x_v).transpose()
- dat = r.list(x=array(x_vals1), y=y_vals)
- set_default_mode(NO_CONVERSION)
- red = r.lm(r("y ~ x"), data= dat) #Reduced model
- set_default_mode(BASIC_CONVERSION)
- summary = r.summary(red)
- redr2 = summary.get('r.squared','NA')
- try:
- rcve = (float(fullr2)-float(redr2))/float(fullr2)
- except:
- rcve = 'NA'
- col_str = ""
- for col in cols:
- col_str = col_str + str(int(x_cols[int(col)]) + 1) + " "
- col_str.strip()
- rcve_col_str = ""
- for col in s:
- if col not in cols:
- rcve_col_str = rcve_col_str + str(int(x_cols[int(col)]) + 1) + " "
- rcve_col_str.strip()
- if len(cols) == len(s): #full model
- rcve_col_str = "-"
- rcve = "-"
- try:
- redr2 = "%.4f" % (float(redr2))
- except:
- pass
- try:
- rcve = "%.4f" % (float(rcve))
- except:
- pass
- print >> fout, "%s\t%s\t%s\t%s" % ( col_str, redr2, rcve_col_str, rcve )
diff --git a/tools/regVariation/rcve.xml b/tools/regVariation/rcve.xml
deleted file mode 100644
index 6807d9ac549..00000000000
--- a/tools/regVariation/rcve.xml
+++ /dev/null
@@ -1,70 +0,0 @@
-
-
-
- rcve.py
- $input1
- $response_col
- $predictor_cols
- $out_file1
- 1>/dev/null
-
-
-
-
-
-
-
-
-
-
-
-
- rpy
-
-
-
-
-
-
-
-
-
-
-
-
-
-.. class:: infomark
-
-**TIP:** If your data is not TAB delimited, use *Edit Datasets->Convert characters*
-
------
-
-.. class:: infomark
-
-**What it does**
-
-This tool computes the RCVE (Relative Contribution to Variance) for all possible variable subsets using the following formula:
-
-**RCVE(i) = [R-sq (full: 1,2,..,i..,p-1) - R-sq(without i: 1,2,...,p-1)] / R-sq (full: 1,2,..,i..,p-1)**,
-which denotes the case where the 'i'th predictor is dropped.
-
-
-In general,
-**RCVE(X+) = [R-sq (full: {X,X+}) - R-sq(reduced: {X})] / R-sq (full: {X,X+})**,
-where,
-
-- {X,X+} denotes the set of all predictors,
-- X+ is the set of predictors for which we compute RCVE (and therefore drop from the full model to obtain a reduced one),
-- {X} is the set of the predictors that are left in the reduced model after excluding {X+}
-
-
-The 4 columns in the output are described below:
-
-- Column 1 (Model): denotes the variables present in the model ({X})
-- Column 2 (R-sq): denotes the R-squared value corresponding to the model in Column 1
-- Column 3 (RCVE_Terms): denotes the variable/s for which RCVE is computed ({X+}). These are the variables that are absent in the reduced model in Column 1. A '-' in this column indicates that the model in Column 1 is the Full model.
-- Column 4 (RCVE): denotes the RCVE value corresponding to the variable/s in Column 3. A '-' in this column indicates that the model in Column 1 is the Full model.
-
-
-
-
diff --git a/tools/regVariation/substitution_rates.py b/tools/regVariation/substitution_rates.py
deleted file mode 100644
index 2fe03d7210c..00000000000
--- a/tools/regVariation/substitution_rates.py
+++ /dev/null
@@ -1,123 +0,0 @@
-#!/usr/bin/env python
-#guruprasad Ananda
-"""
-Estimates substitution rates from pairwise alignments using JC69 model.
-"""
-
-from galaxy import eggs
-from galaxy.tools.util.galaxyops import *
-from galaxy.tools.util import maf_utilities
-import bx.align.maf
-import fileinput
-import sys
-
-def stop_err(msg):
- sys.stderr.write(msg)
- sys.exit()
-
-
-if len(sys.argv) < 3:
- stop_err("Incorrect number of arguments.")
-
-inp_file = sys.argv[1]
-out_file = sys.argv[2]
-fout = open(out_file, 'w')
-int_file = sys.argv[3]
-if int_file != "None": #The user has specified an interval file
- dbkey_i = sys.argv[4]
- chr_col_i, start_col_i, end_col_i, strand_col_i = parse_cols_arg( sys.argv[5] )
-
-
-def rateEstimator(block):
- global alignlen, mismatches
-
- src1 = block.components[0].src
- sequence1 = block.components[0].text
- start1 = block.components[0].start
- end1 = block.components[0].end
- len1 = int(end1)-int(start1)
- len1_withgap = len(sequence1)
- mismatch = 0.0
-
- for seq in range (1, len(block.components)):
- src2 = block.components[seq].src
- sequence2 = block.components[seq].text
- start2 = block.components[seq].start
- end2 = block.components[seq].end
- len2 = int(end2)-int(start2)
- for nt in range(len1_withgap):
- if sequence1[nt] not in '-#$^*?' and sequence2[nt] not in '-#$^*?': # Not a gap or masked character
- if sequence1[nt].upper() != sequence2[nt].upper():
- mismatch += 1
-
- if int_file == "None":
- p = mismatch/min(len1, len2)
- print >> fout, "%s\t%s\t%s\t%s\t%s\t%s\t%d\t%d\t%.4f" % ( src1, start1, end1, src2, start2, end2, min(len1, len2), mismatch, p )
- else:
- mismatches += mismatch
- alignlen += min(len1, len2)
-
-
-def main():
- skipped = 0
- not_pairwise = 0
-
- if int_file == "None":
- try:
- maf_reader = bx.align.maf.Reader( open(inp_file, 'r') )
- except:
- stop_err("Your MAF file appears to be malformed.")
- print >> fout, "#Seq1\tStart1\tEnd1\tSeq2\tStart2\tEnd2\tL\tN\tp"
- for block in maf_reader:
- if len(block.components) != 2:
- not_pairwise += 1
- continue
- try:
- rateEstimator(block)
- except:
- skipped += 1
- else:
- index, index_filename = maf_utilities.build_maf_index( inp_file, species = [dbkey_i] )
- if index is None:
- print >> sys.stderr, "Your MAF file appears to be malformed."
- sys.exit()
- win = NiceReaderWrapper( fileinput.FileInput( int_file ),
- chrom_col=chr_col_i,
- start_col=start_col_i,
- end_col=end_col_i,
- strand_col=strand_col_i,
- fix_strand=True)
- species = None
- mincols = 0
- global alignlen, mismatches
-
- for interval in win:
- alignlen = 0
- mismatches = 0.0
- src = "%s.%s" % ( dbkey_i, interval.chrom )
- for block in maf_utilities.get_chopped_blocks_for_region( index, src, interval, species, mincols ):
- if len(block.components) != 2:
- not_pairwise += 1
- continue
- try:
- rateEstimator(block)
- except:
- skipped += 1
- if alignlen:
- p = mismatches/alignlen
- else:
- p = 'NA'
- interval.fields.append(str(alignlen))
- interval.fields.append(str(mismatches))
- interval.fields.append(str(p))
- print >> fout, "\t".join(interval.fields)
- #num_blocks += 1
-
- if not_pairwise:
- print "Skipped %d non-pairwise blocks" % (not_pairwise)
- if skipped:
- print "Skipped %d blocks as invalid" % (skipped)
-
-
-if __name__ == "__main__":
- main()
diff --git a/tools/regVariation/substitution_rates.xml b/tools/regVariation/substitution_rates.xml
deleted file mode 100644
index a4201a675d9..00000000000
--- a/tools/regVariation/substitution_rates.xml
+++ /dev/null
@@ -1,61 +0,0 @@
-
- for non-coding regions
-
- substitution_rates.py
- $input
- $out_file1
- #if $region.type == "win":
- ${region.input2} ${region.input2.dbkey} ${region.input2.metadata.chromCol},$region.input2.metadata.startCol,$region.input2.metadata.endCol,$region.input2.metadata.strandCol
- #else:
- "None"
- #end if
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-.. class:: infomark
-
-**What it does**
-
-This tool takes a pairwise MAF file as input and estimates substitution rate according to Jukes-Cantor JC69 model. The 3 new columns appended to the output are explained below:
-
-- L: number of nucleotides compared
-- N: number of different nucleotides
-- p = N/L
-
------
-
-.. class:: warningmark
-
-**Note**
-
-Any block/s not containing exactly two sequences, will be omitted.
-
-
-
\ No newline at end of file
diff --git a/tools/regVariation/substitutions.py b/tools/regVariation/substitutions.py
deleted file mode 100644
index 7774ca5e274..00000000000
--- a/tools/regVariation/substitutions.py
+++ /dev/null
@@ -1,85 +0,0 @@
-#!/usr/bin/env python
-#Guruprasad ANanda
-"""
-Fetches substitutions from pairwise alignments.
-"""
-
-from galaxy import eggs
-
-from galaxy.tools.util import maf_utilities
-
-import bx.align.maf
-import sys
-
-def stop_err(msg):
- sys.stderr.write(msg)
- sys.exit()
-
-
-if len(sys.argv) < 3:
- stop_err("Incorrect number of arguments.")
-
-inp_file = sys.argv[1]
-out_file = sys.argv[2]
-fout = open(out_file, 'w')
-
-def fetchSubs(block):
- src1 = block.components[0].src
- sequence1 = block.components[0].text
- start1 = block.components[0].start
- end1 = block.components[0].end
- len1_withgap = len(sequence1)
-
- for seq in range (1, len(block.components)):
- src2 = block.components[seq].src
- sequence2 = block.components[seq].text
- start2 = block.components[seq].start
- end2 = block.components[seq].end
- sub_begin = None
- sub_end = None
- begin = False
-
- for nt in range(len1_withgap):
- if sequence1[nt] not in '-#$^*?' and sequence2[nt] not in '-#$^*?': # Not a gap or masked character
- if sequence1[nt].upper() != sequence2[nt].upper():
- if not(begin):
- sub_begin = nt
- begin = True
- sub_end = nt
- else:
- if begin:
- print >> fout, "%s\t%s\t%s" % ( src1, start1+sub_begin-sequence1[0:sub_begin].count('-'), start1+sub_end-sequence1[0:sub_end].count('-') )
- print >> fout, "%s\t%s\t%s" % ( src2, start2+sub_begin-sequence2[0:sub_begin].count('-'), start2+sub_end-sequence2[0:sub_end].count('-') )
- begin = False
- else:
- if begin:
- print >> fout, "%s\t%s\t%s" % ( src1, start1+sub_begin-sequence1[0:sub_begin].count('-'), end1+sub_end-sequence1[0:sub_end].count('-') )
- print >> fout, "%s\t%s\t%s" % ( src2, start2+sub_begin-sequence2[0:sub_begin].count('-'), end2+sub_end-sequence2[0:sub_end].count('-') )
- begin = False
-
-
-def main():
- skipped = 0
- not_pairwise = 0
- try:
- maf_reader = bx.align.maf.Reader( open(inp_file, 'r') )
- except:
- stop_err("Your MAF file appears to be malformed.")
- print >> fout, "#Chr\tStart\tEnd"
- for block in maf_reader:
- if len(block.components) != 2:
- not_pairwise += 1
- continue
- try:
- fetchSubs(block)
- except:
- skipped += 1
-
- if not_pairwise:
- print "Skipped %d non-pairwise blocks" % (not_pairwise)
- if skipped:
- print "Skipped %d blocks" % (skipped)
-
-
-if __name__ == "__main__":
- main()
diff --git a/tools/regVariation/substitutions.xml b/tools/regVariation/substitutions.xml
deleted file mode 100644
index d982db8725d..00000000000
--- a/tools/regVariation/substitutions.xml
+++ /dev/null
@@ -1,38 +0,0 @@
-
- from pairwise alignments
-
- substitutions.py
- $input
- $out_file1
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-.. class:: infomark
-
-**What it does**
-
-This tool takes a pairwise MAF file as input and fetches substitutions per alignment block.
-
------
-
-.. class:: warningmark
-
-**Note**
-
-Any block/s not containing exactly two sequences, will be omitted.
-
-
-
\ No newline at end of file
diff --git a/tools/regVariation/windowSplitter.py b/tools/regVariation/windowSplitter.py
deleted file mode 100644
index f1d9896d8ae..00000000000
--- a/tools/regVariation/windowSplitter.py
+++ /dev/null
@@ -1,84 +0,0 @@
-#!/usr/bin/env python
-
-"""
-Split into windows.
-
-usage: %prog input size out_file
- -l, --cols=N,N,N,N: Columns for chrom, start, end, strand in file
-"""
-
-import sys
-
-from galaxy import eggs
-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.cols )
- inp_file, winsize, out_file, makesliding, offset = args
- winsize = int(winsize)
- offset = int(offset)
- makesliding = int(makesliding)
- except:
- stop_err( "Data issue, click the pencil icon in the history item to correct the metadata attributes of the input dataset." )
-
- fo = open(out_file,'w')
-
- skipped_lines = 0
- first_invalid_line = 0
- invalid_line = None
- if offset == 0:
- makesliding = 0
-
- for i, line in enumerate( file( inp_file ) ):
- line = line.strip()
- if line and line[0:1] != "#":
- try:
- elems = line.split('\t')
- start = int(elems[start_col_1])
- end = int(elems[end_col_1])
- if makesliding == 0:
- numwin = (end - start)/winsize
- else:
- numwin = (end - start)/offset
- if numwin > 0:
- for win in range(numwin):
- elems_1 = elems
- elems_1[start_col_1] = str(start)
- elems_1[end_col_1] = str(start + winsize)
- fo.write( "%s\n" % '\t'.join( elems_1 ) )
- if makesliding == 0:
- start = start + winsize
- else:
- start = start + offset
- if start+winsize > end:
- break
- except:
- skipped_lines += 1
- if not invalid_line:
- first_invalid_line = i + 1
- invalid_line = line
-
- fo.close()
-
- if makesliding == 1:
- print 'Window size=%d, Sliding=Yes, Offset=%d' % ( winsize, offset )
- else:
- print 'Window size=%d, Sliding=No' % (winsize)
- if skipped_lines > 0:
- print 'Skipped %d invalid lines starting with #%d: "%s"' % ( skipped_lines, first_invalid_line, invalid_line )
-
-
-if __name__ == "__main__":
- main()
diff --git a/tools/regVariation/windowSplitter.xml b/tools/regVariation/windowSplitter.xml
deleted file mode 100644
index 4caaf96edc6..00000000000
--- a/tools/regVariation/windowSplitter.xml
+++ /dev/null
@@ -1,104 +0,0 @@
-
-
- windowSplitter.py $input $size $out_file1 ${wintype.choice} ${wintype.offset} -l ${input.metadata.chromCol},${input.metadata.startCol},${input.metadata.endCol},${input.metadata.strandCol}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-.. class:: infomark
-
-**What it does**
-
-This tool splits the intervals in the input file into smaller intervals based on the specified window-size and window type.
-
------
-
-.. class:: warningmark
-
-**Note**
-
-The positions at the end of the input interval which do not fit into the last window or a new window of required size, will be omitted from the output.
-
------
-
-.. class:: infomark
-
-**About formats**
-
-**BED format** Browser Extensible Data format was designed at UCSC for displaying data tracks in the Genome Browser. It has three required fields and several additional optional ones:
-
-The first three BED fields (required) are::
-
- 1. chrom - The name of the chromosome (e.g. chr1, chrY_random).
- 2. chromStart - The starting position in the chromosome. (The first base in a chromosome is numbered 0.)
- 3. chromEnd - The ending position in the chromosome, plus 1 (i.e., a half-open interval).
-
-The additional BED fields (optional) are::
-
- 4. name - The name of the BED line.
- 5. score - A score between 0 and 1000.
- 6. strand - Defines the strand - either '+' or '-'.
- 7. thickStart - The starting position where the feature is drawn thickly at the Genome Browser.
- 8. thickEnd - The ending position where the feature is drawn thickly at the Genome Browser.
- 9. reserved - This should always be set to zero.
- 10. blockCount - The number of blocks (exons) in the BED line.
- 11. blockSizes - A comma-separated list of the block sizes. The number of items in this list should correspond to blockCount.
- 12. blockStarts - A comma-separated list of block starts. All of the blockStart positions should be calculated relative to chromStart. The number of items in this list should correspond to blockCount.
- 13. expCount - The number of experiments.
- 14. expIds - A comma-separated list of experiment ids. The number of items in this list should correspond to expCount.
- 15. expScores - A comma-separated list of experiment scores. All of the expScores should be relative to expIds. The number of items in this list should correspond to expCount.
-
------
-
-**Example**
-
-- For the following dataset::
-
- chr22 1000 4700 NM_174568 0 +
-
-- running this tool with **Window size as 1000**, will return::
-
- chr22 1000 2000 NM_174568 0 +
- chr22 2000 3000 NM_174568 0 +
- chr22 3000 4000 NM_174568 0 +
-
-- running this tool to make **Sliding windows** of **size 1000** and **offset 500**, will return::
-
- chr22 1000 2000 NM_174568 0 +
- chr22 1500 2500 NM_174568 0 +
- chr22 2000 3000 NM_174568 0 +
- chr22 2500 3500 NM_174568 0 +
- chr22 3000 4000 NM_174568 0 +
- chr22 3500 4500 NM_174568 0 +
-
-
-
-
-
\ No newline at end of file