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

Annotation of /FigKernelPackages/ComparedRegions.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3