From 0177cd59cc1ddeb00a31edfdea6cf1bda8a3d867 Mon Sep 17 00:00:00 2001 From: Daniel Blankenberg Date: Thu, 8 Feb 2007 21:24:39 +0000 Subject: [PATCH] Add tool to allow users to limit a MAF to a set of species. --- tool_conf.xml.sample | 1 + tools/filters/maf/maf_limit_to_species.py | 47 +++++++++++++++++++ tools/filters/maf/maf_limit_to_species.xml | 43 +++++++++++++++++ .../filters/maf/maf_limit_to_species_code.py | 47 +++++++++++++++++++ 4 files changed, 138 insertions(+) create mode 100644 tools/filters/maf/maf_limit_to_species.py create mode 100644 tools/filters/maf/maf_limit_to_species.xml create mode 100644 tools/filters/maf/maf_limit_to_species_code.py diff --git a/tool_conf.xml.sample b/tool_conf.xml.sample index fa0b78799b9..23ca3586de9 100644 --- a/tool_conf.xml.sample +++ b/tool_conf.xml.sample @@ -62,6 +62,7 @@ +
diff --git a/tools/filters/maf/maf_limit_to_species.py b/tools/filters/maf/maf_limit_to_species.py new file mode 100644 index 00000000000..339c5e2a00a --- /dev/null +++ b/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() diff --git a/tools/filters/maf/maf_limit_to_species.xml b/tools/filters/maf/maf_limit_to_species.xml new file mode 100644 index 00000000000..66b124c0089 --- /dev/null +++ b/tools/filters/maf/maf_limit_to_species.xml @@ -0,0 +1,43 @@ + + Remove undesired species from a MAF file + maf_limit_to_species.py $species $input1 $out_file1 $allow_partial $min_species + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +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. + + + + diff --git a/tools/filters/maf/maf_limit_to_species_code.py b/tools/filters/maf/maf_limit_to_species_code.py new file mode 100644 index 00000000000..993138ede48 --- /dev/null +++ b/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 [("You must wait for the MAF file to be created before you can limit species.",'None',True)]