[Bio] / FigKernelScripts / close_connections.pl Repository:
ViewVC logotype

Annotation of /FigKernelScripts/close_connections.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : overbeek 1.1 #!/usr/bin/env /Volumes/seagate/FIGdisk/env/mac/bin/perl
2 :    
3 :     BEGIN {
4 :     @INC = qw(
5 :     /Volumes/seagate/FIGdisk/dist/releases/ross/mac/lib
6 :     /Volumes/seagate/FIGdisk/dist/releases/ross/mac/lib/FigKernelPackages
7 :     /Volumes/seagate/FIGdisk/dist/releases/ross/mac/lib/WebApplication
8 :     /Volumes/seagate/FIGdisk/dist/releases/ross/mac/lib/FortyEight
9 :     /Volumes/seagate/FIGdisk/dist/releases/ross/mac/lib/PPO
10 :     /Volumes/seagate/FIGdisk/dist/ross/mac/lib
11 :     /Volumes/seagate/FIGdisk/dist/ross/mac/lib/FigKernelPackages
12 :     /Volumes/seagate/FIGdisk/env/mac/lib/perl5/5.8.7/darwin-2level
13 :     /Volumes/seagate/FIGdisk/env/mac/lib/perl5/5.8.7
14 :     /Volumes/seagate/FIGdisk/env/mac/lib/perl5/site_perl/5.8.7/darwin-2level
15 :     /Volumes/seagate/FIGdisk/env/mac/lib/perl5/site_perl/5.8.7
16 :     /Volumes/seagate/FIGdisk/env/mac/lib/perl5/site_perl
17 :     .
18 :     /Volumes/seagate/FIGdisk/config
19 :    
20 :     );
21 :     }
22 :     use Data::Dumper;
23 :     use Carp;
24 :     use FIG_Config;
25 :     $ENV{'BLASTMAT'} = "/Volumes/seagate/FIGdisk/BLASTMAT";
26 :     $ENV{'FIG_HOME'} = "/Volumes/seagate/FIGdisk";
27 :     # end of tool_hdr
28 :     ########################################################################
29 :    
30 :     use FIG;
31 :     my $fig = new FIG;
32 :    
33 :     use strict;
34 :    
35 :     # usage: close_connections [MinNormSc] < SimsFile > Connections
36 :    
37 :     # a minimum normalized bit score (bit score divided by the max of the lengths of the input seqs)
38 :     # is used to filter similarities. We use a value of 1.3 as a default, but I have not put
39 :     # any real effort into finding what an optimal value might be.
40 :     my $min_norm_sc = (@ARGV > 0) ? $ARGV[0] : 1.3;
41 :    
42 :    
43 :     my(%genomes,%sims);
44 :     while (defined($_ = <STDIN>))
45 :     {
46 :     chop;
47 :     my($id1,$id2,undef,undef,undef,undef,undef,undef,$ln1,$ln2,$bsc) = split(/\t/,$_);
48 :     $genomes{&FIG::genome_of($id1)} = 1;
49 :     $genomes{&FIG::genome_of($id2)} = 1;
50 :    
51 :     my $norm_bsc = sprintf("%.3f",$bsc / &FIG::max($ln1,$ln2));
52 :     if ($norm_bsc >= $min_norm_sc)
53 :     {
54 :     push(@{$sims{$id1}},[$id2,$norm_bsc]);
55 :     push(@{$sims{$id2}},[$id1,$norm_bsc]);
56 :     }
57 :     }
58 :    
59 :     my @genomes = sort { a <=> $b } keys(%genomes);
60 :     if (@genomes != 2)
61 :     {
62 :     print &Dumper(\%genomes);
63 :     die "We should have just 2 genomes";
64 :     }
65 :    
66 :     my($peg,$i);
67 :     foreach $peg (keys(%sims))
68 :     {
69 :     my @tmp = sort { $b->[1] <=> $a->[1] } @{$sims{$peg}};
70 :     for ($i=0; ($i < @tmp) && ($tmp[$i]->[1] >= ($tmp[0]->[1] - 0.1)); $i++) {}
71 :     $#tmp = $i-1;
72 :     $sims{$peg} = [@tmp];
73 :     }
74 :    
75 :     # We compute two hashes, one for each genome. Each hash contains an entry for each PEG
76 :     # in that genome. The value for a given PEG is [contig,N] where N increases by 1 as you
77 :     # go through the PEGs on a contig. Thus, two PEGs are adjacent iff they have the same
78 :     # contigs and their N values differ by 1.
79 :     #
80 :    
81 :     my $locs1 = &get_locations($fig,$genomes[0]);
82 :     my $locs2 = &get_locations($fig,$genomes[1]);
83 :    
84 :     my @pegs1 = sort { ($locs1->{$a}->[0] cmp $locs1->{$b}->[0]) or
85 :     ($locs1->{$a}->[1] <=> $locs1->{$b}->[1]) }
86 :     keys(%$locs1);
87 :    
88 :     $i = 0;
89 :     my($x,$j,$start,$peg1B,$peg2B,$peg2P);
90 :     while ($i < @pegs1)
91 :     {
92 :     if ($peg2B = &clear_bbh($pegs1[$i],\%sims))
93 :     {
94 :     $peg1B = $pegs1[$i];
95 :     $start = $i;
96 :     $i++;
97 :     my $last = [$peg1B,$peg2B];
98 :     my @matches = ($last);
99 :     my $got_end = 0;
100 :     while (($i < @pegs1) && (! $got_end) && ($peg2P = &poss($fig,$last,$pegs1[$i],\%sims,$locs1,$locs2)))
101 :     {
102 :     $last = [$pegs1[$i],$peg2P];
103 :     push(@matches,$last);
104 :     if (&are_clear_bbhs($pegs1[$i],$peg2P,\%sims))
105 :     {
106 :     $got_end = 1;
107 :     }
108 :     else
109 :     {
110 :     $i++;
111 :     }
112 :     }
113 :    
114 :     if ($got_end && (@matches <= 10))
115 :     {
116 :     my $pair;
117 :     foreach $pair (@matches)
118 :     {
119 :     print join("\t",@$pair),"\n";
120 :     }
121 :     }
122 :     }
123 :     else
124 :     {
125 :     $i++;
126 :     }
127 :     }
128 :    
129 :     sub clear_bbh {
130 :     my($peg1,$sims) = @_;
131 :     my($x,$y,$peg2);
132 :    
133 :     if (($x = $sims->{$peg1}) && (@$x == 1) &&
134 :     ($peg2 = $x->[0]->[0]) &&
135 :     ($y = $sims->{$peg2}) && (@$y == 1) &&
136 :     ($y->[0]->[0] eq $peg1))
137 :     {
138 :     return $peg2;
139 :     }
140 :     else
141 :     {
142 :     return undef;
143 :     }
144 :     }
145 :    
146 :     sub poss {
147 :     my($fig,$last,$peg1b,$sims,$locs1,$locs2) = @_;
148 :     my($x,$y,$i,$peg2b);
149 :    
150 :     my($peg1a,$peg2a) = @$last;
151 :     if ($x = $sims->{$peg1b})
152 :     {
153 :     for ($i=0; ($i < @$x) && (! &adjacent($peg2a,$x->[$i]->[0],$locs2)); $i++) {}
154 :     if ($i < @$x)
155 :     {
156 :     $peg2b = $x->[$i]->[0];
157 :     if ($y = $sims->{$peg2b})
158 :     {
159 :     for ($i=0; ($i < @$y) && ($y->[$i]->[0] ne $peg1b); $i++) {}
160 :     if ($i < @$y)
161 :     {
162 :     return $peg2b;
163 :     }
164 :     }
165 :     }
166 :     }
167 :     return undef;
168 :     }
169 :    
170 :     sub adjacent {
171 :     my($peg1,$peg2,$locs) = @_;
172 :     my($x,$y);
173 :    
174 :     return (($x = $locs->{$peg1}) && ($y = $locs->{$peg2}) && ($x->[0] eq $y->[0]) &&
175 :     (abs($x->[1] - $y->[1]) == 1));
176 :     }
177 :    
178 :     sub are_clear_bbhs {
179 :     my($peg1,$peg2,$sims) = @_;
180 :     my($x,$y);
181 :    
182 :     return (($x = $sims->{$peg1}) && (@$x == 1) && ($x->[0]->[0] eq $peg2) &&
183 :     ($y = $sims->{$peg2}) && (@$y == 1) && ($y->[0]->[0] eq $peg1));
184 :     }
185 :    
186 :     sub get_locations {
187 :     my($fig,$genome,$sims) = @_;
188 :     my($peg,$i);
189 :    
190 :     my $locs = {};
191 :     my @pegs_and_locs = ();
192 :     foreach $peg ($fig->all_features($genome,"peg"))
193 :     {
194 :     my @loc = $fig->feature_location($peg);
195 :     if ((@loc == 1) && ($loc[0] =~ /^(\S+)_(\d+)_(\d+)$/))
196 :     {
197 :     push(@pegs_and_locs,[$peg,$1,int(($2+$3)/2)]);
198 :     }
199 :     }
200 :     @pegs_and_locs = sort { ($a->[1] cmp $b->[1]) or ($a->[2] <=> $b->[2]) } @pegs_and_locs;
201 :     for ($i=0; ($i < @pegs_and_locs); $i++)
202 :     {
203 :     $locs->{$pegs_and_locs[$i]->[0]} = [$pegs_and_locs[$i]->[1],$i];
204 :     }
205 :     return $locs;
206 :     }
207 :    

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3