如何從Perl中的DNA序列提取起始和終止密碼子？

Question

我有一個代碼，試圖確定給定DNA序列的起始密碼子和終止密碼子的位置。 我們定義起始密碼子爲一個ATG序列和端密碼子爲TGA，TAA，TAG序列。

我遇到的問題是下面的代碼僅適用於前兩個序列（DM208659和AF038953），但其餘部分不適用。

我的方法出現了什麼問題？

該代碼可以從here複製粘貼。

 #!/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操作中直接使用ernatively，你可以本地化$[ 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

當if (/tga|taa|tag/g)未能找到結尾密碼子時，它永遠不會突破你的內部循環。它會一直保持匹配/atg/g，永遠不會前進。你可以強行從內環彈出：

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

Answer 3

這一切都取決於你是否想生成可能重疊的序列。例如，序列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的主要全局匹配不包含序列（如上所述），因此我們重新開始搜索以前搜索停止的開始位置，即在我們找到的開始位置之後。這樣我們do包括重疊的序列。