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