Enable 'extract genomic DNA' tool to accept and produce GFF files and added functional tests for this feature.

lib/galaxy/tools/util/gff_util.py

@@ -6,23 +6,37 @@ from bx.intervals.io import NiceReaderWrapper, GenomicInterval
 
 class GFFReaderWrapper( NiceReaderWrapper ):
     """
-    Reader wrapper converts GFF format--starting and ending coordinates are 1-based, closed--to the 'traditional' interval format--0 based, 
-    half-open. This is useful when using GFF files as inputs to tools that expect traditional interval format.
+    Reader wrapper converts GFF format--starting and ending coordinates are 1-based, closed--to the 
+    'traditional'/BED interval format--0 based, half-open. This is useful when using GFF files as inputs 
+    to tools that expect traditional interval format.
     """
     def parse_row( self, line ):
-        interval = GenomicInterval( self, line.split( "\t" ), self.chrom_col, self.start_col, self.end_col, self.strand_col, self.default_strand, fix_strand=self.fix_strand )
-        # Change from 1-based to 0-based format.
-        interval.start -= 1
-        # Add 1 to end to move from closed to open format for end coordinate.
-        interval.end += 1
+        interval = GenomicInterval( self, line.split( "\t" ), self.chrom_col, self.start_col, self.end_col, \
+                                    self.strand_col, self.default_strand, fix_strand=self.fix_strand )
+        interval = convert_gff_coords_to_bed( interval )
         return interval
         
-def convert_to_gff_coordinates( interval ):
+def convert_bed_coords_to_gff( interval ):
     """
-    Converts a GenomicInterval's coordinates to GFF format.
+    Converts an interval object's coordinates from BED format to GFF format. Accepted object types include 
+    GenomicInterval and list (where the first element in the list is the interval's start, and the second 
+    element is the interval's end).
     """
     if type( interval ) is GenomicInterval:
         interval.start += 1
-        interval.end -= 1
-        return interval
+    elif type ( interval ) is list:
+        interval[ 0 ] += 1
     return interval
+    
+def convert_gff_coords_to_bed( interval ):
+    """
+    Converts an interval object's coordinates from GFF format to BED format. Accepted object types include 
+    GenomicInterval and list (where the first element in the list is the interval's start, and the second 
+    element is the interval's end).
+    """
+    if type( interval ) is GenomicInterval:
+        interval.start -= 1
+    elif type ( interval ) is list:
+        interval[ 0 ] -= 1
+    return interval
+    

tools/extract/extract_genomic_dna.py

@@ -5,6 +5,7 @@ usage: %prog $input $out_file1
     -d, --dbkey=N: Genome build of input file
     -o, --output_format=N: the data type of the output file
     -g, --GALAXY_DATA_INDEX_DIR=N: the directory containing alignseq.loc
