parcours d'un tableau et suppression d'élements

Version imprimable

Voir 40 message(s) de cette discussion en une page

Bonjour à tous,
J'ai un problème que je ne parviens pas à résoudre.
J'ai un fichier qui ressemble à ceci :
Code:

1 2 3 4 5 6 7 8 9 10 11 AAAACCCGAAACCGTTTAGCTGT AAAAACCCGAAACCGTTTAGCTGT AAAACCCGAAACCGTTTAGCTGTT AAAACCCGAAACCGTTTAGCTGTTC AAAACCCGAAACCGTTTAGCTGTTT AAAACCCGAAACCGTTTAGCTGTTCT AGAAAAACCCGAAACCGTTTAGCTGT AAAACCCGAAACCGTTTAGCTGTTCTC AAAACCCGAAACCGTTTAGCTGTTCTCC AAAACCCGAAACCGTTTAGCTGTTCTCCG ATAACCCGAAACCGTTTAGCTGTTCTCCG
Tous les séquences sont rangées par taille et par ordre alphabétique dans mon fichier.

Je souhaite supprimer la redondance en supprimant toutes les séquences contenant entièrement une autre séquences (le début doit être identique )par exemple :
Par exemple, pour la première séquence :

Code:

AAAACCCGAAACCGTTTAGCTGT

j'ai toutes ces séquences :
Code:

1 2 3 4 5 6 7 AAAACCCGAAACCGTTTAGCTGTT AAAACCCGAAACCGTTTAGCTGTTC AAAACCCGAAACCGTTTAGCTGTTT AAAACCCGAAACCGTTTAGCTGTTCT AAAACCCGAAACCGTTTAGCTGTTCTC AAAACCCGAAACCGTTTAGCTGTTCTCC AAAACCCGAAACCGTTTAGCTGTTCTCCG
Mais je ne souhaite garder que

Code:

AAAACCCGAAACCGTTTAGCTGT

puisque toutes les autres séquences commencent par cette séquence.
A la fin, j'aimerai avoir :
Code:

