Avec Rust : scores de similitude d'une séquence biologique contre une banque

Avec Rust : scores de similitude d’une séquence biologique contre une banque

Article mis en ligne le 25 octobre 2021

dernière modification le 27 octobre 2021

par Laurent Bloch

L’article précédent nous avait permis de lire une banque de séquences biologiques au format FASTA, en introduisant à cette occasion l’instruction match et quelques subtilités de la possession, du prêt et de l’emprunt de valeurs. Celui d’avant nous avait donné l’algorithme de Needleman et Wunsch pour aligner deux séquences et calculer leur score de similitude. Il ne nous reste plus qu’à combiner les deux pour obtenir tous les scores de la séquence d’intérêt contre celles de la banque. Ce travail de construction se fait aisément grâce aux caractéristiques de Rust : on s’est donné du mal pour avoir un programme qui compile, mais une fois que le compilateur a accepté notre code, il n’y a pas de mauvaises surprises à l’exécution.

Je retiens comme score final la valeur de la dernière case de la matrice (« en bas à droite ») : je n’ignore pas que d’autres choix sont loisibles.

Certaines lignes du programme sont placées en commentaire pour améliorer la concision et la lisibilité des résultats, mais il suffit de les « décommenter » pour obtenir le détail des matrices d’alignement.

Voici :

// src/sequences_matrix/build_sequences_matrix.rs :

// This module was inspired by Vincent Esche's Seal crate,
// but simplified and much more basic, without mmap and so on.
// For pedagogic use.

pub mod build_sequences_matrix {

    use simple_matrix::Matrix;
    use std::str;
    use std::char;
    
    pub fn build_matrix(sequence1: &(String, Vec<u8>),
			sequence2: &(String, Vec<u8>),
			match_bonus: f32,
			gap_penalty: f32) {
//			-> Matrix<f32> {
//	print_seq(&sequence1);
//	print_seq(&sequence2);

	let l_seq1: usize = (sequence1.1).len();
	let l_seq2: usize = (sequence2.1).len();

	println!("Longueur première séquence : {} ", l_seq1);
	println!("Longueur seconde séquence : {} ", l_seq2);
	
	let mut the_mat: Matrix::<f32> = Matrix::new(l_seq2+1, l_seq1+1);

	init_matrix(&mut the_mat, l_seq2+1, l_seq1+1, 0.0);

	nw_matrix(&mut the_mat, l_seq2+1, l_seq1+1, match_bonus, gap_penalty, &sequence1.1, &sequence2.1);

	print_ident(&sequence1);
	print_ident(&sequence2);

//	print_matrix(&the_mat, &sequence2.1, l_seq2+1, l_seq1+1);

	print_score(&the_mat, l_seq2+1, l_seq1+1);
	    
    }

    fn nw_matrix(the_mat: &mut Matrix::<f32>,
		     lin: usize,
		     col: usize,
		     match_bonus: f32,
		     gap_penalty: f32,
		     seq1: &Vec<u8>,
		     seq2: &Vec<u8>) {
	for j in 1..col {
	    the_mat.set(0, j, gap_penalty * j as f32) ;
	}
	let mut score: f32 = 0.0;
	for i in 1..lin {
	    the_mat.set(i, 0, gap_penalty * i as f32) ;
	    for j in 1..col {
		if seq1[j-1] == seq2[i-1] {
		    score = match_bonus} else {
		    score = 0.0}
		the_mat.set(i, j, max3(the_mat.get(i-1,j-1).unwrap()
				       + score,
				       the_mat.get(i-1,j).unwrap()
				       + gap_penalty,
				       the_mat.get(i,j-1).unwrap()
				       + gap_penalty));
	    }
	}
    }

    fn max3(v1: f32, v2: f32, v3: f32) -> f32 {
	let tmp = f32::max(v2, v3);
	if v1 > tmp {
	    return v1 } else {
	    return tmp };
    }
    
    fn init_matrix(the_mat: &mut Matrix::<f32>, lin: usize, col: usize, val: f32) {
	for i in 0..lin {
	    for j in 0..col {
		the_mat.set(i, j, val) ;
	    }
	}
    }

// "print_seq" affiche une séquence selon différents formats.
    pub fn print_seq(sequence: &(String, Vec<u8>)) {
		println!("Ident : {:?}", sequence.0);
//		println!("Séquence : {:?}", sequence.1);
		let sequence_str = str::from_utf8(&sequence.1).unwrap().to_string();
		println!("Séquence : {}", &sequence_str);
    }

    fn print_vector(the_vec: &Vec<u8>) {
	let vec_str = str::from_utf8(the_vec).unwrap().to_string();
	print!("{} ", "   ");
	for c in vec_str.chars() {
	    print!("{} ", c);
	}
	print!("{}", "   \n");
    }
    
    fn print_matrix(the_mat: &Matrix::<f32>, seq2: &Vec<u8>, lin: usize, col: usize) {
	for i in 0..lin {
	    if i > 0 {print!("{} ", char::from(seq2[i-1]))} else
	    {print!("{} ", " ")};
	    for j in 0..col {
		print!("{} ", the_mat.get(i, j).unwrap());
	    }
	    print!("{}", "\n")
	}
    }

