[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.6 - (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 : olson 1.4 use ExpressionDir;
7 : overbeek 1.1
8 : olson 1.4 my($expr_dir);
9 : overbeek 1.5 my $usage = "usage: make_coreg_conjectures_based_on_subsys expr_dir [subsystems-dir]";
10 : overbeek 1.1 (
11 : olson 1.4 ($expr_dir = shift @ARGV)
12 : overbeek 1.1 )
13 :     || die $usage;
14 :    
15 : overbeek 1.5 my %sets;
16 :     if ((@ARGV > 0) && (-s "$ARGV[0]/bindings"))
17 :     {
18 :     foreach $_ (`cat $ARGV[0]/bindings`)
19 :     {
20 :     chop;
21 :     my($ss,$role,$peg) = split(/\t/,$_);
22 :     push(@{$sets{$ss}},$peg);
23 :     }
24 :     }
25 :    
26 : olson 1.4 my $exprO = ExpressionDir->new($expr_dir);
27 :     my $rawF = $exprO->expr_dir . "/rma_normalized.tab";
28 :    
29 : overbeek 1.2 my %bad_subsystems;
30 :    
31 : olson 1.4 #my $sapO = SAPserver->new;
32 :    
33 : overbeek 1.1 my $corrH = &get_corr($rawF);
34 :    
35 : olson 1.4 #my $genomeH = $sapO->genomes_to_subsystems( -ids => [$genome] );
36 :     #my @subs = map { ($_->[1] =~ /^\*?(0|-1)$/) ? () : $_->[0] } @{$genomeH->{$genome}};
37 :     #@subs = grep { ! $bad_subsystems{$_} } @subs;
38 : overbeek 1.2
39 : olson 1.4 #my $subH = $sapO->ids_in_subsystems( -subsystems => \@subs,
40 :     # -genome => $genome);
41 :     my $subH = $exprO->ids_in_subsystems();
42 : overbeek 1.6 my @subs = sort (%sets ? keys(%sets) : keys %$subH);
43 : olson 1.4
44 : overbeek 1.2 my %bad;
45 : olson 1.4
46 : overbeek 1.1 foreach my $sub (@subs)
47 :     {
48 : overbeek 1.2 my %pegs;
49 : overbeek 1.5 my @pegs;
50 : overbeek 1.6 if (%sets)
51 : overbeek 1.5 {
52 :     @pegs = sort { &SeedUtils::by_fig_id($a,$b) }@{$sets{$sub}};
53 :     }
54 :     else
55 : overbeek 1.1 {
56 : overbeek 1.5 my $sub_entry = $subH->{$sub};
57 :     @pegs = ();
58 :     foreach my $role (keys(%$sub_entry))
59 :     {
60 :     my $pegs = $sub_entry->{$role};
61 :     foreach $_ (@$pegs) { $pegs{$_} = 1 }
62 :     }
63 :     @pegs = sort { &SeedUtils::by_fig_id($a,$b) } keys(%pegs);
64 : overbeek 1.2 }
65 :    
66 :     my @sets = grep { @$_ > 1 } split_on_pc(\@pegs,$corrH);
67 :    
68 :     if (@sets > ((@pegs + 2) / 3))
69 :     {
70 :     $bad{$sub} = 1;
71 :     # print STDERR &Dumper([$sub,\@sets]);
72 :     }
73 :     else
74 :     {
75 :     foreach my $set (@sets)
76 :     {
77 :     if (@$set > 1)
78 :     {
79 :     print join(",",@$set),"\tInSubsystem:$sub\n";
80 :     }
81 :     }
82 : overbeek 1.1 }
83 : overbeek 1.2 }
84 :     foreach $_ (keys(%bad))
85 :     {
86 :     print STDERR "bad subsystem\t$_\n";
87 :     }
88 :    
89 :     sub split_on_pc {
90 :     my($pegs,$corrH) = @_;
91 :    
92 :     my @sets = ();
93 :     my %used;
94 :     my $i;
95 :     for ($i=0; ($i < (@$pegs - 1)); $i++)
96 : overbeek 1.1 {
97 : overbeek 1.2 if (! $used{$pegs->[$i]})
98 :     {
99 :     my @poss = ($pegs->[$i]);
100 :     my $j;
101 :     for ($j=$i+1; ($j < @$pegs); $j++)
102 :     {
103 :     if (&corr($pegs->[$j],\@poss,$corrH))
104 :     {
105 :     push(@poss,$pegs->[$j]);
106 :     $used{$pegs->[$j]} = 1;
107 :     }
108 :     }
109 :     push(@sets,\@poss);
110 :     }
111 : overbeek 1.1 }
112 : overbeek 1.2 return @sets;
113 : overbeek 1.1 }
114 :    
115 :     sub get_corr {
116 :     my($rawFF) = @_;
117 :    
118 :     my %gene_to_values;
119 :     open(RAW,"<$rawF") || die "could not open $rawF";
120 :     while (<RAW>)
121 :     {
122 :     chomp;
123 :     my ($gene_id, @gxp_values) = split("\t");
124 :     $gene_to_values{$gene_id} = \@gxp_values;
125 :     }
126 :     close(RAW);
127 :     return \%gene_to_values;
128 :     }
129 :    
130 :     sub compute_pc
131 :     {
132 :     my ($gene_ids, $gxp_hash) = @_;
133 :     my %values = ();
134 :    
135 :     for (my $i = 0; $i < @$gene_ids-1; $i++)
136 :     {
137 :     my $stat = Statistics::Descriptive::Full->new();
138 :     $stat->add_data(@{$gxp_hash->{$gene_ids->[$i]}});
139 :    
140 :     for (my $j = $i+1; $j < @$gene_ids; $j++)
141 :     {
142 :     my ($q, $m, $r, $err) = $stat->least_squares_fit(@{$gxp_hash->{$gene_ids->[$j]}});
143 :     $values{$gene_ids->[$i]}->{$gene_ids->[$j]} = $r;
144 :     $values{$gene_ids->[$j]}->{$gene_ids->[$i]} = $r;
145 :     }
146 :     }
147 :    
148 :     return \%values;
149 :     }
150 :    
151 : overbeek 1.2 sub corr {
152 :     my($peg1,$cluster,$corrH) = @_;
153 :    
154 :     my $hash = &compute_pc([$peg1,@$cluster],$corrH);
155 :    
156 :     # print STDERR &Dumper($peg1,$cluster);
157 :    
158 :     my $sum = 0;
159 :     foreach my $peg2 (@$cluster)
160 :     {
161 :     my $v = $hash->{$peg1}->{$peg2};
162 : overbeek 1.3 if ((! defined($v)) || ($v < 0.4)) { return 0 }
163 : overbeek 1.2 $sum += $v;
164 :     }
165 : overbeek 1.3 return (($sum / @$cluster) >= 0.7);
166 : overbeek 1.2 }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3