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

Diff of /FigKernelScripts/connections_between.pl

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

revision 1.7, Wed Dec 17 16:19:51 2008 UTC revision 1.8, Sun Jan 11 13:43:37 2009 UTC
# Line 1  Line 1 
1    ########################################################################
2  #  #
3  # Copyright (c) 2003-2006 University of Chicago and Fellowship  # Copyright (c) 2003-2006 University of Chicago and Fellowship
4  # for Interpretations of Genomes. All Rights Reserved.  # for Interpretations of Genomes. All Rights Reserved.
# Line 23  Line 24 
24    
25  my $fig;  my $fig;
26    
27  my($simsD,$iden_diff,$neigh_diff,$req_neigh);  my($simsD,$iden_diff);
28    
29  my $usage = "usage: connections_between [-orgdir orgdir] Sims IdenDiff NeighDiff  > connections";  
30    my $usage = "usage: connections_between [-orgdir orgdir]  [-cluster N1] [-insub N2] Sims IdenDiff > connections";
31    
32  my $orgdir;  my $orgdir;
33    my $must_cluster = 0;
34    my $must_be_in_sub = 0;
35    
36  while ((@ARGV > 0) && ($ARGV[0] =~ /^-/))  while ((@ARGV > 0) && ($ARGV[0] =~ /^-/))
37  {  {
38      my $arg = shift @ARGV;      my $arg = shift @ARGV;
# Line 35  Line 40 
40      {      {
41          $orgdir = shift @ARGV;          $orgdir = shift @ARGV;
42      }      }
43        elsif ($arg =~ /-cluster/)
44        {
45            $must_cluster = shift @ARGV;
46        }
47        elsif ($arg =~ /-insub/)
48        {
49            $must_be_in_sub = shift @ARGV;
50        }
51      else      else
52      {      {
53          die "Unknown option $arg\n";          die "Unknown option $arg\n";
# Line 52  Line 65 
65    
66  (  (
67   ($simsD      = shift @ARGV) &&   ($simsD      = shift @ARGV) &&
68   ($iden_diff  = shift @ARGV) &&   ($iden_diff  = shift @ARGV)
  ($neigh_diff = shift @ARGV)  
69  )  )
70      || die $usage;      || die $usage;
71    
# Line 167  Line 179 
179                  my $conn;                  my $conn;
180                  if ($peg2 = $bbhs->{$peg1})                  if ($peg2 = $bbhs->{$peg1})
181                  {                  {
182                        if (! $must_cluster)
183                        {
184                      print TMP "$peg1\t$peg2\n";                      print TMP "$peg1\t$peg2\n";
185                  }                  }
186                  elsif ($conn = $sims{$genomes}->{$peg1})                      else
187                  {                  {
188                      my(@hits);                          my $clustered = &clustered($peg1,$peg2,\%sims,$genome1,$genome2,$bbhs);
189                      my $neigh1 = &neighbors_of($peg1,$neighbors);                          if (@$clustered >= $must_cluster)
   
                     my $tuple;  
                     foreach $tuple (@$conn)  
190                      {                      {
191                          my($iden,$peg2) = @$tuple;                              my $common_subs;
192                          if ($iden >= ($conn->[0]->[0] - $iden_diff))                              if ((! $must_be_in_sub) || (@{($common_subs = &insub($fig,$clustered))} >= $must_be_in_sub))
193                          {                          {
194                              my $neigh2 = &neighbors_of($peg2,$neighbors);                                  print TMP "$peg1\t$peg2\n";
195                              push(@hits,[$peg2,&matches($neigh1,$neigh2,$bbhs)]);                                  print STDERR "connection\t$peg1\t$peg2\t";
196                                    foreach my $x (@$clustered) { print STDERR join(",",@$x),";" }   print STDERR "\t";
197                                    foreach my $x (@$common_subs) { print STDERR join(",",@$x),";" } print STDERR "\n";
198                          }                          }
199                      }                      }
                     @hits = sort { $b->[1] <=> $a->[1] } @hits;  
                     if ((@hits == 1) || ((@hits > 1) && (($hits[0]->[1] - $hits[1]->[1]) > $neigh_diff)))  
                     {  
                         print TMP "$peg1\t$hits[0]->[0]\n";  
200                      }                      }
201                  }                  }
202              }              }
# Line 220  Line 229 
229          print STDERR "$_\n";          print STDERR "$_\n";
230      }      }
231  }  }
232  #unlink($tmpF);  unlink($tmpF);
   
 sub matches {  
     my($n1,$n2,$bbhs) = @_;  
     my($peg,%want);  
233    
234      my $n = 0;  sub clustered {
235      foreach $peg (@$n1)      my($peg1,$peg2,$sims,$genome1,$genome2,$bbhs) = @_;
236        my $neigh1 = &neighbors_of($peg1,$neighbors);
237        my $neigh2 = &neighbors_of($peg2,$neighbors);
238        my %neigh2H = map { $_ => 1 } @$neigh2;
239        my $matches = [];
240        my ($peg_in1,$peg_in2);
241        foreach $peg_in1 (@$neigh1)
242      {      {
243          if ($_ = $bbhs->{$peg})          $peg_in2 = $bbhs->{$peg_in1};
244            if ($neigh2H{$peg_in2})
245          {          {
246              $want{$_} = 1;              push(@$matches,[$peg_in1,$peg_in2]);
247            }
248          }          }
249        return $matches;
250      }      }
251    
252      foreach $peg (@$n2)  sub insub {
253        my($fig,$clustered) = @_;
254    
255        my $n = [];
256        my $tuple;
257        foreach $tuple (@$clustered)
258        {
259            my($peg1,$peg2) = @$tuple;
260            my $func1 = $fig->function_of($peg1);
261            my $func2 = $fig->function_of($peg2);
262            if ($func1 eq $func2)
263      {      {
264          if ($want{$peg})              my @tmp = grep { $fig->usable_subsystem($_) } $fig->peg_to_subsystems($peg1);
265                if (@tmp > 0)
266          {          {
267              $n++;                  @tmp = grep { $fig->usable_subsystem($_) } $fig->peg_to_subsystems($peg2);
268                    if (@tmp > 0)
269                    {
270                        push(@$n,[$peg1,$peg2,$func1]);
271                    }
272                }
273          }          }
274      }      }
275      return $n;      return $n;

Legend:
Removed from v.1.7  
changed lines
  Added in v.1.8

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3