I wanted to write some code which reads in a FASTA file and one hot encodes the sequence which is consequentially saved to a file. A FASTA file is a text based file format commonly used in bioinformatics. Here's an example:
>chr1 ACGTCTCGTGC NNNNNNNNNNN ACGTTGGGGGG >chr2 AGGGTCTGGGT NNNNNNNNNNN It has a header string with > and then some genetic sequence following. Every character in that sequence represents a nucleotide. For example "A" stands for Adenine and "N" for unknown. The number of possible characters in the sequence is limited, but there are some encodings allowing different representations. But this is not that important for me, since I will probably filter beforehand anyway.
The output should be a text file with the one hot encoded sequence and the header for identification. Depending on how one would choose the encoding this could look something like this for the first three nucleotides(ACG):
>chr1 0,0,0,1 0,1,0,0 0,0,1,0 I wrote some Rust code together to achieve this and it works, but it seems to be pretty slow. Maybe I should also mention at this point that the FASTA files can get really big (many GB) in no time. I tried my code (cargo build --release) with an ~60MB .fa file and took about ~1m.30s and some generic Python code I wrote together in 10 minutes took only ~40s.
Now I wanted to know why my Rust Code is slower then it's Python equivalent. Did I write it inn an idiomatic Rust way? I do not want to store the data in something like HDF5. I just want to understand Rust better. Especially I/O heavy text processing tasks, since these are common in bioinformatics.
Here's the code I wrote:
use std::fs::File; use std::io::{BufReader, Write}; use std::io::prelude::*; use std::env; use std::fs::OpenOptions; fn main() -> std::io::Result<()>{ // gather cmd arguments let args: Vec<String> = env::args().collect(); let filename: &String = &args[1]; let outfile: &String = &args[2]; // Open the file to read from let file = File::open(&filename).unwrap(); // Create a buffered reader -> instead of reading everything to ram let reader = BufReader::new(file); let lines = reader.lines(); // this allows to open a file in file append mode to write (or create it) let mut to_write_file = OpenOptions::new() .append(true) .create(true) // Optionally create .open(outfile) .expect("Unable to open file"); // iterate over every line for l in lines { // check if we get a string back (Result<String, Error>) if let Ok(line_string) = l { // check for header lines if line_string.contains("chr"){ println!("Processing: {}",line_string); // write to file so we can keep track of chroms writeln!(to_write_file, "{}",line_string).unwrap(); } else { // iterate over nucleotides for character in line_string.chars() { // one hot encode the nucleotide let nucleotide: String = one_hot_encode_string(String::from(character)); writeln!(to_write_file,"{}", nucleotide)?; } } } } Ok(()) } fn one_hot_encode_string(nucleotide: String) -> String { // One hot encodes a Nucleotide: String let encoding: String; // create return variable // OHE the nucleotide match nucleotide.as_str() { "A" => encoding = String::from("0,0,0,1"), "G" => encoding = String::from("0,0,1,0"), "C" => encoding = String::from("0,1,0,0"), "T" => encoding = String::from("1,0,0,0"), "N" => encoding = String::from("0,0,0,0"), _ => encoding = String::from("NA"), } encoding } Her's also the "naive" python code:
import sys # one hot encode nulceotide def ohe(nucleotide: str) -> str: if nucleotide == "A": return("0,0,0,1") if nucleotide == "G": return("0,0,1,0") if nucleotide == "T": return("0,1,0,0") if nucleotide == "C": return("1,0,0,0") if nucleotide == "N": return("0,0,0,0") else: return("NA") def main() -> None: # get arguments input: str = sys.argv[1] output: str = sys.argv[2] # open in file with open(input) as in_f: # open out file with open(output, 'a') as out_f: for line in in_f: # get header if line.startswith(">"): out_f.write(f"{line.rstrip()}\n") else: # iterate over every nucleotide for character in line.rstrip(): # one hot encode character # print(ohe(character)) out_f.write(f"{ohe(character)}\n") if __name__ == "__main__": main() 
