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

Diff of /FigKernelScripts/make_coreg_conjectures_based_on_subsys.pl

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

revision 1.1, Tue Apr 20 15:43:06 2010 UTC revision 1.2, Tue Apr 20 18:32:51 2010 UTC
# Line 5  Line 5 
5  use Statistics::Descriptive;  use Statistics::Descriptive;
6    
7  my($genome,$rawF);  my($genome,$rawF);
8  my $usage = "usage: make_coreg_conjectures_based_on_subsys Genome RawDataTab";  my $usage = "usage: make_coreg_conjectures_based_on_subsys Genome RawDataTab BadSubsystems";
9  (  (
10   ($genome    = shift @ARGV) &&   ($genome    = shift @ARGV) &&
11   ($rawF      = shift @ARGV)   ($rawF      = shift @ARGV)
12  )  )
13      || die $usage;      || die $usage;
14    
15    my %bad_subsystems;
16    if (@ARGV > 0)
17    {
18        map { $_ =~ /^(\S.*\S)/; $1 => 1 } `cat $ARGV[0]`;
19    }
20    
21  my $sapO    = SAPserver->new;  my $sapO    = SAPserver->new;
22  my $corrH = &get_corr($rawF);  my $corrH = &get_corr($rawF);
23    
24  my $genomeH = $sapO->genomes_to_subsystems( -ids => [$genome] );  my $genomeH = $sapO->genomes_to_subsystems( -ids => [$genome] );
25  my @subs = map { ($_->[1] =~ /^\*?(0|-1)$/) ? () : $_->[0] } @{$genomeH->{$genome}};  my @subs = map { ($_->[1] =~ /^\*?(0|-1)$/) ? () : $_->[0] } @{$genomeH->{$genome}};
26    @subs = grep { ! $bad_subsystems{$_} } @subs;
27    
28  my $subH = $sapO->ids_in_subsystems( -subsystems => \@subs,  my $subH = $sapO->ids_in_subsystems( -subsystems => \@subs,
29                                       -genome     => $genome);                                       -genome     => $genome);
30    my %bad;
31  foreach my $sub (@subs)  foreach my $sub (@subs)
32  {  {
33        my %pegs;
34      my $sub_entry = $subH->{$sub};      my $sub_entry = $subH->{$sub};
35      my @pegs = ();      my @pegs = ();
36      foreach my $role (keys(%$sub_entry))      foreach my $role (keys(%$sub_entry))
37      {      {
38          my $pegs = $sub_entry->{$role};          my $pegs = $sub_entry->{$role};
39          push(@pegs,@$pegs);          foreach $_ (@$pegs) { $pegs{$_} = 1 }
40        }
41        @pegs = sort { &SeedUtils::by_fig_id($a,$b) } keys(%pegs);
42    
43        my @sets = grep { @$_ > 1 } split_on_pc(\@pegs,$corrH);
44    
45        if (@sets > ((@pegs + 2) / 3))
46        {
47            $bad{$sub} = 1;
48    #       print STDERR &Dumper([$sub,\@sets]);
49      }      }
50      @pegs = sort { &SeedUtils::by_fig_id($a,$b) } @pegs;      else
     if (@pegs > 1)  
51      {      {
52          print join(",",@pegs),"\tInSubsystem:$sub\n";          foreach my $set (@sets)
53            {
54                if (@$set > 1)
55                {
56                    print join(",",@$set),"\tInSubsystem:$sub\n";
57                }
58            }
59      }      }
60  }  }
61    foreach $_ (keys(%bad))
62    {
63        print STDERR "bad subsystem\t$_\n";
64    }
65    
66    sub split_on_pc {
67        my($pegs,$corrH) = @_;
68    
69        my @sets = ();
70        my %used;
71        my $i;
72        for ($i=0; ($i < (@$pegs - 1)); $i++)
73        {
74            if (! $used{$pegs->[$i]})
75            {
76                my @poss = ($pegs->[$i]);
77                my $j;
78                for ($j=$i+1; ($j < @$pegs); $j++)
79                {
80                    if (&corr($pegs->[$j],\@poss,$corrH))
81                    {
82                        push(@poss,$pegs->[$j]);
83                        $used{$pegs->[$j]} = 1;
84                    }
85                }
86                push(@sets,\@poss);
87            }
88        }
89        return @sets;
90    }
91    
92  sub get_corr {  sub get_corr {
93      my($rawFF) = @_;      my($rawFF) = @_;
# Line 71  Line 125 
125      return \%values;      return \%values;
126  }  }
127    
128    sub corr {
129        my($peg1,$cluster,$corrH) = @_;
130    
131        my $hash = &compute_pc([$peg1,@$cluster],$corrH);
132    
133    #   print STDERR &Dumper($peg1,$cluster);
134    
135        my $sum = 0;
136        foreach my $peg2 (@$cluster)
137        {
138            my $v = $hash->{$peg1}->{$peg2};
139            if ((! defined($v)) || ($v < 0.3)) { return 0 }
140            $sum += $v;
141        }
142        return (($sum / @$cluster) >= 0.4);
143    }

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3