Discussion :

[Isabella83]

Bonjour à tous.
Quand je lance le script ci dessous sur un fichier test (4,0K), tout se passe bien. En revanche, lorsque je le lance sur mon fichier de 102M, il m'est indiqué après quelques minutes "processus arrêté", est ce que cela vient de mon script ?

Je cherche en fait à supprimer le premier nucleotide puis le second puis le troisième etc jusqu'a ce que ma sequence ait une longueur > à l'entier $int , par exemple, ici $int = 20 :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
 
>seq
CATCTGGTATCTGTTGGGTACATCC 
>seq_nb1
ATCTGGTATCTGTTGGGTACATCC 
>seq_nb2
TCTGGTATCTGTTGGGTACATCC 
>seq_nb3
CTGGTATCTGTTGGGTACATCC 
>seq_nb4
TGGTATCTGTTGGGTACATCC 
>seq_nb5
GGTATCTGTTGGGTACATCC

j'ai donc écrit ce script:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
#!/usr/bin/perl
use strict;
use warnings;
use Carp qw(confess);
use Getopt::Long;
use Bio::SeqIO;
 
my ($fasta_file,$out_file);
my ($int,$cut,$len);
my (%hash,%dico);
 
GetOptions("fasta=s" => \$fasta_file,"seed=i" => \$int ,"out=s" => \$out_file); 
 
 
#open (OUT, ">"."$input"."_"."$number".".fasta") or die "Cannot open outfile!";
 
 
my $in = Bio::SeqIO->new( -file => $fasta_file, '-format' => 'Fasta' );
open(my $out,'>',$out_file) or die "$out_file : $!\n\n";
 
while ( my $seq = $in->next_seq()){
	my $id = $seq->primary_id ;
	my $sequence = $seq->seq ;
	$len = $seq->length();
	my $count = 0;
	while($len >= $int){
		$cut = substr($sequence,$count);
		push(@{$dico{$id}->{'sequence'}}, $cut);
		$len--;
		$count++;
 
 
	}
 
}
 
foreach my $data (keys(%dico)){
	my $i=0;
	foreach my $seq (@{$dico{$data}->{'sequence'}}){
		print {$out} "$data"."_nb$i"."\n$seq\n";
		$i++;
	}
}

Je souhaiterai savoir s'il y a unproblème dans mon script ?
Merci d'avance,
Bonne journée à tous

[Philou67430]

J'avoue ne pas comprendre ton besoin réel (et je ne te demande pas, pour l'instant, de le rendre plus clair).
Cependant, j'ai l'impression que l'on peut remplacer ceci:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
	my $count = 0;
	while($len >= $int){
		$cut = substr($sequence,$count);
		push(@{$dico{$id}->{'sequence'}}, $cut);
		$len--;
		$count++;
	}

par

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

  push @{$dico{$id}->{sequence}}, split //, $sequence;

Bien entendu, ça ne répond pas à ta question (et ça ne la résoudra sans doute pas).
Je ne vois pas d'autre idée que l'erreur puisse venir de Bio::SeqIO (sur une séquence particulière peut-être).
Ton fichier fasta n'est pas vérolé ?

[Isabella83]

Merci pour ta réponse,
Non, mon fichier fasta est ok !
Ca ne peut pas venir d'un problème de mémoire ou autre du fait que mon fichier contienne des millions de sequences ?

[Philou67430]

Si bien sur, c'est tout à fait possible.
Je suis en train d'essayer d'installer Bio::IOSeq pour tester avec un énorme fichier de 235Mo (que j'ai généré avec l'uniligne suivant). Mais j'ai beaucoup de problème pour installer ce module, donc je n'aurai pas de résultat probant avant lundi.

Pour ceux qui voudraient t'aider, voici le script uniligne permettant de générer un fichier fasta aléatoire très gros.

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

$ perl -E 'foreach my $seq (0 .. 500000) { say ">seq_nb$seq"; say join "", map { [qw(A T C G)]->[rand(3)] } 0 .. rand(1000); }' > seq2.fasta

(attention aux quotes si vous utilisez Windows).

On peut faire varier la taille en jouant sur 500000 (nombre de séquences) et 1000 (taille max d'une séquence).

[Isabella83]

La taille maximale de mes séquences est de 29 nt

[stoyak]

Bonjour,

Il est probable que cela provienne de la façon dont tu parses ton fichier. En effet, le module Bio :: Seq met l'ensemble des séquences en mémoire. Vu la taille de ton fichier, il est bien possible que tu dépasses les capacités mémoire.

[Isabella83]

D'accord donc le mieux serait de parcourir avec quelque chose du genre :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
 
 
while (<$in>){
 
	if ($_=~/>/){
		$title=$_;
		$title=~s/>//;
	}
	else{
		$seq{$_}=$title;
	}
}

et ensuite, faire le travail sur ce hash ?

[Isabella83]

Merci beaucoup Stoyak ! Effectivement cela devait probablement venir du module Bio :: Seq, car en faisant comme ci dessous,en meme pas 1 minutes j'ai mon fichier de sortie !!

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
while (<$fasta>){
	chomp;
	if ($_=~/>/){
		$title=$_;
		$title=~s/>//;
	}
	else{
		$len=length($_);
		my $count = 0;
		while($len >= $int){
			$cut = substr($_,$count);
			push(@{$dico{$title}->{'sequence'}}, $cut);
			$len--;
			$count++;
 
 
		}
 
	}
}
 
print "Done 1 !\n";
 
foreach my $data (keys(%dico)){
	my $i=0;
	foreach my $seq (@{$dico{$data}->{'sequence'}}){
		print {$out} ">$data"."_nb$i"."\n$seq\n";
		$i++;
	}
}

[stoyak]

C'est parfait !

Un petit tag résolu peut-être ?!