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

Annotation of /FigKernelScripts/bayes_expected_start_scores.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : overbeek 1.1 # -*- perl -*-
2 :    
3 :     use FIG;
4 :    
5 :     use Memoize;
6 :     memoize('binom');
7 :    
8 :     $0 =~ m/^([^\/]+)$/;
9 :     $usage = "$1 avg_peg_len avg_leader_len < scored.start.objects > re-scored.start.objects";
10 :    
11 :     ((($avg_peg_len, $avg_leader_len) = @ARGV) == 2)
12 :     || die "usage: $usage";
13 :    
14 :     $nu = 2;
15 :     $p_peg = $nu / ($nu + $avg_peg_len/3);
16 :     $p_leader = 1 / (1 + $avg_leader_len/3);
17 :    
18 :     $short_call_weight = -0.50;
19 :     $long_call_weight = -0.25;
20 :    
21 :     $/ = "///\n";
22 :    
23 :     $tick = $tock = 0;
24 :     while (defined($entry = <STDIN>))
25 :     {
26 :     chomp $entry;
27 :    
28 :     if ($tick >= 100) { print STDERR "."; $tick = 0; $tock++; } else { $tick++; }
29 :     if ($tock >= 50) { print STDERR "\n"; $tock = 0; }
30 :    
31 :     @calls = split m{//\n}, $entry;
32 :     $header = shift @calls;
33 :     unless ($calls[-1]) { pop @calls; }
34 :    
35 :     %header_fields = map { m/^([^=]+)=([^=]+)$/ } split /\n/, $header;
36 :     $orf_beg = $header_fields{BEGIN};
37 :     $orf_end = $header_fields{END};
38 :     $orf_len = (1 + abs($orf_end - $orf_beg)) / 3;
39 :    
40 :     undef %scores;
41 :     foreach $call (@calls)
42 :     {
43 :     $fields = {};
44 :     %$fields = map { m/^([^=]+)=([^=]+)$/; $1 => $2 } split /\n/, $call;
45 :    
46 :     $calls{$fields->{START}} = $fields;
47 :     $scores{$fields->{START}} = $fields->{SCORE};
48 :     }
49 :     @candidates = sort { $a <=> $b } keys %scores;
50 :    
51 :     $normalization = 0;
52 :     foreach $candidate (@candidates)
53 :     {
54 :     $called_len = $orf_len - ($candidate-1)/3;
55 :    
56 :     $prob_leader = $p_leader * (1-$p_leader)**(($candidate-1)/3);
57 :     $prob_peg = &binom(($called_len + $nu - 1), $called_len)
58 :     * $p_peg**$nu * (1-$p_peg)**$called_len;
59 :    
60 :     $unnormed_prob{$candidate} = $prob_leader * $prob_peg;
61 :    
62 :     $normalization += $unnormed_prob{$candidate};
63 :     }
64 :    
65 :     $best_pos = 0;
66 :     $best_expect = -9.99999999e99;
67 :     foreach $called_pos (@candidates)
68 :     {
69 :     $expectation = 0;
70 :     foreach $candidate (@candidates)
71 :     {
72 :     $unnormed_prob = $unnormed_prob{$candidate};
73 :    
74 :     if ($called_pos < $candidate)
75 :     {
76 :     $expectation += $unnormed_prob * $long_call_weight * $scores{$candidate};
77 :     }
78 :     elsif ($called_pos > $candidate)
79 :     {
80 :     $expectation += $unnormed_prob * $short_call_weight * $scores{$candidate};
81 :     }
82 :     else # ($called_pos == $candidate)
83 :     {
84 :     $expectation += $unnormed_prob * $scores{$candidate};
85 :     }
86 :     }
87 :     $expectation /= $normalization;
88 :    
89 :     if ($expectation > $best_expect)
90 :     {
91 :     $best_pos = $called_pos;
92 :     $best_expect = $expectation;
93 :     }
94 :    
95 :     $calls{$called_pos}->{RAW_SCORE} = $calls{$called_pos}->{SCORE};
96 :     $calls{$called_pos}->{SCORE} = $expectation;
97 :     }
98 :    
99 :     if ($orf_beg < $orf_end)
100 :     {
101 :     $new_start = $orf_beg + $best_pos - 1;
102 :     }
103 :     else
104 :     {
105 :     $new_start = $orf_beg - $best_pos + 1;
106 :     }
107 :    
108 :     $header_fields{NEW_START_POS} = $best_pos;
109 :     $header_fields{NEW_START_LOC} = "$header_fields{CONTIG}\_$new_start\_$orf_end";
110 :    
111 :     while (($k, $v) = each %header_fields) { print "$k=$v\n"; }
112 :     print "//\n";
113 :     foreach $candidate (sort { $calls{$b}->{SCORE} <=> $calls{$a}->{SCORE} } keys %calls)
114 :     {
115 :     $call = $calls{$candidate};
116 :     while (($k, $v) = each %$call) { print "$k=$v\n"; }
117 :     print "//\n";
118 :     }
119 :     print "///\n";
120 :     }
121 :     print STDERR "\n\n";
122 :    
123 :    
124 :     sub binom
125 :     {
126 :     my ($n, $k) = @_;
127 :     my ($sign, $binom);
128 :    
129 :     $sign = 1;
130 :     if ($n < 0) { $n = ($k - $n - 1); $sign = -1 if ($k % 2); }
131 :     if ($k > $n/2) { $k = ($n - $k); }
132 :     # print STDERR "n = $n, k = $k\n";
133 :    
134 :     return 0 if ($k < 0);
135 :     return 1 if ($k == 0);
136 :    
137 :     $binom = &binom($n-1, $k) + &binom($n-1, $k-1);
138 :    
139 :     return ($sign * $binom);
140 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3