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

Diff of /FigKernelScripts/split_sequences_into_sets.pl

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

revision 1.2, Wed Aug 19 16:11:28 2009 UTC revision 1.3, Sun Aug 23 18:40:43 2009 UTC
# Line 53  Line 53 
53  }  }
54    
55  my $min_size = (@seqs > 20) ? (@seqs / 10) : (@seqs / 2);  my $min_size = (@seqs > 20) ? (@seqs / 10) : (@seqs / 2);
56  my @all = &ProtSims::blastP(\@seqs,\@seqs,$min_size);  my @all = map { my($id1,$id2,undef,undef,undef,undef,$b1,$e1,$b2,$e2,$psc,undef,$ln1,$ln2) = @$_;
57  # &print_sims(\@all);                  [$id1,$id2,(($e1+1)-$b1)/$ln1,(($e2+1)-$b2)/$ln2,$psc]
58                  } &ProtSims::blastP(\@seqs,\@seqs,$min_size);
59    
 # $_ = @all; print STDERR "computed $_ similarities\n";  
60    
61  my @sims = grep { ($_->id1 ne $_->id2) &&  my @sims = grep { my($id1,$id2,$frac1,$frac2,$psc) = @$_;
62                    ($_->psc <= $sim_cutoff) &&                    ($id1 ne $id2) && ($frac1 >= $frac_cov) && ($frac2 >= $frac_cov) } @all;
                   ((($_->e1 - $_->b1) / $_->ln1) >= $frac_cov) &&  
                   ((($_->e2 - $_->b2) / $_->ln2) >= $frac_cov)  
                 } @all;  
   
 # print STDERR "=======\n";  
 # &print_sims(\@sims);  
63    
64  my(%seen,%conn);  my(%seen,%conn);
65  foreach my $sim (@sims)  foreach my $sim (@sims)
66  {  {
67      if (! $seen{"$sim->id1\t$sim->id2"})  # pick only the best similarity between two ids      my($id1,$id2,$frac1,$frac2,$psc) = @$sim;
68      {      my $key = join("\t",sort ($id1,$id2));
69          push(@{$conn{$sim->id1}},$sim);      if (! $seen{$key})
70          $seen{"$sim->id1\t$sim->id2"} = 1;      {
71            $seen{$key} = 1;
72            push(@{$conn{$id1}},$id2);
73            push(@{$conn{$id2}},$id1);
74      }      }
75  }  }
76  undef %seen;  undef %seen;
77  # print STDERR "built connections\n";  my @ids = sort { my $x = $conn{$b} ? @{$conn{$b}} : 0;
78                     my $y = $conn{$a} ? @{$conn{$a}} : 0;
79  my @ids = sort { &cov($conn{$b}) <=> &cov($conn{$a}) }                   ($x <=> $y)
80            keys(%seq_of);                 } keys(%seq_of);
 # sort by sum of coverage  
81    
82  my $n = 1;  my $n = 1;
83  my @sets = ();  my @sets = ();
# Line 109  Line 105 
105      }      }
106  }  }
107  close(SZ);  close(SZ);
 # print STDERR "built set sizes\n";  
108    
109  my @distances = &construct_distance_matrix(\@sets,\@all); # entries are [n1,n2,sc]  my @distances = &construct_distance_matrix(\@sets,\@all); # entries are [n1,n2,sc]
110  open(DIST,">$output/distance.matrix\n") || die "could not open $output/distance.matrix";  open(DIST,">$output/distance.matrix\n") || die "could not open $output/distance.matrix";
# Line 118  Line 113 
113      print DIST join("\t",@$tuple),"\n";      print DIST join("\t",@$tuple),"\n";
114  }  }
115  close(DIST);  close(DIST);
 # print STDERR "built distance matrix\n";  
   
 sub print_sims {  
     my($sims) = @_;  
   
     foreach my $sim (@$sims)  
     {  
         print STDERR join("\t",($sim->id1,$sim->id2,$sim->iden,$sim->b1,$sim->e1,$sim->b2,$sim->e2,$sim->psc,$sim->ln1,$sim->ln2)),"\n";  
     }  
 }  
   
 sub cov {  
     my($sims) = @_;  
     my($x,$n);  
   
     $n = 0;  
     foreach $x (@$sims)  
     {  
         $n += ($x->e1 - $x->b1);  
     }  
     return $n;  
 }  
116    
117  sub run {  sub run {
118      my($cmd) = @_;      my($cmd) = @_;
# Line 159  Line 132 
132      {      {
133          my $x = $conn->{$cluster->[$i]};          my $x = $conn->{$cluster->[$i]};
134          $i++;          $i++;
135          foreach my $peg (map { $_->id2} @$x)          foreach my $peg (@$x)
136          {          {
137              if (! $in{$peg})              if (! $in{$peg})
138              {              {
# Line 188  Line 161 
161      my %best;      my %best;
162      foreach my $sim (@$sims)      foreach my $sim (@$sims)
163      {      {
164          my $id1 = $sim->id1;          my($id1,$id2,undef,undef,$sc) = @$sim;
         my $id2 = $sim->id2;  
         my $sc = $sim->psc;  
165          if ($in{$id1} && $in{$id2} && ($in{$id1} != $in{$id2}))          if ($in{$id1} && $in{$id2} && ($in{$id1} != $in{$id2}))
166          {          {
167              my $key = join("\t",($in{$id1},$in{$id2}));              my $key = join("\t",($in{$id1},$in{$id2}));
168              if ((! ($sofar = $best{$key})) || ($sofar > $sc))              if ((! ($sofar = $best{$key})) || ($sofar < $sc))
169              {              {
170                  $best{$key} = $sc;                  $best{$key} = $sc;
171              }              }

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3