[2] | 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() |
---|