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 :

Configurer des séquences biologiques pour les aligner
Article mis en ligne le 3 juillet 2021
dernière modification le 17 août 2021

par Laurent Bloch

Pour programmer un algorithme génial qui méritera le prix Turing, il faudra souvent avoir au préalable acquis les données, généralement depuis un ou plusieurs fichiers, et les avoir mises sous la forme appropriée. C’est un travail peu glorieux, mais ce n’est pas forcément très facile, cela prend du temps et des lignes de code, et oblige à se faire des idées sur le système d’exploitation. C’est l’objet du présent article.

Codage, caractères, glyphes

Cet article vient à la suite de Premiers programmes en Rust, de Programmation Rust, suite et de Codage de séquences biologiques avec Rust. Je ne suis pas encore arrivé au cœur du sujet, mais au seuil : les programmes de cet épisode peuvent lire un fichier au format FASTA qui contient une (unique) séquence, vérifier que la première ligne est bien conforme au format, l’extraire pour éventuellement documenter la séquence, ensuite retirer du texte proprement dit de la séquence les caractères de saut de ligne (LF, Line Feed) afin d’obtenir une séquence propre, prête à subir les exactions de Messieurs Needleman et Wunsch, ou, si l’on préfère, de Messieurs Smith et Waterman, ce que je prévois de faire dans un prochain épisode.

Pour suivre les conseils de lecteurs compétents, j’ai renoncé au type chaîne de caractères : en effet les caractères de Rust obéissent à la norme UTF-8, qui leur permet d’occuper jusqu’à quatre octets, ce qui est encombrant, mais surtout cette occupation est de taille variable, un caractère peut occuper un, deux, trois ou quatre octets, selon sa signification [1]. Pour éviter ces caractères de taille variable, j’ai décidé d’utiliser pour coder les nucléotides le type unsigned 8, qui correspond à un octet. Je sais qu’il est possible d’utiliser des codages encore plus condensés, mais après tout les algorithmes envisagés valent aussi pour les acides aminés, alors vivent les octets.

Martin Larralde m’a conseillé « de créer un wrapper de Vec, et d’implémenter [moi-même] les traits Display et Debug, de manière à pouvoir visualiser la séquence sous forme de texte (ce qui se fait de façon triviale avec std::str::from_utf8_unchecked , si on respecte l’invariant que le Vec ne contient que des caractères alphanumériques). » Bon, j’ai bien utilisé std::str::from_utf8_unchecked, mais le wrapper attendra que je me sois un peu perfectionné en Rust.

Typage et dépendances

La grande originalité de Rust, et sa qualité exclusive, est son modèle de mémoire, qui garantit à la compilation l’absence de débordement de buffer ou de toute autre zone de mémoire, et aussi le fait qu’un seul sous-programme puisse, à un instant donné, modifier la valeur d’une variable : toute valeur appartient à une variable et une seule, qui est limitée à une portée donnée par la syntaxe. Si un sous-programme passe une de ses variables à un autre sous-programme, ce dernier en devient le propriétaire, et le premier ne peut plus modifier sa valeur.

Ces caractéristiques font la sûreté du langage, mais elles ont un prix. Pour la cinématique des données, il faut y réfléchir à deux fois avant de passer une variable en argument à un autre sous-programme, parce qu’après on ne pourra plus y toucher ; on peut souvent s’en tirer avec des références, comme en Ada. C prétend avoir un typage fort : cette assertion est spécieuse, parce que si le typage est fort, il peut facilement être contourné pour provoquer les failles de sécurité qui font la joie des pirates. Rust ne plaisante pas avec le typage, on n’y coupe pas. Certaines méthodes de Vec<u8> (données sur le tas) renvoient des tranches (slices) de type array [u8] (données sur la pile), il faut leur appliquer la méthode to_vec() (conversion d’array en Vec) pour les utiliser, ainsi :

  1. ident = str::from_utf8(fasta_sequence.split_at(index).0).unwrap().to_string();
  2. sequence_nuc = (fasta_sequence.split_at(index).1).to_vec();

Télécharger

La méthode fasta_sequence.split_at(index) renvoie un tuple (type d’objet commode et indulgent que Rust fournit au programmeur pour le consoler de ses souffrances) de deux éléments de type array [u8] dont le premier élément contient les premiers éléments du vecteur fasta_sequence jusqu’à l’indice index (non compris), et le second les éléments depuis index (compris) jusqu’à fasta_sequence.len() (la longueur totale du vecteur), non compris. Finalement c’est bien conçu, mais ces deux lignes m’ont fait transpirer avant que j’aie compris comment cela marchait.

On remarque que fasta_sequence.split_at(index).0 désigne le premier élement du tuple, et fasta_sequence.split_at(index).1 le second.

Incidemment, cette notation des méthodes, empruntée à la programmation par objets, est bien commode, mais Rust n’est pas un langage à objets, avec toutes ces histoires d’héritage multiple et de surcharge d’opérations qui ont finalement introduit plus de confusion qu’autre chose.

Le programme

Voici le corps du module fasta_files_mgt. Reste à écrire le module sequences_matrix, qui soumettra nos deux séquences à l’algorithme de Messieurs Needleman et Wunsch, ou à celui, cousin, de Messieurs Smith et Waterman. Les autres fichiers, inchangés, sont ceux de l’article Codage de séquences biologiques avec Rust, où l’on trouvera aussi des séquences au format FASTA pour essayer.

Les pages Web dont je me suis inspiré :

Lire un fichier d’octets

