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

Annotation of /FigKernelScripts/close_connections.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3