Beefy Boxes and Bandwidth Generously Provided by pair Networks
The stupid question is the question not asked
 
PerlMonks  

pattern finding algorithm

by kdt2006 (Initiate)
on Jul 31, 2007 at 23:13 UTC ( #629941=perlquestion: print w/replies, xml ) Need Help??

kdt2006 has asked for the wisdom of the Perl Monks concerning the following question:

Hi, I'm checking to see if you guys may be able to help me with an algorithm for finding patterns. I have around 2000 short sequences (of length 9) that are aligned. I want to be able to extract all common patterns on the same positions and report the number of occurrences. For example in the following: 1.ACGCATTCA, 2.ACTGGATAC, 3.TCAGCCATC I would like the following output (where a full stop represents any character): (AC....T..) 2 occurrences -pattern between sequence 1 and 2. (.C.G....C) 2 occurrences -pattern between sequence 2 and 3. (.C.......) 2 occurrences -pattern between sequence 1 and 3. As you can see, the way that I am planning on doing this now requires sum(n-1...1) comparisons. Is there a more efficient way of doing this with less comparisons? If there are any existing algorithms that you think may be better suited, please let me know. Thanks

Replies are listed 'Best First'.
Re: pattern finding algorithm
by GrandFather (Saint) on Jul 31, 2007 at 23:46 UTC

    If you really want to look for any individual common characters between all possible pairs of 2000 sequences then you just gotta do the work, and there is a lot of work to do!

    It may help to tell us why you want to do that. There may be a better solution to your problem than the brute force search implied so far.

    You may find this code interesting however:

    use strict; use warnings; my @sequences = qw(ACGCATTCA ACTGGATAC TCAGCCATC); my %matches; for my $outer (0 .. $#sequences - 1) { for my $inner ($outer + 1 .. $#sequences) { my $mask = $sequences[$outer] ^ $sequences[$inner]; next if index ($mask, "\0") == -1; # No matching characters $mask =~ tr/\0/\xff/c; $mask |= $sequences[$outer]; $mask =~ tr/\xff/./; push @{$matches{$mask}}, [$outer + 1, $inner + 1]; } } for my $match (sort keys %matches) { print "$match pattern between ", join (', ', map {"$_->[0] and $_->[1]"} @{$matches{$match}}), "\n"; }

    Prints:

    .C....... pattern between 1 and 3 .C.G....C pattern between 2 and 3 AC....T.. pattern between 1 and 2

    DWIM is Perl's answer to Gödel
Re: pattern finding algorithm
by Fletch (Chancellor) on Jul 31, 2007 at 23:48 UTC

    Given the sample data in question, I'm going to hazard a guess that you might want to look at bioperl and see if there's not something off the shelf that does what you're looking for.

Re: pattern finding algorithm
by dogz007 (Scribe) on Aug 01, 2007 at 07:18 UTC
    Let's explore the overall efficiency. There are nine positions and four different letters (A,C,G,T). This makes 4^9 = 262,144 unique permutations. Now examine the form of a pattern. Each position in the pattern can either be a match (which you've represented above as the actual matching letter) or a no match (which you've represented with a full stop. So every pattern can have two options (match or no match) for each of the nine positions, making only 2^9-1 = 511 unique permutations for patterns (subtract one for the case of no matches at all). Taking advantage of this fact may significantly reduce the number of comparisons necessary to find all possible matches.

    To generalize this method, let's assume that there are n sequences of m positions. Checking each pattern against all of the sequences would require (2^m-1)*sum(1..n-1) comparisons. Below is my solution to the problem using this method.

    use strict; use List::MoreUtils qw(uniq); # initialize sequence data and emtpy match set my $seq = [ [ qw(A C G C A T T C A) ], [ qw(A C T G G A T A C) ], [ qw(T C A G C C A T C) ] ]; my $n = scalar(@$seq); # num of sequences my $m = scalar(@{$seq->[0]}); # num of positions my $match = {}; # first iterate over all the possible patterns for the sequence data for my $pat (1 .. 2**$m-1) { # now create the pattern by converting the pattern number to # binary and then to a slice with which to compare sequences. # pretty slick, if you ask me. my @pattern = reverse split //, binary($pat); my @slice; for (0 .. $m-1) { push @slice, $_ if $pattern[$_] } # now iterate over every sequence. check to see if it matches any # of the others based on the given pattern. for my $i (0 .. $n-2) { for my $j ( $i+1 .. $n-1 ) { if ( join("", @{$seq->[$i]}[@slice]) eq join("", @{$seq->[$j]}[@slice]) ) { # construct a hash key and store the matches my $key = join "", map { $pattern[$_] ? $seq->[$i][$_] : "." } (0 .. $m-1); push @{$match->{$key}}, $i,$j ; } } } } # now strip out the duplicate sequence numbers for each match (an # unfortunate step I had trouble avoiding), then display the matches for my $key (keys %$match) { @{$match->{$key}} = uniq @{$match->{$key}}; print $key , " occurs in " , join(",", @{$match->{$key}}) , "\n"; } # a short decimal to binary converter I borrowed from Mark Dominus sub binary { my ($n) = @_; return $n if $n == 0 || $n == 1; my $k = int($n/2); my $b = $n % 2; my $E = binary($k); return $E . $b; }

    Prints the following matches (note there are more than previously thought):
    ...G....C occurs in 1,2 .C......C occurs in 1,2 A.....T.. occurs in 0,1 .C.G....C occurs in 1,2 ........C occurs in 1,2 ...G..... occurs in 1,2 ......T.. occurs in 0,1 .C.G..... occurs in 1,2 A........ occurs in 0,1 AC....... occurs in 0,1 .C....T.. occurs in 0,1 AC....T.. occurs in 0,1 .C....... occurs in 0,1,2
      I haven't tried to understand every line, but where exactly is the performance gain compared to the naive algorithm? Instead of comparing m positions you compare 2^(m-1) patterns?
      for my $pat (1 .. 2**$m-1) { ... for my $i (0 .. $n-2) { for my $j ( $i+1 .. $n-1 ) { ... } } }
        Let me see if I can explain, limal. There are n sequences, each of length m. A pattern may consist of a match (1) or no-match (0) in each position of the sequence, such as 100101001. Think of it as a match "template". There are 2^m-1 patterns (disregarding the case of all zeros). To find all the matches, use each match template to compare every sequence to every other sequence, for a total of (2^m-1)*sum(1..n-1) comparisons.

        Now let's see how many it would take to find all of them using the "naive" algorithm. The algorithm used by GrandFather makes only sum(1..n-1) comparisons. But consider a hypothetical set of 6 sequences. At the end of the naive algorithm, sequence 5 might match 6 in some way. However, this match may also occur in earlier sequences. So to find these would require that this specific match be checked against all the other sequences. Therefore, to find all possible matches would require (n-1)*sum(1..n-1) comparisons.

        Depending on the sizes of m and n, either method may be faster. For the case at hand (m=9,n=2000), the naive algorithm requires 1999*sum(1..1999)=3,996,001,000 comparisons, while my pattern template algorithm only requires (2^9-1)*sum(1..1999)=1,021,489,000 comparisons. This makes my algorithm faster (meaning fewer comparisons) by a factor of (n-1)/(2^m-1)=3.9119 . As an added bonus, it also finds all of the submatches, although kdt2006 has expressed that they are not needed.
      Wow! This is my first post and I have to say that I really am impressed with the breadth of knowledge you guys have. The code that I have already is very similar to Godfathers (although written in python). I have used a bit of perl in the past, and I noticed you fellas were pretty good at algorithms so I decided to post here.

      Thanks to the contribution by all of you, I was surprised to see Limal's link to the WebLogo as I plan on doing something similar to it for this project.

      Many thanks to your idea, converting to binary is not something I would have thought of - I guess you can always learn something new.

      The actual sequences being used will be protein sequences (so 20 possibilities each character instead of 4).

      The purpose of this is to try to find position specific patterns in the sequence to try to estimate the binding properties of the sequences (I have additional binding score data).

      In your example, I would like to find a way to display only maximal patterns. For example, in your results between sequence 1 and 2 there are 7 sub patterns (as these can all be said to be sub patterns of (.C.G....C) on line 4). When sub patterns occur on same instances, I would like to only return the maximal pattern. Thanks

        Haven't done something similar yet, but my first idea: Calculate a distance matrix with one of the standard distance measures for protein data for the sequences with known binding scores. Then use a simple cluster algorithm like UPGMA to find clusters with good scores and visualize these sequences in weblogo. If you want to predict binding scores, this could at least help to find some descriptors for machine learning tools.

        sorry, just to let you know the previous post was mine. Forgot to sign in :O

      dogz007:

      Hmmm ... Your output doesn't limit the results to the "best" matches. Specifically, line 4 shows the "best" match between 1 & 2, but you're also reporting subfragments that match (lines 1, 2, 5, 6 & 8). Similarly, the "best" match between 0 & 1 is in the penultimate line.

      ...roboticus

Re: pattern finding algorithm
by salva (Canon) on Aug 01, 2007 at 06:36 UTC
    If you just need to find "good" patterns in your data (and for "good" I mean, patterns very specific and that match a high number of entries), I think that, funnily, a genetic algorithm could produce good results.
Re: pattern finding algorithm
by lima1 (Curate) on Aug 01, 2007 at 07:13 UTC
Re: pattern finding algorithm
by dogz007 (Scribe) on Aug 01, 2007 at 12:41 UTC
    kdt2006, please correct if I'm wrong, but I read your last reply to mean that you want the maximal match in each subset of sequences. For example, I added a new sequence to my code and got the following results:

    my $seq = [ [ qw(A C G C A T T C A) ], [ qw(A C T G G A T A C) ], [ qw(T C A G C C A T C) ], [ qw(C T G G A T C G C) ] ];

    Outputs:
    ..G.AT... occurs in 0,3 ...G....C occurs in 1,2,3 .C......C occurs in 1,2 A.....T.. occurs in 0,1 ........C occurs in 1,2,3 ...G..... occurs in 1,2,3 A........ occurs in 0,1 ....AT... occurs in 0,3 AC....... occurs in 0,1 ..G.A.... occurs in 0,3 ..G..T... occurs in 0,3 .C.G....C occurs in 1,2 ..G...... occurs in 0,3 ......T.. occurs in 0,1 .C.G..... occurs in 1,2 .C....T.. occurs in 0,1 AC....T.. occurs in 0,1 .....T... occurs in 0,3 ....A.... occurs in 0,3 .C....... occurs in 0,1,2

    Notice above how the maximal match between 1 and 2 is .C.G....C, while a submatch of this, ...G....C, is the maximal match between 1, 2, and 3. From what I understand you would want both of these, even though one is a submatch of the other. You would not, however, want ........C or ...G..... because these are submatches of the maximal match for 1,2,3. I still believe this is a better way to find all possible matches. Perhaps the way to best suit it to your needs would be to implement a filter that finds the maximal match for each sequence subset. The following update to my original code does the trick. (Note that true and uniq require List::MoreUtils.

    my %maxMatch; for my $key (keys %$match) { @{$match->{$key}} = uniq @{$match->{$key}}; my $seqSet = join ",", @{$match->{$key}}; unless ($maxMatch{$seqSet}) { $maxMatch{$seqSet} = $key; next; } my $old = true { /[^\.]/ } split //, $maxMatch{$seqSet}; my $new = true { /[^\.]/ } split //, $key; $maxMatch{$seqSet} = $key if $new > $old; } for my $key (keys %maxMatch) { print $maxMatch{$key}, " occurs in $key \n"; }

    Prints (for the four sequence set above):
    ..G.AT... occurs in 0,3 .C.G....C occurs in 1,2 ...G....C occurs in 1,2,3 .C....... occurs in 0,1,2 AC....T.. occurs in 0,1

    Update: I was flipping through my Perl in a Nutshell text and discovered that grep could (and probably should) be used in place of true in the code above. I tested both and they work the same, although grep may run a little faster since it's a built-in command.

      hi dogz07,

      Many thanks, this is what I was looking for.

      Thanks everyone for all the help. Limal, I will take a look at your suggestion.

      Many thanks again

      hello dogz, i tried using the messaging but i couldnt see a text box to type in. ok im a little new here but i was wondering if you could help me in find a pattern in some longer codes that deal with both numbers and letters if you can please email me at tottinum1@gmail.com

Log In?
Username:
Password:

What's my password?
Create A New User
Node Status?
node history
Node Type: perlquestion [id://629941]
Approved by GrandFather
help
Chatterbox?
and the web crawler heard nothing...

How do I use this? | Other CB clients
Other Users?
Others imbibing at the Monastery: (4)
As of 2021-06-24 23:09 GMT
Sections?
Information?
Find Nodes?
Leftovers?
    Voting Booth?
    What does the "s" stand for in "perls"? (Whence perls)












    Results (133 votes). Check out past polls.

    Notices?