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

Dear all,

I have a sequence:

GAATGTTTTAGCAATCTCTTTCTGTCATGAATCCATGGCAGTGACCATACTAATGGTGACTGCCATTGAT
+GGAGGGAGACACA

[download]

and I want to be able to find this word in the sequence:

CTGGATAAGAATGTTTTAGCAATCTCTT

[download]
Problem is, the word is found at the beginning of the sequence, and in a truncated form (first 8 letters of word missing).

I was considering looping through successively smaller substrings of the word until it is found. I'd have to do it from either end of the word too, when trying to match an array of sequences.

My question is, is there a simpler way of using regexp to say find partial match?

I tried googling this but couldn't find the right combination of words (truncated/partial/overlap etc.) so any links to a solution would be greatly appreciated.

Many thanks
Sam

20061019 Janitored by Corion: Added code tags around sequences to enable wrapping

Replies are listed 'Best First'.
Re: Searching for a word that may only exist in part
by GrandFather (Saint) on Oct 18, 2006 at 23:09 UTC

    This may point you in a helpful direction:

    use strict;
use warnings;

my $seq = 'GAATGTTTTAGCAATCTCTTTCTGTC' .
    'ATGAATCCATGGCAGTGACCATACTAAT' .
    'GGTGACTGCCATTGATGGAGGGAGACACA';
    
my $word = 'CTGGATAAGAATGTTTTAGCAATCTCTT';

print slidingMatch($seq, $word);

