Re: pattern finding algorithm

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;
}
[download]

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
[download]

Comment on Re: pattern finding algorithm Select or Download Code

Replies are listed 'Best First'.
Re^2: pattern finding algorithm by lima1 (Curate) on Aug 01, 2007 at 07:42 UTC
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 ) { ... } } }` [download]	[reply] [d/l]
Re^3: pattern finding algorithm by dogz007 (Scribe) on Aug 01, 2007 at 13:59 UTC
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 1999sum(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.	[reply]
Re^4: pattern finding algorithm by lima1 (Curate) on Aug 01, 2007 at 14:53 UTC
Maybe I just need some coffee, but GrandFathers algorithm generates the list of all maximum patterns including the sequence ids. In your example, the hash would contain something like `$matches{$pattern} = [ [1,3], ... [ 5, 6 ], ];` [download] after sum(1..n-1) comparsions. Isn't that what kdt2006 asked for?	[reply] [d/l]
Re^2: pattern finding algorithm by Anonymous Monk on Aug 01, 2007 at 09:06 UTC
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	[reply]
Re^3: pattern finding algorithm by lima1 (Curate) on Aug 01, 2007 at 09:30 UTC
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.	[reply]
Re^3: pattern finding algorithm by kdt2006 (Initiate) on Aug 01, 2007 at 09:10 UTC
sorry, just to let you know the previous post was mine. Forgot to sign in :O	[reply]
Re^2: pattern finding algorithm by roboticus (Chancellor) on Aug 01, 2007 at 10:57 UTC
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	[reply]