Get kürzeste und längste Sequenz in Datei

Question

Ich versuche, die kürzeste und längste Sequenz in einer Datei, die mehrere genbankähnliche Einträge enthält. Beispiel für die Datei:

LOCUS  NM_182854    2912 bp mRNA linear PRI 20-APR-2016 
DEFINITION Homo sapiens mRNA. 
ACCESSION NM_182854 
SOURCE  Homo sapiens (human) 
    ORGANISM Homo sapiens 
      Eukaryota; Metazoa; Chordata; Craniata; Vertebrata; Euteleostomi; 
      Mammalia; Eutheria; Euarchontoglires; Primates; Haplorrhini; 
      Catarrhini; Hominidae; Homo. 

ORIGIN  
     1 gggcgatcag aagcaggtca cacagcctgt ttcctgtttt caaacgggga acttagaaag 
     61 tggcagcccc tcggcttgtc gccggagctg agaaccaaga gctcgaaggg gccatatgac 
     // 

LOCUS  NM_001323410   6992 bp mRNA linear PRI 20-APR-2016 
DEFINITION Homo sapiens mRNA. 
ACCESSION NM_001323410 
SOURCE  Homo sapiens (human) 
    ORGANISM Homo sapiens 
      Eukaryota; Metazoa; Chordata; Craniata; Vertebrata; Euteleostomi; 
      Mammalia; Eutheria; Euarchontoglires; Primates; Haplorrhini; 
      Catarrhini; Hominidae; Homo. 

ORIGIN  
     1 actacttccg gcttccccgc cccgccccgt ccccgggcgt ctccattttg gtctcaggtg 
     61 tggactcggc aagaaccagc gcaagaggga agcagagtta tagctacccc ggc 
     // 

Ich mag die Zugangsnummer, die Art des Organismus, aus der kürzesten Sequenz und die längste Sequenz so weit

mein Code drucken:

#!/usr/bin/perl 

use strict; 
use warnings; 

print "enter file path\n"; 

while (my $line = <>){ 
    chomp $line; 
    my @record = ($line); 

    foreach my $file(@record){ 
    open(IN, "$file") or die "\n error opening file \n;/\n"; 

    $/="//"; 

    while (my $line = <IN>){ 
     my @gb_seq = split ("ORIGIN", $line); 
     my $definition = $gb_seq[0]; 
     my $sequence = $gb_seq[1]; 

     $definition =~ m/ORGANISM[\s\t]+(.+)[\n\s\t]+/; 
     my $organism = $1; 

     if ($definition =~ m/ACCESSION[\s\t]+(\D\D_\d\d\d\d\d\d(\d*))[\n\s\t]+/){ 
     my $accession = $1; 

      $sequence =~ s/\d//g; 
      $sequence =~ s/[\n\s\t]//g; 
      my $size = length($sequence); 
      my @sorted_keys = sort { $a <=> $b } keys my %size; 
      my $shortest = $sorted_keys[0]; 
      my $longest = $sorted_keys[-1]; 

      print "this is the shortest: $accession $organism size: $shortest\n"; 
      print "this is the longest: $accession $organism size: $longest\n"; 
    } 
    }}} 
    exit; 

Ich dachte daran, die Länge in einen Hash zu setzen, um den kürzesten und den längsten zu bekommen, aber da stimmt etwas nicht. Ich bekomme folgende Fehler:

Use of uninitialized value $organism in concatenation (.) or string at test.pl line 39, <IN> chunk 1 
Use of uninitialized value $shortest in concatenation (.) or string at test.pl line 39, <IN> chunk 1. 
Use of uninitialized value $longest in concatenation (.) or string at test.pl line 40, <IN> chunk 1. 

Welchen Teil soll ich ändern? Danke

Answer 1

Sie sagen, dass Sie zwei Stücke von Daten wollen - der Beitritt und der Organismus - für die längste und kürzeste Sequenz. Dies bedeutet, dass Ihre Hash-Werte zwei Elemente speichern müssen. Wenn Sie '//' als Datensatztrennzeichen verwenden, wird das '//' immer noch am Ende jedes Datensatzes angezeigt. Wenn Sie Whitespace und Ziffern aus Ihrer Sequenz herausfiltern, bleibt am Ende immer noch '//' übrig. Als ich den Code durch den Debugger laufen ließ, stellte ich fest, dass die Längen aus diesem Grund alle um 2 waren.

Ein paar andere Dinge:

