Supprimer la redondance dans des sequences

**Isabella83** · 14/12/2012, 15h01

Bonjour à tous,
Je suis confrontée à un problème que je n'arrive pas à résoudre.
J'ai un fichier qui ressemble à cela :

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
>seq1
ACTTTCCACAACGATGGAAGATGATGA
>seq2
ACTTTCCACAACGATGGAAGATGATGAA
>seq3
ATTCCACAACGATGGAAGATGATGAA
>seq4
CTTTCCACAACGATGGAAGATGATGAA
>seq5
NTCCACAACGATGGAAGATGATGAAGA
>seq6
TACTTTCCACAACGATGGAAGATGATGA
>seq7
TACTTTCCACAACGATGGAAGATGATGAA
>seq8
TCCACAACGATGGAAGATGATGA
>seq9
AAAGAAGAAATTGAATAAATATATGTC
>seq10
AAAGAAGAAATTGAATAAATATATGT
>seq11
AAAGAAGAAATTGAATAAATATATGTCA
>seq12
AAAGAAGAAATTGAATAAATATAT
>seq13
AAAGAAGAAATTGAATAAATATATG
>seq14
AAAGAAGAAATTGAATAAATATA

Je souhaiterai trouver la séquence la plus petite que toutes ces séquences partagent entre elle (en gras dans mon exemple), et donc la seule séquence que j'aimerai récupérer est la seq 8 et seq14
>seq8
TCCACAACGATGGAAGATGATGA
>seq14
AAAGAAGAAATTGAATAAATATA

Quelqu'un aurait une idée ?

**Gardyen** · 14/12/2012, 16h39

hum petit pseudo algo:

classement des séquences par taille de la plus petite à la plus grande
je stocke la première séquence, forcément choisie puisque la plus petite
boucle de parcours des autres séquences:
- si la séquence contient l'une des séquences stockées (entièrement et exactement), l'enlever
- sinon la stocker

qu'en dis-tu ?

**Isabella83** · 14/12/2012, 16h58

Alors j'ai fait quelque choses comme ceci mais qui ne fonctionne pas :

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
#!/usr/bin/perl
use strict;
use warnings;
use Carp qw(confess);
use Getopt::Long;
use List::MoreUtils qw(uniq);
use Bio::SeqIO;
 
my %hash;
my ($fasta_file,$out_file);
GetOptions("fasta=s" => \$fasta_file,"out=s" => \$out_file);
open(my $out,'>',$out_file) or die "$out_file : $!\n\n";
my $in = Bio::SeqIO->new( -file => $fasta_file, '-format' => 'Fasta' );
 
# read file in hash
while ( my $seq = $in->next_seq()){
	my $id = $seq->id() ;
	my $sequence = $seq->seq ;
	$hash{$sequence}=$id;
}
 
my @all_seq=keys (%hash); # @seq contient toutes les sequences
my @uniqueseq;
my $find=0;
 
foreach my $seqs (@all_seq){
	$find=0;
	my $seq=uc($seqs);  #uppercase
	foreach (@uniqueseq){
		if ($_=~/$seq/){ # si ma sequence contient 
			$_=$seq; #remplace avec la plus petite sequence
			$find=1;
		}
	}
	if ($find==0){
		push @uniqueseq,$seq;
	}
}
 
foreach (@uniqueseq){
	print ">$hash{$_}\n$_\n";
}

je vais essayer de traduire en perl l'algo que tu m'as donné

**Isabella83** · 14/12/2012, 17h05

j'ai fait comme ceci :

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
44
45
46
#!/usr/bin/perl
use strict;
use warnings;
use Carp qw(confess);
use Getopt::Long;
use List::MoreUtils qw(uniq);
use Bio::SeqIO;
 
my %hash;
my ($fasta_file,$out_file);
GetOptions("fasta=s" => \$fasta_file,"out=s" => \$out_file);
open(my $out,'>',$out_file) or die "$out_file : $!\n\n";
my $in = Bio::SeqIO->new( -file => $fasta_file, '-format' => 'Fasta' );
 
# read file in hash
while ( my $seq = $in->next_seq()){
	my $id = $seq->id() ;
	my $sequence = $seq->seq ;
	$hash{$sequence}=$id;
}
 
my @all_seq=keys (%hash); # @seq contient toutes les sequences
my @sort_all_seq = sort @all_seq;
 
my @uniqueseq;
my $find=0;
 
