Wednesday, October 03, 2007

Nützliche Variablen: $/

Im Forum tauchen immer wieder ähnliche Fragen auf und heute hatte ich wieder ne Mail im Postfach, bei denen es sich die Schreiber unnötig schwer machen. Da ich lange Zeit in dieser Ecke gearbeitet habe und die heutige Mail damit zu tun hatte möchte ich eine nützliche Spezialvariable in Perl mit einer FASTA-Datei erläutern.

FASTA-Dateien sind Dateien in einem bestimmten Format, in denen Sequenzdaten gespeichert sind - also mit Bioinformatik zu tun haben. Perl ist in der Bioinformatik eine weit verbreitete Programmiersprache und wenn man sich der Mächtigkeit von Perl bedient versteht man auch warum...

Diese FASTA-Dateien werden häufig eingelesen. Dabei möchte man meistens komplette Blöcke in ein Array schreiben. Also nicht Zeilenweise wie in den meisten Fällen in denen eine Datei eingelesen wird, sondern in benutzerdefinierten Blöcken.

Eine FASTA-Beispieldatei:
>gi|5524211|gb|AAD44166.1| cytochrome b [Elephas maximus maximus]
LCLYTHIGRNIYYGSYLYSETWNTGIMLLLITMATAFMGYVLPWGQMSFWGATVITNLFSAIPYIGTNLV
EWIWGGFSVDKATLNRFFAFHFILPFTMVALAGVHLTFLHETGSNNPLGLTSDSDKIPFHPYYTIKDFLG
LLILILLLLLLALLSPDMLGDPDNHMPADPLNTPLHIKPEWYFLFAYAILRSVPNKLGGVLALFLSIVIL
GLMPFLHTSKHRSMMLRPLSQALFWTLTMDLLTLTWIGSQPVEYPYTIIGQMASILYFSIILAFLPIAGX
IENY
>gi|1234567|gb|AAD44133.8| cytochrome a [Elephas maximus maximus]
LCLYTHIGRNIYYGSYLYSETWNTGIMLLLITMATAFMGYVLPWGQMSFWGATVITNLFSAIPYIGTNLV
EWIWGGFSVDKATLNRFFAFHFILPFTMVALAHTHLTFLHETGSNNPLGLTSDSDKIPFHPYYTIKDFLG
LLILILLLLLLALLSPDMLGDPDNHMPGLMITADSHIKPEWYFLFAYAILRSVPNKLGGVLALFLSIVIL
GLMPFLHTSKHRSMMLRPLSQALFWTLTMDLLTLTWIGSQPVEYPYTIIGQMASILYFSIILAFLPIAGX
IENY
In vielen Programmen sieht man dann so etwas:
#!/usr/bin/perl

my $fasta_file = '/pfad/zur/datei.fasta';
my @sequences;
my $sequence;
my $counter = 0;

open(FILEHANDLE, "<$fasta_file") or die $!; while( ){
if( $_ =~ /^>/ and $counter > 1 ){
push(@sequences, $sequence);
$sequence = $_;
}
$sequence .= $_;
$counter++;
}
push @sequences, $sequence if $sequence;
close(FILEHANDLE);

foreach(@sequences){
print "Sequenz: $_\n\n";
}

Das ist aber reichlich umständlich. Erstens gibt es schon Module wie etwas BioPerl oder Bio::FASTASequence oder wenn man es doch selbst machen will, sollte man die Sachen von Perl nutzen.

Es gibt eine Spezialvariable $/ (siehe auch perldoc perlvar), die das "Ende des einzulesenden Blocks" festlegt. Im Normalfall - wenn man also nichts verändert hat, ist die Variable auf "\n" gesetzt. Damit wird festgelegt, dass immer bis zu einem "\n" eingelesen wird.

Hier ist das zeilenweise einlesen aber eher ungeschickt. Wenn man sich die FASTA-Datei anschaut, sieht man, dass jeder Eintrag mit ">" an einem Zeilenanfang beginnt.

Also kann man $/ auf "\n>" setzen. Damit ist *ein* einzulesender Block durch ein Zeilenende-Zeichen und einem folgenden ">" begrenzt.

Somit kann man folgendes machen:
#!/usr/bin/perl

my $fasta_file = '/pfad/zur/datei.fasta';
my @sequences;

{
local $/ = "\n>";
open my $fh, '<', $fasta_file or die $!;
while(my $sequence = <$fh> ){
chomp $sequence;
$sequence = '>' . $sequence unless $sequence =~ /^>/;
push @sequences, $sequence;
}
close $fh;
}

foreach(@sequences){
print "Sequenz: $_\n\n";
}

Noch deutlicher wird es bei Dateien, bei denen die interessanten Blöcke durch eine Leerzeile begrenzt sind. Zum Beispiel:
Header1
Info1 zu Header1
Zeile2 der Info1

Header2
Info1 zu Header2
Zeile2 der Info1


Hier kann man $/ auf "\n\n" setzen und schon bekommt man bei jedem Durchlauf der while-Schleife beim Einlesen einen kompletten interessanten Block geliefert.

Aber man sollte aufpassen: Die Änderung von $/ sollte man auf einen möglichst kleinen Block beschränken. Deswegen habe ich oben auch um das Einlesen die geschweiften Klammern gepackt. Aber es gibt ein Problem wenn man *kein* local verwendet und innerhalb eines Programms erst eine FASTA-Datei dann eine andere Datei einlesen will...

chomp arbeitet übrigens auch mit $/, deswegen macht es in meinem Beispielcode auch das richtige und entfert das "\n>" bei jedem Eintrag...

No comments: