Mercurial > repos > xuebing > sharplabtool

comparison tools/encode/split_by_partitions.py @ 0:9071e359b9a3

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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 #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()