Как извлечь начальный и конечный кодоны из последовательностей ДНК в Perl?

Question

У меня есть код ниже, который пытается определить положение начального и конечного кодонов данных последовательностей ДНК. Мы определяем стартовый кодон как a ATG последовательность и конечный кодон как TGA, TAA, TAG последовательности.

У меня проблема в том, что приведенный ниже код работает только для первых двух последовательностей (DM208659 и AF038953), но не для остальных.

Что не так с моим подходом ниже?

Этот код можно скопировать с здесь .

      #!/usr/bin/perl -w


while (<DATA>) {
    chomp;
    print "$_\n";
    my ($id,$rna_sq) = split(/\s+/,$_);

    local $_ = $rna_sq;
    while (/atg/g) {
        my $start = pos() - 2;

        if (/tga|taa|tag/g) {

            my $stop    = pos();
            my $gene    = substr( $_, $start - 1, $stop - $start + 1 ),$/;
            my $genelen = length($gene);
            my $ct      = "$id $start $stop $gene $genelen";
            print "\t$ct\n";

        }

    }

}

__DATA__
DM208659    gtgggcctcaaatgtggagcactattctgatgtccaagtggaaagtgctgcgacatttgagcgtcac
AF038953    gatcccagacctcggcttgcagtagtgttagactgaagataaagtaagtgctgtttgggctaacaggatctcctcttgcagtctgcagcccaggacgctgattccagcagcgccttaccgcgcagcccgaagattcactatggtgaaaatcgccttcaatacccctaccgccgtgcaaaaggaggaggcgcggcaagacgtggaggccctcctgagccgcacggtcagaactcagatactgaccggcaaggagctccgagttgccacccaggaaaaagagggctcctctgggagatgtatgcttactctcttaggcctttcattcatcttggcaggacttattgttggtggagcctgcatttacaagtacttcatgcccaagagcaccatttaccgtggagagatgtgcttttttgattctgaggatcctgcaaattcccttcgtggaggagagcctaacttcctgcctgtgactgaggaggctgacattcgtgaggatgacaacattgcaatcattgatgtgcctgtccccagtttctctgatagtgaccctgcagcaattattcatgactttgaaaagggaatgactgcttacctggacttgttgctggggaactgctatctgatgcccctcaatacttctattgttatgcctccaaaaaatctggtagagctctttggcaaactggcgagtggcagatatctgcctcaaacttatgtggttcgagaagacctagttgctgtggaggaaattcgtgatgttagtaaccttggcatctttatttaccaactttgcaataacagaaagtccttccgccttcgtcgcagagacctcttgctgggtttcaacaaacgtgccattgataaatgctggaagattagacacttccccaacgaatttattgttgagaccaagatctgtcaagagtaagaggcaacagatagagtgtccttggtaataagaagtcagagatttacaatatgactttaacattaaggtttatgggatactcaagatatttactcatgcatttactctattgcttatgccgtaaaaaaaaaaaaaaaaaaaaaaaaaaaaa
BC021011    ggggagtccggggcggcgcctggaggcggagccgcccgctgggctaaatggggcagaggccgggaggggtgggggttccccgcgccgcagccatggagcagcttcgcgccgccgcccgtctgcagattgttctg
DM208660    gggatactcaaaatgggggcgctttcctttttgtctgtactgggaagtgcttcgattttggggtgtccc
AF038954    ggacccaagggggccttcgaggtgccttaggccgcttgccttgctctcagaatcgctgccgccatggctagtcagtctcaggggattcagcagctgctgcaggccgagaagcgggcagccgagaaggtgtccgaggcccgcaaaagaaagaaccggaggctgaagcaggccaaagaagaagctcaggctgaaattgaacagtaccgcctgcagagggagaaagaattcaaggccaaggaagctgcggcattgggatcccgtggcagttgcagcactgaagtggagaaggagacccaggagaagatgaccatcctccagacatacttccggcagaacagggatgaagtcttggacaacctcttggcttttgtctgtgacattcggccagaaatccatgaaaactaccgcataaatggatagaagagagaagcacctgtgctgtggagtggcattttagatgccctcacgaatatggaagcttagcacagctctagttacattcttaggagatggccattaaattatttccatatattataagagaggtccttccactttttggagagtagccaatctagctttttggtaacagacttagaaattagcaaagatgtccagctttttaccacagattcctgagggattttagatgggtaaatagagtcagactttgaccaggttttgggcaaagcacatgtatatcagtgtggacttttcctttcttagatctagtttaaaaaaaaaaaccccttaccattctttgaagaaaggaggggattaaataattttttcccctaacactttcttgaaggtcaggggctttatctatgaaaagttagtaaatagttctttgtaacctgtgtgaagcagcagccagccttaaagtagtccattcttgctaatggttagaacagtgaatactagtggaattgtttgggctgcttttagtttctcttaatcaaaattactagatgatagaattcaagaacttgttacatgtattacttggtgtatcgataatcatttaaaagtaaagactctgtcatgcaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa

