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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3