comment on

In bioinformatics, it is often important to determine if a particular sequence is in part of the DNA just sequenced. However, since DNA sequencers are imperfect devices, exact matches to the sequence being searched for is rare. Tools for similar sorts of problems are available, noteably BLAST, but it can be over kill to run the BLAST, parse the results, and then make a desision based on the results. So I wrode a script that will take a sequence being searched for (typically 20 to 30 bases long) and write a new perl script that contains a series of regexes that will allow exactly one mismatch, deletion or insertion.

Here is some example input sequences:

>test1
ATGCATGCATGCATGC
>test2
ATGCNATGCNATGCNATGCN
[download]

which produces this output:

#!/usr/bin/perl -w
$seq = $ARGV[0];
if (
# test1
  $seq =~ /\w{0,2}TGCATGCATGCATGC/ or
  $seq =~ /A\w{0,2}GCATGCATGCATGC/ or
  $seq =~ /AT\w{0,2}CATGCATGCATGC/ or
  $seq =~ /ATG\w{0,2}ATGCATGCATGC/ or
  $seq =~ /ATGC\w{0,2}TGCATGCATGC/ or
  $seq =~ /ATGCA\w{0,2}GCATGCATGC/ or
  $seq =~ /ATGCAT\w{0,2}CATGCATGC/ or
  $seq =~ /ATGCATG\w{0,2}ATGCATGC/ or
  $seq =~ /ATGCATGC\w{0,2}TGCATGC/ or
  $seq =~ /ATGCATGCA\w{0,2}GCATGC/ or
  $seq =~ /ATGCATGCAT\w{0,2}CATGC/ or
  $seq =~ /ATGCATGCATG\w{0,2}ATGC/ or
  $seq =~ /ATGCATGCATGC\w{0,2}TGC/ or
  $seq =~ /ATGCATGCATGCA\w{0,2}GC/ or
  $seq =~ /ATGCATGCATGCAT\w{0,2}C/ or
  $seq =~ /ATGCATGCATGCATG\w{0,2}/ or
  $seq =~ /ATGCATGCATGCATGC\w{0,2}/ or
  $seq =~ /ZZZ/) {
    $match=$&;
    $name = "test1";
    $matchlength=length($match);
    $firstpos=index($seq,$match);
    print "$name, $match, $matchlength, $firstpos \n";
    print "Got it!\n";
} elsif (
# test2
  $seq =~ /\w{0,2}TGC[ATGCN]ATGC[ATGCN]ATGC[ATGCN]ATGC[ATGCN]/ or
  $seq =~ /A\w{0,2}GC[ATGCN]ATGC[ATGCN]ATGC[ATGCN]ATGC[ATGCN]/ or
  $seq =~ /AT\w{0,2}C[ATGCN]ATGC[ATGCN]ATGC[ATGCN]ATGC[ATGCN]/ or
  $seq =~ /ATG\w{0,2}[ATGCN]ATGC[ATGCN]ATGC[ATGCN]ATGC[ATGCN]/ or
  $seq =~ /ATGC\w{0,2}ATGC[ATGCN]ATGC[ATGCN]ATGC[ATGCN]/ or
  $seq =~ /ATGC[ATGCN]\w{0,2}TGC[ATGCN]ATGC[ATGCN]ATGC[ATGCN]/ or
  $seq =~ /ATGC[ATGCN]A\w{0,2}GC[ATGCN]ATGC[ATGCN]ATGC[ATGCN]/ or
  $seq =~ /ATGC[ATGCN]AT\w{0,2}C[ATGCN]ATGC[ATGCN]ATGC[ATGCN]/ or
  $seq =~ /ATGC[ATGCN]ATG\w{0,2}[ATGCN]ATGC[ATGCN]ATGC[ATGCN]/ or
  $seq =~ /ATGC[ATGCN]ATGC\w{0,2}ATGC[ATGCN]ATGC[ATGCN]/ or
  $seq =~ /ATGC[ATGCN]ATGC[ATGCN]\w{0,2}TGC[ATGCN]ATGC[ATGCN]/ or
  $seq =~ /ATGC[ATGCN]ATGC[ATGCN]A\w{0,2}GC[ATGCN]ATGC[ATGCN]/ or
  $seq =~ /ATGC[ATGCN]ATGC[ATGCN]AT\w{0,2}C[ATGCN]ATGC[ATGCN]/ or
  $seq =~ /ATGC[ATGCN]ATGC[ATGCN]ATG\w{0,2}[ATGCN]ATGC[ATGCN]/ or
  $seq =~ /ATGC[ATGCN]ATGC[ATGCN]ATGC\w{0,2}ATGC[ATGCN]/ or
  $seq =~ /ATGC[ATGCN]ATGC[ATGCN]ATGC[ATGCN]\w{0,2}TGC[ATGCN]/ or
  $seq =~ /ATGC[ATGCN]ATGC[ATGCN]ATGC[ATGCN]A\w{0,2}GC[ATGCN]/ or
  $seq =~ /ATGC[ATGCN]ATGC[ATGCN]ATGC[ATGCN]AT\w{0,2}C[ATGCN]/ or
  $seq =~ /ATGC[ATGCN]ATGC[ATGCN]ATGC[ATGCN]ATG\w{0,2}[ATGCN]/ or
  $seq =~ /ATGC[ATGCN]ATGC[ATGCN]ATGC[ATGCN]ATGC\w{0,2}/ or
  $seq =~ /ATGC[ATGCN]ATGC[ATGCN]ATGC[ATGCN]ATGC[ATGCN]\w{0,2}/ or
  $seq =~ /ZZZ/) {
    $match=$&;
    $name = "test2";
    $matchlength=length($match);
    $firstpos=index($seq,$match);
    print "$name, $match, $matchlength, $firstpos \n";
    print "Got it!\n";
} elsif (
"false"){
}
[download]

