[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.1, Fri Jun 15 21:24:49 2007 UTC revision 1.2, Tue Jun 19 18:52:40 2007 UTC
# Line 55  Line 55 
55    return $data;    return $data;
56  }  }
57    
58    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  # returns the n closest pegs, sorted by taxonomy  # returns the n closest pegs, sorted by taxonomy
263  sub get_closest_pegs {  sub get_closest_pegs {
264    my ($fig, $peg, $is_sprout, $n) = @_;    my ($fig, $peg, $is_sprout, $n) = @_;

Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3