Нахождение самой большой открытой рамки для чтения с помощью Perl

Question

Я создал программу, которая может читать последовательности ДНК, которые могут генерировать комплементарную цепь, которая в дальнейшем транслируется в мРНК.Однако мне нужно найти максимально длинную открытую рамку считывания для этой ДНК.Я что-то закодировал, но когда он распечатывает заявление, я не получаю ответа.Помощь?

Это то, что у меня есть;

# Search for the longest open reading frame for this DNA.
print "\nHere is the largest ORF, from 5' to 3':\n" ;
local $_ = $RNA_seq ;
while ( / AUG /g ) {
    my $start = pos () - 2 ;
    if ( / UGA|UAA|UAG /g ) {
        my $stop = pos ;
        $gene = substr ( $_ , $start - 1 , $stop - $start + 1 ), $/ ;
        print "$gene" ;
    }
}

# The next set of commands translates the ORF found above for an amino acid seq.
print "\nThe largest reading Frame is:\t\t\t" . $protein { "gene" } . "\n" ;
sub translate {
    my ( $gene , $reading_frame ) = @_ ;
    my %protein = ();
    for ( $i = $reading_frame ; $i < length ( $gene ); $i += 3 ) {
        $codon = substr ( $gene , $i , 3 );
        $amino_acid = translate_codon( $codon );
        $protein { $amino_acid }++;
        $protein { "gene" } .= $amino_acid ;
    }
    return %protein ;
}

sub translate_codon {
if ( $_ [ 0 ] =~ / GC[AGCU] /i )             { return A;} # Alanine;
if ( $_ [ 0 ] =~ / UGC|UGU /i )              { return C;} # Cysteine
if ( $_ [ 0 ] =~ / GAC|GAU /i )              { return D;} # Aspartic Acid;
if ( $_ [ 0 ] =~ / GAA|GAG /i )              { return Q;} # Glutamine;
if ( $_ [ 0 ] =~ / UUC|UUU /i )              { return F;} # Phenylalanine;
if ( $_ [ 0 ] =~ / GG[AGCU] /i )             { return G;} # Glycine;
if ( $_ [ 0 ] =~ / CAC|CAU /i )              { return His;} # Histine (start codon);
if ( $_ [ 0 ] =~ / AU[AUC] /i )              { return I;} # Isoleucine;
if ( $_ [ 0 ] =~ / AAA|AAG /i )              { return K;} # Lysine;
if ( $_ [ 0 ] =~ / UUA|UUG|CU[AGCU] /i )     { return Leu;} # Leucine;
if ( $_ [ 0 ] =~ / AUG /i )                  { return M;} # Methionine;
if ( $_ [ 0 ] =~ / AAC|AAU /i )              { return N;} # Asparagine;
if ( $_ [ 0 ] =~ / CC[AGCU] /i )             { return P;} # Proline;
if ( $_ [ 0 ] =~ / CAA|CAG /i )              { return G;} # Glutamine;
if ( $_ [ 0 ] =~ / AGA|AGG|CG[AGCU] /i )     { return R;} # Arginine;
if ( $_ [ 0 ] =~ / AGC|AGU|UC[AGCU] /i )     { return S;} # Serine;
if ( $_ [ 0 ] =~ / AC[AGCU] /i )             { return T;} # Threonine;
if ( $_ [ 0 ] =~ / GU[AGCU] /i )             { return V;} # Valine;
if ( $_ [ 0 ] =~ / UGG /i )                  { return W;} # Tryptophan;
if ( $_ [ 0 ] =~ / UAC|UAU /i )              { return Y;} # Tyrosine;
if ( $_ [ 0 ] =~ / UAA|UGA|UAG /i )          { return "***" ;} # Stop Codons;
}

Я что-то упустил?

Answer 1

Поставьте

use strict;
use warnings;

в начале вашего кода: это поможет обнаружить проблемы

, как правило, с

return A;

вы не возвращает строку, но дескриптор файла, если вы хотите вернуть строку

return 'A';

Answer 2

Если вы хотите поместить пробелы в регулярные выражения для удобства чтения, вам нужно использовать модификатор регулярного выражения x как:

if ( $_ [ 0 ] =~ / GC[AGCU] /xi )

Тогда, если вам нужно представить пробел (чего здесь не следует), просто используйте \s В качестве отступления вместо использования $_ [0] рассмотрим:

sub translate_codon {
    my $codon = shift;
    return q(A) if $codon =~ m/ GC[AGCU] /xi;
    return q(F) if $codon =~ m/ UU[CU]   /xi;
}

Answer 3

Вы уже проверили BioPerl ?Я не проверял, но, возможно, он уже включает в себя функции, которые вам требуются.Или вы сейчас работаете над заданием, для которого вам нужно написать программу самостоятельно?

РЕДАКТИРОВАТЬ

Я не совсем уверен в первом разделе кодаВы отправили.Например, в ваших регулярных выражениях есть пробелы.Включает ли строка, которую вы пытаетесь сопоставить, эти пробелы или все кодоны записаны вместе, как в AUGCCGGAUGA?В последнем случае просто не будет совпадения, даже если присутствуют искомые кодоны (я, вероятно, говорю вам то, что вы знаете, хотя ... Я просто хотел указать на это, на всякий случай :)).

Кроме того, каков код вашей pos функции?

Еще одна вещь: вам не нужно устанавливать $_, вы можете просто сопоставить $RNA_seq сшаблон, вот так:

if ($RNA_seq =~ m/UGA/) { # ...

РЕДАКТИРОВАТЬ 2

Я подумал немного больше о первом разделе, я думаю, что использование функции index будетздесь рекомендуется, так как это сразу дает вам позицию в последовательности:

#!/usr/bin/perl

use strict;
use warnings;
use List::Util qw( min );

my $string = 'UGAAUGGGCCAUAUUUGACUGAGUUUCAGAUGCCAUUGGCGAUUAG';
# the genes:     *-------------*           *---------------*

my $prev = -1;
my @genes;

while (1) {
   my $start = index($string, 'AUG', $prev+1);
   my $stop  = min grep $_ > -1, (index($string, 'UGA', $start+1), index($string, 'UAA', $start+1),
                                  index($string, 'UAG', $start+1));

   # exit the loop if index is -1, i.e. if there was no more match
   last if ($start < 0 or $stop < 0); 

   # using +1 so that 'AUGA' will not count as a gene
   if ($stop > $start+1) { 
      push @genes, substr($string, $start, $stop);
   }

   $prev = $stop; # I'm assuming that no second AUG will come between AUG and a stop codon
}

print "@genes\n";

Это дает вам

AUGGGCCAUAUUUGA AUGCCAUUGGCGAUUAG

Я бы сказал, что, возможно, потребуется некоторое уточнение, но я надеюсь, что это будетбудь полезным все равно.

Answer 4

Есть несколько приложений для ab-initio предсказаний гена . Если вы все еще хотите закодировать его с нуля, я рекомендую изучить динамическое программирование и алгоритм Смита-Уотермана для локального выравнивания последовательностей . Однако обратите внимание, что для эукариотных организмов вы также должны учитывать сплайсинг РНК .