diff --git a/lib/galaxy/datatypes/converters/maf_to_fasta_converter.py b/lib/galaxy/datatypes/converters/maf_to_fasta_converter.py index a8c8067730a..a6401d90cf0 100644 --- a/lib/galaxy/datatypes/converters/maf_to_fasta_converter.py +++ b/lib/galaxy/datatypes/converters/maf_to_fasta_converter.py @@ -14,12 +14,15 @@ def __main__(): input_name = sys.argv.pop(1) out = open( output_name, 'w' ) count = 0 - for count, maf in enumerate( bx.align.maf.Reader( open( input_name, 'r' ) ) ): - for c in maf.components: - spec, chrom = bx.align.maf.src_split( c.src ) - if not spec or not chrom: - spec = chrom = c.src - out.write( "%s\n" % maf_utilities.get_fasta_header( c, suffix = "%s_%i" % ( spec, count ) ) ) + for count, block in enumerate( bx.align.maf.Reader( open( input_name, 'r' ) ) ): + spec_counts = {} + for c in block.components: + spec, chrom = maf_utilities.src_split( c.src ) + if spec not in spec_counts: + spec_counts[ spec ] = 0 + else: + spec_counts[ spec ] += 1 + out.write( "%s\n" % maf_utilities.get_fasta_header( c, { 'block_index' : count, 'species' : spec, 'sequence_index' : spec_counts[ spec ] }, suffix = "%s_%i_%i" % ( spec, count, spec_counts[ spec ] ) ) ) out.write( "%s\n" % c.text ) out.write( "\n" ) out.close() @@ -27,3 +30,12 @@ def __main__(): if __name__ == "__main__": __main__() + + for component in block.components: + spec, chrom = maf_utilities.src_split( component.src ) + if spec not in spec_counts: + spec_counts[ spec ] = 0 + else: + spec_counts[ spec ] += 1 + file_out.write( "%s\n" % maf_utilities.get_fasta_header( component, { 'block_index' : block_num, 'species' : spec, 'sequence_index' : spec_counts[ spec ] }, suffix = "%s_%i_%i" % ( spec, block_num, spec_counts[ spec ] ) ) ) + file_out.write( "%s\n" % component.text ) diff --git a/lib/galaxy/datatypes/converters/maf_to_fasta_converter.xml b/lib/galaxy/datatypes/converters/maf_to_fasta_converter.xml index d4798724d90..643a94ef75c 100644 --- a/lib/galaxy/datatypes/converters/maf_to_fasta_converter.xml +++ b/lib/galaxy/datatypes/converters/maf_to_fasta_converter.xml @@ -1,4 +1,4 @@ - + maf_to_fasta_converter.py $output1 $input1 diff --git a/lib/galaxy/datatypes/converters/maf_to_interval_converter.py b/lib/galaxy/datatypes/converters/maf_to_interval_converter.py index bde1c741e61..375b82a72c9 100644 --- a/lib/galaxy/datatypes/converters/maf_to_interval_converter.py +++ b/lib/galaxy/datatypes/converters/maf_to_interval_converter.py @@ -4,7 +4,8 @@ import sys from galaxy import eggs import pkg_resources; pkg_resources.require( "bx-python" ) -import bx.align.maf +import bx.align.maf +from galaxy.tools.util import maf_utilities assert sys.version_info[:2] >= ( 2, 4 ) @@ -17,15 +18,15 @@ def __main__(): #write interval header line out.write( "#chrom\tstart\tend\tstrand\n" ) try: - for maf in bx.align.maf.Reader( open(input_name, 'r') ): - c = maf.get_component_by_src_start(species) - if c is not None: - out.write( "%s\t%i\t%i\t%s\n" % (bx.align.src_split(c.src)[-1], c.get_forward_strand_start(), c.get_forward_strand_end(), c.strand) ) - count += 1 + for block in bx.align.maf.Reader( open( input_name, 'r' ) ): + for c in maf_utilities.iter_components_by_src_start( block, species ): + if c is not None: + out.write( "%s\t%i\t%i\t%s\n" % ( bx.align.src_split( c.src )[-1], c.get_forward_strand_start(), c.get_forward_strand_end(), c.strand ) ) + count += 1 except Exception, e: print >> sys.stderr, "There was a problem processing your input: %s" % e out.close() - print "%i MAF blocks converted to Genomic Intervals for species %s." % (count, species) + print "%i MAF blocks converted to Genomic Intervals for species %s." % ( count, species ) if __name__ == "__main__": __main__() diff --git a/lib/galaxy/datatypes/converters/maf_to_interval_converter.xml b/lib/galaxy/datatypes/converters/maf_to_interval_converter.xml index d644e146ee0..64ba291f047 100644 --- a/lib/galaxy/datatypes/converters/maf_to_interval_converter.xml +++ b/lib/galaxy/datatypes/converters/maf_to_interval_converter.xml @@ -1,4 +1,4 @@ - + maf_to_interval_converter.py $output1 $input1 ${input1.metadata.dbkey} diff --git a/lib/galaxy/datatypes/sequence.py b/lib/galaxy/datatypes/sequence.py index b68031d8cf5..d9cb1078555 100644 --- a/lib/galaxy/datatypes/sequence.py +++ b/lib/galaxy/datatypes/sequence.py @@ -22,7 +22,7 @@ class Sequence( data.Text ): pass class Alignment( Sequence ): - """Class describing an alignmnet""" + """Class describing an alignment""" """Add metadata elements""" MetadataElement( name="species", desc="Species", default=[], param=metadata.SelectParameter, multiple=True, readonly=True, no_value=None ) @@ -316,6 +316,78 @@ try: import bx.align.maf except: pass +#trying to import maf_utilities here throws an ImportError due to a circular import between jobs and tools: +#from galaxy.tools.util.maf_utilities import build_maf_index_species_chromosomes +#Traceback (most recent call last): +# File "./scripts/paster.py", line 27, in +# command.run() +# File "build/bdist.solaris-2.11-i86pc/egg/paste/script/command.py", line 78, in run +# File "build/bdist.solaris-2.11-i86pc/egg/paste/script/command.py", line 117, in invoke +# File "build/bdist.solaris-2.11-i86pc/egg/paste/script/command.py", line 212, in run +# File "build/bdist.solaris-2.11-i86pc/egg/paste/script/serve.py", line 227, in command +# File "build/bdist.solaris-2.11-i86pc/egg/paste/script/serve.py", line 250, in loadapp +# File "build/bdist.solaris-2.11-i86pc/egg/paste/deploy/loadwsgi.py", line 193, in loadapp +# File "build/bdist.solaris-2.11-i86pc/egg/paste/deploy/loadwsgi.py", line 213, in loadobj +# File "build/bdist.solaris-2.11-i86pc/egg/paste/deploy/loadwsgi.py", line 237, in loadcontext +# File "build/bdist.solaris-2.11-i86pc/egg/paste/deploy/loadwsgi.py", line 267, in _loadconfig +# File "build/bdist.solaris-2.11-i86pc/egg/paste/deploy/loadwsgi.py", line 397, in get_context +# File "build/bdist.solaris-2.11-i86pc/egg/paste/deploy/loadwsgi.py", line 439, in _context_from_explicit +# File "build/bdist.solaris-2.11-i86pc/egg/paste/deploy/loadwsgi.py", line 18, in import_string +# File "/afs/bx.psu.edu/home/dan/galaxy/central/lib/pkg_resources.py", line 1912, in load +# entry = __import__(self.module_name, globals(),globals(), ['__name__']) +# File "/afs/bx.psu.edu/home/dan/galaxy/central/lib/galaxy/web/buildapp.py", line 18, in +# from galaxy import config, jobs, util, tools +# File "/afs/bx.psu.edu/home/dan/galaxy/central/lib/galaxy/jobs/__init__.py", line 3, in +# from galaxy import util, model +# File "/afs/bx.psu.edu/home/dan/galaxy/central/lib/galaxy/model/__init__.py", line 13, in +# import galaxy.datatypes.registry +# File "/afs/bx.psu.edu/home/dan/galaxy/central/lib/galaxy/datatypes/registry.py", line 6, in +# import data, tabular, interval, images, sequence, qualityscore, genetics, xml, coverage, tracks, chrominfo +# File "/afs/bx.psu.edu/home/dan/galaxy/central/lib/galaxy/datatypes/sequence.py", line 344, in +# from galaxy.tools.util.maf_utilities import build_maf_index_species_chromosomes +# File "/afs/bx.psu.edu/home/dan/galaxy/central/lib/galaxy/tools/__init__.py", line 15, in +# from galaxy import util, jobs, model +#ImportError: cannot import name jobs +#so we'll copy and paste for now...terribly icky +#*** ANYCHANGE TO THIS METHOD HERE OR IN maf_utilities MUST BE PROPAGATED *** +def COPIED_build_maf_index_species_chromosomes( filename, index_species = None ): + species = [] + species_chromosomes = {} + indexes = bx.interval_index_file.Indexes() + try: + maf_reader = bx.align.maf.Reader( open( filename ) ) + while True: + pos = maf_reader.file.tell() + block = maf_reader.next() + if block is None: break + for c in block.components: + spec = c.src + chrom = None + if "." in spec: + spec, chrom = spec.split( ".", 1 ) + if spec not in species: + species.append( spec ) + species_chromosomes[spec] = [] + if chrom and chrom not in species_chromosomes[spec]: + species_chromosomes[spec].append( chrom ) + if index_species is None or spec in index_species: + forward_strand_start = c.forward_strand_start + forward_strand_end = c.forward_strand_end + try: + forward_strand_start = int( forward_strand_start ) + forward_strand_end = int( forward_strand_end ) + except ValueError: + continue #start and end are not integers, can't add component to index, goto next component + #this likely only occurs when parse_e_rows is True? + #could a species exist as only e rows? should the + if forward_strand_end > forward_strand_start: + #require positive length; i.e. certain lines have start = end = 0 and cannot be indexed + indexes.add( c.src, forward_strand_start, forward_strand_end, pos, max=c.src_size ) + except Exception, e: + #most likely a bad MAF + log.debug( 'Building MAF index on %s failed: %s' % ( filename, e ) ) + return ( None, [], {} ) + return ( indexes, species, species_chromosomes ) class Maf( Alignment ): """Class describing a Maf alignment""" @@ -333,38 +405,8 @@ class Maf( Alignment ): Parses and sets species, chromosomes, index from MAF file. """ #these metadata values are not accessable by users, always overwrite + indexes, species, species_chromosomes = COPIED_build_maf_index_species_chromosomes( dataset.file_name ) - try: - maf_reader = bx.align.maf.Reader( open( dataset.file_name ) ) - except: - return #not a maf file - species = [] - species_chromosomes = {} - indexes = bx.interval_index_file.Indexes() - while True: - pos = maf_reader.file.tell() - block = maf_reader.next() - if block is None: break - for c in block.components: - spec = c.src - chrom = None - if "." in spec: - spec, chrom = spec.split( ".", 1 ) - if spec not in species: - species.append(spec) - species_chromosomes[spec] = [] - if chrom and chrom not in species_chromosomes[spec]: - species_chromosomes[spec].append( chrom ) - forward_strand_start = c.forward_strand_start - forward_strand_end = c.forward_strand_end - try: - forward_strand_start = int( forward_strand_start ) - forward_strand_end = int( forward_strand_end ) - except ValueError: - continue #start and end are not integers, can't add component to index, goto next component - if forward_strand_end > forward_strand_start: - #require positive length; i.e. certain lines have start = end = 0 and cannot be indexed - indexes.add( c.src, forward_strand_start, forward_strand_end, pos, max=c.src_size ) dataset.metadata.species = species #only overwrite the contents if our newly determined chromosomes don't match stored chrom_file = dataset.metadata.species_chromosomes diff --git a/lib/galaxy/tools/util/maf_utilities.py b/lib/galaxy/tools/util/maf_utilities.py index 9243c031306..09ef49d0f8b 100644 --- a/lib/galaxy/tools/util/maf_utilities.py +++ b/lib/galaxy/tools/util/maf_utilities.py @@ -7,10 +7,41 @@ import pkg_resources; pkg_resources.require( "bx-python" ) import bx.align.maf import bx.intervals import bx.interval_index_file -import sys, os, string, tempfile +import sys, os, string, tempfile +import logging +from copy import deepcopy assert sys.version_info[:2] >= ( 2, 4 ) - + +log = logging.getLogger(__name__) + + +GAP_CHARS = [ '-' ] +SRC_SPLIT_CHAR = '.' + +def src_split( src ): + spec, chrom = bx.align.maf.src_split( src ) + if None in [ spec, chrom ]: + spec = chrom = src + return spec, chrom + +def src_merge( spec, chrom, contig = None ): + if None in [ spec, chrom ]: + spec = chrom = spec or chrom + return bx.align.maf.src_merge( spec, chrom, contig ) + +def get_species_in_block( block ): + species = [] + for c in block.components: + spec, chrom = src_split( c.src ) + if spec not in species: + species.append( spec ) + return species + +def tool_fail( msg = "Unknown Error" ): + print >> sys.stderr, "Fatal Error: %s" % msg + sys.exit() + #an object corresponding to a reference layered alignment class RegionAlignment( object ): @@ -153,69 +184,187 @@ def open_or_build_maf_index( maf_file, index_filename, species = None ): except: return build_maf_index( maf_file, species = species ) +#*** ANYCHANGE TO THIS METHOD HERE OR IN galaxy.datatypes.sequences MUST BE PROPAGATED *** +def build_maf_index_species_chromosomes( filename, index_species = None ): + species = [] + species_chromosomes = {} + indexes = bx.interval_index_file.Indexes() + try: + maf_reader = bx.align.maf.Reader( open( filename ) ) + while True: + pos = maf_reader.file.tell() + block = maf_reader.next() + if block is None: break + for c in block.components: + spec = c.src + chrom = None + if "." in spec: + spec, chrom = spec.split( ".", 1 ) + if spec not in species: + species.append( spec ) + species_chromosomes[spec] = [] + if chrom and chrom not in species_chromosomes[spec]: + species_chromosomes[spec].append( chrom ) + if index_species is None or spec in index_species: + forward_strand_start = c.forward_strand_start + forward_strand_end = c.forward_strand_end + try: + forward_strand_start = int( forward_strand_start ) + forward_strand_end = int( forward_strand_end ) + except ValueError: + continue #start and end are not integers, can't add component to index, goto next component + #this likely only occurs when parse_e_rows is True? + #could a species exist as only e rows? should the + if forward_strand_end > forward_strand_start: + #require positive length; i.e. certain lines have start = end = 0 and cannot be indexed + indexes.add( c.src, forward_strand_start, forward_strand_end, pos, max=c.src_size ) + except Exception, e: + #most likely a bad MAF + log.debug( 'Building MAF index on %s failed: %s' % ( filename, e ) ) + return ( None, [], {} ) + return ( indexes, species, species_chromosomes ) #builds and returns ( index, index_filename ) for specified maf_file def build_maf_index( maf_file, species = None ): - indexes = bx.interval_index_file.Indexes() - try: - maf_reader = bx.align.maf.Reader( open( maf_file ) ) - # Need to be a bit tricky in our iteration here to get the 'tells' right - while True: - pos = maf_reader.file.tell() - block = maf_reader.next() - if block is None: break - for c in block.components: - if species is not None and c.src.split( "." )[0] not in species: - continue - indexes.add( c.src, c.forward_strand_start, c.forward_strand_end, pos ) + indexes, found_species, species_chromosomes = build_maf_index_species_chromosomes( maf_file, species ) + if indexes is not None: fd, index_filename = tempfile.mkstemp() out = os.fdopen( fd, 'w' ) indexes.write( out ) out.close() - return ( bx.align.maf.Indexed( maf_file, index_filename = index_filename, keep_open = True, parse_e_rows = False ), index_filename ) - except: - return ( None, None ) - -def chop_block_by_region( block, src, region, species = None, mincols = 0, force_strand = None ): - ref = block.get_component_by_src( src ) - #We want our block coordinates to be from positive strand - if ref.strand == "-": - block = block.reverse_complement() - ref = block.get_component_by_src( src ) + return ( bx.align.maf.Indexed( maf_file, index_filename = index_filename, keep_open = True, parse_e_rows = False ), index_filename ) + return ( None, None ) + +def component_overlaps_region( c, region ): + if c is None: return False + start, end = c.get_forward_strand_start(), c.get_forward_strand_end() + if region.start >= end or region.end <= start: + return False + return True + +def chop_block_by_region( block, src, region, species = None, mincols = 0 ): + # This chopping method was designed to maintain consistency with how start/end padding gaps have been working in Galaxy thus far: + # behavior as seen when forcing blocks to be '+' relative to src sequence (ref) and using block.slice_by_component( ref, slice_start, slice_end ) + # whether-or-not this is the 'correct' behavior is questionable, but this will at least maintain consistency + # comments welcome + slice_start = block.text_size #max for the min() + slice_end = 0 #min for the max() + old_score = block.score #save old score for later use + # We no longer assume only one occurance of src per block, so we need to check them all + for c in iter_components_by_src( block, src ): + if component_overlaps_region( c, region ): + if c.text is not None: + rev_strand = False + if c.strand == "-": + #We want our coord_to_col coordinates to be returned from positive stranded component + rev_strand = True + c = c.reverse_complement() + start = max( region.start, c.start ) + end = min( region.end, c.end ) + start = c.coord_to_col( start ) + end = c.coord_to_col( end ) + if rev_strand: + #need to orient slice coordinates to the original block direction + slice_len = end - start + end = len( c.text ) - start + start = end - slice_len + slice_start = min( start, slice_start ) + slice_end = max( end, slice_end ) + + if slice_start < slice_end: + block = block.slice( slice_start, slice_end ) + if block.text_size > mincols: + # restore old score, may not be accurate, but it is better than 0 for everything? + block.score = old_score + if species is not None: + block = block.limit_to_species( species ) + block.remove_all_gap_columns() + return block + return None - #save old score here for later use - old_score = block.score - slice_start = max( region.start, ref.start ) - slice_end = min( region.end, ref.end ) - - #slice block by reference species at determined limits - block = block.slice_by_component( ref, slice_start, slice_end ) - - if block.text_size > mincols: - if ( force_strand is None and region.strand != ref.strand ) or ( force_strand is not None and force_strand != ref.strand ): - block = block.reverse_complement() - # restore old score, may not be accurate, but it is better than 0 for everything - block.score = old_score - if species is not None: - block = block.limit_to_species( species ) - block.remove_all_gap_columns() - return block - return None +def orient_block_by_region( block, src, region, force_strand = None ): + #loop through components matching src, + #make sure each of these components overlap region + #cache strand for each of overlaping regions + #if force_strand / region.strand not in strand cache, reverse complement + ### we could have 2 sequences with same src, overlapping region, on different strands, this would cause no reverse_complementing + strands = [ c.strand for c in iter_components_by_src( block, src ) if component_overlaps_region( c, region ) ] + if strands and ( force_strand is None and region.strand not in strands ) or ( force_strand is not None and force_strand not in strands ): + block = block.reverse_complement() + return block + +def get_oriented_chopped_blocks_for_region( index, src, region, species = None, mincols = 0, force_strand = None ): + for block, idx, offset in get_oriented_chopped_blocks_with_index_offset_for_region( index, src, region, species, mincols, force_strand ): + yield block +def get_oriented_chopped_blocks_with_index_offset_for_region( index, src, region, species = None, mincols = 0, force_strand = None ): + for block, idx, offset in get_chopped_blocks_with_index_offset_for_region( index, src, region, species, mincols ): + yield orient_block_by_region( block, src, region, force_strand ), idx, offset + +#split a block with multiple occurances of src into one block per src +def iter_blocks_split_by_src( block, src ): + for src_c in iter_components_by_src( block, src ): + new_block = bx.align.Alignment( score=block.score, attributes=deepcopy( block.attributes ) ) + new_block.text_size = block.text_size + for c in block.components: + if c == src_c or c.src != src: + new_block.add_component( deepcopy( c ) ) #components have reference to alignment, dont want to loose reference to original alignment block in original components + yield new_block + +#split a block into multiple blocks with all combinations of a species appearing only once per block +def iter_blocks_split_by_species( block, species = None ): + def __split_components_by_species( components_by_species, new_block ): + if components_by_species: + #more species with components to add to this block + components_by_species = deepcopy( components_by_species ) + spec_comps = components_by_species.pop( 0 ) + for c in spec_comps: + newer_block = deepcopy( new_block ) + newer_block.add_component( deepcopy( c ) ) + for value in __split_components_by_species( components_by_species, newer_block ): + yield value + else: + #no more components to add, yield this block + yield new_block + + #divide components by species + spec_dict = {} + if not species: + species = [] + for c in block.components: + spec, chrom = src_split( c.src ) + if spec not in spec_dict: + spec_dict[ spec ] = [] + species.append( spec ) + spec_dict[ spec ].append( c ) + else: + for spec in species: + spec_dict[ spec ] = [] + for c in iter_components_by_src_start( block, spec ): + spec_dict[ spec ].append( c ) + + empty_block = bx.align.Alignment( score=block.score, attributes=deepcopy( block.attributes ) ) #should we copy attributes? + empty_block.text_size = block.text_size + #call recursive function to split into each combo of spec/blocks + for value in __split_components_by_species( spec_dict.values(), empty_block ): + sort_block_components_by_block( value, block ) #restore original component order + yield value + + #generator yielding only chopped and valid blocks for a specified region -def get_chopped_blocks_for_region( index, src, region, species = None, mincols = 0, force_strand = None ): - for block, idx, offset in get_chopped_blocks_with_index_offset_for_region( index, src, region, species, mincols, force_strand ): +def get_chopped_blocks_for_region( index, src, region, species = None, mincols = 0 ): + for block, idx, offset in get_chopped_blocks_with_index_offset_for_region( index, src, region, species, mincols ): yield block -def get_chopped_blocks_with_index_offset_for_region( index, src, region, species = None, mincols = 0, force_strand = None ): +def get_chopped_blocks_with_index_offset_for_region( index, src, region, species = None, mincols = 0 ): for block, idx, offset in index.get_as_iterator_with_index_and_offset( src, region.start, region.end ): - block = chop_block_by_region( block, src, region, species, mincols, force_strand ) + block = chop_block_by_region( block, src, region, species, mincols ) if block is not None: yield block, idx, offset #returns a filled region alignment for specified regions -def get_region_alignment( index, primary_species, chrom, start, end, strand = '+', species = None, mincols = 0 ): +def get_region_alignment( index, primary_species, chrom, start, end, strand = '+', species = None, mincols = 0, overwrite_with_gaps = True ): if species is not None: alignment = RegionAlignment( end - start, species ) else: alignment = RegionAlignment( end - start, primary_species ) - return fill_region_alignment( alignment, index, primary_species, chrom, start, end, strand, species, mincols ) + return fill_region_alignment( alignment, index, primary_species, chrom, start, end, strand, species, mincols, overwrite_with_gaps ) #reduces a block to only positions exisiting in the src provided def reduce_block_by_primary_genome( block, species, chromosome, region_start ): @@ -237,14 +386,12 @@ def reduce_block_by_primary_genome( block, species, chromosome, region_start ): return ( start_offset, species_texts ) #fills a region alignment -def fill_region_alignment( alignment, index, primary_species, chrom, start, end, strand = '+', species = None, mincols = 0 ): +def fill_region_alignment( alignment, index, primary_species, chrom, start, end, strand = '+', species = None, mincols = 0, overwrite_with_gaps = True ): region = bx.intervals.Interval( start, end ) region.chrom = chrom region.strand = strand primary_src = "%s.%s" % ( primary_species, chrom ) - - #Order blocks overlaping this position by score, lowest first blocks = [] for block, idx, offset in index.get_as_iterator_with_index_and_offset( primary_src, start, end ): @@ -255,28 +402,40 @@ def fill_region_alignment( alignment, index, primary_species, chrom, start, end, break else: blocks.append( ( score, idx, offset ) ) - + + gap_chars_tuple = tuple( GAP_CHARS ) + gap_chars_str = ''.join( GAP_CHARS ) #Loop through ordered blocks and layer by increasing score - for block_dict in blocks: - block = chop_block_by_region( block_dict[1].get_at_offset( block_dict[2] ), primary_src, region, species, mincols, strand ) - if block is None: continue - start_offset, species_texts = reduce_block_by_primary_genome( block, primary_species, chrom, start ) - for spec, text in species_texts.items(): - try: - alignment.set_range( start_offset, spec, text ) - except: - #species/sequence for species does not exist - pass - + for block_dict in blocks: for block in iter_blocks_split_by_species( block_dict[1].get_at_offset( block_dict[2] ) ): #need to handle each occurance of sequence in block seperately + if component_overlaps_region( block.get_component_by_src( primary_src ), region ): + block = chop_block_by_region( block, primary_src, region, species, mincols ) #chop block + block = orient_block_by_region( block, primary_src, region ) #orient block + start_offset, species_texts = reduce_block_by_primary_genome( block, primary_species, chrom, start ) + for spec, text in species_texts.items(): + #we should trim gaps from both sides, since these are not positions in this species genome (sequence) + text = text.rstrip( gap_chars_str ) + gap_offset = 0 + while text.startswith( gap_chars_tuple ): + gap_offset += 1 + text = text[1:] + if not text: + break + if text: + if overwrite_with_gaps: + alignment.set_range( start_offset + gap_offset, spec, text ) + else: + for i, char in enumerate( text ): + if char not in GAP_CHARS: + alignment.set_position( start_offset + gap_offset + i, spec, char ) return alignment #returns a filled spliced region alignment for specified region with start and end lists -def get_spliced_region_alignment( index, primary_species, chrom, starts, ends, strand = '+', species = None, mincols = 0 ): +def get_spliced_region_alignment( index, primary_species, chrom, starts, ends, strand = '+', species = None, mincols = 0, overwrite_with_gaps = True ): #create spliced alignment object if species is not None: alignment = SplicedAlignment( starts, ends, species ) else: alignment = SplicedAlignment( starts, ends, [primary_species] ) for exon in alignment.exons: - fill_region_alignment( exon, index, primary_species, chrom, exon.start, exon.end, strand, species, mincols) + fill_region_alignment( exon, index, primary_species, chrom, exon.start, exon.end, strand, species, mincols, overwrite_with_gaps ) return alignment #loop through string array, only return non-commented lines @@ -319,29 +478,36 @@ def get_starts_ends_fields_from_gene_bed( line ): starts.append( start ) ends.append( end ) return ( starts, ends, fields ) - -def get_species_in_maf( maf_filename ): - try: - species={} - - file_in = open( maf_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() - return species - except: - return [] + +def iter_components_by_src( block, src ): + for c in block.components: + if c.src == src: + yield c + +def get_components_by_src( block, src ): + return [ value for value in iter_components_by_src( block, src ) ] + +def iter_components_by_src_start( block, src ): + for c in block.components: + if c.src.startswith( src ): + yield c + +def get_components_by_src_start( block, src ): + return [ value for value in iter_components_by_src_start( block, src ) ] + +def sort_block_components_by_block( block1, block2 ): + #orders the components in block1 by the index of the component in block2 + #block1 must be a subset of block2 + #occurs in-place + return block1.components.sort( cmp = lambda x, y: block2.components.index( x ) - block2.components.index( y ) ) + +def get_species_in_maf( maf_filename ): + species = [] + for block in maf.Reader( open( maf_filename ) ): + for spec in get_species_in_block( block ): + if spec not in species: + species.append( spec ) + return species def parse_species_option( species ): if species: diff --git a/tool_conf.xml.sample b/tool_conf.xml.sample index 5b5a0d1fddf..791327cad18 100644 --- a/tool_conf.xml.sample +++ b/tool_conf.xml.sample @@ -90,6 +90,7 @@
+ diff --git a/tools/annotation_profiler/annotation_profiler_for_interval.py b/tools/annotation_profiler/annotation_profiler_for_interval.py index cd280d2e0e1..1220f61ad39 100644 --- a/tools/annotation_profiler/annotation_profiler_for_interval.py +++ b/tools/annotation_profiler/annotation_profiler_for_interval.py @@ -22,10 +22,12 @@ class CachedRangesInFile: fmt_size = struct.calcsize( fmt ) def __init__( self, filename ): self.file_size = os.stat( filename ).st_size - self.file = open( filename, 'rb' ) + self.file = open( filename, 'rb' ) + self.filename = filename self.length = int( self.file_size / self.fmt_size / 2 ) self._cached_ranges = [ None for i in xrange( self.length ) ] def __getitem__( self, i ): + old_i = i if self._cached_ranges[i] is not None: return self._cached_ranges[i] if i < 0: i = self.length + i @@ -35,7 +37,13 @@ class CachedRangesInFile: start = struct.unpack( self.fmt, self.file.read( self.fmt_size ) )[0] end = struct.unpack( self.fmt, self.file.read( self.fmt_size ) )[0] except Exception, e: - raise IndexError, e + print 'filename', self.filename + print 'len', len( self ) + print 'fmtsize', self.fmt_size + print 'index', i + print 'old i', old_i + print 'offset', offset + raise IndexError( str( e ) ) self._cached_ranges[i] = ( start, end ) return start, end def __len__( self ): @@ -141,20 +149,24 @@ class TableCoverageSummary: self.chromosome_coverage[chrom] = bx.bitset.BitSet( chrom_length ) self.chromosome_coverage[chrom].set_range( region_start, region_length ) - for table_name, coverage, regions in self.coverage_reader.iter_table_coverage_regions_by_region( chrom, region_start, region_end ): - if table_name not in self.table_coverage: - self.table_coverage[table_name] = 0 - self.table_chromosome_size[table_name] = {} - self.table_regions_overlaped_count[table_name] = 0 - self.interval_table_overlap_count[table_name] = 0 - self.table_chromosome_count[table_name] = {} - if chrom not in self.table_chromosome_size[table_name]: - self.table_chromosome_size[table_name][chrom] = self.coverage_reader._coverage[table_name][chrom]._total_coverage - self.table_chromosome_count[table_name][chrom] = len( self.coverage_reader._coverage[table_name][chrom]._coverage ) - self.table_coverage[table_name] += coverage - if coverage: - self.interval_table_overlap_count[table_name] += 1 - self.table_regions_overlaped_count[table_name] += regions + try: + for table_name, coverage, regions in self.coverage_reader.iter_table_coverage_regions_by_region( chrom, region_start, region_end ): + if table_name not in self.table_coverage: + self.table_coverage[table_name] = 0 + self.table_chromosome_size[table_name] = {} + self.table_regions_overlaped_count[table_name] = 0 + self.interval_table_overlap_count[table_name] = 0 + self.table_chromosome_count[table_name] = {} + if chrom not in self.table_chromosome_size[table_name]: + self.table_chromosome_size[table_name][chrom] = self.coverage_reader._coverage[table_name][chrom]._total_coverage + self.table_chromosome_count[table_name][chrom] = len( self.coverage_reader._coverage[table_name][chrom]._coverage ) + self.table_coverage[table_name] += coverage + if coverage: + self.interval_table_overlap_count[table_name] += 1 + self.table_regions_overlaped_count[table_name] += regions + except Exception, e: + print "chrom:%s, start:%s, end%s:." % ( chrom, start, end ) + raise e def iter_table_coverage( self ): def get_nr_coverage(): #returns non-redundant coverage, where user's input intervals have been collapse to resolve overlaps diff --git a/tools/maf/genebed_maf_to_fasta.xml b/tools/maf/genebed_maf_to_fasta.xml index e6897ec8ff5..a458ea80c19 100644 --- a/tools/maf/genebed_maf_to_fasta.xml +++ b/tools/maf/genebed_maf_to_fasta.xml @@ -1,8 +1,8 @@ - + given a set of coding exon intervals #if $maf_source_type.maf_source == "user":#interval_maf_to_merged_fasta.py --dbkey=$dbkey --species=$maf_source_type.species --mafSource=$maf_source_type.maf_file --mafIndex=$maf_source_type.maf_file.metadata.maf_index --interval_file=$input1 --output_file=$out_file1 --mafSourceType=$maf_source_type.maf_source --geneBED --mafIndexFileDir=${GALAXY_DATA_INDEX_DIR} #else:#interval_maf_to_merged_fasta.py --dbkey=$dbkey --species=$maf_source_type.species --mafSource=$maf_source_type.maf_identifier --interval_file=$input1 --output_file=$out_file1 --mafSourceType=$maf_source_type.maf_source --geneBED --mafIndexFileDir=${GALAXY_DATA_INDEX_DIR} -#end if +#end if# --overwrite_with_gaps=$overwrite_with_gaps @@ -29,42 +29,48 @@ - - - - - - + + + + + + - - - - - + + + + + - + + + + + - + in aligning species - + + - + + diff --git a/tools/maf/interval2maf.py b/tools/maf/interval2maf.py index 5baaaa90ccc..130f2ee3555 100755 --- a/tools/maf/interval2maf.py +++ b/tools/maf/interval2maf.py @@ -20,6 +20,8 @@ usage: %prog maf_file [options] -i, --interval_file=i: Input interval file -o, --output_file=o: Output MAF file -p, --species=p: Species to include in output + -P, --split_blocks_by_species=P: Split blocks by species + -r, --remove_all_gap_columns=r: Remove all Gap columns -l, --indexLocation=l: Override default maf_index.loc file -z, --mafIndexFile=z: Directory of local maf index file ( maf_index.loc or maf_pairwise.loc ) """ @@ -45,25 +47,21 @@ def __main__(): if options.dbkey: dbkey = options.dbkey else: dbkey = None if dbkey in [None, "?"]: - print >>sys.stderr, "You must specify a proper build in order to extract alignments. You can specify your genome build by clicking on the pencil icon associated with your interval file." - sys.exit() + maf_utilities.tool_fail( "You must specify a proper build in order to extract alignments. You can specify your genome build by clicking on the pencil icon associated with your interval file." ) species = maf_utilities.parse_species_option( options.species ) if options.chromCol: chromCol = int( options.chromCol ) - 1 else: - print >>sys.stderr, "Chromosome column not set, click the pencil icon in the history item to set the metadata attributes." - sys.exit() + maf_utilities.tool_fail( "Chromosome column not set, click the pencil icon in the history item to set the metadata attributes." ) if options.startCol: startCol = int( options.startCol ) - 1 else: - print >>sys.stderr, "Start column not set, click the pencil icon in the history item to set the metadata attributes." - sys.exit() + maf_utilities.tool_fail( "Start column not set, click the pencil icon in the history item to set the metadata attributes." ) if options.endCol: endCol = int( options.endCol ) - 1 else: - print >>sys.stderr, "End column not set, click the pencil icon in the history item to set the metadata attributes." - sys.exit() + maf_utilities.tool_fail( "End column not set, click the pencil icon in the history item to set the metadata attributes." ) if options.strandCol: strandCol = int( options.strandCol ) - 1 else: @@ -71,13 +69,17 @@ def __main__(): if options.interval_file: interval_file = options.interval_file else: - print >>sys.stderr, "Input interval file has not been specified." - sys.exit() + maf_utilities.tool_fail( "Input interval file has not been specified." ) if options.output_file: output_file = options.output_file else: - print >>sys.stderr, "Output file has not been specified." - sys.exit() + maf_utilities.tool_fail( "Output file has not been specified." ) + + split_blocks_by_species = remove_all_gap_columns = False + if options.split_blocks_by_species and options.split_blocks_by_species == 'split_blocks_by_species': + split_blocks_by_species = True + if options.remove_all_gap_columns and options.remove_all_gap_columns == 'remove_all_gap_columns': + remove_all_gap_columns = True #Finish parsing command line #Open indexed access to MAFs @@ -87,16 +89,13 @@ def __main__(): else: index = maf_utilities.maf_index_by_uid( options.mafType, options.mafIndexFile ) if index is None: - print >> sys.stderr, "The MAF source specified (%s) appears to be invalid." % ( options.mafType ) - sys.exit() + maf_utilities.tool_fail( "The MAF source specified (%s) appears to be invalid." % ( options.mafType ) ) elif options.mafFile: index, index_filename = maf_utilities.open_or_build_maf_index( options.mafFile, options.mafIndex, species = [dbkey] ) if index is None: - print >> sys.stderr, "Your MAF file appears to be malformed." - sys.exit() + maf_utilities.tool_fail( "Your MAF file appears to be malformed." ) else: - print >>sys.stderr, "Desired source MAF type has not been specified." - sys.exit() + maf_utilities.tool_fail( "Desired source MAF type has not been specified." ) #Create MAF writter out = bx.align.maf.Writer( open(output_file, "w") ) @@ -105,10 +104,20 @@ def __main__(): num_blocks = 0 num_regions = None for num_regions, region in enumerate( bx.intervals.io.NiceReaderWrapper( open( interval_file, 'r' ), chrom_col = chromCol, start_col = startCol, end_col = endCol, strand_col = strandCol, fix_strand = True, return_header = False, return_comments = False ) ): - src = "%s.%s" % ( dbkey, region.chrom ) - for block in maf_utilities.get_chopped_blocks_for_region( index, src, region, species, mincols ): - out.write( block ) - num_blocks += 1 + src = maf_utilities.src_merge( dbkey, region.chrom ) + for block in index.get_as_iterator( src, region.start, region.end ): + if split_blocks_by_species: + blocks = [ new_block for new_block in maf_utilities.iter_blocks_split_by_species( block ) if maf_utilities.component_overlaps_region( new_block.get_component_by_src_start( dbkey ), region ) ] + else: + blocks = [ block ] + for block in blocks: + block = maf_utilities.chop_block_by_region( block, src, region ) + if block is not None: + block = maf_utilities.orient_block_by_region( block, src, region ) + if remove_all_gap_columns: + block.remove_all_gap_columns() + out.write( block ) + num_blocks += 1 #Close output MAF out.close() diff --git a/tools/maf/interval2maf.xml b/tools/maf/interval2maf.xml index 2cd164ec12a..5c82d50a5a3 100644 --- a/tools/maf/interval2maf.xml +++ b/tools/maf/interval2maf.xml @@ -3,6 +3,9 @@ #if $maf_source_type.maf_source == "user":#interval2maf.py --dbkey=${input1.dbkey} --chromCol=${input1.metadata.chromCol} --startCol=${input1.metadata.startCol} --endCol=${input1.metadata.endCol} --strandCol=${input1.metadata.strandCol} --mafFile=$maf_source_type.mafFile --mafIndex=$maf_source_type.mafFile.metadata.maf_index --interval_file=$input1 --output_file=$out_file1 --mafIndexFile=${GALAXY_DATA_INDEX_DIR}/maf_index.loc --species=$maf_source_type.species #else:#interval2maf.py --dbkey=${input1.dbkey} --chromCol=${input1.metadata.chromCol} --startCol=${input1.metadata.startCol} --endCol=${input1.metadata.endCol} --strandCol=${input1.metadata.strandCol} --mafType=$maf_source_type.mafType --interval_file=$input1 --output_file=$out_file1 --mafIndexFile=${GALAXY_DATA_INDEX_DIR}/maf_index.loc --species=$maf_source_type.species + #end if + --split_blocks_by_species=$split_blocks_by_species_selector.split_blocks_by_species + #if $split_blocks_by_species_selector.split_blocks_by_species == "split_blocks_by_species":# --remove_all_gap_columns=$split_blocks_by_species_selector.remove_all_gap_columns #end if @@ -48,7 +51,22 @@ - + + + + + + + + + + + + + + + + @@ -59,6 +77,7 @@ + @@ -66,6 +85,7 @@ + diff --git a/tools/maf/interval2maf_pairwise.xml b/tools/maf/interval2maf_pairwise.xml index a97401ae68f..b4270840c3e 100644 --- a/tools/maf/interval2maf_pairwise.xml +++ b/tools/maf/interval2maf_pairwise.xml @@ -1,4 +1,4 @@ - + given a set of genomic intervals interval2maf.py --dbkey=${input1.dbkey} --chromCol=${input1.metadata.chromCol} --startCol=${input1.metadata.startCol} --endCol=${input1.metadata.endCol} --strandCol=${input1.metadata.strandCol} --mafType=$mafType --interval_file=$input1 --output_file=$out_file1 --indexLocation=${GALAXY_DATA_INDEX_DIR}/maf_pairwise.loc diff --git a/tools/maf/interval_maf_to_merged_fasta.py b/tools/maf/interval_maf_to_merged_fasta.py index e7de389d3e3..56edbd265e6 100644 --- a/tools/maf/interval_maf_to_merged_fasta.py +++ b/tools/maf/interval_maf_to_merged_fasta.py @@ -19,6 +19,7 @@ usage: %prog maf_file [options] -i, --interval_file=i: Input interval file -o, --output_file=o: Output MAF file -p, --species=p: Species to include in output + -O, --overwrite_with_gaps=O: Overwrite bases found in a lower-scoring block with gaps interior to the sequence for a species. -z, --mafIndexFileDir=z: Directory of local maf_index.loc file usage: %prog dbkey_of_BED comma_separated_list_of_additional_dbkeys_to_extract comma_separated_list_of_indexed_maf_files input_gene_bed_file output_fasta_file cached|user GALAXY_DATA_INDEX_DIR @@ -93,6 +94,11 @@ def __main__(): strand_col = int( options.strandCol ) - 1 mafIndexFile = "%s/maf_index.loc" % options.mafIndexFileDir + + overwrite_with_gaps = True + if options.overwrite_with_gaps and options.overwrite_with_gaps.lower() == 'false': + overwrite_with_gaps = False + #Finish parsing command line #get index for mafs based on type @@ -127,7 +133,7 @@ def __main__(): try: starts, ends, fields = maf_utilities.get_starts_ends_fields_from_gene_bed( line ) #create spliced alignment object - alignment = maf_utilities.get_spliced_region_alignment( index, primary_species, fields[0], starts, ends, strand = '+', species = species, mincols = mincols ) + alignment = maf_utilities.get_spliced_region_alignment( index, primary_species, fields[0], starts, ends, strand = '+', species = species, mincols = mincols, overwrite_with_gaps = overwrite_with_gaps ) primary_name = secondary_name = fields[3] alignment_strand = fields[5] except Exception, e: @@ -136,7 +142,7 @@ def __main__(): else: #Process as standard intervals try: #create spliced alignment object - alignment = maf_utilities.get_region_alignment( index, primary_species, line.chrom, line.start, line.end, strand = '+', species = species, mincols = mincols ) + alignment = maf_utilities.get_region_alignment( index, primary_species, line.chrom, line.start, line.end, strand = '+', species = species, mincols = mincols, overwrite_with_gaps = overwrite_with_gaps ) primary_name = "%s(%s):%s-%s" % ( line.chrom, line.strand, line.start, line.end ) secondary_name = "" alignment_strand = line.strand diff --git a/tools/maf/interval_maf_to_merged_fasta.xml b/tools/maf/interval_maf_to_merged_fasta.xml index facfe79d458..d4e2d7cebe6 100644 --- a/tools/maf/interval_maf_to_merged_fasta.xml +++ b/tools/maf/interval_maf_to_merged_fasta.xml @@ -1,8 +1,8 @@ - + given a set of genomic intervals #if $maf_source_type.maf_source == "user":#interval_maf_to_merged_fasta.py --dbkey=$dbkey --species=$maf_source_type.species --mafSource=$maf_source_type.maf_file --mafIndex=$maf_source_type.maf_file.metadata.maf_index --interval_file=$input1 --output_file=$out_file1 --chromCol=${input1.metadata.chromCol} --startCol=${input1.metadata.startCol} --endCol=${input1.metadata.endCol} --strandCol=${input1.metadata.strandCol} --mafSourceType=$maf_source_type.maf_source --mafIndexFileDir=${GALAXY_DATA_INDEX_DIR} #else:#interval_maf_to_merged_fasta.py --dbkey=$dbkey --species=$maf_source_type.species --mafSource=$maf_source_type.maf_identifier --interval_file=$input1 --output_file=$out_file1 --chromCol=${input1.metadata.chromCol} --startCol=${input1.metadata.startCol} --endCol=${input1.metadata.endCol} --strandCol=${input1.metadata.strandCol} --mafSourceType=$maf_source_type.maf_source --mafIndexFileDir=${GALAXY_DATA_INDEX_DIR} -#end if +#end if# --overwrite_with_gaps=$overwrite_with_gaps @@ -30,25 +30,29 @@ - - - - - - + + + + + + - - - - - + + + + + + + + + @@ -59,21 +63,24 @@ - + + - + + - + + diff --git a/tools/maf/maf_filter.py b/tools/maf/maf_filter.py index ea9b8334a9a..c223f2912ff 100644 --- a/tools/maf/maf_filter.py +++ b/tools/maf/maf_filter.py @@ -46,7 +46,7 @@ def main(): i = 0 blocks_kept = 0 for i, maf_block in enumerate( maf_reader ): - if min_size <= maf_block.components[0].size <= max_size: + if min_size <= maf_block.text_size <= max_size: local = {'maf_block':maf_block, 'ret_val':False} execfile( script_file, {}, local ) if local['ret_val']: diff --git a/tools/maf/maf_limit_size.py b/tools/maf/maf_limit_size.py index 91739306f24..2b0b33227f1 100644 --- a/tools/maf/maf_limit_size.py +++ b/tools/maf/maf_limit_size.py @@ -28,7 +28,7 @@ def __main__(): blocks_kept = 0 i = 0 for i, m in enumerate( maf_reader ): - if min_size <= m.components[0].size <= max_size: + if min_size <= m.text_size <= max_size: maf_writer.write( m ) blocks_kept += 1 print 'Kept %s of %s blocks (%.2f%%).' % ( blocks_kept, i + 1, float( blocks_kept ) / float( i + 1 ) * 100.0 ) diff --git a/tools/maf/maf_limit_size.xml b/tools/maf/maf_limit_size.xml index e2d36cab1db..5c5a8d85e6a 100644 --- a/tools/maf/maf_limit_size.xml +++ b/tools/maf/maf_limit_size.xml @@ -1,4 +1,4 @@ - + by Size maf_limit_size.py $input1 $out_file1 $min_size $max_size diff --git a/tools/maf/maf_limit_to_species.py b/tools/maf/maf_limit_to_species.py index a13a760ec6e..781b6038b44 100644 --- a/tools/maf/maf_limit_to_species.py +++ b/tools/maf/maf_limit_to_species.py @@ -11,13 +11,18 @@ usage: %prog species,species2,... input_maf output_maf allow_partial min_species from galaxy import eggs import pkg_resources; pkg_resources.require( "bx-python" ) import bx.align.maf +from galaxy.tools.util import maf_utilities import sys assert sys.version_info[:2] >= ( 2, 4 ) def main(): - species = sys.argv[1].split( ',' ) + species = maf_utilities.parse_species_option( sys.argv[1] ) + if species: + spec_len = len( species ) + else: + spec_len = 0 try: maf_reader = bx.align.maf.Reader( open( sys.argv[2],'r' ) ) maf_writer = bx.align.maf.Writer( open( sys.argv[3],'w' ) ) @@ -30,10 +35,11 @@ def main(): maf_blocks_kept = 0 for m in maf_reader: - if species != ['None']: + if species: m = m.limit_to_species( species ) m.remove_all_gap_columns() - if ( species == ['None'] or allow_partial or len( m.components ) == len( species ) ) and len( m.components ) > min_species_per_block: + spec_in_block_len = len( maf_utilities.get_species_in_block( m ) ) + if ( not species or allow_partial or spec_in_block_len == spec_len ) and spec_in_block_len > min_species_per_block: maf_writer.write( m ) maf_blocks_kept += 1 diff --git a/tools/maf/maf_split_by_species.py b/tools/maf/maf_split_by_species.py new file mode 100644 index 00000000000..076e644ce22 --- /dev/null +++ b/tools/maf/maf_split_by_species.py @@ -0,0 +1,44 @@ +#!/usr/bin/env python + +""" +Read a maf and split blocks by unique species combinations +""" +import sys +from galaxy import eggs +import pkg_resources; pkg_resources.require( "bx-python" ) +from bx.align import maf +from galaxy.tools.util import maf_utilities +from galaxy.util import string_as_bool + +assert sys.version_info[:2] >= ( 2, 4 ) + +def __main__(): + try: + maf_reader = maf.Reader( open( sys.argv[1] ) ) + except Exception, e: + maf_utilities.tool_fail( "Error opening MAF: %s" % e ) + try: + out = maf.Writer( open( sys.argv[2], "w") ) + except Exception, e: + maf_utilities.tool_fail( "Error opening file for output: %s" % e ) + try: + collapse_columns = string_as_bool( sys.argv[3] ) + except Exception, e: + maf_utilities.tool_fail( "Error determining collapse columns value: %s" % e ) + + start_count = 0 + end_count = 0 + for start_count, start_block in enumerate( maf_reader ): + for block in maf_utilities.iter_blocks_split_by_species( start_block ): + if collapse_columns: + block.remove_all_gap_columns() + out.write( block ) + end_count += 1 + out.close() + + if end_count: + print "%i alignment blocks created from %i original blocks." % ( end_count, start_count + 1 ) + else: + print "No alignment blocks were created." + +if __name__ == "__main__": __main__() diff --git a/tools/maf/maf_split_by_species.xml b/tools/maf/maf_split_by_species.xml new file mode 100644 index 00000000000..4a594da62c6 --- /dev/null +++ b/tools/maf/maf_split_by_species.xml @@ -0,0 +1,216 @@ + + by Species + maf_split_by_species.py $input1 $out_file1 $collapse_columns + + + + + + + + + + + + + + + + + + + + + + + + +**What it does** + +This tool examines each MAF block for multiple occurrences of a species in a single block. When this occurs, a block is split into multiple blocks where every combination of one sequence per species per block is represented. + +The interface for this tool has two inputs: + + * **MAF file to split**. Choose multiple alignments from history to be split by species. + * **Collapse empty alignment columns**. Should alignment columns containing only gaps in the new blocks be removed. + +----- + +**Example 1**: **Collapse empty alignment columns is Yes**: + +For the following alignment:: + + ##maf version=1 + a score=2047408.0 + s species1.chr1 147984545 85 + 245522847 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCTTCGGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTTGTCCTCAG + s species1.chr1 147984545 83 - 245522847 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCTT--GTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTTGTCCTCAG + s species1.chr1 147984645 79 + 245522847 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCTTCGGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTT------AG + s species1.chr1 147984645 79 - 245522847 ATGGCGTCGGCCTCCTCCGGGCCGTCGTC---GGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTTGTC---AG + s species2.chr1 129723125 85 + 229575298 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCTTCGGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTCGTCCTCAG + s species2.chr1 129723125 83 - 229575298 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCT--GGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTCGTCCTCAG + s species2.chr1 129723925 79 + 229575298 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCTTCGGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTC------AG + s species3.chr3 68255714 76 - 258222147 ATGGCGTCCGCCTCCTCAGGGCCAGCGGC---GGCGGGGTTTTCACCCCTTGATTCCGGGGTCCCTGCCGGTACCGC------AG + +the tool will create **a single** history item containing 12 alignment blocks (notice that no columns contain only gaps):: + + ##maf version=1 + a score=2047408.0 + s species1.chr1 147984545 85 + 245522847 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCTTCGGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTTGTCCTCAG + s species2.chr1 129723125 85 + 229575298 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCTTCGGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTCGTCCTCAG + s species3.chr3 68255714 76 - 258222147 ATGGCGTCCGCCTCCTCAGGGCCAGCGGC---GGCGGGGTTTTCACCCCTTGATTCCGGGGTCCCTGCCGGTACCGC------AG + + a score=2047408.0 + s species1.chr1 147984545 83 - 245522847 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCTT--GTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTTGTCCTCAG + s species2.chr1 129723125 85 + 229575298 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCTTCGGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTCGTCCTCAG + s species3.chr3 68255714 76 - 258222147 ATGGCGTCCGCCTCCTCAGGGCCAGCGGC---GGCGGGGTTTTCACCCCTTGATTCCGGGGTCCCTGCCGGTACCGC------AG + + a score=2047408.0 + s species1.chr1 147984645 79 + 245522847 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCTTCGGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTT------AG + s species2.chr1 129723125 85 + 229575298 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCTTCGGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTCGTCCTCAG + s species3.chr3 68255714 76 - 258222147 ATGGCGTCCGCCTCCTCAGGGCCAGCGGC---GGCGGGGTTTTCACCCCTTGATTCCGGGGTCCCTGCCGGTACCGC------AG + + a score=2047408.0 + s species1.chr1 147984645 79 - 245522847 ATGGCGTCGGCCTCCTCCGGGCCGTCGTC---GGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTTGTC---AG + s species2.chr1 129723125 85 + 229575298 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCTTCGGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTCGTCCTCAG + s species3.chr3 68255714 76 - 258222147 ATGGCGTCCGCCTCCTCAGGGCCAGCGGC---GGCGGGGTTTTCACCCCTTGATTCCGGGGTCCCTGCCGGTACCGC------AG + + a score=2047408.0 + s species1.chr1 147984545 85 + 245522847 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCTTCGGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTTGTCCTCAG + s species2.chr1 129723125 83 - 229575298 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCT--GGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTCGTCCTCAG + s species3.chr3 68255714 76 - 258222147 ATGGCGTCCGCCTCCTCAGGGCCAGCGGC---GGCGGGGTTTTCACCCCTTGATTCCGGGGTCCCTGCCGGTACCGC------AG + + a score=2047408.0 + s species1.chr1 147984545 83 - 245522847 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCTT-GTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTTGTCCTCAG + s species2.chr1 129723125 83 - 229575298 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCT-GGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTCGTCCTCAG + s species3.chr3 68255714 76 - 258222147 ATGGCGTCCGCCTCCTCAGGGCCAGCGGC--GGCGGGGTTTTCACCCCTTGATTCCGGGGTCCCTGCCGGTACCGC------AG + + a score=2047408.0 + s species1.chr1 147984645 79 + 245522847 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCTTCGGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTT------AG + s species2.chr1 129723125 83 - 229575298 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCT--GGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTCGTCCTCAG + s species3.chr3 68255714 76 - 258222147 ATGGCGTCCGCCTCCTCAGGGCCAGCGGC---GGCGGGGTTTTCACCCCTTGATTCCGGGGTCCCTGCCGGTACCGC------AG + + a score=2047408.0 + s species1.chr1 147984645 79 - 245522847 ATGGCGTCGGCCTCCTCCGGGCCGTCGTC-GGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTTGTC---AG + s species2.chr1 129723125 83 - 229575298 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCTGGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTCGTCCTCAG + s species3.chr3 68255714 76 - 258222147 ATGGCGTCCGCCTCCTCAGGGCCAGCGGC-GGCGGGGTTTTCACCCCTTGATTCCGGGGTCCCTGCCGGTACCGC------AG + + a score=2047408.0 + s species1.chr1 147984545 85 + 245522847 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCTTCGGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTTGTCCTCAG + s species2.chr1 129723925 79 + 229575298 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCTTCGGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTC------AG + s species3.chr3 68255714 76 - 258222147 ATGGCGTCCGCCTCCTCAGGGCCAGCGGC---GGCGGGGTTTTCACCCCTTGATTCCGGGGTCCCTGCCGGTACCGC------AG + + a score=2047408.0 + s species1.chr1 147984545 83 - 245522847 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCTT--GTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTTGTCCTCAG + s species2.chr1 129723925 79 + 229575298 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCTTCGGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTC------AG + s species3.chr3 68255714 76 - 258222147 ATGGCGTCCGCCTCCTCAGGGCCAGCGGC---GGCGGGGTTTTCACCCCTTGATTCCGGGGTCCCTGCCGGTACCGC------AG + + a score=2047408.0 + s species1.chr1 147984645 79 + 245522847 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCTTCGGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTTAG + s species2.chr1 129723925 79 + 229575298 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCTTCGGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTCAG + s species3.chr3 68255714 76 - 258222147 ATGGCGTCCGCCTCCTCAGGGCCAGCGGC---GGCGGGGTTTTCACCCCTTGATTCCGGGGTCCCTGCCGGTACCGCAG + + a score=2047408.0 + s species1.chr1 147984645 79 - 245522847 ATGGCGTCGGCCTCCTCCGGGCCGTCGTC---GGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTTGTCAG + s species2.chr1 129723925 79 + 229575298 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCTTCGGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTC---AG + s species3.chr3 68255714 76 - 258222147 ATGGCGTCCGCCTCCTCAGGGCCAGCGGC---GGCGGGGTTTTCACCCCTTGATTCCGGGGTCCCTGCCGGTACCGC---AG + +----- + +**Example 1**: **Collapse empty alignment columns is Yes**: + +For the following alignment:: + + ##maf version=1 + a score=2047408.0 + s species1.chr1 147984545 85 + 245522847 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCTTCGGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTTGTCCTCAG + s species1.chr1 147984545 83 - 245522847 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCTT--GTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTTGTCCTCAG + s species1.chr1 147984645 79 + 245522847 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCTTCGGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTT------AG + s species1.chr1 147984645 79 - 245522847 ATGGCGTCGGCCTCCTCCGGGCCGTCGTC---GGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTTGTC---AG + s species2.chr1 129723125 85 + 229575298 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCTTCGGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTCGTCCTCAG + s species2.chr1 129723125 83 - 229575298 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCT--GGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTCGTCCTCAG + s species2.chr1 129723925 79 + 229575298 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCTTCGGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTC------AG + s species3.chr3 68255714 76 - 258222147 ATGGCGTCCGCCTCCTCAGGGCCAGCGGC---GGCGGGGTTTTCACCCCTTGATTCCGGGGTCCCTGCCGGTACCGC------AG + +the tool will create **a single** history item containing 12 alignment blocks (notice that some columns contain only gaps):: + + ##maf version=1 + a score=2047408.0 + s species1.chr1 147984545 85 + 245522847 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCTTCGGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTTGTCCTCAG + s species2.chr1 129723125 85 + 229575298 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCTTCGGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTCGTCCTCAG + s species3.chr3 68255714 76 - 258222147 ATGGCGTCCGCCTCCTCAGGGCCAGCGGC---GGCGGGGTTTTCACCCCTTGATTCCGGGGTCCCTGCCGGTACCGC------AG + + a score=2047408.0 + s species1.chr1 147984545 83 - 245522847 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCTT--GTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTTGTCCTCAG + s species2.chr1 129723125 85 + 229575298 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCTTCGGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTCGTCCTCAG + s species3.chr3 68255714 76 - 258222147 ATGGCGTCCGCCTCCTCAGGGCCAGCGGC---GGCGGGGTTTTCACCCCTTGATTCCGGGGTCCCTGCCGGTACCGC------AG + + a score=2047408.0 + s species1.chr1 147984645 79 + 245522847 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCTTCGGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTT------AG + s species2.chr1 129723125 85 + 229575298 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCTTCGGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTCGTCCTCAG + s species3.chr3 68255714 76 - 258222147 ATGGCGTCCGCCTCCTCAGGGCCAGCGGC---GGCGGGGTTTTCACCCCTTGATTCCGGGGTCCCTGCCGGTACCGC------AG + + a score=2047408.0 + s species1.chr1 147984645 79 - 245522847 ATGGCGTCGGCCTCCTCCGGGCCGTCGTC---GGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTTGTC---AG + s species2.chr1 129723125 85 + 229575298 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCTTCGGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTCGTCCTCAG + s species3.chr3 68255714 76 - 258222147 ATGGCGTCCGCCTCCTCAGGGCCAGCGGC---GGCGGGGTTTTCACCCCTTGATTCCGGGGTCCCTGCCGGTACCGC------AG + + a score=2047408.0 + s species1.chr1 147984545 85 + 245522847 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCTTCGGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTTGTCCTCAG + s species2.chr1 129723125 83 - 229575298 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCT--GGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTCGTCCTCAG + s species3.chr3 68255714 76 - 258222147 ATGGCGTCCGCCTCCTCAGGGCCAGCGGC---GGCGGGGTTTTCACCCCTTGATTCCGGGGTCCCTGCCGGTACCGC------AG + + a score=2047408.0 + s species1.chr1 147984545 83 - 245522847 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCTT--GTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTTGTCCTCAG + s species2.chr1 129723125 83 - 229575298 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCT--GGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTCGTCCTCAG + s species3.chr3 68255714 76 - 258222147 ATGGCGTCCGCCTCCTCAGGGCCAGCGGC---GGCGGGGTTTTCACCCCTTGATTCCGGGGTCCCTGCCGGTACCGC------AG + + a score=2047408.0 + s species1.chr1 147984645 79 + 245522847 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCTTCGGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTT------AG + s species2.chr1 129723125 83 - 229575298 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCT--GGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTCGTCCTCAG + s species3.chr3 68255714 76 - 258222147 ATGGCGTCCGCCTCCTCAGGGCCAGCGGC---GGCGGGGTTTTCACCCCTTGATTCCGGGGTCCCTGCCGGTACCGC------AG + + a score=2047408.0 + s species1.chr1 147984645 79 - 245522847 ATGGCGTCGGCCTCCTCCGGGCCGTCGTC---GGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTTGTC---AG + s species2.chr1 129723125 83 - 229575298 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCT--GGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTCGTCCTCAG + s species3.chr3 68255714 76 - 258222147 ATGGCGTCCGCCTCCTCAGGGCCAGCGGC---GGCGGGGTTTTCACCCCTTGATTCCGGGGTCCCTGCCGGTACCGC------AG + + a score=2047408.0 + s species1.chr1 147984545 85 + 245522847 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCTTCGGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTTGTCCTCAG + s species2.chr1 129723925 79 + 229575298 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCTTCGGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTC------AG + s species3.chr3 68255714 76 - 258222147 ATGGCGTCCGCCTCCTCAGGGCCAGCGGC---GGCGGGGTTTTCACCCCTTGATTCCGGGGTCCCTGCCGGTACCGC------AG + + a score=2047408.0 + s species1.chr1 147984545 83 - 245522847 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCTT--GTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTTGTCCTCAG + s species2.chr1 129723925 79 + 229575298 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCTTCGGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTC------AG + s species3.chr3 68255714 76 - 258222147 ATGGCGTCCGCCTCCTCAGGGCCAGCGGC---GGCGGGGTTTTCACCCCTTGATTCCGGGGTCCCTGCCGGTACCGC------AG + + a score=2047408.0 + s species1.chr1 147984645 79 + 245522847 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCTTCGGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTT------AG + s species2.chr1 129723925 79 + 229575298 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCTTCGGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTC------AG + s species3.chr3 68255714 76 - 258222147 ATGGCGTCCGCCTCCTCAGGGCCAGCGGC---GGCGGGGTTTTCACCCCTTGATTCCGGGGTCCCTGCCGGTACCGC------AG + + a score=2047408.0 + s species1.chr1 147984645 79 - 245522847 ATGGCGTCGGCCTCCTCCGGGCCGTCGTC---GGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTTGTC---AG + s species2.chr1 129723925 79 + 229575298 ATGGCGTCGGCCTCCTCCGGGCCGTCGTCTTCGGTCGGTTTTTCATCCTTTGATCCCGCGGTCCCTTCCTGTACCTC------AG + s species3.chr3 68255714 76 - 258222147 ATGGCGTCCGCCTCCTCAGGGCCAGCGGC---GGCGGGGTTTTCACCCCTTGATTCCGGGGTCCCTGCCGGTACCGC------AG + +------- + +.. class:: infomark + +**About formats** + +**MAF format** multiple alignment format file. This format stores multiple alignments at the DNA level between entire genomes. + + - The .maf format is line-oriented. Each multiple alignment ends with a blank line. + - Each sequence in an alignment is on a single line. + - Lines starting with # are considered to be comments. + - Each multiple alignment is in a separate paragraph that begins with an "a" line and contains an "s" line for each sequence in the multiple alignment. + - Some MAF files may contain two optional line types: + + - An "i" line containing information about what is in the aligned species DNA before and after the immediately preceding "s" line; + - An "e" line containing information about the size of the gap between the alignments that span the current block. + + + + diff --git a/tools/maf/maf_stats.py b/tools/maf/maf_stats.py index 8f38e9bcbb7..0a5aa61c507 100644 --- a/tools/maf/maf_stats.py +++ b/tools/maf/maf_stats.py @@ -64,16 +64,19 @@ def __main__(): total_length += region_length coverage = { dbkey: BitSet( region_length ) } - for block in maf_utilities.get_chopped_blocks_for_region( index, src, region, force_strand='+' ): - #make sure all species are known - for c in block.components: - spec = c.src.split( '.' )[0] + + for block in index.get_as_iterator( src, region.start, region.end ): + for spec in maf_utilities.get_species_in_block( block ): if spec not in coverage: coverage[spec] = BitSet( region_length ) - start_offset, alignment = maf_utilities.reduce_block_by_primary_genome( block, dbkey, region.chrom, region.start ) - for i in range( len( alignment[dbkey] ) ): - for spec, text in alignment.items(): - if text[i] != '-': - coverage[spec].set( start_offset + i ) + for block in maf_utilities.iter_blocks_split_by_species( block ): + if maf_utilities.component_overlaps_region( block.get_component_by_src( src ), region ): + #need to chop and orient the block + block = maf_utilities.orient_block_by_region( maf_utilities.chop_block_by_region( block, src, region ), src, region, force_strand = '+' ) + start_offset, alignment = maf_utilities.reduce_block_by_primary_genome( block, dbkey, region.chrom, region.start ) + for i in range( len( alignment[dbkey] ) ): + for spec, text in alignment.items(): + if text[i] != '-': + coverage[spec].set( start_offset + i ) if summary: #record summary for key in coverage.keys(): diff --git a/tools/maf/maf_stats.xml b/tools/maf/maf_stats.xml index fe89b43f291..59d745fe8fb 100644 --- a/tools/maf/maf_stats.xml +++ b/tools/maf/maf_stats.xml @@ -1,4 +1,4 @@ - + Alignment coverage information maf_stats.py diff --git a/tools/maf/maf_to_fasta.xml b/tools/maf/maf_to_fasta.xml index 841ac024820..3f6ef2138ff 100644 --- a/tools/maf/maf_to_fasta.xml +++ b/tools/maf/maf_to_fasta.xml @@ -1,4 +1,4 @@ - + Converts a MAF formated file to FASTA format #if $fasta_target_type.fasta_type == "multiple":#maf_to_fasta_multiple_sets.py $input1 $out_file1 $fasta_target_type.species $fasta_target_type.complete_blocks @@ -47,7 +47,7 @@ - + diff --git a/tools/maf/maf_to_fasta_concat.py b/tools/maf/maf_to_fasta_concat.py index 3aac8f6b8ff..25665b9d17e 100755 --- a/tools/maf/maf_to_fasta_concat.py +++ b/tools/maf/maf_to_fasta_concat.py @@ -1,7 +1,7 @@ #!/usr/bin/env python """ -Read a maf and print the text as a fasta file, concatenating blocks +Read a maf and output a single block fasta file, concatenating blocks usage %prog species1,species2 maf_file out_file """ @@ -15,28 +15,43 @@ from galaxy.tools.util import maf_utilities assert sys.version_info[:2] >= ( 2, 4 ) def __main__(): - print "Restricted to species:", sys.argv[1] + try: + species = maf_utilities.parse_species_option( sys.argv[1] ) + except Exception, e: + maf_utilities.tool_fail( "Error determining species value: %s" % e ) + try: + input_filename = sys.argv[2] + except Exception, e: + maf_utilities.tool_fail( "Error reading MAF filename: %s" % e ) + try: + file_out = open( sys.argv[3], 'w' ) + except Exception, e: + maf_utilities.tool_fail( "Error opening file for output: %s" % e ) - texts = {} + if species: + print "Restricted to species: %s" % ', '.join( species ) + else: + print "Not restricted to species." - input_filename = sys.argv[2] - output_filename = sys.argv[3] - species = sys.argv[1].split( ',' ) + if not species: + try: + species = maf_utilities.get_species_in_maf( input_filename ) + except Exception, e: + maf_utilities.tool_fail( "Error determining species in input MAF: %s" % e ) - if "None" in species: - species = maf_utilities.get_species_in_maf( input_filename ) - - file_out = open( output_filename, 'w' ) for spec in species: file_out.write( ">" + spec + "\n" ) try: - for block in maf.Reader( open( input_filename, 'r' ) ): - component = block.get_component_by_src_start( spec ) - if component: file_out.write( component.text ) - else: file_out.write( "-" * block.text_size ) - except: - print >>sys.stderr, "Your MAF file appears to be malformed." - sys.exit() + for start_block in maf.Reader( open( input_filename, 'r' ) ): + for block in maf_utilities.iter_blocks_split_by_species( start_block ): + block.remove_all_gap_columns() #remove extra gaps + component = block.get_component_by_src_start( spec ) #blocks only have one occurrence of a particular species, so this is safe + if component: + file_out.write( component.text ) + else: + file_out.write( "-" * block.text_size ) + except Exception, e: + maf_utilities.tool_fail( "Your MAF file appears to be malformed: %s" % e ) file_out.write( "\n" ) file_out.close() diff --git a/tools/maf/maf_to_fasta_multiple_sets.py b/tools/maf/maf_to_fasta_multiple_sets.py index d48b75041cc..b8b7b93ee13 100755 --- a/tools/maf/maf_to_fasta_multiple_sets.py +++ b/tools/maf/maf_to_fasta_multiple_sets.py @@ -1,7 +1,7 @@ #!/usr/bin/env python """ -Read a maf and print the text as a fasta file. +Read a maf and output a multiple block fasta file. """ #Dan Blankenberg import sys @@ -13,35 +13,46 @@ from galaxy.tools.util import maf_utilities assert sys.version_info[:2] >= ( 2, 4 ) def __main__(): - print "Restricted to species:", sys.argv[3] - - input_filename = sys.argv[1] - output_filename = sys.argv[2] - species = sys.argv[3].split( ',' ) - partial = sys.argv[4] - num_species = len( species ) - - file_in = open( input_filename, 'r' ) try: - maf_reader = maf.Reader( file_in ) - - file_out = open( output_filename, 'w' ) - - for block_num, block in enumerate( maf_reader ): - if "None" not in species: - block = block.limit_to_species( species ) - if len( block.components ) < num_species and partial == "partial_disallowed": continue - for component in block.components: - spec, chrom = maf.src_split( component.src ) - if not spec or not chrom: - spec = chrom = component.src - file_out.write( "%s\n" % maf_utilities.get_fasta_header( component, suffix = "%s_%i" % ( spec, block_num ) ) ) - file_out.write( "%s\n" % component.text ) - file_out.write( "\n" ) - file_in.close() + maf_reader = maf.Reader( open( sys.argv[1] ) ) except Exception, e: - print >>sys.stderr, "Your MAF file appears to be malformed:", e - sys.exit() + maf_utilities.tool_fail( "Error opening input MAF: %s" % e ) + try: + file_out = open( sys.argv[2], 'w' ) + except Exception, e: + maf_utilities.tool_fail( "Error opening file for output: %s" % e ) + try: + species = maf_utilities.parse_species_option( sys.argv[3] ) + if species: + num_species = len( species ) + else: + num_species = 0 + except Exception, e: + maf_utilities.tool_fail( "Error determining species value: %s" % e ) + try: + partial = sys.argv[4] + except Exception, e: + maf_utilities.tool_fail( "Error determining keep partial value: %s" % e ) + + if species: + print "Restricted to species: %s" % ', '.join( species ) + else: + print "Not restricted to species." + + for block_num, block in enumerate( maf_reader ): + if species: + block = block.limit_to_species( species ) + if len( maf_utilities.get_species_in_block( block ) ) < num_species and partial == "partial_disallowed": continue + spec_counts = {} + for component in block.components: + spec, chrom = maf_utilities.src_split( component.src ) + if spec not in spec_counts: + spec_counts[ spec ] = 0 + else: + spec_counts[ spec ] += 1 + file_out.write( "%s\n" % maf_utilities.get_fasta_header( component, { 'block_index' : block_num, 'species' : spec, 'sequence_index' : spec_counts[ spec ] }, suffix = "%s_%i_%i" % ( spec, block_num, spec_counts[ spec ] ) ) ) + file_out.write( "%s\n" % component.text ) + file_out.write( "\n" ) file_out.close() if __name__ == "__main__": __main__()