[Bio] / FigKernelPackages / ComparedRegions.pm Repository:
ViewVC logotype

Annotation of /FigKernelPackages/ComparedRegions.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.3 - (view) (download) (as text)

1 : paczian 1.1 package ComparedRegions;
2 :    
3 :     1;
4 :    
5 :     use strict;
6 :     use warnings;
7 :    
8 :     use FIG;
9 :     use FIGV;
10 : parrello 1.3 use Tracer;
11 : paczian 1.1
12 :     sub get_compared_regions {
13 :     my ($params) = @_;
14 :    
15 :     # check parameters
16 :     my $fig = $params->{fig} || return undef;
17 :     my $peg = $params->{id} || return undef;
18 :     my $is_sprout = $params->{is_sprout} || 0;
19 :     my $region_size = $params->{region_size} || 16000;
20 :     my $number_of_genomes = $params->{number_of_genomes} || 5;
21 :    
22 :     # initialize data variable
23 :     my $data = [];
24 :    
25 :     # get the n pegs closest to the one passed
26 :     my @closest_pegs = &get_closest_pegs($fig, $peg, $is_sprout, $number_of_genomes);
27 :     unshift(@closest_pegs, $peg);
28 :    
29 :     # iterate over the returned pegs
30 :     foreach my $peg (@closest_pegs) {
31 :     my $loc = $fig->feature_location($peg);
32 :     my ($contig,$beg,$end) = $fig->boundaries_of($loc);
33 :     if ($contig && $beg && $end) {
34 :     my $mid = int(($beg + $end) / 2);
35 :     my $min = int($mid - ($region_size / 2));
36 :     my $max = int($mid + ($region_size / 2));
37 :     my $features = [];
38 :     my $feature_ids;
39 :     ($feature_ids, undef, undef) = $fig->genes_in_region($fig->genome_of($peg),$contig,$min,$max);
40 :     foreach my $fid (@$feature_ids) {
41 :     my $floc = $fig->feature_location($fid);
42 :     my ($contig1,$beg1,$end1) = $fig->boundaries_of($fig->feature_location($fid));
43 :     $beg1 = &in_bounds($min,$max,$beg1);
44 :     $end1 = &in_bounds($min,$max,$end1);
45 :     push(@$features, { 'start' => $beg1,
46 :     'stop' => $end1,
47 :     'id' => $fid });
48 :     }
49 :     push(@$data, { 'features' => $features,
50 :     'contig' => $contig,
51 :     'offset' => $mid - 8000 });
52 :     }
53 :     }
54 :    
55 :     # return the data
56 :     return $data;
57 :     }
58 :    
59 : dsouza 1.2 sub get_pattern {
60 :     my($params) = @_;
61 :     # get pattern
62 :    
63 :     my $pattern = '';
64 :     my $line;
65 :    
66 :     my $file = $params->{'file'} || return undef;
67 :    
68 :     if ( $file =~ /^([a-z0-9]+)$/ ) {
69 :     $file = 'tmp_' . $file . '.fasta';
70 :     } elsif ( $file =~ /^(tmp_[a-z0-9]+)\.cache$/ ) {
71 :     $file = $1 . '.fasta';
72 :     }
73 :    
74 :     # add path -- need to check if temp directory is correct
75 :     $file = "$FIG_Config::temp/$file";
76 :     (-e $file) || return undef;
77 :    
78 :     open(PAT, "<$file") or die "could not open file '$file': $!";
79 : parrello 1.3 while ( defined($line = <PAT>) ) {
80 : dsouza 1.2 chomp $line;
81 :     $pattern .= $line;
82 :     }
83 :     close(PAT) or die "could not close file '$file': $!";
84 :    
85 :     return $pattern;
86 :     }
87 :    
88 :     sub get_marked_locations {
89 :     my($params) = @_;
90 :     # get Bruce's hit locations
91 :    
92 :     my $fig = $params->{'fig'} || return undef;
93 :     # my $file = $params->{'file_name'} || return undef;
94 :     my $file = $params->{'file'} || 'tmp_19857906a0d1ea9aa7ce1a91a77d11b9.cache';
95 :    
96 :     if ( $file =~ /^[a-z0-9]+$/ ) {
97 :     $file = 'tmp_' . $file . '.cache';
98 :     }
99 :     # add path -- need to check if temp directory is correct
100 :     $file = "$FIG_Config::temp/$file";
101 :     (-e $file) || return undef;
102 :    
103 :     my($line, %locations);
104 :     open(HITS, "<$file") or die "could not open file '$file': $!";
105 :     $line = <HITS>; # column header line
106 :     while ( defined($line = <HITS>) ) {
107 :     my($hit) = split(/\s+/, $line);
108 :    
109 :     if ( $hit =~ /^fig\|(\S+)\.peg\.(\d+)_(\d+)\+(\d+)$/ ) {
110 :     # output from protein sequence scan
111 :     my($org_id, $peg_num, $offset, $ln) = ($1, $2, $3, $4);
112 :    
113 :     my $fid = 'fig|' . $org_id . '.peg.' . $peg_num;
114 :     my $loc = $fig->feature_location($fid);
115 :     my($contig,$beg,$end) = $fig->boundaries_of($loc);
116 :    
117 :     my($hit_beg, $hit_end);
118 :    
119 :     if ( $beg <= $end ) {
120 :     $hit_beg = $beg + ($offset * 3);
121 :     $hit_end = $hit_beg + ($ln * 3);
122 :     } else {
123 :     $hit_beg = $beg - ($offset * 3);
124 :     $hit_end = $hit_beg - ($ln * 3);
125 :     }
126 :    
127 :     push @{$locations{$org_id}{$contig}}, [$hit_beg, $hit_end];
128 :     } else {
129 :     # output from DNA scan
130 :     # this will need to be cleaned up -- check Bruce's file format
131 :     my($org_id, $loc) = split(":", $hit);
132 :     my($contig, $beg, $end) = $fig->boundaries_of($fig->feature_location($loc));
133 :     push @{$locations{$org_id}{$contig}}, [$beg, $end];
134 :     }
135 :     }
136 :     close(HITS) or die "could not close file '$file': $!";
137 :     return \%locations;
138 :     }
139 :    
140 :     sub get_hits_in_regions {
141 :     my($fig, $peg_regions, $hit_loc) = @_;
142 :     # find scan_for_matches hits (from $hit_loc) which overlap the regions
143 :     # containing the closest pegs (from $peg_regions)
144 :    
145 :     my %regions;
146 :    
147 :     # should probably sort hits & regions and perform search intelligently
148 :     foreach my $org_id ( keys %$peg_regions ) {
149 :     foreach my $contig ( keys %{$peg_regions->{$org_id}} ) {
150 :     if ( exists $hit_loc->{$org_id} and exists $hit_loc->{$org_id}{$contig} ) {
151 :     foreach my $rec ( @{$peg_regions->{$org_id}{$contig}} ) {
152 :     my($reg_beg, $reg_end, $peg) = @$rec;
153 :     my $overlapping_hits = [];
154 :     foreach my $hit ( @{$hit_loc->{$org_id}{$contig}} ) {
155 :     my($hit_beg, $hit_end) = @$hit;
156 :     if ( $fig->between($reg_beg, $hit_beg, $reg_end) or $fig->between($reg_beg, $hit_end, $reg_end) ) {
157 :     push @$overlapping_hits, [$hit_beg, $hit_end];
158 :     }
159 :     }
160 :    
161 :     if ( @$overlapping_hits ) {
162 :     push @{$regions{$org_id}{$contig}}, {'id' => $peg,
163 :     'start' => $reg_beg,
164 :     'stop' => $reg_end,
165 :     'hits' => $overlapping_hits};
166 :     }
167 :     }
168 :     }
169 :     }
170 :     }
171 :    
172 :     return \%regions;
173 :     }
174 :    
175 :     sub get_features_in_regions {
176 :     my($fig, $regions) = @_;
177 :     # add features in region to data structure containing closest pegs, with hits in defined window
178 :     foreach my $org_id ( keys %$regions ) {
179 :     foreach my $contig ( keys %{$regions->{$org_id}} ) {
180 :     foreach my $region ( @{$regions->{$org_id}{$contig}} ) {
181 :     my $peg = $region->{'id'};
182 :     my $min = $region->{'start'};
183 :     my $max = $region->{'stop'};
184 :    
185 :     my $features = [];
186 :     my $feature_ids;
187 :     ($feature_ids, undef, undef) = $fig->genes_in_region($fig->genome_of($peg),$contig,$min,$max);
188 :     foreach my $fid (@$feature_ids) {
189 :     my $floc = $fig->feature_location($fid);
190 :     my ($contig1,$beg1,$end1) = $fig->boundaries_of($fig->feature_location($fid));
191 :     $beg1 = &in_bounds($min,$max,$beg1);
192 :     $end1 = &in_bounds($min,$max,$end1);
193 :     push(@$features, { 'start' => $beg1,
194 :     'stop' => $end1,
195 :     'id' => $fid });
196 :     }
197 :    
198 :     $region->{'features'} = $features;
199 :     }
200 :    
201 :     }
202 :     }
203 :    
204 :     return $regions;
205 :     }
206 :    
207 :     sub expand_region {
208 :     my($loc, $region_size) = @_;
209 :     my $half_region = int($region_size/2);
210 :    
211 :     foreach my $org_id ( keys %$loc ) {
212 :     foreach my $contig ( keys %{$loc->{$org_id}} ) {
213 :     foreach my $rec ( @{$loc->{$org_id}{$contig}} ) {
214 :     my($beg, $end, $peg) = @$rec;
215 :     my $mid = int(($beg + $end)/2);
216 :     $rec = [($mid - $half_region), ($mid + $half_region), $peg];
217 :     }
218 :     }
219 :     }
220 :    
221 :     return $loc;
222 :     }
223 :    
224 :     sub get_closest_pegs_with_locs {
225 :     my ($params) = @_;
226 :     # return the locations of the closest pegs
227 :    
228 :     # check parameters
229 :     my $fig = $params->{fig} || return undef;
230 :     my $peg = $params->{id} || return undef;
231 :     my $is_sprout = $params->{is_sprout} || 0;
232 :     my $number_of_genomes = $params->{number_of_genomes} || 50000; # return ALL closest pegs
233 :    
234 :     # initialize data variable
235 :     my $locations = {};
236 :    
237 :     # get the n pegs closest to the one passed
238 :     my @closest_pegs = &get_closest_pegs($fig, $peg, $is_sprout, $number_of_genomes);
239 :     unshift(@closest_pegs, $peg);
240 :    
241 :     # iterate over the returned pegs
242 :     foreach my $peg (@closest_pegs) {
243 :     my $loc = $fig->feature_location($peg);
244 :     my ($contig,$beg,$end) = $fig->boundaries_of($loc);
245 :     if ($contig && $beg && $end) {
246 :     my $org_id;
247 :     if ( $peg =~ /^fig\|([^\|]+)\|/ ) {
248 :     $org_id = $1;
249 :     } else {
250 :     my $org_name = $fig->org_of($peg);
251 :     $org_id = $fig->orgid_of_orgname($org_name);
252 :     }
253 :    
254 :     # add peg and location
255 :     push @{ $locations->{$org_id}{$contig} }, [$beg, $end, $peg];
256 :     }
257 :     }
258 :    
259 :     # return the locations of the closest pegs
260 :     return $locations;
261 :     }
262 :    
263 : parrello 1.3
264 : paczian 1.1 # returns the n closest pegs, sorted by taxonomy
265 :     sub get_closest_pegs {
266 :     my ($fig, $peg, $is_sprout, $n) = @_;
267 :    
268 :     my($id2,$d,$peg2,$i);
269 :    
270 :     my @closest;
271 :     if ($is_sprout) {
272 :     @closest = map { $_->[0] } sort { $a->[1] <=> $b->[1] } $fig->bbhs($peg, 1.0e-10);
273 :     } else {
274 :     @closest = map { $id2 = $_->id2; ($id2 =~ /^fig\|/) ? $id2 : () } $fig->sims($peg,&FIG::max(20,$n*4),1.0e-20,"fig",&FIG::max(20,$n*4));
275 :     }
276 :    
277 :     if (@closest >= ($n-1)) {
278 :     $#closest = $n-2 ;
279 :     }
280 :     my %closest = map { $_ => 1 } @closest;
281 :    
282 :     # there are dragons flying around...
283 :     my @pinned_to = grep { ($_ ne $peg) && (! $closest{$_}) } $fig->in_pch_pin_with($peg);
284 :     my $g1 = $fig->genome_of($peg);
285 :     @pinned_to = map {$_->[1] } sort { $a->[0] <=> $b->[0] } map { $peg2 = $_; $d = $fig->crude_estimate_of_distance($g1,$fig->genome_of($peg2)); [$d,$peg2] } @pinned_to;
286 :    
287 :     if (@closest == ($n-1)) {
288 :     $#closest = ($n - 2) - &FIG::min(scalar @pinned_to,int($n/2));
289 :     for ($i=0; ($i < @pinned_to) && (@closest < ($n-1)); $i++) {
290 :     if (! $closest{$pinned_to[$i]}) {
291 :     $closest{$pinned_to[$i]} = 1;
292 :     push(@closest,$pinned_to[$i]);
293 :     }
294 :     }
295 :     }
296 :    
297 :     if ($fig->possibly_truncated($peg)) {
298 :     push(@closest, &possible_extensions($fig, $peg, \@closest));
299 :     }
300 :     @closest = $fig->sort_fids_by_taxonomy(@closest);
301 :    
302 :     return @closest;
303 :    
304 :     }
305 :    
306 :     sub in_bounds {
307 :     my($min,$max,$x) = @_;
308 :    
309 :     if ($x < $min) { return $min }
310 :     elsif ($x > $max) { return $max }
311 :     else { return $x }
312 :     }
313 :    
314 :     sub possible_extensions {
315 :     my($fig, $peg,$closest_pegs) = @_;
316 :     my($g,$sim,$id2,$peg1,%poss);
317 :    
318 :     $g = &FIG::genome_of($peg);
319 :    
320 :     foreach $peg1 (@$closest_pegs) {
321 :     if ($g ne &genome_of($peg1)) {
322 :     foreach $sim ($fig->sims($peg1,500,1.0e-5,"all")) {
323 :     $id2 = $sim->id2;
324 :     if (($id2 ne $peg) && ($id2 =~ /^fig\|$g\./) && $fig->possibly_truncated($id2)) {
325 :     $poss{$id2} = 1;
326 :     }
327 :     }
328 :     }
329 :     }
330 :     return keys(%poss);
331 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3