Perl script to create a consensus of nucleotide sequences !

Neel — Fri, 12 Oct 2018 10:01:22 -0500

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

BOL: Perl script to create a consensus of nucleotide sequences !

Perl script to create a consensus of nucleotide sequences !