sub slidingMatch {
    my ($seq, $word) = @_;

    return "Complete match at $-[1]\n" if $seq =~ /($word)/;
    
    my @bestMatch;
    
    # Try matching start
    my $str = substr $seq, 0, 3;
    my $wordLen = length $word;
    
    while ($word =~ /$str/g) {
        my $tailLen = length ($word) - (my $start = $-[0]);
        my $subWord = substr ($word, $wordLen - $tailLen);
        
        next if $seq !~ m/^$subWord/;
        @bestMatch = ($start, $tailLen);
        last;
    }
    
    # Try matching end
    my $rword = reverse $word;
    my $rseq = reverse $seq;
    $str = substr $rseq, 0, 3;
    
    while ($rword =~ /$str/g) {
        my $tailLen = length ($rword) - (my $start = $-[0]);
        
        last if @bestMatch && $tailLen < $bestMatch[1];

        my $subWord = substr ($rword, $wordLen - $tailLen);
        
        next if $rseq !~ m/^$subWord/;
        @bestMatch = (length ($seq) - $tailLen, $tailLen);
        last;
    }

    return "Partial match found at $bestMatch[0] for $bestMatch[1] cha
+racters\n" if @bestMatch;
    return "No match found for $word in $seq\n";
}

    [download]

    Prints:

    Partial match found at 8 for 20 characters

    [download]

    Update: fixed bug where worse end match replaced start match.

    Update: fixed tail match position calculation bug (see BrowserUk's reply).

    DWIM is Perl's answer to Gödel
[reply]
[d/l]
[select]

      You still seem to have a bug occasionally. You report an offset of 75, but I make it to be 59? 

      #! perl -slw
use strict;

sub slidingMatch {
    my ($seq, $word) = @_;

    return "Complete match at $-[1]\n" if $seq =~ /($word)/;

    my @bestMatch;

    # Try matching start
    my $str = substr $seq, 0, 3;
    my $wordLen = length $word;

    while ($word =~ /$str/g) {
        my $tailLen = length ($word) - (my $start = $-[0]);
        my $subWord = substr ($word, $wordLen - $tailLen);

        next if $seq !~ m/^$subWord/;
        @bestMatch = ($start, $tailLen);
        last;
    }

    # Try matching end
    my $rword = reverse $word;
    my $rseq = reverse $seq;
    $str = substr $rseq, 0, 3;

    while ($rword =~ /$str/g) {
        my $tailLen = length ($rword) - (my $start = $-[0]);

        last if @bestMatch && $tailLen < $bestMatch[1];

        my $subWord = substr ($rword, $wordLen - $tailLen);

        next if $rseq !~ m/^$subWord/;
        @bestMatch = (length ($seq) - $start, $tailLen);
        last;
    }

    return "Partial match found at $bestMatch[0] for $bestMatch[1] cha
+racters\n" if @bestMatch;
    return "No match found for $word in $seq\n";
}

print slidingMatch 
#012345678 012345678 012345678 012345678 012345678 012345678 
'GAATGTTTTAGCAATCTCTTTCTGTCATGAATCCATGGCAGTGACCATACTAATGGTGACTGCCATTGA
+TGGAGGGAGACACA',
                                                           'CTGCCATTGA
+TGGAGGGAGACACAACGTACGT';
__END__
c:\test>junk
Partial match found at 75 for 24 characters

      [download]
      Examine what is said, not who speaks -- Silence betokens consent -- Love the truth but pardon error.
      Lingua non convalesco, consenesco et abolesco. -- Rule 1 has a caveat! -- Who broke the cabal?
      "Science is about questioning the status quo. Questioning authority".
      In the absence of evidence, opinion is indistinguishable from prejudice.
[reply]
[d/l]
Re: Searching for a word that may only exist in part
by BrowserUk (Patriarch) on Oct 19, 2006 at 05:29 UTC

    This seems to work. Run it on a terminal 100+ wide to see the output properly; the wrapping screws it up. 

    #! perl -slw
use strict;

sub search {
    my( $haystack, $needle, $pos ) = @_;
    my $l = length $needle;
    my( @res1, @res2 );

    return( $pos-1, 0, $l ) if $pos = 1+index $haystack, $needle;

    substr( $haystack, 0, $_ ) eq substr( $needle, -$_ )
        and @res1 = ( 0, $l-$_, $_ ), last
        for reverse 1 .. $l-1;

    substr( $haystack, $_-$l ) eq substr( $needle, 0, $l-$_ )
        and @res2 = ( length( $haystack )-( $l-$_), -$_, $l-$_  ), las
+t
        for 1 .. $l-1;

    return       unless @res1 or @res2;
    return @res1 unless @res2;
    return @res2 unless @res1;
    return $res1[2] > $res2[2] ? @res1 : @res2;
}

my $haystack = 'GAATGTTTTAGCAATCTCTTTCTGTCATGAATCCATGGCAGTGACCATACTAAT
+GGTGACTGCCATTGATGGAGGGAGACACA';
my @needles = (
    'NOT FOUND',
    'CTGGATAAGAATGTTTTAGCAATCTCTT',
    'CTGCCATTGATGGAGGGAGACACAACGTACGT',
    'CATGAATCCATGGCAGTGACCATACTAATGGTGACTG',
    'AGGGAGACACAGAATGTTTTAG',
    'AGGGAGACACAGAATGTTTTAGC',
);

$, = ' ';
for my $needle ( @needles ) {
    my( $hOffset, $nOffset, $l ) = search $haystack, $needle;
    print "No match found in '$haystack'\n for '$needle'\n"
        and next unless defined $hOffset;

    printf "%s%s\n%s%s\n\n",
        ' 'x$nOffset, $haystack,
        ' 'x$hOffset, $needle;
}
__END__
c:\test>579215
No match found in 'GAATGTTTTAGCAATCTCTTTCTGTCATGAATCCATGGCAGTGACCATACT
+AATGGTGACTGCCATTGATGGAGGGAGACACA'
 for 'NOT FOUND'

        GAATGTTTTAGCAATCTCTTTCTGTCATGAATCCATGGCAGTGACCATACTAATGGTGACTG
+CCATTGATGGAGGGAGACACA
CTGGATAAGAATGTTTTAGCAATCTCTT

GAATGTTTTAGCAATCTCTTTCTGTCATGAATCCATGGCAGTGACCATACTAATGGTGACTGCCATTGAT
+GGAGGGAGACACA
                                                           CTGCCATTGAT
+GGAGGGAGACACAACGTACGT

GAATGTTTTAGCAATCTCTTTCTGTCATGAATCCATGGCAGTGACCATACTAATGGTGACTGCCATTGAT
+GGAGGGAGACACA
                         CATGAATCCATGGCAGTGACCATACTAATGGTGACTG

GAATGTTTTAGCAATCTCTTTCTGTCATGAATCCATGGCAGTGACCATACTAATGGTGACTGCCATTGAT
+GGAGGGAGACACA
                                                                      
+  AGGGAGACACAGAATGTTTTAG

           GAATGTTTTAGCAATCTCTTTCTGTCATGAATCCATGGCAGTGACCATACTAATGGTGA
+CTGCCATTGATGGAGGGAGACACA
AGGGAGACACAGAATGTTTTAGC

    [download]
    Examine what is said, not who speaks -- Silence betokens consent -- Love the truth but pardon error.
    Lingua non convalesco, consenesco et abolesco. -- Rule 1 has a caveat! -- Who broke the cabal?
    "Science is about questioning the status quo. Questioning authority".
    In the absence of evidence, opinion is indistinguishable from prejudice.
[reply]
[d/l]
      Dear BrowserUk,

      I ran your search subroutine on some hard data, three 'needles', and three different 'haystacks'. I realised that you didn't set a lower limit of consecutive matches, as you'll see from the output below, it even matches when just the first letter of the needle matches the last letter of the haystack.

      I was going to try and edit your code myself, but the collection of return statements had me confused, so I was wondering, what would you do to your code to set that lower limit, which should probably be a higher number, like 6, now that I think about it.

      Cheers
      Sam 

      TTGTCAGCGAAAAAAATTAAAGCGCAAGATTGTTGGTTTTTGCGTGATGGTGACCGGGCAGCCTAAAGGC
+TATCCTTAACCAGGGAGCTGATT
      GCGAAAAAAATTAAAGCGCAAGATTGTTGGTTTTTGCGTGATGGTGAC
 
                                              TTGTCAGCGAAAAAAATTAAAGCG
+CAAGATTGTTGGTTTTTGCGTGATGGTGACCGGGCAGCCTAAAGGCTATCCTTAACCAGGGAGCTGATT
CGCGACAACCGGAATATGAAAGCAAAGCGCAGCGTCTGAATAACGTTT

                  GAATGTTTTAGCAATCTCTTTCTGTCATGAATCCATGGCAGTGACCATACTA
+ATGGTGACTGCCATTGATGGAGGGAGACACAGTGCACTGGCAAACTCACAC
CATTACATTGCTGGATAAGAATGTTTTAGCAATCTCTTTCTGTCATGA
 
GAATGTTTTAGCAATCTCTTTCTGTCATGAATCCATGGCAGTGACCATACTAATGGTGACTGCCATTGAT
+GGAGGGAGACACAGTGCACTGGCAAACTCACAC
                                                                      
+                                CGCGACAACCGGAATATGAAAGCAAAGCGCAGCGTCT
+GAATAACGTTT
 
TAATCAAAACCAATAAACACGAAATAATCCCCATGCCGGTGAAGAAGGGGCGTGACTTTAGCGAAATGTT
+GCCGTCGCGACAACCGGAATATGAAAGCAAAGCGCAGCGTCTGAATAACGTTTATGCTGAAAGCGGATG
+AATAAGGAGATGCG
                                                                      
+                                                                     
+           GCGAAAAAAATTAAAGCGCAAGATTGTTGGTTTTTGCGTGATGGTGAC
 
TAATCAAAACCAATAAACACGAAATAATCCCCATGCCGGTGAAGAAGGGGCGTGACTTTAGCGAAATGTT
+GCCGTCGCGACAACCGGAATATGAAAGCAAAGCGCAGCGTCTGAATAACGTTTATGCTGAAAGCGGATG
+AATAAGGAGATGCG
                                                                      
+     CGCGACAACCGGAATATGAAAGCAAAGCGCAGCGTCTGAATAACGTTT

      [download]
      Update: Removed redundant lines
[reply]
[d/l]

        All that is needed is a test in each loop to quit when the partial match attempt is less than whatever minimum you want to specify. Rather than duplicate the test in both loops, I've amalgamated them into a single loop which allows for a single test, with the nice side-effect of being 25% more efficient than the original above.

        I've also added some explanitary comments and (I think) improved the variable names. Let me know if anything is unclear.

        I've posted the output of a run with the minimum match set to 3, to test the breakpoint, but set it to whatever value makes sense for your purposes. 

        #! perl -slw
use strict;
use constant MINMATCH => 3;

## returns a list of 3 integers.
## offset into the haystack
## offset into the needle ( can be negative, subtract from length(need
+le) )
## length of match.
sub search {
    my( $haystack, $needle, $pos ) = @_;
    my( @start, @end );
    my $l = length $needle;

    ## A full match found
    return( $pos-1, 0, $l ) if $pos = 1+index $haystack, $needle;

    ## iterate, and attempt matches at both ends
    for( 1 .. $l-1 ) {
        my $r = $l - $_;   ## reverse offset

        last if $r < MINMATCH; ## quit if we've reached the minimum ma
+tch length

        ## try a partial match at the beginning
        substr( $haystack, $_-$l ) eq substr( $needle, 0, $l-$_ )
            and @start = ( length( $haystack )-( $l-$_), -$_, $l-$_  )
+;

        ## try a partial match at the end
        substr( $haystack, 0, $r ) eq substr( $needle, -$r )
            and @end = ( 0, $_, $r );

        ## Quit if we got either.
        last if @start or @end;
    }
    return        unless @start or @end;  ## No match
    return @start unless @end;            ## No partial at the end, st
+art is best
    return @end   unless @start;          ## No partial at the start, 
+end is best
    return $start[2] >= $end[2] ? @start : @end; ## Got both, longest 
+(or start if equal) is best
}

