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

Annotation of /FigKernelPackages/ComparedRegions.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (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 :     # returns the n closest pegs, sorted by taxonomy
59 :     sub get_closest_pegs {
60 :     my ($fig, $peg, $is_sprout, $n) = @_;
61 :    
62 :     my($id2,$d,$peg2,$i);
63 :    
64 :     my @closest;
65 :     if ($is_sprout) {
66 :     @closest = map { $_->[0] } sort { $a->[1] <=> $b->[1] } $fig->bbhs($peg, 1.0e-10);
67 :     } else {
68 :     @closest = map { $id2 = $_->id2; ($id2 =~ /^fig\|/) ? $id2 : () } $fig->sims($peg,&FIG::max(20,$n*4),1.0e-20,"fig",&FIG::max(20,$n*4));
69 :     }
70 :    
71 :     if (@closest >= ($n-1)) {
72 :     $#closest = $n-2 ;
73 :     }
74 :     my %closest = map { $_ => 1 } @closest;
75 :    
76 :     # there are dragons flying around...
77 :     my @pinned_to = grep { ($_ ne $peg) && (! $closest{$_}) } $fig->in_pch_pin_with($peg);
78 :     my $g1 = $fig->genome_of($peg);
79 :     @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;
80 :    
81 :     if (@closest == ($n-1)) {
82 :     $#closest = ($n - 2) - &FIG::min(scalar @pinned_to,int($n/2));
83 :     for ($i=0; ($i < @pinned_to) && (@closest < ($n-1)); $i++) {
84 :     if (! $closest{$pinned_to[$i]}) {
85 :     $closest{$pinned_to[$i]} = 1;
86 :     push(@closest,$pinned_to[$i]);
87 :     }
88 :     }
89 :     }
90 :    
91 :     if ($fig->possibly_truncated($peg)) {
92 :     push(@closest, &possible_extensions($fig, $peg, \@closest));
93 :     }
94 :     @closest = $fig->sort_fids_by_taxonomy(@closest);
95 :    
96 :     return @closest;
97 :    
98 :     }
99 :    
100 :     sub in_bounds {
101 :     my($min,$max,$x) = @_;
102 :    
103 :     if ($x < $min) { return $min }
104 :     elsif ($x > $max) { return $max }
105 :     else { return $x }
106 :     }
107 :    
108 :     sub possible_extensions {
109 :     my($fig, $peg,$closest_pegs) = @_;
110 :     my($g,$sim,$id2,$peg1,%poss);
111 :    
112 :     $g = &FIG::genome_of($peg);
113 :    
114 :     foreach $peg1 (@$closest_pegs) {
115 :     if ($g ne &genome_of($peg1)) {
116 :     foreach $sim ($fig->sims($peg1,500,1.0e-5,"all")) {
117 :     $id2 = $sim->id2;
118 :     if (($id2 ne $peg) && ($id2 =~ /^fig\|$g\./) && $fig->possibly_truncated($id2)) {
119 :     $poss{$id2} = 1;
120 :     }
121 :     }
122 :     }
123 :     }
124 :     return keys(%poss);
125 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3