1
\$\begingroup\$

I'm learning Rust, and I have a simple program that I hope to use as a learning exercise. My goal here is to get a better idea of the proper way of doing things.

I'm trying to use the fastq crate to do some manipulation of fastq files - these are simple text files that are standard for storing high throughput [DNA/RNA] sequencing data. What I would like to do is modify a subset of the entries, and write the modified entries to a file.

Right now I have an initial proof of concept below that works to modify the first byte of each of the quality values. I want to determine if I'm doing things sensibly or not. My general thought process:

  1. use fastq::RefRecord::to_owned_record() to create a copy of the Record that I'll have ownership of
  2. Create a bytes::BytesMut from the quality field byte slice
  3. Modify the bytes of the BytesMut
  4. Set the quality field in my OwnedRecord from the modified BytesMut
  5. Verify that it's been set

Does the process above and the code below seem sensible? Are there things you'd do differently, and why? Since I'm so new to Rust, I feel that I'm probably looking at things in a naive and suboptimal way. Let me know your general thoughts.

Dependencies:

  • bytes
  • fastq

Here is my example code:

use bytes::BytesMut;
use fastq::{Parser, parse_path, Record, RefRecord};
use std::env::args;

const READS: &str = r#"@read1/ENST00000266263.10;mate1:84-183;mate2:264-363
GACAGCCAGGGGCCAGCGGGTGGCAGTGCCCAGGACATAGAGAGAGGCAGCACACACGCGGTTGATGGTGAAGCCCGGAATGGCCACAGAGGCTAGAGCC
+
IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII
@read2/ENST00000266263.10;mate1:163-262;mate2:283-382
GATGCCATTGACAAAGGCAAGAAGGCTGGAGAGGTGCCCAGCCCTGAAGCAGGCCGCAGCGCCAGGGTGACTGTGGCTGTGGTGGACACCTTTGTATGGC
+
IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII
@read3/ENST00000266263.10;mate1:86-185;mate2:265-364
GGACAGCCAGGGGCCAGCGGGTGGCAGTGCCCAGGACATAGAGAGAGGCAGCANACACACGGTTGATGGTGAAGCCCGGAATGGCCACAGAGGCTAGAGC
+
IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII!IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII
@read4/ENST00000266263.10;mate1:297-396;mate2:401-500
CAGGAGGAGCTGGGCTTCCCCACTGTTAGGTAGAGCTTGCGCAGGCTGGAGTCCAGGAGGAAATCCACCGACCTGTCAATGGGGTGGATAATGATGGGGA
+
IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII
"#;

fn main() {

    let mut total: usize = 0;

    let parser = Parser::new(READS.as_bytes());
    parser.each(|record| {
        println!("Before mod");
        println!("{}", String::from_utf8_lossy(record.qual()));
        let mut owned_rec = RefRecord::to_owned_record(&record);
        let mut curr_bytes = BytesMut::from(owned_rec.qual());
        curr_bytes[0] = b'$';
        owned_rec.qual = curr_bytes.to_vec();
        println!("After mod");
        println!("{}", String::from_utf8_lossy(owned_rec.qual()));
        total += 1;
        true
    }).expect("Invalid fastq file");
    println!("{}", total);
}

In addition, I'm hoping to move this process into a function, so I'm wondering about the best procedure for doing that. I'm thinking of a function modify_qual that could perform this for me. What would your function signature look like, for something like this? Maybe something that consumes the record and returns the modified record? Something else? I understand the concepts of mutability, passing by reference on the surface level, but don't have a good mental model of when to do certain things.

\$\endgroup\$

1 Answer 1

1
\$\begingroup\$

this feels like a very general question so it is a bit hard to find an answer for. However I have some leads for you. The fastq crate has not been updated for 4 years, so I would advise against relying on that. The bio-rs crate is a very nice replacement, with a lot more functionality besides fastq files. I am unsure why you would need the BytesMut instead of relying on Vec<u8>. Besides that, the code you show does not do anything apart from changing the score of the first Amino Acid which makes it hard to envision what you want from a function. But as a general rule you want to 'hit and forget' whit these (possibly extremely large) fastq files. So that means that you read one record at a time, do you calculations, save whatever information you need in another vector/file and move on to the next. So in general you will not need to return anything from the inner function.

use bio::io::fastq;
use bio::io::fastq::FastqRead;
use bio::io::fastq::Record;

...

fn main() {
    // Note: structure derived from the example given by `bio-rs`.
    let mut total = 0;
    let mut reader = fastq::Reader::new(READS.as_bytes());
    let mut record = fastq::Record::new();

    reader
        .read(&mut record)
        .expect("Failed to parse the first record");
    while !record.is_empty() {
        record
            .check()
            .unwrap_or_else(|_| panic!("I got a rubbish record! Index: {}", total));
        total += 1;

        change_record(&record);

        reader
            .read(&mut record)
            .unwrap_or_else(|_| panic!("Failed to parse record! Index: {}", total));
    }
    println!("Total sequences: {}", total);
}

// Note: the the parameter `record` is passed by reference and nothing is returned.
fn change_record(record: &Record) {
    let mut qual = record.qual().to_vec();
    qual[0] = b'$';
    let new_record = Record::with_attrs(record.id(), record.desc(), record.seq(), &qual);
}
```
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.