foreach my $seqs (@sort_all_seq){
	$find=0;
	my $seq=uc($seqs);  #uppercase
	foreach (@uniqueseq){
		if ($_=~/$seq/){ # si ma sequence contient 
			$_=$seq; #remplace avec la plus petite sequence
			$find=1;
		}
	}
	if ($find==0){
		push @uniqueseq,$seq;
	}
}
 
my @final = uniq @uniqueseq;
 
foreach (@final){
	print ">$hash{$_}\n$_\n";
}

Ca fonctionne sur mon fichier exemple sur la sequence 8, mais pas pour la 14 ...
J'ai modifié comme ceci, ca à l'air de fonctionner à présent :

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
44
45
46
47
48
49
#!/usr/bin/perl
use strict;
use warnings;
use Carp qw(confess);
use Getopt::Long;
use List::MoreUtils qw(uniq);
use Bio::SeqIO;
 
my %hash;
my ($fasta_file,$out_file);
GetOptions("fasta=s" => \$fasta_file,"out=s" => \$out_file);
open(my $out,'>',$out_file) or die "$out_file : $!\n\n";
my $in = Bio::SeqIO->new( -file => $fasta_file, '-format' => 'Fasta' );
 
# read file in hash
while ( my $seq = $in->next_seq()){
	my $id = $seq->id() ;
	my $sequence = $seq->seq ;
	$hash{$sequence}=$id;
}
 
my @all_seq=keys (%hash); # @seq contient toutes les sequences
my @sort_all_seq = sort @all_seq;
 
my @uniqueseq;
my $find=0;
 
foreach my $seqs (@sort_all_seq){
	$find=0;
	my $seq=uc($seqs);  #uppercase
	foreach (@uniqueseq){
		if ($_=~/$seq/){ # si ma sequence contient 
			$_=$seq; #remplace avec la plus petite sequence
			$find=1;
		}
		if ($seq=~/$_/){ # si ma sequence contient 
			$find=1;
		}
	}
	if ($find==0){
		push @uniqueseq,$seq;
	}
}
 
my @final = uniq @uniqueseq;
 
foreach (@final){
	print ">$hash{$_}\n$_\n";
}

**Gardyen** · 14/12/2012, 17h25

je voyais plus simple

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 List::MoreUtils qw(uniq);
use Bio::SeqIO;
 
my %hash;
my ($fasta_file,$out_file);
GetOptions("fasta=s" => \$fasta_file,"out=s" => \$out_file);
open(my $out,'>',$out_file) or die "$out_file : $!\n\n";
my $in = Bio::SeqIO->new( -file => $fasta_file, '-format' => 'Fasta' );
 
# read file in hash
my $smallest_sequence;
while ( my $seq = $in->next_seq()){
	my $id = $seq->id();
	my $sequence = uc($seq->seq) ;
	$hash{$sequence}=$id;
 
	if (!$smallest_sequence) {
		$smallest_sequence = $sequence;
	}
 
	$smallest_sequence = $sequence unless length($smallest_sequence) < length($sequence);
}
 
my @all_seq=sort {length($a) <=> length($b)} keys (%hash); # @seq contient toutes les sequences
my @uniqueseq = ($smallest_sequence);
 
foreach my $seqs (@all_seq){
	my $seq=uc($seqs);  #uppercase
 
	my @found = grep {$seq =~ /$_/} @uniqueseq;
	if (!@found) {
		push @uniqueseq, $seq;
	}
}
 
foreach (@uniqueseq){
	print ">$hash{$_}\n$_\n";
}

le classement par ordre de taille permet d'être sûr qu'aucune séquence plus courte que celles dans unique_seq puisse exister avec la même séquence cible

**Isabella83** · 14/12/2012, 17h28

ah oui effectivement, ce sera surement moins long que ce que j'ai fait !
Merci beaucoup !

**Isabella83** · 04/02/2013, 09h41

Bonjour,
Je réécris dans ce post car je réutilise ce script, mais cela fait maintenant 4 jours qu'il tourne sur 2196958 séquences, sur un serveur ( disques : 6To brut, mémoire vive : 16Go).
J'aurai voulu savoir si cela était normal, ou si j'attendais pour rien ?
Comment puis je "optimiser" ce script pour que le résultat me soit retourné plus rapidement ?
je ne connais que Perl comme langage de programmation ... Faudrait il que je passe dans un autre langage de prog ?
Merci d'avance !

**trex7g2** · 04/02/2013, 10h38

Bonjour,

En théorie si ton script fonctionne sur un plus petit jeu de données, il n'y a aucune raison que tu te retrouves dans une boucle infinie.

Une lecture rapide me fait croire que tu fais d'assez gros tableaux... Je ne serais donc pas surpris que cela prenne du temps (et cela quelque soit ton langage).

Quand je sais qu'un script fa être long, je n'hésite pas a imprimer dans la console ou dans un fichier séparé des indicateurs de mon avancement : Par exemple, si tu es dans une boucle, tu peux mettre un compteur et faire un print tout les 10 000 ou 100 000 occurences :

Code perl :

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

1
2
3
4
					$i++;
					if ($i%100000==0){
						print "$i\n";
					}

Si tu as plusieurs étapes, tu peux tout simplement indiquer "début de la deuxieme étape" etc....

**Isabella83** · 04/02/2013, 10h44

oui justement , j'ai mit des print, juste après le hash ( il s'affiche au bout de quelques secondes), et ensuite, avant l'affichage (qui ne s'est toujours pas affiché) depuis 4 jours !
N'y a t'il pas une autre solution que d'utiliser des tableaux ?

**trex7g2** · 04/02/2013, 10h50

Il y a des chances que tu passes/aies passé beaucoup de temps dans le sort... Trier 2 millions de valeurs c'est assez long. Et cette fonction est en C donc impossible d'améliorer/optimiser cette étape.

Le grep aussi peut être long si tu as beaucoup de séquences uniques.

Tu as regardé ce qui est limitant de la mémoire ou du CPU ? Tu as vérifié que ton job était bien en train de tourner ? C'est peut être ton serveur qui sature sous le nombre de jobs aussi.

**Isabella83** · 04/02/2013, 10h52

Non je ne sais pas du tout comment on fait cela ...

**Isabella83** · 04/02/2013, 10h59

Est ce que je ne pourrai pas faire un sort avec une commande bash, creer un fichier trier et ensuite lire ce fichier et faire la suite du programme ? je ne sais pas si le sort unix est plus rapide ...

**trex7g2** · 04/02/2013, 10h59

Si tu es sur linux le plus simple c'est l'interface htop (je vais me faire tirer les oreilles par les pros la

) :

Dans ta console, lances la commande htop. Tu vas obtenir une interface graphique qui te permettra de suivre ton CPU et ta mémoire. htop a l'avantage d'être assez simples a utiliser. Tu peux te servir des flèches de ton clavier pour te diriger dans la liste des logiciels lancés. Plus d'info la (sourceforge). Il y a pleins d'explications sur le net à htop que tu devrais trouver assez facilement

.

**Isabella83** · 04/02/2013, 11h12

Alors j'ai fait 'htop' :
pour le script que j'ai lancé,
PRI : 20
VIRT 917M
RES : 884M
SHR : 1932
S : R
CPU : 99%
MEM : 11%

Je n'y comprends pas grand chose !

Sous tes conseils, j'ai rajouté un print après le sort, le sort se fait très rapidement, je pense que ce qui est long c'est le parcours de tous les tableaux ...

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
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
#!/usr/bin/perl
#	Input file: a fasta file
#	Output file: a unique fasta file
#	Command line: perl test.pl infile.fasta
 
use strict;
use warnings;
use List::MoreUtils qw(uniq);
#read the file into a hash
my %hash;
my $title;
my $infile=shift or die "give me a infile\n";
open (IN,"$infile");
	while (<IN>){
		$_=~s/\n//;
		$_=~s/\r//;
		if ($_=~/>/){
			$title=$_;
			$title=~s/>//;
		}
		else{
			$hash{$_}=$title;
		}
	}
 
print "hash ok \n";	
 
close IN;
#remove the abundant sequences
 
 
	my @all_seq=keys (%hash); # @seq contient toutes les sequences
	print "sort ok \n";
	my @sort_all_seq = sort @all_seq;
	my @uniqueseq;
	my $find=0;
	foreach my $seqs (@sort_all_seq){
		$find=0;
		my $seq=uc($seqs); #uppercase (retourne la chaine en majuscule)
		foreach (@uniqueseq){
			if ($_=~/$seq/){ # si ma sequence contient 
			$_=$seq; #remplace avec la plus petite sequence
			$find=1;
		}
		if ($seq=~/$_/){ 
			$find=1;
		}
	}
	if ($find==0){
		push @uniqueseq,$seq;
	}
	}
 
print "writing ....\n";
 
