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