    fn print_score(the_mat: &Matrix::<f32>, lin: usize, col: usize) {
	println!("Score de similarité : {} ", the_mat.get(lin-1, col-1).unwrap());
	print!("{}", "\n")
    }

    fn print_ident(sequence: &(String, Vec<u8>)) {
	println!("Ident : {:?}", sequence.0);
    }
}

Et le programme de lecture des séquences, adapté de l’article précédent :

// src/fasta_files_mgt/fasta_multiple_cmp.rs :

// https://linuxfr.org/forums/programmationautre/posts/rust-lire-des-donnees-de-type-i8-depuis-un-fichier
// 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/
// https://docs.rs/simple-matrix/0.1.2/simple_matrix/

pub mod fasta_multiple_cmp {

    use std::env;
    use std::fs;
    use std::fs::File;
    use std::io;
    use std::io::Read;
    use std::io::{prelude::*, BufReader};
    use std::io::Lines;
    use std::fs::Metadata;
    use std::str;

    use crate::sequences_matrix::build_sequences_matrix::build_sequences_matrix::print_seq;
    use crate::sequences_matrix::build_sequences_matrix::build_sequences_matrix::build_matrix;
    
    pub struct Config {
	pub query_filename: String,
	pub bank_filename: String,
	pub match_bonus: f32,
	pub gap_penalty: f32
    }
    
    impl Config {
	pub fn new(args: &[String]) -> Config {
	    if args.len() < 5 {
		panic!("pas assez d'arguments");
	    }
	    let query_filename = args[1].clone();
	    let bank_filename = args[2].clone();
	    let match_bonus: f32 = args[3].parse()
		.expect("Ce n'est pas un nombre !");
	    let gap_penalty: f32 = args[4].parse()
		.expect("Ce n'est pas un nombre !");
		
	    Config {query_filename, bank_filename, match_bonus, gap_penalty}
	}
    }
	
    pub fn get_filenames() {
	let args: Vec<String> = env::args().collect();
	let config = Config::new(&args);
	
	println!("Alignement de {} avec {} \n", config.query_filename, config.bank_filename);
	
        let f_query = fasta_open_file(config.query_filename);
        let f_bank = fasta_open_file(config.bank_filename);

	read_sequences(f_query,
		       f_bank,
		       config.match_bonus,
		       config.gap_penalty);

    }

    fn fasta_open_file(filename: String) -> (File) {
        let f = File::open(filename).expect("Fichier non trouvé !");
        f
    }

    fn get_sequence<B: BufRead>(count: &mut u8, ident: &mut String, lines: &mut Lines<B>)
				   -> (String, Vec<u8>) {
	let mut sequence: (String, Vec<u8>) = (String::new(), vec![]);
	let mut sequence_nuc: Vec<u8> = vec![];
	
	for line in lines {
	    let the_line = line.unwrap();
	    if the_line.len() > 0 {
		let first = &the_line[0..1];
		match first {
		    first if first == ">" => {
			if *count == 0 {
			    *ident = the_line.clone();
			    *count += 1;
			} else {
			    sequence = (ident.to_string(), sequence_nuc.clone());
			    println!("Numéro : {}", count);
			    *ident = the_line.clone();
			    sequence_nuc = vec![];
			    *count += 1;
			    return sequence;
			}
		    }
		    first if first != ">" => {
			sequence_nuc.extend(the_line.as_bytes())}
		    &_ => {}
		}
	    }
	}
	sequence = (ident.to_string(), sequence_nuc.clone());
	println!("Numéro : {}", count);
	sequence
    }


//    fn read_sequences(f: File) {
//	let fb = BufReader::new(&f);
//	let mut lines = fb.lines();
//	let mut count: u8 = 0;
//	let mut ident = String::new();
//	loop {
//	    let sequence = get_sequence(&mut count, &mut ident, &mut lines);
//	    if sequence.1.len() == 0 {
//		break} else {
//		print_seq(&sequence);
//	    }
//	}
//  }

    fn read_sequences(f_query: File,
		      f_bank: File,
		      match_bonus: f32,
		      gap_penalty: f32) {
	let fq = BufReader::new(&f_query);
	let mut fq_iter = fq.lines();
	let mut count: u8 = 0;
	let mut ident = String::new();
	let query_sequence = get_sequence(&mut count, &mut ident, &mut fq_iter);
	print_seq(&query_sequence);

	let fb = BufReader::new(&f_bank);
	let mut fb_iter = fb.lines();
	count = 0;
	loop {
	    let bank_sequence = get_sequence(&mut count, &mut ident, &mut fb_iter);
	    if bank_sequence.1.len() == 0 {
		break} else {
//		print_seq(&bank_sequence);
		build_matrix(&query_sequence,
			     &bank_sequence,
			     match_bonus,
			     gap_penalty);
	    }
	}
    }
}

Dans la même rubrique

Votre inscription a été enregistrée avec succès !

Au menu

Dans la même rubrique