Discussion :

[Isabella83]

Bonjour à tous,
Je viens sur le forum C car j'ai un problème.
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
>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

Code :

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

1
2
3
4
>seq8
TCCACAACGATGGAAGATGATGA
>seq14
AAAGAAGAAATTGAATAAATATA

J'ai ce script perl qui fonctionne très bien sur un petit jeu de données:

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,">tmp.fa");
foreach (@final){
	print OUT ">$hash{$_}\n$_\n";
}
close OUT;

Mais malheureusement, j'ai un fichier avec plus de 2 millions de séquences, et cela fait 4 jours que mon script toune sur un serveur ( disques : 6To brut, mémoire vive : 16Go) ... je trouve que cela est très long, et je voulais savoir si vous pensiez qu'en langage C , cela pouvait aller plus vite ?
Pb : je ne sais pas coder en C !

[Bktero]

Pb : je ne sais pas coder en C !

En partant du principe que le script en Perl est bien écrit, avec les bonnes optimisations algorithmiques, alors la réponse est probablement "non". En effet, il faudra sacrément optimiser le code C pour le faire aller vraiment plus vite et si tu es débutante, ça risque d'être compliqué.

[phi1981]

2 millions de séquences .... même en comptant 64 octets par séquence, ça ne fait jamais que 128 Mo de données.
4 jours pour les traiter me parait beaucoup. Ton script ne tournerait il tout simplement pas en boucle ?

De toute façons, avant de savoir si ça irait mieux en C, tu dois savoir où ton programme perd du temps. S'il s'agit de l'accès disque, réecrire ton programme en C ne va pas te faire gagner grand chose.
Si c'est le traitement en lui même qui est gourmand, peut être utiliser un algorithme différent aurait un sens.

[Obsidian]

Bonsoir,

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

Si tu es sûre qu'au moins une des séquences de ton fichier est contenue en entier dans toutes les autres, il te suffit de les trier par taille, voire simplement examiner la taille de chacune d'elle sans les classer, et de ne conserver que celles qui ont la taille minimale. De là, tu recherches uniquement ces chaînes au sein des autres. Tu passes alors d'une complexité exponentielle à une complexité linéaire, et tu gagneras beaucoup de temps. Cela dit, il y a peu de chances pour que ce soit réellement ce que tu cherches à faire.

S'il s'agit de trouver un motif commun, c'est plus difficile. Et dans ce cas, il faut rechercher non pas la plus courte séquence commune mais la plus longue. Sinon, la plus courte séquence serait forcément « A » , « C », « G » ou « T » et ne contiendrait qu'une paire de bases.

Si, en plus, tes séquences peuvent contenir des insertions, des délétions ou des SNPs, vois du côté de Blast. Cela dit, j'ai quitté le milieu de la bio-info il y a un moment et je ne sais plus quel est l'état de l'art en la matière.

Mais malheureusement, j'ai un fichier avec plus de 2 millions de séquences, et cela fait 4 jours que mon script toune sur un serveur ( disques : 6To brut, mémoire vive : 16Go) ... je trouve que cela est très long, et je voulais savoir si vous pensiez qu'en langage C , cela pouvait aller plus vite ? Pb : je ne sais pas coder en C !

Il y a des choses que tu pourrais effectivement faire directement, comme mapper des fichiers en mémoire mais il y a probablement des optimisations à faire au préalable sur ton programme. Ça dépend de ta problématique exacte.

En partant du principe que le script en Perl est bien écrit, avec les bonnes optimisations algorithmiques, alors la réponse est probablement "non". En effet, il faudra sacrément optimiser le code C pour le faire aller vraiment plus vite et si tu es débutante, ça risque d'être compliqué.