my @final = uniq @uniqueseq;
#outout the final result
open (OUT,">SF_file.fa");
foreach (@final){
	print OUT ">$hash{$_}\n$_\n";
}
close OUT;

**trex7g2** · 04/02/2013, 11h29

A l'oeil tu satures complètement ta vitesse de calcul, ce qui est plutôt bien (je préfère saturer le CPU à la mémoire... )

tu devrais aussi mettre un print après ton sort de la ligne 34;

Ensuite il faut voir ce que fait ton script : Supposons que le premier million de séquences soient uniques, de la séquence 1 000 001 à 2 000 000 tu vas devoir faire 1 000 000 de greps (séquences uniques) sur 1 000 000 de séquences... Sachant que grep est une fonction (très bien optimisée cela va de soit) qui demande un peu de temps de calcul...

**Isabella83** · 04/02/2013, 11h33

Oui, le "print sort" après la ligne 34 est vraiment très rapide (quelques secondes).
Donc je ne sais pas quoi penser ... comment optimiser car 4 jours sans aucun résultat c'est vraiment long je trouve

**trex7g2** · 04/02/2013, 11h44

Alors plus qu'une seule possibilité... Trop de greps! (si toutes tes séquences sont uniques tu arrives a plusieurs milliers de milliards de greps! si tu supposes que chacun prend 1/100eme de secondes, je te laisse faire le calcul du temps...)

tu peux eventuellement diviser ton fichier en 100, faire ton script sur tous, puis réunir ton fichier résultat et refaire ton script... Comme ça ton nombre de calcul cahngera de dimension

**Isabella83** · 04/02/2013, 11h49

oui j'avais penser à diviser mon fichier original en plusieurs fichier et faire un script bash pour appliquer ce script à tous ces fichiers temporaires ... Tu penses que cela diminuerai radicalement le temps de calcul ?

**Gardyen** · 04/02/2013, 12h14

on peut aussi ajouter quelques améliorations:
par exemple tu classes par ordre alphabétique tes séquences:

Code :

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

my @sort_all_seq = sort @all_seq;

, ce qui je l'avoue me semble inutile

un classement par taille réduirait les calculs puisqu'on est certain que les séquences suivantes sont forcément plus grandes que celles déjà traitées. du coup pas besoin de remplacer les séquences.

Code :

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

1
2
3
4
5
foreach (@uniqueseq){
			if ($_=~/$seq/){ # si ma sequence contient 
			$_=$seq; #remplace avec la plus petite sequence
			$find=1;
		}

ici tu vas parcourir toutes tes séquences même si tu en trouves une. un last après le $find=1; permettrait de réduire le nombre de séquences parcourues.

**Isabella83** · 04/02/2013, 12h17

ok comme ceci :

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
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
#!/usr/bin/perl
#	Input file: a fasta file
#	Output file: a unique fasta file
#	Command line: perl test.pl infile.fasta
 
use strict;
use warnings;
use List::MoreUtils qw(uniq);
#read the file into a hash
my %hash;
my $title;
my $infile=shift or die "give me a infile\n";
open (IN,"$infile");
	while (<IN>){
		$_=~s/\n//;
		$_=~s/\r//;
		if ($_=~/>/){
			$title=$_;
			$title=~s/>//;
		}
		else{
			$hash{$_}=$title;
		}
	}
 
print "hash ok \n";	
 
close IN;
#remove the abundant sequences
 
	my @sort_all_seq = sort {length($a) <=> length($b)} keys (%hash);
	print "sort ok \n";
	my @uniqueseq;
	my $find=0;
	foreach my $seqs (@sort_all_seq){
		$find=0;
		my $seq=uc($seqs); #uppercase (retourne la chaine en majuscule)
		foreach (@uniqueseq){
			if ($_=~/$seq/){ # si ma sequence contient 
			$_=$seq; #remplace avec la plus petite sequence
			$find=1;
			last;
		}
		if ($seq=~/$_/){ 
			$find=1;
                        # ajout d'un last ici ?
		}
	}
	if ($find==0){
		push @uniqueseq,$seq;
	}
	}
 
print "writing ....\n";
 
my @final = uniq @uniqueseq;
#outout the final result
open (OUT,">SF_banque_embryon_02h.fa");
foreach (@final){
	print OUT ">$hash{$_}\n$_\n";
}
close OUT;

dois je rajouter un "last là ou j'ai indiqué dans mon code ?

Supprimer la redondance dans des sequences

Bioinformatique Perl

Discussions similaires

Partager

Partager