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

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

import galaxy-central

行番号 
1#!/usr/bin/env python
2#Dan Blankenberg
3#%prog bounding_region_file mask_intervals_file intervals_to_mimic_file out_file mask_chr mask_start mask_end interval_chr interval_start interval_end interval_strand use_mask allow_strand_overlaps
4import sys, random
5from copy import deepcopy
6from galaxy import eggs
7import pkg_resources
8pkg_resources.require( "bx-python" )
9import bx.intervals.io
10import bx.intervals.intersection
11import psyco_full
12
13assert sys.version_info[:2] >= ( 2, 4 )
14
15max_iters = 5
16
17def stop_err( msg ):
18    sys.stderr.write( msg )
19    sys.exit()
20
21#Try to add a random region
22def add_random_region( mimic_region, bound, exist_regions, plus_mask, minus_mask, overlaps ):
23    region_length, region_strand = mimic_region
24    plus_count = plus_mask.count_range()
25    minus_count = minus_mask.count_range()
26    gaps = []
27
28    if region_strand == "-":
29        gaps = minus_mask.get_gaps( region_length )
30    else:
31        gaps = plus_mask.get_gaps( region_length )
32   
33    while True:
34        try:
35            gap_length, gap_start, gap_end = gaps.pop( random.randint( 0, len( gaps ) - 1 ) )
36        except:
37            break
38        try:
39            start = random.randint( bound.start + gap_start, bound.start + gap_end - region_length - 1 )
40        except ValueError, ve:
41            stop_err( "Exception thrown generating random start value: %s" %str( ve ) )
42
43        end = start + region_length
44        try_plus_mask = plus_mask.copy()
45        try_minus_mask = minus_mask.copy()
46       
47        if region_strand == "-":
48            try_minus_mask.set_range( start - bound.start, end - bound.start )
49        else:
50            try_plus_mask.set_range( start - bound.start, end - bound.start )
51       
52        rand_region = bx.intervals.io.GenomicInterval( None, [bound.chrom, start, end, region_strand], 0, 1, 2, 3, "+", fix_strand=True )
53       
54        if try_plus_mask.count_range() == plus_count + region_length or try_minus_mask.count_range() == minus_count + region_length:
55            if overlaps in ["strand", "all"]: #overlaps allowed across strands
56                exist_regions.append( rand_region )
57                if overlaps == "strand":
58                    return exist_regions, True, try_plus_mask, try_minus_mask
59                else: #overlaps allowed everywhere
60                    return exist_regions, True, plus_mask, minus_mask
61            else: #no overlapping anywhere
62                exist_regions.append( rand_region )
63                if region_strand == "-":
64                    return exist_regions, True, try_minus_mask.copy(), try_minus_mask
65                else:
66                    return exist_regions, True, try_plus_mask, try_plus_mask.copy()
67    return exist_regions, False, plus_mask, minus_mask
68
69def main():
70    includes_strand = False
71    region_uid = sys.argv[1]
72    mask_fname = sys.argv[2]
73    intervals_fname = sys.argv[3]
74    out_fname = sys.argv[4]
75    try:
76        mask_chr = int( sys.argv[5] ) - 1
77    except:
78        stop_err( "'%s' is an invalid chrom column for 'Intervals to Mask' dataset, click the pencil icon in the history item to edit column settings." % str( sys.argv[5] ) )
79    try:
80        mask_start = int( sys.argv[6] ) - 1
81    except:
82        stop_err( "'%s' is an invalid start column for 'Intervals to Mask' dataset, click the pencil icon in the history item to edit column settings." % str( sys.argv[6] ) )
83    try:
84        mask_end = int( sys.argv[7] ) - 1
85    except:
86        stop_err( "'%s' is an invalid end column for 'Intervals to Mask' dataset, click the pencil icon in the history item to edit column settings." % str( sys.argv[7] ) )
87    try:
88        interval_chr = int( sys.argv[8] ) - 1
89    except:
90        stop_err( "'%s' is an invalid chrom column for 'File to Mimick' dataset, click the pencil icon in the history item to edit column settings." % str( sys.argv[8] ) )
91    try:
92        interval_start = int( sys.argv[9] ) - 1
93    except:
94        stop_err( "'%s' is an invalid start column for 'File to Mimick' dataset, click the pencil icon in the history item to edit column settings." % str( sys.argv[9] ) )
95    try:
96        interval_end = int( sys.argv[10] ) - 1
97    except:
98        stop_err( "'%s' is an invalid end column for 'File to Mimick' dataset, click the pencil icon in the history item to edit column settings." % str( sys.argv[10] ) )
99    try:
100        interval_strand = int( sys.argv[11] ) - 1
101        includes_strand = True
102    except:
103        interval_strand = -1
104    if includes_strand:
105        use_mask = sys.argv[12]
106        overlaps = sys.argv[13]
107    else:
108        use_mask = sys.argv[11]
109        overlaps = sys.argv[12]
110    available_regions = {}
111    loc_file = "%s/regions.loc" % sys.argv[-1]
112   
113    for i, line in enumerate( file( loc_file ) ):
114        line = line.rstrip( '\r\n' )
115        if line and not line.startswith( '#' ):
116            fields = line.split( '\t' )
117            #read each line, if not enough fields, go to next line
118            try:
119                build = fields[0]
120                uid = fields[1]
121                description = fields[2]
122                filepath = fields[3]
123                available_regions[uid] = filepath
124            except:
125                continue
126
127    if region_uid not in available_regions:
128        stop_err( "Region '%s' is invalid." % region_uid )
129    region_fname = available_regions[region_uid].strip()
130
131    #set up bounding regions to hold random intervals
132    bounds = []
133    for bound in bx.intervals.io.NiceReaderWrapper( open( region_fname, 'r' ), chrom_col=0, start_col=1, end_col=2, fix_strand=True, return_header=False, return_comments=False ):
134        bounds.append( bound )
135    #set up length and number of regions to mimic
136    regions = [ [] for i in range( len( bounds ) ) ]
137
138    for region in bx.intervals.io.NiceReaderWrapper( open( intervals_fname, 'r' ), chrom_col=interval_chr, start_col=interval_start, end_col=interval_end, strand_col=interval_strand, fix_strand=True, return_header=False, return_comments=False ):
139        #loop through bounds, find first proper bounds then add
140        #if an interval crosses bounds, it will be added to the first bound
141        for i in range( len( bounds ) ):
142            if bounds[i].chrom != region.chrom:
143                continue
144            intersecter = bx.intervals.intersection.Intersecter()
145            intersecter.add_interval( bounds[i] )
146            if len( intersecter.find( region.start, region.end ) ) > 0:
147                regions[i].append( ( region.end - region.start, region.strand ) ) #add region to proper bound and go to next region
148                break
149    for region in regions:
150        region.sort()
151        region.reverse()
152   
153    #read mask file
154    mask = []
155    if use_mask != "no_mask":
156        for region in bx.intervals.io.NiceReaderWrapper( open( mask_fname, 'r' ), chrom_col=mask_chr, start_col=mask_start, end_col=mask_end, fix_strand=True, return_header=False, return_comments=False ):
157            mask.append( region )
158
159    try:
160        out_file = open ( out_fname, "w" )
161    except:
162        stop_err( "Error opening output file '%s'." % out_fname )
163
164    i = 0
165    i_iters = 0
166    region_count = 0
167    best_regions = []
168    num_fail = 0
169    while i < len( bounds ):
170        i_iters += 1
171        #order regions to mimic
172        regions_to_mimic = regions[i][0:]
173        if len( regions_to_mimic ) < 1: #if no regions to mimic, skip
174            i += 1
175            i_iters = 0
176            continue
177        #set up region mask
178        plus_mask = Region( bounds[i].end - bounds[i].start )
179        for region in mask:
180            if region.chrom != bounds[i].chrom: continue
181            mask_start = region.start - bounds[i].start
182            mask_end = region.end - bounds[i].start
183            if mask_start >= 0 and mask_end > 0:
184                plus_mask.set_range( mask_start, mask_end )
185        minus_mask = plus_mask.copy()
186        random_regions = []
187        num_added = 0
188        for j in range( len( regions[i] ) ):
189            random_regions, added, plus_mask, minus_mask = add_random_region( regions_to_mimic[j], bounds[i], random_regions, plus_mask, minus_mask, overlaps )
190            if added:
191                num_added += 1
192        if num_added == len( regions_to_mimic ) or i_iters >= max_iters:
193            if len( best_regions ) > len( random_regions ):
194                random_regions = best_regions.copy()
195            num_fail += ( len( regions_to_mimic ) - len( random_regions ) )
196            i_iters = 0
197            best_regions = []
198            for region in random_regions:
199                print >>out_file, "%s\t%d\t%d\t%s\t%s\t%s" % ( region.chrom, region.start, region.end, "region_" + str( region_count ), "0", region.strand )
200                region_count += 1
201        else:
202            i -= 1
203            if len( best_regions ) < len( random_regions ):
204                best_regions = random_regions[:]
205        i+=1
206   
207    out_file.close()
208    if num_fail:
209        print "After %i iterations, %i regions could not be added." % (max_iters, num_fail)
210        if use_mask == "use_mask":
211            print "The mask you have provided may be too restrictive."
212
213class Region( list ):
214    """
215    A list for on/off regions
216    """
217    def __init__( self, size=0 ):
218        for i in range( size ):
219            self.append( False )
220    def copy( self ):
221        return deepcopy( self )
222    def set_range( self, start=0, end=None ):
223        if start < 0:
224            start = 0
225        if ( not end and end != 0 ) or end > len( self ):
226            end = len( self )
227        for i in range( start, end ):
228            self[i]=True
229    def count_range( self, start=0, end=None ):
230        if start < 0:
231            start = 0
232        if ( not end and end != 0 ) or end > len( self ):
233            end = len( self )
234        return self[start:end].count( True )
235    def get_gaps( self, min_size = 0 ):
236        gaps = []
237        start = end = 0
238        while True:
239            try:
240                start = self[end:].index( False ) + end
241            except:
242                break
243            try:
244                end = self[start:].index( True ) + start
245            except:
246                end = len( self )
247            if end > start and end - start >= min_size:
248                gaps.append( ( end - start, start, end ) )
249        gaps.sort()
250        gaps.reverse()
251        return gaps
252
253if __name__ == "__main__": main()
Note: リポジトリブラウザについてのヘルプは TracBrowser を参照してください。