Mercurial > repos > devteam > flanking_features

annotate flanking_features.py @ 1:850c05b9af00 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload commit 33927a87ba2eee9bf0ecdd376a66241b17b3d734
author devteam
date Tue, 13 Oct 2015 12:50:14 -0400
parents e928e029f6eb
children 94248d5b9b8b
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
1 #!/usr/bin/env python
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
2 #By: Guruprasad Ananda
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
3 """
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
4 Fetch closest up/downstream interval from features corresponding to every interval in primary
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
5
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
6 usage: %prog primary_file features_file out_file direction
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
7 -1, --cols1=N,N,N,N: Columns for start, end, strand in first file
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
8 -2, --cols2=N,N,N,N: Columns for start, end, strand in second file
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
9 -G, --gff1: input 1 is GFF format, meaning start and end coordinates are 1-based, closed interval
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
10 -H, --gff2: input 2 is GFF format, meaning start and end coordinates are 1-based, closed interval
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
11 """
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
12
1
850c05b9af00 planemo upload commit 33927a87ba2eee9bf0ecdd376a66241b17b3d734
devteam
parents: 0
diff changeset
13 import fileinput
850c05b9af00 planemo upload commit 33927a87ba2eee9bf0ecdd376a66241b17b3d734
devteam
parents: 0
diff changeset
14 import sys
0
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
15 from bx.cookbook import doc_optparse
1
850c05b9af00 planemo upload commit 33927a87ba2eee9bf0ecdd376a66241b17b3d734
devteam
parents: 0
diff changeset
16 from bx.intervals.io import Comment, GenomicInterval, Header, NiceReaderWrapper
0
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
17 from bx.intervals.operations import quicksect
1
850c05b9af00 planemo upload commit 33927a87ba2eee9bf0ecdd376a66241b17b3d734
devteam
parents: 0
diff changeset
18 from bx.tabular.io import ParseError
850c05b9af00 planemo upload commit 33927a87ba2eee9bf0ecdd376a66241b17b3d734
devteam
parents: 0
diff changeset
19 from galaxy.tools.util.galaxyops import fail, parse_cols_arg, skipped
850c05b9af00 planemo upload commit 33927a87ba2eee9bf0ecdd376a66241b17b3d734
devteam
parents: 0
diff changeset
20 from utils.gff_util import convert_bed_coords_to_gff, GFFIntervalToBEDReaderWrapper
0
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
21
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
22 assert sys.version_info[:2] >= ( 2, 4 )
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
23
1
850c05b9af00 planemo upload commit 33927a87ba2eee9bf0ecdd376a66241b17b3d734
devteam
parents: 0
diff changeset
24
850c05b9af00 planemo upload commit 33927a87ba2eee9bf0ecdd376a66241b17b3d734
devteam
parents: 0
diff changeset
25 def get_closest_feature(node, direction, threshold_up, threshold_down, report_func_up, report_func_down):
0
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
26 #direction=1 for +ve strand upstream and -ve strand downstream cases; and it is 0 for +ve strand downstream and -ve strand upstream cases
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
27 #threhold_Up is equal to the interval start for +ve strand, and interval end for -ve strand
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
28 #threhold_down is equal to the interval end for +ve strand, and interval start for -ve strand
1
850c05b9af00 planemo upload commit 33927a87ba2eee9bf0ecdd376a66241b17b3d734
devteam
parents: 0
diff changeset
29 if direction == 1:
0
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
30 if node.maxend <= threshold_up:
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
31 if node.end == node.maxend:
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
32 report_func_up(node)
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
33 elif node.right and node.left:
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
34 if node.right.maxend == node.maxend:
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
35 get_closest_feature(node.right, direction, threshold_up, threshold_down, report_func_up, report_func_down)
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
36 elif node.left.maxend == node.maxend:
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
37 get_closest_feature(node.left, direction, threshold_up, threshold_down, report_func_up, report_func_down)
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
38 elif node.right and node.right.maxend == node.maxend:
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
39 get_closest_feature(node.right, direction, threshold_up, threshold_down, report_func_up, report_func_down)
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
40 elif node.left and node.left.maxend == node.maxend:
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
41 get_closest_feature(node.left, direction, threshold_up, threshold_down, report_func_up, report_func_down)
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
42 elif node.minend <= threshold_up:
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
43 if node.end <= threshold_up:
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
44 report_func_up(node)
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
45 if node.left and node.right:
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
46 if node.right.minend <= threshold_up:
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
47 get_closest_feature(node.right, direction, threshold_up, threshold_down, report_func_up, report_func_down)
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
48 if node.left.minend <= threshold_up:
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
49 get_closest_feature(node.left, direction, threshold_up, threshold_down, report_func_up, report_func_down)
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
50 elif node.left:
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
51 if node.left.minend <= threshold_up:
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
52 get_closest_feature(node.left, direction, threshold_up, threshold_down, report_func_up, report_func_down)
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
53 elif node.right:
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
54 if node.right.minend <= threshold_up:
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
55 get_closest_feature(node.right, direction, threshold_up, threshold_down, report_func_up, report_func_down)
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
56 elif direction == 0:
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
57 if node.start > threshold_down:
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
58 report_func_down(node)
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
59 if node.left:
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
60 get_closest_feature(node.left, direction, threshold_up, threshold_down, report_func_up, report_func_down)
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
61 else:
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
62 if node.right:
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
63 get_closest_feature(node.right, direction, threshold_up, threshold_down, report_func_up, report_func_down)
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
64
1
850c05b9af00 planemo upload commit 33927a87ba2eee9bf0ecdd376a66241b17b3d734
devteam
parents: 0
diff changeset
65
0
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
66 def proximal_region_finder(readers, region, comments=True):
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
67 """
1
850c05b9af00 planemo upload commit 33927a87ba2eee9bf0ecdd376a66241b17b3d734
devteam
parents: 0
diff changeset
68 Returns an iterator that yields elements of the form [ <original_interval>, <closest_feature> ].
850c05b9af00 planemo upload commit 33927a87ba2eee9bf0ecdd376a66241b17b3d734
devteam
parents: 0
diff changeset
69 Intervals are GenomicInterval objects.
0
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
70 """
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
71 primary = readers[0]
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
72 features = readers[1]
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
73 either = False
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
74 if region == 'Upstream':
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
75 up, down = True, False
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
76 elif region == 'Downstream':
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
77 up, down = False, True
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
78 else:
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
79 up, down = True, True
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
80 if region == 'Either':
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
81 either = True
1
850c05b9af00 planemo upload commit 33927a87ba2eee9bf0ecdd376a66241b17b3d734
devteam
parents: 0
diff changeset
82
0
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
83 # Read features into memory:
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
84 rightTree = quicksect.IntervalTree()
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
85 for item in features:
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
86 if type( item ) is GenomicInterval:
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
87 rightTree.insert( item, features.linenum, item )
1
850c05b9af00 planemo upload commit 33927a87ba2eee9bf0ecdd376a66241b17b3d734
devteam
parents: 0
diff changeset
88
0
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
89 for interval in primary:
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
90 if type( interval ) is Header:
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
91 yield interval
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
92 if type( interval ) is Comment and comments:
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
93 yield interval
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
94 elif type( interval ) == GenomicInterval:
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
95 chrom = interval.chrom
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
96 start = int(interval.start)
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
97 end = int(interval.end)
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
98 strand = interval.strand
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
99 if chrom not in rightTree.chroms:
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
100 continue
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
101 else:
1
850c05b9af00 planemo upload commit 33927a87ba2eee9bf0ecdd376a66241b17b3d734
devteam
parents: 0
diff changeset
102 root = rightTree.chroms[chrom] # root node for the chrom tree
0
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
103 result_up = []
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
104 result_down = []
1
850c05b9af00 planemo upload commit 33927a87ba2eee9bf0ecdd376a66241b17b3d734
devteam
parents: 0
diff changeset
105 if (strand == '+' and up) or (strand == '-' and down):
0
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
106 #upstream +ve strand and downstream -ve strand cases
1
850c05b9af00 planemo upload commit 33927a87ba2eee9bf0ecdd376a66241b17b3d734
devteam
parents: 0
diff changeset
107 get_closest_feature(root, 1, start, None, lambda node: result_up.append( node ), None)
850c05b9af00 planemo upload commit 33927a87ba2eee9bf0ecdd376a66241b17b3d734
devteam
parents: 0
diff changeset
108
0
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
109 if (strand == '+' and down) or (strand == '-' and up):
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
110 #downstream +ve strand and upstream -ve strand case
1
850c05b9af00 planemo upload commit 33927a87ba2eee9bf0ecdd376a66241b17b3d734
devteam
parents: 0
diff changeset
111 get_closest_feature(root, 0, None, end - 1, None, lambda node: result_down.append( node ))
850c05b9af00 planemo upload commit 33927a87ba2eee9bf0ecdd376a66241b17b3d734
devteam
parents: 0
diff changeset
112
0
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
113 if result_up:
1
850c05b9af00 planemo upload commit 33927a87ba2eee9bf0ecdd376a66241b17b3d734
devteam
parents: 0
diff changeset
114 if len(result_up) > 1: # The results_up list has a list of intervals upstream to the given interval.
0
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
115 ends = []
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
116 for n in result_up:
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
117 ends.append(n.end)
1
850c05b9af00 planemo upload commit 33927a87ba2eee9bf0ecdd376a66241b17b3d734
devteam
parents: 0
diff changeset
118 res_ind = ends.index(max(ends)) # fetch the index of the closest interval i.e. the interval with the max end from the results_up list
0
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
119 else:
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
120 res_ind = 0
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
121 if not(either):
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
122 yield [ interval, result_up[res_ind].other ]
1
850c05b9af00 planemo upload commit 33927a87ba2eee9bf0ecdd376a66241b17b3d734
devteam
parents: 0
diff changeset
123
850c05b9af00 planemo upload commit 33927a87ba2eee9bf0ecdd376a66241b17b3d734
devteam
parents: 0
diff changeset
124 if result_down:
0
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
125 if not(either):
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
126 #The last element of result_down will be the closest element to the given interval
1
850c05b9af00 planemo upload commit 33927a87ba2eee9bf0ecdd376a66241b17b3d734
devteam
parents: 0
diff changeset
127 yield [ interval, result_down[-1].other ]
850c05b9af00 planemo upload commit 33927a87ba2eee9bf0ecdd376a66241b17b3d734
devteam
parents: 0
diff changeset
128
0
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
129 if either and (result_up or result_down):
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
130 iter_val = []
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
131 if result_up and result_down:
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
132 if abs(start - int(result_up[res_ind].end)) <= abs(end - int(result_down[-1].start)):
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
133 iter_val = [ interval, result_up[res_ind].other ]
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
134 else:
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
135 #The last element of result_down will be the closest element to the given interval
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
136 iter_val = [ interval, result_down[-1].other ]
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
137 elif result_up:
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
138 iter_val = [ interval, result_up[res_ind].other ]
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
139 elif result_down:
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
140 #The last element of result_down will be the closest element to the given interval
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
141 iter_val = [ interval, result_down[-1].other ]
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
142 yield iter_val
1
850c05b9af00 planemo upload commit 33927a87ba2eee9bf0ecdd376a66241b17b3d734
devteam
parents: 0
diff changeset
143
850c05b9af00 planemo upload commit 33927a87ba2eee9bf0ecdd376a66241b17b3d734
devteam
parents: 0
diff changeset
144
0
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
145 def main():
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
146 options, args = doc_optparse.parse( __doc__ )
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
147 try:
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
148 chr_col_1, start_col_1, end_col_1, strand_col_1 = parse_cols_arg( options.cols1 )
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
149 chr_col_2, start_col_2, end_col_2, strand_col_2 = parse_cols_arg( options.cols2 )
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
150 in1_gff_format = bool( options.gff1 )
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
151 in2_gff_format = bool( options.gff2 )
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
152 in_fname, in2_fname, out_fname, direction = args
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
153 except:
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
154 doc_optparse.exception()
1
850c05b9af00 planemo upload commit 33927a87ba2eee9bf0ecdd376a66241b17b3d734
devteam
parents: 0
diff changeset
155
0
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
156 # Set readers to handle either GFF or default format.
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
157 if in1_gff_format:
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
158 in1_reader_wrapper = GFFIntervalToBEDReaderWrapper
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
159 else:
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
160 in1_reader_wrapper = NiceReaderWrapper
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
161 if in2_gff_format:
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
162 in2_reader_wrapper = GFFIntervalToBEDReaderWrapper
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
163 else:
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
164 in2_reader_wrapper = NiceReaderWrapper
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
165
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
166 g1 = in1_reader_wrapper( fileinput.FileInput( in_fname ),
1
850c05b9af00 planemo upload commit 33927a87ba2eee9bf0ecdd376a66241b17b3d734
devteam
parents: 0
diff changeset
167 chrom_col=chr_col_1,
850c05b9af00 planemo upload commit 33927a87ba2eee9bf0ecdd376a66241b17b3d734
devteam
parents: 0
diff changeset
168 start_col=start_col_1,
850c05b9af00 planemo upload commit 33927a87ba2eee9bf0ecdd376a66241b17b3d734
devteam
parents: 0
diff changeset
169 end_col=end_col_1,
850c05b9af00 planemo upload commit 33927a87ba2eee9bf0ecdd376a66241b17b3d734
devteam
parents: 0
diff changeset
170 strand_col=strand_col_1,
850c05b9af00 planemo upload commit 33927a87ba2eee9bf0ecdd376a66241b17b3d734
devteam
parents: 0
diff changeset
171 fix_strand=True )
0
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
172 g2 = in2_reader_wrapper( fileinput.FileInput( in2_fname ),
1
850c05b9af00 planemo upload commit 33927a87ba2eee9bf0ecdd376a66241b17b3d734
devteam
parents: 0
diff changeset
173 chrom_col=chr_col_2,
850c05b9af00 planemo upload commit 33927a87ba2eee9bf0ecdd376a66241b17b3d734
devteam
parents: 0
diff changeset
174 start_col=start_col_2,
850c05b9af00 planemo upload commit 33927a87ba2eee9bf0ecdd376a66241b17b3d734
devteam
parents: 0
diff changeset
175 end_col=end_col_2,
850c05b9af00 planemo upload commit 33927a87ba2eee9bf0ecdd376a66241b17b3d734
devteam
parents: 0
diff changeset
176 strand_col=strand_col_2,
850c05b9af00 planemo upload commit 33927a87ba2eee9bf0ecdd376a66241b17b3d734
devteam
parents: 0
diff changeset
177 fix_strand=True )
0
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
178
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
179 # Find flanking features.
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
180 out_file = open( out_fname, "w" )
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
181 try:
1
850c05b9af00 planemo upload commit 33927a87ba2eee9bf0ecdd376a66241b17b3d734
devteam
parents: 0
diff changeset
182 for result in proximal_region_finder([g1, g2], direction):
0
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
183 if type( result ) is list:
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
184 line, closest_feature = result
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
185 # Need to join outputs differently depending on file types.
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
186 if in1_gff_format:
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
187 # Output is GFF with added attribute 'closest feature.'
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
188
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
189 # Invervals are in BED coordinates; need to convert to GFF.
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
190 line = convert_bed_coords_to_gff( line )
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
191 closest_feature = convert_bed_coords_to_gff( closest_feature )
1
850c05b9af00 planemo upload commit 33927a87ba2eee9bf0ecdd376a66241b17b3d734
devteam
parents: 0
diff changeset
192
0
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
193 # Replace double quotes with single quotes in closest feature's attributes.
1
850c05b9af00 planemo upload commit 33927a87ba2eee9bf0ecdd376a66241b17b3d734
devteam
parents: 0
diff changeset
194 out_file.write( "%s closest_feature \"%s\" \n" %
850c05b9af00 planemo upload commit 33927a87ba2eee9bf0ecdd376a66241b17b3d734
devteam
parents: 0
diff changeset
195 ( "\t".join( line.fields ),
0
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
196 "\t".join( closest_feature.fields ).replace( "\"", "\\\"" )
1
850c05b9af00 planemo upload commit 33927a87ba2eee9bf0ecdd376a66241b17b3d734
devteam
parents: 0
diff changeset
197 ) )
0
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
198 else:
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
199 # Output is BED + closest feature fields.
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
200 output_line_fields = []
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
201 output_line_fields.extend( line.fields )
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
202 output_line_fields.extend( closest_feature.fields )
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
203 out_file.write( "%s\n" % ( "\t".join( output_line_fields ) ) )
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
204 else:
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
205 out_file.write( "%s\n" % result )
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
206 except ParseError, exc:
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
207 fail( "Invalid file format: %s" % str( exc ) )
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
208
1
850c05b9af00 planemo upload commit 33927a87ba2eee9bf0ecdd376a66241b17b3d734
devteam
parents: 0
diff changeset
209 print "Direction: %s" % (direction)
0
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
210 if g1.skipped > 0:
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
211 print skipped( g1, filedesc=" of 1st dataset" )
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
212 if g2.skipped > 0:
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
213 print skipped( g2, filedesc=" of 2nd dataset" )
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
214
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
215 if __name__ == "__main__":
e928e029f6eb Imported from capsule None
devteam
parents:
diff changeset
216 main()