our @haystacks = (
    'TTGTCAGCGAAAAAAATTAAAGCGCAAGATTGTTGGTTTTTGCGTGATGGTGACCGGGCAGCCTA
+AAGGCTATCCTTAACCAGGGAGCTGATT',
    'GAATGTTTTAGCAATCTCTTTCTGTCATGAATCCATGGCAGTGACCATACTAATGGTGACTGCCA
+TTGATGGAGGGAGACACAGTGCACTGGCAAACTCACAC',
    'TAATCAAAACCAATAAACACGAAATAATCCCCATGCCGGTGAAGAAGGGGCGTGACTTTAGCGAA
+ATGTTGCCGTCGCGACAACCGGAATATGAAAGCAAAGCGCAGCGTCTGAATAACGTTTATGCTGAAAGC
+GGATGAATAAGGAGATGCG',
);
our @needles = (
    'GCGAAAAAAATTAAAGCGCAAGATTGTTGGTTTTTGCGTGATGGTGAC',
    'CATTACATTGCTGGATAAGAATGTTTTAGCAATCTCTTTCTGTCATGA',
    'CGCGACAACCGGAATATGAAAGCAAAGCGCAGCGTCTGAATAACGTTT',
);

$, = ' ';
for my $needle ( @needles ) {
    for my $haystack ( @haystacks ) {
        my( $hOffset, $nOffset, $l ) = search $haystack, $needle;
        if( defined $hOffset ) {
            print "$hOffset, $nOffset, $l";
            printf "%s%s\n%s%s\n\n",
                ' 'x$nOffset, $haystack,
                ' 'x$hOffset, $needle;
        }
        else {
            print "No match found in '$haystack'\n for '$needle'\n";
        }
    }
}

