Neel
Bioinformatics Beginner
Our Sponsors
Perl script to create a consensus of nucleotide sequences !
By Neel 2233 days ago
use strict;
use warnings;
my @instances = qw ( AAAAA ATCGA ATAAA );
my @instances2 = qw ( AAAAA AACGA ATAAA AGAAA AGAAA);
print consensus(@instances),"\n"; # ATAAA
print consensus(@instances2),"\n"; # ATAAA
exit;
sub consensus{
my @mi = @_;
chomp(@mi);
my $motif_count=0;
my @words =();
my %H = ( A=>[], T=>[], C=>[], G=>[] );
s/\s//g for @mi;
my ($w) = sort {$b <=> $a} map {length} @mi; # set w to the length of the longest element
foreach my $j ( 0 .. $w-1 ){
# Initialize the base counts.
my %h = ( a=>0, t=>0, c=>0, g=>0 );
my @mi_letters = map { [split '', uc $_] } @mi;
foreach my $j ( 0 .. $w-1 ){
$H{ $_->[$j] }->[$j]++ for @mi_letters;
}
push @{$H{ uc $_ }}, $h{$_} for keys %h; # example: push @{$H{G}}, $g;
}
my @cons = ();
my %prefOrder = ( A=>1, T=>2, C=>3, G=>4 );
foreach my $B ( 0 .. $w-1 ){
push @cons, [ sort { ($H{$b}->[$B]||0) <=> ($H{$a}->[$B]||0) || $prefOrder{$b} <=> $prefOrder{$a} } qw/A T G C/ ]->[0];
}
return @cons;
}
#reference https://www.perlmonks.org/bare/?node_id=500962