and finally, the code that writes the code:

#!/usr/bin/perl -w
use strict;

my $FILE_IN  = $ARGV[0];
my $FILE_OUT = "patternmatcher.pl";
open FILE_IN,  $FILE_IN;
open FILE_OUT, ">$FILE_OUT";

print FILE_OUT "\#!/usr/bin/perl -w\n";
print FILE_OUT "\$seq = \$ARGV[0]\;\n";
print FILE_OUT "if (\n";

my $wild = "\\w\{0,2\}";
my $onN  = "X";
my $onN2 = "\[ATGCN\]";
my $name;
while (<FILE_IN>) {
    if ( $_ =~ /^\>(.+)$/ ) {
        $name = $1;
    }
    elsif ( $_ =~ /^(\w+)$/ ) {
        my $seq = $1;
        my $len = length $seq;
        print FILE_OUT "\# $name\n";

        for my $i ( 0 .. $len ) { #go through $seq base by base, repla
+cing one character with the wild card
            my $newseq = $seq;
            substr( $newseq, $i, 1 ) = $wild;
            $newseq =~ s/N/$onN/go;  # these to lines are to allow Ns 
+to match any character
            $newseq =~ s/X/$onN2/go; #
            print FILE_OUT "  \$seq =~ \/$newseq\/ or \n";
        }
        print FILE_OUT "  \$seq =~ \/ZZZ\/) \{\n"; #impossible match t
+o close the if
        print FILE_OUT "    \$match=\$\&\;\n";
        print FILE_OUT "    \$name = \"$name\"\;\n";
        print FILE_OUT "    \$matchlength=length\(\$match\)\;\n";
        print FILE_OUT "    \$firstpos=index\(\$seq,\$match)\;\n";
        print FILE_OUT
          "    print \"\$name, \$match, \$matchlength, \$firstpos \\n\
+"\;\n";
        print FILE_OUT "    print \"Got it!\\n\"\;\n";
        print FILE_OUT "\} elsif \(\n";

    }
}
print FILE_OUT "\"false\"\)\{\n"; # false value to end the series of e
+lsifs

print FILE_OUT "\}\n";
close FILE_IN;
close FILE_OUT;

system( "chmod", "u+x", $FILE_OUT );
[download]

In reply to An imperfect pattern matcher writer by scain

Are you posting in the right place? Check out Where do I post X? to know for sure.
Posts may use any of the Perl Monks Approved HTML tags. Currently these include the following:
<code> <a> <b> <big> <blockquote> <br /> <dd> <dl> <dt> <em> <font> <h1> <h2> <h3> <h4> <h5> <h6> <hr /> <i> <li> <nbsp> <ol> <p> <small> <strike> <strong> <sub> <sup> <table> <td> <th> <tr> <tt> <u> <ul>
Snippets of code should be wrapped in <code> tags not <pre> tags. In fact, <pre> tags should generally be avoided. If they must be used, extreme care should be taken to ensure that their contents do not have long lines (<70 chars), in order to prevent horizontal scrolling (and possible janitor intervention).
Want more info? How to link or How to display code and escape characters are good places to start.


There's more than one way to do things
	PerlMonks