Ça dépend : il n'y a pas que l'algorithmique. Il faut également tenir compte des capacités de la machine quand on fait du traitement à grand échelle. Les bio-informaticiens au taquet dans les deux disciplines sont rares et c'est normal. On trouve par conséquent beaucoup de scripts Perl ou Python qui ne sont pas forcément sales en eux-mêmes, mais qui pourraient être beaucoup efficaces.

Un cas typique que j'aime à raconter, autour de 2004 ou 2005 : Un collègue avait écrit un petit script Perl pour faire un reverse-complement d'un fichier de séquences, qui consistait donc un un fichier texte avec des lignes de « ACGT » mais dans lequel chaque ligne pouvait faire éventuellement plusieurs mégas. Le gars avait donc fait un « while (<>) » pour faire une boucle ligne-à-ligne, chacune contenue à tour de rôle dans une variable, puis utilisait « chop » pour extraire le dernier caractère. Il faisait une suite de « if … else » pour en trouver le complément puis ajoutait ce caractère à une autre variable, initialement vide.

Sur le plan sémantique, c'était très propre : « chop » servait exactement à ce que l'on voulait faire, indépendamment du contexte et sans avoir à faire de mesure ou d'indexation préalable. Seulement, en interne, sa chaîne fonctionnait comme en C : le Perl ne semblait pas maintenir la taille de la chaîne en interne et à chaque appel, il fallait reparcourir en entier les méga-octets qui pouvaient la constituer pour trouver le marqueur de fin de chaîne, et pareil pour ajouter le nouveau caractère à la chaîne de sortie, sans compter que tout cela était interprété.

Au bout de deux heures, son script n'avait pas terminé ! Juste pour la curiosité, j'avais ré-écrit le même programme en assembleur, en utilisant le mapping de fichier au niveau du système plutôt qu'une lecture ordinaire, et en utilisant un table d'équivalence pour convertir les bases, même si cela ne concernait que quatre caractères. À l'exécution, le programme a converti le fichier entier en 300 millisecondes.

[Bktero]

Je sais pas si c'est toi qui m'a mis -1 sur mon précédent message, mais vous êtes potentiellement 2 à désapprouver ma réponse. Je vais donc préciser un peu la chose ^^

Il faut distinguer deux cas :

1. L'algorithme est bon, mais tu sais où Perl pêche et comment tu pourrais l'améliorer en C. Dans un cas comme celui-ci, bien sûr, ça vaut la peine d'essayer et ça pourrait être plus rapide !
2. Dans le cas contraire, et notamment sans connaissance suffisamment poussées en C pour savoir comment recoder de manière plus efficace, je vois difficilement cette tentative couronnée de succès.

En fait, ton exemple me corrobore presque ! (oui, j'aime retourner les situations ). Tu savais ce qu'il y avait comme soucis avec le code Perl et tu as réécris ça en..... assembleur ! Pardon, mais ce n'est pas à la portée du premier venu et tu n'es pas parti à l'aveuglette ! Tu cites Python, c'était d'ailleurs un point que je trouvais positif quand j'ai réfléchi à "cela vaut-il la peine d'apprendre Python ?". Tu peux écrire rapidement du code Python et si certaines fonctions sont trop lentes, tu peux les réécrire en C et les intégrer dans le programme Python. Mais cela sous-entend bien que tu sais où sont les problèmes et que tu as une idée de comment les corriger en C.

Globalement, je suis d'accord avec toi : il est probable qu'on puisse diminuer le temps de traitement mais ça dépend (de la machine, des faiblesses de Perl sur un point particulier, et aussi des capacités à recoder ça en C ^^)

[Obsidian]

Bonjour,

Je sais pas si c'est toi qui m'a mis -1 sur mon précédent message

Non, ce n'est pas moi.

mais vous êtes potentiellement 2 à désapprouver ma réponse.

En fait, ça peut être n'importe qui. Il y a beaucoup plus de lecteurs que de postants sur une discussion (115 affichages pour celle-ci). Il arrive aussi que certains se trompent de vote, même si c'est moins fréquent.