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

Diff of /FigKernelScripts/clustered_hypotheticals.pl

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

revision 1.1, Sun Dec 27 01:48:45 2009 UTC revision 1.2, Sun Dec 27 16:06:53 2009 UTC
# Line 7  Line 7 
7  my $sapObject = SAPserver->new();  my $sapObject = SAPserver->new();
8    
9  my $genomeHash = $sapObject->all_genomes(-complete => 1);  my $genomeHash = $sapObject->all_genomes(-complete => 1);
10  for my $genome (sort { $a <=> $b } keys %$genomeHash) {  my @completeG = keys(%$genomeHash);
11    my $isP = $sapObject->is_prokaryotic(-ids => \@completeG);
12    my @complete_proks = grep { $isP->{$_} } @completeG;
13    my %seen;
14    
15    foreach my $genome (sort { $a <=> $b } @complete_proks)
16    {
17      print STDERR "processing $genome\n";      print STDERR "processing $genome\n";
18    
19      my $genomeName = $genomeHash->{$genome};      my $genomeName = $genomeHash->{$genome};
20      my $geneHash = $sapObject->feature_assignments(-genome => $genome,      my $geneHash = $sapObject->feature_assignments(-genome => $genome,
21                                                     -type => 'peg');                                                     -type => 'peg');
22      my @hypotheticalGenes = grep { &SeedUtils::hypo($geneHash->{$_}) } sort { &SeedUtils::by_fig_id($a,$b) } keys(%$geneHash);      my @hypotheticalGenes = sort { &SeedUtils::by_fig_id($a,$b) }
23                                grep { &SeedUtils::hypo($geneHash->{$_}) }
24                                grep { ! $seen{$_} }
25                                keys(%$geneHash);
26    
27      my $couplingHash = $sapObject->conserved_in_neighborhood(-ids => \@hypotheticalGenes);      my $couplingHash = $sapObject->conserved_in_neighborhood(-ids => \@hypotheticalGenes);
28        my @all_coupled_ids = map { map { $_->[1] } @{$couplingHash->{$_}}} keys(%$couplingHash);
29        my $subHash         = $sapObject->ids_to_subsystems(-ids => \@all_coupled_ids);
30    
31      for my $gene (@hypotheticalGenes) {      foreach my $gene (@hypotheticalGenes)
32        {
33          my $couplingList = $couplingHash->{$gene};          my $couplingList = $couplingHash->{$gene};
34            if (defined $couplingList)
35            {
36                my ($bestCoupler, $bestScore, $bestFunction) = (undef, 0, '');
37                my $best_pairset;
38    
39          if (defined $couplingList) {              foreach  my $coupling (@$couplingList)
40                {
41                    my ($score, $coupler, $function, $pairset) = @$coupling;
42    
43              my $subHash = $sapObject->ids_to_subsystems(-ids => [ map { $_->[1]} @$couplingList ]);                  if ($subHash->{$coupler} && $score > $bestScore)
44              my ($bestCoupler, $bestScore, $bestFunction) = (undef, 0, '');                  {
             my %seen;  
             for my $coupling (@$couplingList) {  
                 my ($score, $coupler, $function, $pid) = @$coupling;  
                 if ( ! $seen{$pid}) {  
                         $seen{$pid} = 1;  
                         if ($subHash->{$coupler} && $score > $bestScore) {  
45                              $bestCoupler = $coupler;                              $bestCoupler = $coupler;
46                              $bestScore = $score;                              $bestScore = $score;
47                              $bestFunction = $function;                      $bestFunction = $function ? $function : '';
48                        $best_pairset = $pairset;
49                          }                          }
50                  }                  }
51              }  
52              if (defined $bestCoupler) {              if (defined $bestCoupler)
53                {
54                  print join("\t", $gene, $geneHash->{$gene}, $genome, $genomeName,                  print join("\t", $gene, $geneHash->{$gene}, $genome, $genomeName,
55                                   $bestCoupler, $bestScore, $bestFunction) . "\n";                                   $bestCoupler, $bestScore, $bestFunction) . "\n";
56                    &set_seen_for_pairset($gene,$best_pairset,\%seen,$sapObject);
57              }              }
58          }          }
59      }      }
60  }  }
61    
62    sub set_seen_for_pairset {
63        my($peg,$pairset,$seen,$sapObject) = @_;
64    
65        my $pairsetsH = $sapObject->pairsets( -ids => [$pairset] );
66        if ($_ = $pairsetsH->{$pairset})
67        {
68            my($sc,$pairs) = @$_;
69            my $got = 0;
70            my @set =   map { if ($_->[0] eq $peg) { $got = 1 }; $_->[0] } @$pairs;
71            if (! $got)
72            {
73                @set = map { $_->[1] } @$pairs;
74            }
75            foreach $_ (@set)
76            {
77                $seen->{$_} = 1;
78            }
79        }
80    }
81    
82    

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3