debruijn/
summarizer.rs

1use bimap::BiMap;
2use clap::ValueEnum;
3use serde::{de::DeserializeOwned, Deserialize, Serialize};
4use statrs::distribution::{ContinuousCDF, Normal, StudentsT};
5use summarydata_derive::SummaryData;
6use crate::{BaseQuality, EdgeMap, EdgeMult, Exts, Kmer, KmerDataItem, Tags};
7use std::{cmp::min_by, collections::HashMap, error::Error, fmt::{Debug, Display}, mem};
8
9/// inner type for [`Tags`] and group markers
10#[cfg(not(feature = "sample128"))]
11pub type Marker = u64;
12
13/// inner type for [`Tags`] and group markers
14#[cfg(feature = "sample128")]
15pub type Marker = u128;
16
17/// type for IDs (e.g. gene IDs)
18#[cfg(not(feature = "id4b"))]
19pub type ID = u16;
20
21/// type for IDs (e.g. gene IDs)
22#[cfg(feature = "id4b")]
23pub type ID = u32;
24
25/// type for tags
26pub type Tag = u8;
27
28#[derive(Debug, PartialEq, Eq, PartialOrd, Ord, Clone, Copy, Serialize, Deserialize, Hash)]
29/// type for IDs and tags together
30pub struct IDTag {
31    id: ID,
32    tag: Tag
33}
34
35impl IDTag {
36    pub fn new(id: ID, tag: Tag) -> IDTag {
37        IDTag { id, tag }
38    }
39
40    pub fn tag(&self) -> Tag {
41        self.tag
42    }
43}
44
45#[derive(Debug, Clone, Serialize, Deserialize, PartialEq, Eq)]
46/// translate tags and IDs (e.g. into sample labels and gene names)
47pub struct Translator {
48    ids: Option<BiMap<String, ID>>,
49    tags: Option<BiMap<String,Tag>>,
50}
51
52impl Translator {
53    /// make a new [`Translator`] for tags and IDs
54    pub fn new(ids: BiMap<String, ID>, tags: BiMap<String,Tag>) -> Translator {
55        Translator { ids: Some(ids), tags: Some(tags) }
56    }
57
58    /// make an empty [`Translator`]
59    pub fn empty() -> Translator {
60        Translator { ids: None, tags: None }
61    }
62
63    /// make a new [`Translator`] for tags
64    pub fn new_tag_translator(hashed_tags: BiMap<String, Tag>) -> Translator {
65        Translator { ids: None, tags: Some(hashed_tags) }
66    }
67
68    /// make a new [`Translator`] for IDs
69    pub fn new_id_translator(hashed_ids: BiMap<String, ID>) -> Translator {
70        Translator { ids: Some(hashed_ids), tags: None }
71    }
72
73    /// get the tag translator, returns None if the `Translator` does not contain a tag translator
74    pub fn tag_translator(&self) -> &Option<BiMap<String, Tag>> {
75        &self.tags
76    }
77
78    /// get the tag translator, returns None if the `Translator` does not contain a tag translator
79    pub fn id_translator(&self) -> &Option<BiMap<String, ID>> {
80        &self.ids
81    }
82
83    /// get a mutable reference to the tag translator, returns None if the `Translator` does not contain a tag translator
84    pub fn mut_id_translator(&mut self) -> &mut Option<BiMap<String, ID>> {
85        &mut self.ids
86    }
87
88    /// dissolve the `Translator` into its underlying [`BiMap`]s
89    pub fn dissolve(self) -> (Option<BiMap<String, ID>>, Option<BiMap<String,Tag>>) {
90        (self.ids, self.tags)
91    }
92}
93
94fn id_format(ids: &[ID], translator: &Translator, id_group_translator: Option<&HashMap<ID, ID>>) -> String {
95    if let Some(id_gr_tr) = id_group_translator {
96        // translate the ids (genes) into id groups (orthogroups)
97        let mut t_ids = ids
98            .iter()
99            .map(|id| id_gr_tr.get(id).unwrap_or_else(|| panic!("ID does not exist - ids {:?}", ids)))
100            .collect::<Vec<_>>();
101        t_ids.sort();
102        t_ids.dedup();
103        format!("{:?}", t_ids)
104    } else if let Some(id_translator) = translator.id_translator() {
105        // translate the ids (genes) into their names
106        let t_ids = ids
107            .iter()
108            .map(|id| id_translator.get_by_right(id).unwrap_or_else(|| panic!("ID does not exist - ids {:?}", ids)))
109            .collect::<Vec<_>>();
110        format!("{:?}", t_ids)
111    } else {
112        // do not translate
113        format!("{:?}", ids)
114    }
115}
116
117#[derive(Debug, PartialEq)]
118struct NotEnoughSamplesError {}
119
120impl Error for NotEnoughSamplesError {}
121
122impl Display for NotEnoughSamplesError {
123    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
124        write!(f, "not enough samples were supplied to perform a statistical test")
125    }
126}
127
128/// Configuration for summary processes. It it used to filter the k-mers during
129/// and after graph construction and to make statistical analyses based on k-mer
130/// occurrence in the sample groups. 
131/// 
132/// For any statistical analyses regarding the sample groups, the `SummaryData` 
133/// requires a [`SampleInfo`].
134/// 
135/// The available options are:
136/// - filter the k-mers by 
137///     - their number of occurrences
138///     - their quality based on the phred scores from the reads
139///     - their p-values regarding sample groups
140///     - if they occurr in at least a specific fraction of one or both of
141///       the sample groups 
142/// - filter the k-mer occurrences by the quality - should this disconnect the 
143///   k-mer, alll occurrences will be used
144/// - set the number of occurrences which is stored with the k-mers to be rounded
145///   to a number of significant digits
146/// - choose a statistical test on which the p-value calculation is based, by default 
147///   this is set to Welch's t-test; Student's t-test and the Mann-Whitney U test 
148///   are also available
149/// 
150/// By default, all filter options are turned off.
151/// 
152/// ```
153/// use debruijn::summarizer::{GroupFrac, SummaryConfig, SampleInfo, StatTest};
154/// use debruijn::BaseQuality;
155/// 
156/// let sample_info = SampleInfo::new(0b1100, 0b0011, vec![100, 100, 100, 100]);
157/// let summary_config = SummaryConfig::new(sample_info.clone())
158///     .with_min_kmer_obs(2)
159///     .with_min_quality(BaseQuality::Medium)
160///     .with_max_p(Some(0.05)) // can also be none which can avoid p-value calculations and save time
161///     .with_group_frac(GroupFrac::One, 0.3);
162/// 
163/// let summary_config2 = SummaryConfig::new(sample_info.clone())
164///     .with_min_kmer_obs(4)
165///     .with_min_quality_for_edge(BaseQuality::Marginal);
166/// 
167/// let summary_config3 = SummaryConfig::new(sample_info)
168///     .with_significant(Some(5))
169///     .with_stat_test(StatTest::StudentsTTest);
170/// ```
171/// 
172/// To filter an already constructed graph, please use the original `SummaryConfig`
173/// and change the settings in plase with the `set_...` methods. This is so the program
174/// knows that some settings have been changed and corresponding values have
175/// to be re-calculated.
176/// 
177/// ```
178/// use debruijn::summarizer::{SummaryConfig, SampleInfo, StatTest};
179/// use debruijn::BaseQuality;
180/// 
181/// let sample_info = SampleInfo::new(0b1100, 0b0011, vec![100, 100, 100, 100]);
182/// let mut summary_config = SummaryConfig::new(sample_info);
183/// 
184/// // ... construct the graph
185/// 
186/// summary_config.set_max_p(Some(0.05));
187/// summary_config.set_stat_test(StatTest::StudentsTTest);
188/// 
189/// // ...
190/// ```
191/// 
192/// By setting the `with_min_kmer_obs` to 0, k-mers for which all occurences were filtered
193/// out with `with_min_quality_for_edge`, can still be included, with empty k-mer data.
194#[derive(Debug, Serialize, Deserialize, Clone, PartialEq)]
195pub struct SummaryConfig {
196    min_kmer_obs: usize,
197    significant: Option<u32>,
198    group_frac: GroupFrac,
199    frac_cutoff: f32,
200    sample_info: SampleInfo,
201    max_p: Option<f32>,
202    stat_test: StatTest,
203    stat_test_changed: bool,
204    min_quality: BaseQuality,
205    min_quality_for_edge: BaseQuality
206}
207
208impl SummaryConfig {
209    /// make a new `SummaryConfig`
210    /// 
211    /// arguments: 
212    /// * `sample_info`: a [`SampleInfo`] with information about the sample groups,
213    ///   which is required for any statistical analysis
214    pub fn new(sample_info: SampleInfo) -> Self {
215        SummaryConfig::empty().with_sample_info(sample_info)
216    }
217
218    /// make an empty `SummaryConfig`. A proper [`SampleInfo`] is required for any
219    /// statistical analysis regarding sample groups.
220    pub fn empty() -> Self {
221        SummaryConfig { 
222            min_kmer_obs: 1, 
223            significant: None, 
224            group_frac: GroupFrac::None, 
225            frac_cutoff: 0., 
226            sample_info: SampleInfo::empty(), 
227            max_p: None, 
228            stat_test: StatTest::WelchsTTest, 
229            stat_test_changed: false,
230            min_quality: BaseQuality::NoCall,
231            min_quality_for_edge: BaseQuality::NoCall,
232        }
233    }
234
235    /// produce a new `SummaryConfig` which will filter k-mers by their number 
236    /// of observations
237    pub fn with_min_kmer_obs(&self, min_kmer_obs: usize) -> Self {
238        let mut config = self.clone();
239        config.min_kmer_obs = min_kmer_obs;
240        config
241    }
242
243    /// modify the number of k-mers observations required for each k-mer to be 
244    /// included in the graph
245    pub fn set_min_kmer_obs(&mut self, min_kmer_obs: usize) {
246        self.min_kmer_obs = min_kmer_obs;
247    }
248
249    /// produce a new `SummaryConfig` which will round the number of observations
250    /// of the k-mer to `significant_digits`
251    pub fn with_significant(&self, significant_digits: Option<u32>) -> Self {
252        let mut config = self.clone();
253        config.significant = significant_digits;
254        config
255    }
256
257    /// modify the number signigicant digits the number of observations stored 
258    /// with the k-mer will be rounded to
259    pub fn set_significant(&mut self, significant_digits: Option<u32>) {
260        self.significant = significant_digits;
261    }
262
263    /// produce a new `SummaryConfig` which will require the k-mer to be ovserved
264    /// in at least a fraction of `frac_cutoff` of either one, both or none of the
265    /// sample groups
266    pub fn with_group_frac(&self, group_frac: GroupFrac, frac_cutoff: f32) -> Self {
267        let mut config = self.clone();
268        config.group_frac = group_frac;
269        config.frac_cutoff = frac_cutoff;
270        config
271    }
272
273    /// modify the group fract settings which require the k-mer to be ovserved
274    /// in at least a fraction of `frac_cutoff` of either one, both or none of the
275    /// sample groups
276    pub fn set_group_frac(&mut self, group_frac: GroupFrac, frac_cutoff: f32) {
277        self.group_frac = group_frac;
278        self.frac_cutoff = frac_cutoff;
279    }
280
281    /// produce a new `SummaryConfig` which contains the given [`SampleInfo`]
282    fn with_sample_info(self, sample_info: SampleInfo) -> Self {
283        let mut config = self.clone();
284        config.sample_info = sample_info;
285        config
286    }
287
288    /// produce a new `SummaryConfig` which will filter k-mers by their p-value
289    /// regarding occurrence in the sample groups
290    pub fn  with_max_p(&self, max_p: Option<f32>) -> Self {
291        let mut config = self.clone();
292        config.max_p = max_p;
293        config
294    }
295
296    /// modify the maximum p-value a k-mer is allowed to have
297    pub fn  set_max_p(&mut self, max_p: Option<f32>) {
298        self.max_p = max_p;
299    }
300
301    /// produce a new `SummaryConfig` which will calculate the p-values based on
302    /// the given statistical test
303    pub fn with_stat_test(&self, stat_test: StatTest) -> Self {
304        let mut config = self.clone();
305        config.stat_test = stat_test;
306        config
307    }
308
309    /// modify the statistical test, which is used to calculate p-values
310    /// regarding observations in the two sample groups
311    pub fn set_stat_test(&mut self, stat_test: StatTest) {
312        if stat_test != self.stat_test { self.stat_test_changed = true }
313        self.stat_test = stat_test;
314    }
315
316    /// produce a new `SummaryConfig` which will filter k-mers by their quality
317    pub fn with_min_quality(&self, min_quality: BaseQuality) -> Self {
318        let mut config = self.clone();
319        config.min_quality = min_quality;
320        config
321    }
322
323    /// modify the minimum quality required for each k-mer to be 
324    /// included in the graph
325    pub fn set_min_quality(&mut self, min_quality: BaseQuality) {
326        self.min_quality = min_quality;
327    }
328
329    /// produce a new `SummaryConfig` which will keep k-mers from being counted 
330    /// if it's quality is too low - should this disconnect the k-mer it will 
331    /// still be counted
332    pub fn with_min_quality_for_edge(&self, min_quality_for_edge: BaseQuality) -> Self {
333        let mut config = self.clone();
334        config.min_quality_for_edge = min_quality_for_edge;
335        config
336    }
337
338    /// modify the minimum quality required for each k-mer observation for its
339    /// edges to be counted
340    pub fn set_min_quality_for_edge(&mut self, min_quality_for_edge: BaseQuality) {
341        self.min_quality_for_edge = min_quality_for_edge;
342    }
343
344    /// get the binary encoded group affiliation of tags
345    pub fn get_markers(&self) -> (Marker, Marker) {
346        self.sample_info.get_markers()
347    }
348
349    /// get the [`SampleInfo`] stored in the `SummaryConfig`
350    pub fn sample_info(&self) -> &SampleInfo {
351        &self.sample_info
352    }
353}
354
355/// In how many of the two sample groups does a specified percentage of the samples 
356/// have to be present
357#[derive(Copy, Clone, PartialEq, PartialOrd, ValueEnum, Debug, Serialize, Deserialize)]
358#[serde(rename_all = "kebab-case")]
359pub enum GroupFrac {
360    None, 
361    One, 
362    Both,
363}
364
365impl std::fmt::Display for GroupFrac {
366    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
367        match &self {
368            GroupFrac::None => write!(f, "none"),
369            GroupFrac::One => write!(f, "one"),
370            GroupFrac::Both => write!(f, "both")            
371        }
372    }
373}
374
375/// Statistical test for calculation of p-values
376#[derive(Copy, Clone, PartialEq, PartialOrd, ValueEnum, Debug, Serialize, Deserialize)]
377#[serde(rename_all = "kebab-case")]
378pub enum StatTest {
379    StudentsTTest,
380    WelchsTTest,
381    UTest,
382}
383
384impl std::fmt::Display for StatTest {
385    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
386        match self {
387            Self::StudentsTTest => write!(f, "students-t-test"),
388            Self::WelchsTTest => write!(f, "welchs-t-test"),
389            Self::UTest => write!(f, "u-test"),
390            
391        }
392    }
393    
394}
395
396/// contains information about the samples required for graph construction
397/// 
398/// ### Example:
399/// 
400/// - Sample IDs in group 1: 0, 1, 2
401/// - Sample IDs in group 2: 3, 4, 5, 6
402/// 
403/// ```
404/// use debruijn::summarizer::{SampleInfo, Marker};
405/// 
406/// let marker0: Marker = 0b0000111; // = 7
407/// let marker1: Marker = 0b1111000; // = 120
408/// 
409/// let sample_kmers = vec![1232, 12323, 24342, 24234, 345456, 21234, 546456];
410/// assert_eq!(marker0.count_ones() + marker1.count_ones(), sample_kmers.len() as u32);
411/// 
412/// let sample_info = SampleInfo::new(marker0, marker1, sample_kmers);
413/// 
414/// ```
415/// 
416#[derive(Clone, Serialize, Deserialize, Debug, PartialEq, Eq, PartialOrd, Ord)]
417pub struct SampleInfo {
418    marker0: Marker,
419    marker1: Marker,
420    count0: u8,
421    count1: u8,
422    sample_kmers: Vec<u64>,
423}
424
425impl SampleInfo {
426    /// make a new [`SampleInfo`]
427    /// 
428    /// ### Arguments
429    /// * `marker0`: a [`M`] which binary-encodes the affiliation of the tags to a group
430    /// * `marker1`: same as `marker0`, for a second group
431    /// * `sample_kmers`: a [`Vec<u64>`] containing numbers of non-unique k-mers for each sample
432    ///   at the index of the sample-id
433    pub fn new(marker0: Marker, marker1: Marker, sample_kmers: Vec<u64>) -> Self {
434        let count0 = marker0.count_ones() as u8;
435        let count1 = marker1.count_ones() as u8;
436        assert_eq!(count0+count1, sample_kmers.len() as u8);
437
438        SampleInfo { marker0, marker1, count0, count1, sample_kmers }
439    }
440
441    /// make a new, empty [`SampleInfo`]
442    pub fn empty() -> Self {
443        SampleInfo { marker0: 0, marker1: 0, count0: 0, count1: 0, sample_kmers: Vec::new() }
444    }
445
446    /// get the binary encoded group affiliation of tags
447    pub fn get_markers(&self) -> (Marker, Marker) {
448        (self.marker0, self.marker1)
449    }
450}
451
452/// count the ocurrences of the tags
453fn tag_counter(tag_vec: &[Tag]) -> Vec<u32> {
454    let mut tag_counter = 1;
455    let mut tag_counts: Vec<u32> = Vec::new();
456
457    // count the occurences of the labels
458    for i in 1..tag_vec.len() {
459        if tag_vec[i] == tag_vec[i-1] {
460            tag_counter += 1;
461        } else {
462            tag_counts.push(tag_counter);
463            tag_counter = 1;
464        }
465    }
466    tag_counts.push(tag_counter);
467    tag_counts.shrink_to_fit();
468
469    tag_counts
470}
471
472#[derive(Debug)]
473struct TagSummary {
474    all_exts: Exts,
475    tag_vec: Vec<Tag>,
476    tag_counts: Vec<u32>,
477    sum: u32,
478    edge_mults: EdgeMult,
479    highest_quality: Option<BaseQuality>,
480}
481
482/// summarize the k-mers, exts and labels, also include an [`EdgeMult`]
483fn summarize_tags<K: Kmer, F: Iterator<Item = KmerDataItem<K, Tag>>>(items: F) -> TagSummary {
484    let mut all_exts = Exts::empty();
485    let mut tag_vec: Vec<Tag> = Vec::with_capacity(items.size_hint().0);
486    let mut edge_mults = EdgeMult::new();
487    let mut highest_quality = None;
488
489    let mut nobs = 0;
490    for item in items {
491        tag_vec.push(item.data); 
492        all_exts = all_exts.add(item.exts);
493        edge_mults.add_exts(item.exts);
494        nobs += 1;
495
496        if let Some(q) = item.quality {
497            if let Some(hq) = highest_quality {
498                if q > hq { highest_quality = Some(q) }
499            } else {
500                highest_quality = Some(q)
501            }
502        }
503    }
504
505    assert_eq!(all_exts, edge_mults.exts());
506
507    tag_vec.sort();
508
509    let tag_counts = tag_counter(&tag_vec);
510
511    tag_vec.dedup();
512
513    TagSummary { all_exts, tag_vec, tag_counts, sum: nobs, edge_mults, highest_quality }
514}
515
516fn summarize_tags_edge_q<K: Kmer, F: Iterator<Item = KmerDataItem<K, Tag>>>(items: F, config: &SummaryConfig) 
517-> TagSummary
518{    
519    // filter the k-mer occurences by their quality -> only use exts and data from k-mers with good enough quality
520    let items_filtered = items.filter(|item| 
521        match item.quality {
522            None => true,
523            Some(q) => q >= config.min_quality_for_edge
524        }
525    );
526    
527    summarize_tags(items_filtered)
528}
529
530#[derive(Debug)]
531struct IDTagSummary {
532    all_exts: Exts,
533    tag_vec: Vec<Tag>,
534    tag_counts: Vec<u32>,
535    sum: u32,
536    id_vec: Vec<ID>,
537    edge_mults: EdgeMult,
538    highest_quality: Option<BaseQuality>,
539}
540
541/// summarize the k-mers, exts and labels
542fn summarize_tags_ids<K: Kmer, F: Iterator<Item = KmerDataItem<K, IDTag>>>(items: F) -> IDTagSummary {
543    let mut all_exts = Exts::empty();
544    let mut tag_vec = Vec::with_capacity(items.size_hint().0);
545    let mut id_vec = Vec::new();
546    let mut edge_mults = EdgeMult::new();
547    let mut highest_quality = None;
548
549    let mut nobs = 0;
550    for item in items {
551        tag_vec.push(item.data.tag); 
552        id_vec.push(item.data.id);
553        all_exts = all_exts.add(item.exts);
554        edge_mults.add_exts(item.exts);
555        nobs += 1;
556
557        if let Some(q) = item.quality {
558            if let Some(hq) = highest_quality {
559                if q > hq { highest_quality = Some(q) }
560            } else {
561                highest_quality = Some(q)
562            }
563        }
564    }
565
566    tag_vec.sort();
567    id_vec.sort();
568
569    let tag_counts = tag_counter(&tag_vec);
570
571    tag_vec.dedup();
572    id_vec.dedup();
573    id_vec.shrink_to_fit();
574
575    IDTagSummary {all_exts, tag_vec, tag_counts, sum: nobs, id_vec, edge_mults, highest_quality}
576}
577
578fn summarize_tags_ids_edge_q<K: Kmer, F: Iterator<Item = KmerDataItem<K, IDTag>>>(items: F, config: &SummaryConfig) 
579-> IDTagSummary
580{
581    // filter the k-mer occurences by their quality -> only use exts and data from k-mers with good enough quality
582    let items_filtered = items.filter(|item|
583        match item.quality {
584            None => true,
585            Some(q) => q >= config.min_quality_for_edge
586        }
587    );
588
589    summarize_tags_ids(items_filtered)
590}
591
592/// round an unsigned integer to the specified amount of digits,
593/// if the integer is shorter than the number if digits, it returns the original integer
594pub fn round_digits(number: u32, digits: u32) -> u32 {
595    let length = (number as f32).log10() as u32 + 1;
596    if digits > length { return number }
597    let empty = length - digits;
598    ((number as f32/ 10i32.pow(empty) as f32).round() * 10i32.pow(empty) as f32) as u32
599}
600
601// check if the k-mer is valid according to the GroupFrac rule and its n obs
602fn valid_counts(tags: Tags, nobs: Option<u32>, config: &SummaryConfig) -> bool {
603    let nobs_valid = match nobs {
604        Some(n) => n as usize >= config.min_kmer_obs,
605        None => true
606    };
607
608    match config.group_frac {
609        GroupFrac::None => nobs_valid,
610        GroupFrac::Both => {
611            // get amount of labels in tags from each group
612            let dist0= tags.bit_and_dist(config.sample_info.marker0);
613            let dist1= tags.bit_and_dist(config.sample_info.marker1);
614    
615            assert_eq!(dist0 + dist1, tags.to_tag_vec().len());
616    
617            // valid if:
618            // - n obs >= min obs AND
619            // - observed in at least one third of samples in both groups
620            nobs_valid
621                && (dist0 as f32 / config.sample_info.count0 as f32 >= config.frac_cutoff) 
622                && (dist1 as f32 / config.sample_info.count1 as f32 >= config.frac_cutoff)
623        },
624        GroupFrac::One => {
625            // get amount of labels in tags from each group
626            let dist0= tags.bit_and_dist(config.sample_info.marker0);
627            let dist1= tags.bit_and_dist(config.sample_info.marker1);
628
629    
630            assert_eq!(dist0 + dist1, tags.to_tag_vec().len());
631    
632            // valid if:
633            // - n obs >= min obs AND
634            // - observed in at least one third of samples in one group
635            nobs_valid
636                && ((dist0 as f32 / config.sample_info.count0 as f32 >= config.frac_cutoff) 
637                    | (dist1 as f32 / config.sample_info.count1 as f32 >= config.frac_cutoff))
638        }
639    }
640}
641
642enum PInfo<'a> {
643    PValue { p: f32 },
644    Calculate { tag_vec: &'a [Tag], tag_counts: &'a Vec<u32> }
645}
646
647fn valid_p(p_info: PInfo, config: &SummaryConfig) -> bool {
648    match config.max_p {
649        Some(max_p) => {
650            match p_info {
651                PInfo::PValue { p } => p <= max_p,
652                PInfo::Calculate { tag_vec, tag_counts } => {
653                    match p_value(tag_vec, tag_counts, config) {
654                        Ok(p) => p <= max_p,
655                        Err(_) => true
656                    } 
657                }
658            }
659        },
660        None => true
661    }
662}
663
664fn p_value(tag_vec: &[Tag], tag_counts: &[u32], config: &SummaryConfig) -> Result<f32, NotEnoughSamplesError> {
665    match config.stat_test {
666        StatTest::StudentsTTest => students_t_test(tag_vec, tag_counts, &config.sample_info),
667        StatTest::WelchsTTest => welchs_t_test(tag_vec, tag_counts, &config.sample_info),
668        StatTest::UTest => u_test(tag_vec, tag_counts, &config.sample_info),
669    }
670}
671
672// perform a student's t-test
673fn students_t_test(tag_vec: &[Tag], tag_counts: &[u32], sample_info: &SampleInfo) -> Result<f32, NotEnoughSamplesError> {
674    let n0 = sample_info.count0 as f64;
675    let n1 = sample_info.count1 as f64;
676
677    if (n0 < 2.) | (n1 < 2.) { return Err(NotEnoughSamplesError {})}
678
679    let mut counts_g0 = Vec::new();
680    let mut counts_g1 = Vec::new();
681
682    let (m0, m1) = sample_info.get_markers();
683
684    for (label, count) in tag_vec.iter().zip(tag_counts) {
685        let bin_rep = (2 as Marker).pow(*label as u32);
686        let norm = *count as f64 / sample_info.sample_kmers[*label as usize] as f64;
687        if (m0 & bin_rep) > 0 { counts_g0.push(norm); }
688        if (m1 & bin_rep) > 0 { counts_g1.push(norm); }
689    }
690
691    let mean0 = counts_g0.iter().sum::<f64>() / n0;
692    let mean1 = counts_g1.iter().sum::<f64>() / n1;
693
694    let df = n0 + n1 - 2.;
695
696    let var0 = (counts_g0.iter().map(|count| (*count - mean0).powi(2)).sum::<f64>() + (n0 - counts_g0.len() as f64) * mean0.powi(2)) / (n0 - 1.);
697    let var1 = (counts_g1.iter().map(|count| (*count - mean1).powi(2)).sum::<f64>() + (n1 - counts_g1.len() as f64) * mean1.powi(2)) / (n1 - 1.);
698
699    let s = ((1./n0 + 1./n1) * ((n0 - 1.) * var0 + (n1 - 1.) * var1) / df).sqrt();
700
701    let t = (mean0 - mean1) / s;
702
703    let t_dist = StudentsT::new(0.0, 1.0, df).expect("error creating student dist: check if you have enough samples (at least 3)");
704
705    let p_value = 2. * (1. - t_dist.cdf(t.abs())) as f32;
706
707    Ok(p_value)
708}
709
710// perform a welch's t-test
711fn welchs_t_test(tag_vec: &[Tag], tag_counts: &[u32], sample_info: &SampleInfo) -> Result<f32, NotEnoughSamplesError> {
712    let n0 = sample_info.count0 as f64;
713    let n1 = sample_info.count1 as f64;
714
715    if (n0 < 2.) | (n1 < 2.) { return Err(NotEnoughSamplesError {})}
716
717    let mut counts_g0 = Vec::new();
718    let mut counts_g1 = Vec::new();
719
720    let (m0, m1) = sample_info.get_markers();
721
722    for (label, count) in tag_vec.iter().zip(tag_counts) {
723        let bin_rep = (2 as Marker).pow(*label as u32);
724        let norm = *count as f64 / sample_info.sample_kmers[*label as usize] as f64;
725        if (m0 & bin_rep) > 0 { counts_g0.push(norm); }
726        if (m1 & bin_rep) > 0 { counts_g1.push(norm); }
727    }
728
729    let mean0 = counts_g0.iter().sum::<f64>() / n0;
730    let mean1 = counts_g1.iter().sum::<f64>() / n1;
731
732    let s0 = (counts_g0.iter().map(|count| (*count - mean0).powi(2)).sum::<f64>() + (n0 - counts_g0.len() as f64) * mean0.powi(2)) / (n0 - 1.);
733    let s1 = (counts_g1.iter().map(|count| (*count - mean1).powi(2)).sum::<f64>() + (n1 - counts_g1.len() as f64) * mean1.powi(2)) / (n1 - 1.);
734
735    let s0 = s0.sqrt();
736    let s1 = s1.sqrt();
737
738    let t = (mean0 - mean1) / (s0.powi(2)/n0 + s1.powi(2)/n1).sqrt();
739
740    let df0 = n0 - 1.;
741    let df1 = n1 - 1.;
742
743    let df = ((s0.powi(2) / n0 + s1.powi(2) / n1).powi(2) / 
744        (s0.powi(4) / (n0.powi(2) * df0) 
745            + s1.powi(4) / (n1.powi(2) * df1))).floor();
746
747    let t_dist = StudentsT::new(0.0, 1.0, df).expect("error creating student dist: check if you have enough samples");
748
749    let p_value = 2. * (1. - t_dist.cdf(t.abs())) as f32;
750
751    Ok(p_value)
752}
753
754// perform a mann-whitney-u-test
755fn u_test(tag_vec: &[Tag], tag_counts: &[u32], sample_info: &SampleInfo) -> Result<f32, NotEnoughSamplesError> {
756
757    let n0 = sample_info.count0 as f64;
758    let n1 = sample_info.count1 as f64;
759
760    if (n0 < 2.) | (n1 < 2.) { return Err(NotEnoughSamplesError {})}
761
762
763    let mut counts_g0 = Vec::new();
764    let mut counts_g1 = Vec::new();
765
766    let (m0, m1) = sample_info.get_markers();
767
768    for (label, count) in tag_vec.iter().zip(tag_counts) {
769        let bin_rep = (2 as Marker).pow(*label as u32);
770        let norm = *count as f64 / sample_info.sample_kmers[*label as usize] as f64;
771        if (m0 & bin_rep) > 0 { counts_g0.push(norm); }
772        if (m1 & bin_rep) > 0 { counts_g1.push(norm); }
773    }
774
775    let n = n0 + n1;
776    let m = (n0 * n1) / 2.;
777
778    // TODO check if more efficient way possible
779
780    let mut all_counts = vec![(0u8, 0f64); n0 as usize - counts_g0.len()];
781    all_counts.append(&mut vec![(1, 0f64); n1 as usize - counts_g1.len()]);
782
783    all_counts.append(&mut counts_g0.iter().map(|elt| (0, *elt)).collect());
784    all_counts.append(&mut counts_g1.iter().map(|elt| (1, *elt)).collect());
785
786    all_counts.sort_by(|(_, x), (_, y)| x.total_cmp(y));
787
788    let chunked = all_counts
789        .chunk_by(|(_, x), (_, y)| x == y);
790
791    let mut ranks = Vec::new();
792    let mut tie_factor = 0.;
793
794    for chunk in chunked {
795        let rank = (chunk.len() + 1) as f64 / 2. + ranks.len() as f64;
796        ranks.append(&mut chunk.iter().map(|(g, _)| (*g, rank)).collect());
797        if !chunk.is_empty() {
798            tie_factor += (chunk.len().pow(3) - chunk.len()) as f64;
799        }
800    }
801
802    let mut rank_sum0 = 0.;
803    let mut rank_sum1 = 0.;
804    ranks.iter().for_each(|(group, rank)| match group { 
805        0 => rank_sum0 += rank, 
806        1 => rank_sum1 += rank, 
807        _ => panic!("should not happen"),
808    });
809
810    let u0 = rank_sum0 - n0 * (n0 + 1.) / 2.;
811    let u1 = rank_sum1 - n1 * (n1 + 1.) / 2.;
812    
813    let u = min_by(u0, u1, |a, b| a.total_cmp(b));
814
815    let s = ((n0 * n1 * (n + 1.) / 12.) - (n0 * n1 * tie_factor / (12. * n * (n - 1.)))).sqrt();
816    //let s = ((n0 * n1 / 12.) * ((n + 1.) - (tie_factor / n * (n - 1.)))).sqrt();
817    // supposedly same term but returns NaN???
818
819    let z = (u - m) / s;
820
821    let dist = Normal::standard();
822    let p_value =  2. * (1. - dist.cdf(z.abs())) as f32;
823
824    Ok(p_value)
825}
826
827// calculate the log2 of the log change of the two groups
828fn log2_fold_change(tags: Tags, counts: &[u32], sample_info: &SampleInfo) -> f32 {
829    let mut norm_count_g0 = 0.;
830    let mut norm_count_g1 = 0.;
831
832    let (m0, m1) = sample_info.get_markers();
833
834    for (label, count) in tags.to_tag_vec().iter().zip(counts) {
835        let bin_rep = (2 as Marker).pow(*label as u32);
836        // normalize with number of k-mers in the sample
837        let norm = *count as f64 / sample_info.sample_kmers[*label as usize] as f64;
838        if (m0 & bin_rep) > 0 { norm_count_g0 += norm; }
839        if (m1 & bin_rep) > 0 { norm_count_g1 += norm; }
840    }
841
842    // normalize with the number of samples in the group
843    norm_count_g0 /= sample_info.count0 as f64;
844    norm_count_g1 /= sample_info.count1 as f64;
845
846    (norm_count_g0 / norm_count_g1).log2() as f32
847}
848/// Trait for summarizing k-mers, determines the data saved in the graph nodes
849pub trait SummaryData<DI>: Clone + Debug + Send + Sync + PartialEq + Serialize + DeserializeOwned {
850    /// format the node data 
851    fn print(&self, translator: &Translator, config: &SummaryConfig, id_group_translator: Option<&HashMap<ID, ID>>) -> String;
852    /// format the node data in for json
853    fn print_ol(&self, translator: &Translator, config: &SummaryConfig, id_group_translator: Option<&HashMap<ID, ID>>) -> String;
854    /// format the node data in one line
855    fn print_json(&self, translator: &Translator, config: &SummaryConfig, id_group_translator: Option<&HashMap<ID, ID>>) -> String;
856    /// get `Tags` and the overall count, returns `None` if data is insufficient
857    fn tags(&self) -> Option<Tags> { None }
858    /// get the size of the structure, including contents of boxed slices
859    fn mem(&self) -> usize;
860    /// get the number of observations, returns `None` if data is insufficient
861    fn sum(&self) -> Option<u32> { None }
862    /// get the IDs, returns `None` if data is insufficient
863    fn ids(&self) -> Option<&[ID]> { None }
864    /// get the p-value, returns `None` if data is insufficient
865    fn p_value(&self, _config: &SummaryConfig) -> Option<f32> { None }
866    /// get the log2(fold change), returns `None` if data is insufficient
867    fn fold_change(&self, _config: &SummaryConfig) -> Option<f32> { None }
868    /// get the number of samples the sequence was observed in, returns `None` if data is insufficient
869    fn sample_count(&self) -> Option<usize> { None }
870    /// get the coverage of the node edges
871    fn edge_mults(&self) -> Option<&EdgeMult> { None }
872    /// get the quality of the node k-mer
873    fn quality(&self) -> Option<BaseQuality> { None }
874    /// fix the [`EdgeMult`] by removing hanging edges
875    fn fix_edge_data(&mut self, _exts: Exts) { }
876    /// set the edge mults
877    fn set_edge_mults(&mut self, _edge_mults: Option<EdgeMult>) { }
878    /// get a reference to the mapped ids,  returns `None` if data is insufficient
879    fn mapped_ids(&self) -> Option<&[ID]> { None }
880    /// add mapped ids to the node data
881    fn set_mapped_ids(&mut self, _mapped_ids: Box<[ID]>) { }
882    /// get a reference to the IDs mapped to the node edges, returns `None` id data is insuffivient
883    fn mapped_edge_ids(&self) -> Option<&EdgeMap> { None }
884    /// add mapped IDs to the node's edges
885    fn set_mapped_edge_ids(&mut self, _mapped_edge_ids: Option<EdgeMap>) { }
886    /// check if the data can be joined into one
887    fn join_test(&self, other: &Self) -> bool { self == other }
888    /// check if node is valid according to: min kmer obs, group fraction, p-value
889    fn valid(&self, _config: &SummaryConfig) -> bool { true }
890    /// summarize k-mers
891    fn summarize<K: Kmer, F: Iterator<Item = KmerDataItem<K, DI>>>(items: F, config: &SummaryConfig) -> (bool, Exts, Self);
892    /// check summerizer kind
893    fn summarizer() -> Summarizers;
894}
895// TODO: move SummaryData::print functionality to Display trait?
896
897/// Number of observations for the k-mer
898impl SummaryData<Tag> for u32 {
899    fn print(&self, _: &Translator, _: &SummaryConfig, _: Option<&HashMap<ID, ID>>) -> String {
900        format!("sum: {}", self)
901    }
902
903    fn print_ol(&self, _: &Translator, _: &SummaryConfig, _: Option<&HashMap<ID, ID>>) -> String {
904        format!("sum: {}", self)
905    }
906
907    fn print_json(&self, _: &Translator, _: &SummaryConfig, _: Option<&HashMap<ID, ID>>) -> String {
908        format!("\"sum\": {}", self)
909    }
910
911    fn mem(&self) -> usize {
912        mem::size_of::<Self>()
913    }
914
915    fn sum(&self) -> Option<u32> {
916        Some(*self)
917    }
918
919    fn valid(&self, config: &SummaryConfig) -> bool {
920        *self >= config.min_kmer_obs as u32
921    }
922
923    fn summarize<K: Kmer, F: Iterator<Item = KmerDataItem<K, Tag>>>(items: F, config: &SummaryConfig) -> (bool, Exts, Self) {
924        let summary = summarize_tags_edge_q(items, config);
925
926        let valid_p = valid_p(PInfo::Calculate { tag_vec: &summary.tag_vec, tag_counts: &summary.tag_counts}, config);
927        let valid_q = if let Some(q) = summary.highest_quality { q >= config.min_quality } else {true };
928
929        let tags = Tags::from_tag_vec(&summary.tag_vec);
930
931        let valid  = valid_counts(tags, Some(summary.sum), config) && valid_p && valid_q;
932
933        let sum = match config.significant {
934            Some(digits) => round_digits(summary.sum, digits),
935            None => summary.sum  
936        };
937
938        (valid, summary.all_exts, sum)
939    }
940
941    fn summarizer() -> Summarizers {
942        Summarizers::Sum
943    }
944}
945
946/// the samples the k-mer was observed with, stored in a Vec -
947/// unlike [`TagsData`], this can hold 256 uniqe sample IDs, but uses 
948/// more memory (min 3x)
949#[derive(Debug, Clone, PartialEq, Serialize, Deserialize, SummaryData)]
950// aligned would be 16 Bytes, packed would be 12 Bytes
951pub struct TagVecData {
952    tag_vec: Vec<Tag>,
953}
954
955/// the IDs the k-mer was observed with and its number of observations
956/// ID could be gene-, read-, or orthogroup-ID
957#[derive(Debug, Clone, PartialEq, Serialize, Deserialize, SummaryData)]
958// aligned would be 16 Bytes, packed would be 12 Bytes
959pub struct IDData {
960    ids: Box<[ID]>,
961}
962
963/// the IDs the k-mer was observed with and its number of observations
964/// ID could be gene-, read-, or orthogroup-ID
965#[derive(Debug, Clone, PartialEq, Serialize, Deserialize, SummaryData)]
966// aligned would be 16 Bytes, packed would be 12 Bytes
967pub struct IDSumData {
968    ids: Box<[ID]>,
969    sum: u32,
970}
971
972/// the tags the k-mer was observed with
973#[derive(Debug, Clone, PartialEq, Serialize, Deserialize, SummaryData)]
974pub struct TagsData {
975    tags: Tags,
976}
977
978/// the tags the k-mer was observed with and its number of observations
979#[derive(Debug, Clone, PartialEq, Serialize, Deserialize, SummaryData)]
980// aligned would be 16 Bytes, packed would be 12 Bytes
981pub struct TagsSumData {
982    tags: Tags,
983    sum: u32,
984}
985
986/// Implementation of [`SummaryData<Tag>`]
987/// 
988/// Contains the tags the k-mer was observed with, how many times it 
989/// was observed with each label, and how many times it was observed overall
990#[derive(Debug, Clone, PartialEq, Serialize, Deserialize, SummaryData)]
991pub struct TagsCountsSumData {
992    tags: Tags,
993    sum: u32,
994    counts: Box<[u32]>,
995}
996
997/// Implementation of [`SummaryData<Tag>`]
998/// 
999/// Contains the tags the k-mer was observed with and how many times it 
1000/// was observed with each label
1001#[derive(Debug, Clone, PartialEq, Serialize, Deserialize, SummaryData)]
1002pub struct TagsCountsData {
1003    tags: Tags,
1004    counts: Box<[u32]>
1005}
1006
1007impl TagsCountsData {
1008    #[inline]
1009    pub fn sum(&self) -> u32 {
1010        self.counts.iter().sum::<u32>()
1011    }
1012}
1013
1014/// Implementation of [`SummaryData<Tag>`]
1015/// 
1016/// Contains the tags the k-mer was observed with, how many times it 
1017/// was observed with each label, and a p-value
1018#[derive(Debug, Clone, PartialEq, Serialize, Deserialize, SummaryData)]
1019pub struct TagsCountsPData {
1020    tags: Tags,
1021    counts: Box<[u32]>,
1022    p_value: f32
1023}
1024
1025impl TagsCountsPData {
1026    #[inline]
1027    pub fn sum(&self) -> u32 {
1028        self.counts.iter().sum::<u32>()
1029    }
1030}
1031
1032/// Implementation of [`SummaryData<Tag>`]
1033/// 
1034/// Contains the tags the k-mer was observed with, how many times it 
1035/// was observed with each label, and the edge multiplicites/coverage
1036#[derive(Debug, Clone, PartialEq, Serialize, Deserialize, SummaryData)]
1037pub struct TagsCountsEMData {
1038    tags: Tags,
1039    counts: Box<[u32]>,
1040    edge_mults: EdgeMult,
1041}
1042
1043impl TagsCountsEMData {
1044    #[inline]
1045    pub fn sum(&self) -> u32 {
1046        self.counts.iter().sum::<u32>()
1047    }
1048}
1049
1050/// Implementation of [`SummaryData<Tag>`]
1051/// 
1052/// Contains the tags the k-mer was observed with, how many times it 
1053/// was observed with each label, a p-value, and the edge multiplicites/coverage
1054#[derive(Debug, Clone, PartialEq, Serialize, Deserialize, SummaryData)]
1055pub struct TagsCountsPEMData {
1056    tags: Tags,
1057    counts: Box<[u32]>,
1058    p_value: f32,
1059    edge_mults: EdgeMult,
1060}
1061
1062impl TagsCountsPEMData {
1063    #[inline]
1064    pub fn sum(&self) -> u32 {
1065        self.counts.iter().sum::<u32>()
1066    }
1067}
1068
1069// Implementation of [`SummaryData<Tag>`]
1070/// 
1071/// Contains the tags the k-mer was observed with, how many times it 
1072/// was observed with each label, a p-value, and the edge multiplicites/coverage
1073#[derive(Debug, Clone, PartialEq, Serialize, Deserialize, SummaryData)]
1074pub struct TagsCountsPEMQualityData {
1075    tags: Tags,
1076    counts: Box<[u32]>,
1077    p_value: f32,
1078    quality: BaseQuality,
1079    edge_mults: EdgeMult,
1080}
1081
1082impl TagsCountsPEMQualityData {
1083    #[inline]
1084    pub fn sum(&self) -> u32 {
1085        self.counts.iter().sum::<u32>()
1086    }
1087}
1088
1089/// Implementation of [`SummaryData<IDTag>`]
1090/// 
1091/// Contains the tags the k-mer was observed with and how many times it 
1092/// was observed with each label
1093#[derive(Debug, Clone, PartialEq, Serialize, Deserialize, SummaryData)]
1094pub struct IDTagsCountsData {
1095    ids: Box<[ID]>,
1096    tags: Tags,
1097    counts: Box<[u32]>
1098}
1099
1100impl IDTagsCountsData {
1101    #[inline]
1102    pub fn sum(&self) -> u32 {
1103        self.counts.iter().sum::<u32>()
1104    }
1105}
1106
1107/// Implementation of [`SummaryData<Tag>`]
1108/// 
1109/// Contains the tags the k-mer was observed with, how many times it 
1110/// was observed with each label, a p-value, and the edge multiplicites/coverage
1111#[derive(Debug, Clone, PartialEq, Serialize, Deserialize, SummaryData)]
1112pub struct IDTagsCountsPEMData {
1113    tags: Tags,
1114    counts: Box<[u32]>,
1115    ids: Box<[ID]>,
1116    p_value: f32,
1117    edge_mults: EdgeMult,
1118}
1119
1120impl IDTagsCountsPEMData {
1121    #[inline]
1122    pub fn sum(&self) -> u32 {
1123        self.counts.iter().sum::<u32>()
1124    }
1125}
1126
1127/// Implementation of [`SummaryData<Tag>`]
1128/// 
1129/// Contains the IDs the k-mer was observed with and the edge multiplicites/coverage
1130#[derive(Debug, Clone, PartialEq, Serialize, Deserialize, SummaryData)]
1131pub struct IDEMData {
1132    ids: Box<[ID]>,
1133    edge_mults: EdgeMult,
1134}
1135
1136/// Implementation of [`SummaryData<Tag>`]
1137/// 
1138/// Contains the IDs the k-mer was observed with, a placeholder for mapped ids, and edge multiplicites/coverage
1139#[derive(Debug, Clone, PartialEq, Serialize, Deserialize, SummaryData)]
1140pub struct IDMapEMData {
1141    ids: Box<[ID]>,
1142    map_ids: Box<[ID]>,
1143    edge_mults: EdgeMult,
1144}
1145
1146/// Implementation of [`SummaryData<Tag>`]
1147/// 
1148/// Contains the IDs the k-mer was observed with, a placeholder for mapped ids, and edge multiplicites/coverage
1149#[derive(Debug, Clone, PartialEq, Serialize, Deserialize, SummaryData)]
1150pub struct IDMapEMQualityData {
1151    ids: Box<[ID]>,
1152    map_ids: Box<[ID]>,
1153    quality: BaseQuality,
1154    edge_mults: EdgeMult,
1155    
1156}
1157
1158/// Implementation of [`SummaryData<Tag>`]
1159/// 
1160/// Contains the IDs the k-mer was observed with, a placeholder for mapped ids, and edge multiplicites/coverage
1161#[derive(Debug, Clone, PartialEq, Serialize, Deserialize, SummaryData)]
1162pub struct SumMapEMQualityData {
1163    sum: u32,
1164    map_ids: Box<[ID]>,
1165    quality: BaseQuality,
1166    edge_mults: EdgeMult,
1167}
1168
1169/// Implementation of [`SummaryData<Tag>`]
1170/// 
1171/// Contains the IDs the k-mer was observed with, a placeholder for mapped ids, and edge multiplicites/coverage
1172#[derive(Debug, Clone, PartialEq, Serialize, Deserialize, SummaryData)]
1173pub struct MapEMEmapQualityData {
1174    map_ids: Box<[ID]>,
1175    edge_mults: EdgeMult,
1176    edge_maps: EdgeMap,
1177    quality: BaseQuality,
1178}
1179
1180/// Implementation of [`SummaryData<Tag>`]
1181/// 
1182/// Contains how many times the k-mer was observed in each group, only validates count
1183#[derive(Debug, Clone, PartialEq, Serialize, Deserialize, SummaryData)]
1184pub struct GroupCountData {
1185    group1: u32,
1186    group2: u32,
1187}
1188
1189impl GroupCountData {
1190    #[inline(always)]
1191    fn sum(&self) -> u32 {
1192        self.group1 + self.group2
1193    }
1194}
1195
1196/// Implementation of [`SummaryData<Tag>`]
1197/// 
1198/// Contains the relative number of observations for the k-mer (in percent) 
1199/// in group 1 and the absolute overall count, only validates count
1200#[derive(Debug, Clone, PartialEq, Serialize, Deserialize, SummaryData)]
1201pub struct RelCountData {
1202    percent: u32,
1203    sum: u32
1204}
1205
1206#[derive(Serialize, Deserialize, Debug, Clone, Copy, PartialEq, Eq)]
1207pub enum Summarizers {
1208    Sum,
1209    TagVecData,
1210    IDData,
1211    IDSumData,
1212    TagsData,
1213    TagsSumData,
1214    TagsCountsData,
1215    TagsCountsSumData,
1216    TagsCountsPData,
1217    TagsCountsEMData,
1218    TagsCountsPEMData,
1219    TagsCountsPEMQualityData,
1220    IDTagsCountsData,
1221    IDTagsCountsPEMData,
1222    IDEMData,
1223    IDMapEMData,
1224    IDMapEMQualityData,
1225    SumMapEMQualityData,
1226    MapEMEmapQualityData,
1227    GroupCountData,
1228    RelCountData
1229}
1230
1231#[cfg(test)]
1232mod test {
1233    
1234
1235    use bimap::BiMap;
1236    use crate::{Tags, clean_graph::CleanGraph, compression::{ CheckCompress, ScmapCompress, compress_graph, compress_kmers_with_hash}, dna_string::DnaString, filter::filter_kmers, graph::{BaseGraph, Node}, kmer::{Kmer8, Kmer16}, reads::{Reads, ReadsPaired}, summarizer::{self, ID, IDData, IDTag, NotEnoughSamplesError, SampleInfo, SummaryData, TagsData, Translator, id_format, p_value, students_t_test, u_test, valid_p, welchs_t_test}};
1237    use crate::Exts;
1238
1239    use super::{log2_fold_change, round_digits, SummaryConfig, TagsCountsSumData};
1240
1241    #[test]
1242    fn test_p_value() -> Result<(), NotEnoughSamplesError> {
1243
1244        /*
1245        group 1: 111111100000 = 4064
1246        group 2: 000000011111 = 31
1247         */
1248
1249        let sample_kmers = vec![1; 12];
1250        let sample_info = SampleInfo::new(31, 4064, sample_kmers);
1251
1252        let summary_config_w = SummaryConfig::new(sample_info.clone()).with_stat_test(summarizer::StatTest::WelchsTTest);
1253        let summary_config_t = SummaryConfig::new(sample_info.clone()).with_stat_test(summarizer::StatTest::StudentsTTest);
1254        let summary_config_u = SummaryConfig::new(sample_info.clone()).with_stat_test(summarizer::StatTest::UTest);
1255        
1256        let tag_vec = [0, 1, 2, 3, 4, 8];
1257        let tag_counts = vec![1; 6];
1258        
1259        let p = welchs_t_test(&tag_vec, &tag_counts, &sample_info)?;
1260        assert_eq!((p * 1000.).round(), 1.);
1261        let p = students_t_test(&tag_vec, &tag_counts, &sample_info)?;
1262        assert_eq!((p * 10000.).round(), 5.);
1263        let p = u_test(&tag_vec, &tag_counts, &sample_info)?;
1264        assert_eq!((p * 1000.).round(), 5.);
1265
1266        let p: f32 = p_value(&tag_vec, &tag_counts, &summary_config_w)?;
1267        assert_eq!((p * 1000.).round(), 1.);
1268        let p = p_value(&tag_vec, &tag_counts, &summary_config_t)?;
1269        assert_eq!((p * 10000.).round(), 5.);
1270        let p = p_value(&tag_vec, &tag_counts, &summary_config_u)?;
1271        assert_eq!((p * 1000.).round(), 5.);
1272
1273
1274        // test with different kmer counts
1275        let sample_kmers = vec![12, 3345, 3478, 87, 1, 2, 666, 98111, 23982938, 555, 122, 7238];
1276
1277        let sample_info = SampleInfo::new(31, 4064, sample_kmers);
1278        let summary_config_w = SummaryConfig::new(sample_info.clone()).with_stat_test(summarizer::StatTest::WelchsTTest);
1279        let summary_config_t = SummaryConfig::new(sample_info.clone()).with_stat_test(summarizer::StatTest::StudentsTTest);
1280        let summary_config_u = SummaryConfig::new(sample_info.clone()).with_stat_test(summarizer::StatTest::UTest);
1281
1282        let tag_vec = [0, 1, 7, 8, 9, 10];
1283        let tag_counts = vec![1; 6];
1284
1285        let p = welchs_t_test(&tag_vec, &tag_counts, &sample_info)?;
1286        assert_eq!((p * 10000.).round(), 4113.);
1287        let p = students_t_test(&tag_vec, &tag_counts, &sample_info)?;
1288        assert_eq!((p * 10000.).round(), 2955.);
1289        let p = u_test(&tag_vec, &tag_counts, &sample_info)?;
1290        assert_eq!((p * 1000.).round(), 862.);
1291
1292        let p: f32 = p_value(&tag_vec, &tag_counts, &summary_config_w)?;
1293        assert_eq!((p * 10000.).round(), 4113.);
1294        let p: f32 = p_value(&tag_vec, &tag_counts, &summary_config_t)?;
1295        assert_eq!((p * 10000.).round(), 2955.);
1296        let p: f32 = p_value(&tag_vec, &tag_counts, &summary_config_u)?;
1297        assert_eq!((p * 1000.).round(), 862.);
1298
1299
1300        // test: not enough samples for statistical test
1301
1302        /*
1303        group 1: 000000100 = 4
1304        group 2: 000000011 = 3
1305         */
1306
1307        let sample_kmers = vec![3, 3, 3];
1308        let sample_info = SampleInfo::new(0b100, 0b11, sample_kmers);
1309        let summary_config_w = SummaryConfig::new(sample_info.clone()).with_stat_test(summarizer::StatTest::WelchsTTest);
1310        let summary_config_t = SummaryConfig::new(sample_info.clone()).with_stat_test(summarizer::StatTest::StudentsTTest);
1311        let summary_config_u = SummaryConfig::new(sample_info.clone()).with_stat_test(summarizer::StatTest::UTest);
1312
1313        let tag_vec = [0, 1, 2];
1314        let tag_counts = vec![2, 3, 1];
1315
1316        if welchs_t_test(&tag_vec, &tag_counts, &sample_info).is_ok() { panic!("should throw err") }
1317        if students_t_test(&tag_vec, &tag_counts, &sample_info).is_ok() { panic!("should throw err") }
1318        if u_test(&tag_vec, &tag_counts, &sample_info).is_ok() { panic!("should throw err") }
1319
1320        if p_value(&tag_vec, &tag_counts, &summary_config_w).is_ok() { panic!("should throw err") }
1321        if p_value(&tag_vec, &tag_counts, &summary_config_t).is_ok() { panic!("should throw err") }
1322        if p_value(&tag_vec, &tag_counts, &summary_config_u).is_ok() { panic!("should throw err") }
1323
1324        Ok(())
1325    }
1326
1327    #[test]
1328    fn test_valid_p() {
1329        let sample_kmers = vec![1; 12];
1330        let sample_info = SampleInfo::new(31, 4064, sample_kmers);
1331        let summary_config_m = SummaryConfig::new(sample_info.clone()).with_max_p(None);
1332        let summary_config_p = SummaryConfig::new(sample_info.clone()).with_max_p(Some(0.05));
1333
1334        let tag_vec = [0, 1, 2, 3, 4, 8];  // p should be 0.001
1335        let tag_counts = vec![1; 6];
1336        let vp = valid_p(summarizer::PInfo::Calculate { tag_vec: &tag_vec, tag_counts: &tag_counts }, &summary_config_m);
1337        assert!(vp);
1338        let vp = valid_p(summarizer::PInfo::Calculate { tag_vec: &tag_vec, tag_counts: &tag_counts }, &summary_config_p);
1339        assert!(vp);
1340
1341
1342        let tag_vec = [0, 7, 8, 9, 10]; // p should be 0.2238
1343        let tag_counts = vec![1; 5];
1344
1345        let vp = valid_p(summarizer::PInfo::Calculate { tag_vec: &tag_vec, tag_counts: &tag_counts }, &summary_config_m);
1346        assert!(vp);
1347        let vp = valid_p(summarizer::PInfo::Calculate { tag_vec: &tag_vec, tag_counts: &tag_counts }, &summary_config_p);
1348        assert!(!vp);
1349    }
1350
1351    #[test]
1352    #[cfg(not(feature = "sample128"))]
1353    fn test_data_valid() {
1354        use crate::{kmer::Kmer16, test::build_test_graph};
1355
1356        let mut graph: BaseGraph<Kmer8, TagsCountsSumData> = BaseGraph::new(false);
1357
1358        let tags = Tags::from_tag_vec(&vec![0, 2, 6]);
1359        let counts: Box<[u32]> = [1, 3, 5].into();
1360        let sum = counts.iter().sum::<u32>();
1361        graph.add(&DnaString::from_acgt_bytes("AAAAAAAA".as_bytes()), Exts::empty(), TagsCountsSumData { tags, counts, sum });
1362
1363        let tags = Tags::from_tag_vec(&vec![0]);
1364        let counts: Box<[u32]> = [1].into();
1365        let sum = counts.iter().sum::<u32>();
1366        graph.add(&DnaString::from_acgt_bytes("CCCCCCCC".as_bytes()), Exts::empty(), TagsCountsSumData { tags, counts, sum });
1367                
1368        
1369        let graph = graph.finish();
1370
1371        graph.print();
1372
1373        let sample_kmers = vec![123, 234, 12334, 34, 1232, 123, 123, 34];
1374        let sample_info = SampleInfo::new(0b00100101, 0b11011010, sample_kmers);
1375        let config = SummaryConfig::new(sample_info).with_stat_test(summarizer::StatTest::StudentsTTest);
1376
1377        let censor_nodes = CleanGraph::new(|node: &Node<'_, Kmer8, TagsCountsSumData>| !node.data().valid(&config))
1378                    .find_bad_nodes(&graph);
1379        println!("censor nodes: {:?}", censor_nodes);
1380        let filter_graph = compress_graph(false, &ScmapCompress::new(), graph, Some(censor_nodes));
1381
1382        filter_graph.print();
1383
1384        // larger test
1385
1386        let (_, _, ser_graph) = build_test_graph::<Kmer16, TagsCountsSumData, _>();
1387        let (graph, _translator, mut config) = ser_graph.dissolve();
1388
1389        config.set_min_kmer_obs(3);
1390
1391        let node3 = graph.get_node(3);
1392        assert!(!node3.data().valid(&config));
1393
1394        let bad_nodes = graph.find_bad_nodes(|node| node.data().valid(&config));
1395        let bad_node_correct = vec![3, 15, 20, 21, 23, 35, 36, 37, 40, 41, 46, 58, 61, 66, 67, 76, 77, 80, 89, 98, 99];
1396        assert_eq!(bad_nodes, bad_node_correct);
1397        let _filtered_graph = compress_graph(false, &ScmapCompress::new(), graph, Some(bad_nodes));
1398
1399    }
1400
1401    #[test]
1402    fn test_fold_change() {    
1403        let marker0 = 0b0000000001111;
1404        let marker1 = 0b1111111110000;
1405        let sample_kmers = vec![2834, 2343, 12, 1234, 345345, 122, 234, 23455, 231, 2, 3564, 12344, 34555];
1406        let sample_info = SampleInfo::new(marker0, marker1, sample_kmers);
1407        //let summary_config = SummaryConfig::new(1, None, GroupFrac::None, 0.33, sample_info.clone(), None, summarizer::StatTest::WelchsTTest);
1408
1409        let labels = vec![0, 1, 2, 3, 7, 8];
1410        let tags = Tags::from_tag_vec(&labels);
1411        let counts = vec![1, 6, 9, 3, 6, 10];
1412        let fold_change = log2_fold_change(tags, &counts, &sample_info);
1413        assert_eq!(fold_change, 5.286_453_2);
1414
1415        let labels = vec![0, 6, 7, 8, 10, 11];
1416        let tags = Tags::from_tag_vec(&labels);
1417        let counts = vec![12, 3, 7, 1, 22, 6];
1418        let fold_change = log2_fold_change(tags, &counts, &sample_info);
1419        assert_eq!(fold_change, -1.339_324_5);
1420
1421        // x/0 = inf -> log(inf) = inf
1422        let labels = vec![0, 1];
1423        let tags = Tags::from_tag_vec(&labels);
1424        let counts = vec![12, 3];
1425        let fold_change = log2_fold_change(tags, &counts, &sample_info);
1426        assert_eq!(fold_change, f32::INFINITY);
1427
1428        // 0/x = 0 -> log2(0) = -inf
1429        let labels = vec![7, 8];
1430        let tags = Tags::from_tag_vec(&labels);
1431        let counts = vec![12, 3];
1432        let fold_change = log2_fold_change(tags, &counts, &sample_info);
1433        assert_eq!(fold_change, f32::NEG_INFINITY);       
1434    }
1435
1436    #[test]
1437    fn test_round_digits() {
1438        assert_eq!(round_digits(18293092, 4), 18290000);
1439        assert_eq!(round_digits(333, 4), 333);
1440        assert_eq!(round_digits(129552, 2), 130000);
1441        assert_eq!(round_digits(1829399, 6), 1829400);
1442        assert_eq!(round_digits(1829399, 0), 0);
1443    }
1444
1445
1446    #[test]
1447    fn test_id_format() {
1448        let id_translator = [("A", 0), ("B", 1), ("C", 2), ("D", 3), ("E", 4), ("F", 5), ("G", 6)].into_iter().map(|(a, b)| (a.to_string(), b as ID)).collect();
1449        let id_gr_tr = [(0 as ID, 0 as ID), (1, 0), (2, 0), (3, 0), (4, 1), (5, 1), (6, 1)].into_iter().collect();
1450
1451        let translator = Translator::new_id_translator(id_translator);
1452        let e_tr = Translator::new_tag_translator(BiMap::new());
1453
1454        assert_eq!("[\"A\", \"B\", \"C\"]", id_format(&[0, 1, 2], &translator, None));
1455        assert_eq!("[0, 1]", id_format(&[0, 1, 5], &translator, Some(&id_gr_tr)));
1456        assert_eq!("[0, 1, 2]", id_format(&[0, 1, 2], &e_tr, None));
1457    }
1458
1459    #[test]
1460    fn test_summary_config() {
1461        let sample_info = SampleInfo::new(0b11, 0b1100, vec![12, 12, 12, 12]);
1462        let config1 = SummaryConfig::new(sample_info.clone())
1463            .with_group_frac(summarizer::GroupFrac::One, 0.3)
1464            .with_max_p(Some(0.3))
1465            .with_min_kmer_obs(2)
1466            .with_min_quality(crate::BaseQuality::Marginal)
1467            .with_min_quality_for_edge(crate::BaseQuality::Medium)
1468            .with_significant(Some(4))
1469            .with_stat_test(summarizer::StatTest::StudentsTTest);
1470
1471        let config2 = SummaryConfig {
1472            min_kmer_obs: 2,
1473            significant: Some(4),
1474            group_frac: summarizer::GroupFrac::One,
1475            frac_cutoff: 0.3,
1476            sample_info,
1477            max_p: Some(0.3),
1478            stat_test: summarizer::StatTest::StudentsTTest,
1479            stat_test_changed: false,
1480            min_quality: crate::BaseQuality::Marginal,
1481            min_quality_for_edge: crate::BaseQuality::Medium
1482        };
1483
1484        assert_eq!(config1, config2)
1485    }
1486
1487    #[test]
1488    fn test_summarize_edge_q() {
1489        let read1 = "AGCTAGCGATGCTAGCTAGCATCGTAGCTAGCAAGCTGATCAAGTCGATGCTGACTGATGCTAGCTGACTGATCGATGCTAGCTGATC";
1490        let qual1 = "CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC";
1491        let read2 = "GAAATCGTAGCTGAAAAATCGTAGCTGATCGTAGCTGATCTATTCTAGCTGATCAAGTCGATCGGAGGGGTTTCGGAGTTTCGGGATTCGTAT";
1492        let qual2 = "CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC-CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC";
1493
1494        /*
1495        should produce three sequences, not connected to each other:
1496        -> the first read as a whole and the second read, with the false k-mer missing and producing a gap
1497        seq: AGCTAGCGATGCTAGCTAGCATCGTAGCTAGCAAGCTGATCAAGTCGATGCTGACTGATGCTAGCTGACTGATCGATGCTAGCTGATC, node: Node { id:0, Exts: |, L:[] R:[], Seq: 88, Data: TagsData { tags: [0] } }
1498        seq: GAAATCGTAGCTGAAAAATCGTAGCTGATCGTAGCTGATCTATTCTAGCTGATCAAGTCGA, node: Node { id:1, Exts: |, L:[] R:[], Seq: 61, Data: TagsData { tags: [1] } }
1499        seq: GCTGATCAAGTCGATCGGAGGGGTTTCGGAGTTTCGGGATTCGTAT, node: Node { id:2, Exts: |, L:[] R:[], Seq: 46, Data: TagsData { tags: [1] } }
1500        [0]
1501        [1]
1502        [2]
1503        */
1504
1505        // without IDs
1506
1507        let mut reads = Reads::new_with_quality(crate::reads::Strandedness::Forward);
1508        reads.add_read(DnaString::from_acgt_bytes(read1.as_bytes()), None, 0, Some(qual1.as_bytes()));
1509        reads.add_read(DnaString::from_acgt_bytes(read2.as_bytes()), None, 1, Some(qual2.as_bytes()));
1510
1511        let reads_paired = ReadsPaired::Unpaired { reads };
1512
1513        let sample_info = SampleInfo::new(0b10, 0b01, vec![73, 78]);
1514        let summary_config = SummaryConfig::new(sample_info)
1515            .with_min_quality_for_edge(crate::BaseQuality::Medium);
1516
1517        let (kmers, _) = filter_kmers::<TagsData, Kmer16, _>(&reads_paired, &summary_config, false, 1., false);
1518
1519        let comp_spec = CheckCompress::new(|a: TagsData, _b| a, |a, b| a.join_test(b));
1520        let mut graph = compress_kmers_with_hash(true, &comp_spec, kmers, false, false).finish();
1521        graph.fix_exts(None);
1522
1523        for node_id in 0..graph.len() {
1524            let node = graph.get_node(node_id);
1525            println!{"seq: {:?}, node: {:?}", node.sequence(), node}
1526        }
1527
1528        let expected_components = [vec![0], vec![1], vec![2]];
1529
1530        for (i, component) in graph.iter_components().enumerate() {
1531            println!("{:?}", component);
1532            assert_eq!(component, expected_components[i]);
1533        }
1534
1535        // with IDs
1536
1537        let mut reads = Reads::new_with_quality(crate::reads::Strandedness::Forward);
1538        reads.add_read(DnaString::from_acgt_bytes(read1.as_bytes()), None, IDTag::new(0, 0), Some(qual1.as_bytes()));
1539        reads.add_read(DnaString::from_acgt_bytes(read2.as_bytes()), None, IDTag::new(1, 1), Some(qual2.as_bytes()));
1540
1541        let reads_paired = ReadsPaired::Unpaired { reads };
1542
1543        let (kmers, _) = filter_kmers::<IDData, Kmer16, _>(&reads_paired, &summary_config, false, 1., false);
1544
1545        let comp_spec = CheckCompress::new(|a: IDData, _b| a, |a, b| a.join_test(b));
1546        let mut graph = compress_kmers_with_hash(true, &comp_spec, kmers, false, false).finish();
1547        graph.fix_exts(None);
1548
1549        for node_id in 0..graph.len() {
1550            let node = graph.get_node(node_id);
1551            println!{"seq: {:?}, node: {:?}", node.sequence(), node}
1552        }
1553
1554        let expected_components = [vec![0], vec![1], vec![2]];
1555
1556        for (i, component) in graph.iter_components().enumerate() {
1557            println!("{:?}", component);
1558            assert_eq!(component, expected_components[i]);
1559        }
1560
1561    }
1562}
debruijn/summarizer.rs

debruijn/
summarizer.rs