+    -G, --gff: input and output file, when it is interval, coordinates are treated as GFF format (1-based, half-open) rather than 'traditional' 0-based, closed format.
 """
 from galaxy import eggs
 import pkg_resources
@@ -14,6 +15,7 @@ from bx.cookbook import doc_optparse
 import bx.seq.nib
 import bx.seq.twobit
 from galaxy.tools.util.galaxyops import *
+from galaxy.tools.util.gff_util import *
 
 assert sys.version_info[:2] >= ( 2, 4 )
     
@@ -50,6 +52,7 @@ def __main__():
         chrom_col, start_col, end_col, strand_col = parse_cols_arg( options.cols )
         dbkey = options.dbkey
         output_format = options.output_format
+        gff_format = options.gff
         GALAXY_DATA_INDEX_DIR = options.GALAXY_DATA_INDEX_DIR
         input_filename, output_filename = args
     except:
@@ -80,6 +83,8 @@ def __main__():
                 chrom = fields[chrom_col]
                 start = int( fields[start_col] )
                 end = int( fields[end_col] )
+                if gff_format:
+                    start, end = convert_gff_coords_to_bed( [start, end] )
                 if includes_strand_col:
                     strand = fields[strand_col]
             except:
@@ -162,7 +167,11 @@ def __main__():
                     c = b
             else: # output_format == "interval"
                 meta_data = "\t".join( fields )
-                fout.write( "%s\t%s\n" % ( meta_data, str( sequence ) ) )
+                if gff_format:
+                    format_str = "%s seq \"%s\";\n"
+                else:
+                    format_str = "%s\t%s\n"
+                fout.write( format_str % ( meta_data, str( sequence ) ) )
 
     fout.close()
 

tools/extract/extract_genomic_dna.xml

@@ -1,20 +1,27 @@
 <tool id="Extract genomic DNA 1" name="Extract Genomic DNA" version="2.2.1">
   <description>using coordinates from assembled/unassembled genomes</description>
-  <command interpreter="python">extract_genomic_dna.py $input $out_file1 -1 ${input.metadata.chromCol},${input.metadata.startCol},${input.metadata.endCol},${input.metadata.strandCol} -d $dbkey -o $out_format -g ${GALAXY_DATA_INDEX_DIR}</command>
+  <command interpreter="python">
+      extract_genomic_dna.py $input $out_file1 -d $dbkey -o $out_format -g ${GALAXY_DATA_INDEX_DIR}
+      #if isinstance( $input.datatype, $__app__.datatypes_registry.get_datatype_by_extension('gff').__class__):
+        -1 1,4,5,7 --gff
+      #else:
+        -1 ${input.metadata.chromCol},${input.metadata.startCol},${input.metadata.endCol},${input.metadata.strandCol}
+      #end if
+  </command>
   <inputs>
-    <param format="interval" name="input" type="data" label="Fetch sequences corresponding to Query">
-      <validator type="unspecified_build" />
-      <validator type="dataset_metadata_in_file" filename="alignseq.loc" metadata_name="dbkey" metadata_column="1" message="Sequences are not currently available for the specified build." line_startswith="seq" />
+    <param format="interval,gff" name="input" type="data" label="Fetch sequences corresponding to Query">
+        <validator type="unspecified_build" />
+        <validator type="dataset_metadata_in_file" filename="alignseq.loc" metadata_name="dbkey" metadata_column="1" message="Sequences are not currently available for the specified build." line_startswith="seq" />
     </param>
 	<param name="out_format" type="select" label="Output data type">
-	<option value="fasta">FASTA</option>
-	<option value="interval">Interval</option>
+    	<option value="fasta">FASTA</option>
+    	<option value="interval">Interval</option>
 	</param>
   </inputs>
   <outputs>
-    <data format="fasta" name="out_file1" metadata_source="input">
+    <data format="input" name="out_file1" metadata_source="input">
       <change_format>
-        <when input="out_format" value="interval" format="interval" />
+        <when input="out_format" value="fasta" format="fasta" />
       </change_format>
     </data>
   </outputs>
@@ -34,6 +41,17 @@
       <param name="out_format" value="interval"/>
       <output name="out_file1" file="extract_genomic_dna_out3.interval" />
     </test>
+    <!-- Test GFF file support. -->
+    <test>
+      <param name="input" value="gff_filtering_out1.gff" dbkey="mm9" ftype="gff" />
+      <param name="out_format" value="interval"/>
+      <output name="out_file1" file="extract_genomic_dna_out4.gff" />
+    </test>
+    <test>
+      <param name="input" value="gff_filtering_out1.gff" dbkey="mm9" ftype="gff" />
+      <param name="out_format" value="fasta"/>
+      <output name="out_file1" file="extract_genomic_dna_out5.fasta" />
+    </test>
   </tests>
   <help>
 
@@ -90,7 +108,7 @@ Extracting sequences with **FASTA** output data type returns::
     CACCAAAACCCTCATCAAGACAATTGTCACCAGGATCAATGACATTTCAC
     ACACG
 
-Extrracting sequences with **Interval** output data type returns::
+Extracting sequences with **Interval** output data type returns::
 
     chr7    127475281       127475310       NM_000230       0       +       GTAGGAATCGCAGCGCCAGCGGTTGCAAG
     chr7    127485994       127486166       NM_000230       0       +       GCCCAAGAAGCCCATCCTGGGAAGGAAAATGCATTGGGGAACCCTGTGCGGATTCTTGTGGCTTTGGCCCTATCTTTTCTATGTCCAAGCTGTGCCCATCCAAAAAGTCCAAGATGACACCAAAACCCTCATCAAGACAATTGTCACCAGGATCAATGACATTTCACACACG

tools/new_operations/gops_intersect.py

@@ -70,7 +70,7 @@ def main():
         for line in intersect( [g1,g2], pieces=pieces, mincols=mincols ):
             if type( line ) == GenomicInterval:
                 if in1_gff_format:
-                    line = convert_to_gff_coordinates( line )
+                    line = convert_bed_coords_to_gff( line )
                 out_file.write( "%s\n" % "\t".join( line.fields ) )
             else:
                 out_file.write( "%s\n" % line )

tools/new_operations/gops_subtract.py

@@ -71,7 +71,7 @@ def main():
         for line in subtract( [g1,g2], pieces=pieces, mincols=mincols ):
             if type( line ) is GenomicInterval:
                 if in1_gff_format:
-                    line = convert_to_gff_coordinates( line )
+                    line = convert_bed_coords_to_gff( line )
                 out_file.write( "%s\n" % "\t".join( line.fields ) )
             else:
                 out_file.write( "%s\n" % line )