|
0
|
1 #!/usr/bin/env python
|
|
|
2
|
|
|
3 """
|
|
|
4 Read a maf file and write out a new maf with only blocks having the
|
|
|
5 required species, after dropping any other species and removing
|
|
|
6 columns containing only gaps.
|
|
|
7
|
|
|
8 usage: %prog species,species2,... input_maf output_maf allow_partial min_species_per_block
|
|
|
9 """
|
|
|
10 #Dan Blankenberg
|
|
|
11 from galaxy import eggs
|
|
|
12 import pkg_resources; pkg_resources.require( "bx-python" )
|
|
|
13 import bx.align.maf
|
|
|
14 from galaxy.tools.util import maf_utilities
|
|
|
15 import sys
|
|
|
16
|
|
|
17 assert sys.version_info[:2] >= ( 2, 4 )
|
|
|
18
|
|
|
19 def main():
|
|
|
20
|
|
|
21 species = maf_utilities.parse_species_option( sys.argv[1] )
|
|
|
22 if species:
|
|
|
23 spec_len = len( species )
|
|
|
24 else:
|
|
|
25 spec_len = 0
|
|
|
26 try:
|
|
|
27 maf_reader = bx.align.maf.Reader( open( sys.argv[2],'r' ) )
|
|
|
28 maf_writer = bx.align.maf.Writer( open( sys.argv[3],'w' ) )
|
|
|
29 except:
|
|
|
30 print >>sys.stderr, "Your MAF file appears to be malformed."
|
|
|
31 sys.exit()
|
|
|
32 allow_partial = False
|
|
|
33 if int( sys.argv[4] ): allow_partial = True
|
|
|
34 min_species_per_block = int( sys.argv[5] )
|
|
|
35
|
|
|
36 maf_blocks_kept = 0
|
|
|
37 for m in maf_reader:
|
|
|
38 if species:
|
|
|
39 m = m.limit_to_species( species )
|
|
|
40 m.remove_all_gap_columns()
|
|
|
41 spec_in_block_len = len( maf_utilities.get_species_in_block( m ) )
|
|
|
42 if ( not species or allow_partial or spec_in_block_len == spec_len ) and spec_in_block_len > min_species_per_block:
|
|
|
43 maf_writer.write( m )
|
|
|
44 maf_blocks_kept += 1
|
|
|
45
|
|
|
46 maf_reader.close()
|
|
|
47 maf_writer.close()
|
|
|
48
|
|
|
49 print "Restricted to species: %s." % ", ".join( species )
|
|
|
50 print "%i MAF blocks have been kept." % maf_blocks_kept
|
|
|
51
|
|
|
52 if __name__ == "__main__":
|
|
|
53 main()
|
