data-beans 0.7.1

Sparse genomics data backends, QC, algorithms, and simulation
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122

use super::*;

fn names(v: &[&str]) -> Vec<Box<str>> {
    v.iter().map(|s| Box::from(*s)).collect()
}

#[test]
fn gene_index_resolves_hgnc_renames_both_ways() {
    // A matrix on the CURRENT HGNC release; a marker panel on an older one. Without the
    // alias tier every one of these silently drops out of the panel.
    let dict = names(&[
        "H4C3", "H2BC21", "H1-0", "MARCHF2", "SEPTIN7", "MTARC1", "CD8A",
    ]);
    let idx = GeneIndex::build(&dict);
    assert_eq!(idx.match_gene("HIST1H4C"), Some(0));
    assert_eq!(idx.match_gene("HIST2H2BE"), Some(1));
    assert_eq!(idx.match_gene("H1F0"), Some(2));
    assert_eq!(idx.match_gene("MARCH2"), Some(3), "rule-based MARCH family");
    assert_eq!(idx.match_gene("SEPT7"), Some(4), "rule-based SEPT family");
    assert_eq!(
        idx.match_gene("MARC1"),
        Some(5),
        "MARC1 is MTARC1, not a MARCH-family member"
    );

    // ...and the reverse: an old-release matrix, a current panel.
    let old = names(&["HIST1H4C", "MARCH2", "SEPT7", "CD8A"]);
    let idx_old = GeneIndex::build(&old);
    assert_eq!(idx_old.match_gene("H4C3"), Some(0));
    assert_eq!(idx_old.match_gene("MARCHF2"), Some(1));
    assert_eq!(idx_old.match_gene("SEPTIN7"), Some(2));

    // The alias tier resolves through the faba `{gene}/count/{track}` row keys and the
    // `ENSG…_SYMBOL` form too, since it runs on the extracted symbol.
    let faba = names(&["ENSG00000197061_H4C3/count/spliced", "CD8A/count/spliced"]);
    let idx_faba = GeneIndex::build(&faba);
    assert_eq!(idx_faba.match_gene("HIST1H4C"), Some(0));

    // A gene with no alias must not be coerced into one.
    assert_eq!(idx.match_gene("ZZZ9"), None);
}

#[test]
fn gene_index_tiers_and_idf() {
    let dict = names(&["ENSG00000153563_CD8A", "MS4A1", "A_FOO", "FOO"]);
    let idx = GeneIndex::build(&dict);

    // exact (case-insensitive)
    assert_eq!(idx.match_gene("ms4a1"), Some(1));
    // symbol = last `_`-segment
    assert_eq!(idx.match_gene("CD8A"), Some(0));
    // exact preferred over an earlier-indexed suffix match (the refinement)
    assert_eq!(idx.match_gene("FOO"), Some(3));
    // unmatched
    assert_eq!(idx.match_gene("ZZZ9"), None);

    // bare ENSG query resolves to the `ENSG…_SYMBOL` row (and the combined
    // form resolves too) — HGNC / ENSG / ENSG_HGNC reconcile both ways.
    assert_eq!(idx.match_gene("ENSG00000153563"), Some(0));
    assert_eq!(idx.match_gene("ENSG00000153563_CD8A"), Some(0));

    // dict keyed by bare ENSG: a combined query still resolves via the
    // leading-id tier.
    let dict2 = names(&["ENSG00000153563", "MS4A1"]);
    let idx2 = GeneIndex::build(&dict2);
    assert_eq!(idx2.match_gene("ENSG00000153563_CD8A"), Some(0));
    assert_eq!(idx2.match_gene("ensg00000153563"), Some(0));

    // IDF: ubiquitous (df == C) → 0; exclusive → ln(C)
    assert_eq!(idf_weight(4, 4), 0.0);
    assert!(idf_weight(4, 1) > 1.38 && idf_weight(4, 1) < 1.39);
}

#[test]
fn a_row_is_named_by_id_alone_when_its_name_adds_nothing() {
    assert_eq!(id_name("FID1", "GENE1"), Box::from("FID1_GENE1"));
    assert_eq!(id_name("FID2", "FID2"), Box::from("FID2"));
    assert_eq!(id_name("FID3", ""), Box::from("FID3"));
}

#[test]
fn loci_match_loci_strictly_by_key() {
    let dict = names(&["X:0-100", "chr1:0-100", "GENE1", "ENSG000_GENE2"]);
    let idx = GeneIndex::build(&dict);
    assert_eq!(idx.match_gene("chrX:0-100"), Some(0));
    assert_eq!(idx.match_gene("CHR1:0-100"), Some(1));
    assert_eq!(idx.match_gene("x:0-100"), None, "chromosome case is kept");
    assert_eq!(idx.match_gene("X:0-101"), None);
    assert_eq!(idx.match_gene("GENE2"), Some(3));
    // A gene query never lands on a locus row through the fallback tiers.
    assert_eq!(idx.match_gene("X"), None);
    assert_eq!(idx.match_gene("chr1"), None);
    // A locus query still reaches a row whose `/`-core is that locus.
    let idx = GeneIndex::build(&names(&["GENE1", "chr2:5-9/count/spliced"]));
    assert_eq!(idx.match_gene("chr2:5-9"), Some(1));
    assert_eq!(idx.match_gene("chr2:5-9/count/spliced"), Some(1));
}

#[test]
fn peak_rows_import_in_colon_form_and_gene_rows_stay() {
    let mut ids = names(&["chr1-100-200", "GENE_1_2", "chrUn_CTG1v1_5_10", "chr2:1-9"]);
    let mut nm = names(&["chr1-100-200", "GENE_1_2", "", "chr2:1-9"]);
    let types = names(&["Peaks", "Gene Expression", "Peaks", "Peaks"]);
    colon_peak_names(&mut ids, &mut nm, Some(&types));
    assert_eq!(
        ids,
        names(&["chr1:100-200", "GENE_1_2", "chrUn_CTG1v1:5-10", "chr2:1-9"])
    );
    assert_eq!(nm, names(&["chr1:100-200", "GENE_1_2", "", "chr2:1-9"]));
}

#[test]
fn an_untyped_list_converts_only_when_every_row_is_an_interval() {
    let mut ids = names(&["chr1-100-200", "chr2_5_9"]);
    let mut nm = names(&["", ""]);
    assert_eq!(colon_peak_names(&mut ids, &mut nm, None), 2);
    assert_eq!(ids, names(&["chr1:100-200", "chr2:5-9"]));
    let mut ids = names(&["chr1-100-200", "GENE1"]);
    let mut nm = names(&["", ""]);
    assert_eq!(colon_peak_names(&mut ids, &mut nm, None), 0);
    assert_eq!(ids, names(&["chr1-100-200", "GENE1"]));
}