annotate maxent_score5.pl @ 12:6e7e036c13ed

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