Jeux Tecleux Novembre 2009
mou...@igbmc.u-strasbg.fr
Le Jeux Tecleux de ce mois-ci porte encore sur la génomique ;-)
Je soumets à votre sagacité l'écriture d'un code optimisé concernant
la comparaison entre 2 séquences de protéines:
Seq1 ACD...FFFI.M
| | ||
Seq2 A.....FYFI.K
(note: les points ne sont pas des points de suspension)
Il s'agit de déterminer le pourcentage d'identité qu'il y a entre des
séquences de protéine.
Pour cela on compte le nombre de lettres identiques qui se font face
pour 2 sequences donnees. Dans notre exemple ce nombre vaut 4.
Et on divise par le nombre de lettres de la plus petite des deux
sequences. Dans notre exemple le nombre de lettres vaut 6 (Seq2).
Donc le pourcentage vaut p = 4/6 = 0.6666
Ca parait simple, MAIS :
- les sequences font entre 1000 et 10 000 caracteres de long
- on a en general entre 300 et 500 sequences, il faut le pourcentage
d'identite de tous les couples, soit par exemple, (300*299/2) paires
de sequences ...
L'aspect vitesse devient crucial !
Pour amorcer, je fournis une proc qui genere une liste de sequence
(aleatoire ici) qui peut s'apparenter a un cas reel ....
Que le meilleur gagne !
Luc Moulinier
(Merci a GS pour le coup de main !)
proc Data {} {
set ll [list A C D E F G H I K L M N P Q R S T V W
Y . . . . . . . . . . . . . . . . . . . . .]
for {set n 0} {$n < 100} {incr n} {
set seq ""
for {set i 0} {$i < 1000} {incr i} {
append seq [lindex $ll [expr {round(40*rand())}]]
}
lappend Lseq $seq
}
return $Lseq
}
Kroc
proc Compare {a b} {
set A [split $a {}]
set B [split $b {}]
if {[llength $A] > [llength $B]} {
set max [llength $A]
} else {
set max [llength $B]
}
set nA 0
set nB 0
set ni 0
for {set i 0} {$i < $max} {incr i} {
if {[string is upper -strict [lindex $A $i]]} {
incr nA
}
if {[string is upper -strict [lindex $B $i]]} {
incr nB
}
if {[lindex $A $i] eq [lindex $B $i]} {
incr ni
}
}
if {$nA > $nB} {
return "$max éléments, $ni identités / $nB\
séquences = [expr {1.0 * $ni / $nB}]"
} else {
return "$max éléments, $ni identités / $nA\
séquences = [expr {1.0 * $ni / $nA}]"
}
}
Chez moi c'est assez rapide :
% set A [Data] ; set B [Data] ; Compare $A $B
100099 éléments, 27488 identités / 48666 séquences =
0.5648296552007561
% time {Compare $A $B} 100
199254.1 microseconds per iteration
--
David Zolli
Kroc
proc Compare {a b} {
set A [split $a {}]
set B [split $b {}]
if {[llength $A] > [llength $B]} {
set max [llength $A]
} else {
set max [llength $B]
}
set nA 0
set nB 0
set ni 0
for {set i 0} {$i < $max} {incr i} {
if {[string is upper -strict [lindex $A $i]]} {
incr nA
}
if {[string is upper -strict [lindex $B $i]]} {
incr nB
}
if {[lindex $A $i] ne "." && [lindex $A $i] eq [lindex $B $i]} {
incr ni
}
}
if {$nA > $nB} {
return "$max éléments, $ni identités / $nB\
séquences = [expr {1.0 * $ni / $nB}]"
} else {
return "$max éléments, $ni identités / $nA\
séquences = [expr {1.0 * $ni / $nA}]"
}
}
Du coup, c'est un poil plus lent :
% Compare $A $B
100099 éléments, 1323 identités / 48666 séquences =
0.027185303908272715
34% time {Compare $A $B} 100
216356.32 microseconds per iteration
--
David Zolli
Kroc
proc Compare {a b} {
set max 0
set nA 0
set nB 0
set ni 0
foreach A [split $a {}] B [split $b {}] {
incr max
if {[string is upper -strict $A]} {
incr nA
}
if {[string is upper -strict $B]} {
incr nB
}
if {$A ne "." && $A eq $B} {
incr ni
}
}
if {$nA > $nB} {
return "$max éléments, $ni identités / $nB\
séquences = [expr {1.0 * $ni / $nB}]"
} else {
return "$max éléments, $ni identités / $nA\
séquences = [expr {1.0 * $ni / $nA}]"
}
}
% Compare $A $B
100099 éléments, 1323 identités / 48666 séquences =
0.027185303908272715
% time {Compare $A $B} 100
220807.61 microseconds per iteration
--
David Zolli
Eric Hassold
> Voil� ma proposition :
> ...
> Chez moi c'est assez rapide :
>
> % set A [Data] ; set B [Data] ; Compare $A $B
> 0.5648296552007561
>
> % time {Compare $A $B} 100
> 199254.1 microseconds per iteration
Oui, mais 200ms fois (n*(n-1)/2), avec n=100, ca fait tout de suite 15
minutes. Et si, comme Luc le decrit, le probleme typique est plutot avec
n=500, et des sequences de plusieurs milliers de symboles, alors la, ca
commence a faire vraiment long. D'ou... meme si ca fait mal de le dire :
quand "la vitesse compte", alors il faut se resoudre a abandonner Tcl
et basculer en langage compile, au moins pour le goulot d'etranglement
de tout le calcul.
Hors sujet (donc hors-jeux) puisque pas en Tcl au sens strict, mais
juste pour se faire une idee, la meme solution (avec le meme algorithme
elementaire), mais implemente avec Odyce (ou critcl). Teste avec eTcl,
mais un interprete Tcl autre, avec critcl disponible devrait marcher aussi.
Le gain est vraiment ... enorme (vitesse d'execution de 40x a 100x plus
rapide pour la version Odyce que pour la version pure Tcl). Dans le cas
"extreme" de 500 sequences de 10000 caracteres, on est de l'ordre d'une
vingtaine de secondes sur un PC tres standard (100 fois plus, ca
commence a faire long).
Eric
-----
Eric Hassold
Evolane - http://www.evolane.com/
-----
# Teste en utilisant l'emulation critcl dans Odyce
package require critcl
::critcl::ccode {
static double
SeqMatch(char *seq1, char *seq2, int l)
{
int l1, l2;
int r;
int i;
double ratio;
l1=0;
l2=0;
r=0;
for (i=0;i<l;i++) {
if (seq1[i]!='.') l1++;
if (seq2[i]!='.') l2++;
if (seq1[i]!='.' && seq1[i]==seq2[i]) {
r++;
}
}
l=((l1<l2) ? l1 : l2);
if (l<=0) {
ratio=0.0;
} else {
ratio = ((double) r)/l;
}
return ratio;
}
}
# Une fonction Tcl pour comparer deux sequences
::critcl::ccommand GMatchC {dummy interp objc objv} {
int n;
Tcl_Obj *resultPtr;
char *seq1, *seq2;
int l1,l2, l;
int r;
int i;
double ratio;
if (objc!=3) {
Tcl_WrongNumArgs(interp,1,objv,"seq1 seq2");
return TCL_ERROR;
}
/* Use Byte array instead of strings */
seq1=Tcl_GetByteArrayFromObj(objv[1],&l1);
seq2=Tcl_GetByteArrayFromObj(objv[2],&l2);
if (l1!=l2 || l1<=0) {
return TCL_ERROR;
}
Tcl_SetObjResult(interp,Tcl_NewDoubleObj(SeqMatch(seq1,seq2,l1)));
return TCL_OK;
}
# Une fonction Tcl pour comparer toutes les paires de sequences
::critcl::ccommand GMatchAllC {dummy interp objc objv} {
int n;
Tcl_Obj *resultPtr;
char *seq1, *seq2;
int l1,l2, l;
int r;
int i,j;
int nbseqs;
Tcl_Obj *matrix;
Tcl_Obj *row;
double ratio;
int pobjc;
Tcl_Obj **pobjv;
if (objc!=2) {
Tcl_WrongNumArgs(interp,1,objv,"lsequences");
return TCL_ERROR;
}
if (Tcl_ListObjGetElements(interp,objv[1],&pobjc,&pobjv)!=TCL_OK) {
return TCL_ERROR;
}
nbseqs=pobjc;
matrix = Tcl_NewListObj(0,NULL);
for (i=0;i<nbseqs-1;i++) {
seq1=Tcl_GetByteArrayFromObj(pobjv[i],&l1);
row = Tcl_NewListObj(0,NULL);
for (j=i+1;j<nbseqs;j++) {
seq2=Tcl_GetByteArrayFromObj(pobjv[j],&l2);
if (l1!=l2 || l1<=0) {
ratio=0;
} else {
ratio=SeqMatch(seq1,seq2,l1);
}
Tcl_ListObjAppendElement(interp,row,Tcl_NewDoubleObj(ratio));
}
Tcl_ListObjAppendElement(interp,matrix,row);
}
Tcl_SetObjResult(interp,matrix);
return TCL_OK;
}
# Generation des donnees de test
proc GTest {{nbseqs 100} {seqlen 1000}} {
set ll [list A C D E F G H I K L M N P Q R S T V W Y . . . . . . . .
. . . . . . . . . . . . .]
for {set n 0} {$n < $nbseqs} {incr n} {
set seq ""
for {set i 0} {$i < $seqlen} {incr i} {
append seq [lindex $ll [expr {round(40*rand())}]]
}
lappend Lseq $seq
}
return $Lseq
}
# Calcul de la "longueur" d'une sequence
proc GLength {s} {
return [string length [string map [list "." ""] $s]]
}
proc GMatchPost {seq1 len1 seq2 len2} {
set l [expr {($len1<$len2) ? $len1 : $len2}]
if {$l<=0} {
return 0
}
set r 0
foreach c1 [split $seq1 ""] c2 [split $seq2 ""] {
if {$c1 eq $c2 && $c1 ne "."} {
incr r
}
}
return [expr {($r*1.0)/$l}]
}
proc GMatch {seq1 seq2} {
GMatchPost $seq1 [GLength $seq1] $seq2 [GLength $seq2]
}
proc GMatchAll {lseqs} {
# Pre-compute length of each sequence
set l [list]
foreach seq $lseqs {
lappend l [GLength $seq]
}
# Find matching sequences for each pair
set nbseqs [llength $lseqs]
set c [list]
for {set i 0} {$i<$nbseqs-1} {incr i} {
set s1 [lindex $lseqs $i]
set l1 [lindex $l $i]
for {set j [expr {$i+1}]} {$j<$nbseqs} {incr j} {
set s2 [lindex $lseqs $j]
set l2 [lindex $l $j]
lappend c [GMatchPost $s1 $l1 $s2 $l2]
}
}
return $c
}
puts "Preparing test data"
# set lseqs [GTest 400 5000]
set lseqs [GTest 100 1000]
# Pure Tcl version
set t1 [time {GMatchAll $lseqs}]
set ms1 [expr {[lindex $t1 0]/1000}]
puts "PureTcl: $ms1 milliseconds"
# C version
set t2 [time {GMatchAllC $lseqs}]
set ms2 [expr {[lindex $t2 0]/1000}]
puts "Odyce: $ms2 milliseconds"
if {$ms2<$ms1} {
puts "Odyce version is [expr {$ms1/$ms2}] times faster than Pure Tcl"
}
Kroc
proc Compare {sa sb} {
set max [llength [set A [split $sa {}]]]
set lb [llength [set B [split $sb {}]]]
if {$lb > $max} {set max $lb}
for {set i 0} {$i < $max} {incr i} {
if {[set a [lindex $A $i]] ne "."} {incr nA}
if {[set b [lindex $B $i]] ne "."} {incr nB}
if {$a ne "." && $a eq $b} {incr ni}
}
return [expr {$nA>$nB?1.0*$ni/$nB:1.0*$ni/$nA}]
}
% time {Compare $A $B} 100
72978.15 microseconds per iteration
--
David Zolli
Kroc
> % time {Compare $A $B} 100
> 72978.15 microseconds per iteration
Je ne suis pas peu fier : ma version Tcl est plus rapide que celle
d'Éric :
% time {GMatchAll [list $A $B]} 100
73419.53 microseconds per iteration
Notez l'écart monstrueux ! ;^)
--
David Zolli
mou...@igbmc.u-strasbg.fr
Je suis moi-mpeme pas mécontent de constater que je n'étais pas loin
non plus de la solution de David, à quelques broutilles près.
J'en profite pour vous poser quelques questions :
- il semble selon vos scripts, que vous favorisiez "for {...} " a
"foreach", car for serait plus rapide. ou peut-on trouver ce genre
d'info ?
- Si j'utilise critcl, y a t il un moyen simple de savoir si gcc est
installe sur la becane ? Cette appli est destine aux biologistes, et
il est hors de question de leur demander d'installer gcc/mingw ou
autre ...
Merci encore !
Kroc
wrote:
> Merci les goutous !! je n'en espérais pas tant !
> Je suis moi-mpeme pas mécontent de constater que je n'étais pas loin
> non plus de la solution de David, à quelques broutilles près.
Bien Luc ! La Force est avec toi ! ;^)
> J'en profite pour vous poser quelques questions :
> - il semble selon vos scripts, que vous favorisiez "for {...} " a
> "foreach", car for serait plus rapide. ou peut-on trouver ce genre
> d'info ?
Des infos piquées aux gens du TCT lorsqu'on bavarde avec eux, et
beaucoup de tests.
> - Si j'utilise critcl, y a t il un moyen simple de savoir si gcc est
> installe sur la becane ? Cette appli est destine aux biologistes, et
> il est hors de question de leur demander d'installer gcc/mingw ou
> autre ...
Le principe génial d'Odyce c'est que le compilateur (basé sur TinyCC)
est intégré à l'extension Tcl. Du coup, ça marche partout et tout le
temps ! C'est chouette, non ?
--
David Zolli