Skip to main content

histogram/
histogram.rs

1#![allow(unused)]
2
3use helicase::input::*;
4use helicase::*;
5
6const CONFIG: Config = ParserOptions::default()
7    .ignore_headers()
8    .dna_columnar()
9    // don't stop the iterator at the end of a record
10    .return_record(false)
11    .config();
12
13fn main() {
14    let path = std::env::args().nth(1).expect("No input file given");
15
16    // create a parser with the desired options
17    let mut parser = FastxParser::<CONFIG>::from_file(&path).expect("Cannot open file");
18    let mut num_a = 0;
19    let mut num_c = 0;
20    let mut num_t = 0;
21    let mut num_g = 0;
22
23    while let Some(_event) = parser.next() {
24        let cdna = parser.get_dna_columnar();
25        let (high_bits, hi) = cdna.high_bits();
26        let (low_bits, lo) = cdna.low_bits();
27        let rem = cdna.len() % 64;
28        for (&hi, &lo) in high_bits.iter().zip(low_bits) {
29            num_a += (!hi & !lo).count_ones() as usize;
30            num_c += (!hi & lo).count_ones() as usize;
31            num_t += (hi & !lo).count_ones() as usize;
32            num_g += (hi & lo).count_ones() as usize;
33        }
34        if rem > 0 {
35            num_a += (!hi & !lo).count_ones() as usize - (64 - rem);
36            num_c += (!hi & lo).count_ones() as usize;
37            num_t += (hi & !lo).count_ones() as usize;
38            num_g += (hi & lo).count_ones() as usize;
39        }
40    }
41
42    println!("a: {num_a}");
43    println!("c: {num_c}");
44    println!("t: {num_t}");
45    println!("g: {num_g}");
46}