gistfile1.pl

use 5.12.4;
use Storable;
use Data::Dump qw/dump/;
use G;

my $f_energy = retrieve('./energy.store');
my $gb = load ('ecoli', 'no msg');
my $ribosome;

open my $RIBO, '<' ,$ARGV[0] || die;
while (<$RIBO>) {
        s/[\r\n]+\z//;
        my ($pos, $r_value) = split /\t/;
        $ribosome->{direct}->[$pos] = $r_value;
}

my @count;

for my $cds ( $gb->feature('CDS') ) {
        next unless $gb->{$cds}->{direction} eq 'direct';
        my $start = $gb->{$cds}->{start};
        my $end = $gb->{$cds}->{end};

        my @array;
        for my $i ($start .. $end){
                push @array, $ribosome->{direct}->[$i] if defined $ribosome->{direct}->[$i];
        }
        my $mean = mean @array;
        for my $i ($start .. $end){
                if (defined $ribosome->{direct}->[$i] && $ribosome->{direct}->[$i] >= (10 * $mean)) {
                        for my $j (0 .. 19) {
                                my @peakpos = grep {$_ <= -4 } @{$f_energy->{direct}}[$i - $j .. $i];
                                $count[$j]++ if scalar @peakpos;
                        }
                }
        }
}

for my $i (0 .. $#count) {
        say join "\t", ($i,$count[$i]);
}