Dear Monks,

I have 2 DNA sequences and want to perform a dotplot using a sliding window of a particular size with particular mismatch allowed. The following script produces a text matrix in which "1" stands for a hit.

The problem is, that this code works very slowly with long DNA sequences (>10000). I already use my four CPU cores for calculation but I've no idea how I could further speed it up. Any suggestions are welcome!

### CREATE TWO SAMPLE DNA SEQUENCES ###
@nucleotides=('A','T','G','C');
foreach(1..1000)
    {
    $seq1.=$nucleotides[int(rand(3))];
    }
foreach(1..1000)
    {
    $seq2.=$nucleotides[int(rand(3))];
    }


### SETTINGS FOR THE DOTPLOT ###
$window_size=5;
$max_mismatch=1;


@seq1=split('',$seq1);
foreach(1..$window_size)
    {
    shift@seq1;
    }
@seq2=split('',$seq2);
foreach(1..$window_size)
    {
    shift@seq2;
    }

open(OUT,">ID_matrix.txt");
$number_of_windows_1=((length$seq1)-$window_size)+1;
$number_of_windows_2=((length$seq2)-$window_size)+1;

$time_start=time;
foreach$window_no(0..$number_of_windows_1-1)
    {
    @seq2_temp=@seq2;
    if($window_no==0)
        {
        $current_window=substr($seq1,0,$window_size);
        @current_window=split('',$current_window);
        }
    else
        {
        shift@current_window;
        $next_character=shift@seq1;
        push(@current_window,$next_character);
        }
    foreach$query_no(0..$number_of_windows_2-1)
        {
        if($query_no==0)
            {
            $query_window=substr($seq2,0,$window_size);
            @query_window=split('',$query_window);
            }
        else
            {
            shift@query_window;
            $next_character=shift@seq2_temp;
            push(@query_window,$next_character);
            }
        $count_matches=0;
        foreach(0..$window_size-1)
            {
            if($current_window[$_]eq$query_window[$_])
                {
                $count_matches++;
                last if($_-$count_matches>$max_mismatch);
                }
            }
        if($count_matches>=$window_size-$max_mismatch)
            {
            print OUT "1";
            }
        else
            {
            print OUT "0";
            }
        }
    print OUT"\n";
    }
$time_end=time;
$time_used=$time_end-$time_start;
close OUT;
print"Time used: $time_used seconds.\n";
system("pause");
exit;

[download]

In reply to Speed up DNA dotplot by Microcebus

Title:
Use:  <p> text here (a paragraph) </p>
and:  <code> code here </code>
to format your post, it's "PerlMonks-approved HTML":


  • Posts are HTML formatted. Put <p> </p> tags around your paragraphs. Put <code> </code> tags around your code and data!
  • Titles consisting of a single word are discouraged, and in most cases are disallowed outright.
  • Read Where should I post X? if you're not absolutely sure you're posting in the right place.
  • Please read these before you post! —
  • Posts may use any of the Perl Monks Approved HTML tags:
    a, abbr, b, big, blockquote, br, caption, center, col, colgroup, dd, del, details, div, dl, dt, em, font, h1, h2, h3, h4, h5, h6, hr, i, ins, li, ol, p, pre, readmore, small, span, spoiler, strike, strong, sub, summary, sup, table, tbody, td, tfoot, th, thead, tr, tt, u, ul, wbr
  • You may need to use entities for some characters, as follows. (Exception: Within code tags, you can put the characters literally.)
            For:     Use:
    & &amp;
    < &lt;
    > &gt;
    [ &#91;
    ] &#93;
  • Link using PerlMonks shortcuts! What shortcuts can I use for linking?
  • See Writeup Formatting Tips and other pages linked from there for more info.