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

Annotation of /FigKernelScripts/ex_make_on_off_calls.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : overbeek 1.1 use Data::Dumper;
2 :     use File::Basename;
3 :     use Getopt::Long;
4 :     use SeedEnv;
5 :     use strict;
6 :    
7 :     my $usage = "usage: ex_make_on_off_calls -d DataDir\n";
8 :    
9 :     my $dataD;
10 :     my $rc = GetOptions('d=s' => \$dataD);
11 :    
12 :     if ((! $rc) || (! -d $dataD))
13 :     {
14 :     print STDERR $usage; exit ;
15 :     }
16 :    
17 :     my $initial_calls = &make_initial_calls($dataD);
18 :     open(CALLS,">$dataD/peg.on.off.calls")
19 :     || die " could not open $dataD/initial.peg.on.off.calls";
20 :     foreach my $exp (sort keys(%$initial_calls))
21 :     {
22 :     my $locusH = $initial_calls->{$exp};
23 :     foreach my $locus (keys(%$locusH))
24 :     {
25 :     my $v = $locusH->{$locus};
26 :     print CALLS join("\t",($exp,$locus,$v)),"\n";
27 :     }
28 :     }
29 :     close(CALLS);
30 :    
31 :     sub make_initial_calls {
32 :     my($dataD) = @_;
33 :    
34 :     my %roles_always_on = map { chop; ($_ => 1) } `cat $dataD/roles.always.on`;
35 :     my %pegs_always_on = map { chop;
36 : overbeek 1.2 my(undef,$locus,undef,$func) = split(/\t/,$_);
37 : overbeek 1.1 $roles_always_on{$func} ? ($locus => 1) : () }
38 :     `cat $dataD/aliases.with.function`;
39 :     my %exp_peg_val;
40 :     open(EPV,"<$dataD/gene.locus.exp.val.tuples")
41 :     || die "could not open $dataD/gene.locus.exp.val.tuples";
42 :     while (defined($_ = <EPV>))
43 :     {
44 :     chop;
45 :     my($gene,$locus,$exp,$v) = split(/\t/,$_);
46 :     $exp_peg_val{$exp}->{$locus} = $v;
47 :     }
48 :     close(EPV);
49 :     my %on_for_exp;
50 :     my %off_for_exp;
51 :     my %diff_by_exp;
52 :     foreach my $exp (keys(%exp_peg_val))
53 :     {
54 :     my $locusH = $exp_peg_val{$exp};
55 :     my @on_values = sort { $a <=> $b }
56 :     map { $pegs_always_on{$_} ? $locusH->{$_} : () } keys(%$locusH);
57 :     my $I_10 = int(@on_values * 0.1);
58 :     my $on_threshold = $on_values[$I_10];
59 :     $on_for_exp{$exp} = $on_threshold;
60 :     my @not_on_values = sort { $a <=> $b }
61 :     map { ((! $pegs_always_on{$_} ) && ($locusH->{$_} < $on_threshold)) ? $locusH->{$_} : () }
62 :     keys(%$locusH);
63 :     my $I_80 = int(@not_on_values * 0.8);
64 :     my $off_threshold = $not_on_values[$I_80];
65 :     $off_for_exp{$exp} = $off_threshold;
66 :     $diff_by_exp{$exp} = $on_threshold - $off_threshold;
67 :     }
68 :     my @diffs = sort { $a <=> $b } map { $diff_by_exp{$_} } keys(%diff_by_exp);
69 :     my $I_25 = int(@diffs * 0.25);
70 :     my $diff_threshold = $diffs[$I_25];
71 :     foreach my $exp (keys(%off_for_exp))
72 :     {
73 :     if (($on_for_exp{$exp} - $off_for_exp{$exp}) < $diff_threshold)
74 :     {
75 :     $off_for_exp{$exp} = $on_for_exp{$exp} - $diff_threshold;
76 :     }
77 :     }
78 :     my %peg_on_off_by_exp;
79 :     foreach my $exp (keys(%exp_peg_val))
80 :     {
81 :     my $on_threshold = $on_for_exp{$exp};
82 :     my $off_threshold = $off_for_exp{$exp};
83 :    
84 :     my $locusH = $exp_peg_val{$exp};
85 :     foreach my $peg (keys(%$locusH))
86 :     {
87 :     my $v = $locusH->{$peg};
88 :     if ($v >= $on_threshold)
89 :     {
90 :     $peg_on_off_by_exp{$exp}->{$peg} = 1;
91 :     }
92 :     elsif ($v <= $off_threshold)
93 :     {
94 :     $peg_on_off_by_exp{$exp}->{$peg} = -1;
95 :     }
96 :     else
97 :     {
98 :     $peg_on_off_by_exp{$exp}->{$peg} = 0;
99 :     }
100 :     }
101 :     }
102 :     return \%peg_on_off_by_exp;
103 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3