Реализация этого алгоритма с меньшим количеством строк кода на 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)

РЕДАКТИРОВАТЬ 3: Вот лучшее предложение алгоритма:

osl.iu.edu / ~ chemuell / проекты / bioinf / dotplot.ppt

Любая идея улучшить код в соответствии с этим лучшим алгоритмом приветствуется.

Answer 1

Во-первых, внутренний цикл for совершенно не нужен. Избавление от него ускорит ваш код.

Во-вторых, ваш код не работает на $ строгость, отличную от 0.

Fix:

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

Кстати, то, что мы пытаемся достичь, удивительно похоже на поиск самой длинной общей подстроки .

Обновление Чтобы получить $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 ...
}