Lire et décrire un fichier

  1. // src/fasta_files_mgt/fasta_open_read.rs :
  2.  
  3. // https://linuxfr.org/forums/programmationautre/posts/rust-lire-des-donnees-de-type-i8-depuis-un-fichier
  4. // https://www.it-swarm-fr.com/fr/file-io/quelle-est-la-maniere-de-facto-de-lire-et-decrire-des-fichiers-dans-rust-1.x/1054845808/
  5. // https://docs.rs/simple-matrix/0.1.2/simple_matrix/
  6.  
  7. pub mod fasta_open_read {
  8.  
  9.     use std::env;
  10.     use std::fs;
  11.     use std::fs::File;
  12.     use std::io::Read;
  13.     use std::fs::Metadata;
  14.     use std::str;
  15.  
  16.     use crate::sequences_matrix::build_sequences_matrix::build_sequences_matrix::build_matrix;
  17.    
  18.     pub struct Config {
  19.         pub filename1: String,
  20.         pub filename2: String,
  21.         pub match_bonus: f32,
  22.         pub gap_penalty: f32
  23.     }
  24.    
  25.     impl Config {
  26.         pub fn new(args: &[String]) -> Config {
  27.             if args.len() < 5 {
  28.                 panic!("pas assez d'arguments");
  29.             }
  30.             let filename1 = args[1].clone();
  31.             let filename2 = args[2].clone();
  32.             let match_bonus: f32 = args[3].parse()
  33.                 .expect("Ce n'est pas un nombre !");
  34.             let gap_penalty: f32 = args[4].parse()
  35.                 .expect("Ce n'est pas un nombre !");
  36.                
  37.             Config {filename1, filename2, match_bonus, gap_penalty}
  38.         }
  39.     }
  40.        
  41.     pub fn get_filenames() {
  42.         let args: Vec<String> = env::args().collect();
  43.         let config = Config::new(&args);
  44.        
  45.         println!("Alignement de {} avec {} \n", config.filename1, config.filename2);
  46.        
  47.         let mut data1 = fasta_open_file(config.filename1);
  48.         let mut data2 = fasta_open_file(config.filename2);
  49.  
  50.         let sequences = read_seq(&mut data1, &mut data2);
  51.         build_matrix(sequences.0,
  52.                      sequences.1,
  53.                      config.match_bonus,
  54.                      config.gap_penalty);
  55.     }
  56.  
  57.     fn fasta_open_file(filename: String) -> (File, Metadata) {
  58.         let metadata = fs::metadata(&filename).expect("Fichier non trouvé.");
  59.         let f = File::open(filename).expect("Fichier non trouvé !");
  60.         (f, metadata)
  61.     }
  62.  
  63. // "read_seq" lit chacun des deux fichiers au format FASTA et ses métadonnées.
  64.     fn read_seq(data1: &mut (File, Metadata),
  65.                 data2: &mut (File, Metadata))
  66.                 -> ((String, Vec<u8>), (String, Vec<u8>)) {
  67.         let mut seq1 = vec![0; data1.1.len() as usize];
  68.         let mut seq2 = vec![0; data2.1.len() as usize];
  69.         data1.0.read(&mut seq1).expect("débordement de buffer");
  70.         data2.0.read(&mut seq2).expect("débordement de buffer");
  71.        
  72.         let seq1_str = str::from_utf8(&seq1).unwrap().to_string();
  73.         let seq2_str = str::from_utf8(&seq2).unwrap().to_string();
  74.        
  75.         println!("1er caractère : {:?}", &seq1[0]);
  76.         println!("1er caractère : {:?}", &seq1_str.bytes().nth(0));
  77.         let sequence1 = split_seq(seq1);
  78.         let sequence2 = split_seq(seq2);
  79.  
  80.         (sequence1, sequence2)
  81.     }
  82.  
  83. // D'un fichier au format FASTA, extraire d'une part la première ligne, de commentaire,
  84. // identifiée par le caractère ">" en première position, dont le texte documentera la
  85. // séquence sous le nom "ident", d'autre part les lignes suivantes, qui contiennent la
  86. // séquence proprement dite des nucléotides (ou des acides aminés), sous le nom
  87. // "sequence_nuc". Après traitement de "sequence_nuc" par la fonction "remove_LF",
  88. // la fonction "split_seq" renvoie le tuple "(ident, sequence_clean)".
  89.     fn split_seq(fasta_sequence: Vec<u8>) -> (String, Vec<u8>) {
  90.             let mut ident: String = "abcd".to_string();
  91.             let mut sequence_nuc: Vec<u8> = vec![3, 4, 5];
  92.             if fasta_sequence[0] == 62 {
  93.                 let index = fasta_sequence.iter().position(|x| *x == 10).unwrap();
  94.                 println!("Position du premier LF : {}", index);
  95.                 ident = str::from_utf8(fasta_sequence.split_at(index).0).unwrap().to_string();
  96.                 sequence_nuc = (fasta_sequence.split_at(index).1).to_vec();
  97.             }
  98.         let sequence_clean: Vec<u8> = remove_LF(&sequence_nuc);
  99.         (ident, sequence_clean)
  100.     }
  101.  
  102. // Cette fonction "remove_LF" reçoit les lignes de nucléotides ou d'acides aminés
  103. // et en retire les caractères LF (10 selon le code Ascii) pour former une seule ligne.
  104.     fn remove_LF(sequence_nuc: &Vec<u8>) -> Vec<u8> {
  105.         let mut sequence_clean: Vec<u8> = (&sequence_nuc).to_vec();
  106.         for i in (0..sequence_clean.len()).rev() {
  107.             if sequence_clean[i] == 10 {
  108.                 sequence_clean.remove(i);
  109.             }
  110.         }
  111.         sequence_clean
  112.     }
  113. }

Télécharger