Mon article précédent consacré à Rust donne une version de l’algorithme classique de Needleman et Wunsch pour l’alignement de séquences biologiques : chacune des séquences est enregistrée dans un fichier FASTA, le programme lit les deux fichiers et calcule le score de similitude en fonction des coûts de gap et de mutation passés en paramètres, ainsi que l’alignement.

Fort bien, mais ce n’est pas ce que les biologistes font le plus souvent : en général ils ont plutôt une séquence qui les intéresse, dite séquence query (désolé pour l’anglais, mais c’est le terme employé communément), et ils veulent l’aligner avec une collection d’autres séquences, par exemple pour sélectionner les cinq ou les dix qui rendent le meilleur score, et qui sont probablement les plus similaires à la séquence query, selon les paramètres utilisés.

Il faut donc d’une part lire un fichier qui contient la séquence query, comme dans l’article précédent, d’autre part lire un fichier qui contient de multiples séquences, toujours au format FASTA (si les séquences sont dans une banque d’un format différent, il sera loisible de les en extraire au format FASTA en utilisant Biopython, c’est un exercice auquel cette bibliothèque et ce langage excellent et on aurait tort de s’en priver). Voici un exemple d’une telle banque de séquences, ce sont des séquences extraites de génomes d’orchidées.

Donc, lire un fichier qui contient non plus une séquence, mais plusieurs : facile, Rust fournit un itérateur sur les lignes de fichier, on saute les lignes vides, on repère les lignes de commentaire (identifiées par leur premier caractère « > ») pour les garder comme identifiants de la séquence et on recolle les lignes de nucléotides (ou d’acides aminés) les unes aux autres par la méthode suivante :

  1. sequence_nuc.extend(the_line.as_bytes())

Facile ? C’était sans compter avec les règles de possession, de prêt et d’emprunt de Rust, plus quelques petits problèmes de typage. Ainsi, non seulement une fonction perd la possession d’une variable mutable dès qu’elle la passe en argument à une autre fonction (ça j’étais au courant et généralement on peut s’en sortir avec un passage par référence), mais, tout en restant dans la même fonction, si la variable a été utilisée dans une boucle, elle n’est plus disponible au tour suivant de l’itération, et si la variable est un itérateur on a une impression de sorcellerie. J’ai beaucoup ramé, finalement Jmb sur StackOverflow m’a donné une solution simple (après que d’autres généreux contributeurs m’eussent donné des solutions très biscornues, ce qui prouve que je ne suis pas le seul à trouver cela compliqué).

Bref, voilà le code, d’abord main.rs :

  1. extern crate fasta_sequences_read;
  2.  
  3. use fasta_sequences_read::fasta_sequences_read::open_sequences_file;
  4.  
  5. fn main() {
  6.     open_sequences_file();
  7. }

Télécharger

Puis mod.rs :

  1. pub mod fasta_sequences_read;

Et enfin lib.rs, qui contient la substance :

  1. pub mod fasta_sequences_read {
  2.  
  3.     use std::env;
  4.     use std::fs::File;
  5.     use std::io::{BufRead, BufReader};
  6.     use std::str;
  7.     use std::io::Lines;
  8.  
  9.     pub fn open_sequences_file() {
  10.         let args: Vec<String> = env::args().collect();
  11.         let bank_filename = args[1].clone();
  12.         let f_bank = File::open(bank_filename).expect("Fichier non trouvé !");
  13.        
  14.         read_sequences(f_bank);
  15.     }
  16.  
  17.     fn print_seq(sequence: &(String, Vec<u8>)) {
  18.         println!("Ident : {}", sequence.0);
  19.         let sequence_str = str::from_utf8(&sequence.1).unwrap().to_string();
  20.         println!("Séquence : {}", &sequence_str);
  21.     }
  22.  
  23.     fn read_sequences(f: File) {
  24.         let fb = BufReader::new(&f);
  25.         let mut lines = fb.lines();
  26.         let mut count: u8 = 0;
  27.         let mut ident = String::new();
  28.         loop {
  29.             let sequence = get_sequence(&mut count, &mut ident, &mut lines);
  30.             if sequence.1.len() == 0 {
  31.                 break} else {
  32.                 print_seq(&sequence);
  33.             }
  34.         }
  35.     }
  36.  
  37.     fn get_sequence<B: BufRead>(count: &mut u8, ident: &mut String, lines: &mut Lines<B>)
  38.                                    -> (String, Vec<u8>) {
  39.         let mut sequence: (String, Vec<u8>) = (String::new(), vec![]);
  40.         let mut sequence_nuc: Vec<u8> = vec![];
  41.        
  42.         for line in lines {
  43.             let the_line = line.unwrap();
  44.             if the_line.len() > 0 {
  45.                 let first = &the_line[0..1];
  46.                 match first {
  47.                     first if first == ">" => {
  48.                         if *count == 0 {
  49.                             *ident = the_line.clone();
  50.                             *count += 1;
  51.                         } else {
  52.                             sequence = (ident.to_string(), sequence_nuc.clone());
  53.                             println!("Numéro : {}", count);
  54.                             *ident = the_line.clone();
  55.                             sequence_nuc = vec![];
  56.                             *count += 1;
  57.                             return sequence;
  58.                         }
  59.                     }
  60.                     first if first != ">" => {
  61.                         sequence_nuc.extend(the_line.as_bytes())}
  62.                     &_ => {}
  63.                 }
  64.             }
  65.         }
  66.         sequence = (ident.to_string(), sequence_nuc.clone());
  67.         println!("Numéro : {}", count);
  68.         sequence
  69.     }
  70. }

Télécharger

Et l’alignement de séquences, me direz-vous ? À chaque jour suffit sa peine, ce sera pour le prochain article. Cela dit, vous avez tous les éléments, ne reste qu’à les assembler, ce qui devrait vous être possible.