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

Diff of /FigKernelPackages/ComparedRegions.pm

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.3, Tue Jun 19 21:12:53 2007 UTC revision 1.4, Thu Jun 28 16:51:39 2007 UTC
# Line 56  Line 56 
56    return $data;    return $data;
57  }  }
58    
59  sub get_pattern {  sub get_scan_for_matches_pattern {
60      my($params) = @_;      my($params) = @_;
61      # get pattern      # get scan_for_matches pattern from the 'fasta' file
   
     my $pattern = '';  
     my $line;  
62    
63        # parse input file name
64      my $file = $params->{'file'} || return undef;      my $file = $params->{'file'} || return undef;
65        if ($file =~ /^([a-z0-9]+)$/ or $file =~ /^tmp_([a-z0-9]+)\.cache$/) {
66      if ( $file =~ /^([a-z0-9]+)$/ ) {          $file = 'tmp_' . $1 . '.fasta';
67          $file = 'tmp_' . $file . '.fasta';      } else {
68      } elsif ( $file =~ /^(tmp_[a-z0-9]+)\.cache$/ ) {          return undef;
         $file = $1 . '.fasta';  
69      }      }
70    
71      # add path -- need to check if temp directory is correct      # add path
72      $file    = "$FIG_Config::temp/$file";      $file    = "$FIG_Config::temp/$file";
73        # check if file exists
74      (-e $file) || return undef;      (-e $file) || return undef;
75    
76        my $pattern = '';
77        my $line;
78    
79      open(PAT, "<$file") or die "could not open file '$file': $!";      open(PAT, "<$file") or die "could not open file '$file': $!";
80      while ( defined($line = <PAT>) ) {      while ( defined($line = <PAT>) ) {
81          chomp $line;          #chomp $line;
82          $pattern .= $line;          $pattern .= $line;
83      }      }
84      close(PAT) or die "could not close file '$file': $!";      close(PAT) or die "could not close file '$file': $!";
# Line 85  Line 86 
86      return $pattern;      return $pattern;
87  }  }
88    
89  sub get_marked_locations {  sub get_scan_for_matches_hits {
90      my($params) = @_;      my($params) = @_;
91      # get Bruce's hit locations      # get Bruce's hit locations
92    
93      my $fig  = $params->{'fig'} || return undef;      my $fig  = $params->{'fig'} || return undef;
94      # my $file = $params->{'file_name'} || return undef;      my $file = $params->{'file'} || return undef;
     my $file = $params->{'file'} || 'tmp_19857906a0d1ea9aa7ce1a91a77d11b9.cache';  
95    
96        # parse input file name
97      if ( $file =~ /^[a-z0-9]+$/ ) {      if ( $file =~ /^[a-z0-9]+$/ ) {
98          $file = 'tmp_' . $file . '.cache';          $file = 'tmp_' . $file . '.cache';
99      }      }
100      # add path -- need to check if temp directory is correct  
101        # add path
102      $file    = "$FIG_Config::temp/$file";      $file    = "$FIG_Config::temp/$file";
103        # check if file exists
104      (-e $file) || return undef;      (-e $file) || return undef;
105    
106      my($line, %locations);      my($line, %locations);
107      open(HITS, "<$file") or die "could not open file '$file': $!";      open(HITS, "<$file") or die "could not open file '$file': $!";
108      $line = <HITS>;  # column header line      $line = <HITS>;  # ignore column header line
109      while ( defined($line = <HITS>) ) {      while ( defined($line = <HITS>) ) {
110          my($hit) = split(/\s+/, $line);          my($hit) = split(/\s+/, $line);
111    
# Line 112  Line 115 
115    
116              my $fid = 'fig|' . $org_id . '.peg.' . $peg_num;              my $fid = 'fig|' . $org_id . '.peg.' . $peg_num;
117              my $loc = $fig->feature_location($fid);              my $loc = $fig->feature_location($fid);
118    
119                if (not $fig->is_deleted_fid($fid))
120                {
121              my($contig,$beg,$end) = $fig->boundaries_of($loc);              my($contig,$beg,$end) = $fig->boundaries_of($loc);
122    
123              my($hit_beg, $hit_end);              my($hit_beg, $hit_end);
# Line 125  Line 131 
131              }              }
132    
133              push @{$locations{$org_id}{$contig}}, [$hit_beg, $hit_end];              push @{$locations{$org_id}{$contig}}, [$hit_beg, $hit_end];
134                }
135          } else {          } else {
136              # output from DNA scan              # output from DNA scan
137              # this will need to be cleaned up -- check Bruce's file format              # this will need to be cleaned up depending on the format Bruce puts the output in
138              my($org_id, $loc) = split(":", $hit);              my($org_id, $loc) = split(":", $hit);
139              my($contig, $beg, $end) = $fig->boundaries_of($fig->feature_location($loc));              my($contig, $beg, $end) = $fig->boundaries_of($fig->feature_location($loc));
140              push @{$locations{$org_id}{$contig}}, [$beg, $end];              push @{$locations{$org_id}{$contig}}, [$beg, $end];
# Line 137  Line 144 
144      return \%locations;      return \%locations;
145  }  }
146    
147  sub get_hits_in_regions {  sub add_hits_in_regions {
148      my($fig, $peg_regions, $hit_loc) = @_;      my($fig, $regions, $patscan_hits) = @_;
149      # find scan_for_matches hits (from $hit_loc) which overlap the regions      # find scan_for_matches hits (from $patscan_hits) which overlap the regions
150      # containing the closest pegs (from $peg_regions)      # containing the closest pegs (from $peg_regions)
151    
152      my %regions;      # should probably sort regions and perform search intelligently
153        foreach my $org_id (keys %$regions) {
154      # should probably sort hits & regions and perform search intelligently          foreach my $contig (keys %{$regions->{$org_id}}) {
155      foreach my $org_id ( keys %$peg_regions ) {              if (exists $patscan_hits->{$org_id} and exists $patscan_hits->{$org_id}{$contig}) {
156          foreach my $contig ( keys %{$peg_regions->{$org_id}} ) {                  foreach my $region (@{$regions->{$org_id}{$contig}}) {
157              if ( exists $hit_loc->{$org_id} and exists $hit_loc->{$org_id}{$contig} ) {                      my($reg_beg, $reg_end) = ($region->{'reg_beg'}, $region->{'reg_end'});
158                  foreach my $rec ( @{$peg_regions->{$org_id}{$contig}} ) {                      # build up list of hits which overlap the region
                     my($reg_beg, $reg_end, $peg) = @$rec;  
159                      my $overlapping_hits = [];                      my $overlapping_hits = [];
160                      foreach my $hit ( @{$hit_loc->{$org_id}{$contig}} ) {                      foreach my $hit (@{$patscan_hits->{$org_id}{$contig}}) {
161                          my($hit_beg, $hit_end) = @$hit;                          my($hit_beg, $hit_end) = @$hit;
162                          if ( $fig->between($reg_beg, $hit_beg, $reg_end) or $fig->between($reg_beg, $hit_end, $reg_end) ) {                          if ( $fig->between($reg_beg, $hit_beg, $reg_end) or $fig->between($reg_beg, $hit_end, $reg_end) ) {
163                              push @$overlapping_hits, [$hit_beg, $hit_end];                              push @$overlapping_hits, [$hit_beg, $hit_end];
# Line 159  Line 165 
165                      }                      }
166    
167                      if ( @$overlapping_hits ) {                      if ( @$overlapping_hits ) {
168                          push @{$regions{$org_id}{$contig}}, {'id'    => $peg,                          # add the overlapping hits in the region
169                                                               'start' => $reg_beg,                          $region->{'hits'} = $overlapping_hits;
                                                              'stop'  => $reg_end,  
                                                              'hits'  => $overlapping_hits};  
170                      }                      }
171                  }                  }
172              }              }
173          }          }
174      }      }
175    
176      return \%regions;      return $regions;
177  }  }
178    
179  sub get_features_in_regions {  sub add_features_in_regions {
180      my($fig, $regions) = @_;      my($fig, $regions) = @_;
181      # add features in region to data structure containing closest pegs, with hits in defined window      # add features in region to data structure containing closest pegs
182      foreach my $org_id ( keys %$regions ) {      foreach my $org_id ( keys %$regions ) {
183          foreach my $contig ( keys %{$regions->{$org_id}} ) {          foreach my $contig ( keys %{$regions->{$org_id}} ) {
184              foreach my $region ( @{$regions->{$org_id}{$contig}} ) {              foreach my $region ( @{$regions->{$org_id}{$contig}} ) {
185                  my $peg = $region->{'id'};                  # get closest peg and boundaries of region
186                  my $min = $region->{'start'};                  my $peg = $region->{'peg'};
187                  my $max = $region->{'stop'};                  my $min = $region->{'reg_beg'};
188                    my $max = $region->{'reg_end'};
189    
190                    # get list of features which fall in this region
191                    my($feature_ids) = $fig->genes_in_region($fig->genome_of($peg),$contig,$min,$max);
192    
193                    # make up a list of all these features
194                  my $features = [];                  my $features = [];
                 my $feature_ids;  
                 ($feature_ids, undef, undef) = $fig->genes_in_region($fig->genome_of($peg),$contig,$min,$max);  
195                  foreach my $fid (@$feature_ids) {                  foreach my $fid (@$feature_ids) {
196                      my $floc = $fig->feature_location($fid);                      my $floc = $fig->feature_location($fid);
197                      my ($contig1,$beg1,$end1) = $fig->boundaries_of($fig->feature_location($fid));                      my ($contig1,$beg1,$end1) = $fig->boundaries_of($fig->feature_location($fid));
198                      $beg1 = &in_bounds($min,$max,$beg1);                      $beg1 = &in_bounds($min,$max,$beg1);
199                      $end1 = &in_bounds($min,$max,$end1);                      $end1 = &in_bounds($min,$max,$end1);
200                      push(@$features, { 'start' => $beg1,                      push(@$features, { 'feat_beg' => $beg1,
201                                         'stop' => $end1,                                         'feat_end' => $end1,
202                                         'id' => $fid });                                         'feat_id'  => $fid });
203                  }                  }
204    
205                    # add the feature data to the region
206                  $region->{'features'} = $features;                  $region->{'features'} = $features;
207              }              }
208            }
209        }
210    
211        return $regions;
212    }
213    
214    sub sort_regions {
215        my($fig, $regions) = @_;
216    
217        # sort regions based on location on the contig, regardless of strand
218        # i.e. sort on the 'left-most' location of the region on the contig
219        foreach my $org_id (keys %$regions) {
220            foreach my $contig (keys %{$regions->{$org_id}}) {
221                my $contig_regions = $regions->{$org_id}{$contig};
222    
223                if (@$contig_regions > 1) {
224                    my @sorted_regions = map {$_->[0]}
225                                           sort {$a->[1] <=> $b->[1]}
226                                             map {[$_, $fig->min($_->{'reg_beg'}, $_->{'reg_end'})]}
227                                               @$contig_regions;
228    
229                    $regions->{$org_id}{$contig} = \@sorted_regions;
230                }
231          }          }
232      }      }
233    
234      return $regions;      return $regions;
235  }  }
236    
237  sub expand_region {  sub collapse_regions {
238      my($loc, $region_size) = @_;      my($fig, $regions, $window_size) = @_;
239      my $half_region = int($region_size/2);      # collapse regions which display regions which show basically the same information
240        # if two (or more) regions:
241        # a. overlap by more than half the region size, and
242        # b. the pattern match is located entirely within the overlap,
243        # only one of the regions will be retained for display.
244        # since the regions are centred on the 'close' peg, this should give a reasonable
245        # neighborhood.
246        # if regions get collapsed, the region displayed is the 'left-most' on the contig.
247    
248        my $half_region = 0.5 * $window_size;
249    
250      foreach my $org_id ( keys %$loc ) {      foreach my $org_id (keys %$regions) {
251          foreach my $contig ( keys %{$loc->{$org_id}} ) {          foreach my $contig (keys %{$regions->{$org_id}}) {
252              foreach my $rec ( @{$loc->{$org_id}{$contig}} ) {              my $contig_regions = $regions->{$org_id}{$contig};
253                  my($beg, $end, $peg) = @$rec;  
254                  my $mid = int(($beg + $end)/2);              if (@$contig_regions > 1) {
255                  $rec = [($mid - $half_region), ($mid + $half_region), $peg];                  # hash for index of regions which should not be displayed
256                    my %discard;
257                    for (my $i = 0; $i < @$contig_regions; $i++) {
258                        for (my $j = ($i+1); $j < @$contig_regions; $j++) {
259                            my $overlap = &regions_overlap($contig_regions->[$i], $contig_regions->[$j]);
260    
261                            if ($overlap and
262                                ($overlap > $half_region) and
263                                &all_hits_overlap_region($contig_regions->[$i], $contig_regions->[$j]{'hits'}))
264                            {
265                                # region i display contains all relevant details from region j,
266                                # region j can be discarded.
267    
268                                $discard{$j} = 1;
269                            }
270                            else
271                            {
272                                # region i display does not contain all relevant details from region j,
273                                # region j needs to be kept.
274                                # next value of $i needs to be $j
275                                $i = $j - 1;
276                                # set $j so that we exit second loop
277                                $j = @$contig_regions;
278              }              }
279          }          }
280      }      }
281    
282      return $loc;                  my @keep = grep {not exists $discard{$_}} (0..$#{$contig_regions});
283                    $regions->{$org_id}{$contig} = [@$contig_regions[@keep]];
284                }
285            }
286        }
287        return $regions;
288    }
289    
290    sub all_hits_overlap_region {
291        my($region, $hits) = @_;
292        # Check whether ALL hits overlap the region, if yes return 1, if no return 0
293        # Return 1 if no hits present
294        # This will have problems for hits which are longer than the region displayed
295    
296        my($reg_beg, $reg_end) = ($region->{'reg_beg'}, $region->{'reg_end'});
297    
298        foreach my $hit (@$hits) {
299            my($hit_beg, $hit_end) = sort {$a <=> $b} @$hit;
300            # get length of hit
301            my $hit_ln = $hit_end - $hit_beg + 1;
302            if (&overlap($reg_beg, $reg_end, $hit_beg, $hit_end) < $hit_ln) {
303                # overlap is less than length of hit
304                return 0;
305            }
306        }
307    
308        # all hits overlap the region
309        return 1;
310    }
311    
312    sub regions_overlap {
313        my($r1, $r2) = @_;
314        # return overlap in bp between two regions
315        return &overlap($r1->{'reg_beg'}, $r1->{'reg_end'}, $r2->{'reg_beg'}, $r2->{'reg_end'});
316    }
317    
318    sub overlap {
319        my($x1, $x2, $y1, $y2) = @_;
320        # return 1 if regions [$x1,$x2] and [$y1,$y2] overlap, otherwise return 0
321    
322        my $overlap = 0;
323        ($x1, $x2) = sort {$a <=> $b} ($x1, $x2);
324        ($y1, $y2) = sort {$a <=> $b} ($y1, $y2);
325    
326        if (($x2 < $y1) or ($y2 < $x1)) {
327            $overlap = 0;
328        } elsif ($x1 <= $y1) {
329            if ($y2 <= $x2) {
330                $overlap = $y2 - $y1 + 1;
331            } else {
332                $overlap = $x2 - $y1 + 1;
333            }
334        } elsif ($y1 <= $x1) {
335            if ($x2 <= $y2) {
336                $overlap = $x2 - $x1 + 1;
337            } else {
338                $overlap = $y2 - $x1 + 1;
339            }
340  }  }
341    
342  sub get_closest_pegs_with_locs {      return $overlap;
343    }
344    
345    sub closest_peg_regions {
346      my ($params) = @_;      my ($params) = @_;
     # return the locations of the closest pegs  
347    
348      # check parameters      # check parameters
349      my $fig = $params->{fig} || return undef;      my $fig          = $params->{'fig'} || return undef;
350      my $peg = $params->{id}  || return undef;      my $closest_pegs = $params->{'closest_pegs'} || return undef;
351      my $is_sprout = $params->{is_sprout} || 0;      my $window_size  = $params->{'window_size'} || return undef;
352      my $number_of_genomes = $params->{number_of_genomes} || 50000; # return ALL closest pegs      my $half_region  = int($window_size/2);
353    
354      # initialize data variable      # initialize data variable
355      my $locations = {};      my $regions = {};
   
     # get the n pegs closest to the one passed  
     my @closest_pegs = &get_closest_pegs($fig, $peg, $is_sprout, $number_of_genomes);  
     unshift(@closest_pegs, $peg);  
356    
357      # iterate over the returned pegs      # iterate over the pegs
358      foreach my $peg (@closest_pegs) {      foreach my $peg (@$closest_pegs) {
359            # get location of peg
360          my $loc = $fig->feature_location($peg);          my $loc = $fig->feature_location($peg);
361          my ($contig,$beg,$end) = $fig->boundaries_of($loc);          my($peg_contig,$peg_beg,$peg_end) = $fig->boundaries_of($loc);
362          if ($contig && $beg && $end) {          if ($peg_contig && $peg_beg && $peg_end) {
363                # get organism ID
364              my $org_id;              my $org_id;
365              if ( $peg =~ /^fig\|([^\|]+)\|/ ) {              if ( $peg =~ /^fig\|([^\|]+)\|/ ) {
366                  $org_id = $1;                  $org_id = $1;
# Line 251  Line 369 
369                  $org_id      = $fig->orgid_of_orgname($org_name);                  $org_id      = $fig->orgid_of_orgname($org_name);
370              }              }
371    
372              # add peg and location              # find mid-point of peg
373              push @{ $locations->{$org_id}{$contig} }, [$beg, $end, $peg];              my $mid = int(($peg_beg + $peg_end)/2);
374    
375                # add peg, peg location and region location to hash
376                push @{ $regions->{$org_id}{$peg_contig} }, {'peg'     => $peg,
377                                                             'peg_beg' => $peg_beg,
378                                                             'peg_end' => $peg_end,
379                                                             'reg_beg' => ($mid - $half_region),
380                                                             'reg_end' => ($mid + $half_region)};
381          }          }
382      }      }
383    
384      # return the locations of the closest pegs      # return regions data
385      return $locations;      return $regions;
386  }  }
387    
   
 # returns the n closest pegs, sorted by taxonomy  
388  sub get_closest_pegs {  sub get_closest_pegs {
389    my ($fig, $peg, $is_sprout, $n) = @_;      my ($params) = @_;
390        # returns the n closest pegs, sorted by taxonomy
391    
392        # check parameters
393        my $fig       = $params->{'fig'} || return undef;
394        my $peg       = $params->{'id'}  || return undef;
395        my $is_sprout = $params->{'is_sprout'} || 0;
396        my $n         = $params->{'number_of_genomes'} || 50000; # return ALL closest pegs
397    
398        # get the n pegs closest to the one passed
399    
400    my($id2,$d,$peg2,$i);    my($id2,$d,$peg2,$i);
401    
# Line 299  Line 431 
431    }    }
432    @closest = $fig->sort_fids_by_taxonomy(@closest);    @closest = $fig->sort_fids_by_taxonomy(@closest);
433    
434    return @closest;      unshift(@closest, $peg);
435        return \@closest;
436  }  }
437    
438  sub in_bounds {  sub in_bounds {

Legend:
Removed from v.1.3  
changed lines
  Added in v.1.4

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3