Traitement fichier 2 go

**JP2mars** · 13/10/2009, 09h30

Bonjour,

je dois traiter de très gros fichiers de séquences, ~ 2go. Fichiers issus de séquençage.
Je veux compter le nombre de A,T,C,G,N pour chaque cycle, c'est à dire à chaque position, (exemple 36) et cela pour chaque tile (100 tiles en tout). Je prend en compte le fait qu'il puisse y avoir jusqu'à 8 fichiers.

voilà un extrait de mon code :

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
 
 
my %nbycycle;
 
while(<ENTREE>){
 
		my @ligne = split /\t/, $_;
 
		my @sequence = split //, $ligne[8];
		my $i=1;
 
		foreach my $seq (@sequence){
 
			$nbycycle{$ligne[2]}{$ligne[3]}{$i}{$seq}++; #{fichier}{tile}{cycle}{base}
			$i++;
		}
 
	}

Le résultat est ok mais je voudrais que ce soit plus rapide (~25 min pour l'instant pour un fichier). Donc y aurai t'il un moyen d'améliorer l'algo ou tester avec un autre langage ou autre chose.

Merci.

**djibril** · 13/10/2009, 10h18

montre nous quelques lignes de ton fichier.

**JP2mars** · 13/10/2009, 10h26

HWI-EAS387 90810 5 1 10 89 0 1 GGGGCCAGGGAGTGAGATCAGGAAGGGGACAGGCAA `````aa_`aa^X`a`\a_XXU_^a`aUGXYJ]aaa chr16.fa 85230887 F 36 117 Y
HWI-EAS387 90810 5 1 10 1835 0 1 TGGGAGGCCCAGGACACATGCCAGAGGAGGGATTAG a__Ra]OO[I_P\GGN_aR^^`^`\\_``^WW\P\[ chr3.fa 128241989 F 3T9C22 81 N
HWI-EAS387 90810 5 1 10 198 0 1 CACTTTCCCCTGTCTTTCAGCTCTGTCCGTGTGCAG a`aXZ`a]baa]`[aba`aaaa^NYXaYFVDNO^^V chr22.fa 16170918 F 36 114

Y

Voilà un exemple, ce qui m'intéresse c'est la 4ème colonne et la 9ème biensur.

**50Nio** · 13/10/2009, 14h42

Bonjour,
Globalement, tu gagnera peu ou rien en écrivant dans un autre langage, Perl est compilé puis exécuté en tant que code machine, donc tu pers un peu de temps au lancement du script, pas à l'exécution. Dans ton cas, si c'est pour gagner 1ms...

Sur un fichier de 1,8Giga avec ton code ça me prend 2mn.
Sortir les "my" de la boucle ne change rien (on est toujours à 2mn +/- 2s).

En utilisant des tableaux plutôt que des hash on retombe aussi sur la même vitesse de traitement.

Les idées simples ne donnent rien, c'est ptêt une question à poser au niveau algorithme...

**JP2mars** · 13/10/2009, 16h05

Envoyé par 50Nio

Bonjour,
Globalement, tu gagnera peu ou rien en écrivant dans un autre langage, Perl est compilé puis exécuté en tant que code machine, donc tu pers un peu de temps au lancement du script, pas à l'exécution. Dans ton cas, si c'est pour gagner 1ms...

Sur un fichier de 1,8Giga avec ton code ça me prend 2mn.
Sortir les "my" de la boucle ne change rien (on est toujours à 2mn +/- 2s).

En utilisant des tableaux plutôt que des hash on retombe aussi sur la même vitesse de traitement.

Les idées simples ne donnent rien, c'est ptêt une question à poser au niveau algorithme...

2 min ??
Moi il met 20min (au lieu de 25 comme j'ai dit au début, j'ai fait du vide dans le script).
Mais j'ai testé seulement la partie que je vous ai montré et elle me prend ~20min.
Il faut dire que je teste ça sur mon pc de bureau qui a 4 Mo de ram.
En routine il tournera sur un serveur de calcul beaucoup plus costaud.

Merci de vos réponses.

**djibril** · 13/10/2009, 16h17

Envoyé par JP2mars

4 Mo de ram.

, ça existe encore ? Tu m'étonnes qu'il mette 20 minutes

**JP2mars** · 13/10/2009, 17h25

Envoyé par djibril

, ça existe encore ? Tu m'étonnes qu'il mette 20 minutes

Go, autant pour moi.

**50Nio** · 14/10/2009, 08h27

Pour les 2mn, oui ça met 2mn, mais sur un serveur (32Go de RAM, 8CPUs ...).

Après, c'est pas le temps qui est important, c'est la variation avec diverses méthodes, et là je n'ai pas réussi à avoir un gain significatif.

**JP2mars** · 14/10/2009, 09h14

Merci pour toutes ces réponses !
Sachant que le serveur a des propirétés proches de celui de 50Nio, je vais le faire tester.
J'ai pas accès à ce serveur

.
Si ça prend alors 2 min, ça sera bon.

Encore merci.

**JP2mars** · 28/10/2009, 14h31

Bonjour,

je suis de retour avec mon problème d'optimisation.
Le test a été fait fait sur un serveur plus puissant que mon pc (il vaut mieux pour le serveur)
mais peu de gain niveau temps. C'est pas moi qui est testé, j'ai juste eu les retours.

J'ai fait des motifs dans mon code et de près de 25 min je suis passé à 15 min (même fichier même pc). J'ai modifié des bouts de code autre que celui montré. Mais il faut que ça aille encore plus vite (pas pour moi, moi ça me va

).

J'ai testé la mise en mémoire sur des plus petit fichier (car le serveur ayant plus de mémoire les gros fichiers pourraient être traité différemment) mais rien de gagner.

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
 
my (%N,%A,%C,%G,%T);
 
while(<ENTREE>){
 
		chomp($_);	
		my @ligne = split /\t/, $_;
		my $sequence = $ligne[8];
 
 
 
		my $e = 0;
		my $n = 0;
		while (1) {
			$e = index($sequence, "A", $e);
			last if ($e == -1);
			$e++;				
			$A{$ligne[3]}{$e}++;
 
		}
		$e = 0;
		$n = 0;
		while (1) {
			$e = index($sequence, "C", $e);
			last if ($e == -1);
			$e++;				
			$C{$ligne[3]}{$e}++;
 
		}
		 $e = 0;
		 $n = 0;
		while (1) {
			$e = index($sequence, "G", $e);
			last if ($e == -1);
			$e++;				
			$G{$ligne[3]}{$e}++;
 
		}
		 $e = 0;
		 $n = 0;
		while (1) {
			$e = index($sequence, "T", $e);
			last if ($e == -1);
			$e++;				
			$T{$ligne[3]}{$e}++;
 
		}
		 $e = 0;
		 $n = 0;
		while (1) {
			$e = index($sequence, "N", $e);
			last if ($e == -1);
			$e++;				
			$N{$ligne[3]}{$e}++;
 
		}
 
 
	}
 
	print "fichier $file ok\n";
	close(ENTREE);

Donc si vous avez une idée dévastatrice , mis à part avoir un pros de 100Go

. je suis preneur.

Merci encore de prendre le temps de repondre.

**djibril** · 28/10/2009, 15h54

Montre nous un exemple de ligne de ton fichier et un exemple de résultat que tu souhaites.

**JP2mars** · 28/10/2009, 16h36

Code :

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

1
2
3
HWI-EAS387	90810	5	100	9	883	0	1	AAGCGGTCCATCCTCACTCAGGCTCTAACTTCTTTC	BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB	NM											N
HWI-EAS387	90810	5	100	10	428	0	1	AGCCCTCCATCAATATGGGGGAACTACACTTTAGAT	^Xa\Tab_H]O[J^K`QGT^NHM[^\\LY_BBBBBB	NM											N
HWI-EAS387	90810	5	100	10	420	0	1	GGGGAGGCAAGTATCTCTGTATGGAGAGCTGGTCTC	^QZMYHTZ`X\JHO^^^[I]WOHZQYIQ_JYOX``b	NM											N

Voilà un exemple, ce qui m'intéresse c'est la 4ème colonne et la 9ème biensur.

Résultat :

Code :

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

1
2
3
4
5
6
7
8
9
10
Tile	Cycle	A	C	G	T	N
1	1	35879	26146	25766	35577	370
1	2	37203	23508	26328	35994	705
...
1	36	35592	25316	23325	37755	1750
2	1	34777	24895	26191	35063	345
...
2	36	35930	24834	24191	34910	1406
...
100  36     .....

Merci

PS: les symbols à la suite de la séquence c'est du code. Chaque signe correspond au score qualité.
Je devrais traiter ça en plus après.

**djibril** · 28/10/2009, 17h04

Bon je ne comprends pas ton résultat. Je pensais que tu souhaitais juste compter le nombre d'acides nucléique pour chaque ligne de ton fichier (chaque séquence)!! Or là, si je prends :

Code :

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

1
2
Tile	Cycle	A	C	G	T	N
1	1	35879	26146	25766	35577	370

Qu'est ce que ça veut dire ?

**djibril** · 28/10/2009, 17h13

Ton code est résumé à ceci :

Code :

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

1
2
3
4
5
6
7
8
9
10
11
12
13
my ( %N, %A, %C, %G, %T );
while ( my $ligne = <ENTREE> ) {
 
  chomp($ligne);
  my ( undef, undef, undef, $TOTO, undef, undef, undef, undef, $sequence ) = split( /\t/, $ligne );
  my $A{$TOTO} = () = $sequence =~ m{A}g;
  my $C{$TOTO} = () = $sequence =~ m{C}g;
  my $G{$TOTO} = () = $sequence =~ m{G}g;
  my $T{$TOTO} = () = $sequence =~ m{T}g;
  my $N{$TOTO} = () = $sequence =~ m{N}g;
}
close(ENTREE);
print "fichier $file ok\n";

Est ce que ça te convient, ai je compris ton souci ?

**JP2mars** · 28/10/2009, 17h19

Je vais tester.

Les données sont issus d'un séquenceur solexa. Le but est de vérifier pour chaque cycle le nombre de A,C,G,T,N.
Au sein d'une même séquence ça ne m'intéresse pas en plus j'ai plein d'exemple à ce sujet et je me serait pas permi de poster sinon.
Donc pour mes 15 Millions de séquences il faut que pour chaque position (1 à 36) je compte le nombre de nucléotides.

Si je trompe le code que tu me donne c'est pour compter le nombre de ACTG pour une même séquence, non?

Merci encore !

**JP2mars** · 28/10/2009, 17h21

Envoyé par djibril

Bon je ne comprends pas ton résultat. Je pensais que tu souhaitais juste compter le nombre d'acides nucléique pour chaque ligne de ton fichier (chaque séquence)!! Or là, si je prends :

Code :

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

1
2
Tile	Cycle	A	C	G	T	N
1	1	35879	26146	25766	35577	370

Qu'est ce que ça veut dire ?

Une piste sur un séquenceur est divisé en 100 cases appelé Tile.
Pour chacune de ces Tiles je dois pour chaque cycle compter le nombre de ACTGN.
La position 1 de la séquence correspond au premier cycle.

En tout j'ai 100*36 ligne pour lesquels j'ai mon nombre de ACTGN.

**djibril** · 28/10/2009, 17h22

Envoyé par JP2mars

Je vais tester.

Les données sont issus d'un séquenceur solexa. Le but est de vérifier pour chaque cycle le nombre de A,C,G,T,N.
Au sein d'une même séquence ça ne m'intéresse pas en plus j'ai plein d'exemple à ce sujet et je me serait pas permi de poster sinon.
Donc pour mes 15 Millions de séquences il faut que pour chaque position (1 à 36) je compte le nombre de nucléotides.

Ca veut que si ta séquence fait 40 acides nucléiques, tu t'en fou des 4 derniers ? Si oui, est ce que ça veut dire que pour toute séquence, tu souhaites compter le nombre de ATCGN qu'il y a pour la sous séquence allant de la position 1 à 36 ? as tu un exemple ?

Envoyé par JP2mars

Si je trompe le code que tu me donne c'est pour compter le nombre de ACTG pour une même séquence, non?

oui

**JP2mars** · 28/10/2009, 17h25

Envoyé par djibril

Ca veut que si ta séquence fait 40 acides nucléiques, tu t'en fou des 4 derniers ? Si oui, est ce que ça veut dire que pour toute séquence, tu souhaites compter le nombre de ATCGN qu'il y a pour la sous séquence allant de la position 1 à 36 ? as tu un exemple ?

oui

Un run de séquençage dure 36 cycles, donc mes séquences feront toujours 36 de longueurs.

Sachant que bientôt ce sera 72 puis 100. Mais dans les cas toutes les séquences feront la même taille.

**djibril** · 28/10/2009, 17h31

On ne se comprends peut être pas.
Voici ce que tu nous as donné comme exemple de fichier :

Code :

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

1
2
3
HWI-EAS387	90810	5	100	9	883	0	1	AAGCGGTCCATCCTCACTCAGGCTCTAACTTCTTTC	BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB	NM											N
HWI-EAS387	90810	5	100	10	428	0	1	AGCCCTCCATCAATATGGGGGAACTACACTTTAGAT	^Xa\Tab_H]O[J^K`QGT^NHM[^\\LY_BBBBBB	NM											N
HWI-EAS387	90810	5	100	10	420	0	1	GGGGAGGCAAGTATCTCTGTATGGAGAGCTGGTCTC	^QZMYHTZ`X\JHO^^^[I]WOHZQYIQ_JYOX``b	NM											N

Tes séquences dans ton fichier font toutes 36 acides nucléiques, Ok pas de soucis.

si je prends cette exemples : AAGCGGTCCATCCTCACTCAGGCTCTAACTTCTTTC

voici ce qu'il y a :

A : 7/36 (19%)
T : 13/36 (36%)
C : 5/36 (14%)
G : 11/36 (31%)

Mais toi, ça ne te suffit pas car ce résultat correspond au cycle 36, c'est ça ?

Peux tu nous donner un exemple de résultat avec cette séquence, histoire qu'on se comprenne mieux.

Merci

**JP2mars** · 28/10/2009, 17h49

Code :

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

1
2
3
4
5
HWI-EAS387	90810	1	1	10	1773	0	1	ACAGAGTGCCATGCTGGGGCTGTGCCTTCATGGGAG	``a[a``a``aa__a\_a_`aZXRW[GYV\]VOYPQ	chr21.fa		43058980	R	T35	80						Y
HWI-EAS387	90810	1	1	10	754	0	1	GACCTGCCCACCCTACTAACACCCCGGTACCCTCCC	``aaa_aaaaaaaaaaaaaaaaaaaaa`aaaaaaaa	chr13.fa		24191148	F	36	118						Y
HWI-EAS387	90810	1	1	10	107	0	1	GGGCAGACTTTGCAGGGTTGTTTCAGGCTAAATTGG	aabK_S_b`b_I[aaab`bba]aabaaba`[G\aZZ	NM											N
HWI-EAS387	90810	1	1	10	707	0	1	ACGGGCTAGGAAACAACACCAGTAACCCGCTGTTGG	IIZDZ[HWQG^^X`BBBBBBBBBBBBBBBBBBBBBB	NM											N
HWI-EAS387	90810	1	1	10	1975	0	1	GATTGACAATGATTTTTATTATGATGTGACTTCTCA	a`baaaaaaba`aabaaaaa`a`aa__^[`aa`a`_	NM

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
 
Tile    Cycle   A       C       G       T       N
1       1       2       0       3       0       0
1       2       2       2       1       0       0
1       3       1       1       2       1       0
1       4       0       2       2       1       0
1       5       2       0       2       1       0
1       6       1       1       3       0       0
1       7       1       2       0       2       0
1       8       2       2       1       0       0
1       9       1       2       1       1       0
1       10      1       1       1       2       0
1       11      2       1       1       1       0
1       12      2       1       1       1       0
1       13      1       2       1       1       0
1       14      1       2       0       2       0
1       15      2       0       1       2       0
1       16      1       1       2       1       0
1       17      0       1       2       2       0
1       18      3       0       1       1       0
1       19      1       1       1       2       0
1       20      0       3       1       1       0
1       21      3       0       0       2       0
1       22      0       1       2       2       0
1       23      0       1       1       3       0
1       24      2       2       1       0       0
1       25      2       2       0       1       0
1       26      0       2       3       0       0
1       27      0       1       2       2       0
1       28      0       2       1       2       0
1       29      2       1       1       1       0
1       30      2       3       0       0       0
1       31      1       1       0       3       0
1       32      1       1       2       1       0
1       33      0       1       1       3       0
1       34      0       1       1       3       0
1       35      1       2       2       0       0
1       36      1       1       3       0       0

A la première position j'ai donc 2 A et 3 G.
A la seconde j'ai 2 A, 2 C, et 1 G.

L'exemple que tu me donnes correspond seulement à une séquence.
Ce que je souhaite c'est que pour toutes les séquences de mon fichiers, à une position donnée (qui correspond donc à un cycle) je compte le nombre de ACTGN (N non trouvé). Et cela en prenant en compte le numéro du tile de la colonne 4.

J'espère être plus clair.

Traitement fichier 2 go

Bioinformatique Perl

Discussions similaires

Partager

Partager