root/galaxy-central/tools/encode/split_by_partitions.py

リビジョン 2, 4.5 KB (コミッタ: hatakeyama, 14 年 前)

import galaxy-central

  • 属性 svn:executable の設定値 *
行番号 
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
6from __future__ import division
7
8import sys
9from galaxy import eggs
10import pkg_resources; pkg_resources.require( "bx-python" )
11from bx.bitset import *
12from bx.bitset_builders import *
13
14assert sys.version_info[:2] >= ( 2, 4 )
15
16def stop_err( msg ):
17    sys.stderr.write( msg )
18    sys.exit()
19
20def 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
125if __name__ == "__main__": main()
Note: リポジトリブラウザについてのヘルプは TracBrowser を参照してください。