Because it's weekend. It is fully functional. I hope you can follow the code, otherwise I'm happy to explain... (Fasta parser is based on Fasta-Parser.pl)
#!/usr/bin/env perl
use warnings;
use strict;
# Usage:
# perl dinuc.pl FILE.FA
my @dipeps_of_interest = qw(
GA
AL
MM
DE
DV
VD
DW
QD
SD
DD
ED
DY
VE
EN
II
KE
NV
VP
FV
SS
WK
KK
);
open(FASTA, $ARGV[0]) or die $!;
while( my $fa = next_fasta_record(\*FASTA) ){
my %fa = fasta_record2hash($fa);
my %dipeps = count_dipeps($fa{seq});
# or
# my %dipeps = count_dipeps_fast($fa{seq});
# report dipeps counts
print "$fa{id}\n";
my $sum = 0;
my $abs = 0;
foreach (@dipeps_of_interest) {
if (exists $dipeps{$_}) {
$sum+= $dipeps{$_};
print "$_\t$dipeps{$_}\n";
}else {
$abs++;
print "$_\t0\n";
}
}
print "sum\t$sum\n";
print "abs\t$abs\n";
}
close FASTA;
# count dipeps of sequence in a hash
sub count_dipeps{
my $seq = shift;
my %dipeps;
while($seq =~ /(?=(\w\w))/g){ # read to chars at a time
$dipeps{$1}++;
}
return %dipeps;
}
# count dipeps of sequence in a hash
# faster than regex and works on seqs >36kbp
sub count_dipeps_fast{
my $seq = shift;
my %dipeps;
my $length = length($seq)-2;
foreach(unpack("(A2X1)".$length."A2", $seq)){
$dipeps{$_}++
}
return %dipeps;
}
# read a single fasta record from a file handle into a string
sub next_fasta_record{
my $fh = shift;
local $/="\n>"; # change record separator form "\n" to "\n>"
my $fa = <$fh>; # read a fasta record
return unless defined($fa); # end of file
chomp($fa); # remove "\n>" from end of record
$fa =~ s/^>//; # (stupid way to) fix first record in file
return '>'.$fa;
}
# split a single fasta record into a hash
sub fasta_record2hash{
my ($head, $seq) = split("\n", $_[0], 2); # split header and seq
my ($id, $desc) = split(/\s/, $head, 2); # split id and desc
$seq =~ tr/\n//d; # remove newlines
my %fa = (
id => $id,
desc => $desc,
seq => $seq,
);
return %fa;
}
What do you mean by "
sum result of all the sequences"