|
0
|
1 #!/usr/bin/env python
|
|
|
2 #Original script from /home/james/work/encode/feature_partitions/split_by_partitions.py
|
|
|
3
|
|
|
4 #Usage: python(2.4) split_by_partitions.py partition_index in_file out_file chrCol startCol endCol strandCol
|
|
|
5
|
|
|
6 from __future__ import division
|
|
|
7
|
|
|
8 import sys
|
|
|
9 from galaxy import eggs
|
|
|
10 import pkg_resources; pkg_resources.require( "bx-python" )
|
|
|
11 from bx.bitset import *
|
|
|
12 from bx.bitset_builders import *
|
|
|
13
|
|
|
14 assert sys.version_info[:2] >= ( 2, 4 )
|
|
|
15
|
|
|
16 def stop_err( msg ):
|
|
|
17 sys.stderr.write( msg )
|
|
|
18 sys.exit()
|
|
|
19
|
|
|
20 def main():
|
|
|
21 GALAXY_DATA_INDEX_DIR = sys.argv[1]
|
|
|
22 partition_index = '%s/encode_feature_partitions/partition_list.txt' % GALAXY_DATA_INDEX_DIR
|
|
|
23 partition_offset = "%s/encode_feature_partitions/" % GALAXY_DATA_INDEX_DIR
|
|
|
24
|
|
|
25 warnings = []
|
|
|
26
|
|
|
27 # Load up the partitions
|
|
|
28 partitions = list()
|
|
|
29 try:
|
|
|
30 for line in open( partition_index ):
|
|
|
31 name, score, filename = line.split()
|
|
|
32 partitions.append( ( name, score, binned_bitsets_from_file( open( partition_offset+filename ) ) ) )
|
|
|
33 except:
|
|
|
34 stop_err( "Error loading partitioning dataset." )
|
|
|
35
|
|
|
36 try:
|
|
|
37 in_file = open( sys.argv[2] )
|
|
|
38 except:
|
|
|
39 stop_err( "Bad input data." )
|
|
|
40
|
|
|
41 try:
|
|
|
42 out_file = open( sys.argv[3], "w" )
|
|
|
43 except:
|
|
|
44 stop_err( "Bad output file." )
|
|
|
45
|
|
|
46 try:
|
|
|
47 chrCol = int( sys.argv[4] ) - 1
|
|
|
48 except:
|
|
|
49 stop_err( "Bad chr column: %s" % ( str( sys.argv[4] ) ) )
|
|
|
50 try:
|
|
|
51 startCol = int( sys.argv[5] ) - 1
|
|
|
52 except:
|
|
|
53 stop_err( "Bad start column: %s" % ( str( sys.argv[5] ) ) )
|
|
|
54 try:
|
|
|
55 endCol = int( sys.argv[6] ) - 1
|
|
|
56 except:
|
|
|
57 stop_err( "Bad end column: %s" % ( str( sys.argv[6] ) ) )
|
|
|
58 try:
|
|
|
59 strandCol = int( sys.argv[7] )-1
|
|
|
60 except:
|
|
|
61 strandCol = -1
|
|
|
62
|
|
|
63 line_count = 0
|
|
|
64 skipped_lines = 0
|
|
|
65 first_invalid_line = None
|
|
|
66 invalid_line = ''
|
|
|
67 try:
|
|
|
68 for line in in_file:
|
|
|
69 line_count += 1
|
|
|
70 line = line.rstrip( '\r\n' )
|
|
|
71 if line and not line.startswith( '#' ):
|
|
|
72 fields = line.split( '\t' )
|
|
|
73 try:
|
|
|
74 chr, start, end = fields[chrCol], int( fields[startCol] ), int( fields[endCol] )
|
|
|
75 except:
|
|
|
76 skipped_lines += 1
|
|
|
77 if first_invalid_line is None:
|
|
|
78 first_invalid_line = line_count
|
|
|
79 invalid_line = line
|
|
|
80 continue
|
|
|
81 label = "input_line_" + str( line_count ) #if input file type was known to be bed, then could guess at label column
|
|
|
82
|
|
|
83 if strandCol < 0:
|
|
|
84 strand = "+"
|
|
|
85 else:
|
|
|
86 try:
|
|
|
87 strand = fields[strandCol]
|
|
|
88 except:
|
|
|
89 strand = "+"
|
|
|
90
|
|
|
91 # Find which partition it overlaps
|
|
|
92 overlap = 0
|
|
|
93 for name, score, bb in partitions:
|
|
|
94 # Is there at least 1bp overlap?
|
|
|
95 if chr in bb:
|
|
|
96 overlap = bb[chr].count_range( start, end-start )
|
|
|
97 if overlap > 0:
|
|
|
98 break
|
|
|
99 else:
|
|
|
100 # No overlap with any partition? For now throw this since the
|
|
|
101 # partitions tile the encode regions completely, indicate an interval
|
|
|
102 # that does not even overlap an encode region
|
|
|
103 warning = "warning: Interval (%s, %d, %d) does not overlap any partition" % ( chr, start, end ) + ", line[" + str( line_count ) + "]. "
|
|
|
104 warnings.append( warning )
|
|
|
105 name = "no_overlap"
|
|
|
106 score = 0
|
|
|
107 # Annotate with the name of the partition
|
|
|
108 frac_overlap = overlap / ( end-start )
|
|
|
109 # BED6 plus?
|
|
|
110 print >>out_file, "%s\t%d\t%d\t%s\t%s\t%s\t%s\t%0.4f" % ( chr, start, end, label, score, strand, name, frac_overlap )
|
|
|
111 except:
|
|
|
112 out_file.close()
|
|
|
113 in_file.close()
|
|
|
114 stop_err( "Unknown error while processing line # %d: %s" % ( line_count, line ) )
|
|
|
115 out_file.close()
|
|
|
116 in_file.close()
|
|
|
117
|
|
|
118 if warnings:
|
|
|
119 warn_msg = "This tool is useful on ENCODE regions only, %d warnings, 1st is: " % len( warnings )
|
|
|
120 warn_msg += warnings[0]
|
|
|
121 print warn_msg
|
|
|
122 if skipped_lines:
|
|
|
123 print "Skipped %d invalid lines starting at line # %d: %s" % ( skipped_lines, first_invalid_line, invalid_line )
|
|
|
124
|
|
|
125 if __name__ == "__main__": main()
|
