This doesn't change your functionality at all (i don't think), but does shorten all the code up and hopefully makes it much easier to modify from here.. Basically it makes use of hashes instead of having to declare a bunch of variables (e.g. $a, $g, $t, $c) so that everything is easier to work with.

One thing i didn't understand was the $w usage -- i assumed in my code that you want $w to be the length of the longest item, but your code was setting it to be the length of the last item...
use strict; use warnings; my @instances = qw ( AAAAA ATCGA ATAAA ); my @instances2 = qw ( AAAAA AACGA ATAAA ATAAA ); 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 lengt +h of the longest element foreach my $j ( 0 .. $w-1 ){ # Initialize the base counts. my %h = ( a=>0, t=>0, c=>0, g=>0 ); foreach my $mi ( @mi ) { $h{ lc substr($mi,$j,1) }++; # example: $h{g}++ +; } 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] <=> $H{$a}->[$B] || $prefOrder +{$b} <=> $prefOrder{$a} } qw/A T G C/ ]->[0]; } return @cons; }
Update: more efficient way to set %H would be:
my %H = ( A=>[], T=>[], C=>[], G=>[] ); my @mi_letters = map { [split '', uc $_] } @mi; foreach my $j ( 0 .. $w-1 ){ $H{ $_->[$j] }->[$j]++ for @mi_letters; } # note: to avoid warnings, change the sort further down to include t +his (the '||0'): # ($H{$b}->[$B]||0) <=> ($H{$a}->[$B]||0)

In reply to Re: Finding Consensus String by davidrw
in thread Finding Consensus String by monkfan

Title:
Use:  <p> text here (a paragraph) </p>
and:  <code> code here </code>
to format your post, it's "PerlMonks-approved HTML":



  • Posts are HTML formatted. Put <p> </p> tags around your paragraphs. Put <code> </code> tags around your code and data!
  • Titles consisting of a single word are discouraged, and in most cases are disallowed outright.
  • Read Where should I post X? if you're not absolutely sure you're posting in the right place.
  • Please read these before you post! —
  • Posts may use any of the Perl Monks Approved HTML tags:
    a, abbr, b, big, blockquote, br, caption, center, col, colgroup, dd, del, details, div, dl, dt, em, font, h1, h2, h3, h4, h5, h6, hr, i, ins, li, ol, p, pre, readmore, small, span, spoiler, strike, strong, sub, summary, sup, table, tbody, td, tfoot, th, thead, tr, tt, u, ul, wbr
  • You may need to use entities for some characters, as follows. (Exception: Within code tags, you can put the characters literally.)
            For:     Use:
    & &amp;
    < &lt;
    > &gt;
    [ &#91;
    ] &#93;
  • Link using PerlMonks shortcuts! What shortcuts can I use for linking?
  • See Writeup Formatting Tips and other pages linked from there for more info.