Add tool to allow users to limit a MAF to a set of species.

tool_conf.xml.sample

@@ -62,6 +62,7 @@
     <tool file="extract/extractAxt_wrapper.xml" />
     <tool file="extract/interval2maf.xml" />
     <tool file="extract/user_interval2maf.xml" />
+    <tool file="filters/maf/maf_limit_to_species.xml"/>
   </section>
   <section name="Alignment Viewers" id="align_view">
     <tool file="visualization/GMAJ.xml" />

tools/filters/maf/maf_limit_to_species.py

@@ -0,0 +1,47 @@
+#!/usr/bin/env python2.3
+
+"""
+Read a maf file and write out a new maf with only blocks having the 
+required species, after dropping any other species and removing
+columns containing only gaps.
+
+usage: %prog species,species2,... input_maf output_maf allow_partial min_species_per_block
+"""
+
+import pkg_resources; pkg_resources.require( "bx-python" )
+import bx.align.maf
+import sys
+
+
+def main():
+
+    species = sys.argv[1].split( ',' )
+
+    maf_reader = bx.align.maf.Reader( open(sys.argv[2],'r') )
+    maf_writer = bx.align.maf.Writer( open(sys.argv[3],'w') )
+    allow_partial = False
+    if int(sys.argv[4]): allow_partial = True
+    min_species_per_block = int(sys.argv[5])
+    
+    maf_blocks_kept = 0
+
+    for m in maf_reader:        
+        new_components = []    
+        for comp in m.components:
+            if comp.src.split( '.' )[0] in species:
+                new_components.append( comp )
+        m.components = new_components
+        m.remove_all_gap_columns()
+        output_block = False
+        if allow_partial: output_block = True
+        elif len( m.components ) == len(species): output_block = True
+        if output_block and len( m.components ) > min_species_per_block:
+            maf_writer.write( m )
+            maf_blocks_kept += 1
+        
+    maf_reader.close()
+    maf_writer.close()
+    print "%i MAF blocks have been kept." % maf_blocks_kept
+
+if __name__ == "__main__": 
+    main()

tools/filters/maf/maf_limit_to_species.xml

@@ -0,0 +1,43 @@
+<tool id="MAF_Limit_To_Species1" name="Maf Limit to Species">
+    <description>Remove undesired species from a MAF file</description>
+    <command interpreter="python2.4">maf_limit_to_species.py $species $input1 $out_file1 $allow_partial $min_species</command>
+    <inputs>
+        <page>
+            <param format="maf" name="input1" type="data" label="MAF file"/>
+            <param name="allow_partial" type="select">
+              <label>Exclude blocks which have missing species</label>
+              <option value="1">No</option>
+              <option value="0">Yes</option>
+            </param>
+            <param name="min_species" type="select">
+              <label>Exclude blocks which have only one species</label>
+              <option value="1">Yes</option>
+              <option value="0">No</option>
+            </param>
+        </page>
+        <page>
+            <param name="species" type="select" label="Species to keep" dynamic_options="get_available_species( input1.file_name )" display="checkboxes" multiple="true"/>
+        </page>
+    </inputs>
+    <outputs>
+        <data format="maf" name="out_file1" />
+    </outputs>
+    <tests>
+      <test>
+        <param name="input1" value="4.maf"/>
+        <param name="species" values="bosTau2,canFam2,hg17,panTro1,rheMac2,rn3"/>
+        <param name="allow_partial" value="0"/>
+        <param name="min_species" value="0"/>
+        <output name="out_file1" file="cf_maf_limit_to_species.dat"/>
+      </test>
+    </tests>
+    <help>
+This tool allows the user to remove any undesired species from a MAF file. Columns which contain only gaps are removed.
+
+The user can choose to keep or discard blocks which do not have all the desired species. 
+
+The user can also choose to include blocks which contain only a single species.
+    </help>
+    <code file="maf_limit_to_species_code.py"/>
+</tool>
+

tools/filters/maf/maf_limit_to_species_code.py

@@ -0,0 +1,47 @@
+import pkg_resources; pkg_resources.require( "bx-python" )
+from bx.align import maf
+# No initialization required.
+
+#return lists of species available, showing gapped and ungapped base counts
+def get_available_species( input_filename ):
+    try:
+        rval = []
+        species={}
+        
+        file_in = open(input_filename, 'r')
+        maf_reader = maf.Reader( file_in )
+        
+        for i, m in enumerate( maf_reader ):
+            l = m.components
+            for c in l:
+                spec,chrom = maf.src_split( c.src )
+                if not spec or not chrom:
+                    spec = chrom = c.src
+                species[spec]=spec
+        
+        file_in.close()
+        
+        species = species.keys()
+        species.sort()
+
+        file_in = open(input_filename, 'r')
+        maf_reader = maf.Reader( file_in )
+        
+        species_sequence={}
+        for s in species: species_sequence[s] = []
+        
+        for m in maf_reader:
+            for s in species:
+                c = m.get_component_by_src_start( s ) 
+                if c: species_sequence[s].append( c.text )
+                else: species_sequence[s].append( "-" * m.text_size )
+        
+        file_in.close()
+        
+        for spec in species:
+            species_sequence[spec] = "".join(species_sequence[spec])
+            rval.append( (spec + ": "+str(len(species_sequence[spec]) - species_sequence[spec].count("-"))+" nongap, "+ str(len(species_sequence[spec])) + " total bases",spec,False) )
+                
+        return rval
+    except:
+        return [("<B>You must wait for the MAF file to be created before you can limit species.</B>",'None',True)]