Beefy Boxes and Bandwidth Generously Provided by pair Networks
Just another Perl shrine

Re: Using grep in a scalar context

by perlhappy (Novice)
on Feb 06, 2013 at 16:12 UTC ( #1017464=note: print w/replies, xml ) Need Help??

in reply to Using grep in a scalar context

Hi newbie1991,

I have a good idea of what you want to do and the data your dealing with (I do quite a lot of bioinformatics based work).

Anyway, I've written two scripts for you to look at. I've kept the code simple and commented so you should be ok with it. Over a chromosome, this mightnt be as fast as it could be but should be ok.

Firstly, you wont want to split the sequence into an array, unless you are absolutely sure you arent going to miss out a count on an odd number occurrance of the acid eg FAAAD would be split into FA-AA-D? or ?F-AA-AD so you would only count AA once whereas its actually got 2 pairs.

This is the first script. This will find only AA pairs and count them. The sequence is ASDTDAAFRASEQSAAAFDG (its in the code) so the number of AA's should be 3.

#!usr/bin/perl -w use strict; my $string = "ASDTDAAFRASEQSAAAFDG"; my $counter = 0; for(my $i = 0; $i < length($string)-1; $i++){ my $amino = substr($string, $i, 2); $counter++ if($amino eq "AA"); #print "$amino\n"; } print "number of AA matches = $counter\n";

This is the second script, its a bit more complicated. Instead of only counting the number of AA's it will count all pairs. It creates these on the fly and if it encounters one it has already created it just increments the value. Just to note for the 22 possible amino acids the number of possible pairs will be much higher (484 I think)

#!usr/bin/perl -w use strict; my $string = "ASDTDAAFRASEQSAAAFDG"; my %acids; for(my $i = 0; $i < length($string)-1; $i++){ my $amino = substr($string, $i, 2); if(exists $acids{$amino}){ $acids{$amino}++; }else{ $acids{$amino} = 1; } #print "$amino\n"; } print "These are the occurrence of all amino acid pairs including over +laps as separate counts\n"; foreach my $amino(keys %acids){ print "$amino\t$acids{$amino}\n"; }

I hope either of these do what you want

Replies are listed 'Best First'.
Re^2: Using grep in a scalar context
by AnomalousMonk (Canon) on Feb 07, 2013 at 16:38 UTC

    I have no bioinformatic background, but I'd like to offer a couple of comments on your code, specifically the version that counts overlapping letter pairs (would 'digrams' be an appropriate term for these?).

    my %acids; for(my $i = 0; $i < length($string)-1; $i++){ my $amino = substr($string, $i, 2); if(exists $acids{$amino}){ $acids{$amino}++; }else{ $acids{$amino} = 1; } #print "$amino\n"; }

    Because it is not necessary to check for the existence of a hash key before incrementing its value (due to autovivification), the body of this for-loop can be reduced to a single statement:
        ++$acids{ substr $string, $i, 2 }
    This will almost certainly yield a speed benefit.

    Alternatively, in 5.10+ versions of Perl, the entire for-loop can be replaced by a single regex (tested):
        $string =~ m{ (?= (..) (?{ ++$pairs2{$^N} }) (*FAIL)) }xms;
    This may or may not increase speed; you will have to Benchmark this for yourself. The alternate regex
        m{ (?= .. (?{ ++$pairs2{${^MATCH}} }) (*FAIL)) }xmsp
    also works (note the additional  /p regex modifier) and may be slightly faster because no capturing group is used. Again, Benchmark-ing will tell the tale.

    >perl -wMstrict -le "use Test::More tests => 2; use Data::Dump; ;; my $string = 'ABCCCDEAB'; ;; my %pairs1; $pairs1{$_}++ for $string =~ /(?=(..))/g; ;; local our %pairs2; $string =~ m{ (?= .. (?{ ++$pairs2{${^MATCH}} }) (*FAIL)) }xmsp; ;; my %pairs3; for (my $i = 0; $i < length($string) - 1; ++$i) { ++$pairs3{ substr $string, $i, 2 } } ;; dd \%pairs1, \%pairs2, \%pairs3; is_deeply \%pairs1, \%pairs2, '1 & 2, same results'; is_deeply \%pairs1, \%pairs3, '1 & 3, same results'; " 1..2 ( { AB => 2, BC => 1, CC => 2, CD => 1, DE => 1, EA => 1 }, { AB => 2, BC => 1, CC => 2, CD => 1, DE => 1, EA => 1 }, { AB => 2, BC => 1, CC => 2, CD => 1, DE => 1, EA => 1 }, ) ok 1 - 1 & 2, same results ok 2 - 1 & 3, same results
Re^2: Using grep in a scalar context
by newbie1991 (Acolyte) on Feb 12, 2013 at 13:58 UTC
    So I went back to the start after posting this and essentially came up with your first option, I'm glad you recommended it as well because I was worried that it was childish/too roundabout. Thanks so much everyone. :)

Log In?

What's my password?
Create A New User
Node Status?
node history
Node Type: note [id://1017464]
and all is quiet...

How do I use this? | Other CB clients
Other Users?
Others exploiting the Monastery: (3)
As of 2017-01-19 01:50 GMT
Find Nodes?
    Voting Booth?
    Do you watch meteor showers?

    Results (167 votes). Check out past polls.