__END__
C:\test>579215-2
6, 0, 48
TTGTCAGCGAAAAAAATTAAAGCGCAAGATTGTTGGTTTTTGCGTGATGGTGACCGGGCAGCCTAAAGGC
+TATCCTTAACCAGGGAGCTGATT
      GCGAAAAAAATTAAAGCGCAAGATTGTTGGTTTTTGCGTGATGGTGAC

No match found in 'GAATGTTTTAGCAATCTCTTTCTGTCATGAATCCATGGCAGTGACCATACT
+AATGGTGACTGCCATTGATGGAGGGAGACACAGTGCACTGGCAAACTCACAC'
 for 'GCGAAAAAAATTAAAGCGCAAGATTGTTGGTTTTTGCGTGATGGTGAC'

150, -45, 3
TAATCAAAACCAATAAACACGAAATAATCCCCATGCCGGTGAAGAAGGGGCGTGACTTTAGCGAAATGTT
+GCCGTCGCGACAACCGGAATATGAAAGCAAAGCGCAGCGTCTGAATAACGTTTATGCTGAAAGCGGATG
+AATAAGGAGATGCG
                                                                      
+                                                                     
+           GCGAAAAAAATTAAAGCGCAAGATTGTTGGTTTTTGCGTGATGGTGAC

No match found in 'TTGTCAGCGAAAAAAATTAAAGCGCAAGATTGTTGGTTTTTGCGTGATGGT
+GACCGGGCAGCCTAAAGGCTATCCTTAACCAGGGAGCTGATT'
 for 'CATTACATTGCTGGATAAGAATGTTTTAGCAATCTCTTTCTGTCATGA'

0, 18, 30
                  GAATGTTTTAGCAATCTCTTTCTGTCATGAATCCATGGCAGTGACCATACTA
+ATGGTGACTGCCATTGATGGAGGGAGACACAGTGCACTGGCAAACTCACAC
CATTACATTGCTGGATAAGAATGTTTTAGCAATCTCTTTCTGTCATGA

