Perl: grep in gx() aufrufen

Computer: Software & Programmierung Perl & CGI
#1

Hallo

Ich habe gerade erst mit Perl angefangen. Daher kenne ich mich leider noch nicht wirklich aus.

Ich habe folgenden (noch unfertigen) Code geschrieben. Jedoch funktioniert er noch nicht wie gewünscht.

$data_file=„demo“;
open(DAT, $data_file) or die(„Could not open file!“;
while(){
#Nimmt nur Zeilen auf die mit „>“ anfangen
next if !($_=~/^>/);
push @seq_names,$_;
}
close(DAT);

for($k=0;$kU93944.1_238-348 in U93944.1/238-348 um
$tmp = substr($seq_names[$k],1);
$tmp =~ tr !_!/!;
chomp($tmp);
@testen = qx(grep „’$tmp\|SS_cons’“ FILE | grep -A1 „$tmp“);
}

Leider funktioniert der grep Aufruf nicht. Wenn ich ihn per echo ausgeben lasse, dann steht er eigentlich richtig da. Denn er soll lauten:

grep ‚U93944.1/238-348|SS_cons‘ FILE | grep -A1 U93944.1/238-348

Zur Erklärung.
Das FILE hat ungefähr 100 000 Zeilen und hat folgende Form:

irgendwelche zeilen
U93944.1/238-348
a9a344.3/28-38
E251.8/124-238
usw.
SS_cons

irgendwelche zeilen

U93944.1/238-348
a9a344.3/28-38
E251.8/124-238
usw.
SS_cons

usw.

Irgendwann kommt ein Block von Daten wo U93944.1/238-348 nicht enthalten ist. In dem Fall interessiere ich mich dann für SS_cons auch NICHT!

Deswegen filtere ich mir mit grep alle Zeilen mit U93944.1/238-348 und SS_cons aus.

Dies würde dann so aussehen:

U93944.1/238-348
SS_cons
U93944.1/238-348
SS_cons
SS_cons
SS_cons
SS_cons

Die letzten SS_cons interessieren mich nicht. Deswegen filtere ich noch 1 Zeile nach dem letzten U93944.1/238-348.

Somit sieht es dann so aus:

U93944.1/238-348
SS_cons
U93944.1/238-348
SS_cons

Jedoch weiß man nicht welchen String (also U93944.1/238-348) man sucht. Es können auch mehrere sein, die in der $data_file stehen.

Das Endergebnis soll dann so aussehen:

U93944.1/238-348
SS_cons
U93944.1/238-348
SS_cons
A93723.1/0-28
SS_cons
A93723.1/0-28
SS_cons
A93723.1/0-28
SS_cons
S203.3/29-189
SS_cons

Das Ergebnis muss ich dann noch weiter verarbeiten.
Ich hoffe, ich habe nicht all zu verwirrend geschrieben :smiley:.

Ich glaube außerdem das ich in der letzten for-schleife irgendwie noch ein push benutzen muss damit das ganze Ergebnis in einem Array steht.

Für Hilfe bedanke ich mich bereits im voraus :wink:

Mfg
Dennis

#2

@testen = qx(grep „’$tmp\|SS_cons’“ FILE | grep -A1
„$tmp“);

sind da ein paar anführungsstriche zu viel?

Leider funktioniert der grep Aufruf nicht. Wenn ich ihn per
echo ausgeben lasse, dann steht er eigentlich richtig da. Denn
er soll lauten:

grep ‚U93944.1/238-348|SS_cons‘ FILE | grep -A1
U93944.1/238-348

da überlistest du dich vielleicht selbst. wenn du einen string an qx()
übergibst und ihn vorher testen willst, testest du ihn nicht mit echo,
sondern vorher in perl. das echo gibt z.b. anführungszeichen nicht (alle)
mit aus.
und wenn du das in perl machst, siehst du, dass du auf „’…’“ matchst, also
doppelte plus einfache anführungszeichen.

#3

@testen = qx(grep „’$tmp\|SS_cons’“ FILE | grep -A1
„$tmp“);

sind da ein paar anführungsstriche zu viel?

Naja, es müssen soweit ich das festgestellt habe so viele sein.
Damit ich eine Variable benutzen kann muss ich sie in " " setzen, also muss ich „$tmp“ schreiben damit da der Wert eingesetzt wird.
Damit ich einen grep Aufruf mit ODER machen kann muss ich dies in der Konsole auch mit ’ ’ machen. Sprich grep ‚wort1|wort2‘ FILE.
Ansonsten geht es ja schon in einer normalen Konsole nicht.
Ich kann mich natürlich auch irren.
Aber wenn ich

@testen = qx(grep „$tmp\|SS_cons“ FILE | grep -A1 „$tmp“);

mache gehts nicht.
Wenn ich

@testen = qx(grep ‚$tmp\|SS_cons‘ FILE | grep -A1 „$tmp“);

mache gehts auch nicht.

Und wenn ich ohne anführungsstriche mache sagt mir die Konsole:
„SS_cons: not found“ und „hängt“ sich sozusagen auf indem er nicht weiter mach und ich muss mit STR+c das Programm abbrechen.

Mfg
Dennis

#4

Ok, ich habs umgestellt.

Wenn man

grep -e Wort1 -e Wort2 … -e WortN FILE

schreibt, dann ist es das gleiche wie

grep ‚Wort1|Wort2|…|WortN‘ File

Somit umgehe ich das Problem.

Meine Lösung sieht nun so aus:

@testen = qx(grep -e „$tmp“ -e SS_cons FILE | grep -A1 „$tmp“);

Gruß
Dennis

#5

Hallo Dennis

@testen = qx(grep -e „$tmp“ -e SS_cons FILE | grep -A1

Dein Problem ist ein *typischer* Anwendungsfall
für Perl. Die Einbeziehung von externem grep halte
ich für kontraproduktiv, das kann man viel besser
direkt in Perl machen, sprich: das ganze ist am
Ende ein kleines schnelles Perlprogramm mit 4-5
Zeilen.

Kannst Du mal bitte ein konkretes Beispiel
für Deine beiden Dateien irgendwo temporär
als zip hochladen (zum testen). Ich vermute,
man kann bereits die 100Klines-Datei in Records,
getrennt an ‚SS_cons‘ einlesen.

 [pseudo]

 ...
 $tmp = $seq\_names[$k]
 $/='SS\_CONS' # input record separator
 while( $record = ) {
 print grep $record =~ /$\_/sg for @seq\_names
 }
 ...

 [/pseudo]

Das würde wahrscheinlich für 100Klines
in Sekundenbruchteilen durchlaufen.

Grüße

CMБ

#6

Natürlich wäre ich über eine effizientere Methode sehr dankbar. Denn ich muss ziemlich viele solcher Daten erstellen…

Auf www.wehrle.it/test.tar.gz habe ich die 2 Testdaten hochgeladen.
Man muss in der demo.ref.fa die Namen rausfiltern und dann in der test.data die ganze Zeile rausholen wo der Name steht und die darunterliegende SS_cons.

Ich kenne mich mit Linux nicht wirklich aus und das war der einzige Ansatz der mir dann halt eingefallen ist. Ich hatte schon öfters gelesen, das man das Perl grep benutzen soll. Aber dort wurde immer gesagt, man soll den Inhalt von der Datei in ein Array stecken und dann da mit grep durchgehen. Jedoch kam mir das nicht sehr effizient vor wenn ich ein Array mit 100 000 Datensätze habe…

In der Datei „ergebnis“ sieht man wie es dann hinterher aussieht (nach meiner langsamen Methode). Jedoch ist das noch nicht das Endergebnis. Diese Datei muss ich dann noch weiter verarbeiten und bestimmte zeichen rauslöschen falls sie an der gleichen Position eines anderen Strings vorkommen. Dafür habe ich mir aber schon ne for-Schleife geschrieben (aber auch noch nicht ganz fertig).

Vielen Dank für die Hilfe!

Gruß
Dennis

[Bei dieser Antwort wurde das Vollzitat nachträglich automatisiert entfernt]

#7

Naja, es müssen soweit ich das festgestellt habe so viele
sein.
Damit ich eine Variable benutzen kann muss ich sie in " "
setzen, also muss ich „$tmp“ schreiben damit da der Wert
eingesetzt wird.

nein, da verwechselst du was. qx() interpoliert automatisch, d.h.
die shell kriegt von $tmp gar nix mit, sie kriegt schon den inhalt
dieser perl-variable.

#8

Hallo Dennis,

Natürlich wäre ich über eine effizientere Methode sehr
dankbar. Denn ich muss ziemlich viele solcher Daten
erstellen…

Es gibt zwei prinzipielle Möglichkeiten, je nach
dem, wie *viele* Sequenzschlüssel in der ‚demo.ref.fa‘-
Datei stehen.

  • (1) Einen Suchausdruck aus *allen* keys basteln, also
    (AF158725.1_345-497|AF444327.1_152-305|…)
    und über das komplette Datenfile laufen lassen,

  • (2) Das Datenfile in „Paragraphen/Records“ (Abschnitte mit
    einer Leerzeile dazwischen) einlesen und für jeden Record
    die Suche der in einem Feld gespeicherten Sequenzschlüssel
    durchführen.

Die Methode „zwei“ ist straightforward (s.u.), die Methode
„eins“ ist etwas abstrakter.

Man muss in der demo.ref.fa die Namen rausfiltern und dann in
der test.data die ganze Zeile rausholen wo der Name steht und
die darunterliegende SS_cons.

OK, ich habe festgestellt, daß man in den Schlüsselworten die
‚_‘ durch ‚/‘ ersetzen muß. Meine Variante für (2) sähe so aus:

use strict;
use warnings;

my ($fn\_keys, $fn\_data) = qw(demo.ref.fa test.data); # Dateinamen
my %erg; # suchergebnis merken
# (1) Steuerdatei oeffnen und einlesen
open my $fhk,'(\S+)/g } ; # '\_' zu '/' wandeln
close $fhk;

local $/=''; # In 'paragraph mode' (Block\n\n) lesen
open my $fhd, ') { # Einen 'Block' lesen
 for my $s (@seqn) { # Im 'Block' auf alle keywords testen
 push @{ $erg{$s} }, "$1\n$2\n" # und gefundene Daten sammeln in $erg{$\_}
 if $r =~ /^($s[^\n]+).+?(\#=GC SS\_cons[^\n]+)/ms # \*WENN\* o.k.
 } 
}
close $fhd;

print map "$\_\t" . @{$erg{$\_}} . " hits!\n", sort keys %erg; # Hits zählen oder
print '='x79 . "\n";
print map join('', @{$erg{$\_}}), sort keys %erg; # Daten weiterverarbeiten

Dabei liegen die Resultate nach der Suche im Hash
%erg vor, dessen Struktur wäre also dann:

%erg = ('AF158725.1/345-497' =\> 
 [ 'AF158725.1/345-497 AUCU....\n #=GC SS\_cons \> ...',
 ...
 ],
 'Y07976.1/124-271' =\> 
 [
 '...',
 '...'
 ],
 ...
 );

Damit lassen sich die Daten gleich weiterverarbeiten.
Nun wäre noch interessant zu wissen, was genau mit den
Daten weiter passieren soll.

Grüße

CMБ

#9

Anmerkung:

In den Keys:

AF158725.1/345-497, AF444327.1_152/305

kommen Zeichen vor, die in regulären Ausdrücken
eine besondere bedeutung haben. Hier muß man also
aufpassen! In Perl benutzt man zur Maskierung
den Quotemeta-Operator (quotemeta($str)) oder
in Strings bzw. regulären Ausdrücken die
\Q$str\E - escape-Sequenz. Im konkreten Fall
hätte das aber zufälligerweise keine Auswirkung :wink:

if $r =~ /^($s[^\n]+).+?(#=GC SS_cons[^\n]+)/ms #

wird z.B:

 if $r =~ /^(\Q$s\E[^\n]+).+?(\#=GC SS\_cons[^\n]+)/ms #

Ausserdem kann man das Hilfsarray @seq ganz weglassen,
also etwa so hier:

use strict;
use warnings;

my ($fn\_keys, $fn\_data) = qw(demo.ref.fa test.data); # Dateinamen

open my $fhk,'my %erg = map{ y{\_}{/}; /^\>(\S+)/ ? ($1, []) : () } ; 
close $fhk;

local $/ = ''; # In 'paragraph mode' (Block\n\n) lesen
open my $fhd, ') { # Einen 'Block' lesen
**while(my ($key,$val) = each %erg)** { # Im 'Block' auf alle keywords testen
 if($record =~ /^(\Q$key\E[^\n]+).+?(\#=GC SS\_cons[^\n]+)/ms) {
 push @$val, "$1\n$2\n" # und gefundene Daten sammeln in $erg{$\_}
 }
 } 
}
close $fhd;

print map "$\_\t" . @{$erg{$\_}} . " hits!\n", sort keys %erg; # Hits zählen oder
print '='x79 . "\n";
print map join('', @{$erg{$\_}}), sort keys %erg; # Daten weiterverarbeiten

Grüße

CMБ

#10

Damit lassen sich die Daten gleich weiterverarbeiten.
Nun wäre noch interessant zu wissen, was genau mit den
Daten weiter passieren soll.

Grüße

CMБ

Hallo

Erstmal danke für die Mühe!!

Wie gesagt, ich fange grade erst querbeet an mit Perl zu arbeiten. Daher verstehe ich den ganzen Code noch nicht. Ich werd mich aber reinarbeiten.

Nun zu deiner Frage.
Also ich habe viele Ordner und Unterordner in denen solche demo.ref.fa (vielemehr DATEINAME.ref.fa) stehen. In den jeweiligen Ordner muss dann das Endergebnis dann auch abgespeichert werden. In jeder Datei stehen dann eigentlich 2-15 solcher Namen die rausgefiltert werden müssen.

Auf http://www.wehrle.it/testdaten habe ich ein Beispiel, wie das noch zu bearbeitende Ergebnis aussehen soll.
Hierbei gilt zu beachten das die Daten so wie sie aus „test.data“ rausgenommen werden noch dahingehend weiter bearbeitet werden müssen, das wenn bei beiden Sequenzen an der gleichen Position ein „.“ ist (egal was an der Position bei SS_cons steht), das diese Position gelöscht werden muss, sowohl in der Sequenz selber als auch in der SS_cons dann.

Also Beispielsweise:
>SEQUENZ1
ACG.GUCG…AAUCCG.AUUCCC
…>>
^ ^
>SEQUENZ2
UGC.UCCGG.GGUUCAAGC.GU.
…>>
^ ^
Hier müsste jetzt (von der 1. Sequenz ausgehend) der 1. und der 3. „.“ (siehe ^) gelöscht werden (sowie auch die Position bei SS_cons).

Endergebnis:
>SEQUENZ1
ACGGUCG.AAUCCG.AUUCCC
…>>

>SEQUENZ2
UGCUCCGGGGUUCAAGC.GU.
…>>

Dieses Endergebnis muss dann in DATEINAME.ref.fa.cons gespeichert werden.
Die Ordnerstruktur sieht so aus, dass das Perlscript in einem Ordner aufgerufen wird in denen viele verschiedene Unterordner sind. In diesen verschiedenen Unterordnern liegen dann u.a. die DATEINAME.ref.fa (natürlich mehrere).

Vielen Dank für die Hilfe!!

Gruß
Dennis

#11

Hallo Dennis

Also ich habe viele Ordner und Unterordner in denen solche
demo.ref.fa (vielemehr DATEINAME.ref.fa) stehen. In den
jeweiligen Ordner muss dann das Endergebnis dann auch
abgespeichert werden. In jeder Datei stehen dann eigentlich
2-15 solcher Namen die rausgefiltert werden müssen.

OK, heißt das, Du hast in einem Verzeichnis eine
100K-lines ‚test.data‘ und und einer Hierarchie von
Unterverzeichnissen jeweils ein *.ref.fa pro Verzeichnis,
welche sich dann auf das eine ‚test.data‘ im
Hauptverzeichnis beziehen? Wieviele Unterverzeichnisse
wären das denn so?

Dann würde man, um eine Laufzeit im unteren Sekundenbereich
beizubehalten, wahrscheinlich die ‚test.data‘-Datei Datei
einmal einlesen und ‚record‘-Weise im Hauptspeicher halten.

Die Ordnerstruktur sieht so aus, dass das Perlscript in einem
Ordner aufgerufen wird in denen viele verschiedene Unterordner
sind. In diesen verschiedenen Unterordnern liegen dann u.a.
die DATEINAME.ref.fa (natürlich mehrere).

Wie gesagt, es wäre interessant zu wissen ob sich die
.fa-Dateien immer aus dieselbe Sequenzdatei beziehen.
Und wie viele .fa-Dateien pro Unterverzeichnis vorkommen.

Die rekusive Suche durch die Verzeichnisse ist recht
einfach, das würde so aussehen:

use strict;
use warnings;
use File::Find;

my ($fn\_data, $startdir) = @ARGV; 
$startdir = '.' unless defined $startdir;

my @filelist = ();
scancontrol($startdir) or die "keine .ref.fa files";

#
# hier kommt das eigentliche Programm hin, das
# wird eine Schleife ueber @filelist
#

sub scancontrol {
 my $from\_where = shift;
 print "searching directories below '$from\_where':\n";
 find( { wanted =\> \&wanted, # sub for processing
 follow =\> 1, # follow symlinks
 no\_chdir =\> 1, # stay in original directory
 follow\_skip =\> 2 # silently ignore duplicated files
 },
 $from\_where
 );
 return scalar @filelist 
}

sub wanted {
 if( -f && /\.ref\.fa$/ ) {
 push @filelist, $\_ ;
 print "$\_\n"
 }
}

Obiges Programmfragment würde Dir schon
alle vollständigen Pfadnamen der ‚.ref.fa‘ -
Dateien liefern. Jetzt müßte man einfügen,
auf welchen .data-Files diese arbeiten sollen.

Siehe: http://search.cpan.org/~nwclark/perl-5.8.5/lib/File/…

Grüße

CMБ

#12

Hallo

Nochmals vielen Dank für die Mühe!!

Also es gibt nur EINE .data File (heißt in wirklichkeit Rfamm.seed).
Die Ordnerstruktur sieht so aus:

Hauptordner:
Rfam.seed

  1. Unterordner: -> vielleicht ca. 80 .ref.fa Dateien (aber auch noch andere Dateien)
  2. Unterordner: -> ebenfalls viele .ref.fa Dateien und andere
    usw.

Also es geht nicht weiter runter als ein Ordner.
Im Hauptverzeichnis liegt also die Rfam.seed und dann gibts da halt viele andere Ordner in denen dann viele .ref.fa Dateien liegen.
Wie gesagt für jede einzelne .ref.fa Datei muss eine .ref.fa.cons Datei erstellt werden.

Aber generell könnte man natürlich den Dateiinhalt von der Rfam.seed in einem Record lassen und darin dann navigieren. Denn es soll ja eigentlich für alle .ref.fa die .cons Datei erstellt werden. Es wäre vielleicht nur von Vorteil, wenn man vorher ne sicherheitsabfrage machen würde das falls für die Datei schon eine .cons Datei existiert, dass man dann die nicht noch mal erstellt. Also das ist vor allem dann halt wichtig falls neue .ref.fa Dateien hinzukommen und nicht noch mal alles neu berechnet werden soll.

Vielen Dank für deine Hilfe!!

Gruß
Dennis

[Bei dieser Antwort wurde das Vollzitat nachträglich automatisiert entfernt]

#13

Hallo Dennis

gestern hatte ich schwer zu tun, da kam ich zu nix :wink:

Also es gibt nur EINE .data File (heißt in wirklichkeit
Rfamm.seed).

OK

Die Ordnerstruktur sieht so aus:
Hauptordner:
Rfam.seed

  1. Unterordner: -> vielleicht ca. 80 .ref.fa Dateien (aber
    auch noch andere Dateien)
  2. Unterordner: -> ebenfalls viele .ref.fa Dateien und
    andere
    usw.

Also es geht nicht weiter runter als ein Ordner.
Im Hauptverzeichnis liegt also die Rfam.seed und dann gibts da
halt viele andere Ordner in denen dann viele .ref.fa Dateien
liegen.
Wie gesagt für jede einzelne .ref.fa Datei muss eine
.ref.fa.cons Datei erstellt werden.

Soweit klar. Nur eine Frage: Du liest ja die Sequenzen
aus den .ref.fa-Dateien ein (nach &gt:wink:. Kann es auch vor-
kommen, daß in unterschiedlichen .ref.fa-Dateien die
gleichen Sequenzen vorkommen? Daß z.B. die
Sequenz X53538.1_2090-2241 sowohl in
/unterordner1/x1.ref.fa
als auch in
/unterordner5/x22.ref.fa
vorkommt? Oder geht das nicht. *Wenn*
nämlich jede Sequenz nur einmal vorkommen
kann, wäre das Programm kürzer (man würde die
Sequenz als Hash-Schlüssel für den Pfad und alles
andere nehmen) :wink:

Es wäre vielleicht nur von Vorteil, wenn man
vorher ne sicherheitsabfrage machen würde das falls für die
Datei schon eine .cons Datei existiert, dass man dann die
nicht noch mal erstellt. Also das ist vor allem dann halt
wichtig falls neue .ref.fa Dateien hinzukommen und nicht noch
mal alles neu berechnet werden soll.

Das ist kein Problem:

if(-f 'x22.ref.fa.cons') {
 warn "! x22.ref.fa.cons schon da !\n"
}
else {
 open ..., 'x22.ref.fa.cons' ...
 usw.
}

Grüße

CMБ

#14

Hallo,

hierzu noch eine Anmerkung/Frage:

Hierbei gilt zu beachten das die Daten so wie sie aus
„test.data“ rausgenommen werden noch dahingehend weiter
bearbeitet werden müssen, das wenn bei beiden Sequenzen an der
gleichen Position ein „.“ ist (egal was an der Position bei
SS_cons steht), das diese Position gelöscht werden muss,
sowohl in der Sequenz selber als auch in der SS_cons dann.

Also Beispielsweise:
>SEQUENZ1
ACG.GUCG…AAUCCG.AUUCCC
…>>
^ ^
>SEQUENZ2
UGC.UCCGG.GGUUCAAGC.GU.
…>>
^ ^
Hier müsste jetzt (von der 1. Sequenz ausgehend) der 1. und
der 3. „.“ (siehe ^) gelöscht werden (sowie auch die Position
bei SS_cons).

Endergebnis:
>SEQUENZ1
ACGGUCG.AAUCCG.AUUCCC
…>>

>SEQUENZ2
UGCUCCGGGGUUCAAGC.GU.
…>>

Was bedeutet hier: „wenn bei beiden Sequenzen“?
In jeder .ref.fa-Datei sind ja jeweils 5-15 Sequenzen
(wie Du schreibst) und von jeder Sequenz gibt es
ja im ‚rfam.seed‘ auch wieder etliche Einträge.

Was also meinst Du mit „bei beiden Sequenzen“?

Grüße

CMБ

#15

Soweit klar. Nur eine Frage: Du liest ja die Sequenzen

Hallo

Und ein weiteres mal vielen Dank!!

Leider kann man das nicht so sagen. Ich kann dir leider versprechen das innerhalb eines Unterordners die Sequenz öfters vorkommt!

Also Seq1 wir in der einen Datei mit Seq2 verglichen und in der nächsten Datei wird Seq1 mit Seq3 verglichen usw.

Gruß
Dennis

#16

Hallo,

Leider kann man das nicht so sagen. Ich kann dir leider
versprechen das innerhalb eines Unterordners die Sequenz
öfters vorkommt!

Also Seq1 wir in der einen Datei mit Seq2 verglichen und in
der nächsten Datei wird Seq1 mit Seq3 verglichen usw.

Aha. OK, aber in jeder ref.fa-Datei stehen
ja immer *viele* Sequenzbezeichner (oder
immern nur zwei?). Was wird womit verglichen,
wenn in einer .ref.fa-Datei 10 Sequenzbezeichner
stehen?

Grüße

CMБ

#17

Hallo

Sorry, ich hatte gestern viel zu tun, deswegen erst jetzt meine Antwort.

Also in einer .ref.fa Datei stehen immer mind. 2 Sequenzen. Es können aber natürlich auch mehr sein! Also da kann ich leider keine feste Anzahl sagen, bzw. es sollte halt für beliebige Anzahl funktionieren. Eigentlich dürften aber auch nicht mehr als 15 solcher Sequenzen drin sein in einer Datei.

Ich habe gemerkt, das du bei deinem Re8 zwei mal geantwortet hattest und ich dabei eine Antwort überlesen hatte. Da hast du ja wegen den Sequenzen gefragt.

Also es ist eigentlich ganz einfach. Wenn beispielsweise in einer .ref.fa Datei 5 Sequenzen drin sind, dann holt man sich ja die Namen aus der .ref.fa Datei. Anhand der Namen werden dann ja die Sequenzen und die SS_cons „Strings“ aus der Rfam.seed rausgeholt. Dann muss man jetzt einfach hergehen und in der ersten Sequenz schauen wo überall ein „.“ ist (also nicht in der SS_cons, sondern in der Sequenz). Jetzt muss man einfach bei den restlichen 4 Sequenzen schauen ob an der gleichen Position auch ein „.“ ist. Wenn ja, dann muss in allen 5 Sequenzen und in allen SS_cons genau die Position gelöscht werden.

Beispiel für 3 Sequenzen:
>Seq1
A A A A . A A A . A
>Seq2
A A A A . A A A A A
>Seq3
A A A A . A A A . A
1 2 3 4 5 6 7 8 9 10 (Position)

Hier müsste die 5. Position in jeder Sequenz und die 5. Position in jeder SS_cons einfach rausgelöscht werden, da alle Sequenzen an der 5. Position ein „.“ haben. Die 9. Position muss allerdings NICHT gelöscht werden, da Seq2 an Position 9 kein „.“ besitzt. Somit passiert an der Stelle nichts und es muss in der 1. Sequenz weiter nach „.“ geschaut werden.

Ergebnis:
Seq1
A A A A A A A . A
Seq2
A A A A A A A A A
Seq3
A A A A A A A . A
1 2 3 4 5 6 7 8 9 (Position)

Also am einfachsten ist es jeder Position der ersten Sequenz durchzugehen und dann halt schauen wo so ein „.“ ist und dann muss in allen anderen geschaut werden ob an dieser Position auch so ein „.“ ist.

Aha. OK, aber in jeder ref.fa-Datei stehen
ja immer *viele* Sequenzbezeichner (oder
immern nur zwei?). Was wird womit verglichen,
wenn in einer .ref.fa-Datei 10 Sequenzbezeichner
stehen?

Grüße

CMБ

Man muss also um die „.“ rauszulöschen nicht jede Sequenz mit jeder vergleichen, sondern nur die 1. mit allen anderen.

Gruß
Dennis

#18

Hallo

Vielen Dank für deine bisherige Hilfe!!!
Dank deiner Hilfe bin ich schon fast fertig.

Ich habe (vorerst) nur noch ein letztes Problem.

Das Ergebnis sieht noch so aus:

U133691/6624-6776     GACUGCGUCAGCCAGC
#=GC S_cons                …>. (-1.0)
U482281/7-166             CAAUCUCGAUGAAGGCCGCAGC
#=GC S_cons                …>… (-1.0)
U482281/7-166             A…A
#=GC S_cons                … (-1.0)

Aussehen soll es jedoch so:

>U133691/6624-6776
GACUGCGUCAGCCAGC
…>. (-1.0)

(Leere Zeile, also \n)

>U482281/7-166
CAAUCUCGAUGAAGGCCGCAGC
…>… (-1.0)
A…A
… (-1.0)

Also ich denke das es irgendwie mit dem split() Operator gehen könnte. Man
müsste erstmal die Sequenz von der SS_cons trennen. Danach müsste man noch mit
substr() irgendwie von hinten durchgehen bis zum ersten mal ein Leerzeichen
gefunden wird.

Aber wie es genau geht hab ich noch nicht rausgefunden. Werde aber morgen mal
weiter versuchen.

Mit freundlichen Grüßen
Dennis Wehrle

Bisheriger Code:

use strict;
use warnings;
use File::Find;

my ($fn_data, $startdir) = @ARGV;
$startdir = ‚.‘ unless defined $startdir;

my @filelist = ();
scancontrol($startdir) or die „keine .ref.fa files“;

while(@filelist){
job(pop(@filelist));
}

sub scancontrol {
my $from_where = shift;
print „durchsuche das Verzeichnis ‚$from_where‘:\n“;
find( { wanted => &wanted, # sub for processing
follow => 1, # follow symlinks
no_chdir => 1, # stay in original directory
follow_skip => 2 # silently ignore duplicated files
},
$from_where
);
return scalar @filelist
}

sub wanted {
if( -f && /.ref.fa$/ ) {
push @filelist, $_ ;
#print „$_\n“
}
}

sub job{
my $fn_keys = $_[0];

if(-f „$fn_keys.cons“){
warn „! $fn_keys.cons existiert schon !\n“
}
else{
open my $fhk,’(\S+)/ ? ($1, []) : () } ;
close $fhk;

local $/ = ‚‘; # In ‚paragraph mode‘ (Block\n\n) lesen
open my $fhd, ') { # Einen ‚Block‘ lesen
while(my ($key,$val) = each %erg) { # Im ‚Block‘ auf alle keywords testen
if($record =~ /^(\Q$key\E[^\n]+).+?(#=GC SS_cons[^\n]+)/ms) {
push @$val, „$1\n$2 (-1.0)\n“ # und gefundene Daten sammeln in
$erg{$_}
}
}
}
close $fhd;

Namen der Sequenzen

my @schluessel = keys(%erg);

Anzahl der Sequenzen

my $anzahl = keys(%erg);

Laenge der einzelnen Teilsequenzen

my @laenge;

Anzahl der Teilsequenzen

my $hits = @{$erg{$schluessel[0]}};

Boolean-Wert ob Spalte gelöscht werden soll

my $loesche=1;

Speichert die verschiedenen Längen der Teilsequenzen.

Die Daten liegen folgendermaßen vor: "Sequenzname Sequenz #=GC SS_cons

Konsensusstrukur (-1.0)"

Daher muss zunächst " (-1.0)" (=7 Zeichen) abgezogen und der Rest durch 2

geteilt werden. Somit ergibt sich die Länge für „Sequenzname Sequenz“
for(my $j=0;$j’, „$fn_keys.cons“ or die „cant’t open $fn_data: $!“;
print $fhc map join(’’, @{$erg{$_}}), sort keys %erg; # Daten
weiterverarbeiten
print „$fn_keys.cons wurde erstellt!\n“;
close $fhc;
}
}