1. Wenn regexs, verwenden Sie ‚erweiterten Modus‘, /x, so dass Sie Leerzeichen, um readabillity einschließen
2. Sie ein erfolgreiches Spiel vermuten, wenn Sie $definition raben - besser um deine Regexs zu testen und zuzuweisen, bei Übereinstimmung,
3. Anstatt die Länge im Hash zu speichern (und die Sequenz selbst zu verlieren), könnten Sie auch die Sequenz speichern und die Längen später berechnen;
4. umbenannt ich die Variable $line-$chunk, da sie mehrere Linien
5. All das Zeug enthält mit der Berechnung der kürzesten und längsten und Drucken der resuts braucht, um aus der Schleife zu tun. An seiner Stelle muss lediglich ein Eintrag in den Hash vorgenommen werden. Wie oben beschrieben, müssen die Hash-Werte ein Array mit zwei Werten sein - dem Beitritt und dem Organismus.
6. Sie entfernen Ziffern aus der Sequenz in einem Befehl und dann Whitespace aus der Sequenz in einem anderen - könnte sie auch beide zusammenbringen. Wenn wir gerade dabei sind, können wir auch die '/' am Ende des Datensatzes entfernen.

Angesichts der Mods oben, bekomme ich;

use v5.14; 
use warnings; 

print "Enter file path: "; 
chomp(my $filename = <>); 
open(IN, $filename) or die "\n error opening file \n;/\n"; 

$/ = "//" ; 

my %organisms ; 
while (my $chunk = <IN>) { 
    next if $chunk =~ /^\s*\n\s*$/ ; 
    my ($definition , $sequence) = split "ORIGIN", $chunk ; 

    my $organism ; 
    $definition =~ m/ ORGANISM [\s\t]+ (.+) [\n\s\t]+ /x 
     ? $organism = $1 
     : die "Couldnt find ORGANISM line" ; 

    my $accession ; 
    $definition =~ m/ ACCESSION [\s\t]+ (\D\D _ \d{6} (\d*)) [\n\s\t]+ /x 
     ? $accession = $1 
     : die "Cant find ACCESSION line" ; 

    $sequence =~ s/[\d\n\s\t\/]//g; 
    $organisms{ $sequence } = [ $accession , $organism ] ; 
} 


my @sorted_keys = sort { length $a <=> length $b } keys %organisms ; 
my $shortest = $sorted_keys[0]; 
my $longest = $sorted_keys[-1]; 

say "this is the shortest: ", $organisms{$shortest}->[0], 
         ", ", $organisms{$shortest}->[1], 
        " size: ", length $shortest, "\n", 
       " sequence: ", $shortest ; 

say "this is the longest: ", $organisms{$longest}->[0], 
         ", ", $organisms{$longest}->[1], 
        " size: ", length $longest, "\n", 
       " sequence: ", $longest ; 

exit; 

wenn es auf Ihren Daten lief, produziert es;

$ ./sequence.pl 
Enter file path: data.txt 
this is the shortest: NM_001323410, Homo sapiens size: 113 
sequence: actacttccggcttccccgccccgccccgtccccgggcgtctccattttggtctcaggtgtggactcggcaagaaccagcgcaagagggaagcagagttatagctaccccggc 
this is the longest: NM_182854, Homo sapiens size: 120 
sequence: gggcgatcagaagcaggtcacacagcctgtttcctgttttcaaacggggaacttagaaagtggcagcccctcggcttgtcgccggagctgagaaccaagagctcgaaggggccatatgac 

UPDATE Das Problem mit dem obigen Code ist, dass, wenn die gleiche Sequenz in zwei Stücken erscheint, dann in der Hash-und verloren überschrieben Daten werden werden. Unten ist eine aktualisierte Version, die Daten in einem Array von Arrays speichert, die das Problem advoid sind.Es erzeugt genau die gleiche Ausgabe:

Answer 2

Wir müssen Einträge extremer Länge finden, während wir in der Lage sind, den Datensatz zu identifizieren, zu dem sie gehören. Lesen von Aufzeichnungen von // ist wieder eine nette Idee. Dann ist jeder Datensatz jedoch eine Zeichenkette, und es ist schwieriger, die Sequenz direkt daraus zu ziehen, als sie zuerst in Zeilen zu zerlegen. So können wir auch Zeile für Zeile gehen, da es klare Markierungen für alles gibt, was benötigt wird.

Eine Auswahl der Datenstruktur ist wichtig und hängt vom Zweck ab. Hier habe ich Daten so organisieren, dass es einfach ist, mit zu arbeiten, in einen Hash mit Elementen

%block = ('accession' => { 'type' => type, 'sequence' => sequence }, ...) 

Die Suche durchführen, sobald die Daten gelesen werden, in würde stark durch die Organisation dies durch ‚Sequenz‘ unterstützt werden (statt durch "Beitritt"), aber das würde die Zusammenarbeit sehr erschweren. Ich nehme an, dass dies am Ende für mehr verwendet werden kann, und dass ein kleiner Geschwindigkeitsverlust nicht von Bedeutung ist. Wenn das einzige Ziel hier wäre, die spezifische Frage mit optimaler Leistung zu beantworten, wären andere Ansätze geeigneter. Kommentare folgen dem Code.

use warnings; 
use strict; 
use feature qw(say); 

my $file = 'data_seqs.txt'; 
open my $fh, '<', $file or die "Can't open $file -- $!"; 

# Hash, helper variables, flag (inside a sequence?), sequence-end marker 
my (%block, $accession, $sequence); 
my $is_seq = 0; 
my $end_marker = qr(\s*//); # marks end of sequence: // 

while (my $line = <$fh>) 
{ 
    chomp($line); 
    next if $line =~ /^\s*$/;  # skip empty lines 

    if ($line =~ /$end_marker/) { # done with the sequence 
     $is_seq = 0; 
     $sequence = ''; 
     next; 
    } 

    if ($line =~ /^\s*ACCESSION\s*(\w+)/) { 
     $accession = $1; 
    } 
    elsif ($line =~ /^\s*ORGANISM\s*(.+)/) { 
     $block{$accession}{'type'} = $1; 
    } 
    elsif ($line =~ /^\s*ORIGIN/) { # start sequence on next line 
     $is_seq = 1; 
    } 
    elsif ($is_seq) {    # read (and add to) sequence 
     if ($line =~ /^\s*\d+\s*(.*)/) { 
      $block{$accession}{'sequence'} .= $1; 
     } 
     else { warn "Not sequence? Line: $line " } 
    } 
} 

# Identify keys for max and min lenght. Initialize with any keys 
my ($max, $min) = keys %block; 

foreach my $acc (keys %block) 
{ 
    my $current_len = length($block{$acc}{'sequence'}); 
    if ($current_len > length($block{$max}{'sequence'})) { 
     $max = $acc; 
    } 
    if ($current_len < length($block{$min}{'sequence'})) { 
     $min = $acc; 
    } 
} 

say "Maximum length sequence: ACCESSION: $max, ORGANISM: " . $block{$max}{'type'}; 
say "Minimum length sequence: ACCESSION: $min, ORGANISM: " . $block{$min}{'type'}; 

use Data::Dumper; 
print Dumper(\%block); 

Diese Drucke (Muldenkipper der Ausdruck weggelassen)

 
Maximum length sequence: ACCESSION: NM_182854, ORGANISM Homo sapiens 
Minimum length sequence: ACCESSION: NM_001323410, ORGANISM Homo sapiens 

einen Kommentar über die Effizienz der Suche

Ein weit verbreiteter Ansatz ein Reverse-Lookup-Hash zuerst bauen würde, dann eine Bibliothek verwenden, sagen aus List::Utils, um Max und Min zu finden, dann schaue nach, wo sie hingehören. Dafür müssen wir den Lookup-Hash erstellen und wir würden die Bibliothek zweimal verwenden, während das Durchsuchen von Hand, wie oben beschrieben, dazu führt, dass man die Struktur durchläuft und auch einfacher ist. Eine andere Option wäre, Hash-Top-Level-Schlüssel als Sequenzen zu haben und dann direkt Max und Min zu finden. Ein solches Hash wäre jedoch wesentlich schwieriger zu handhaben. Ein weiterer Ansatz wäre, Daten in einer Struktur zu organisieren, die ein effizienteres Abrufen dieser spezifischen Informationen ermöglicht, die wahrscheinlich auf Arrays basieren.

Allerdings scheint der Effizienzgewinn den großen Komfortverlust nicht zu rechtfertigen. Wenn sich die Geschwindigkeit als Problem herausstellt, sollte dies berücksichtigt werden.

Wenn Sie mit mehreren Dateien arbeiten müssen, ändern Sie einfach die Schleife zu while (<>) und senden Sie sie in der Befehlszeile. Alle Zeilen von allen werden dann Zeile für Zeile gelesen und der Code bleibt gleich.

Es kann sein, dass ich einige Begriffe missverstanden habe. Ich entferne keine leeren Stellen aus der "Sequenz" und verwende Wörter in der ersten Zeile nur für "Typ", nur um ein paar Kandidaten zu nennen. Diese sind einfach zu justieren, lass es mich wissen.