debruijn/
lib.rs

1// Copyright 2017 10x Genomics
2
3//! # debruijn: a De Bruijn graph library for DNA seqeunces in Rust.
4//! This library provides tools for efficient construction DeBruijn graphs (dBG)
5//! from DNA sequences, tracking arbitrary metadata associated with kmers in the
6//! graph, and performing path-compression of unbranched graph paths to improve
7//! speed and reduce memory consumption.
8//!
9//! Most applications of `debruijn` will follow this general workflow:
10//! 1. You generate a set of sequences to make a dBG from.
11//! 2. You pass those sequences to the `filter_kmers` function, which converts the sequences into kmers, while tracking 'metadata' about each kmer in a very customizable way. The metadata could be read count, a set of colors, a set of read counts split by haplotype, a UMI count, etc.
12//! 3. The the library will convert the kmers to a compressed dBG. You can also customize the rules for how to compress the dBG and how to 'combine' the per-kmer metadata.
13//!
14//! Then you can use the final compressed dBG how you like. There are some methods for simplifying and re-building the  graph, but those could be developed more.
15//!
16//! ## Examples
17//! - [Local phased SV assembly tool in our Long Ranger package](https://github.com/10XGenomics/longranger/blob/master/lib/pvc/src/asm_caller.rs#L205)
18//! - [Single-cell VDJ assember](https://github.com/10XGenomics/cellranger/blob/master/lib/rust/vdj_asm/src/asm.rs#L191)
19//! - [Build a colored, compressed dBG of a transcriptome reference](https://github.com/10XGenomics/rust-pseudoaligner/blob/master/src/build_index.rs#L40)
20//!
21//! All the data structures in debruijn-rs are specialized to the 4 base DNA alphabet,
22//! and use 2-bit packed encoding of base-pairs into integer types, and efficient methods for
23//! reverse complement, enumerating kmers from longer sequences, and transfering data between
24//! sequences.
25//!
26//! ## Encodings
27//! Most methods for ingesting sequence data into the library have a form named 'bytes',
28//! which expects bases encoded as the integers 0,1,2,3, and a separate form names 'ascii',
29//! which expects bases encoded as the ASCII letters A,C,G,T.
30
31use bimap::BiMap;
32use clap::ValueEnum;
33use serde_derive::{Deserialize, Serialize};
34use summarizer::Marker;
35use std::fmt::{self, Debug, Display};
36use std::hash::Hash;
37use std::marker::PhantomData;
38use std::{array, mem};
39use std::ops::Range;
40
41use crate::compression::{CheckCompress, compress_kmers_with_hash};
42use crate::dna_string::DnaString;
43use crate::filter::filter_kmers;
44use crate::reads::{ReadData, Reads, ReadsPaired};
45use crate::serde::{SerGraph, SerKmers, SerReads};
46use crate::summarizer::{ID, SampleInfo, SummaryConfig, SummaryData, Tag, Translator};
47
48pub mod clean_graph;
49pub mod compression;
50pub mod dna_string;
51pub mod reads;
52pub mod filter;
53pub mod summarizer;
54pub mod graph;
55pub mod kmer;
56pub mod msp;
57pub mod neighbors;
58pub mod vmer;
59pub mod fastq;
60pub mod colors;
61pub mod serde;
62
63const BUF: usize = 64*1024;
64const BUCKETS: usize = 256;
65const ALPHABET_SIZE: usize = 4;
66const PROGRESS_STYLE: &str = "{msg} [{elapsed_precise}] {bar:60.cyan/blue} ({pos}/{len})";
67
68#[cfg(any(target_arch = "x86", target_arch = "x86_64"))]
69mod bitops_avx2;
70
71#[cfg(test)]
72pub mod test;
73
74/// Convert a 2-bit representation of a base to a char
75#[inline]
76pub fn bits_to_ascii(c: u8) -> u8 {
77    match c {
78        0u8 => b'A',
79        1u8 => b'C',
80        2u8 => b'G',
81        3u8 => b'T',
82        _ => b'X',
83    }
84}
85
86/// Convert an ASCII-encoded DNA base to a 2-bit representation,
87/// transforming bytes outside of ACGTacgt to A
88#[inline]
89pub fn base_to_bits(c: u8) -> u8 {
90    match c {
91        b'A' | b'a' => 0u8,
92        b'C' | b'c' => 1u8,
93        b'G' | b'g' => 2u8,
94        b'T' | b't' => 3u8,
95        _ => 0u8,
96    }
97}
98
99/// Convert an ASCII-encoded DNA base to a 2-bit representation,
100/// second value is `false` if the base was ambiguous
101#[inline]
102pub fn base_to_bits_checked(c: u8) -> (u8, bool) {
103    match c {
104        b'A' | b'a' => (0u8, true),
105        b'C' | b'c' => (1u8, true),
106        b'G' | b'g' => (2u8, true),
107        b'T' | b't' => (3u8, true),
108        _ => (0u8, false)
109    }
110}
111
112#[inline]
113pub fn dna_only_base_to_bits(c: u8) -> Option<u8> {
114    match c {
115        b'A' | b'a' => Some(0u8),
116        b'C' | b'c' => Some(1u8),
117        b'G' | b'g' => Some(2u8),
118        b'T' | b't' => Some(3u8),
119        _ => None,
120    }
121}
122
123/// Convert an ASCII-encoded DNA base to a 2-bit representation
124#[inline]
125pub fn is_valid_base(c: u8) -> bool {
126    matches!(c, b'A' | b'C' | b'G' | b'T' | b'a' | b'c' | b'g' | b't')
127}
128
129/// Convert a 2-bit representation of a base to a char
130#[inline]
131pub fn bits_to_base(c: u8) -> char {
132    match c {
133        0u8 => 'A',
134        1u8 => 'C',
135        2u8 => 'G',
136        3u8 => 'T',
137        _ => 'X',
138    }
139}
140
141/// The complement of a 2-bit encoded base
142#[inline(always)]
143pub fn complement(base: u8) -> u8 {
144    (!base) & 0x3u8
145}
146
147/// Trait for interacting with DNA sequences
148pub trait Mer: Sized + fmt::Debug {
149    /// Length of DNA sequence
150    fn len(&self) -> usize;
151
152    /// True if the sequence is empty.
153    fn is_empty(&self) -> bool;
154
155    /// Get 2-bit encoded base at position `pos`
156    fn get(&self, pos: usize) -> u8;
157
158    /// Set base at `pos` to 2-bit encoded base `val`
159    fn set_mut(&mut self, pos: usize, val: u8);
160
161    /// Set `nbases` positions in the sequence, starting at `pos`.
162    /// Values must  be packed into the upper-most bits of `value`.
163    fn set_slice_mut(&mut self, pos: usize, nbases: usize, value: u64);
164
165    /// Return a new object containing the reverse complement of the sequence
166    fn rc(&self) -> Self;
167
168    /// Iterate over the bases in the sequence
169    fn iter(&'_ self) -> MerIter<'_, Self> {
170        MerIter {
171            sequence: self,
172            i: 0,
173        }
174    }
175
176    /// Count the number of A/T bases in the kmer
177    fn at_count(&self) -> u32 {
178        let mut count = 0;
179        for i in 0..self.len() {
180            let base = self.get(i);
181            if base == 0 || base == 3 {
182                count += 1;
183            }
184        }
185        count
186    }
187
188    /// Count the number of G/C bases in the kmer
189    fn gc_count(&self) -> u32 {
190        let mut count = 0;
191        for i in 0..self.len() {
192            let base = self.get(i);
193            if base == 1 || base == 2 {
194                count += 1;
195            }
196        }
197        count
198    }
199}
200
201/// Iterator over bases of a DNA sequence (bases will be unpacked into bytes).
202pub struct MerIter<'a, M: 'a + Mer> {
203    sequence: &'a M,
204    i: usize,
205}
206
207impl<'a, M: 'a + Mer> Iterator for MerIter<'a, M> {
208    type Item = u8;
209
210    fn next(&mut self) -> Option<u8> {
211        if self.i < self.sequence.len() {
212            let value = self.sequence.get(self.i);
213            self.i += 1;
214            Some(value)
215        } else {
216            None
217        }
218    }
219}
220
221/// Encapsulates a Kmer sequence with statically known K.
222pub trait Kmer: Mer + Sized + Copy + PartialEq + PartialOrd + Eq + Ord + Hash {
223    /// Create a Kmer initialized to all A's
224    fn empty() -> Self;
225
226    /// K value for this concrete type.
227    fn k() -> usize;
228
229    /// Return the rank of this kmer in an lexicographic ordering of all kmers
230    /// E.g. 'AAAA' -> 0, 'AAAT' -> 1, etc. This will panic if K > 32.
231    fn to_u64(&self) -> u64;
232
233    /// Construct a kmer from the given lexicographic rank of the kmer.
234    /// If K > 32, the leads bases will be A's.
235    fn from_u64(value: u64) -> Self;
236
237    // Compute Hamming distance between self and other
238    fn hamming_dist(&self, other: Self) -> u32;
239
240    /// Add the base `v` to the left side of the sequence, and remove the rightmost base
241    fn extend_left(&self, v: u8) -> Self;
242
243    /// Add the base `v` to the right side of the sequence, and remove the leftmost base
244    fn extend_right(&self, v: u8) -> Self;
245
246    /// Add the base `v` to the side of sequence given by `dir`, and remove a base at the opposite side
247    fn extend(&self, v: u8, dir: Dir) -> Self {
248        match dir {
249            Dir::Left => self.extend_left(v),
250            Dir::Right => self.extend_right(v),
251        }
252    }
253
254    /// Generate all the extension of this sequence given by `exts` in direction `Dir`
255    fn get_extensions(&self, exts: Exts, dir: Dir) -> Vec<Self> {
256        let ext_bases = exts.get(dir);
257        ext_bases.iter().map(|b| self.extend(*b, dir)).collect()
258    }
259
260    /// Return the minimum of the kmer and it's reverse complement, and a flag indicating if sequence was flipped
261    fn min_rc_flip(&self) -> (Self, bool) {
262        let rc = self.rc();
263        //println!("kmer flip: self: {:?}, rc: {:?}, t/f: {}", self, rc, (*self < rc));
264        if *self < rc {
265            (*self, false)
266        } else {
267            (rc, true)
268        }
269    }
270
271    /// Return the minimum of the kmer and it's reverse complement
272    fn min_rc(&self) -> Self {
273        let rc = self.rc();
274        if *self < rc {
275            *self
276        } else {
277            rc
278        }
279    }
280
281    /// Test if this Kmer and it's reverse complement are the same
282    fn is_palindrome(&self) -> bool {
283        self.len().is_multiple_of(2) && *self == self.rc()
284    }
285
286    /// Create a Kmer from the first K bytes of `bytes`, which must be encoded as the integers 0-4.
287    fn from_bytes(bytes: &[u8]) -> Self {
288        if bytes.len() < Self::k() {
289            panic!("bytes not long enough to form kmer")
290        }
291
292        let mut k0 = Self::empty();
293
294        for (i, b) in bytes.iter().take(Self::k()).enumerate() {
295            k0.set_mut(i, *b)
296        }
297
298        k0
299    }
300
301    /// Create a Kmer from the first K bytes of `bytes`, which must be encoded as ASCII letters A,C,G, or T.
302    fn from_ascii(bytes: &[u8]) -> Self {
303        if bytes.len() < Self::k() {
304            panic!("bytes not long enough to form kmer")
305        }
306
307        let mut k0 = Self::empty();
308
309        for (i, b) in bytes.iter().take(Self::k()).enumerate() {
310            k0.set_mut(i, base_to_bits(*b))
311        }
312
313        k0
314    }
315
316    /// Return String containing Kmer sequence
317    fn to_string(&self) -> String {
318        let mut s = String::with_capacity(self.len());
319        for pos in 0..self.len() {
320            s.push(bits_to_base(self.get(pos)))
321        }
322        s
323    }
324
325    /// Generate vector of all kmers contained in `str` encoded as 0-4.
326    fn kmers_from_bytes(str: &[u8]) -> Vec<Self> {
327        if str.len() < Self::k() {
328            return Vec::default();
329        }
330        let mut k0 = Self::empty();
331        for (i, v) in str.iter().take(Self::k()).enumerate() {
332            k0.set_mut(i, *v);
333        }
334
335        let mut r = Vec::with_capacity(str.len() - Self::k() + 1);
336        r.push(k0);
337
338        for v in str.iter().skip(Self::k()) {
339            k0 = k0.extend_right(*v);
340            r.push(k0);
341        }
342
343        r
344    }
345
346    /// Generate vector of all kmers contained in `str`, encoded as ASCII ACGT.
347    fn kmers_from_ascii(str: &[u8]) -> Vec<Self> {
348        if str.len() < Self::k() {
349            return Vec::default();
350        }
351        let mut k0 = Self::empty();
352        for (i, b) in str.iter().take(Self::k()).enumerate() {
353            k0.set_mut(i, base_to_bits(*b));
354        }
355
356        let mut r = Vec::with_capacity(str.len() - Self::k() + 1);
357        r.push(k0);
358
359        for v in str.iter().skip(Self::k()) {
360            k0 = k0.extend_right(base_to_bits(*v));
361            r.push(k0);
362        }
363
364        r
365    }
366
367    fn has_low_complexity(&self) -> bool {
368        let a = Self::from_u64(0);
369        let c = Self::from_u64((0..(Self::k()*2)).filter(|&x| (x % 2 == 0) | (x == 0)).map(|x| 2usize.pow(x as u32)).sum::<usize>() as u64);
370        let g = Self::from_u64((0..(Self::k()*2)).filter(|&x| x % 2 != 0).map(|x| 2u64.pow(x as u32)).sum::<u64>());
371        let t = Self::from_u64(2u64.pow((Self::k()*2) as u32) - 1);
372
373        (self == &a) | (self == &c) | (self == &g) | (self == &t)
374    }
375}
376
377/// An immutable interface to a Mer sequence.
378pub trait MerImmut: Mer + Clone {
379    fn set(&self, pos: usize, val: u8) -> Self {
380        let mut new = self.clone();
381        new.set_mut(pos, val);
382        new
383    }
384
385    fn set_slice(&self, pos: usize, nbases: usize, bits: u64) -> Self {
386        let mut new = self.clone();
387        new.set_slice_mut(pos, nbases, bits);
388        new
389    }
390}
391
392impl<T> MerImmut for T where T: Mer + Clone {}
393
394/// A DNA sequence with run-time variable length, up to a statically known maximum length
395pub trait Vmer: Mer + PartialEq + Eq {
396    /// Create a new sequence with length `len`, initialized to all A's
397    fn new(len: usize) -> Self;
398
399    /// Maximum sequence length that can be stored in this type
400    fn max_len() -> usize;
401
402    /// Create a Vmer from a sequence of bytes
403    fn from_slice(seq: &[u8]) -> Self {
404        let mut vmer = Self::new(seq.len());
405        for (i, v) in seq.iter().enumerate() {
406            vmer.set_mut(i, *v);
407        }
408
409        vmer
410    }
411
412    /// Efficiently extract a Kmer from the sequence
413    fn get_kmer<K: Kmer>(&self, pos: usize) -> K;
414
415    /// Get the first Kmer from the sequence
416    fn first_kmer<K: Kmer>(&self) -> K {
417        self.get_kmer(0)
418    }
419
420    /// Get the last kmer in the sequence
421    fn last_kmer<K: Kmer>(&self) -> K {
422        self.get_kmer(self.len() - K::k())
423    }
424
425    /// Get the terminal kmer of the sequence, on the both side of the sequence
426    fn both_term_kmer<K: Kmer>(&self) -> (K, K) {
427        (self.first_kmer(), self.last_kmer())
428    }
429
430    /// Get the terminal kmer of the sequence, on the side of the sequence given by dir
431    fn term_kmer<K: Kmer>(&self, dir: Dir) -> K {
432        match dir {
433            Dir::Left => self.first_kmer(),
434            Dir::Right => self.last_kmer(),
435        }
436    }
437
438    /// Iterate over the kmers in the sequence
439    fn iter_kmers<K: Kmer>(&self) -> KmerIter<'_, K, Self> {
440        let kmer = if self.len() >= K::k() {
441            self.first_kmer()
442        } else {
443            // Default kmer, will not be used
444            K::empty()
445        };
446
447        KmerIter {
448            bases: self,
449            kmer,
450            pos: K::k(),
451        }
452    }
453
454    /// Iterate over the kmers and their extensions, given the extensions of the whole sequence
455    fn iter_kmer_exts<K: Kmer>(&self, seq_exts: Exts) -> KmerExtsIter<'_, K, Self> {
456        let kmer = if self.len() >= K::k() {
457            self.first_kmer()
458        } else {
459            // Default kmer, will not be used
460            K::empty()
461        };
462
463        KmerExtsIter {
464            bases: self,
465            exts: seq_exts,
466            kmer,
467            pos: K::k(),
468        }
469    }
470}
471
472#[derive(Debug, Clone, Copy)]
473pub struct KmerDataItem<K: Kmer, DI> {
474    pub kmer: K,
475    pub exts: Exts,
476    pub data: DI,
477    pub quality: Option<BaseQuality>
478}
479
480impl<K: Kmer, DI> KmerDataItem<K, DI> {
481    pub fn new(kmer: K, exts: Exts, data: DI, quality: Option<BaseQuality>) -> KmerDataItem<K, DI> {
482        KmerDataItem { kmer, exts, data, quality }
483    }
484}
485
486/// A newtype wrapper around a `Vec<u8>` with implementations
487/// of the `Mer` and `Vmer` traits.
488#[derive(Debug, Clone, Eq, PartialEq, Ord, PartialOrd)]
489pub struct DnaBytes(pub Vec<u8>);
490
491impl Mer for DnaBytes {
492    fn len(&self) -> usize {
493        self.0.len()
494    }
495
496    fn is_empty(&self) -> bool {
497        self.0.is_empty()
498    }
499
500    fn get(&self, pos: usize) -> u8 {
501        self.0[pos]
502    }
503
504    /// Set base at `pos` to 2-bit encoded base `val`
505    fn set_mut(&mut self, pos: usize, val: u8) {
506        self.0[pos] = val
507    }
508
509    /// Set `nbases` positions in the sequence, starting at `pos`.
510    /// Values must  be packed into the upper-most bits of `value`.
511    fn set_slice_mut(&mut self, _pos: usize, _nbases: usize, _value: u64) {
512        unimplemented!();
513        //for i in pos .. (pos + nbases) {
514        //
515        //}
516    }
517
518    /// Return a new object containing the reverse complement of the sequence
519    fn rc(&self) -> Self {
520        unimplemented!();
521    }
522}
523
524impl Vmer for DnaBytes {
525    /// Create a new sequence with length `len`, initialized to all A's
526    fn new(len: usize) -> Self {
527        DnaBytes(vec![0u8; len])
528    }
529
530    /// Maximum sequence length that can be stored in this type
531    fn max_len() -> usize {
532        1 << 48
533    }
534
535    /// Efficiently extract a Kmer from the sequence
536    fn get_kmer<K: Kmer>(&self, pos: usize) -> K {
537        K::from_bytes(&self.0[pos..pos + K::k()])
538    }
539}
540
541/// A newtype wrapper around a `&[u8]` with implementations
542/// of the `Mer` and `Vmer` traits.
543#[derive(Debug, Eq, PartialEq, Ord, PartialOrd)]
544pub struct DnaSlice<'a>(pub &'a [u8]);
545
546impl Mer for DnaSlice<'_> {
547    fn len(&self) -> usize {
548        self.0.len()
549    }
550
551    fn is_empty(&self) -> bool {
552        self.0.is_empty()
553    }
554
555    fn get(&self, pos: usize) -> u8 {
556        self.0[pos]
557    }
558
559    /// Set base at `pos` to 2-bit encoded base `val`
560    fn set_mut(&mut self, _pos: usize, _val: u8) {
561        unimplemented!()
562    }
563
564    /// Set `nbases` positions in the sequence, starting at `pos`.
565    /// Values must  be packed into the upper-most bits of `value`.
566    fn set_slice_mut(&mut self, _pos: usize, _nbases: usize, _value: u64) {
567        unimplemented!();
568        //for i in pos .. (pos + nbases) {
569        //
570        //}
571    }
572
573    /// Return a new object containing the reverse complement of the sequence
574    fn rc(&self) -> Self {
575        unimplemented!();
576    }
577}
578
579impl Vmer for DnaSlice<'_> {
580    /// Create a new sequence with length `len`, initialized to all A's
581    fn new(_len: usize) -> Self {
582        unimplemented!();
583    }
584
585    /// Maximum sequence length that can be stored in this type
586    fn max_len() -> usize {
587        1 << 48
588    }
589
590    /// Efficiently extract a Kmer from the sequence
591    fn get_kmer<K: Kmer>(&self, pos: usize) -> K {
592        K::from_bytes(&self.0[pos..pos + K::k()])
593    }
594}
595
596/// Direction of motion in a DeBruijn graph
597#[derive(Copy, Clone, Debug, Serialize, Deserialize, PartialEq)]
598pub enum Dir {
599    Left,
600    Right,
601}
602
603impl Dir {
604    /// Return a fresh Dir with the opposite direction
605    pub fn flip(&self) -> Dir {
606        match *self {
607            Dir::Left => Dir::Right,
608            Dir::Right => Dir::Left,
609        }
610    }
611
612    /// Return a fresh Dir opposite direction if do_flip == True
613    pub fn cond_flip(&self, do_flip: bool) -> Dir {
614        if do_flip {
615            self.flip()
616        } else {
617            *self
618        }
619    }
620
621    /// Pick between two alternatives, depending on the direction
622    pub fn pick<T>(&self, if_left: T, if_right: T) -> T {
623        match self {
624            Dir::Left => if_left,
625            Dir::Right => if_right,
626        }
627    }
628
629    /// get the index of the base in dir for [`Exts`] and [`EdgeMult`]
630    fn index(&self, base: u8) -> u8 {
631        match self {
632            Self::Right => ALPHABET_SIZE as u8 - 1 - base,
633            Self::Left => 2 * ALPHABET_SIZE as u8 - 1 - base,
634        }
635    }
636
637    /// get the index range of the dir for [`Exts`] and [`EdgeMult`]
638    fn index_range(&self) -> Range<usize>{
639        match self {
640            Self::Right => 0..ALPHABET_SIZE,
641            Self::Left => ALPHABET_SIZE..(2 * ALPHABET_SIZE),
642        }
643    }
644}
645
646/// Store single-base extensions for a DNA Debruijn graph.
647///
648/// 8 bits, 4 higher order ones represent extensions to the right, 4 lower order ones
649/// represent extensions to the left. For each direction the bits (from lower order
650/// to higher order) represent whether there exists an extension with each of the
651/// letters A, C, G, T. So overall the bits are:
652///  right   left
653/// T G C A T G C A
654#[derive(Eq, PartialEq, Copy, Clone, Ord, PartialOrd, Hash, Serialize, Deserialize)]
655pub struct Exts {
656    pub val: u8,
657}
658
659impl Exts {
660    pub fn new(val: u8) -> Self {
661        Exts { val }
662    }
663
664    pub fn empty() -> Exts {
665        Exts { val: 0u8 }
666    }
667
668    pub fn from_single_dirs(left: Exts, right: Exts) -> Exts {
669        Exts {
670            val: (right.val << 4) | (left.val & 0xf),
671        }
672    }
673
674    pub fn merge(left: Exts, right: Exts) -> Exts {
675        Exts {
676            val: left.val & 0x0f | right.val & 0xf0,
677        }
678    }
679
680    pub fn add(&self, v: Exts) -> Exts {
681        Exts {
682            val: self.val | v.val,
683        }
684    }
685
686    /// subtract an Exts from an Exts
687    pub fn subtract(&self, v: Exts) -> Exts {
688        Exts { val: self.val & !v.val }
689    }
690
691    pub fn set(&self, dir: Dir, pos: u8) -> Exts {
692        let shift = pos
693            + match dir {
694                Dir::Right => 4,
695                Dir::Left => 0,
696            };
697
698        let new_val = self.val | (1u8 << shift);
699        Exts { val: new_val }
700    }
701
702    pub fn remove(&self, dir: Dir, pos: u8) -> Exts {
703        let shift = pos
704            + match dir {
705                Dir::Right => 4,
706                Dir::Left => 0,
707            };
708
709        let new_val = self.val & !(1u8 << shift);
710        Exts { val: new_val }
711    }
712
713    #[inline]
714    fn dir_bits(&self, dir: Dir) -> u8 {
715        match dir {
716            Dir::Right => self.val >> 4,
717            Dir::Left => self.val & 0xf,
718        }
719    }
720
721    pub fn get(&self, dir: Dir) -> Vec<u8> {
722        let bits = self.dir_bits(dir);
723        let mut v = Vec::with_capacity(4);
724        for i in 0..4 {
725            if bits & (1 << i) > 0 {
726                v.push(i);
727            }
728        }
729
730        v
731    }
732
733    pub fn has_ext(&self, dir: Dir, base: u8) -> bool {
734        let bits = self.dir_bits(dir);
735        (bits & (1 << base)) > 0
736    }
737
738    pub fn from_slice_bounds(src: &[u8], start: usize, length: usize) -> Exts {
739        let l_extend = if start > 0 {
740            1u8 << (src[start - 1])
741        } else {
742            0u8
743        };
744        let r_extend = if start + length < src.len() {
745            1u8 << src[start + length]
746        } else {
747            0u8
748        };
749
750        Exts {
751            val: (r_extend << 4) | l_extend,
752        }
753    }
754
755    pub fn from_dna_string(src: &dna_string::DnaString, start: usize, length: usize) -> Exts {
756        let l_extend = if start > 0 {
757            1u8 << (src.get(start - 1))
758        } else {
759            0u8
760        };
761        let r_extend = if start + length < src.len() {
762            1u8 << src.get(start + length)
763        } else {
764            0u8
765        };
766
767        Exts {
768            val: (r_extend << 4) | l_extend,
769        }
770    }
771
772    pub fn num_exts_l(&self) -> u8 {
773        self.num_ext_dir(Dir::Left)
774    }
775
776    pub fn num_exts_r(&self) -> u8 {
777        self.num_ext_dir(Dir::Right)
778    }
779
780    pub fn num_ext_dir(&self, dir: Dir) -> u8 {
781        let e = self.dir_bits(dir);
782        (e & 1u8) + ((e & 2u8) >> 1) + ((e & 4u8) >> 2) + ((e & 8u8) >> 3)
783    }
784
785    pub fn mk_left(base: u8) -> Exts {
786        Exts::empty().set(Dir::Left, base)
787    }
788
789    pub fn mk_right(base: u8) -> Exts {
790        Exts::empty().set(Dir::Right, base)
791    }
792
793    pub fn mk(left_base: u8, right_base: u8) -> Exts {
794        Exts::merge(Exts::mk_left(left_base), Exts::mk_right(right_base))
795    }
796
797    pub fn get_unique_extension(&self, dir: Dir) -> Option<u8> {
798        if self.num_ext_dir(dir) != 1 {
799            None
800        } else {
801            let e = self.dir_bits(dir);
802            for i in 0..4 {
803                if (e & (1 << i)) > 0 {
804                    return Some(i);
805                }
806            }
807
808            None
809        }
810    }
811
812    pub fn single_dir(&self, dir: Dir) -> Exts {
813        match dir {
814            Dir::Right => Exts { val: self.val >> 4 },
815            Dir::Left => Exts {
816                val: self.val & 0xfu8,
817            },
818        }
819    }
820
821    /// Complement the extension bases for each direction
822    pub fn complement(&self) -> Exts {
823        let v = self.val;
824
825        // swap bits
826        let mut r = (v & 0x55u8) << 1 | ((v >> 1) & 0x55u8);
827
828        // swap pairs
829        r = (r & 0x33u8) << 2 | ((r >> 2) & 0x33u8);
830        Exts { val: r }
831    }
832
833    pub fn reverse(&self) -> Exts {
834        let v = self.val;
835        let r = (v & 0xf) << 4 | (v >> 4);
836        Exts { val: r }
837    }
838
839    pub fn rc(&self) -> Exts {
840        self.reverse().complement()
841    }
842}
843
844impl fmt::Debug for Exts {
845    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
846        let mut s = String::new();
847
848        for b in self.get(Dir::Left) {
849            s.push(bits_to_base(b));
850        }
851        s.push('|');
852
853        for b in self.get(Dir::Right) {
854            s.push(bits_to_base(b));
855        }
856
857        write!(f, "{}", s)
858    }
859}
860
861/// Iterate over the `Kmer`s of a DNA sequence efficiently
862pub struct KmerIter<'a, K: Kmer, D>
863where
864    D: 'a,
865{
866    bases: &'a D,
867    kmer: K,
868    pos: usize,
869}
870
871impl<K: Kmer, D: Mer> Iterator for KmerIter<'_, K, D> {
872    type Item = K;
873
874    #[inline]
875    fn next(&mut self) -> Option<K> {
876        if self.pos <= self.bases.len() {
877            let retval = self.kmer;
878
879            if self.pos < self.bases.len() {
880                self.kmer = self.kmer.extend_right(self.bases.get(self.pos));
881            }
882
883            self.pos += 1;
884            Some(retval)
885        } else {
886            None
887        }
888    }
889}
890
891/// Iterate over the `(Kmer, Exts)` tuples of a sequence and it's extensions efficiently
892pub struct KmerExtsIter<'a, K: Kmer, D>
893where
894    D: 'a,
895{
896    bases: &'a D,
897    exts: Exts,
898    kmer: K,
899    pos: usize,
900}
901
902impl<K: Kmer, D: Mer> Iterator for KmerExtsIter<'_, K, D> {
903    type Item = (K, Exts);
904
905    fn next(&mut self) -> Option<(K, Exts)> {
906        if self.pos <= self.bases.len() {
907            let next_base = if self.pos < self.bases.len() {
908                self.bases.get(self.pos)
909            } else {
910                0u8
911            };
912
913            let cur_left = if self.pos == K::k() {
914                self.exts
915            } else {
916                Exts::mk_left(self.bases.get(self.pos - K::k() - 1))
917            };
918
919            let cur_right = if self.pos < self.bases.len() {
920                Exts::mk_right(next_base)
921            } else {
922                self.exts
923            };
924
925            let cur_exts = Exts::merge(cur_left, cur_right);
926
927            let retval = self.kmer;
928            self.kmer = self.kmer.extend_right(next_base);
929            self.pos += 1;
930            Some((retval, cur_exts))
931        } else {
932            None
933        }
934    }
935}
936
937
938/// Compress up to 64 tags to one `u64` (8 bytes) (or up to 128 tags to a  
939/// `u128`(16 bytes) with the feature `sample128` enabled)
940#[derive(Clone, PartialEq, Copy, Serialize, Deserialize)]
941pub struct Tags {
942    pub val: Marker,
943}
944
945impl Tags {
946
947    /// Make a new Tags from a `u64` (or `u128` with the feature `sample128` enabled)
948    pub fn new(val: Marker) -> Self {
949        Tags { val }
950    }
951
952    /// get number of labels saved in the Tags
953    pub fn len(&self) -> usize {
954        self.val.count_ones() as usize
955    }
956
957    /// check if the tags are empty
958    pub fn is_empty(&self) -> bool {
959        self.val == 0
960    }
961
962    /// encodes a sorted (!) Vec<Tag> and encodes it as a u64
963    pub fn from_tag_vec(vec: &Vec<Tag>) -> Self {
964        let mut x = 0;
965        
966        // if the vector is empty, return an empty Tags
967        if vec.is_empty() { return Tags { val: 0 } }
968
969        // panic if Tags would overflow
970        if ( *vec.last().expect("vector empty when it shouldn't be") ) / 8 as Tag >= mem::size_of::<Tags>() as Tag { 
971            panic!("too many tags - maximum number of supported tags is 64 by default, 128 with compile flag / feature '--feature sample128'") 
972        }
973        
974        // iterate backwards over all elements of the vector
975        for i in (1..vec.len()).rev() {
976            x += 1;
977            x <<= vec[i] - vec[i-1];
978        }
979
980        x += 1;
981        x <<= vec[0];
982
983        Tags { val: x }
984    }
985
986    // turn Tags into Vec<Tag>
987    pub fn to_tag_vec(&self) -> Vec<Tag> {
988        let mut x = self.val;
989        let mut vec: Vec<Tag> = Vec::new();
990
991        // do bit-wise right shifts trough u64
992        // each time first digit is 1 (is an odd number), push i to vec
993        for i in 0..(mem::size_of::<Tags>()*8) as Tag {
994            if !x.is_multiple_of(2) {
995                vec.push(i)
996            }
997            x >>= 1;
998        }
999
1000        vec
1001    }
1002
1003    // directly translate Tags to Vec<&str>
1004    // str_map is translatror BiMap between Tag and &str 
1005    pub fn to_string_vec<'a>(&'a self, str_map: &'a BiMap<String, Tag>) -> Vec<&'a str> {
1006        let mut x = self.val;
1007        let mut vec: Vec<&str> = Vec::with_capacity(x.count_ones() as usize);
1008
1009        // iterate through bits of the u64
1010        for i in 0..(mem::size_of::<Tags>()*8) as Tag {
1011            // check if odd number: current first bit is 1
1012            if !x.is_multiple_of(2) {
1013                match str_map.get_by_right(&{ i }) {
1014                    Some(label) => vec.push(label),
1015                    None => panic!("tried to access label that does not exist!"),
1016                }
1017            }
1018            // shift the u64 bitise to rotate though it
1019            x >>= 1;
1020        }
1021        vec    
1022    }
1023
1024
1025    /// compares the value of the tags with another value (marker) with a bit-wise and,
1026    /// returns true if the result is greater than 0:
1027    /// `00101 & 01000 -> false`
1028    /// `00101 & 00100 -> true`
1029    pub fn bit_and(&self, marker: Marker) -> bool {
1030        (self.val & marker) > 0
1031    }
1032
1033    /// compares the value of the tags with another value (marker) with a bit-wise and
1034    /// counts the overlaps:
1035    /// `00101 & 01000 -> 0`
1036    /// `00101 & 00100 -> 1`
1037    /// `00101 & 00101 -> 2`
1038    pub fn bit_and_dist(&self, marker: Marker) -> usize {
1039        (self.val & marker).count_ones() as usize
1040    }
1041
1042    /// get an iterator over the tags in the [`Tags`]
1043    pub fn iter(&self) -> TagsIterator {
1044        TagsIterator::new(*self)
1045    }
1046
1047    /// get the memory of the [`Tags`] (depends on activated features)
1048    pub fn mem(&self) -> usize {
1049        mem::size_of::<Marker>()
1050    }
1051}
1052
1053impl fmt::Debug for Tags {
1054    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
1055        write!(f, "{:?}", self.to_tag_vec())
1056    }
1057}
1058
1059pub struct TagsIterator {
1060    tags: Tags,
1061    i: Tag
1062}
1063
1064impl TagsIterator  {
1065    fn new(tags: Tags) -> TagsIterator {
1066        TagsIterator {tags, i: 0}
1067    }
1068}
1069
1070impl Iterator for TagsIterator {
1071    type Item = Tag;
1072
1073    fn next(&mut self) -> Option<Self::Item> {
1074        loop {
1075            if self.i as usize == mem::size_of::<Tags>()*8 { return None }
1076            let result = !self.tags.val.is_multiple_of(2);
1077            self.tags.val >>= 1;
1078            self.i += 1;
1079            if result { return Some(self.i - 1); }
1080        }
1081    }
1082}
1083
1084pub struct TagsFormatter<'a> {
1085    tags: Tags,
1086    translator: &'a Translator
1087}
1088
1089impl<'a> TagsFormatter<'a> {
1090    pub fn new(tags: Tags, translator: &'a Translator) -> TagsFormatter<'a> {
1091        TagsFormatter {
1092            tags,
1093            translator
1094        }
1095    }
1096}
1097
1098impl fmt::Display for TagsFormatter<'_> {
1099    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
1100        if let Some(tag_translator) = self.translator.tag_translator() {
1101            let tag_vec = self.tags.to_string_vec(tag_translator);
1102
1103            writeln!(f, "samples:")?;
1104
1105            for label in tag_vec.into_iter() {
1106                writeln!(f, "{}", label)?
1107            }
1108        } else {
1109            write!(f, "samples: {:?}", self.tags.to_tag_vec())?
1110        }
1111
1112        Ok(())
1113    }
1114}
1115
1116pub struct TagsCountsFormatter<'a> {
1117    tags: Tags,
1118    counts: &'a [u32],
1119    translator: &'a Translator
1120}
1121
1122impl<'a> TagsCountsFormatter<'a> {
1123    pub fn new(tags: Tags, counts: &'a [u32], translator: &'a Translator) -> TagsCountsFormatter<'a> {
1124        TagsCountsFormatter {
1125            tags,
1126            counts,
1127            translator
1128        }
1129    }
1130}
1131
1132impl fmt::Display for TagsCountsFormatter<'_> {
1133    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
1134        writeln!(f, "{:<20} - counts", "samples")?;
1135
1136        if let Some(tag_translator) = self.translator.tag_translator() {
1137            let label_vec = self.tags.to_string_vec(tag_translator);
1138
1139            for (label, count) in label_vec.into_iter().zip(self.counts) {
1140                writeln!(f, "{:<20} - {}", label, count)?
1141            }
1142        } else {
1143            let tag_vec = self.tags.to_tag_vec();
1144
1145            for (tag, count) in tag_vec.into_iter().zip(self.counts) {
1146                writeln!(f, "{:<20} - {}", tag, count)?
1147            }
1148        }
1149
1150
1151        Ok(())
1152    }
1153}
1154
1155
1156// would be more intuitive with left and right switched but Exts were built this way
1157/// multiplicities or coverage for each of the 8 possible edges
1158/// indices: 
1159/// 0: T right
1160/// 1: G right
1161/// 2: C right
1162/// 3: A right
1163/// 4: T left
1164/// 5: G left
1165/// 6: C left
1166/// 7: A left
1167#[derive(PartialEq, PartialOrd, Eq, Ord, Serialize, Deserialize, Clone)]
1168pub struct EdgeMult {
1169    edge_mults: [u32; 2*ALPHABET_SIZE],
1170}
1171
1172impl Debug for EdgeMult {
1173    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
1174        let edge_f = ["A:", ", C:", ", G:", ", T:", " | A:", ", C:", ", G:", ", T:"];
1175        for (ef, em) in edge_f.iter().zip(self.edge_mults.iter().rev()) {
1176             write!(f, "{} {}", ef, em)?
1177        }
1178        Ok(())
1179    }
1180}
1181
1182impl Display for EdgeMult {
1183    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
1184        let base = ["A", "C", "G", "T"];
1185        for (i, b) in (0..ALPHABET_SIZE).rev().zip(base) {
1186            writeln!(f, "{}: {} | {}", 
1187                b, 
1188                self.edge_mults[i + ALPHABET_SIZE], 
1189                self.edge_mults[i]
1190            )?
1191        }
1192
1193        Ok(())
1194    }
1195}
1196
1197impl EdgeMult {
1198    /// a new, empty `EdgeMult`
1199    pub fn new() -> Self {
1200        EdgeMult { edge_mults: [0;  2 * ALPHABET_SIZE] }
1201    }
1202
1203    /// a new `EdgeMult` with values
1204    pub fn new_from(edge_mults: [u32; 2*ALPHABET_SIZE]) -> Self {
1205        EdgeMult { edge_mults }
1206    }
1207
1208    /// add a count to the an edge
1209    pub fn add(&mut self, base: u8, dir: Dir, count: u32) {
1210        self.edge_mults[dir.index(base) as usize] += count
1211    }
1212
1213    /// remove an edge from the `EdgeMult`
1214    pub fn remove(&mut self, base: u8, dir: Dir) {
1215        self.edge_mults[dir.index(base) as usize] = 0;
1216    }
1217
1218    /// add an [`Exts`] to the `EdgeMult`
1219    pub fn add_exts(&mut self, exts: Exts) {
1220        let mut exts = exts.val;
1221        for index in (0..(2 * ALPHABET_SIZE)).rev() {
1222            if !exts.is_multiple_of(2) {
1223                self.edge_mults[index] += 1
1224            }
1225            exts >>= 1;
1226        }
1227    }
1228
1229    /// get the edge multiplicities as an array 
1230    pub fn edge_mults(&self) -> [u32; 2*ALPHABET_SIZE] {
1231        self.edge_mults
1232    }
1233
1234    /// get the edge multiplicities to the right of the node
1235    pub fn right(&self) -> &[u32] {
1236        &self.edge_mults[(Dir::Right).index_range()]
1237    }
1238
1239    /// get the edge multiplicities to the left of the node
1240    pub fn left(&self) -> &[u32] {
1241        &self.edge_mults[(Dir::Left).index_range()]
1242    }
1243
1244    /// get the sum of all edges of the node
1245    pub fn sum(&self) -> u32 {
1246        self.edge_mults.iter().sum::<u32>()
1247    }
1248
1249    /// get the multiplicity of a certain edge
1250    pub fn edge_mult(&self, base: u8, dir: Dir) -> u32 {
1251        self.edge_mults[dir.index(base) as usize]
1252    }
1253
1254    /// get the [`Exts`] corresponding to the `EdgeMult`
1255    pub fn exts(&self) -> Exts {
1256        let mut exts_val = 0u8;
1257        for (i, edge) in self.edge_mults.iter().rev().enumerate() {
1258            if *edge > 0 { exts_val += 2u8.pow(i as u32) }
1259        }
1260
1261        Exts::new(exts_val)
1262    }
1263
1264    /// clean the edge mults by removing edge counts that led to filtered kmers
1265    /// based on a correct [`Exts`]
1266    pub fn clean_edges(&mut self, exts: Exts) {
1267        let mut exts = exts.val;
1268        for index in (0..(2 * ALPHABET_SIZE)).rev() {
1269            if exts.is_multiple_of(2) {
1270                self.edge_mults[index] = 0;
1271            }
1272            exts >>= 1;
1273        }
1274        
1275    }
1276
1277    pub fn rc(&mut self) {
1278        self.edge_mults.reverse();
1279    }
1280
1281    pub fn combine(left: &EdgeMult, right: &EdgeMult) -> EdgeMult {
1282        let mut combined = [0u32; 2*ALPHABET_SIZE];
1283        (0..ALPHABET_SIZE).for_each(|i| combined[i] = right.edge_mults[i]);
1284        (ALPHABET_SIZE..(2*ALPHABET_SIZE)).for_each(|i| combined[i] = left.edge_mults[i]);
1285
1286        EdgeMult::new_from(combined)
1287    }
1288
1289    pub fn from_single_dirs(left: &Option<SingleDirEdgeMult>, right: &Option<SingleDirEdgeMult>) -> Option<EdgeMult> {
1290        if let Some(l_em) = left {
1291            if let Some(r_em) = right {
1292                let mut combined = [0u32; 2*ALPHABET_SIZE];
1293                (0..ALPHABET_SIZE).for_each(|i| combined[i] = r_em.edge_mults[i]);
1294                (0..ALPHABET_SIZE).for_each(|i| combined[i + ALPHABET_SIZE] = l_em.edge_mults[i]);
1295        
1296                return Some(EdgeMult::new_from(combined))
1297            }
1298        }
1299
1300        None
1301    }
1302
1303    pub fn single_dir(&self, dir: Dir) -> SingleDirEdgeMult {
1304        match dir {
1305            Dir::Left => SingleDirEdgeMult::new(self.edge_mults[ALPHABET_SIZE..(2*ALPHABET_SIZE)]
1306                .try_into().expect("Error: slice has incorrect length")),
1307            Dir::Right => SingleDirEdgeMult::new(self.edge_mults[0..ALPHABET_SIZE]
1308                .try_into().expect("Error: slice has incorrect length")),
1309        }
1310    }
1311}
1312
1313impl Default for EdgeMult {
1314    fn default() -> Self {
1315        Self::new()
1316    }
1317}
1318
1319#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord)]
1320pub struct SingleDirEdgeMult {
1321    edge_mults: [u32; ALPHABET_SIZE]
1322}
1323
1324impl SingleDirEdgeMult {
1325    pub fn new(edge_mults: [u32; ALPHABET_SIZE]) -> Self {
1326        SingleDirEdgeMult { edge_mults }
1327    }
1328
1329    pub fn complement(&self) -> Self {
1330        let mut reverse = self.edge_mults;
1331        reverse.reverse();
1332        SingleDirEdgeMult::new(reverse)
1333    }
1334
1335    /// get the multiplicity of a certain edge
1336    pub fn edge_mult(&self, base: u8) -> u32 {
1337        self.edge_mults[(ALPHABET_SIZE as u8 - 1 - base) as usize]
1338    }
1339}
1340
1341// would be more intuitive with left and right switched but Exts were built this way
1342/// mapped transcript/gene/chromosome IDs for each of the 8 possible edges
1343/// indices: 
1344/// 0: T right
1345/// 1: G right
1346/// 2: C right
1347/// 3: A right
1348/// 4: T left
1349/// 5: G left
1350/// 6: C left
1351/// 7: A left
1352#[derive(PartialEq, PartialOrd, Eq, Ord, Serialize, Deserialize, Clone)]
1353pub struct EdgeMap {
1354    edge_maps: [Box<[ID]>; 2*ALPHABET_SIZE],
1355}
1356
1357impl EdgeMap {
1358    /// create a new [`EdgeMap`] by supplying the underlying mapped IDs.
1359    pub fn new(edge_maps: [Box<[ID]>; 2*ALPHABET_SIZE]) -> EdgeMap {
1360        EdgeMap { edge_maps }
1361    }
1362
1363    /// get the IDs mapped to the specified edge
1364    fn edge_map(&self, base: u8, dir: Dir) -> &[ID] {
1365        &self.edge_maps[dir.index(base) as usize]
1366    }
1367
1368    /// set the IDs mapped to the edge at the specified index
1369    fn set_edge_map_at_index(&mut self, edge_map: Box<[ID]>, index: usize) {
1370        self.edge_maps[index] = edge_map;
1371    }
1372
1373    /// add an ID to the edge at the specified index
1374    fn add_id_to_edge_map_at_index(&mut self, id: ID, index: usize) {
1375        let mut em = self.edge_maps[index].to_vec();
1376        em.push(id);
1377
1378        self.set_edge_map_at_index(em.into(), index);
1379    }
1380
1381    /// add an ID at an [`Exts`] to the [`EdgeMap`]
1382    pub fn add_id(&mut self, exts: Exts, id: ID) {
1383        let mut exts = exts.val;
1384        for index in (0..(2 * ALPHABET_SIZE)).rev() {
1385            if !exts.is_multiple_of(2) {
1386                self.add_id_to_edge_map_at_index(id, index);
1387            }
1388            exts >>= 1;
1389        }
1390    }
1391
1392    /// returns true if the [`EdgeMap`] does not contain any IDs
1393    pub fn is_empty(&self) -> bool {
1394        let mut empty = true;
1395
1396        for emap in self.edge_maps.iter() {
1397            if !emap.is_empty() { empty = false }
1398        }
1399
1400        empty
1401    }
1402
1403    /// returns the heap memory used by the [`EdgeMap`] - the stack memory size
1404    /// is always 128 bytes (8 edges * (8 byte pointer + 8 byte length))
1405    pub fn mem_heap(&self) -> usize {
1406        let mut heap = 0;
1407
1408        for emap in self.edge_maps.iter() {
1409            heap += mem::size_of_val(&**emap);
1410        }
1411
1412        heap
1413    }
1414
1415    /// create an [`EdgeMap`] from two [`SingleDirEdgeMap`]s
1416    pub fn from_single_dirs(left: &Option<SingleDirEdgeMap>, right: &Option<SingleDirEdgeMap>) -> Option<EdgeMap> {
1417        if let Some(l_em) = left {
1418            if let Some(r_em) = right {
1419                let mut combined = EdgeMap::default().edge_maps;
1420                (0..ALPHABET_SIZE).for_each(|i| combined[i] = r_em.edge_maps[i].clone());
1421                (0..ALPHABET_SIZE).for_each(|i| combined[i + ALPHABET_SIZE] = l_em.edge_maps[i].clone());
1422        
1423                return Some(EdgeMap::new(combined))
1424            }
1425        }
1426
1427        None
1428    }
1429
1430    /// get the [`SingleDirEdgeMap`] in the specified [`Dir`]
1431    pub fn single_dir(&self, dir: Dir) -> SingleDirEdgeMap {
1432        let singe_dir: &[Box<[ID]>; 4] = match dir {
1433            Dir::Left => self.edge_maps[ALPHABET_SIZE..(2*ALPHABET_SIZE)].try_into().expect("Error: slice has incorrect length"),
1434            Dir::Right => self.edge_maps[0..ALPHABET_SIZE].try_into().expect("Error: slice has incorrect length"),
1435        };
1436
1437        SingleDirEdgeMap::new(singe_dir.clone())
1438    }
1439
1440    /// clean the edge maps by removing IDs mapped to edges that led to filtered kmers,
1441    /// based on a correct [`Exts`]
1442    pub fn clean_edges(&mut self, exts: Exts) {
1443        let mut exts = exts.val;
1444        for index in (0..(2 * ALPHABET_SIZE)).rev() {
1445            if exts.is_multiple_of(2) {
1446                self.edge_maps[index] = [].into();
1447            }
1448            exts >>= 1;
1449        }
1450        
1451    }
1452}
1453
1454impl Default for EdgeMap {
1455    fn default() -> Self {
1456        let edge_maps = array::from_fn(|_n| Vec::new().into());
1457        Self { edge_maps }
1458    }
1459}
1460
1461impl Debug for EdgeMap {
1462    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
1463        let edge_f = ["A:", ", C:", ", G:", ", T:", " | A:", ", C:", ", G:", ", T:"];
1464        for (ef, em) in edge_f.iter().zip(self.edge_maps.iter().rev()) {
1465             write!(f, "{} {:?}", ef, em)?
1466        }
1467        Ok(())
1468    }
1469}
1470
1471impl Display for EdgeMap {
1472    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
1473        let base = ["A", "C", "G", "T"];
1474        for (i, b) in (0..ALPHABET_SIZE).rev().zip(base) {
1475            writeln!(f, "{}: {:?} | {:?}", 
1476                b, 
1477                self.edge_maps[i + ALPHABET_SIZE], 
1478                self.edge_maps[i]
1479            )?
1480        }
1481
1482        Ok(())
1483    }
1484}
1485
1486/// mapped transcript/gene/chromosome IDs for each of the 4 possible edges
1487/// indices in one direction: 
1488/// 
1489/// 0: T 
1490/// 1: G 
1491/// 2: C 
1492/// 3: A 
1493#[derive(Debug, Clone, PartialEq, Eq, PartialOrd, Ord)]
1494pub struct SingleDirEdgeMap {
1495    edge_maps: [Box<[ID]>; ALPHABET_SIZE]
1496}
1497
1498impl SingleDirEdgeMap {
1499    pub fn new(edge_maps: [Box<[ID]>; ALPHABET_SIZE]) -> Self {
1500        SingleDirEdgeMap { edge_maps }
1501    }
1502
1503    pub fn complement(&self) -> Self {
1504        let mut reverse = self.edge_maps.clone();
1505        reverse.reverse();
1506        SingleDirEdgeMap::new(reverse)
1507    }
1508
1509    /// get the IDs mapped to the specified edge
1510    pub fn edge_map(&self, base: u8) -> &[ID] {
1511        &self.edge_maps[(ALPHABET_SIZE as u8 - 1 - base) as usize]
1512    }
1513}
1514
1515// TODO add methods
1516#[derive(Debug, Serialize, Deserialize)]
1517pub struct Label {
1518    group: char,
1519    sample_label: String
1520}
1521
1522/// category for the 
1523#[derive(Debug, Deserialize, Serialize, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, ValueEnum)]
1524#[serde(rename_all = "kebab-case")]
1525pub enum BaseQuality {
1526    NoCall,
1527    Marginal,
1528    Medium,
1529    High
1530}
1531
1532impl BaseQuality {
1533    fn from_u64(quality: u64) -> BaseQuality {
1534        match quality {
1535            0 => Self::NoCall,
1536            1 => Self::Marginal,
1537            2 => Self::Medium,
1538            3 => Self::High,
1539            _ => panic!("invalid base quality value")
1540        }
1541    }
1542
1543    fn as_char(&self) -> char {
1544        match self {
1545            Self::NoCall => '#',
1546            Self::Marginal => '-',
1547            Self::Medium => ';',
1548            Self::High => 'C',
1549        }
1550    }
1551}
1552
1553impl fmt::Display for BaseQuality {
1554    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
1555        match self {
1556            Self::NoCall => write!(f, "no-call"),
1557            Self::Marginal => write!(f, "marginal"),
1558            Self::Medium => write!(f, "medium"),
1559            Self::High => write!(f, "high"),
1560        }
1561    }
1562}
1563
1564#[derive(Debug, Deserialize, Serialize, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)]
1565pub struct QualityBins {
1566    marginal_top_cutoff: u8,
1567    high_bottom_cutoff: u8,
1568}
1569
1570impl Default for QualityBins  {
1571    fn default() -> Self {
1572        Self { marginal_top_cutoff: 15, high_bottom_cutoff: 30 }
1573    }
1574}
1575
1576impl QualityBins {
1577    pub fn new(marginal_top_cutoff: u8, high_bottom_cutoff: u8) -> QualityBins {
1578        Self { marginal_top_cutoff, high_bottom_cutoff }
1579    }
1580
1581    fn base_quality(&self, score: u8) -> BaseQuality {
1582        if score <= 2 {
1583            BaseQuality::NoCall
1584        } else if score < self.marginal_top_cutoff {
1585            BaseQuality::Marginal
1586        } else if score > self.high_bottom_cutoff {
1587            BaseQuality::High
1588        } else {
1589            BaseQuality::Medium
1590        }
1591    }
1592
1593    fn base_quality_from_ascii_bytes(&self, char: u8) -> BaseQuality {
1594        let score = char - 33;
1595        self.base_quality(score)
1596    }
1597}
1598
1599#[derive(Deserialize, Serialize, Clone, PartialEq, PartialOrd)]
1600pub struct QualityVec {
1601    storage: Vec<BaseQuality>
1602}
1603
1604impl QualityVec {
1605    fn from_vec(quality_vec: Vec<BaseQuality>) -> QualityVec {
1606        QualityVec { storage: quality_vec }
1607    }
1608
1609    fn from_ascii_bytes(quality_scores: &[u8], quality_bins: QualityBins) -> QualityVec {
1610        let mut vec = Vec::new();
1611        for score in quality_scores {
1612            vec.push(quality_bins.base_quality_from_ascii_bytes(*score));
1613        }
1614
1615        QualityVec { storage: vec }
1616    }
1617
1618    fn iter_k_lowest_q<K: Kmer>(&'_ self) -> KLowestQualityIter<'_, K> {
1619        KLowestQualityIter { 
1620            quality_vec: self, 
1621            start_pos: 0, 
1622            phantom_data: PhantomData,
1623        }
1624    }
1625
1626    fn len(&self) -> usize {
1627        self.storage.len()
1628    }
1629}
1630
1631impl fmt::Debug for QualityVec {
1632    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
1633        for q in self.storage.iter() {
1634            write!(f, "{}", q.as_char())?;
1635        }
1636
1637        Ok(())
1638    }
1639}
1640
1641pub struct KLowestQualityIter<'a, K: Kmer> {
1642    quality_vec: &'a QualityVec,
1643    start_pos: usize,
1644    phantom_data: PhantomData<K>
1645}
1646
1647impl<K: Kmer> Iterator for KLowestQualityIter<'_, K> {
1648    type Item = BaseQuality;
1649
1650    fn next(&mut self) -> Option<Self::Item> {
1651        let end_pos = self.start_pos + K::k();
1652        if end_pos <= self.quality_vec.len() {
1653            let range = self.start_pos..end_pos;
1654
1655            let quality = self.quality_vec.storage[range]
1656                .iter()
1657                .min()
1658                .expect("missing base quality");
1659
1660            self.start_pos += 1;
1661            Some(*quality)
1662        } else {
1663            None
1664        }
1665    }
1666}
1667
1668/// add the alignment buffer to a structure
1669/// - `size_heap`: contents of boxes/vectors -> are stored separately and do not go into alignment calculation
1670pub fn size_aligned(size_stack: usize, size_heap: usize, align: usize) -> usize {
1671    let empty = size_stack % align;
1672    let buffer = match empty {
1673        0 => 0,
1674        _ => align - empty
1675    };
1676
1677    buffer + size_heap + size_stack
1678}
1679
1680pub fn build_test_graph<K, SD, DI>() -> (SerReads<DI>, SerKmers<K, SD>, SerGraph<K, SD>)
1681where
1682    K: Kmer +  Send + Sync,
1683    SD: SummaryData<DI>,
1684    DI: ReadData
1685{
1686    /*
1687    transcrips
1688    GCAGCTAGCTAGCGCGACTACGATCGTAGCGCAGCGAGCAGGGGGGGGGATAGCTGTCGCGGGGACGTATTATTATTAAAATTGCGGCGCGAGCTATTCGAGCGGAGCGAGCGACAGGAGCGGAGTTTGCGGTACGGGATTTTCGGATATCGGC
1689    GCGATTATTTTGCGGGGGATTTTCGGTAGCGACTGGGGGGGGGTATCGATCGTGACAGCTTTCGACTGGGAGCGCAGCTAGGCAGGACGCATTAATTATATATCATTATTTTTTTCTATAAAAAAAAAAGAGCTAGCGATCGACGCGATCGAC
1690    TATATTATCGGCTGAGCGAGCGGGGGCAGCTATATTACGCGATAAAGAGCCCCCCGAGGCGAGGCGGACTTACGTAGCGCAGGCACCATGACGAGCTAGCAGTCAGTCGTAGCGATCA
1691    GCTAGCTAGCTGACTACGATCGACGGGGAGCATTAATTAGAAAAAAGAGAGAGACAGCTTTCGACTGGGAGCGCAGCTAGGCAGGACGCATTACTATCTATTATTATATATCATTATTTGCGATTGGGGTGCTAGCATGCGT
1692    */
1693
1694    let raw_reads = [
1695        ["GCAGCTAGCTAGCGCGACTACGATCGTAGCGCAGCGAGCAGGGGGGGGGA", "gene1", "sample1", "CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC;;CCCCCCCCCCCCCCC"],
1696        ["CTAGCGCGACTACGATCGTAGCGCAGCGAGCAGGGGGGGGGATAGCTGTC", "gene1", "sample1", "CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC"],
1697        ["GCAGCGAGCAGGGGGGGGGATAGCTGTCGCGGGGACGTATTATTATTAAA", "gene1", "sample1", "CCCCCCCCCCCCCCCCCCCCCC-CCCCCCCCCCCCCCCCCCCCCCCCCCC"],
1698        ["CGGGGACGTATTATTATTAAAATTGCGGCGCGAGCTATTCGAGCGGAGCG", "gene1", "sample1", "CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC"],
1699        ["TATTCGAGCGGAGCGAGCGACAGGAGAGGAGTTTGCGGTACGGGATTTTC", "gene1", "sample1", "CCCCCCCCCCCCCCCCCCCCCCCCCC--CCCCCCCCCCCCCCCCCCCCCC"],
1700        ["CGGAGCGAGCGACAGGAGCGGAGTTTGCGGTACGGGATTTTCGGATATCG", "gene1", "sample1", "CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC"],
1701        ["GAGCGAGCGACAGGAGCGGAGTTTGCGGTACGGGATTTTCGGATATCGGC", "gene1", "sample1", "CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC"],
1702        ["GCAGCTAGCTAGCGCGACTACGATCGTAGCGCAGCGAGCAGGGGGGGGGA", "gene1", "sample1", "CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC-CCCCCCCCCCCCCC"],
1703        ["TATTATTAAAATTGCGGCGCGAGCTATTCGAGCGGAGCGAGCGACAGGAG", "gene1", "sample1", "CCCCCCCCCCCCCCCCCC-CC--CCCCCCCCCCCCCCCCCCCCCCCCCCC"],
1704        ["ACTACGATCGTAGCGCAGCGAGCAGGGGGGGGGATAGCTGTCGCGGGGAC", "gene1", "sample1", "CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC"],
1705        ["AGCGGAGCGAGCGACAGGAGCGGAGTTTGCGGTACGGGATTTTCGGATAT", "gene1", "sample1", "CCCCCCCCCCCCCCCCCCCCC--CCCC--CCCCCCCCCCCCCCCCCCCCC"],
1706        ["TAGCTAGCGCGACTACGATCGTAGCGCAGCGAGCAGGGGGGGGGATAGCT", "gene1", "sample1", "CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC"],
1707        ["TACGATCGTAGCGCAGCGAGCAGGGGGGGGGATAGCTGTCGCGGGGACGT", "gene1", "sample1", "CCCCCCCCCCCCCCCCC--CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC"],
1708        ["ATTGCGGCGCGAGCTATTCGAGCGGAGCGAGCGACAGGAGCGGAGTTTGC", "gene1", "sample1", "CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC"],
1709        ["CAGCTAGCTAGCGCGACTACGATCGTAGCGCAGCGAGCAGGGGGGGGGAT", "gene1", "sample1", "CCCCCCCCCCCCCCCCCCCCCCCC-CCCCCCCCCCCCCCCCCCCCCCCCC"],
1710        ["GCGATTATTTTGCGGGGGATTTTCGGTAGCGACTGGGGGGGGGTATCGAT", "gene2", "sample2", "CCCCCCCC---CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC"],
1711        ["GGGGATTTTCGGTAGCGACTGGGGGGGGGTATCGATCGTGACAGCTTTCG", "gene2", "sample2", "CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC;CCCCCCCCCCCC"],
1712        ["TTTCGACTGGGAGCGCAGCTAGGCAGGACGCATTACTATCTATTATTATA", "gene2", "sample2", "CCCCCCCCCCCCCCCCCCCCCCCCCCCCCC.-CCCCCCCCCCCCCCCCCC"],
1713        ["CAGGACGCATTACTATCTATTATTATATATCATTATTTTTTTCTATAAAA", "gene2", "sample2", "CCCCCCCCCCCCCCCC---CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC"],
1714        ["TTTCGGTAGCGACTGGGGGGGGGTATCGATCGTGACAGCTTTCGACTGGG", "gene2", "sample2", "CCCCCCCCCCCCCC-CCCC-CCCCCCCCCCCCCCCCCCCCCCCCCCCCCC"],
1715        ["AGCTAGGCAGGACGCATTACTATCTATTATTATATATCATTATTTTTTTC", "gene2", "sample2", "CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC"],
1716        ["AGGACGCATTACTATCTATTATTATATATCATTATTTTTTTCTATAAAAA", "gene2", "sample2", "CCCCCCCCCCCCCCCCCCCCCCCCCCCCC--;CCCCCCCCCCCCCCCCCC"],
1717        ["GGGATTTTCGGTAGCGACTGGGGGGGGGTATCGATCGTGACAGCTTTCGA", "gene2", "sample2", "CCCCCCCCCCCC---CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC"],
1718        ["CTTTCGACTGGGAGCGCAGCTAGGCAGGACGCATTACTATCTATTATTAT", "gene2", "sample2", "CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC"],
1719        ["CGACTGGGGGGGGGTATCGATCGTGACAGCTTTCGACTGGGAGCGCAGCT", "gene2", "sample2", "CCCCCCCCCCCCCCCCCCCCCCCCCCCC-CCCCCCCCCCCCCCCCCCCCC"],
1720        ["ATTATATATCATTATTTTTTTCTATAAAAAAAAAAGAGCTAGCGATCGAC", "gene2", "sample2", "CCCCCCCCCCCCCCCCCCC-CCCCCCCCCCCCCCCCCCCCCCCCCCCCCC"],
1721        ["TAATTATATATCATTATTTTTTTCTATAAAAAAAAAAGAGCTAGCGATCG", "gene2", "sample2", "CCCCCC-CCCCCCCCCCCCCCCCCCCCCCCC--CCCCCCCCCCCCCCCCC"],
1722        ["ATCATTATTTTTTTCTATAAAAAAAAAAGAGCTAGCGATCGACGCGATCG", "gene2", "sample2", "CCCCCCCCCCCCCCCC;CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC"],
1723        ["CGTGACAGCTTTCGACTGGGAGCGCAGCTAGGCAGGACGCATTAATTATA", "gene2", "sample2", "CCCCCCCCCCCCCCCCCC;;;CCCCCCCCCCCCCCCCCCCCCCCCCCCCC"],
1724        ["GGGGGTATCGATCGTGACAGCTTTCGACTGGGAGCGCAGCTAGGCAGGAC", "gene2", "sample2", "CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC"],
1725        ["TATATTATCGGCTGAGCGAGCGGGGGGAGCTATATTACGCGATAAAGAGC", "gene3", "sample3", "CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC"],
1726        ["AGCGAGCGGGGGCAGCTATATTACGCGATAAAGAGCCCCCCGAGGCGAGG", "gene3", "sample3", "CCCCCCCCCCCCCCCCCCCCCCCCC---CCCCCCCCCCCCCCCCCCCCCC"],
1727        ["CTATATTACGCGATAAAGAGCCCCCCGAGGCGAGGCGGACTTACGTAGCG", "gene3", "sample3", "CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC"],
1728        ["CGCGATAAAGAGCCCCCCGAGGCGAGGCGGACTTACGTAGCGCAGGCACC", "gene3", "sample3", "CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC"],
1729        ["CTTACGTAGCGCAGGCACCATGACGAGCTAGCAGTCAGTCGTAGCGATCA", "gene3", "sample3", "CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC"],
1730        ["GCGATAAAGAGCCCCCCGAGGCGAGGCGGACTTACGTAGCGCAGGCACCA", "gene3", "sample3", "CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC"],
1731        ["TATATTACGCGATAAAGAGCCCCCCGAGGCGAGGCGGACTTACGTAGCGC", "gene3", "sample3", "CCCCCCCCCCCCCCCCCC-----CCCCCCCCCCCCCCCCCCCCCCCCCCC"],
1732        ["GAGCGGGGGCAGCTATATTACGCGATAAAGAGCCCCCCGAGGCGAGGCGG", "gene3", "sample3", "CCCCCCCCCCCCCCCCCCCCCCC;;;CCCCC;;CCCCCCCCCCCCCCCCC"],
1733        ["AGAGCCCCCCGAGGCGAGGCGGACTTACGTAGCGCAGGCACCATGACGAG", "gene3", "sample3", "CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC"],
1734        ["GGCGGACTTACGTAGCGCAGGCACCATGACGAGCTAGCAGTCAGTCGTAG", "gene3", "sample3", "CCCCCCCCCCCCC---CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC"],
1735        ["ACTTACGTAGCGCAGGCACCATGACGAGCTAGCAGTCAGTCGTAGCGATC", "gene3", "sample3", "CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC"],
1736        ["GGGGCAGCTATATTACGCGATAAAGAGCCCCCCGAGGCGAGGCGGACTTA", "gene3", "sample3", "CCCCCCCCCCCCCCCCCCCCCCCCCCCC---CC-CCCCCCCCCCCCCCCC"],
1737        ["CCGAGGCGAGGCGGACTTACGTAGCGCAGGCACCATGACGAGCTAGCAGT", "gene3", "sample3", "CCCCCCCCCCCC-CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC"],
1738        ["GGCGAGGCGGACTTACGTAGCGCAGGCACCATGACGAGCTAGCAGTCAGT", "gene3", "sample3", "CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC"],
1739        ["GAGCCCCCCGAGGCGAGGCGGACTTACGTAGCGCAGGCACCATGACGAGC", "gene3", "sample3", "CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC"],
1740        ["GCTAGCTAGCTGACTACGATCGACGGGGAGCATTAATTAGAAAAAAGAGA", "gene4", "sample4", "CCCCCCCCCCCCCCCCCC--CCCCCCCCCCCCCCCCCCCCCCCCCCCCCC"],
1741        ["ACTATCTATTATTATATATCATTATTTGCGATTGGGGTGCTAGCATGCGT", "gene4", "sample4", "CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC"],
1742        ["AGGCAGGACGCATTACTATCTATTATTATATATCATTATTTGCGATTGGG", "gene4", "sample4", "CCCCCCCCCCCCCCCCCCCCCCCCCCC-CCCCCCCCCCCCCCCCCCCCCC"],
1743        ["ATCGACGGGGAGCATTAATTAGAAAAAAGAGAGAGACAGCTTTCGACTGG", "gene4", "sample4", "CCCCCCCCCCCCCCCC-CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC"],
1744        ["AGGACGCATTACTATCTATTATTATATATCATTATTTGCGATTGGGGTGC", "gene4", "sample4", "CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC"],
1745        ["GAGCATTAATTAGAAAAAAGAGAGAGACAGCTTTCGACTGGGAGCGCAGC", "gene4", "sample4", "CCCCCCCCCCCCCCCCCCCCCCCCCCC-CCCCCCCCCCCCCCCCCCCCCC"],
1746        ["ACTGGGAGCGCAGCTAGGCAGGACGCATTACTATCTATTATTATATATCA", "gene4", "sample4", "CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC"],
1747        ["GATCGACGGGGAGCATTAATTAGAAAAAAGAGAGAGACAGCTTTCGACTG", "gene4", "sample4", "CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC"],
1748        ["CGCAGCTAGGCAGGACGCATTACTATCTATTATTATATATCATTATTTGC", "gene4", "sample4", "CCCCCCCCCC-CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC"],
1749        ["CGATCGACGGGGAGCATTAATTAGAAAAAAGAGAGAGACAGCTTTCGACT", "gene4", "sample4", "CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC"],
1750        ["ATCGACGGGGAGCATTAATTAGAAAAAAGAGAGAGACAGCTTTCGACTGG", "gene4", "sample4", "CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC--CCCCCCCCCCCCCC"],
1751        ["GCAGCTAGGCAGGACGCATTACTATCTATTATTATATATCATTATTTGCG", "gene4", "sample4", "CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC"],
1752        ["GACGGGGAGCATTAATTAGAAAAAAGAGAGAGACAGCTTTCGACTGGGAG", "gene4", "sample4", "CCCC----CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC"],
1753        ["CGATCGACGGGGAGCATTAATTAGAAAAAAGAGAGAGACAGCTTTCGACT", "gene4", "sample4", "CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC---CCCCCCCCCCCCCC"],
1754        ["TAGGCAGGACGCATTACTATCTATTATTATATATCATTATTTGCGATTGG", "gene4", "sample4", "CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC-CCCCCCCCCCCCCCCCC"]
1755    ];
1756
1757    let mut reads = Reads::new_with_quality(crate::reads::Strandedness::Forward);
1758    let mut id_translator = BiMap::new();
1759    let mut tag_translator = BiMap::new();
1760
1761    for [seq, gene, sample, quality] in raw_reads {
1762        let gene = String::from(gene);
1763        let sample = String::from(sample);
1764
1765        let id = match id_translator.get_by_left(&gene) {
1766            Some(id) => *id,
1767            None =>  {
1768                let new_id = id_translator.len() as ID;
1769                id_translator.insert(gene, new_id);
1770                new_id
1771            }
1772        };
1773
1774        let tag = match tag_translator.get_by_left(&sample) {
1775            Some(tag) => *tag,
1776            None =>  {
1777                let new_tag = tag_translator.len() as Tag;
1778                tag_translator.insert(sample, new_tag);
1779                new_tag
1780            }
1781        };
1782
1783        reads.add_read(DnaString::from_acgt_bytes(seq.as_bytes()), None, DI::new(id, tag), Some(quality.as_bytes()));
1784    }
1785
1786    let translator = Translator::new(id_translator, tag_translator);
1787    let reads_paired = ReadsPaired::Unpaired { reads };
1788    let sample_kmers = reads_paired.tag_kmers_vec(K::k(), 4);
1789    let ser_reads = SerReads::new(reads_paired, translator.clone());
1790    
1791    let sample_info = SampleInfo::new(0b1100, 0b0011, sample_kmers);
1792    let summary_config = SummaryConfig::new(sample_info);
1793
1794    let (kmers, _) = filter_kmers::<SD, K, _>(
1795        ser_reads.reads(), 
1796        &summary_config, 
1797        false, 
1798        5., 
1799        false
1800    );
1801
1802    let ser_kmers = SerKmers::new(kmers.clone(), translator.clone(), summary_config.clone());
1803
1804    let comp_spec = CheckCompress::new(|d: SD, _| d, |d, d1| d.join_test(d1));
1805    let graph = compress_kmers_with_hash(true, &comp_spec, kmers, false, false).finish();
1806    let ser_graph = SerGraph::new(graph, translator, summary_config);
1807
1808    (ser_reads, ser_kmers, ser_graph)
1809}
1810
1811#[cfg(test)]
1812mod tests {
1813    use bimap::BiMap;
1814
1815    use crate::{ALPHABET_SIZE, BaseQuality, Dir, EdgeMap, EdgeMult, Exts, Kmer, QualityBins, Tags, TagsCountsFormatter, TagsFormatter, kmer::{Kmer4, Kmer17}, size_aligned, summarizer::{ID, Marker, Tag, Translator}};
1816
1817    #[test]
1818    fn test_dir_index() {
1819        let bases = [0, 1, 2, 3];
1820
1821        for (i, (dir, base)) in [Dir::Right, Dir::Left].iter().flat_map(|d| bases.into_iter().rev().map(move |b| (d, b))).enumerate() {
1822            assert_eq!(i as u8, dir.index(base))
1823        }
1824
1825        assert_eq!((Dir::Left).index_range(), 4..8);
1826        assert_eq!((Dir::Right).index_range(), 0..4);
1827    }
1828
1829    #[test]
1830    fn test_remove_ext() {
1831        let ext = Exts::new(0b11111011);
1832        assert_eq!(ext.remove(Dir::Right, 0).val, 0b11101011);
1833        assert_eq!(ext.remove(Dir::Left, 1).val, 0b11111001);
1834    }
1835
1836    #[test]
1837    fn test_edge_mult() {
1838        let mut edge_mult = EdgeMult::new();
1839        assert_eq!(edge_mult.edge_mults, [0; 2*ALPHABET_SIZE]);
1840
1841        let exts = Exts::new(0b11111111);
1842        edge_mult.add_exts(exts);
1843        assert_eq!(edge_mult.edge_mults, [1, 1, 1, 1, 1, 1, 1, 1]);
1844
1845        edge_mult.add(0, crate::Dir::Right, 2);
1846        edge_mult.add(2, crate::Dir::Left, 78989);
1847        assert_eq!(edge_mult.edge_mults, [1, 1, 1, 3, 1, 78990, 1, 1]);
1848
1849        let exts = Exts::new(0);
1850        let comp = [1, 1, 1, 3, 1, 78990, 1, 1];
1851        edge_mult.add_exts(exts);
1852        assert_eq!(edge_mult.edge_mults, comp);
1853        assert_eq!(edge_mult.sum(), comp.iter().sum::<u32>());
1854
1855        let exts = Exts::new(0b10101010);
1856        edge_mult.add_exts(exts);
1857        assert_eq!(edge_mult.edge_mults, [2, 1, 2, 3, 2, 78990, 2, 1]);
1858
1859        let exts = Exts::new(0b01010101);
1860        edge_mult.add_exts(exts);
1861        assert_eq!(edge_mult.edge_mults, [2, 2, 2, 4, 2, 78991, 2, 2]);
1862
1863        let clean_exts = Exts::new(0b10010101);
1864        edge_mult.clean_edges(clean_exts);
1865        assert_eq!(edge_mult.edge_mults, [2, 0, 0, 4, 0, 78991, 0, 2]);
1866
1867        let mut em = EdgeMult::new();
1868        let exts = Exts::new(0b00011101);
1869        println!("{:?}", exts);
1870        em.add_exts(exts);
1871        println!("{}", em);
1872        println!("{:?}", em);
1873        assert_eq!(em.exts(), exts);
1874
1875        assert_eq!(em.left(), &[1, 1, 0, 1]);
1876        assert_eq!(em.right(), &[0, 0, 0, 1]);
1877
1878        // single dir em
1879        let sdir_em = em.single_dir(Dir::Left);
1880        assert_eq!(sdir_em.edge_mults, [1, 1, 0, 1]);
1881        let sdir_em = em.single_dir(Dir::Right);
1882        assert_eq!(sdir_em.edge_mults, [0, 0, 0, 1]);
1883        let em_2 = EdgeMult::from_single_dirs(&Some(em.single_dir(Dir::Left)), &Some(em.single_dir(Dir::Right)));
1884        assert_eq!(em_2.unwrap(), em);
1885        assert_eq!(sdir_em.complement().edge_mults, [1, 0, 0, 0]);
1886
1887        // reverse complement
1888        em.rc();
1889        assert_eq!(em.edge_mults, [1, 0, 1, 1, 1, 0, 0, 0]);
1890    }
1891
1892    #[test]
1893    fn test_edge_map() {
1894        let empty_emaps: [Box<[ID]>; 8] = Default::default();
1895        let default_emap = EdgeMap::default();
1896
1897        let mut emap = EdgeMap::new(empty_emaps.clone());
1898
1899        assert!(emap.is_empty());
1900
1901        assert_eq!(emap.edge_maps, empty_emaps);
1902        assert_eq!(emap, default_emap);
1903
1904        let m = vec![1, 2];
1905        emap.set_edge_map_at_index(m.clone().into(), 1); // right G
1906        let r = emap.edge_map(2, Dir::Right);
1907
1908        assert!(!emap.is_empty());
1909
1910        assert_eq!(&m, r);
1911
1912        emap.add_id_to_edge_map_at_index(3, 1);
1913        let r = emap.edge_map(2, Dir::Right);
1914        assert_eq!(&[1, 2, 3], r);
1915
1916        let exp_mem = 3*std::mem::size_of::<ID>();
1917        let r_mem = emap.mem_heap();
1918        assert_eq!(exp_mem, r_mem);
1919
1920        emap.add_id(Exts::new(0b01000010), 4); // right G and left C
1921        assert_eq!(&[1, 2, 3, 4], emap.edge_map(2, Dir::Right));
1922        assert_eq!(&[4], emap.edge_map(1, Dir::Left));
1923
1924        let left = emap.single_dir(Dir::Left);
1925        let right = emap.single_dir(Dir::Right);
1926        let new_emap = EdgeMap::from_single_dirs(&Some(left), &Some(right)).unwrap();
1927        assert_eq!(emap, new_emap);
1928
1929        let mut clean_emap = EdgeMap::default();
1930        clean_emap.add_id_to_edge_map_at_index(4, 6); // left C
1931
1932        emap.clean_edges(Exts::new(0b00000010));
1933
1934        assert_eq!(emap, clean_emap);
1935
1936        assert_eq!("A: [] | []\nC: [4] | []\nG: [] | []\nT: [] | []\n", &format!("{}", emap));
1937        assert_eq!("A: [], C: [4], G: [], T: [] | A: [], C: [], G: [], T: []", &format!("{:?}", emap));
1938
1939    }
1940
1941    #[test]
1942    fn test_base_quality() {
1943        let bq = BaseQuality::from_u64(0);
1944        assert_eq!(&format!("{bq} {}", bq.as_char()), "no-call #");
1945
1946        let bq = BaseQuality::from_u64(1);
1947        assert_eq!(&format!("{bq} {}", bq.as_char()), "marginal -");
1948
1949        let bq = BaseQuality::from_u64(2);
1950        assert_eq!(&format!("{bq} {}", bq.as_char()), "medium ;");
1951
1952        let bq = BaseQuality::from_u64(3);
1953        assert_eq!(&format!("{bq} {}", bq.as_char()), "high C");
1954
1955        let bins = QualityBins::new(15, 30);
1956        assert_eq!(bins, QualityBins::default());
1957    }
1958
1959    #[test]
1960    #[should_panic]
1961    fn test_base_quality_panic() {
1962        let _bq = BaseQuality::from_u64(5);
1963    }
1964
1965    #[test]
1966    fn test_bit_and_dist() {
1967        let marker: Marker = 0b1111000011110000111100001111000011110000111100001111000011110000;
1968        println!("marker:   {:064b}", marker);
1969
1970        let tags = Tags::from_tag_vec(&vec![0, 1, 4]);
1971        println!("tags:     {:064b}", tags.val);
1972        let dist = tags.bit_and_dist(marker);
1973        println!("dist: {}", dist);
1974        assert_eq!(tags.len(), 3);
1975
1976        let tags = Tags::from_tag_vec(&vec![1, 5, 19, 25, 32]);
1977        println!("tags:     {:064b}", tags.val);
1978        let dist = tags.bit_and_dist(marker);
1979        println!("dist: {}", dist);
1980        assert_eq!(tags.len(), 5);
1981
1982        let tags = Tags::from_tag_vec(&vec![0, 1, 2, 3, 4, 5, 6, 7, 63]);
1983        println!("tags:     {:064b}", tags.val);
1984        let dist = tags.bit_and_dist(marker);
1985        println!("dist: {}", dist);
1986        assert_eq!(tags.len(), 9);
1987
1988        let tags = Tags::from_tag_vec(&vec![31]);
1989        println!("tags:     {:064b}", tags.val);
1990        let dist = tags.bit_and_dist(marker);
1991        println!("dist: {}", dist);
1992        assert_eq!(tags.len(), 1);
1993
1994        let tags = Tags::from_tag_vec(&vec![63]);
1995        println!("tags:     {:064b}", tags.val);
1996        let dist = tags.bit_and_dist(marker);
1997        println!("dist: {}", dist);
1998        assert_eq!(tags.len(), 1);
1999    }
2000
2001    #[test]
2002    fn test_tag_formatter() {
2003        let mut tag_translator = BiMap::new();
2004        let samples = vec!["A", "B", "C", "D", "E", "F", "G"];
2005
2006        for (i, label) in samples.into_iter().enumerate() {
2007            tag_translator.insert(label.to_string(), i as Tag);
2008        }
2009
2010        let translator = Translator::new_tag_translator(tag_translator);
2011
2012        let tags = Tags::from_tag_vec(&vec![0, 1, 4]);
2013        let counts = vec![1, 2, 3].into_boxed_slice();
2014        print!("{}", TagsCountsFormatter::new(tags, &counts, &translator));
2015
2016        let tags = Tags::from_tag_vec(&vec![0, 1, 4, 6]);
2017        let counts = vec![1, 2, 3, 0].into_boxed_slice();
2018        print!("{}", TagsCountsFormatter::new(tags, &counts, &translator));
2019
2020        let tags = Tags::from_tag_vec(&vec![0, 1, 4]);
2021        print!("{}", TagsFormatter::new(tags, &translator));
2022
2023        let tags = Tags::from_tag_vec(&vec![0, 1, 4, 6]);
2024        print!("{}", TagsFormatter::new(tags, &translator));
2025
2026    }
2027
2028    #[test]
2029    fn test_iter_tags() {
2030        let tags = Tags::from_tag_vec(&vec![0, 1, 4, 12, 32, 63]);
2031        for tag in tags.iter() {
2032            println!("tag: {tag}")
2033        }
2034    }
2035
2036    #[test]
2037    fn test_kmer_complexity() {
2038        let kmer = Kmer4::empty();
2039        assert!(kmer.has_low_complexity());
2040        let kmer = Kmer4::from_ascii("TTTTTTTT".as_bytes());
2041        assert!(kmer.has_low_complexity());
2042        let kmer = Kmer4::from_ascii("CCCCCCCC".as_bytes());
2043        assert!(kmer.has_low_complexity());
2044        let kmer = Kmer4::from_ascii("GGGGGGGG".as_bytes());
2045        assert!(kmer.has_low_complexity());
2046
2047        let kmer = Kmer4::from_ascii("ACGATCGA".as_bytes());
2048        assert!(!kmer.has_low_complexity());
2049        let kmer = Kmer4::from_ascii("AGCAGCTC".as_bytes());
2050        assert!(!kmer.has_low_complexity());
2051
2052        let kmer = Kmer17::empty();
2053        assert!(kmer.has_low_complexity());
2054        let kmer = Kmer17::from_ascii("TTTTTTTTTTTTTTTTT".as_bytes());
2055        assert!(kmer.has_low_complexity());
2056        let kmer = Kmer17::from_ascii("CCCCCCCCCCCCCCCCC".as_bytes());
2057        assert!(kmer.has_low_complexity());
2058        let kmer = Kmer17::from_ascii("GGGGGGGGGGGGGGGGG".as_bytes());
2059        assert!(kmer.has_low_complexity());
2060
2061        let kmer = Kmer17::from_ascii("ACGATCGAGACTGACTG".as_bytes());
2062        assert!(!kmer.has_low_complexity());
2063        let kmer = Kmer17::from_ascii("AGCAGCTCAGCTAGCTG".as_bytes());
2064        assert!(!kmer.has_low_complexity());
2065    }
2066}
2067
2068
debruijn/lib.rs

debruijn/
lib.rs
