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

Annotation of /FigKernelScripts/make_coreg_conjectures_based_on_subsys.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : overbeek 1.1 use SeedEnv;
2 :     use strict;
3 :     use Carp;
4 :     use Data::Dumper;
5 :     use Statistics::Descriptive;
6 :    
7 :     my($genome,$rawF);
8 : overbeek 1.2 my $usage = "usage: make_coreg_conjectures_based_on_subsys Genome RawDataTab BadSubsystems";
9 : overbeek 1.1 (
10 :     ($genome = shift @ARGV) &&
11 :     ($rawF = shift @ARGV)
12 :     )
13 :     || die $usage;
14 :    
15 : overbeek 1.2 my %bad_subsystems;
16 :     if (@ARGV > 0)
17 :     {
18 :     map { $_ =~ /^(\S.*\S)/; $1 => 1 } `cat $ARGV[0]`;
19 :     }
20 :    
21 : overbeek 1.1 my $sapO = SAPserver->new;
22 :     my $corrH = &get_corr($rawF);
23 :    
24 :     my $genomeH = $sapO->genomes_to_subsystems( -ids => [$genome] );
25 :     my @subs = map { ($_->[1] =~ /^\*?(0|-1)$/) ? () : $_->[0] } @{$genomeH->{$genome}};
26 : overbeek 1.2 @subs = grep { ! $bad_subsystems{$_} } @subs;
27 :    
28 : overbeek 1.1 my $subH = $sapO->ids_in_subsystems( -subsystems => \@subs,
29 :     -genome => $genome);
30 : overbeek 1.2 my %bad;
31 : overbeek 1.1 foreach my $sub (@subs)
32 :     {
33 : overbeek 1.2 my %pegs;
34 : overbeek 1.1 my $sub_entry = $subH->{$sub};
35 :     my @pegs = ();
36 :     foreach my $role (keys(%$sub_entry))
37 :     {
38 :     my $pegs = $sub_entry->{$role};
39 : overbeek 1.2 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 :     else
51 :     {
52 :     foreach my $set (@sets)
53 :     {
54 :     if (@$set > 1)
55 :     {
56 :     print join(",",@$set),"\tInSubsystem:$sub\n";
57 :     }
58 :     }
59 : overbeek 1.1 }
60 : overbeek 1.2 }
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 : overbeek 1.1 {
74 : overbeek 1.2 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 : overbeek 1.1 }
89 : overbeek 1.2 return @sets;
90 : overbeek 1.1 }
91 :    
92 :     sub get_corr {
93 :     my($rawFF) = @_;
94 :    
95 :     my %gene_to_values;
96 :     open(RAW,"<$rawF") || die "could not open $rawF";
97 :     while (<RAW>)
98 :     {
99 :     chomp;
100 :     my ($gene_id, @gxp_values) = split("\t");
101 :     $gene_to_values{$gene_id} = \@gxp_values;
102 :     }
103 :     close(RAW);
104 :     return \%gene_to_values;
105 :     }
106 :    
107 :     sub compute_pc
108 :     {
109 :     my ($gene_ids, $gxp_hash) = @_;
110 :     my %values = ();
111 :    
112 :     for (my $i = 0; $i < @$gene_ids-1; $i++)
113 :     {
114 :     my $stat = Statistics::Descriptive::Full->new();
115 :     $stat->add_data(@{$gxp_hash->{$gene_ids->[$i]}});
116 :    
117 :     for (my $j = $i+1; $j < @$gene_ids; $j++)
118 :     {
119 :     my ($q, $m, $r, $err) = $stat->least_squares_fit(@{$gxp_hash->{$gene_ids->[$j]}});
120 :     $values{$gene_ids->[$i]}->{$gene_ids->[$j]} = $r;
121 :     $values{$gene_ids->[$j]}->{$gene_ids->[$i]} = $r;
122 :     }
123 :     }
124 :    
125 :     return \%values;
126 :     }
127 :    
128 : overbeek 1.2 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 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3