bigsig 0.1.0

Large-scale Sequence Search with BItsliced Genomic Signature Index (BIGSIG)
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

use bit_vec_serde::BitVec;
use xxh3;
use fnv;
use kmer;

pub fn batch_search(
    files: Vec<&str>,
    bigsi_map: &fnv::FnvHashMap<usize, BitVec>,
    colors_accession: &fnv::FnvHashMap<usize, String>,
    _n_ref_kmers: &fnv::FnvHashMap<String, usize>,
    bloom_size: usize,
    num_hash: usize,
    k_size: usize,
    _cov: f64,
) {
    for file in files {
        //only fasta formatted file!
        eprintln!("Counting k-mers, this may take a while!");
        let vec_query = kmer::read_fasta(file.to_owned());
        let kmers_query = kmer::kmerize_vector(&vec_query, k_size, 1);
        eprintln!("{} kmers in query", kmers_query.len());
        if kmers_query.len() == 0 {
            eprintln!("Warning! no kmers in query; maybe your kmer length is larger than your query length?");
        } else {
            let mut kmer_slices = Vec::new();
            for k in kmers_query.keys() {
                for i in 0..num_hash {
                    let bit_index =
                    xxh3::hash64_with_seed(&k.as_bytes(), i as u64) % bloom_size as u64;
                    let bi = bit_index as usize;
                    if !bigsi_map.contains_key(&bi) {
                        break;
                    } else {
                        kmer_slices.push(&bigsi_map[&bi]);
                    }
                }
            }
            if kmer_slices.len() < (num_hash * kmers_query.len()) {
                eprintln!("No perfect hits!");
            } else {
                //bit-wise AND
                let mut first = kmer_slices[0].to_owned();
                for i in 1..(num_hash * kmers_query.len()) {
                    let j = i as usize;
                    first.intersect(&kmer_slices[j]);
                }
                let mut hits = Vec::new();
                for i in 0..first.len() {
                    if first[i] {
                        hits.push(&colors_accession[&i]);
                    }
                }
                eprintln!("{} hits", hits.len());
                for h in &hits {
                    println!("{}\t{}\t{}\t1.00", file, h, kmers_query.len());
                }
            }
        }
    }
}

pub fn batch_search_mf(
    files: Vec<&str>,
    bigsi_map: &fnv::FnvHashMap<usize, BitVec>,
    colors_accession: &fnv::FnvHashMap<usize, String>,
    _n_ref_kmers: &fnv::FnvHashMap<String, usize>,
    bloom_size: usize,
    num_hash: usize,
    k_size: usize,
    _cov: f64,
) {
    for file in files {
        let (labels, sequences) = super::kmer::read_fasta_mf(file.to_owned());
        for (i, label) in labels.iter().enumerate() {
            //only fasta formatted file!
            //eprintln!("Counting k-mers, this may take a while!");
            let kmers = super::kmer::kmerize_string(sequences[i].to_owned(), k_size);
            //let kmers_query = kmer_fa::clean_map(unfiltered, filter);
            match kmers{
                None    => println!("Warning! no kmers in query '{}'; maybe your kmer length is larger than your query length?", label),
                Some(kmers_query) =>
                {eprintln!("{} kmers in query", kmers_query.len());
                let mut kmer_slices = Vec::new();
                for k in kmers_query.keys() {
                    for i in 0..num_hash {
                        let bit_index = xxh3::hash64_with_seed(&k.as_bytes(), i as u64)
                            % bloom_size as u64;
                        let bi = bit_index as usize;
                        if !bigsi_map.contains_key(&bi) {
                            break;
                        } else {
                            kmer_slices.push(&bigsi_map[&bi]);
                        }
                    }
                }
                if kmer_slices.len() < (num_hash * kmers_query.len()) {
                    eprintln!("No perfect hits!");
                } else {
                    //bit-wise AND
                    let mut first = kmer_slices[0].to_owned();
                    for i in 1..(num_hash * kmers_query.len()) {
                        let j = i as usize;
                        first.intersect(&kmer_slices[j]);
                    }
                    let mut hits = Vec::new();
                    for i in 0..first.len() {
                        if first[i] {
                            hits.push(&colors_accession[&i]);
                        }
                    }
                    eprintln!("{} hits", hits.len());
                    for h in &hits {
                        println!("{}\t{}\t{}\t1.00", label.to_string(), h, kmers_query.len());
                    }
                }
                },
            }
        }
    }
}