1 2 3 4 AAAACCCGAAACCGTTTAGCTGT AAAAACCCGAAACCGTTTAGCTGT AGAAAAACCCGAAACCGTTTAGCTGT ATAACCCGAAACCGTTTAGCTGTTCTCCG
J'ai donc écris ce script perl :
Code:

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 my ($file_file,$out_file); my (%hash,%dico); my (@results,@sort_all_seq,@tmp); GetOptions("file=s" => \$file_file); open(my $file,'<',$file_file) or die "$file_file : $!\n\n"; while(<$file>){ chomp; push(@sort_all_seq,$_); } print "seq in array : OK !\n"; my $i=0; foreach my $value(@sort_all_seq){ push(@tmp,$value); my @dix = grep { $sort_all_seq[$_] =~ m/^$value/ } 0..$#sort_all_seq ; push(@tmp,$value) if(scalar(@dix)>1); my $o = 0; foreach my $val (@dix) { splice (@sort_all_seq,$val-$o,1); $o++; } $i++; } print "tmp : @tmp\n";
qui me retourne :
Code:

1 2 3 AAAACCCGAAACCGTTTAGCTGT AAAACCCGAAACCGTTTAGCTGT AGAAAAACCCGAAACCGTTTAGCTGT
Je pense que cela vient du fait que je modifie à chaque fois mon tableau @sort_all_seq.
Auriez vous une idée pour résoudre mon problème ?
Merci d'avance

06/10/2014, 19h18
Lolo78
Hmmm, modifier un tableau (ajouter ou retirer des éléments) alors que l'on est en train de le balayer avec une boucle for est une mauvaise idée.

Lis ton tableau et travaille sur une copie du même tableau. Quelque chose comme ceci:
Code:

1 2 3 4 5 6 my @result = @sort_all_seq; foreach my $value(@sort_all_seq){ @result = grep { not m/^$value/ } @result ; # note comment cette syntaxe est plus simple unshift @result, $value; }
On pourrait ajouter des optimisations, mais l'essentiel y est.

En fait, en y réfléchissant, l'optimisation à laquelle je pensais n'est pas seulement une optimisation, elle est nécessaire au bon fonctionnement de l'algo.

Je te propose ceci:

Code:

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
 
use strict;
use warnings;
use feature qw /say/;
 
my @sort_all_seq = split /\s+/, 
"AAAACCCGAAACCGTTTAGCTGT
AAAAACCCGAAACCGTTTAGCTGT
AAAACCCGAAACCGTTTAGCTGTT
AAAACCCGAAACCGTTTAGCTGTTC
AAAACCCGAAACCGTTTAGCTGTTT
AAAACCCGAAACCGTTTAGCTGTTCT
AGAAAAACCCGAAACCGTTTAGCTGT
AAAACCCGAAACCGTTTAGCTGTTCTC
AAAACCCGAAACCGTTTAGCTGTTCTCC
AAAACCCGAAACCGTTTAGCTGTTCTCCG
ATAACCCGAAACCGTTTAGCTGTTCTCCG";
 
 
my @result = @sort_all_seq;
my $old_val= "W"; # une valeur quelconque qui ne peut être reconnue dans les chaînes de bases
foreach my $value(@sort_all_seq){
     next if $value =~ $old_val; # l'optimisation nécessaire au bon fonctionnement
     @result = grep { not /^$value/ } @result ; # note comment cette syntaxe est plus simple
     unshift @result, $value; # on a éliminé toutes les chaines ressemblant à $value, il faut remettre $value
     $old_val = $value;
}
print "$_\n" for sort {length ($a) <=> length ($b)} @result;

Ce qui imprime:

Code:

1
2
3
4
5
6
$ perl tableau_genes.pl
 
AAAACCCGAAACCGTTTAGCTGT
AAAAACCCGAAACCGTTTAGCTGT
AGAAAAACCCGAAACCGTTTAGCTGT
ATAACCCGAAACCGTTTAGCTGTTCTCCG

@Lolo78, n'y aurait il pas intérêt à trier le tableau auparavant?

Quelque chose comme ça?
Code:

1 2 3 4 5 6 foreach my $val(sort @sort_all_seq) { if ($val !~ /^$old_value/) { print "$val\n"; $old_value = $val; } }

Oui, tu as raison. J'avais compris que le fichier d'origine était trié, mais il est fait trié par taille de séquence. Pour que la variable $old_val puisse pleinement jouer son rôle, il faudrait tries les données par ordre alphabétique. En revanche, comme le tableau @sort_all_seq n'est modifié à aucun moment, il faut le trier une seule fois, avant d'entrer dans la boucle for, et pas dans la bouclo for. Ce qui donne (pas testé):
Code:

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 use strict; use warnings; use feature qw /say/; my @sort_all_seq = sort split /\s+/, "AAAACCCGAAACCGTTTAGCTGT AAAAACCCGAAACCGTTTAGCTGT AAAACCCGAAACCGTTTAGCTGTT AAAACCCGAAACCGTTTAGCTGTTC AAAACCCGAAACCGTTTAGCTGTTT AAAACCCGAAACCGTTTAGCTGTTCT AGAAAAACCCGAAACCGTTTAGCTGT AAAACCCGAAACCGTTTAGCTGTTCTC AAAACCCGAAACCGTTTAGCTGTTCTCC AAAACCCGAAACCGTTTAGCTGTTCTCCG ATAACCCGAAACCGTTTAGCTGTTCTCCG"; my @result = @sort_all_seq; my $old_val= "W"; # une valeur quelconque qui ne peut être reconnue dans les chaînes de bases foreach my $value(@sort_all_seq){ next if $value =~ $old_val; # l'optimisation nécessaire au bon fonctionnement @result = grep { not /^$value/ } @result ; # note comment cette syntaxe est plus simple unshift @result, $value; # on a éliminé toutes les chaines ressemblant à $value, il faut remettre $value $old_val = $value; } print "$_\n" for sort {length ($a) <=> length ($b)} @result;
Mais je pense que ça marchera mieux et sans tri d'aucune sorte si on remplace ces deux lignes:
Code:

1 2 @result = grep { not /^$value/ } @result ; # note comment cette syntaxe est plus simple unshift @result, $value;
par :

Code:

@result = grep { $_ eq $value or not /^$value/ } @result ;

07/10/2014, 10h08
Isabella83

Merci.
En fait, je voulais supprimer les éléments de mon tableau au fur et à mesure car j'ai beaucoup de séquences (4millions de séquences uniques) et je voulais limiter le temps de calcul au maximum. J'avais essayé autre chose et le script à tourner 3 jours sans me donner de résultat ...
Je vais essayer avec votre solution et voir ce que ça donne !
07/10/2014, 13h53
Lolo78

Le tableau @result est épuré au fur et à mesure et l'optimisation utilisant $old_value devrait améliorer sensiblement les choses (mais ça dépend de tes données et de leur ordre).

Je serai curieux de savoir combien de temps ça prendra.

Bonne journée,

Laurent.

Une proposition de script que j'ai testée avec ton jeu d'essai, et qui semble fonctionner (elle n'utilise pas de regexp mais la fonction index, puisque dans ton cas, la simple recherche de sous-chaine est suffisante).

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
#!/usr/bin/perl
 
use strict;
use warnings;
 
use Getopt::Long;
 
my ($file_file);
 
GetOptions("file=s" => \$file_file);
open my $FILE, '<', $file_file or die "$file_file : $!\n\n";
 
chomp(my @sort_all_seq = <$FILE>);
print scalar(@sort_all_seq), " sequences in array : OK !\n";
 
my @root_seq;
while (@sort_all_seq) {
  push @root_seq, my $seq = shift @sort_all_seq;
  @sort_all_seq = grep index($_, $seq) != 0, @sort_all_seq;
}
print "Root sequences :\n", map "$_\n",@root_seq, "\n";

07/10/2014, 23h34
Lolo78

Tu as bien raison, Philou, de proposer d'utiliser la fonction index au lieu d'une regex, quand il s'agit de faire du "exact match" comme ici, c'est nettement plus rapide (même si ce n'est certainement pas cette partie du code qui prend le plus de temps, mais tant qu'à faire, c'est une optimisation qui ne coûte rien, autant ne pas s'en priver). J'y avais aussi pensé, mais elle ne m'avais pas paru si importante tant que j'ignorais le volume de données.

Cela dit, ce qui coûte cher ici, c'est de devoir passer en revue le tableau un très grand nombre de fois. Et ça, je vois difficilement un moyen de faire autrement, compte tenu de la nature du problème. Mais avec un bon algo éliminant les données inutiles rapidement, j'imagine que ça ne devrait pas prendre trois jours ou plus, mais la nature et la répartition exacte des données ont tellement d'importance sur les vraies performances qu'il est impossible de le dire sans avoir testé avec de vraies données.

@ Philou: j'ai tout de même un très gros doute sur ton algo. Modifier (brutalement, en plus, en fait le remplacer) le tableau @sort_all_seq alors que tu es en train de le parcourir avec une boucle while me paraît pour le moins extrêmement périlleux. Il se peut que ça marche, mais j'ai de très très gros doutes, en fait, je n'y crois guère, Je ne vois pas comment Perl pourrait gérer un truc pareil. Mais, bon, on ne sait jamais, si ça marche, j'aurais appris quelque chose de nouveau. Ma solution avec deux tableaux différents a le mérite, sauf erreur, de ne pas prendre ce genre de risque, même si elle est intrinsèquement un peu plus lente (mais pas tant que cela, à mon avis, parce que c'est la taille du tableau @result, épuré progressivement, qui gouverne probablement les performances de l'algo que j'ai proposé, beaucoup plus que celle du tableau en entrée) que la tienne, si la tienne marche.

@ Isabella: il est possible d'imaginer des solutions d'optimisations nettement plus pointues, par exemple en utilisant des hash sur les séquences de base pour éliminer plus rapidement les séquences que tu veux rejeter, ou autres, mais il faudrait que tu fournisses un échantillon nettement plus fourni de données en entrée (je dirais quelques milliers ou même plutôt quelques dizaines de milliers de séquences réelles) pour permettre de tester le fonctionnement et de mesurer les performances de diverses solutions avec des données quasi-réelles. Si tu peux attacher à un post un fichier de test réaliste, on pourra sans doute aller plus loin pour t'aider.

Bonne soirée à tous,
Laurent.
08/10/2014, 09h56
Isabella83

Oui cela fait 1 jour maintenant et toujours pas fini ... Je pense qu'il faut optimiser effectivement mais je ne suis pas une pro pour ça.
Je n'ai pas réussi à charger le fichier en pièce jointe mais il est téléchargeable à l'adresse suivante : http://dl.free.fr/getfile.pl?file=/kJYd6cJb
Merci pour votre aide !

3 poitns me semblent important: séquences triées par ordre alphabétique, triées par taille et comparaison des séquences mais doivent commencer pareil.

avec les 3 points précédents, je te propose d'essayer de traiter directement les lignes (pas testé):
Code:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 open (FH_IN, $file) or die "Impossible d'ouvrir le fichier $file\n"; my $taille = 10000000000; my $reference; my @Seq_conservee; while (<FH_IN>){ my $ligne = $_; $ligne =~/\s+$//; if ( $taille < length($ligne) ){ $reference = $ligne; push( @Seq_conservee, $ligne ); } else{ if ( $ligne !~ /^$reference/ ){ push( @Seq_conservee, $ligne ); } } $taille = length($ligne); }
j'ai mis un push mais tu peux aussi écrire les séquences dans un autre fi

08/10/2014, 10h51
Isabella83

Oui dans mon fichier j'ai trié les séquences par taille et par ordre alpha, de façon à avoir les plus courtes en premiers !
J'essaie ça et je vous tiens au courant.
08/10/2014, 11h36
6ril23

sinon une manière de faire plus bioinfo serait de réaliser un alignement de séquences et de filtrer les résultats (100% d'identité, ...)
08/10/2014, 11h37
Isabella83

Un alignement de séquence sur les séquences ?
08/10/2014, 11h39
6ril23

oui tu alignes les séquences (avec Blast par exemple) avec de paramètres assez restrictif

mais je pense que la méthode précédente si elle fonctionne sera meilleur et plus facile

je corrige mon message precedent

Code:

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
use warnings;
use feature qw /say/;
my @sort_all_seq = split /\s+/, 
"AAAACCCGAAACCGTTTAGCTGT
AAAAACCCGAAACCGTTTAGCTGT
AAAACCCGAAACCGTTTAGCTGTT
AAAACCCGAAACCGTTTAGCTGTTC
AAAACCCGAAACCGTTTAGCTGTTT
AAAACCCGAAACCGTTTAGCTGTTCT
AGAAAAACCCGAAACCGTTTAGCTGT
AAAACCCGAAACCGTTTAGCTGTTCTC
AAAACCCGAAACCGTTTAGCTGTTCTCC
AAAACCCGAAACCGTTTAGCTGTTCTCCG
ATAACCCGAAACCGTTTAGCTGTTCTCCG";
my @result = @sort_all_seq;
            my %tout;
            foreach my $v (@sort_all_seq){
            my $n = length($v); 
            my $s = $v;
            for (my $var = $n; $var > 0; $var--) {
 
                    $s = substr($s,0,$var);
                    if ( (exists( $tout{$s} ))) {
                            my @tob = $tout{$s};
                            push($tob[0], $v);
                            #unshift @sort_all_seq;
                            last;
                    }   
            }   
            if ( (!exists( $tout{$s} ))) {
            my @tab;
            push(@tab, $v);
            $tout{$v} = [@tab];
    }
 
    }   
    foreach my $k (keys(%tout)) {
       my @valeur=$tout{$k};
       my $long = @{$valeur[0]};
       print "$k est la clef et et la taille est  $long\n";
       my $a1 = ${valeur}[0][0];
       my $a2 = ${valeur}[0][1];
       my $a3 = ${valeur}[0][2];
       my $a4 = ${valeur}[0][3];
       if ($a1) {say $a1, " valeur 1";}
       if ($a2) {say $a2, " valeur 2";}
       if ($a3) {say $a3, " valeur 3";}
       if ($a4) {say $a4, " valeur etc";}
       say $long;
       }

affiche
AAAACCCGAAACCGTTTAGCTGT est la clef et et la taille est 8
AAAACCCGAAACCGTTTAGCTGT valeur 1
AAAACCCGAAACCGTTTAGCTGTT valeur 2
AAAACCCGAAACCGTTTAGCTGTTC valeur 3
AAAACCCGAAACCGTTTAGCTGTTT valeur etc
8
AGAAAAACCCGAAACCGTTTAGCTGT est la clef et et la taille est 1
AGAAAAACCCGAAACCGTTTAGCTGT valeur 1
1
AAAAACCCGAAACCGTTTAGCTGT est la clef et et la taille est 1
AAAAACCCGAAACCGTTTAGCTGT valeur 1
1
ATAACCCGAAACCGTTTAGCTGTTCTCCG est la clef et et la taille est 1
ATAACCCGAAACCGTTTAGCTGTTCTCCG valeur 1
1

Je ne l'ai pas encore fait tourner, j'ai une erreur en fait, et j'ai dû faire autre chose entre temps !
Code:

1 2 Use of uninitialized value $reference in regexp compilation at smaller_family_test.pl line 24, <$FILE> line 1 Use of uninitialized value $reference in regexp compilation at smaller_family_test.pl line 24, <$FILE> line 2
J'ai juste modifié :

Code:

$ligne =~/^(\s)*$/;

Car j'avais une erreur

08/10/2014, 12h08
6ril23

Citation:

Envoyé par Isabella83

J'ai juste modifié :

Code:

$ligne =~/^(\s)*$/;

Car j'avais une erreur

en fait il faut modifier la ligne 7 de mon code:
$ligne =~/\s+$//; --> $ligne =~ s/\s+$//;
c'est à dire enlever les éventuels espaces, tabulations en fin de ligne, ainsi que le retour à la ligne

Ta ligne de code ci-dessus peut servir si tu as des lignes vides (cas que je n'ai pas pris en compte)
08/10/2014, 12h10
Isabella83

OK
Mais j'ai toujours l'erreur que je ne comprends pas
08/10/2014, 12h55
Philou67430

Citation:

Envoyé par Lolo78

@ Philou: j'ai tout de même un très gros doute sur ton algo. Modifier (brutalement, en plus, en fait le remplacer) le tableau @sort_all_seq alors que tu es en train de le parcourir avec une boucle while me paraît pour le moins extrêmement périlleux. Il se peut que ça marche, mais j'ai de très très gros doutes, en fait, je n'y crois guère, Je ne vois pas comment Perl pourrait gérer un truc pareil. Mais, bon, on ne sait jamais, si ça marche, j'aurais appris quelque chose de nouveau. Ma solution avec deux tableaux différents a le mérite, sauf erreur, de ne pas prendre ce genre de risque, même si elle est intrinsèquement un peu plus lente (mais pas tant que cela, à mon avis, parce que c'est la taille du tableau @result, épuré progressivement, qui gouverne probablement les performances de l'algo que j'ai proposé, beaucoup plus que celle du tableau en entrée) que la tienne, si la tienne marche.

Je n'ai aucun doute sur le fait que cela fonctionne, car la boucle "while" ne travaille pas sur le tableau, elle ne fait qu'évaluer à chaque boucle si le tableau est vide ou non (c'est une technique que j'ai déjà utilisée par le passé, et peut-être même pour des séquences d'ADN). Ce n'est exactement comme si l'on utilisait un itérateur sur le tableau pendant qu'on le modifie.

Je vais lancer un test avec le fichier d'Isabella.

Voir 40 message(s) de cette discussion en une page