annotate splicesitescore/score5.pl @ 8:ec97befe722b

Uploaded
author xuebing
date Sun, 01 Apr 2012 08:14:45 -0400
parents 3fd430270c0c
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
22486256ed0d Uploaded
xuebing
parents:
diff changeset
1 use strict;
22486256ed0d Uploaded
xuebing
parents:
diff changeset
2
22486256ed0d Uploaded
xuebing
parents:
diff changeset
3
22486256ed0d Uploaded
xuebing
parents:
diff changeset
4 my $inputfile = $ARGV[0];
22486256ed0d Uploaded
xuebing
parents:
diff changeset
5 my $usemaxent = 1;
22486256ed0d Uploaded
xuebing
parents:
diff changeset
6
8
ec97befe722b Uploaded
xuebing
parents: 7
diff changeset
7 my $modelpath = "/home/wuxbl/rowley/tools/maxentscan/";
0
22486256ed0d Uploaded
xuebing
parents:
diff changeset
8 my %me2x5 = &makescorematrix($modelpath.'me2x5');
22486256ed0d Uploaded
xuebing
parents:
diff changeset
9 my %seq = &makesequencematrix($modelpath.'splicemodels/splice5sequences');
22486256ed0d Uploaded
xuebing
parents:
diff changeset
10
22486256ed0d Uploaded
xuebing
parents:
diff changeset
11 my %bgd;
22486256ed0d Uploaded
xuebing
parents:
diff changeset
12 $bgd{'A'} = 0.27;
22486256ed0d Uploaded
xuebing
parents:
diff changeset
13 $bgd{'C'} = 0.23;
22486256ed0d Uploaded
xuebing
parents:
diff changeset
14 $bgd{'G'} = 0.23;
22486256ed0d Uploaded
xuebing
parents:
diff changeset
15 $bgd{'T'} = 0.27;
22486256ed0d Uploaded
xuebing
parents:
diff changeset
16
22486256ed0d Uploaded
xuebing
parents:
diff changeset
17
22486256ed0d Uploaded
xuebing
parents:
diff changeset
18
22486256ed0d Uploaded
xuebing
parents:
diff changeset
19 open (FILE,"<$inputfile") || die "can't open!\n";
22486256ed0d Uploaded
xuebing
parents:
diff changeset
20
22486256ed0d Uploaded
xuebing
parents:
diff changeset
21 while(<FILE>) {
22486256ed0d Uploaded
xuebing
parents:
diff changeset
22 chomp;
22486256ed0d Uploaded
xuebing
parents:
diff changeset
23 if (/^\s*$/) { #discard blank lines;
22486256ed0d Uploaded
xuebing
parents:
diff changeset
24 next;
22486256ed0d Uploaded
xuebing
parents:
diff changeset
25 }
22486256ed0d Uploaded
xuebing
parents:
diff changeset
26 elsif (/^>/) { #discard comment lines;
22486256ed0d Uploaded
xuebing
parents:
diff changeset
27 print $_."\t";
22486256ed0d Uploaded
xuebing
parents:
diff changeset
28 next;
22486256ed0d Uploaded
xuebing
parents:
diff changeset
29 }
22486256ed0d Uploaded
xuebing
parents:
diff changeset
30 elsif (/[NQWERYUIOPLKJHFDSZXVBM]/) {
22486256ed0d Uploaded
xuebing
parents:
diff changeset
31 next;
22486256ed0d Uploaded
xuebing
parents:
diff changeset
32 }
22486256ed0d Uploaded
xuebing
parents:
diff changeset
33 else {
22486256ed0d Uploaded
xuebing
parents:
diff changeset
34 $_ =~ s/\cM//g; #gets rid of carriage return
22486256ed0d Uploaded
xuebing
parents:
diff changeset
35 my $str = $_;
22486256ed0d Uploaded
xuebing
parents:
diff changeset
36 print $str."\t";
22486256ed0d Uploaded
xuebing
parents:
diff changeset
37 $str = uc($str);
22486256ed0d Uploaded
xuebing
parents:
diff changeset
38 if ($usemaxent) {
22486256ed0d Uploaded
xuebing
parents:
diff changeset
39 print sprintf("%.2f",&log2(&scoreconsensus($str)*$me2x5{$seq{&getrest($str)}}))."\n";
22486256ed0d Uploaded
xuebing
parents:
diff changeset
40 }
22486256ed0d Uploaded
xuebing
parents:
diff changeset
41 }
22486256ed0d Uploaded
xuebing
parents:
diff changeset
42 }
22486256ed0d Uploaded
xuebing
parents:
diff changeset
43
22486256ed0d Uploaded
xuebing
parents:
diff changeset
44
22486256ed0d Uploaded
xuebing
parents:
diff changeset
45 sub makesequencematrix{
22486256ed0d Uploaded
xuebing
parents:
diff changeset
46 my $file = shift;
22486256ed0d Uploaded
xuebing
parents:
diff changeset
47 my %matrix;my $n=0;
22486256ed0d Uploaded
xuebing
parents:
diff changeset
48 open(SCOREF, $file) || die "Can't open $file!\n";
22486256ed0d Uploaded
xuebing
parents:
diff changeset
49 while(<SCOREF>) {
22486256ed0d Uploaded
xuebing
parents:
diff changeset
50 chomp;
22486256ed0d Uploaded
xuebing
parents:
diff changeset
51 $_=~ s/\s//;
22486256ed0d Uploaded
xuebing
parents:
diff changeset
52 $matrix{$_} = $n;
22486256ed0d Uploaded
xuebing
parents:
diff changeset
53 $n++;
22486256ed0d Uploaded
xuebing
parents:
diff changeset
54 }
22486256ed0d Uploaded
xuebing
parents:
diff changeset
55 close(SCOREF);
22486256ed0d Uploaded
xuebing
parents:
diff changeset
56 return %matrix;
22486256ed0d Uploaded
xuebing
parents:
diff changeset
57 }
22486256ed0d Uploaded
xuebing
parents:
diff changeset
58 sub makescorematrix{
22486256ed0d Uploaded
xuebing
parents:
diff changeset
59 my $file = shift;
22486256ed0d Uploaded
xuebing
parents:
diff changeset
60 my %matrix;my $n=0;
22486256ed0d Uploaded
xuebing
parents:
diff changeset
61 open(SCOREF, $file) || die "Can't open $file!\n";
22486256ed0d Uploaded
xuebing
parents:
diff changeset
62 while(<SCOREF>) {
22486256ed0d Uploaded
xuebing
parents:
diff changeset
63 chomp;
22486256ed0d Uploaded
xuebing
parents:
diff changeset
64 $_=~ s/\s//;
22486256ed0d Uploaded
xuebing
parents:
diff changeset
65 $matrix{$n} = $_;
22486256ed0d Uploaded
xuebing
parents:
diff changeset
66 $n++;
22486256ed0d Uploaded
xuebing
parents:
diff changeset
67 }
22486256ed0d Uploaded
xuebing
parents:
diff changeset
68 close(SCOREF);
22486256ed0d Uploaded
xuebing
parents:
diff changeset
69 return %matrix;
22486256ed0d Uploaded
xuebing
parents:
diff changeset
70 }
22486256ed0d Uploaded
xuebing
parents:
diff changeset
71
22486256ed0d Uploaded
xuebing
parents:
diff changeset
72 sub getrest{
22486256ed0d Uploaded
xuebing
parents:
diff changeset
73 my $seq = shift;
22486256ed0d Uploaded
xuebing
parents:
diff changeset
74 my @seqa = split(//,uc($seq));
22486256ed0d Uploaded
xuebing
parents:
diff changeset
75 return $seqa[0].$seqa[1].$seqa[2].$seqa[5].$seqa[6].$seqa[7].$seqa[8];
22486256ed0d Uploaded
xuebing
parents:
diff changeset
76 }
22486256ed0d Uploaded
xuebing
parents:
diff changeset
77 sub scoreconsensus{
22486256ed0d Uploaded
xuebing
parents:
diff changeset
78 my $seq = shift;
22486256ed0d Uploaded
xuebing
parents:
diff changeset
79 my @seqa = split(//,uc($seq));
22486256ed0d Uploaded
xuebing
parents:
diff changeset
80 my %bgd;
22486256ed0d Uploaded
xuebing
parents:
diff changeset
81 $bgd{'A'} = 0.27;
22486256ed0d Uploaded
xuebing
parents:
diff changeset
82 $bgd{'C'} = 0.23;
22486256ed0d Uploaded
xuebing
parents:
diff changeset
83 $bgd{'G'} = 0.23;
22486256ed0d Uploaded
xuebing
parents:
diff changeset
84 $bgd{'T'} = 0.27;
22486256ed0d Uploaded
xuebing
parents:
diff changeset
85 my %cons1;
22486256ed0d Uploaded
xuebing
parents:
diff changeset
86 $cons1{'A'} = 0.004;
22486256ed0d Uploaded
xuebing
parents:
diff changeset
87 $cons1{'C'} = 0.0032;
22486256ed0d Uploaded
xuebing
parents:
diff changeset
88 $cons1{'G'} = 0.9896;
22486256ed0d Uploaded
xuebing
parents:
diff changeset
89 $cons1{'T'} = 0.0032;
22486256ed0d Uploaded
xuebing
parents:
diff changeset
90 my %cons2;
22486256ed0d Uploaded
xuebing
parents:
diff changeset
91 $cons2{'A'} = 0.0034;
22486256ed0d Uploaded
xuebing
parents:
diff changeset
92 $cons2{'C'} = 0.0039;
22486256ed0d Uploaded
xuebing
parents:
diff changeset
93 $cons2{'G'} = 0.0042;
22486256ed0d Uploaded
xuebing
parents:
diff changeset
94 $cons2{'T'} = 0.9884;
22486256ed0d Uploaded
xuebing
parents:
diff changeset
95 my $addscore = $cons1{$seqa[3]}*$cons2{$seqa[4]}/($bgd{$seqa[3]}*$bgd{$seqa[4]});
22486256ed0d Uploaded
xuebing
parents:
diff changeset
96 return $addscore;
22486256ed0d Uploaded
xuebing
parents:
diff changeset
97 }
22486256ed0d Uploaded
xuebing
parents:
diff changeset
98
22486256ed0d Uploaded
xuebing
parents:
diff changeset
99 sub log2{
22486256ed0d Uploaded
xuebing
parents:
diff changeset
100 my ($val) = @_;
22486256ed0d Uploaded
xuebing
parents:
diff changeset
101 return log($val)/log(2);
22486256ed0d Uploaded
xuebing
parents:
diff changeset
102 }