No match found in 'TAATCAAAACCAATAAACACGAAATAATCCCCATGCCGGTGAAGAAGGGGC
+GTGACTTTAGCGAAATGTTGCCGTCGCGACAACCGGAATATGAAAGCAAAGCGCAGCGTCTGAATAACG
+TTTATGCTGAAAGCGGATGAATAAGGAGATGCG'
 for 'CATTACATTGCTGGATAAGAATGTTTTAGCAATCTCTTTCTGTCATGA'

No match found in 'TTGTCAGCGAAAAAAATTAAAGCGCAAGATTGTTGGTTTTTGCGTGATGGT
+GACCGGGCAGCCTAAAGGCTATCCTTAACCAGGGAGCTGATT'
 for 'CGCGACAACCGGAATATGAAAGCAAAGCGCAGCGTCTGAATAACGTTT'

No match found in 'GAATGTTTTAGCAATCTCTTTCTGTCATGAATCCATGGCAGTGACCATACT
+AATGGTGACTGCCATTGATGGAGGGAGACACAGTGCACTGGCAAACTCACAC'
 for 'CGCGACAACCGGAATATGAAAGCAAAGCGCAGCGTCTGAATAACGTTT'

75, 0, 48
TAATCAAAACCAATAAACACGAAATAATCCCCATGCCGGTGAAGAAGGGGCGTGACTTTAGCGAAATGTT
+GCCGTCGCGACAACCGGAATATGAAAGCAAAGCGCAGCGTCTGAATAACGTTTATGCTGAAAGCGGATG
+AATAAGGAGATGCG
                                                                      
+     CGCGACAACCGGAATATGAAAGCAAAGCGCAGCGTCTGAATAACGTTT

        [download]
        Examine what is said, not who speaks -- Silence betokens consent -- Love the truth but pardon error.
        Lingua non convalesco, consenesco et abolesco. -- Rule 1 has a caveat! -- Who broke the cabal?
        "Science is about questioning the status quo. Questioning authority".
        In the absence of evidence, opinion is indistinguishable from prejudice.
[reply]
[d/l]
      Dear BrowserUK, Mreece and Grandfather,

      Many thanks for your wise offerings.

      Sam

[reply]
Re: Searching for a word that may only exist in part
by mreece (Friar) on Oct 19, 2006 at 01:13 UTC
    i don't see any way to do it besides eating letters from the back, then starting over from the front.

    here is one way, though index() would surely be faster than m// here. 

    my $sequence = "GAATGTTTTAGCAATCTCTTTCTGTCATGAATCCATGGCAGTGACCATACTAAT
+GGTGACTGCCATTGATGGAGGGAGACACA";
my $find = "CTGGATAAGAATGTTTTAGCAATCTCTT";

my $found;

MATCH: {
    my $tail = $find;
    while ( length($tail) > 2 and not $found ) {
        ($found) = $sequence =~ /($tail)/  # find match
            or substr( $tail, 0, 1, '');   # or eat first letter

    }
    last MATCH if $found;

    my $head = $find;
    ## can chop first since exact match already failed
    while ( chop $head and length($head) > 2 and not $found ) {
        ($found) = $sequence =~ /($head)/;
    }
}

print "found? $found\n";

    [download]
    updated: to provide better(?) var names
[reply]
[d/l]

      For:

      my $sequence = "111...1111...11";
my $find = "11111";

      [download]

      Prints:

      found? 1111

      [download]

      whereas the OP says "if I do not find the whole word within a sequence and start truncating the word, then it can only match at either end of a sequence and not within".

      DWIM is Perl's answer to Gödel
[reply]
[d/l]
[select]
        oh, i misunderstood, thinking only at the beginning or end of the 'find' sequence (1111 is the beginning (..and end..) of 11111).

        the regexes are fixable (/(^$tail|$tail$)/ and same for $head?) easily enough.. i just wanted an excuse to use chop ;-)

[reply]
[d/l]
[select]
Re: Searching for a word that may only exist in part
by GrandFather (Saint) on Oct 18, 2006 at 22:06 UTC

    Are there constraints on how much can be missing? Do you expect that the part that remains will be unique? Will the part that matches match exactly (it's not a fuzzy match)? If it's not unique I presume that a longest match is a best match. Are there any other factors that may help constrain the search?

    DWIM is Perl's answer to Gödel
[reply]
      Good question, I forgot to add that if I do not find the whole word within a sequence and start truncating the word, then it can only match at either end of a sequence and not within.

      In other words, I am looking for an exact match. The only constraint is that I'd want at least 3 consecutive letters to match.

      Thanks
      Sam

[reply]

Back to Seekers of Perl Wisdom