Discussion :

[trex7g2]

un last en ligne 46 serait une amélioration oui

[Gardyen]

ah il y a un problème d'accolade !
il en manque une à la ligne 47, du coup le if de la ligne 49 est répété pour chaque séquence unique testée !!! au lieu de l'être uniquement si il n'a rien trouvé après avoir parcouru toutes les séquences uniques.

donc tu te retrouves avec une augmentation exponentielle de tes séquences uniques, ce qui doit causer tes problèmes ^^
au temps pour moi j'ai mal lu l'indentation

ensuite vu que tu as trié tes séquences par taille, aucune des séquences testées ne sera contenue dans une séquence unique, à la limite tu pourras avoir 2 séquences identiques en cas de duplicat, mais pas plus => plus besoin de remplacer les séquences.
tu peux donc commenter/effacer le premier if, du coup il faut ajouter un last dans le if restant

[trex7g2]

je viens de relire ton code brièvement et il me semble qu'il y a une erreur :
ligne 39 :

Code :

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

if ($_=~/$seq/)

or si j'ai bien suivi, $_ correspond à ton jeu de séquences uniques. Par conséquent c'est le motif que tu cherches. Tu devrais par conséquent faire le contraire non ?

Code :

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

if ($seq=~/$_/)

[Lolo78]

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.

Je n'ai pas cherché à comprendre le code, mais trier 2 millions d'enregistrements, ce n'est pas si long. Aujourd'hui, au boulot, j'avais 35 millions d'enregistrements, le tri dure moins de 15 minutes. Pour moi, le problème n'est pas là.

[Isabella83]

J'ai essayé de supprimer le premier if, mais cela ne me retourne pas le bon résultat, donc je pense qu'il est indispensable ...
Cela me renvoit seulement la séquence n°7 de mon exemple ...

[djibril]

Bonjour,

Après lecture très rapide, pourquoi ne pas procéder en deux fois ?
C'est-à-dire lire le fichier deux fois.
Une première fois pour trouver les séquences les plus courtes (en stockant les id) et une seconde lecture pour les récupérer et mettre dans un autre fichier.

Cela permet de ne pas stocker en mémoire les séquences car je ne connais pas la taille de tes fichiers.

[Isabella83]

Bonjour,
Je ne suis pas sure de bien comprendre, tu suggères de couper le script en deux ? Mais cela serait plus rapide ?
j'ai plus de 2 millions de séquences à analyser

[djibril]

Heu, j'avais pas bien lu l'énoncé, tu ne connais pas par avance la similarité entre tes séquences ?

Dans ce cas, je ne comprends pas très bien ton exemple ? Tu peux reexpliquer ?

[Isabella83]

Bonjour,
Non je ne connais pas la similarité entre mes séquences;
J'ai un jeu de données, parmi lequel certaines séquences partagent un morceau de séquence en commun comme l'exemple que j'avais donné :

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
>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 les seules séquences que j'aimerai récupérer sont les seq 8 et seq14
>seq8
TCCACAACGATGGAAGATGATGA
>seq14
AAAGAAGAAATTGAATAAATATA

[djibril]

Entre la séquence 8 et 14, y a rien en commun. Donc, comment tu différencies ?
Pour moi, soit tu cherches une similarité entre toutes les séquences, soit on trouvera tout et n'importe quoi.

[Isabella83]

Non il n'y a rien entre la séquence 8 et 14 , mais ce sont les plus petites séquences, la séquence 8 est contenue dans les séquences 1 à 7 et la séquence 14 est contenu dans les séquence 9 à 13 ...

[djibril]

ça j'avais compris !

Je suis juste en train de te demander ton algorithme .
Car de ce que je vois de ton exemple, si ton fichier contient 10 séquences complétement différentes, tu aura donc 10 séquences en résultat, on est bien d'accord ?

Car pour faire ça, il faut comparer les séquences deux par deux, c'est ce que tu fais ?

[stoyak]

Bonjour,

Donc si je comprends bien, tu veux récupérer un nombre minimal de séquences dans ton jeu de données. Et pas

plus petite que toutes ces séquences partagent entre elle

.... Car dans ce cas, et comme le dit djibril, tu ne récupérerais aucune des séquences.

Si j'ai bien compris ce que tu cherches à faire, voici comment je procèderai :
- Lecture ligne à ligne des séquences ;
- Stockage de la première séquence dans un hash $1;
- Passage à la deuxième ligne : comparaison de la séquence $2 avec celle(s) stockée(s) ;
- Si la séquence $2 est plus courte et comprise dans celle stockée, ajout de s2 dans le hash et suppression de $1;
- Si la séquence $1 est plus longue et comprise dans $2, conservation de $1 et passage à la ligne suivante ;
- Si les séquences sont différentes, stockage de $2.

Cela te permet d'effectuer les comparaisons ligne à ligne sans stocker l'ensemble des séquences.

[trex7g2]

si j'ai bien compris dans ton fichier tu as des "familles" de séquences ayant toutes un bot en commun, et tu cherches à n'avoir que la signature de ces familles, c'est à dire la partie commune.

Dans l'algo que l'on tentait de faire avec gardyen, le but était donc de trier les séquences par taille, de manière a avoir au début les séquences les plus courtes, c'est à dire celles qui ne contiennent que la signature de leur famille respective. Puis on vérifie que toutes les autres séquences appartiennent à une famille que l'on connait déjà.

je ne suis pas sur que ce soit beaucoup plus clair

[Gardyen]

j'ai réécrit le premier script que j'avais donné, à la rigueur tu peux changer la ligne avec le grep avec le foreach/last pour éviter de parcourir toutes les séquences uniques

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
#!/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 = uc($seq->seq) ;
	$hash{$sequence}=$id;
 }
 
my @all_seq=sort {length($a) <=> length($b)} keys (%hash); # @seq contient toutes les sequences
my @uniqueseq = ();
 
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";
}

[Isabella83]

Merci Gardyen, j'ai mit le script a tourner, wait and see !
stoyak, je vais essayer d’écrire le script correspondant à l'algo et voir si c'est plus court en terme de temps !