Answer 1

Я исключил использование $_ (Я особенно вздрогнул, когда вы local измерили его - вы сделали это правильно, но зачем заставлять себя беспокоиться, если какая-то другая функция будет сжимать $_, а не использовать $rna_sq который уже доступен?

Кроме того, я исправил $start и $stop, чтобы в строке были индексы, основанные на 0 (что сделало остальную часть математики более простой), и вычислил $genelen рано, чтобы его можно было использовать непосредственно в substr операция. (В качестве альтернативы вы можете локализовать $[ в 1, чтобы использовать индексы массивов на основе 1, см. perldoc perlvar .)

use strict;
use warnings;
while (my $line = <DATA>) {
    chomp $line;
    print "processing $line\n";
    my ($id, $rna_sq) = split(/\s+/, $line);

    while ($rna_sq =~ /atg/g) {
        # $start and $stop are 0-based indexes
        my $start = pos($rna_sq) - 3; # back up to include the start sequence

        # discard remnant if no stop sequence can be found
        last unless $rna_sq =~ /tga|taa|tag/g;

        my $stop    = pos($rna_sq);
        my $genelen = $stop - $start;

        my $gene    = substr($rna_sq, $start, $genelen);
        print "\t" . join(' ', $id, $start+1, $stop, $gene, $genelen) . "\n";
    }
}

Answer 2

Все зависит от того, хотите ли вы генерировать последовательности, которые могут перекрываться. Например, последовательность AF038954 содержит atgaccatcctccagacatacttccggcagaacagggatga, конец которой перекрывается с atgaagtcttggacaacctcttggcttttgtctgtga. Вы хотите сообщить о них обоих?

Если вы не хотите сообщать о последовательностях, которые перекрываются, это очень простая проблема, которую вы можете решить с помощью одного регулярного выражения:

while (<DATA>) {
    chomp;
    print "processing $_\n";
    my ($id, $rna_sq) = split;

    while ($rna_sq =~ /(atg.*?(?:tga|taa|tag))/g) {
      printf "\t%8s %4i %4i %s %i\n",
             $id,
             pos($rna_sq) - length($1) + 1,
             pos($rna_sq),
             $1,
             length($1);
      }
}

Регулярное выражение (atg.*?(?:tga|taa|tag)) соответствует вашему необходимому началу, затем как можно меньше из того, что будет дальше (это ?, чтобы остановить .* быть "жадным"), затем к вашему требуемому концу. Итерирование по нему с помощью цикла while перезапускает после этого совпадения, что отвечает требованию не искать перекрытия.

Если вы хотите хотите, чтобы сообщалось о перекрывающихся последовательностях, вам нужен двухэтапный процесс: найдите начало, найдите конец, а затем найдите другой старт, подбирая с того места, где вы остановились, в поисках начала. последний раз. Но вы все равно можете выполнить более простую работу, используя второе регулярное выражение:

while (<DATA>) {
    chomp;
    print "processing $_\n";
    my ($id, $rna_sq) = split;

    while ($rna_sq =~ /atg/g) {
      if ($' =~ /(.*?(?:tga|taa|tag))/) {
        my $match = "atg$1";
        printf "\t%8s %4i %4i %s %i\n",
               $id,
               pos($rna_sq) - 2,
               pos($rna_sq) - 3 + length($match),
               $match,
               length($match);
      }
    }
}

Здесь мы используем (обычно не рекомендуется) $' специальную переменную, которая содержит содержимое после совпадения. Мы смотрим в этом, чтобы найти конец последовательности и вывести детали. Поскольку наше основное глобальное совпадение с $rna_seq не включает в себя последовательность (как и выше), мы возобновляем поиск с начала, на котором прервался предыдущий поиск, то есть сразу после найденного начала. Таким образом, мы делаем включаем перекрывающиеся последовательности.

Answer 3

Это никогда не выходит из вашего внутреннего цикла, когда if (/tga|taa|tag/g) не может найти конечный кодон. Он постоянно соответствует /atg/g, никогда не продвигаясь дальше. Вы можете принудительно извлечь его из внутреннего цикла:

if (/tga|taa|tag/g) {
    ...
}
else {
    last;
}