PageRenderTime 31ms CodeModel.GetById 14ms app.highlight 13ms RepoModel.GetById 2ms app.codeStats 0ms

/tools/maf/maf_to_fasta_multiple_sets.py

https://bitbucket.org/cistrome/cistrome-harvard/
Python | 58 lines | 47 code | 6 blank | 5 comment | 13 complexity | a1854887c47cbf17d0030a6cc46daed5 MD5 | raw file
 1#!/usr/bin/env python
 2
 3"""
 4Read a maf and output a multiple block fasta file.
 5"""
 6#Dan Blankenberg
 7import sys
 8from galaxy import eggs
 9import pkg_resources; pkg_resources.require( "bx-python" )
10from bx.align import maf
11from galaxy.tools.util import maf_utilities
12
13assert sys.version_info[:2] >= ( 2, 4 )
14
15def __main__():
16    try:
17        maf_reader = maf.Reader( open( sys.argv[1] ) )
18    except Exception, e:
19        maf_utilities.tool_fail( "Error opening input MAF: %s" % e )
20    try:
21        file_out = open( sys.argv[2], 'w' )
22    except Exception, e:
23        maf_utilities.tool_fail( "Error opening file for output: %s" % e )
24    try:
25        species = maf_utilities.parse_species_option( sys.argv[3] )
26        if species:
27            num_species = len( species )
28        else:
29            num_species = 0
30    except Exception, e:
31        maf_utilities.tool_fail( "Error determining species value: %s" % e )
32    try:
33        partial = sys.argv[4]
34    except Exception, e:
35        maf_utilities.tool_fail( "Error determining keep partial value: %s" % e )
36    
37    if species:
38        print "Restricted to species: %s" % ', '.join( species )
39    else:
40        print "Not restricted to species."
41    
42    for block_num, block in enumerate( maf_reader ):
43        if species:
44            block = block.limit_to_species( species )
45            if len( maf_utilities.get_species_in_block( block ) ) < num_species and partial == "partial_disallowed": continue
46        spec_counts = {}
47        for component in block.components:
48            spec, chrom = maf_utilities.src_split( component.src )
49            if spec not in spec_counts:
50                spec_counts[ spec ] = 0
51            else:
52                spec_counts[ spec ] += 1
53            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 ] ) ) )
54            file_out.write( "%s\n" % component.text )
55        file_out.write( "\n" )
56    file_out.close()
57
58if __name__ == "__main__": __main__()