Wie kann ich Start- und Endcodon aus DNA-Sequenzen in Perl extrahieren?

Question

Ich habe einen Code unten, die versuchen, die Position der Start-und End-Codon der gegebenen DNA-Sequenzen zu identifizieren. Wir definieren Startcodon als a ATG Reihenfolge und Endcodon als TGA,TAA,TAG Sequenzen.

Das Problem, das ich habe, ist, dass der nachstehende Code nur für die ersten beiden Sequenzen (DM208659 und AF038953) funktioniert, nicht aber für den Rest.

Was ist an meinem nachstehenden Ansatz falsch?

Dieser Code kann kopiert und eingefügt werden aus aquí .

      #!/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

Ich habe die Verwendung von $_ (Mich schauderte es besonders, als Sie local Sie haben es richtig gemacht, aber warum sollten Sie sich Gedanken darüber machen, ob eine andere Funktion Sie in den Ruin treiben könnte? $_ statt der Verwendung von $rna_sq die bereits verfügbar ist?

Zusätzlich habe ich korrigiert $start y $stop zu 0-basierten Indizes in der Zeichenkette zu machen (was den Rest der Mathematik einfacher machte), und berechnete $genelen früh, so dass es direkt in der Datenbank verwendet werden kann. substr Betrieb. (Alternativ können Sie auch $[ auf 1 setzen, um 1-basierte Array-Indizes zu verwenden, siehe 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

Es geht darum, nie aus der inneren Schleife auszubrechen, wenn die if (/tga|taa|tag/g) ein Endcodon nicht findet. Es passt immer wieder /atg/g wiederholt, ohne weiter voranzukommen. Sie könnten ihn gewaltsam aus der inneren Schleife auswerfen:

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

Answer 3

Es kommt darauf an, ob Sie Sequenzen erzeugen wollen, die sich überschneiden könnten. Die Sequenz AF038954 enthält zum Beispiel atgaccatcctccagacatacttccggcagaacagggatga , dessen Ende sich überschneidet mit atgaagtcttggacaacctcttggcttttgtctgtga . Möchten Sie sie beide melden?

Wenn Sie keine sich überschneidenden Sequenzen melden wollen, ist dies ein sehr einfaches Problem, das Sie mit einer einzigen Regexp lösen können:

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);
      }
}

Der Regexp (atg.*?(?:tga|taa|tag)) Ihrem gewünschten Start entspricht, dann so wenig wie möglich von dem, was danach kommt (das ist die ? zum Anhalten der .* "gierig" zu sein), dann ist das gewünschte Ergebnis erreicht. Die Iteration darüber mit dem while Schleifen-Neustarts nach diese Übereinstimmung, die die Anforderung erfüllt, nicht nach Überschneidungen zu suchen.

Wenn Sie tun Wenn Sie überlappende Sequenzen melden möchten, benötigen Sie einen zweistufigen Prozess: Finden Sie den Anfang, finden Sie das Ende und finden Sie dann einen weiteren Anfang, indem Sie dort weitermachen, wo Sie bei der letzten Suche aufgehört haben. Mit einer zweiten Regexp lässt sich die Aufgabe jedoch einfacher lösen:

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);
      }
    }
}

Hier verwenden wir die (im Allgemeinen nicht empfohlene) $' eine spezielle Variable, die den Inhalt nach der Übereinstimmung enthält. Darin suchen wir nach dem Ende der Sequenz und geben die Details aus. Da unser globaler Hauptabgleich gegen $rna_seq die Sequenz nicht einschließt (wie oben), beginnen wir die Suche an der Stelle neu, an der die vorherige Suche aufgehört hat, d. h. direkt nach dem gefundenen Anfang. Auf diese Weise können wir tun überlappende Sequenzen enthalten.