(module nw:lb (main main) (import nw:matrices) (import nw:chains) (import nw:alignment)) (define (nw-2 s1 s2 match-bonus gap-penalty) (let ((ncol (+ (chain-length s1) 2)) (nlin (+ (chain-length s2) 2))) (let ((T (make-matrix nlin ncol 0))) (matrix-margins T s1 s2) ;; (do ((j 2 (+ j 1))) ((= j ncol)) (matrix-set! T 1 j (* (- j 1) gap-penalty))) (do ((i 2 (+ i 1))) ((= i nlin) T) (matrix-set! T i 1 (* (- i 1) gap-penalty)) (do ((j 2 (+ j 1))) ((= j ncol)) (let ((val (max (+ (matrix-ref T (- i 1) (- j 1)) (if (char=? (matrix-ref T i 0) (matrix-ref T 0 j)) match-bonus 0)) (+ (matrix-ref T (- i 1) j) gap-penalty) (+ (matrix-ref T i (- j 1)) gap-penalty)))) (matrix-set! T i j val))))))) ;; on suppose qu'une séquence est dans un fichier fasta ;; read-fasta lit ce fichier et en rend la première séquence (define (read-fasta port) (let ((titre (read-line port))) (if (or (eof-object? titre) (zero? (string-length titre)) (not (char=? (string-ref titre 0) #\>))) (error 'read-fasta "not a fasta file" port) (let loop ((str "")) (let ((lu (read-line port))) (if (or (eof-object? lu) (char=? #\> (string-ref lu 0))) (cons titre str) (begin (print lu) (loop (string-append str lu))))))))) (define (usage) ; pour corriger une erreur d'invocation (print "nw fichier-1 fichier-2 match-bonus gap-penalty") (exit 1)) (define (main argv) (if (not (= (length argv) 5)) (usage) (let ((f1 (cadr argv)) (f2 (caddr argv)) (match-bonus (string->number (cadddr argv))) (gap-penalty (string->number (cadddr (cdr argv))))) (let* ((s1 (make-chain (cdr (let ((port (open-input-file f1))) (read-fasta port))))) (s2 (make-chain (cdr (let ((port (open-input-file f2))) (read-fasta port))))) (the-score-matrix (nw-2 s1 s2 match-bonus gap-penalty))) (matrix-print the-score-matrix) (matrix-print (alignment the-score-matrix match-bonus gap-penalty))))))