在Perl中使用较少的代码行实现此算法

Question

我想在Perl中实现此算法。 让我们接受的是：

• DNA1 = GACTAGGC
• DNA2 = AGCTAGGA

DNA1的

第一要素是G，我们会发现，如果没有为G在DNA2并将其指向带点。我们继续它直到结束，因此图像显示相同的元素交点为点。

下一步是：连接点。要指向点首先应该在一个小方格的左上角，另一个在右下方（我的意思是线条应该有135度）如果严格度是2，这意味着拒绝从2出现的线和少于2个点（这意味着如果严格度为3，则图像上只会有一条线）。

最后一步是：wordcount。如果wordcount是1（它是一个图像），这意味着比较元素一个接一个。如果是3，则表示将它们中的3个比较。您可以编写WORDCOUNT是1，因为它始终是1

我搜索关于它的程序，这是我所：

$infile1 = "DNA1.txt"; 
$infile2 = "DNA2.txt"; 

$outfile = "plot.txt"; 
$wordsize=0; 
$stringency=0; 

open inf, $infile1 or die "STOP! File $infile1 not found.\n"; 
$sequence1=<inf>; 
chomp $sequence1; 
@seq1=split //,$sequence1; 
close inf; 

open inf, $infile2 or die "STOP! File $infile2 not found.\n"; 
$sequence2=<inf>; 
chomp $sequence2; 
@seq2=split //,$sequence2; 
close inf; 

$Lseq1=$#seq1+1; 
$Lseq2=$#seq2+1; 

open ouf, ">$outfile"; 

for ($i=0;$i<$Lseq1;$i++){ 
print ouf "\n"; 
for ($j=0;$j<$Lseq2;$j++){ 
    $match=0; 
    for ($w=0;$w<=$wordsize;$w++){ 
    if($seq1[$i+$w] eq $seq2[$j+$w]){ 
     $match++; 
    } 
    } 
    if($match > $stringency){ 
    print ouf "1"; 
    } 
    else{ 
    print ouf "0"; 
    } 
} 
} 

您可以检查有关错误，我怎么可以优化我的代码Perl中的代码更少？

PS：我认为可以每次接受$ wordsize等于$严格。

编辑1：我编辑了我的代码，它适用于只是点。

编辑2：算法是这样的：

qseq, sseq = sequences 
win = number of elements to compare for each point 
Strig = number of matches required for a point 

for each q in qseq: 
    for each s in sseq: 
    if CompareWindow(qseq[q:q+win], s[s:s+win], strig): 
     AddDot(q, s) 

EDIT 3：这里是一个更好的算法建议：

osl.iu.edu/~chemuell/projects/bioinf/dotplot.ppt

任何想法，可根据该改进代码更好的算法是可喜的。

Answer 1

首先，最里面的for循环是完全没有必要的。摆脱它会加速你的代码。

其次，您的代码不超过0.1

修复等$严谨的工作：

use strict; 
use warnings; 

my $s1 = 'GACTAGGC'; 
my $s2 = 'AGCTAGGA'; 
my $stringency = 0; 

my @s1 = split //, $s1; 
my @s2 = split //, $s2; 
my @L; 
for my $i (0..$#s1) { 
    for my $j (0..$#s2) { 
     if ($s1[$i] ne $s2[$j]) { 
     $L[$i][$j] = 0; 
     } elsif ($i == 0 || $j == 0) { 
     $L[$i][$j] = 1; 
     } else { 
     $L[$i][$j] = $L[$i-1][$j-1] + 1; 
     } 

     print $L[$i][$j] <= $stringency ? "0" : "1"; 
    } 
    print("\n"); 
} 

现在我们有一个高效的算法，我们可以调整的实施。

use strict; 
use warnings; 

my $s1 = 'GACTAGGC'; 
my $s2 = 'AGCTAGGA'; 
my $stringency = 0; 

my @s1 = split //, $s1; 
my @s2 = split //, $s2; 
my @L = (0) x @s2; 
for my $i (0..$#s1) { 
    for my $j (0..$#s2) { 
     if ($s1[$i] eq $s2[$j]) { 
     ++$L[$j]; 
     } else { 
     $L[$j] = 0; 
     } 

     print $L[$j] <= $stringency ? "0" : "1"; 
    } 

    print("\n"); 
    pop @L; 
    unshift @L, 0; 
} 

如果你想发生了什么更好的主意，更换

print $L[$j] <= $stringency ? "0" : "1"; 

与

print $L[$j]; 

你就会得到这样

01000110 
10001002 
00100000 
00020000 
10003001 
02000410 
01000150 
00200000 

顺便说一句，试图实现的是非常相似的o找到longest common substring。

更新要想从文件$s1和$s2，在一次一行，

open(my $fh1, '<', ...) or die(...); 
open(my $fh2, '<', ...) or die(...); 

for (;;) { 
    my $s1 = <$fh1>; 
    my $s2 = <$fh2>; 

    die("Files have different length\n") 
     if defined($s1) && !defined($s2) 
     || !defined($s1) && defined($s2); 

    last if !defined(($s1); 

    chomp($s1); 
    chomp($s2); 

    ... code to generate graph ... 
}