debruijn/
reads.rs

1use std::collections::HashMap;
2use std::mem::take;
3use std::ops::Range;
4use bimap::BiMap;
5use itertools::Itertools;
6use serde::de::DeserializeOwned;
7use serde_derive::{Deserialize, Serialize};
8use std::fmt::{Debug, Display};
9use std::hash::Hash;
10use std::{mem, str};
11use crate::dna_string::DnaString;
12use crate::summarizer::{IDTag, Tag, ID};
13use crate::{BaseQuality, Exts, Kmer, QualityBins, QualityVec, Vmer, base_to_bits, base_to_bits_checked};
14
15#[derive(Debug, Serialize, Deserialize, PartialEq, Eq, PartialOrd, Ord, Clone, Hash, Copy)]
16pub enum Strandedness {
17    Forward,
18    Reverse,
19    Unstranded
20}
21
22/// a sequencing read and additional information
23#[derive(Debug, PartialEq, )]
24pub struct Read<D> {
25    seq: DnaString,
26    exts: Exts,
27    data: D,
28    strand: Strandedness,
29    quality: Option<QualityVec>
30}
31
32impl<D: Clone + Copy> Read<D> {
33    /// a new `Read`
34    pub fn new(seq: DnaString, exts: Exts, data: D, strand: Strandedness, quality: Option<QualityVec>) -> Read<D> {
35        Read { seq, exts, data, strand, quality }
36    }
37
38    /// the sequence of the `Read`
39    pub fn seq(&self) -> &DnaString {
40        &self.seq
41    }
42
43    /// the [`Exts`] of the `Read`
44    pub fn exts(&self) -> Exts {
45        self.exts
46    }
47
48    /// the data of the `Read`
49    pub fn data(&self) -> D {
50        self.data
51    }
52
53    /// the strandedness of the `Read`
54    pub fn stranded(&self) -> Strandedness {
55        self.strand
56    }
57
58    pub fn iter_kmer_exts_quality<'a, K: Kmer + 'a>(&'a self) -> Box<dyn Iterator<Item = (K, Exts, Option<BaseQuality>)> + 'a> {
59        if let Some(quality) = self.quality.as_ref() {
60            Box::new(self.seq()
61                .iter_kmer_exts::<K>(self.exts)
62                .zip(quality
63                    .iter_k_lowest_q::<K>()
64                    .map(Some)
65                )
66                .map(|((kmer, exts), quality)| (kmer, exts, quality))
67            )
68        } else {
69            Box::new(self.seq()
70                .iter_kmer_exts::<K>(self.exts)
71                .map(|(kmer, exts)| (kmer, exts, None)))
72        }
73    }
74}
75
76/// two paired sequencing reads
77pub struct PairedRead<D> {
78    read1: Read<D>,
79    read2: Read<D>
80}
81
82impl<D> PairedRead<D> {
83    /// a new `PairedRead`
84    pub fn new(read1: Read<D>, read2: Read<D>) -> PairedRead<D> {
85        PairedRead {
86            read1,
87            read2
88        }
89    }
90
91    /// the paired reads
92    pub fn reads(&self) -> (&Read<D>, &Read<D>) {
93        (&self.read1, &self.read2)
94    }
95}
96
97/// Store many DNA sequences together with an Exts and data each compactly packed together
98/// 
99/// #### fields:
100/// 
101/// * `storage`: `Vec` with 2-bit encoded DNA bases of all sequences
102/// * `ends`:  `Vec` with the ends (exclusive) of the separate sequences in the `Reads`
103/// * `exts`: `Option<Vec>` with one Exts for each sequence
104/// * `data`: `Vec` with data for each sequence
105/// * `len`: length of all sequences together
106/// * `stranded`: [`Stranded`] conveying the strandedness and direction of the reads
107#[derive(Ord, PartialOrd, Clone, PartialEq, Eq, Hash, Serialize, Deserialize, Debug)]
108pub struct Reads<D> {
109    storage: Vec<u64>,
110    ends: Vec<usize>,
111    exts: Option<Vec<Exts>>,
112    quality: Option<Vec<u64>>,
113    data: Vec<D>,
114    len: usize,
115    stranded: Strandedness,
116    quality_bins: QualityBins
117
118}
119
120impl<D: Clone + Copy> Reads<D> {
121
122    /// Returns a new `Reads`
123    pub fn new(stranded: Strandedness) -> Self {
124        Reads {
125            storage: Vec::new(),
126            ends: Vec::new(),
127            exts: None,
128            data: Vec::new(),
129            quality: None,
130            len: 0,
131            stranded,
132            quality_bins: QualityBins::default(),
133        }
134    }
135
136    pub fn new_with_quality(stranded: Strandedness) -> Self {
137        Reads {
138            storage: Vec::new(),
139            ends: Vec::new(),
140            exts: None,
141            data: Vec::new(),
142            quality: Some(Vec::new()),
143            len: 0,
144            stranded,
145            quality_bins: QualityBins::default(),
146        }
147    }
148
149    #[inline(always)]
150    pub fn stranded(&self) -> Strandedness {
151        self.stranded
152    }
153
154    #[inline(always)]
155    /// Returns the number of reads stored
156    pub fn n_reads(&self) -> usize {
157        self.ends.len()
158    }
159
160    /// get the memory required for the reads
161    pub fn mem(&self) -> usize {
162        let exts_size = if let Some(e_vec) = self.exts.as_ref() { size_of_val(&**e_vec) } else { 0 };
163        let quality_size_inner = if let Some(q) = self.quality.as_ref() {
164            size_of_val(&*q)
165        } else {
166            0
167        };
168        mem::size_of_val(self) + size_of_val(&*self.storage) + size_of_val(&*self.data) + size_of_val(&*self.ends) + exts_size + quality_size_inner
169    }
170
171    /// set the strandedness and the direction of the reads
172    pub fn set_stranded(&mut self, stranded: Strandedness) {
173        self.stranded = stranded
174    }
175
176    /// set custom quality bins
177    pub fn set_custom_quality_bins(&mut self, quality_bins: QualityBins) {
178        self.quality_bins = quality_bins
179    }
180
181    /// add exts to `Reads` if needed - use after adding sequence and new end
182    fn add_exts(&mut self, exts: Option<Exts>) {
183        match exts {
184            Some(e) => {
185                // only add exts if some exts are not empty
186                match self.exts.as_mut() {
187                    // already has exts, simply append new exts
188                    Some(e_vec) => e_vec.push(e),
189                    // no exts so far
190                    None => {
191                        // check if exts are empty
192                        if e != Exts::empty() {
193                            // if not, add vector of empty exts and then push new exts
194                            self.exts = Some(vec![Exts::empty(); self.n_reads() - 1]);
195                            self.exts.as_mut().unwrap().push(e);
196                        } // else keep no exts
197                    }
198                }
199            }
200            None => {
201                if let Some(e_vec) = self.exts.as_mut() {
202                    // Reads has exts, but no exts here, so add empty exts
203                    e_vec.push(Exts::empty())
204                } // else do nothing
205            }
206        }
207        
208    }
209
210    /// Adds a new read to the `Reads`
211    // maybe push_base until u64 is full and then do extend like in DnaString::extend ? with accellerated mode
212    pub fn add_read<V: Vmer>(&mut self, seq: V, exts: Option<Exts>, data: D, quality_scores: Option<&[u8]>) {
213        // check if we have quality scores
214        if let (Some(q), true) = (quality_scores, self.quality.is_some()) {
215            // we do, add bases and quality
216            assert_eq!(seq.len(), q.len(), "mismatch in read and quality score length");
217            for (base, &score) in seq.iter().zip(q) {
218                self.push_base_and_quality(base, score)
219            }
220        } else if self.quality.is_none() {
221            // we do not, only add base
222            for base in seq.iter() {
223                self.push_base(base);
224            }
225        } else {
226            // mismatch in quality scores, panic
227            panic!("Error: quality scores have to be added to all reads or none")
228        }
229
230        self.ends.push(self.len);
231        self.add_exts(exts);
232        self.data.push(data);
233    }
234
235    /// Transforms a `[(vmer, exts, data)]` into a `Reads` - watch for memory usage
236    // TODO test if memory efficient
237    pub fn from_vmer_vec<V: Vmer, S: IntoIterator<Item=(V, Exts, D)>>(vec_iter: S, stranded: Strandedness) -> Self {
238        let mut reads = Reads::new(stranded);
239        for (vmer, exts, data) in vec_iter {
240            for base in vmer.iter() {
241                reads.push_base(base);
242            }
243            reads.ends.push(reads.len);
244            reads.add_exts(Some(exts));
245            reads.data.push(data);
246        }
247
248        reads.shrink_to_fit();
249        
250        reads
251    }
252
253
254    /// add ASCII encoded bases to the `Reads`
255    /// 
256    /// will transform all ascii characters outside of ACGTacgt into A
257    /// 
258    /// if `Reads` previously contained no exts, no new exts will be added
259    /// see also: [`Reads::add_from_bytes_checked`]
260    pub fn add_from_bytes(&mut self, bytes: &[u8], exts: Option<Exts>, data: D) {
261        
262        // fill the last incomplete u64 block
263        let missing = 32 - (self.len % 32);
264        if missing != 0 {
265            if  missing > bytes.len() {
266                let fill = bytes.iter().map(|c| base_to_bits(*c));
267                self.extend(fill);
268                self.ends.push(self.len);
269                self.add_exts(exts);
270                self.data.push(data);
271                return;
272            } else {
273                let fill = bytes[0..missing].iter().map(|c| base_to_bits(*c));
274                self.extend(fill);
275            }
276        }
277        
278        // Accelerated avx2 mode. Should run on most machines made since 2013.
279        #[cfg(any(target_arch = "x86", target_arch = "x86_64"))]
280        {
281            if is_x86_feature_detected!("avx2") {
282                for chunk in bytes[missing..bytes.len()].chunks(32) {
283                    if chunk.len() == 32 {
284                        let (conv_chunk, _) = unsafe { crate::bitops_avx2::convert_bases(chunk) };
285                        let packed = unsafe { crate::bitops_avx2::pack_32_bases(conv_chunk) };
286                        self.storage.push(packed);
287                        self.len += 32;
288                    } else {
289                        let b = chunk.iter().map(|c| base_to_bits(*c));
290                        self.extend(b);
291                    }
292                }
293                self.ends.push(self.len);
294                self.add_exts(exts);
295                self.data.push(data);
296                
297                return;
298            }
299        }
300
301        let b = bytes.iter().map(|c| base_to_bits(*c));
302        self.extend(b);
303        self.ends.push(self.len);
304        self.add_exts(exts);
305        self.data.push(data);
306        
307        
308    }
309
310    /// add ASCII encoded bases to the Reads
311    /// 
312    /// will return `false` if the bytes contained characters outside of `ACGTacgt`, otherwise return true and add the bases
313    /// see also: [`Reads::add_from_bytes`]
314    pub fn add_from_bytes_checked(&mut self, bytes: &[u8], exts: Option<Exts>, data: D) -> bool {
315
316        let (_, corrects): (Vec<u8>, Vec<bool>) = bytes.iter().map(|c| base_to_bits_checked(*c)).collect();
317        if corrects.iter().contains(&false) { return false }
318
319        
320        // fill the last incomplete u64 block
321        let missing = 32 - (self.len % 32);
322        if missing != 0 {
323            if  missing > bytes.len() {
324                let fill = bytes.iter().map(|c| base_to_bits(*c));
325                self.extend(fill);
326                self.add_exts(exts);
327                self.data.push(data);
328                self.ends.push(self.len);
329                return true;
330            } else {
331                let fill = bytes[0..missing].iter().map(|c| base_to_bits(*c));
332                self.extend(fill);
333            }
334        }
335        
336        // Accelerated avx2 mode. Should run on most machines made since 2013.
337        #[cfg(any(target_arch = "x86", target_arch = "x86_64"))]
338        {
339            if is_x86_feature_detected!("avx2") {
340                for chunk in bytes[missing..bytes.len()].chunks(32) {
341                    if chunk.len() == 32 {
342                        let (conv_chunk, _) = unsafe { crate::bitops_avx2::convert_bases(chunk) };
343                        let packed = unsafe { crate::bitops_avx2::pack_32_bases(conv_chunk) };
344                        self.storage.push(packed);
345                        self.len += 32;
346                    } else {
347                        let b = chunk.iter().map(|c| base_to_bits(*c));
348                        self.extend(b);
349                    }
350                }
351
352                self.ends.push(self.len);
353                self.add_exts(exts);
354                self.data.push(data);
355                
356                return true;
357            }
358        }
359
360        let b = bytes.iter().map(|c| base_to_bits(*c));
361        self.extend(b);
362        self.ends.push(self.len);
363        self.add_exts(exts);
364        self.data.push(data);
365
366        true         
367    }
368
369
370    /// Add new 2-bit encoded base to the `Reads`
371    fn push_base(&mut self, base: u8) {
372        let bit = (self.len % 32) * 2;
373        if bit != 0 {
374            match self.storage.pop() {
375                Some(last) => {
376                    let last = last + ((base as u64) << (64 - bit - 2));
377                    self.storage.push(last);
378                },
379                None => panic!("tried to push base to empty vector (?)")
380            }
381        } else {
382            self.storage.push((base as u64) << 62);
383        }
384        self.len += 1; 
385    }
386
387    /// Simultaniously add new 2-bit encoded base and quality to the `Reads`
388    fn push_base_and_quality(&mut self, base: u8, score: u8) {
389        let Some(quality) = self.quality.as_mut() else { return; };
390        let base_quality = self.quality_bins.base_quality_from_ascii_bytes(score);
391
392        let bit = (self.len % 32) * 2;
393        if bit != 0 {
394            match self.storage.pop() {
395                Some(last) => {
396                    let last = last + ((base as u64) << (64 - bit - 2));
397                    self.storage.push(last);
398                },
399                None => panic!("tried to push base to empty vector (?)")
400            }
401
402            match quality.pop() {
403                Some(last) => {
404                    let last = last + ((base_quality as u64) << (64 - bit - 2));
405                    quality.push(last);
406                },
407                None => panic!("tried to push quality to empty vector (?)")
408            }
409        } else {
410            self.storage.push((base as u64) << 62);
411            quality.push((base_quality as u64) << 62);
412        }
413        self.len += 1; 
414    }
415
416    /// extend the reads' storage by 2-bit encoded bases
417    fn extend(&mut self, mut bytes: impl Iterator<Item = u8>) {
418        // fill the last incomplete u64 block
419        while self.len % 32 != 0 {
420            match bytes.next() {
421                Some(b) => self.push_base(b),
422                None => return,
423            }
424        }
425
426        let mut bytes = bytes.peekable();
427
428        // chunk the remaining items into groups of at most 32 and handle them together
429        while bytes.peek().is_some() {
430            let mut val: u64 = 0;
431            let mut offset = 62;
432            let mut n_added = 0;
433
434            for _ in 0..32 {
435                if let Some(b) = bytes.next() {
436                    assert!(b < 4);
437                    val |= (b as u64) << offset;
438                    offset -= 2;
439                    n_added += 1;
440                } else {
441                    break;
442                }
443            }
444
445            self.storage.push(val);
446            self.len += n_added;
447        }
448    }
449
450    #[inline(always)]
451    fn addr(&self, i: &usize) -> (usize, usize) {
452        (i / 32, (i % 32 ) * 2)
453    }
454
455    /// get the `i`th read in a `Reads`
456    pub fn get_read(&self, i: usize) -> Option<Read<D>> {
457        if i >= self.n_reads() { return None }
458
459        // get the read sequence
460        let mut sequence = DnaString::new();
461        let end = self.ends[i];
462        //let start = if i != 0 { self.ends[i-1] } else { 0 };
463        let start = match i {
464            0 => 0,
465            1.. => self.ends[i-1]
466        };
467
468        for b in start..end {
469            let (block, bit) = self.addr(&b);
470            let base = ((self.storage[block] >> (62 - bit)) & 3u64) as u8;
471            sequence.push(base);
472        }
473
474        // get the quality score sequence
475        let quality = if let Some(quality) = self.quality.as_ref() {
476            let mut base_qualities = Vec::new();
477            for q in start..end {
478                let (block, bit) = self.addr(&q);
479                let base_quality = BaseQuality::from_u64((quality[block] >> (62 - bit)) & 3u64);
480                base_qualities.push(base_quality);
481            }
482            
483            Some(QualityVec::from_vec(base_qualities))
484        } else {
485            None
486        };
487
488        let exts = match self.exts {
489            Some(ref e_vec) => e_vec[i],
490            None => Exts::empty()
491        };
492
493        Some(Read::new(sequence, exts, self.data[i], self.stranded, quality))
494    }
495
496
497    /// shrink the vectors' capacity to fit the length
498    /// 
499    /// use sparsely
500    pub fn shrink_to_fit(&mut self)  {
501        self.storage.shrink_to_fit();
502        self.data.shrink_to_fit();
503        if let Some(e_vec) = self.exts.as_mut() { e_vec.shrink_to_fit(); }
504        self.ends.shrink_to_fit();
505    }
506
507    /// Iterate over the reads as (DnaString, Exts, D).
508    pub fn iter(&self) -> ReadsIter<'_, D> {
509        ReadsIter {
510            reads: self,
511            i: 0,
512            end: self.n_reads(),
513            length: self.n_reads(),
514        }
515    }
516
517    /// Iterate over a range start reads as (DnaString, Exts, D).
518    pub fn partial_iter(&self, range: Range<usize>) -> ReadsIter<'_, D> {
519        assert!(range.end <= self.n_reads());
520        assert!(range.start < self.n_reads());
521        assert!(range.start < range.end);
522        ReadsIter {
523            reads: self,
524            i: range.start,
525            end: range.end,
526            length: (range.end - range.start)
527        }
528    }
529
530    pub fn info(&self) -> String {
531        format!("Reads {{ n reads: {}, stranded: {:?} }}", self.n_reads(), self.stranded)
532    }
533}
534
535impl<D: ReadData> Reads<D> {
536    /// get the number of k-mers for each unique data value
537    pub fn tag_kmers(&self, k: usize) -> HashMap<Tag, usize> {
538        let mut hm = HashMap::new();
539
540        self.iter().for_each(|read| {
541            let kmers = read.seq.len().saturating_sub(k - 1);
542            if let Some(tag) = read.data.get_tag() {
543                if let Some(count) = hm.get_mut(&tag) {
544                    *count += kmers;
545                } else {
546                    hm.insert(tag, kmers);
547                }
548            }
549            
550        });
551
552        hm
553    }
554}
555
556impl<D: Clone + Copy> Default for Reads<D> {
557    fn default() -> Self {
558        Self::new(Strandedness::Unstranded)
559    }
560}
561
562/// Iterator over values of a DnaStringoded sequence (values will be unpacked into bytes).
563pub struct ReadsIter<'a, D> {
564    reads: &'a Reads<D>,
565    i: usize,
566    end: usize,
567    length: usize,
568}
569
570impl<D: Clone + Copy> Iterator for ReadsIter<'_, D> {
571    type Item = Read<D>;
572
573    fn next(&mut self) -> Option<Self::Item> {
574        if (self.i < self.reads.n_reads()) && (self.i < self.end) {
575            let value = self.reads.get_read(self.i);
576            self.i += 1;
577            value
578        } else {
579            None
580        }
581    }
582}
583
584impl<D: Copy> ExactSizeIterator for ReadsIter<'_, D> {
585    fn len(&self) -> usize {
586        self.length
587    }
588}
589
590impl<D: Clone + Copy + Debug> Display for Reads<D> {
591    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
592        let vec: Vec<_> = self.iter().collect();
593        write!(f, "{:?}", vec)
594    }
595}
596
597#[derive(Debug, Clone, Serialize, Deserialize, PartialEq, Eq, PartialOrd, Ord, Hash)]
598pub enum ReadsPaired<D> {
599    Empty,
600    Unpaired { reads: Reads<D> },
601    Paired { paired1: Reads<D>, paired2: Reads<D> },
602    Combined {paired1: Reads<D>, paired2: Reads<D>, unpaired: Reads<D>}
603}
604
605impl<D: Clone + Copy> ReadsPaired<D> {
606    /// return an iterable element (`Vec`) with all contained `Reads`
607    pub fn iterable(&self) -> Vec<&Reads<D>> {
608        match self {
609            Self::Empty => vec![],
610            Self::Unpaired { reads  } => vec![reads],
611            Self::Paired { paired1, paired2 } => vec![paired1, paired2],
612            Self::Combined { paired1, paired2, unpaired } => vec![paired1, paired2, unpaired],
613        }
614    }
615
616    /// the overall number of reads
617    pub fn n_reads(&self) -> usize {
618        match self {
619            Self::Empty => 0,
620            Self::Unpaired { reads  } => reads.n_reads(),
621            Self::Paired { paired1, paired2 } => paired1.n_reads() + paired2.n_reads(),
622            Self::Combined { paired1, paired2, unpaired } => paired1.n_reads() + paired2.n_reads() + unpaired.n_reads(),
623        }
624    }
625
626    /// the number of paired reads
627    pub fn n_read_pairs(&self) -> usize {
628        match self {
629            Self::Empty => 0,
630            Self::Unpaired { reads: _ } => 0,
631            Self::Paired { paired1, paired2 } => {
632                assert_eq!(paired1.n_reads(), paired2.n_reads());
633                paired1.n_reads()
634            },
635            Self::Combined { paired1, paired2, unpaired: _ } => {
636                assert_eq!(paired1.n_reads(), paired2.n_reads());
637                paired1.n_reads()
638            }
639        }
640    }
641
642    /// the number of unpaired reads
643    pub fn n_unpaired_reads(&self) -> usize {
644        match self {
645            Self::Empty => 0,
646            Self::Unpaired { reads } => reads.n_reads(),
647            Self::Paired { paired1: _, paired2: _ } => 0,
648            Self::Combined { paired1: _, paired2: _, unpaired } => {
649                unpaired.n_reads()
650            }
651        }
652    }
653
654    /// get the read with the index `i` from the `ReadsPaired` - with multiple 
655    /// underlying `Reads` its is counted linearly trough paired1, paired2, unpaired
656    pub fn get_read(&self, i: usize) -> Option<Read<D>> {
657        match self {
658            ReadsPaired::Empty => None,
659            ReadsPaired::Unpaired { reads } => reads.get_read(i),
660            ReadsPaired::Paired { paired1, paired2 } => {
661                if i < paired1.n_reads() {
662                    paired1.get_read(i)
663                } else if (i - paired1.n_reads()) < paired2.n_reads() {
664                    paired2.get_read(i -  paired1.n_reads())
665                } else {
666                    None
667                }
668            },
669            ReadsPaired::Combined { paired1, paired2, unpaired } => {
670                if i < paired1.n_reads() {
671                    paired1.get_read(i)
672                } else if (i - paired1.n_reads()) < paired2.n_reads() {
673                    paired2.get_read(i -  paired1.n_reads())
674                } else if (i - (paired1.n_reads() + paired2.n_reads())) < unpaired.n_reads() {
675                    unpaired.get_read(i - (paired1.n_reads() + paired2.n_reads()))
676                } else {
677                    None
678                }
679            },
680        }
681    }
682
683    pub fn get_paired_read(&self, i: usize) -> Option<PairedRead<D>> {
684        match self {
685            ReadsPaired::Empty => None,
686            ReadsPaired::Unpaired { reads: _ } => None,
687            ReadsPaired::Paired { paired1, paired2 } => {
688                if let (Some(read1), Some(read2)) = (paired1.get_read(i), paired2.get_read(i)) {
689                    Some(PairedRead::new(read1, read2))
690                } else {
691                    None
692                }
693            }
694            ReadsPaired::Combined { paired1, paired2, unpaired: _ } => {
695                if let (Some(read1), Some(read2)) = (paired1.get_read(i), paired2.get_read(i)) {
696                    Some(PairedRead::new(read1, read2))
697                } else {
698                    None
699                }
700            },
701        }
702    }
703
704    pub fn mem(&self) -> usize {
705        match self {
706            Self::Empty => 0,
707            Self::Unpaired { reads } => reads.mem(),
708            Self::Paired { paired1, paired2 } => paired1.mem() + paired2.mem(),
709            Self::Combined { paired1, paired2, unpaired } => paired1.mem() + paired2.mem() + unpaired.mem(),
710        }
711    }
712
713    /// transform a tuple of two paired [`Reads`] and one unpaired [`Reads`] into a `ReadsPaired`
714    /// depending on the contents of the [`Reads`]
715    pub fn from_reads((paired1, paired2, unpaired): (Reads<D>, Reads<D>, Reads<D>)) -> Self {
716        // first two elements should be paired reads and thus have same n
717        assert_eq!(paired1.n_reads(), paired2.n_reads(), "Error: R1 read and R2 read counts have to match");
718
719        if (paired1.n_reads() + paired2.n_reads() + unpaired.n_reads()) == 0 {
720            // no reads
721            ReadsPaired::Empty
722        } else if paired1.n_reads() == 0 && unpaired.n_reads() > 0 {
723            // only reads in third element -> unpaired
724            ReadsPaired::Unpaired { reads: unpaired }
725        } else if paired1.n_reads() > 0 && unpaired.n_reads() == 0 {
726            // reads in first and second element -> paired
727            ReadsPaired::Paired { paired1, paired2 }
728        } else if paired1.n_reads() > 0 && unpaired.n_reads() > 0 {
729            // reads in all elements: both paired and unpaired reads
730            ReadsPaired::Combined { paired1, paired2, unpaired }
731        } else {
732            panic!("error in transforming Reads into ReadsPaired")
733        }
734    }
735
736    pub fn iter(&self) -> Box<dyn Iterator<Item = Read<D>> + '_> {
737        match self {
738            ReadsPaired::Empty => panic!("Error: no reads to process"),
739            ReadsPaired::Unpaired { reads } => Box::new(reads.iter()),
740            ReadsPaired::Paired { paired1, paired2 } => Box::new(paired1.iter().chain(paired2.iter())),
741            ReadsPaired::Combined { paired1, paired2, unpaired } => Box::new(paired1.iter().chain(paired2.iter()).chain(unpaired.iter())),
742        }
743    }
744
745    pub fn iter_partial(&self, range: Range<usize>) -> Box<dyn Iterator<Item = Read<D>> + '_> {
746        match self {
747            Self::Empty => panic!("Error: no reads to process"),
748            Self::Unpaired { reads } => Box::new(reads.partial_iter(range)),
749            Self::Paired { paired1, paired2 } => {
750                let n_p1 = paired1.n_reads();
751                if range.start >= n_p1 {
752                    // range is fully in paired2
753                    Box::new(paired2.partial_iter((range.start - n_p1)..(range.end - n_p1)))
754                } else if range.end <= n_p1 {
755                    // range is fully in paired1
756                    Box::new(paired1.partial_iter(range))
757                } else {
758                    // range is both in paired1 and paired2
759                    Box::new(paired1.partial_iter(range.start..n_p1).chain(paired2.partial_iter(0..(range.end - n_p1))))
760                }
761            },
762            Self::Combined { paired1, paired2, unpaired } => {
763                let n_p1 = paired1.n_reads();
764                let n_p2 = paired2.n_reads();
765                let n_p12 = n_p1 + paired2.n_reads();
766                if range.end <= n_p1 {
767                    // range is only in paired1
768                    Box::new(paired1.partial_iter(range))
769                } else if range.end >= n_p1 && range.end <= n_p12 && range.start >= n_p1 && range.start <= n_p12 {
770                    // range is only in paired2
771                    Box::new(paired2.partial_iter((range.start - n_p1)..(range.end - n_p1)))
772                } else if range.start >= n_p12 {
773                    // range is only in unpaired
774                    Box::new(unpaired.partial_iter((range.start - n_p12)..(range.end - n_p12)))
775                } else if range.start <= n_p1 && range.end >= n_p1 && range.end <= n_p12 {
776                    // range is in paired1 and paired2
777                    Box::new(paired1.partial_iter(range.start..n_p1).chain(paired2.partial_iter(0..(range.end - n_p1))))
778                } else if range.start >= n_p1 && range.start <= n_p12 && range.end >= n_p12 {
779                    // range is in paired2 and unpaired
780                    Box::new(paired2.partial_iter((range.start - n_p1)..n_p2).chain(unpaired.partial_iter(0..(range.end - n_p12))))
781                } else {
782                    // range is in paired1, paired2, and in unpaired
783                    Box::new(paired1.partial_iter(range.start..n_p1).chain(paired2.partial_iter(0..n_p2)).chain(unpaired.partial_iter(0..(range.end - n_p12))))
784                }
785            }
786        }
787    }
788
789    /// if the `ReadsPaired` is of `Combined` type, remove the unpaired reads,
790    /// returns the number of reads that were removed
791    pub fn decombine(&mut self) -> usize {
792        if let Self::Combined { paired1, paired2, unpaired } = self {
793            let rm_reads = unpaired.n_reads();
794            *self = ReadsPaired::Paired { paired1: take(paired1), paired2: take(paired2) };
795            rm_reads
796        } else {
797            0
798        }
799    }
800}
801
802impl<DI: ReadData> ReadsPaired<DI> {
803    /// get the number of k-mers for each unique data value
804    pub fn tag_kmers(&self, k: usize) -> HashMap<Tag, usize> {
805        match self {
806            Self::Empty => HashMap::new(),
807            Self::Unpaired { reads } => reads.tag_kmers(k),
808            Self::Paired { paired1, paired2 } => {
809                let mut hm_p1 = paired1.tag_kmers(k);
810                let hm_p2 = paired2.tag_kmers(k);
811
812                // combine the values for underlying Reads
813                hm_p2.into_iter().for_each(|(data, kmers)| {
814                   if let Some(count) = hm_p1.get_mut(&data) {
815                    *count += kmers;
816                   } else {
817                    hm_p1.insert(data, kmers);
818                   }
819                });
820
821                hm_p1
822            },
823            Self::Combined { paired1, paired2, unpaired } => {
824                let mut hm_p1: HashMap<u8, usize> = paired1.tag_kmers(k);
825                let hm_p2 = paired2.tag_kmers(k);
826                let hm_up = unpaired.tag_kmers(k);
827
828                // combine the values for underlying Reads
829                hm_p2.into_iter().for_each(|(data, kmers)| {
830                   if let Some(count) = hm_p1.get_mut(&data) {
831                    *count += kmers;
832                   } else {
833                    hm_p1.insert(data, kmers);
834                   }
835                });
836
837                hm_up.into_iter().for_each(|(data, kmers)| {
838                    if let Some(count) = hm_p1.get_mut(&data) {
839                     *count += kmers;
840                    } else {
841                     hm_p1.insert(data, kmers);
842                    }
843                 });
844
845                hm_p1
846            }
847        }
848    }
849
850    /// return the number of k-mers occuring with each u8-encoded tag, 
851    /// with the tag as the index
852    /// if there are no tags saved in the Readspauired, it returns a vector of the
853    /// length `n_sampeles`, filles with zeroes
854    pub fn tag_kmers_vec(&self, k: usize, n_samples: usize) -> Vec<u64> {
855        let hashed_kmer_counts = self.tag_kmers(k);
856
857        let mut kmer_counts = vec![0; n_samples];
858
859        for (tag, kmer_count) in hashed_kmer_counts {
860            kmer_counts[tag as usize] += kmer_count as u64;
861        }
862
863        kmer_counts
864    }
865}
866
867impl<D: Clone + Copy> Display for ReadsPaired<D> {
868    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
869        match self {
870            Self::Empty => write!(f, "empty ReadsPaired"),
871            Self::Unpaired { reads } => write!(f, "unpaired ReadsPaired: \n{}", reads.info()),
872            Self::Paired { paired1, paired2 } => write!(f, "paired ReadsPaired: \n{}\n{}", paired1.info(), paired2.info()),
873            Self::Combined { paired1, paired2, unpaired } => write!(f, "combined ReadsPaired: \n{}\n{}\n{}", paired1.info(), paired2.info(), unpaired.info()),
874        }
875    }
876}
877
878/// Trait for ReadData, [`ID`]s can only be generated from [Marbel](https://github.com/jlab/marbel) reads
879pub trait ReadData: PartialEq + Hash + serde::Serialize + DeserializeOwned + Debug + Clone + Copy + Eq + Send + Sync + Ord {
880    /// generate a read data from an ID and a tag
881    fn new(id: ID, tag: Tag) -> Self;
882    /// geneate a read data, [`ID`]s and [`IDTag`]s can only be generated from [Marbel](https://github.com/jlab/marbel) reads
883    fn read_data(gene_ids: &mut BiMap<String, ID>, read_name: &[u8], tag: Tag) -> Self;
884    /// if available, get a tag
885    fn get_tag(&self) -> Option<Tag>;
886    /// retrun a ReadDatas enum to check which kind of ReadData is present
887    fn read_datas() -> ReadDatas;
888}
889
890impl ReadData for Tag {
891    fn new(_id: ID, tag: Tag) -> Self {
892        tag
893    }
894
895    fn read_data(_: &mut BiMap<String, ID>, _: &[u8], tag: Tag) -> Self {
896        tag
897    }
898
899    fn get_tag(&self) -> Option<Tag> {
900        Some(*self)
901    }
902
903    fn read_datas() -> ReadDatas {
904        ReadDatas::Tag
905    }
906}
907
908impl ReadData for ID {
909    fn new(id: ID, _tag: Tag) -> Self {
910        id    
911    }
912
913    fn read_data(gene_ids: &mut BiMap<String, ID>, read_name: &[u8], _: Tag) -> Self {
914
915        // read name is e.g. "B7R87_RS28825_2_0/1" -> gene: "B7R87_RS28825"
916        // split at '_' and use first two elements and reconnect with '_'
917        let read_name_sting = str::from_utf8(read_name).expect("error reading read name").to_string();
918        let mut split_iter = read_name_sting.split('_');
919        let mut gene = String::new();
920
921        let Some(gene1) = split_iter.next() else {
922            panic!("no gene names found in reads - only use id summarizers with marbel data - read name: {}", read_name_sting)
923        };
924        gene.push_str(gene1);
925
926        gene.push('_');
927
928        let Some(gene2) = split_iter.next() else {
929            panic!("no gene names found in reads - only use id summarizers with marbel data - read name: {}", read_name_sting)
930        };
931        gene.push_str(gene2);
932
933        // if gene not in gene_ids, then add, else get gene id
934        let new_id = gene_ids.len() as ID;
935        if new_id == ID::MAX { panic!("number of genes has surpassed 65.5k limit of u16") }
936        match gene_ids.get_by_left(&gene) {
937            Some(id) => *id,
938            None => {
939                gene_ids.insert(gene, new_id);
940                new_id
941            },
942        }
943    }
944
945    fn get_tag(&self) -> Option<Tag> {
946        None
947    }
948
949    fn read_datas() -> ReadDatas {
950        ReadDatas::ID
951    }
952}
953
954impl ReadData for IDTag {
955    fn new(id: ID, tag: Tag) -> Self {
956        Self::new(id, tag)
957    }
958
959    fn read_data(gene_ids: &mut BiMap<String, ID>, read_name: &[u8], tag: Tag) -> Self {
960        let id = ID::read_data(gene_ids, read_name, tag);
961        IDTag::new(id, tag)
962    }
963
964    fn get_tag(&self) -> Option<Tag> {
965        Some(self.tag())
966    }
967
968    fn read_datas() -> ReadDatas {
969        ReadDatas::IDTag
970    }
971}
972
973#[derive(PartialEq, Eq, Debug, Clone, Copy, Serialize, Deserialize)]
974pub enum ReadDatas {
975    ID,
976    Tag,
977    IDTag
978}
979
980#[cfg(test)]
981mod tests {
982    use std::{collections::HashMap, time};
983
984    use bimap::BiMap;
985    use itertools::enumerate;
986    use rand::random;
987
988    use crate::{Exts, QualityVec, dna_string::DnaString, reads::{Read, Strandedness}, summarizer::{ID, IDTag, Tag}, test::random_dna};
989    use crate::reads::ReadData;
990    use super::{Reads, ReadsPaired};
991
992    #[test]
993    fn test_add() {
994
995        let fastq = vec![
996            (DnaString::from_acgt_bytes(str::as_bytes("ACGATCGT")), Exts::empty(), 6u8, str::as_bytes("CCC-CC;#")),
997            (DnaString::from_acgt_bytes(str::as_bytes("GGGGGG")), Exts::empty(), 5u8, str::as_bytes("CCCCCC")),
998            (DnaString::from_acgt_bytes(str::as_bytes("TTGGTT")), Exts::empty(), 7u8, str::as_bytes("CCC-CC")),
999            (DnaString::from_acgt_bytes(str::as_bytes("ACCAC")), Exts::empty(), 8u8, str::as_bytes("-CC;#")),
1000            (DnaString::from_acgt_bytes(str::as_bytes("TCCCT")), Exts::empty(), 9u8, str::as_bytes("CCCC;")),
1001            (DnaString::from_acgt_bytes(str::as_bytes("ACCAC")), Exts::empty(), 8u8, str::as_bytes("C-C;C")),
1002            (DnaString::from_acgt_bytes(str::as_bytes("TCCCT")), Exts::empty(), 9u8, str::as_bytes("CCC-C")),
1003        ];
1004
1005
1006        let mut reads = Reads::new_with_quality(Strandedness::Unstranded);
1007        for (read, _, data, quality) in fastq.clone() {
1008            reads.add_read(read, None, data, Some(quality));
1009        }
1010
1011        println!("reads: {:#?}", reads);
1012
1013
1014        /* for no in reads.storage.iter() {
1015            println!("{:#b}", no)
1016        } */
1017
1018        assert_eq!(reads.storage, vec![1791212948343256433, 5140577499666710528]);
1019
1020        for (i, _) in fastq.iter().enumerate() {
1021            //println!("read {}: {:?}", i, reads.get_read(i))
1022            assert_eq!(
1023                Read::new(
1024                    fastq[i].0.clone(), 
1025                    fastq[i].1, 
1026                    fastq[i].2, 
1027                    Strandedness::Unstranded, 
1028                    Some(QualityVec::from_ascii_bytes(fastq[i].3, reads.quality_bins))
1029                ), 
1030                reads.get_read(i).unwrap())
1031        }
1032
1033        for read in reads.iter() {
1034            let seq = read.seq;
1035            println!("{:?}, {}", seq, seq.len())
1036        }
1037        println!();
1038
1039        for read in reads.partial_iter(5..7) {
1040            println!("{:?}", read)
1041        }
1042
1043        println!("memory usage: {}", reads.mem())
1044    }
1045
1046    #[test]
1047    fn test_get_read() {
1048        let mut reads = Reads::new(Strandedness::Unstranded);
1049        //reads.add_read(DnaString::from_acgt_bytes("AGCTAGCTAGC".as_bytes()), Exts::empty(), 67u8);
1050        reads.add_from_bytes("ACGATCGNATGCTAGCTGATCGGCGACGATCGATGCTAGCTGATCGTAGCTGACTGATCGATCG".as_bytes(), None, 67u8);
1051        let read = reads.get_read(0);
1052        println!("{:?}", read);
1053        println!("{:#066b}", reads.storage[0]);
1054        println!("{:?}", reads)
1055    }
1056
1057    #[test]
1058    fn test_add_from_bytes() {
1059        let dna = [
1060            "AAGCGGAGATTATTCACGAGCATCGCGTAC".as_bytes(),
1061            "GATCGATGCATGCTAGA".as_bytes(),
1062            "ACGTAAAAAAAAAATTATATAACGTACGTAAAAAAAAAATTATATAACGTAACGTAAAAAAAAAAATTATAATAACGT".as_bytes(),
1063            "AGCTAGCTAGCTGACTGAGCGACTGA".as_bytes(),
1064            "AGCTAGCTAGCTGACTGAGCGACTGACGGATC".as_bytes(),
1065            "TTTTTTTTTTTTTTTTTTTTTTTT".as_bytes(),
1066            "ACGATCGAATGCTAGCTGATCGGCGACGATCGATGCTAGCTGATCGTAGCTGACTGATCGATCG".as_bytes(),
1067            "ACGATCGATGCTAGCTGATCGGCGACGATCGATGCTAGCTGATCGTAGCTGACTGATCGATCGAAGGGCAGTTAGGCCGTAAGCGCGAT".as_bytes(),
1068        ];
1069
1070        let mut reads: Reads<u8> = Reads::new(Strandedness::Unstranded);
1071        for seq in dna {
1072            reads.add_from_bytes(seq, None, random());
1073
1074        }
1075
1076        for (i, read) in enumerate(reads.iter()) {
1077            let sequence = DnaString::from_acgt_bytes(dna[i]);
1078            assert_eq!(read.seq, sequence);
1079        }
1080    }
1081
1082    #[test]
1083    fn test_add_from_bytes_checked() {
1084        let dna = [
1085            "AAGCGGAGATTATTCACGAGCATCGCGTAC".as_bytes(),
1086            "GATCGATGCATGCTAGA".as_bytes(),
1087            "ACGTAAAAAAAAAATTATATAACGTACGTAAAAAAAAAANTTATATAACGTAACGTAAAAAAAAAAATTATAATAACGT".as_bytes(),
1088            "AGCTAGCTAGCTGACNGAGCGACTGA".as_bytes(),
1089            "AGCTAGCTAGCTGACTGAGCGACTGACGGATC".as_bytes(),
1090            "TTTTTTTTTTTTTTTTTTTTTTTT".as_bytes(),
1091            "ACGATCGAATGCTAGCTGATCGGCGACGATCGATGCTAGCTGATCGTAGCTGACNNNTGATCGATCG".as_bytes(),
1092            "ACGATCGATGCTAGCTGATCGGCGACGATCGATGCTAGCTGATCGTAGCTGACTGATCGATCGAAGGGCAGTTAGGCCGTAAGCGCGAT".as_bytes(),
1093            "A".as_bytes(),
1094            "AAAAN".as_bytes(),
1095            "NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN".as_bytes(),
1096        ];
1097
1098        let mut reads: Reads<u8> = Reads::new(Strandedness::Unstranded);
1099        let mut corrects = Vec::new();
1100        for seq in dna {
1101            corrects.push(reads.add_from_bytes_checked(seq, None, random()));
1102        }
1103
1104        let mut read_counter = 0;
1105
1106        for (i, correct) in enumerate(corrects) {
1107            let sequence = DnaString::from_acgt_bytes_checked(dna[i]);
1108            match correct {
1109                true => {
1110                    let read = reads.get_read(read_counter).unwrap();
1111                    assert_eq!(sequence.unwrap(), read.seq);
1112                    read_counter += 1;
1113
1114                },
1115                false => assert!(sequence.is_err()),
1116            }
1117        }
1118    }
1119
1120    #[test]
1121    fn test_speed_from_bytes() {
1122        let dnas = [
1123            "AAGCGGAGATTATTCACGAGCATCGCGTAC".as_bytes(),
1124            "GATCGATGCATGCTAGA".as_bytes(),
1125            "ACGTAAAAAAAAAATTATATAACGTACGTAAAAAAAAAANTTATATAACGTAACGTAAAAAAAAAAATTATAATAACGT".as_bytes(),
1126            "AGCTAGCTAGCTGACNGAGCGACTGA".as_bytes(),
1127            "AGCTAGCTAGCTGACTGAGCGACTGACGGATC".as_bytes(),
1128            "TTTTTTTTTTTTTTTTTTTTTTTT".as_bytes(),
1129            "ACGATCGAATGCTAGCTGATCGGCGACGATCGATGCTAGCTGATCGTAGCTGACNNNTGATCGATCG".as_bytes(),
1130            "ACGATCGATGCTAGCTGATCGGCGACGATCGATGCTAGCTGATCGTAGCTGACTGATCGATCGAAGGGCAGTTAGGCCGTAAGCGCGAT".as_bytes(),
1131            "A".as_bytes(),
1132            "AAAAN".as_bytes(),
1133            "NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN".as_bytes(),
1134        ];
1135
1136        const REPS: usize = 500000;
1137        /*
1138        test with 5 ml:
1139            through DnaString: 88.09323 s 
1140            direct to Read: 44.41913 s
1141         */
1142
1143
1144        let ds_start= time::Instant::now();
1145        let mut reads: Reads<u8> = Reads::new(Strandedness::Unstranded);
1146        for _i in 0..REPS {
1147            for dna in dnas {
1148                reads.add_read(DnaString::from_acgt_bytes(dna), None, random(), None);
1149            }
1150        }
1151        let ds_finish = ds_start.elapsed();
1152
1153        let r_start= time::Instant::now();
1154        let mut reads: Reads<u8> = Reads::new(Strandedness::Unstranded);
1155        for _i in 0..REPS {
1156            for dna in dnas {
1157                reads.add_from_bytes(dna, None, random());
1158            }
1159        }
1160        let r_finish = r_start.elapsed();
1161
1162        println!("through DnaString: {} s \n direct to Read: {} s", ds_finish.as_secs_f32(), r_finish.as_secs_f32())
1163
1164
1165    }
1166
1167
1168    #[test]
1169    fn test_reads_stranded() {
1170        let mut reads: Reads<u8> = Reads::new(Strandedness::Forward);
1171        assert_eq!(reads.stranded(), Strandedness::Forward);
1172        reads.set_stranded(Strandedness::Reverse);
1173        assert_eq!(reads.stranded(), Strandedness::Reverse);
1174
1175    }
1176
1177    #[test]
1178    fn test_reads_data_kmers() {
1179        let mut reads = Reads::new(Strandedness::Unstranded);
1180        let seqs = [
1181            ("ACGATCGTACGTACGTAGCTAGCTGCTAGCTAGCTGACTGACTGA", 0),
1182            ("CGATGCTATCAGCGAGCGATCGTACGTAGCTACG", 1),
1183            ("CGATCGACGAGCAGCGTATGCTACGAGCTGACGATCTACGA", 2),
1184            ("CACACACGGCATCGATCGAGCAGCATCGACTACGTA", 3),
1185        ];
1186
1187        seqs.iter().for_each(|(read, tag)| reads.add_from_bytes(read.as_bytes(), None, *tag as u8));
1188        let data_kmers = reads.tag_kmers(16);
1189       
1190        let comp_hm: HashMap<u8, usize> = [(0, 30), (1, 19), (2, 26), (3, 21)].into_iter().collect();
1191
1192        assert_eq!(comp_hm, data_kmers);
1193    }
1194
1195    #[test]
1196    fn test_reads_add_exts() {
1197        let mut raw_reads = Vec::new();
1198        for _i in 0..10 {
1199            raw_reads.push((DnaString::from_bytes(&random_dna(100)), Exts::new(rand::random::<u8>()), rand::random::<u8>()));
1200        }
1201        let reads = Reads::from_vmer_vec(raw_reads.clone(), Strandedness::Unstranded);
1202        let new_raw_reads = reads.iter().map(|read| (read.seq, read.exts, read.data)).collect::<Vec<_>>();
1203
1204        assert_eq!(raw_reads, new_raw_reads)
1205    }
1206
1207    #[test]
1208    fn test_reads_info() {
1209        let mut reads = Reads::new(Strandedness::Unstranded);
1210        let seqs = [
1211            ("ACGATCGTACGTACGTAGCTAGCTGCTAGCTAGCTGACTGACTGA", 0),
1212            ("CGATGCTATCAGCGAGCGATCGTACGTAGCTACG", 1),
1213            ("CGATCGACGAGCAGCGTATGCTACGAGCTGACGATCTACGA", 2),
1214            ("CACACACGGCATCGATCGAGCAGCATCGACTACGTA", 3),
1215        ];
1216
1217        seqs.iter().for_each(|(read, tag)| reads.add_from_bytes(read.as_bytes(), None, *tag as u8));
1218
1219        assert_eq!(reads.info(), "Reads { n reads: 4, stranded: Unstranded }".to_string());
1220    }
1221
1222    #[test]
1223    fn test_reads_paired() {
1224        let mut p1 = Reads::new(Strandedness::Unstranded);
1225        let mut p2 = Reads::new(Strandedness::Unstranded);
1226        let mut up = Reads::new(Strandedness::Unstranded);
1227
1228        let reads_p1 = [
1229            "ACGATCGTACGTACGTAGCTAGCTGCTAGCTAGCTGACTGACTGA",
1230            "CGATGCTATCAGCGAGCGATCGTACGTAGCTACG",
1231            "CGATCGACGAGCAGCGTATGCTACGAGCTGACGATCTACGA",
1232            "CACACACGGCATCGATCGAGCAGCATCGACTACGTA",
1233        ];
1234
1235        let reads_p2 = [
1236            "AGCTAGCTAGCTACTGATCGTAGCTAGCTGATCGA",
1237            "AGCGATCGTACGTAGCTAGCTA",
1238            "CGATCGATCGACTAGCGTAGCTGACTGAC",
1239            "CAGATGCTCTGCTGACTGACTGATCGTACTGACTAGCATCTAGC",
1240        ];
1241
1242        let reads_up = [
1243            "CGTACTAGCTGACGTAC",
1244            "CGATGCTAGCTAGCTAGCGATCG",
1245        ];
1246
1247        let tags = (0..4).collect::<Vec<u8>>();
1248
1249        reads_p1.iter().enumerate().for_each(|(i, read)| p1.add_from_bytes(read.as_bytes(), None, tags[i]));
1250        reads_p2.iter().enumerate().for_each(|(i, read)| p2.add_from_bytes(read.as_bytes(), None, tags[i]));
1251        reads_up.iter().enumerate().for_each(|(i, read)| up.add_from_bytes(read.as_bytes(), Some(Exts::new(i as u8)), tags[i]));
1252
1253        let empty: ReadsPaired<u8> = ReadsPaired::from_reads((Reads::new(Strandedness::Unstranded), Reads::new(Strandedness::Unstranded), Reads::new(Strandedness::Unstranded)));
1254        assert_eq!(empty, ReadsPaired::Empty);
1255        assert_eq!(empty.mem(), 0);
1256        assert_eq!(empty.n_reads(), 0);
1257        assert_eq!(empty.iterable(), Vec::<&Reads<u8>>::new());
1258
1259        let unpaired = ReadsPaired::from_reads((Reads::new(Strandedness::Unstranded), Reads::new(Strandedness::Unstranded), up.clone()));
1260        assert_eq!(ReadsPaired::Unpaired { reads: up.clone() }, unpaired);
1261        println!("exts up {:?}", unpaired);
1262        assert_eq!(unpaired.mem(), 172);
1263        assert_eq!(unpaired.n_reads(), 2);
1264        assert_eq!(unpaired.iterable(), vec![&up]);
1265      
1266        let paired = ReadsPaired::from_reads((p1.clone(), p2.clone(), Reads::new(Strandedness::Unstranded)));
1267        assert_eq!(ReadsPaired::Paired { paired1: p1.clone(), paired2: p2.clone() }, paired);        
1268        assert_eq!(paired.mem(), 424);
1269        assert_eq!(paired.n_reads(), 8);
1270        assert_eq!(paired.iterable(), vec![&p1, &p2]);
1271
1272        let combined = ReadsPaired::from_reads((p1.clone(), p2.clone(), up.clone()));
1273        assert_eq!(ReadsPaired::Combined { paired1: p1.clone(), paired2: p2.clone(), unpaired: up.clone() }, combined);
1274        assert_eq!(combined.mem(), 596);
1275        assert_eq!(combined.n_reads(), 10);
1276        assert_eq!(combined.iterable(), vec![&p1, &p2, &up]);
1277
1278
1279        // test iter
1280
1281        assert_eq!(unpaired.iter().collect::<Vec<_>>(), up.iter().collect::<Vec<_>>());
1282        assert_eq!(paired.iter().collect::<Vec<_>>(), p1.iter().chain(p2.iter()).collect::<Vec<_>>());
1283        assert_eq!(combined.iter().collect::<Vec<_>>(), p1.iter().chain(p2.iter()).chain(up.iter()).collect::<Vec<_>>());
1284
1285        // test partial iter
1286
1287        assert_eq!(unpaired.iter_partial(0..1).collect::<Vec<_>>(), up.partial_iter(0..1).collect::<Vec<_>>());
1288
1289        assert_eq!(paired.iter_partial(0..1).collect::<Vec<_>>(), p1.partial_iter(0..1).collect::<Vec<_>>());
1290        assert_eq!(paired.iter_partial(5..7).collect::<Vec<_>>(), p2.partial_iter(1..3).collect::<Vec<_>>());
1291        assert_eq!(paired.iter_partial(1..8).collect::<Vec<_>>(), p1.partial_iter(1..4).chain(p2.partial_iter(0..4)).collect::<Vec<_>>());
1292
1293        assert_eq!(combined.iter_partial(0..1).collect::<Vec<_>>(), p1.partial_iter(0..1).collect::<Vec<_>>());
1294        assert_eq!(combined.iter_partial(5..7).collect::<Vec<_>>(), p2.partial_iter(1..3).collect::<Vec<_>>());
1295        assert_eq!(combined.iter_partial(8..10).collect::<Vec<_>>(), up.partial_iter(0..2).collect::<Vec<_>>());
1296        assert_eq!(combined.iter_partial(1..8).collect::<Vec<_>>(), p1.partial_iter(1..4).chain(p2.partial_iter(0..4)).collect::<Vec<_>>());
1297        assert_eq!(combined.iter_partial(6..9).collect::<Vec<_>>(), p2.partial_iter(2..4).chain(up.partial_iter(0..1)).collect::<Vec<_>>());
1298        assert_eq!(combined.iter_partial(1..9).collect::<Vec<_>>(), p1.partial_iter(1..4).chain(p2.partial_iter(0..4)).chain(up.partial_iter(0..1)).collect::<Vec<_>>());
1299
1300        // test tag kmers (and data_kmers)
1301        assert_eq!(unpaired.tag_kmers_vec(16, 2), vec![2, 8]);
1302        assert_eq!(paired.tag_kmers_vec(16, 4), vec![50, 26, 40, 50]);
1303        assert_eq!(combined.tag_kmers_vec(16, 4), vec![52, 34, 40, 50]);
1304        assert_eq!(empty.tag_kmers_vec(16, 0), Vec::<u64>::new());
1305
1306        // test decombine
1307        let mut paired_dc = paired.clone();
1308        let rm_reads_p = paired_dc.decombine();
1309        let mut combined_dc = combined.clone();
1310        let rm_reads_up = combined_dc.decombine();
1311        assert_eq!(paired_dc, paired);
1312        assert_eq!(rm_reads_p, 0);
1313        assert_eq!(combined_dc, paired);
1314        assert_eq!(rm_reads_up, 2);
1315
1316        // display
1317        assert_eq!(format!("{}", empty), "empty ReadsPaired".to_string());
1318        assert_eq!(format!("{}", unpaired), "unpaired ReadsPaired: 
1319Reads { n reads: 2, stranded: Unstranded }".to_string());
1320        assert_eq!(format!("{}", paired), "paired ReadsPaired: 
1321Reads { n reads: 4, stranded: Unstranded }
1322Reads { n reads: 4, stranded: Unstranded }".to_string());
1323        assert_eq!(format!("{}", combined), "combined ReadsPaired: 
1324Reads { n reads: 4, stranded: Unstranded }
1325Reads { n reads: 4, stranded: Unstranded }
1326Reads { n reads: 2, stranded: Unstranded }".to_string());
1327
1328    }
1329
1330    #[test]
1331    #[should_panic]
1332    fn test_reads_paired_panic() {
1333        let mut p1 = Reads::new(Strandedness::Unstranded);
1334
1335        let reads_p1 = [
1336            "ACGATCGTACGTACGTAGCTAGCTGCTAGCTAGCTGACTGACTGA",
1337            "CGATGCTATCAGCGAGCGATCGTACGTAGCTACG",
1338            "CGATCGACGAGCAGCGTATGCTACGAGCTGACGATCTACGA",
1339            "CACACACGGCATCGATCGAGCAGCATCGACTACGTA",
1340        ];
1341
1342        reads_p1.iter().for_each(|read| p1.add_from_bytes(read.as_bytes(), None, 0u8));
1343
1344        let _ = ReadsPaired::from_reads((p1, Reads::new(Strandedness::Unstranded), Reads::new(Strandedness::Unstranded)));
1345    }
1346
1347    #[test]
1348    fn test_read_data() {
1349        let read_name_1 = "B7R87_RS28825_2_0/1".as_bytes(); // gene B7R87_RS28825
1350        let read_name_2 = "B7R87_RS21825_2_0/1".as_bytes(); // gene B7R87_RS21825
1351
1352        let tag = 0 as Tag;
1353        
1354        let mut ids = BiMap::new();
1355
1356        let id = ID::read_data(&mut ids, read_name_1, tag);
1357        assert_eq!(id, 0);
1358
1359        let id_tag = IDTag::read_data(&mut ids, read_name_1, tag);
1360        assert_eq!(id_tag, IDTag::new(id, tag));
1361        assert_eq!(id_tag.get_tag(), Some(tag));
1362
1363        let id = ID::read_data(&mut ids, read_name_2, tag);
1364        assert_eq!(id, 1);
1365        assert_eq!(id.get_tag(), None);
1366
1367        assert_eq!(Tag::read_data(&mut ids, read_name_1, tag), tag);
1368        assert_eq!(tag.get_tag(), Some(tag));
1369    }
1370
1371    #[test]
1372    #[should_panic]
1373    fn test_read_data_panic() {
1374        let mut ids = BiMap::new();
1375
1376        let read_name_1 = "B7R87_RS28825_2_0/1".as_bytes(); // gene B7R87_RS28825
1377        let read_name_2 = "B7R87_RS21825_2_0/1".as_bytes(); // gene B7R87_RS21825
1378
1379        let _ = ID::read_data(&mut ids, read_name_1, 0);
1380        let _ = ID::read_data(&mut ids, read_name_2, 0);
1381
1382        // trying to add invalid read name -> panic
1383        let _ = ID::read_data(&mut ids, "AAAAAA".as_bytes(), 0);
1384    }
1385}
debruijn/reads.rs

debruijn/
reads.rs
