Mercurial > repos > iuc > maf_stats

comparison maf_stats.py @ 0:5a92f4c476b5 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/interval2maf/ commit f7b04f8c6edca90b1e9bc867a1569996aec76d7a"
author iuc
date Thu, 13 Aug 2020 09:30:47 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:5a92f4c476b5
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 from __future__ import print_function
7
8 import sys
9
10 import bx.intervals.io
11 from bx.bitset import BitSet
12 from galaxy.tools.util import maf_utilities
13
14
15 def __main__():
16 maf_source_type = sys.argv.pop(1)
17 input_maf_filename = sys.argv[1].strip()
18 input_interval_filename = sys.argv[2].strip()
19 output_filename = sys.argv[3].strip()
20 dbkey = sys.argv[4].strip()
21 try:
22 chr_col = int(sys.argv[5].strip()) - 1
23 start_col = int(sys.argv[6].strip()) - 1
24 end_col = int(sys.argv[7].strip()) - 1
25 except Exception:
26 print("You appear to be missing metadata. You can specify your metadata by clicking on the pencil icon associated with your interval file.", file=sys.stderr)
27 sys.exit()
28 summary = sys.argv[8].strip()
29 if summary.lower() == "true":
30 summary = True
31 else:
32 summary = False
33
34 mafIndexFile = "%s/maf_index.loc" % sys.argv[9]
35 try:
36 maf_index_filename = sys.argv[10].strip()
37 except Exception:
38 maf_index_filename = None
39 index = index_filename = None
40 if maf_source_type == "user":
41 # index maf for use here
42 index, index_filename = maf_utilities.open_or_build_maf_index(input_maf_filename, maf_index_filename, species=[dbkey])
43 if index is None:
44 print("Your MAF file appears to be malformed.", file=sys.stderr)
45 sys.exit()
46 elif maf_source_type == "cached":
47 # access existing indexes
48 index = maf_utilities.maf_index_by_uid(input_maf_filename, mafIndexFile)
49 if index is None:
50 print("The MAF source specified (%s) appears to be invalid." % (input_maf_filename), file=sys.stderr)
51 sys.exit()
52 else:
53 print('Invalid source type specified: %s' % maf_source_type, file=sys.stdout)
54 sys.exit()
55
56 out = open(output_filename, 'w')
57
58 num_region = None
59 num_bad_region = 0
60 species_summary = {}
61 total_length = 0
62 # loop through interval file
63 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)):
64 src = "%s.%s" % (dbkey, region.chrom)
65 region_length = region.end - region.start
66 if region_length < 1:
67 num_bad_region += 1
68 continue
69 total_length += region_length
70 coverage = {dbkey: BitSet(region_length)}
71
72 for block in index.get_as_iterator(src, region.start, region.end):
73 for spec in maf_utilities.get_species_in_block(block):
74 if spec not in coverage:
75 coverage[spec] = BitSet(region_length)
76 for block in maf_utilities.iter_blocks_split_by_species(block):
77 if maf_utilities.component_overlaps_region(block.get_component_by_src(src), region):
78 # need to chop and orient the block
79 block = maf_utilities.orient_block_by_region(maf_utilities.chop_block_by_region(block, src, region), src, region, force_strand='+')
80 start_offset, alignment = maf_utilities.reduce_block_by_primary_genome(block, dbkey, region.chrom, region.start)
81 for i in range(len(alignment[dbkey])):
82 for spec, text in alignment.items():
83 if text[i] != '-':
84 coverage[spec].set(start_offset + i)
85 if summary:
86 # record summary
87 for key in coverage.keys():
88 if key not in species_summary:
89 species_summary[key] = 0
90 species_summary[key] = species_summary[key] + coverage[key].count_range()
91 else:
92 # print coverage for interval
93 coverage_sum = coverage[dbkey].count_range()
94 out.write("%s\t%s\t%s\t%s\n" % ("\t".join(region.fields), dbkey, coverage_sum, region_length - coverage_sum))
95 keys = list(coverage.keys())
96 keys.remove(dbkey)
97 keys.sort()
98 for key in keys:
99 coverage_sum = coverage[key].count_range()
100 out.write("%s\t%s\t%s\t%s\n" % ("\t".join(region.fields), key, coverage_sum, region_length - coverage_sum))
101 if summary:
102 out.write("#species\tnucleotides\tcoverage\n")
103 for spec in species_summary:
104 out.write("%s\t%s\t%.4f\n" % (spec, species_summary[spec], float(species_summary[spec]) / total_length))
105 out.close()
106 if num_region is not None:
107 print("%i regions were processed with a total length of %i." % (num_region + 1, total_length))
108 if num_bad_region:
109 print("%i regions were invalid." % (num_bad_region))
110 maf_utilities.remove_temp_index_file(index_filename)
111
112
113 if __name__ == "__main__":
114 __main__()