4

我在下面有一个代码,试图识别给定 DNA 序列的起始和结束密码子的位置。我们将起始密码子定义为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
4

3 回答 3

4

我删除了$_( 当你定义它时我特别不寒而栗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";
    }
}
于 2009-10-13T04:54:27.600 回答
1

if (/tga|taa|tag/g)当找不到结束密码子时,它永远不会脱离你的内循环。它一直在重复匹配/atg/g,从不进一步前进。您可以将其从内循环中强行弹出:

if (/tga|taa|tag/g) {
    ...
}
else {
    last;
}
于 2009-10-13T04:14:56.420 回答
1

这完全取决于您是否要生成可能重叠的序列。例如,序列 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不包括序列(如上所示),所以我们重新开始搜索之前搜索停止的开始,即在我们找到的开始之后。通过这种方式,我们确实包含了重叠序列。

于 2009-10-13T08:57:32.387 回答