Mercurial > repos > xuebing > sharplabtool
comparison tools/maf/maf_stats.py @ 0:9071e359b9a3
Uploaded
author | xuebing |
---|---|
date | Fri, 09 Mar 2012 19:37:19 -0500 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:9071e359b9a3 |
---|---|
1 #!/usr/bin/env python | |
2 #Dan Blankenberg | |
3 """ | |
4 Reads a list of intervals and a maf. Outputs a new set of intervals with statistics appended. | |
5 """ | |
6 | |
7 import sys | |
8 from galaxy import eggs | |
9 import pkg_resources; pkg_resources.require( "bx-python" ) | |
10 import bx.intervals.io | |
11 from bx.bitset import BitSet | |
12 from galaxy.tools.util import maf_utilities | |
13 | |
14 assert sys.version_info[:2] >= ( 2, 4 ) | |
15 | |
16 def __main__(): | |
17 maf_source_type = sys.argv.pop( 1 ) | |
18 input_maf_filename = sys.argv[1].strip() | |
19 input_interval_filename = sys.argv[2].strip() | |
20 output_filename = sys.argv[3].strip() | |
21 dbkey = sys.argv[4].strip() | |
22 try: | |
23 chr_col = int( sys.argv[5].strip() ) - 1 | |
24 start_col = int( sys.argv[6].strip() ) - 1 | |
25 end_col = int( sys.argv[7].strip() ) - 1 | |
26 except: | |
27 print >>sys.stderr, "You appear to be missing metadata. You can specify your metadata by clicking on the pencil icon associated with your interval file." | |
28 sys.exit() | |
29 summary = sys.argv[8].strip() | |
30 if summary.lower() == "true": summary = True | |
31 else: summary = False | |
32 | |
33 mafIndexFile = "%s/maf_index.loc" % sys.argv[9] | |
34 try: | |
35 maf_index_filename = sys.argv[10].strip() | |
36 except: | |
37 maf_index_filename = None | |
38 index = index_filename = None | |
39 if maf_source_type == "user": | |
40 #index maf for use here | |
41 index, index_filename = maf_utilities.open_or_build_maf_index( input_maf_filename, maf_index_filename, species = [dbkey] ) | |
42 if index is None: | |
43 print >>sys.stderr, "Your MAF file appears to be malformed." | |
44 sys.exit() | |
45 elif maf_source_type == "cached": | |
46 #access existing indexes | |
47 index = maf_utilities.maf_index_by_uid( input_maf_filename, mafIndexFile ) | |
48 if index is None: | |
49 print >> sys.stderr, "The MAF source specified (%s) appears to be invalid." % ( input_maf_filename ) | |
50 sys.exit() | |
51 else: | |
52 print >>sys.stdout, 'Invalid source type specified: %s' % maf_source_type | |
53 sys.exit() | |
54 | |
55 out = open(output_filename, 'w') | |
56 | |
57 num_region = None | |
58 species_summary = {} | |
59 total_length = 0 | |
60 #loop through interval file | |
61 for num_region, region in enumerate( bx.intervals.io.NiceReaderWrapper( open( input_interval_filename, 'r' ), chrom_col = chr_col, start_col = start_col, end_col = end_col, fix_strand = True, return_header = False, return_comments = False ) ): | |
62 src = "%s.%s" % ( dbkey, region.chrom ) | |
63 region_length = region.end - region.start | |
64 total_length += region_length | |
65 coverage = { dbkey: BitSet( region_length ) } | |
66 | |
67 | |
68 for block in index.get_as_iterator( src, region.start, region.end ): | |
69 for spec in maf_utilities.get_species_in_block( block ): | |
70 if spec not in coverage: coverage[spec] = BitSet( region_length ) | |
71 for block in maf_utilities.iter_blocks_split_by_species( block ): | |
72 if maf_utilities.component_overlaps_region( block.get_component_by_src( src ), region ): | |
73 #need to chop and orient the block | |
74 block = maf_utilities.orient_block_by_region( maf_utilities.chop_block_by_region( block, src, region ), src, region, force_strand = '+' ) | |
75 start_offset, alignment = maf_utilities.reduce_block_by_primary_genome( block, dbkey, region.chrom, region.start ) | |
76 for i in range( len( alignment[dbkey] ) ): | |
77 for spec, text in alignment.items(): | |
78 if text[i] != '-': | |
79 coverage[spec].set( start_offset + i ) | |
80 if summary: | |
81 #record summary | |
82 for key in coverage.keys(): | |
83 if key not in species_summary: species_summary[key] = 0 | |
84 species_summary[key] = species_summary[key] + coverage[key].count_range() | |
85 else: | |
86 #print coverage for interval | |
87 coverage_sum = coverage[dbkey].count_range() | |
88 out.write( "%s\t%s\t%s\t%s\n" % ( "\t".join( region.fields ), dbkey, coverage_sum, region_length - coverage_sum ) ) | |
89 keys = coverage.keys() | |
90 keys.remove( dbkey ) | |
91 keys.sort() | |
92 for key in keys: | |
93 coverage_sum = coverage[key].count_range() | |
94 out.write( "%s\t%s\t%s\t%s\n" % ( "\t".join( region.fields ), key, coverage_sum, region_length - coverage_sum ) ) | |
95 if summary: | |
96 out.write( "#species\tnucleotides\tcoverage\n" ) | |
97 for spec in species_summary: | |
98 out.write( "%s\t%s\t%.4f\n" % ( spec, species_summary[spec], float( species_summary[spec] ) / total_length ) ) | |
99 out.close() | |
100 if num_region is not None: | |
101 print "%i regions were processed with a total length of %i." % ( num_region + 1, total_length ) | |
102 maf_utilities.remove_temp_index_file( index_filename ) | |
103 | |
104 if __name__ == "__main__": __main__() |