Site WWW de Laurent Bloch
Slogan du site

ISSN 2271-3905
Cliquez ici si vous voulez visiter mon autre site, orienté vers des sujets moins techniques.

Pour recevoir (au plus une fois par semaine) les nouveautés de ce site, indiquez ici votre adresse électronique :

Avec des expressions régulières :
Recherche de motifs dans des séquences biologiques
Article mis en ligne le 7 octobre 2015
dernière modification le 17 octobre 2015

par Laurent Bloch

Un de mes étudiants du Cnam, Paul Bachelerie, m’a soumis un projet qu’il souhaitait résoudre en Scheme, et dont la solution demandait quelques connaissances non abordées en cours. Le laboratoire où il accomplissait un stage lui avait demandé de le faire en Python, et il voulait comparer les deux langages.

Le projet consistait à rechercher, au moyen d’expressions régulières, des motifs dans la séquence d’un génome bactérien enregistré dans un fichier au format Fasta. C’est plus compliqué que la recherche exacte d’un mot, pour laquelle on peut utiliser l’algorithme de Knuth-Morris-Pratt, mais moins compliqué que la recherche de similitudes par optimisation d’alignement selon l’algorithme de Needleman et Wunsch.

Le format Fasta est très simple : une ligne d’en-tête décrit la séquence, dont le texte suit, en lignes de 70 caractères :

Et voici le génome complet ici :

Le génome d’Escherichia coli (E. coli)

Toute la solution consiste à s’affranchir des caractères de fin de ligne (cf. la procédure multiline, et la documentation qui l’explique) et à lire tout le fichier en mémoire (heureusement c’est un génome procaryote).

La procédure multiline invoque la procédure format, utilisée ici pour insérer après chaque caractère du motif l’expression régulière décrite ainsi : ~(\\n?), qui introduit la possibilité d’un saut de ligne à n’importe quel emplacement du motif. Ainsi :

  1. (define seq "ACCATTGCAATATAGA")
  2.  
  3. (format "~(\\n?)" (string->list seq))
  4.  
  5. -> A\n?C\n?C\n?A\n?T\n?T\n?G\n?C\n?A\n?A\n?T\n?A\n?T\n?A\n?G\n?A

Télécharger

Quant à la procédure read-string utilisée par get-data elle lit jusqu’à épuisement le flux de données ouvert par in-port et renvoie le résultat sous forme d’une chaîne de caractères.

Manuel Serrano, auteur de Bigloo, nous promet pour bientôt une version basée sur mmap, qui fera en gros la même chose mais plus efficacement.

Voilà le code (ci-dessous celui du module principal) :

  1. (module pregexp-search
  2.    (export (pregexp-parse the-pattern in-file)))
  3. ;;
  4. (define (pregexp-parse the-pattern in-file)
  5.    (let ((in-port (open-input-file in-file)))
  6.       (read-line in-port) ;; sauter 1ère ligne fasta
  7.       (analyse the-pattern (get-data in-port))
  8.       (close-input-port in-port)))
  9. ;;
  10. (define (multiline str)
  11.    (format "~(\\n?)" (string->list str)))
  12. ;;
  13. (define (get-data in-port)
  14.    (read-string in-port)) ;; Pour lire TOUT le fichier dans une ligne
  15. ;;
  16. (define (put-data start pos-0 end-pos)
  17.    (print "From " (+ start pos-0) " to " (+ start end-pos)))
  18. ;;
  19. (define (analyse pattern sequence)
  20.    (let ((regexp (pregexp (multiline pattern) 'MULTILINE))
  21.           (len (string-length sequence)))
  22.      (let loop ((start 0))
  23.        (let ((result (pregexp-match-positions regexp sequence start len)))
  24.          (if result
  25.              (begin
  26.                 (put-data start (caar result) (cdar result))
  27.                 (loop (+ start (cdar result)))) ) ) )))

Télécharger

Et le module de pilotage de l’ensemble :

  1. (module benchmark-screen
  2.    (main bench)
  3.    (import pregexp-search))
  4. ;;
  5. (define (bench Args)
  6.    (let ((the-pattern (cadr Args))
  7.          (in-file (caddr Args)))
  8.       (newline)
  9.       (print "Résultats Posix Regexp Multiline")
  10.       (multiple-value-bind (res rtime stime utime)
  11.          (time (lambda () (pregexp-parse the-pattern in-file)))
  12.          (print "real: " rtime "ms sys: " stime "ms user: " utime))
  13.       ))

Télécharger

Ce programme s’utilise ainsi :

./bin/benchmark-screen CAATATAGGCATAGCGCACA sequence_E-coli.fasta

et le résultat est ainsi :

Résultats Posix Regexp Multiline
From 138 to 159
real: 230ms sys: 100ms user: 130

Il faut une version récente du compilateur Bigloo