Eh toi ! rev_​comp, tu l'écris comment ?

Écrire un algo­rithme de rev_​comp ou com­plé­ment inverse, on l'a tous déjà fait. Aujourd'hui c'est deve­nu un clas­sique des cours d'algorithmique en étude de bio­in­for­ma­tique, c'est un algo­rithme simple mais qui demande de savoir uti­li­ser les struc­tures de contrôle de base. À la fois un bon exer­cice pra­tique et péda­go­gique, doit-il cepen­dant res­ter implé­men­té comme à nos débuts ? Cette ques­tion m'a tra­ver­sé l'esprit durant un cours où mon implé­men­ta­tion et celles de mes cama­rades étaient toutes dif­fé­rentes avec à chaque fois dif­fé­rents algo­rithmes de base. Ne vou­lant en res­ter là, j'ai donc pas­sé mon week-end (et au final les deux mois sui­vants), à étu­dier quels étaient les algo­rithmes les plus per­for­mants par­mi les­quels ceux que je vais vous pré­sen­ter plus en avant.

Le complément inverse, c'est quoi exactement ?

Nous savons que l'ADN est com­po­sé de deux brins anti­pa­ral­lèles (si tu ne le savais pas je suis déso­lé je viens de te dévoi­ler un élé­ment essen­tiel de tes cours de bio­lo­gie). Pour les 3 sta­tis­ti­ciens, l'informaticien et le phy­si­cien du fond on va résu­mer : l'ADN est un poly­mère com­po­sé de 2 chaînes paral­lèles mais dont l'orientation est oppo­sée. Les deux orien­ta­tions sont nom­mées  5'-3' ou 3'-5' (pour en savoir plus c'est par ici). Les deux poly­mères sont com­plé­men­taires : il y aura tou­jours un A en face d'un T, un C en face d'un G, et inver­se­ment.

Cette com­plé­men­ta­ri­té per­met à par­tir d'un seul brin de recons­truire l'autre que ce soit in vivo lors de la répli­ca­tion, ou in sili­co pour retrou­ver les infor­ma­tions pla­cées sur l'autre brin. Un algo­rithme de com­plé­ment prend en entrée une séquence et retourne en sor­tie sa séquence com­plé­men­taire. Par conven­tion, toutes les séquences sont écrites dans le sens 5'-3'. Par consé­quent une fois le brin com­plé­men­taire recons­ti­tué sous la forme 3'-5' l’algorithme doit l'inverser.

Ain­si, dès qu'on ana­lyse une séquence, on s’intéresse aus­si à son com­plé­ment car il peut conte­nir des infor­ma­tions impor­tantes (un gène, une séquence de régu­la­tion, etc…), c'est donc une famille d'algorithmes essen­tielle en bio­in­for­ma­tique.

L'algorithme du jeune bioinformaticien (noob)

Je pense ne pas me trom­per en disant que votre pre­mier algo­rithme de com­plé­ment inverse res­sem­blait sûre­ment à ceci (écrit en pseu­do-code pour ne frois­ser per­sonne) :

Cet algo­rithme est juste, il don­ne­ra le bon résul­tat, mais n'en existe-t-il pas de plus per­for­mant ? Car oui, si vous appe­lez deux fois la fonc­tion dans l'ensemble de votre pro­gramme, perdre un peu de temps ce n'est pas très grave, mais si vous l’appelez plus de dix mille fois comme quand vous tra­vaillez sur des don­nées NGS ça peut rapi­de­ment deve­nir limi­tant.

Comme cha­cun le sait le tra­vail prin­ci­pal d'un algo­rith­mi­cien est de décou­per son pro­blème jusqu’à ce qu'il devienne ridi­cu­le­ment simple, j'ai donc décou­pé le pro­blème en deux : d'un côté le com­plé­ment et de l'autre l'inversion.

Présentation des algorithmes de complément (bien barbu ou pas)

L'algorithme de com­plé­ment va prendre en entrée un nucléo­tide et ren­voyer son com­plé­ment.

L'algorithme naïf (comp_​naif)

On va reprendre la sous-par­tie com­plé­ment de l'algorithme naïf pré­sen­té plus haut.

Rien de bien com­pli­qué, juste une série de si et on retourne la bonne valeur dès que pos­sible, mais atten­tion, si jamais on nous passe une lettre qui n'est pas un nucléo­tide, on ren­ver­ra sys­té­ma­ti­que­ment un T. Il fau­dra pro­ba­ble­ment adap­ter votre algo­rithme en fonc­tion des nucléo­tides pré­sents dans vos jeux de don­nées, il existe de nom­breux cas où votre séquence ne contien­dra pas que des A, des C, des T ou des G, on peut avoir des X, des N, des I, des B, je vous invite à aller consul­ter l'alphabet IUPAC pour l'ensemble des lettres pos­sibles et leur signi­fi­ca­tion.

L'algorithme de correspondance (comp_​hash)

Cet algo­rithme se base sur une table de cor­res­pon­dance (ou table de hachage) : pour chaque nucléo­tide, on asso­cie via la table de cor­res­pon­dance son com­plé­ment. Sans plus attendre voi­ci l'algorithme.

Petite sub­ti­li­té dans l'exemple ci-des­sus, on déclare la table de cor­res­pon­dance à l'extérieur de la fonc­tion pour évi­ter de la recréer à chaque appel. A part cela, il n'y a pas grand-chose de com­pli­qué, cet algo­rithme est assez simple.

Voi­ci cepen­dant une pré­ci­sion sur le fonc­tion­ne­ment d'une table de cor­res­pon­dance :

• la clef est trans­for­mée en valeur par une fonc­tion de hachage
• si la valeur n'existe pas dans la table, créa­tion de la case et ajout de la valeur asso­ciée
• si la valeur existe déjà dans la table, ges­tion de col­li­sion (qui va exi­ger un grand nombre d'opérations)

Je vous invite à lire l'article wiki­pe­dia sur les tables de hachage pour plus de pré­ci­sion. Ce qu'il est impor­tant de rete­nir, c'est qu'une table de cor­res­pon­dance fonc­tionne comme un tableau, mais est un peu plus longue a l’exécution.

L'algorithme du tableau vide (comp_​tab)

Vous savez que dans votre ordi­na­teur les lettres ne sont que des nombres et qu'il existe dif­fé­rents types d'enco­dage des carac­tères. Mais celui qui est le plus sou­vent uti­li­sé à l'heure actuelle (même s'il est de plus en plus rem­pla­cé par d'autres) est l'enco­dage ASCII.

Grâce à la jolie table on découvre que 'T' majus­cule vaut 84, donc si on demande la case T d'un tableau on demande la 85ᵉ case du tableau (rap­pel : Dans les lan­gages infor­ma­tiques sérieux les tableaux com­mencent à l'indice zéro).

En C et en C++ la tra­duc­tion de la lettre en sa valeur ASCII se fait auto­ma­ti­que­ment, dans d'autres lan­gages (notam­ment Python) vous ris­quez de devoir uti­li­ser des fonc­tions spé­ciales. Le tableau créé va occu­per 85 cases mémoire néces­saires pour sto­cker un char, alors qu'on va n'en uti­li­ser que 4 : l'algorithme n'est pas opti­mal dans son usage de la mémoire.

L'algorithme du polynôme de degré trois (comp_​pol3)

Nous savons que chaque lettre est asso­ciée à un nombre, on peut sup­po­ser qu'il existe une fonc­tion mathé­ma­tique qui asso­cie la valeur de 'A' à la valeur de 'T', la valeur de 'T' à la valeur de 'A', etc… . On cherche donc une fonc­tion qui, pour les nombres [65, 67, 71, 84], asso­cie les nombres [84, 71, 67, 65]. Les mathé­ma­ti­ciens ont un super outil pour trou­ver ce type de fonc­tion, ça s'appelle une inter­po­la­tion de Lagrange, je vous épargne les cal­culs mais le résul­tat cor­res­pond à ce poly­nôme d'ordre 3 :

On peut se dire qu'on a juste à mettre tout ça dans le code sans se poser de ques­tion, mais il y a un petit truc cool. Comp­tons le nombre d'opérations à effec­tuer pour trans­for­mer un nucléo­tide en son com­plé­ment : on a 8 mul­ti­pli­ca­tions et 3 addi­tions. Bien évi­de­ment, il existe une solu­tion pour réduire ce nombre d'opérations, et elle nous vient des vieux mathé­ma­ti­ciens : on l'appellera la méthode de Hor­ner (même si cette appel­la­tion est dis­cu­tée). Elle consiste essen­tiel­le­ment à réécrire le poly­nôme sous une autre forme pour réduire le nombre d’opérations.

On a réduit le nombre d’opérations à 3 mul­ti­pli­ca­tions et 3 addi­tions, c'est peut-être un détail pour vous mais pour votre CPU, ça veut dire beau­coup.

On a fini avec les dif­fé­rents algo­rithmes de com­plé­ment. On a donc 4 algo­rithmes. Voi­ci une petite liste réca­pi­tu­la­tive de leurs avan­tages et leurs incon­vé­nients :

• comp_​naif :
  • Avan­tage : simple, pas d'usage mémoire
  • Incon­vé­nients : la com­plexi­té et lon­gueur du code qui aug­mente beau­coup avec le nombre de lettres à gérer
• comp_​hash :
  • Avan­tage : rien de par­ti­cu­lier
  • Incon­vé­nients : usage mémoire non nul, per­for­mance dépen­dante de la fonc­tion de hachage uti­li­sée
• comp_​tab :
  • Avan­tage : une seule opé­ra­tion quel que soit le nombre de lettres à gérer
  • Incon­vé­nients : usage mémoire impor­tant par rap­port aux autres algo­rithmes (faire atten­tion pour l'embarqué)
• comp_​pol3 :
  • Avan­tage : la classe de faire des maths, aucun usage mémoire
  • Incon­vé­nients : nombre d'opérations constant mais néces­site de recal­cu­ler la fonc­tion si l'on veut sup­por­ter de nou­velles lettres, et le poly­nôme sera de degré nombre de lettres — 1

Les algorithmes d'inversion (il ne faut pas avoir Alzheimer)

Pour que vous com­pre­niez bien cette par­tie il va fal­loir une petite expli­ca­tion sur com­ment fonc­tionne votre ordi­na­teur. Il existe plu­sieurs zones où votre ordi­na­teur stocke les variables qu'il uti­lise :

• le cache CPU, une zone mémoire très petite, où l'on stocke les variables qui sont uti­li­sées par les opé­ra­tions en cours (les 2 nombres de l’addition qu'on effec­tue par exemple)
• la mémoire RAM où sont sto­ckées toutes les don­nées sur les­quelles le pro­gramme tra­vaille
• le disque dur où sont sto­ckées toutes les autres don­nées de façon plus pérenne

C'est un résu­mé très gros­sier, mais ce qu'il est impor­tant de noter c'est que tous ces espaces ont des carac­té­ris­tiques dif­fé­rentes en termes de temps d'accès, le plus rapide étant le cache CPU, sui­vi de la mémoire RAM, puis enfin le disque dur (source). On peut aus­si dire que lire une valeur est plus rapide que l'écrire, du moins en ce qui concerne la mémoire RAM (source). Encore une chose avant d'écrire dans la mémoire, il faut deman­der à avoir accès à cette zone mémoire (pour ne pas piquer celle de quelqu'un d'autre), on dit alors qu'on alloue cette zone mémoire à notre pro­gramme, cette opé­ra­tion prend évi­de­ment du temps.

L'algorithme d'inversion fait un peu plus qu'inverser la séquence, il construit la séquence de qui va conte­nir le com­plé­ment. Nous allons donc consi­dé­rer qu'il existe une fonc­tion nuc2nuc(), qui prend en para­mètre un nucléo­tide et ren­voie son com­plé­ment.

Gestion naïve de la mémoire (rev_​naif)

La ges­tion naïve reprend le code du jeune bio­in­for­ma­ti­cien vus pré­cé­dem­ment mais avec quelques amé­lio­ra­tions.

Si on regarde pré­ci­sé­ment ce qui se passe à chaque tour de boucle, on lit un nucléo­tide et on crée une nou­velle chaîne de carac­tère avec la conca­té­na­tion de son com­plé­ment et la chaîne résul­tat. On peut tra­duire en disant qu'on écrit dans la mémoire RAM n chaînes de carac­tères de taille 1 à n (pour n = taille de notre séquence).

Réservation de l'espace mémoire nécessaire à l'avance (rev_​allocate)

Écrire en mémoire plu­sieurs fois la même infor­ma­tion et une très mau­vaise idée, mais si on alloue direc­te­ment l'ensemble de la chaîne de carac­tères dès le début, on écri­ra qu'une seule fois chaque nucléo­tide.

On alloue l'ensemble de la chaîne, et on ins­crit les com­plé­ments depuis la fin.

Pourquoi réserver de la mémoire ? On en a déjà ! (rev_​switch)

C'est une bonne ques­tion. On a déjà une chaîne de carac­tères qui fait la bonne taille, pour­quoi ré-allouer une chaîne de la même taille ? Uti­li­sons celle à dis­po­si­tion.

Cet algo­rithme prend les deux cases extrêmes, cal­cule leur com­plé­ment et les échange, et conti­nue en pro­gres­sant vers le centre, on n'alloue aucune mémoire, donc on gagne du temps. Mais ce gain de temps a un coût. En effet, on modi­fie la chaîne de carac­tères d'origine, ce qui peut éven­tuel­le­ment poser pro­blème.

Allez, une petite liste réca­pi­tu­la­tive des avan­tages et des incon­vé­nients des algo­rithmes d'inversion :

• rev_​naif :
  • Avan­tage : aucun
  • Incon­vé­nients : réa­lise beau­coup d'écritures
• rev_​allocate :
  • Avan­tage : on écrit le reverse com­plé­ment une seule fois
  • Incon­vé­nients : on alloue une deuxième plage mémoire
• rev_​switch
  • Avan­tage : aucune allo­ca­tion
  • Incon­vé­nient : modi­fie la chaîne de carac­tères d'origine

Nous avons vu tous les algo­rithmes, je vais main­te­nant vous pré­sen­ter les résul­tats que j'ai obte­nus, mais tout d'abord un petit point sur le pro­to­cole de test.

Le protocole de test (vachement original comme titre)

Vous trou­ve­rez l'ensemble du code source et des scripts uti­li­sés pour obte­nir les résul­tats pré­sen­tés (et même plus) dans ce dépôt Github.

J'ai implé­men­té l'ensemble des algo­rithmes de com­plé­ment et d'inversion de façon à pou­voir les com­bi­ner, nous avons donc 12 algo­rithmes dif­fé­rents. Ils ont été implé­men­tés dans le lan­gage C++, et peuvent éven­tuel­le­ment être géné­ra­li­sés à d'autres lan­gages mais il faut res­ter pru­dent notam­ment avec les lan­gages inter­pré­tés (Dino­py­thon, Python, Perl, etc.).

Le jeu de test est com­po­sé de séquences d'ADN géné­rées aléa­toi­re­ment, d'une lon­gueur com­prise entre 10 et 10 000 nucléo­tides, avec des GC% com­pris entre 0 et 100 % avec un pas de 5 %.

Mon pro­gramme prend en entrée deux para­mètres : le nombre de fois qu'on va rap­pe­ler chaque fonc­tion, et la séquence à uti­li­ser. J'ai para­mé­tré les scripts d'appel pour que :

En sor­tie, les pro­grammes me four­nissent la somme des temps de cal­cul de chaque fonc­tion, sous forme de tableau CSV.

J'ai appe­lé ce pro­gramme dix fois pour pou­voir avoir une moyenne pour chaque algo­rithme, pour chaque taille de séquence et pour chaque GC%.

Résultat

Les résul­tats pré­sen­tés ont été obte­nus sur un ordi­na­teur por­table avec un Intel core I3 3110M caden­cé a 2.40 GHtz, et pro­pul­sé par une Fedo­ra 22. Ils sont cohé­rents avec des résul­tats obte­nus sur un autre ordi­na­teur lui aus­si d'architecture x86.

Quel est le meilleur algorithme en fonction du GC% ?

On s'attend à ce que le seul algo­rithme impac­té par les varia­tions du GC% soit l'algorithme naïf. En effet, si la séquence n'est com­po­sée que de A et que le pre­mier if test si le nucléo­tide est un A, alors on effec­tue­ra qu'une seule opé­ra­tion, mais si c'est la der­nière, on effec­tue­ra quatre opé­ra­tions.

Mis à part l'algorithme de hachage, tous les autres ne semblent pas impac­tés par le taux de GC, ce qui contre­dit notre hypo­thèse. Peut-être que l’algorithme de hachage uti­li­sé par la biblio­thèque stan­dard C++ pro­voque une col­li­sion entre deux nucléo­tides.

Évolution par rapport à la taille de la séquence

Le taux de GC n'a pas d’impact impor­tant sur les temps de cal­cul de nos algo­rithmes, mais on se doute que la taille des séquences pose­ra pro­blème. Notam­ment pour l'algorithme naïf d'inversion.

La lon­gueur des séquences en abs­cisse suit une échelle loga­rith­mique, les courbes ont l'allure d'une expo­nen­tielle, on a donc une crois­sance qua­dra­tique des temps de cal­cul, ce qui est cohé­rent avec notre réflexion sur l'algorithme d'inversion naïf.

Évolution pour chaque algorithme d'inversion

Est-ce qu'utiliser un algo­rithme d'inversion qui n'est pas naïf aide à réduire l'augmentation du temps de cal­cul ?

La réponse est évi­dem­ment oui, l'augmentation du temps de cal­cul n'est plus qua­dra­tique, ce qui est une très bonne chose. Ce gra­phique ne nous per­met pas de savoir lequel des deux algo­rithmes est le plus per­for­mant.

Tout d'abord un petit point sur l'allure des courbes : avec l'augmentation de la taille des séquences, les temps de cal­cul dimi­nuent, ceci est lié à notre pro­to­cole de test. On lance plus de fois l'algorithme quand la séquence est plus courte, on peut impu­ter ce phé­no­mène aux nom­breux appels de fonc­tions et d'allocations mémoire.

Logi­que­ment, l'algorithme le plus per­for­mant est celui qui ne fait pas d'allocation mémoire.

Conclusion

Pour conclure sur quel est le meilleur algo­rithme de com­plé­ment inverse, voi­ci un petit réca­pi­tu­la­tif de leurs per­for­mances :

• Algo­rithme de com­plé­ment :
  • En termes de temps de cal­cul :
    • pre­mier comp_​tab
    • second comp_​naif
    • troi­sième comp_​pol3
    • qua­trième comp_​hash
  • En termes de consom­ma­tion mémoire
    • pre­mier comp_​naif et comp_​pol3
    • second comp_​hash
    • troi­sième comp_​tab
• Algo­rithme d'inversion en termes de temps
  • pre­mier rev_​switch
  • second rev_​allocate
  • troi­sième rev_​naif

Si l'on veut l'algorithme le plus per­for­mant, en temps de cal­cul, il faut uti­li­ser comp_​tab et rev_​switch, mais si on ne veut pas allouer de mémoire en plus, il faut uti­li­ser comp_​naif et rev_​switch. Bien évi­de­ment ces choix sont à véri­fier et modu­ler en fonc­tion de votre cas par­ti­cu­lier, mais il me semble que ce n'est pas une mau­vaise chose de s'interroger sur quel algo­rithme on uti­lise  pour ce pro­blème qui paraît simple.

J’espère que cet article vous aura inté­res­sé et que vous serez vigi­lant la pro­chaine fois que vous écri­rez un algo­rithme de com­plé­ment inverse. Si jamais vous trou­viez d'autres façons d'écrire votre com­plé­ment inverse (j'en ait per­son­nel­le­ment déjà trou­vé d'autres, mais je vous en par­le­rais dans un autre billet qui sera encore plus tech­nique que celui là), vous pou­vez les pré­sen­ter dans les com­men­taires et/​ou faire une pull request sur le pro­jet git bench_​rev_​comp. Et s'il existe d'autres pro­blèmes clas­siques de la bio­in­for­ma­tique dont vous sou­hai­te­riez connaître les solu­tions, par­lez en aus­si en com­men­taire, je serai ravi de me pen­cher des­sus.

Remerciement :

Je tiens à remer­cier : Pierre Peter­lon­go, Lucas Bour­neuf, Pierre Vignet, les membres de mon mas­ter qui ont sup­por­té mes mails, je remer­cie très sin­cè­re­ment tous les relec­teurs qui m'ont gran­de­ment aidé dans la rédac­tion de cet article : Mathu­rin, m4rsu, Norore, Hed­jour, lelouar, Sacha Schutz, ZaZo0o, Kum­qua­tum.