debruijn/msp.rs
1// Copyright 2017 10x Genomics
2
3//! Methods for minimum substring partitioning of a DNA string
4//!
5//! simple_scan method is based on:
6//! Li, Yang. "MSPKmerCounter: a fast and memory efficient approach for k-mer counting." arXiv preprint arXiv:1505.06550 (2015).
7
8use crate::DnaSlice;
9use crate::Exts;
10use crate::Kmer;
11use crate::Vmer;
12use std::cmp::min;
13use std::cmp::Ordering;
14use std::iter::Iterator;
15use std::ops::Range;
16
17#[derive(Debug)]
18pub struct MspInterval {
19 bucket: u16,
20 start: u32,
21 len: u16,
22}
23
24impl MspInterval {
25 pub fn new(bucket: u16, start: u32, len: u16) -> MspInterval {
26 MspInterval { bucket, start, len }
27 }
28
29 pub fn start(&self) -> usize {
30 self.start as usize
31 }
32
33 pub fn len(&self) -> usize {
34 self.len as usize
35 }
36
37 pub fn is_empty(&self) -> bool {
38 self.len == 0
39 }
40
41 pub fn end(&self) -> usize {
42 (self.len as u32 + self.start) as usize
43 }
44
45 pub fn range(&self) -> Range<usize> {
46 self.start()..self.start() + self.len()
47 }
48
49 pub fn bucket(&self) -> u16 {
50 self.bucket
51 }
52}
53
54/// Determine MSP substrings of seq, for given k and p.
55/// Returns a vector of tuples indicating the substrings, and the pmer values:
56/// (p-mer value, min p-mer position, start position, end position)
57/// permutation is a permutation of the lexicographically-sorted set of all pmers.
58/// A permutation of pmers sorted by their inverse frequency in the dataset will give the
59/// most even bucketing of MSPs over pmers.
60#[deprecated(note = "Please use the `Scanner` type instead ")]
61pub fn simple_scan<V: Vmer, P: Kmer>(
62 k: usize,
63 seq: &V,
64 permutation: &[usize],
65 rc: bool,
66) -> Vec<MspInterval> {
67 // Can't partition strings shorter than k
68 assert!(seq.len() >= k);
69 assert!(P::k() <= 8);
70 assert!(seq.len() < 1 << 32);
71
72 let score = |pi: &P| {
73 if rc {
74 min(
75 permutation[pi.to_u64() as usize],
76 permutation[pi.rc().to_u64() as usize],
77 )
78 } else {
79 permutation[pi.to_u64() as usize]
80 }
81 };
82
83 let scanner = Scanner::new(seq, score, k);
84 let res = scanner.scan();
85
86 res.into_iter()
87 .map(|slc| MspInterval {
88 bucket: slc.bucket() as u16,
89 start: slc.start,
90 len: slc.len,
91 })
92 .collect()
93}
94
95/// Represents a sequence interval composed of
96/// successive k-mers that share a
97/// common minizer p-mer.
98#[derive(Debug)]
99pub struct MspIntervalP<P> {
100 /// The minimizing p-mer in this interval
101 pub minimizer: P,
102 /// The start of the sequence interval
103 pub start: u32,
104 /// The length of the sequence interval
105 pub len: u16,
106 /// The position of the minimizer
107 pub minimizer_pos: u32,
108}
109
110impl<P: Kmer> MspIntervalP<P> {
111 /// The shard 'bucket' identifer
112 /// to put this sequence in.
113 /// This is the reverse-complement canonicalized
114 /// value of the p-mer.
115 pub fn bucket(&self) -> u64 {
116 self.minimizer.min_rc().to_u64()
117 }
118}
119
120#[derive(PartialEq, Eq, Clone, Copy, Debug)]
121struct MinPos<P> {
122 val: usize,
123 pos: usize,
124 kmer: P,
125}
126
127impl<P: Eq> Ord for MinPos<P> {
128 fn cmp(&self, other: &Self) -> Ordering {
129 let val_cmp = self.val.cmp(&other.val);
130 if val_cmp != Ordering::Equal {
131 return val_cmp;
132 }
133
134 let pos_cmp = self.pos.cmp(&other.pos);
135 match pos_cmp {
136 Ordering::Equal => Ordering::Equal,
137 Ordering::Less => Ordering::Greater,
138 Ordering::Greater => Ordering::Less,
139 }
140 }
141}
142
143impl<P: Eq> PartialOrd for MinPos<P> {
144 fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
145 let val_cmp = self.val.cmp(&other.val);
146 if val_cmp != Ordering::Equal {
147 return Some(val_cmp);
148 }
149
150 let pos_cmp = self.pos.cmp(&other.pos);
151 Some(match pos_cmp {
152 Ordering::Equal => Ordering::Equal,
153 Ordering::Less => Ordering::Greater,
154 Ordering::Greater => Ordering::Less,
155 })
156 }
157}
158
159/// Determine MSP substrings of a sequence, for given k and p.
160/// The `scan()` method Returns a vector of tuples indicating the substrings,
161/// and the p-mer values as a set of `MspIntervalP<P>` values. A user-supplied
162/// score function is used to rank p-mers for the purposes of finding the minimizer.
163/// A permutation is a permutation of the lexicographically-sorted set of all pmers.
164/// A permutation of pmers sorted by their inverse frequency in the dataset will give the
165/// most even bucketing of MSPs over pmers.
166pub struct Scanner<'a, V, F, P> {
167 seq: &'a V,
168 score: F,
169 k: usize,
170 _mp: MinPos<P>,
171}
172
173impl<'a, V, F, P> Scanner<'a, V, F, P>
174where
175 V: Vmer,
176 P: Kmer,
177 F: Fn(&P) -> usize,
178{
179 /// Build a scanner for the `sequence`, using minimizer function `score_func`,
180 /// and kmer length `k`. The p-mer length is set by the type `P`.
181 pub fn new(sequence: &'a V, score_func: F, k: usize) -> Scanner<'a, V, F, P> {
182 Scanner {
183 seq: sequence,
184 score: score_func,
185 k,
186 _mp: MinPos {
187 val: 0,
188 pos: 0,
189 kmer: P::empty(),
190 },
191 }
192 }
193
194 fn mp(&self, pos: usize) -> MinPos<P> {
195 let kmer = self.seq.get_kmer::<P>(pos);
196 let val = (self.score)(&kmer);
197 MinPos { pos, val, kmer }
198 }
199
200 fn incr(&self, mp: &MinPos<P>) -> MinPos<P> {
201 let pos = mp.pos + 1;
202 let kmer = mp.kmer.extend_right(self.seq.get(pos + P::k() - 1));
203 let val = (self.score)(&kmer);
204 MinPos { pos, val, kmer }
205 }
206
207 pub fn scan(&self) -> Vec<MspIntervalP<P>> {
208 // Can't partition strings shorter than k
209 assert!(self.seq.len() >= self.k);
210 assert!(self.seq.len() < 1 << 32);
211
212 let seq = self.seq;
213 let m = seq.len();
214
215 let k = self.k;
216 let p = P::k();
217
218 let find_min = |start, stop| {
219 let mut min_pos = self.mp(start);
220 let mut current = min_pos;
221
222 while current.pos < stop {
223 current = self.incr(¤t);
224 min_pos = min(min_pos, current);
225 }
226
227 min_pos
228 };
229
230 let mut min_positions = Vec::with_capacity(16);
231
232 let mut min_pos = find_min(0, k - p);
233 let mut end_pos = self.mp(k - p);
234
235 min_positions.push((0, min_pos));
236
237 for i in 1..(m - k + 1) {
238 // end_pos always corresponds to i + k - p
239 end_pos = self.incr(&end_pos);
240
241 if i > min_pos.pos {
242 min_pos = find_min(i, i + k - p);
243 min_positions.push((i, min_pos));
244 } else if end_pos.val < min_pos.val {
245 min_pos = end_pos;
246 min_positions.push((i, min_pos));
247 }
248 }
249
250 let mut slices = Vec::with_capacity(min_positions.len());
251
252 // Generate the slices of the final string
253 for p in 0..min_positions.len() - 1 {
254 let (start_pos, min_pos) = min_positions[p];
255 let (next_pos, _) = min_positions[p + 1];
256
257 let interval = MspIntervalP {
258 minimizer: min_pos.kmer,
259 minimizer_pos: min_pos.pos as u32,
260 start: start_pos as u32,
261 len: (next_pos + k - 1 - start_pos) as u16,
262 };
263 slices.push(interval);
264 }
265
266 let (last_pos, min_pos) = min_positions[min_positions.len() - 1];
267 let last_interval = MspIntervalP {
268 minimizer: min_pos.kmer,
269 minimizer_pos: min_pos.pos as u32,
270 start: last_pos as u32,
271 len: (m - last_pos) as u16,
272 };
273 slices.push(last_interval);
274
275 slices
276 }
277}
278
279pub fn msp_sequence<P, V>(
280 k: usize,
281 seq: &[u8],
282 permutation: Option<&[usize]>,
283 rc: bool,
284) -> Vec<(u32, Exts, V)>
285where
286 P: Kmer,
287 V: Vmer,
288{
289 let p = P::k();
290
291 // Make sure the substrings will fit into the Vmers
292 assert!(V::max_len() >= 2 * k - p);
293
294 if seq.len() < k {
295 return Vec::new();
296 }
297
298 let default_perm: Option<Vec<usize>> = match permutation {
299 Some(_) => None,
300 None => Some((0..1 << (2 * p)).collect()),
301 };
302
303 let perm = permutation.unwrap_or_else(|| default_perm.as_ref().unwrap());
304
305 let score = |pi: &P| {
306 if rc {
307 min(perm[pi.to_u64() as usize], perm[pi.rc().to_u64() as usize])
308 } else {
309 perm[pi.to_u64() as usize]
310 }
311 };
312
313 let dna = DnaSlice(seq);
314 let msp_parts = Scanner::new(&dna, score, k).scan();
315 msp_parts
316 .into_iter()
317 .map(|msp| {
318 let v =
319 V::from_slice(&seq[(msp.start as usize)..(msp.start as usize + msp.len as usize)]);
320 let exts = Exts::from_slice_bounds(seq, msp.start as usize, msp.len as usize);
321 (msp.bucket() as u32, exts, v)
322 })
323 .collect()
324}
325
326#[cfg(test)]
327mod tests {
328 use super::*;
329 //use crate::dna_string::DnaString;
330 //use crate::kmer::{Kmer10, Kmer12, Kmer14, Kmer15, Kmer1}
331 use crate::kmer::{Kmer5, Kmer8};
332 use crate::test;
333 use crate::DnaSlice;
334 use std::collections::HashSet;
335 use std::iter::FromIterator;
336
337 fn all_kmers<T>(k: usize, seq: &[T]) -> Vec<&[T]> {
338 (0..(seq.len() - k + 1)).map(|i| &seq[i..i + k]).collect()
339 }
340
341 fn test_all_kmers(k: usize, full_seq: &[u8], slices: Vec<MspInterval>) {
342 let start_kmers = HashSet::from_iter(all_kmers(k, full_seq));
343
344 let mut sliced_kmers = HashSet::new();
345
346 for msp in slices {
347 let slc = &full_seq[(msp.start as usize)..(msp.start as usize + msp.len as usize)];
348 sliced_kmers.extend(all_kmers(k, slc));
349 }
350
351 if start_kmers != sliced_kmers {
352 println!("start kmers: {:?}", start_kmers);
353 println!("sliced kmers: {:?}", sliced_kmers);
354 panic!("kmer sets not equal");
355 }
356 }
357
358 #[test]
359 fn test1() {
360 let v = [1u8, 2u8, 3u8, 4u8, 5u8, 6u8];
361 let s = &v[..];
362 let ak = all_kmers(2, s);
363
364 if ak[0] < ak[1] {
365 println!("sorts!")
366 }
367
368 println!("{:?}", ak);
369
370 let v = [6u8, 5u8, 4u8, 3u8, 2u8, 1u8, 0u8];
371 let s = &v[..];
372 let mut ak = all_kmers(2, s);
373
374 println!("{:?}", ak);
375 ak.sort();
376
377 if ak[0] < ak[1] {
378 println!("sorts!")
379 }
380 println!("{:?}", ak);
381 }
382
383 #[test]
384 fn test_slice() {
385 let p = 8;
386 let permutation: Vec<usize> = (0..(1 << (2 * p))).collect();
387
388 for _ in 0..100 {
389 let k = 50usize;
390 let dna = test::random_dna(150);
391 println!("{:?}", dna);
392 #[allow(deprecated)]
393 let slices = super::simple_scan::<_, Kmer8>(k, &DnaSlice(&dna), &permutation, true);
394 println!(
395 "Made {} slices from dna of length {:?}",
396 slices.len(),
397 dna.len()
398 );
399
400 println!("slices: {:?}", slices);
401 test_all_kmers(k, &dna[..], slices);
402 }
403 }
404/*
405 fn check_msp_slices<P, F>(
406 k: usize,
407 full_seq: &DnaString,
408 slices: Vec<MspIntervalP<P>>,
409 score: F,
410 ) where
411 P: Kmer,
412 F: Fn(&P) -> usize,
413 {
414 // Check all the correctness properties of the slices returned by MSP.
415 // 1. Each pmer is covered by one and only one slice
416 // 2. Slices are at least p and most 2k-p long
417 // 3. The selected pmer exists in the slice and is minimal given the scoring function
418 // 4. No slices can be extended to the right by covering an additional pmer without violating the above constraints
419
420 let p = P::k();
421
422 /*
423 For debugging:
424 println!("k: {}, p: {}", k, p);
425 println!("seq: {:?}", full_seq);
426 println!("slices: {:#?}", slices);
427 */
428
429 // 1. each kmer is covered exactly once.
430 let mut covered = vec![false; full_seq.len() - k + 1];
431
432 for s in &slices {
433 let start = s.start as usize;
434 let end = s.start as usize + s.len as usize - k + 1;
435
436 for (i, c) in covered.iter_mut().take(end).skip(start).enumerate() {
437 assert!(!*c, "at {}\nbase already covered!", i);
438
439 *c = true;
440 }
441 }
442 assert!(covered.iter().all(|x| *x), "a pmer wasn't covered");
443
444 //2. p <= slice.len <= 2k-p
445 for s in &slices {
446 assert!(s.len as usize >= p);
447 assert!(s.len as usize <= 2 * k - p);
448 }
449
450 // 3. pmer exists and is the best possible within the slice
451 for s in &slices {
452 let slice_score = score(&s.minimizer);
453
454 let start = s.start as usize;
455 let end = s.start as usize + s.len as usize - p + 1;
456
457 // check the correct minimizer info is reported
458 assert_eq!(s.minimizer, full_seq.get_kmer(s.minimizer_pos as usize));
459
460 for i in start..end {
461 let pmer: P = full_seq.get_kmer(i);
462 let score = score(&pmer);
463
464 assert!(score >= slice_score, "found better pmer within slice")
465 }
466 }
467
468 // 4. No slices can be extended to the right by covering an additional pmer without violating the above constraints
469 for s in &slices[0..(slices.len() - 1)] {
470 // the next kmer beyond the end of the slice must either:
471 // a. not cover the minimizer
472 // b. contain a better minimizer
473 let next_kmer_pos = s.start as usize + s.len as usize - k + 1;
474 let next_kmer_covers_pmer = next_kmer_pos <= s.minimizer_pos as usize;
475
476 let next_pmer_pos = s.start as usize + s.len as usize - p + 1;
477 let next_pmer: P = full_seq.get_kmer(next_pmer_pos);
478 let next_score = score(&next_pmer);
479
480 let correct_end = next_score < score(&s.minimizer) || !next_kmer_covers_pmer;
481
482 assert!(
483 correct_end,
484 "detected a sequence slice that ended before it should have."
485 )
486 }
487 }
488
489 fn test_new_slicer<P: Kmer>(k: usize) {
490 let p = P::k();
491 if p >= k {
492 return;
493 }
494
495 let score = |p: &P| p.to_u64() as usize;
496
497 for i in 0..(20 * k) {
498 let len = 2 * i;
499 if len < k {
500 continue;
501 }
502
503 let dna = test::random_dna(len);
504 let dna = DnaString::from_bytes(&dna);
505
506 // use lexicographic ordering
507 let scanner = Scanner::new(&dna, score, k);
508 let slices = scanner.scan();
509 check_msp_slices(k, &dna, slices, score);
510
511 // use AT count ordering
512 let at_score = |p: &P| p.at_count() as usize;
513 let scanner = Scanner::new(&dna, at_score, k);
514 let slices = scanner.scan();
515 check_msp_slices(k, &dna, slices, at_score);
516 }
517
518 for i in 0..(4 * k) {
519 let len = i;
520 if len < k {
521 continue;
522 }
523
524 let dna = DnaString::blank(i);
525
526 let scanner = Scanner::new(&dna, score, k);
527 let slices = scanner.scan();
528 check_msp_slices(k, &dna, slices, score);
529 }
530 }
531
532 // TODO check why this test takes so long
533 #[test]
534 fn test_msp_scanner() {
535 for k in 16..64 {
536 test_new_slicer::<Kmer5>(k);
537 test_new_slicer::<Kmer8>(k);
538 test_new_slicer::<Kmer10>(k);
539 test_new_slicer::<Kmer12>(k);
540 test_new_slicer::<Kmer14>(k);
541 test_new_slicer::<Kmer15>(k);
542 test_new_slicer::<Kmer16>(k);
543 }
544 } */
545
546 #[test]
547 fn test_sample() {
548 // for testing MSP on specific error cases
549
550 // let v1 : Vec<u8> = vec![3, 0, 3, 0, 3, 2, 3, 1, 0, 0, 2, 0, 0, 1, 1, 3, 0, 2, 0, 3, 2, 3, 3, 1, 2, 0, 2, 2, 3, 1, 0, 1, 2, 1, 2, 3, 2, 1, 0, 2, 1, 0, 2, 2, 1, 0, 2, 3, 0, 3, 2, 3, 0, 2, 0, 0, 1, 0, 2, 3, 1, 3, 2, 0, 2, 2, 2, 2, 1, 2, 0, 2, 1, 0, 0, 1, 3, 0, 0, 0, 2, 2, 3, 0, 0, 3, 2, 3, 1, 2, 0, 2, 0, 3, 2, 1, 3, 2, 2, 3, 3, 1, 2, 3, 0, 3, 1, 1, 2, 2, 2, 3, 1, 1, 2, 0, 2, 3, 3, 2, 0, 3, 1, 1, 1, 1, 3, 3, 3, 0, 1, 0, 0, 3, 0, 2, 3, 2, 2, 2, 3, 2, 3, 1, 0, 2, 3, 0, 2, 2, 2, 1, 0, 2, 3, 3, 3, 2, 2, 1, 3, 1, 0, 2, 0, 1, 0, 1, 2, 2, 2, 2, 2, 3, 0, 3, 3, 1, 2, 3, 2, 0, 2, 1, 0, 3, 1, 1, 0, 2, 1, 3, 0, 1, 2, 1, 1, 1, 3, 1, 0, 3, 0, 0, 2, 2, 0, 3, 2, 3, 2, 0, 2, 1, 2, 3, 0, 1, 3, 1, 1, 2, 2, 3, 0, 3, 3, 3, 3, 0, 3, 3, 0, 0, 0, 2, 3, 3, 1, 3, 3, 1, 2, 2, 0, 3, 2, 1, 0, 2];
551 // let v2 : Vec<u8> = vec![3, 0, 3, 0, 3, 2, 3, 1, 0, 0, 2, 0, 0, 1, 1, 3, 0, 2, 0, 3, 2, 3, 3, 1, 2, 0, 2, 2, 3, 2, 0, 3, 2, 1, 2, 3, 2, 1, 0, 2, 1, 0, 2, 2, 1, 0, 2, 3, 0, 3, 2, 3, 0, 2, 3, 0, 1, 0, 2, 3, 1, 3, 2, 0, 2, 2, 0, 2, 1, 2, 0, 2, 1, 0, 0, 1, 3, 0, 0, 0, 2, 2, 3, 0, 0, 3, 2, 3, 1, 2, 0, 2, 0, 3, 2, 1, 3, 2, 2, 3, 3, 1, 2, 3, 0, 3, 1, 1, 2, 2, 2, 3, 1, 1, 2, 0, 0, 3, 3, 2, 0, 3, 1, 1, 1, 1, 3, 3, 3, 0, 1, 0, 0, 3, 0, 2, 3, 1, 2, 1, 3, 2, 3, 1, 0, 2, 3, 0, 2, 2, 2, 2, 1, 2, 3, 3, 3, 2, 2, 1, 3, 1, 0, 2, 0, 1, 0, 1, 2, 2, 2, 2, 2, 3, 0, 3, 3, 1, 2, 3, 2, 0, 2, 1, 0, 3, 1, 1, 0, 2, 1, 3, 0, 1, 2, 1, 2, 1, 3, 1, 0, 3, 0, 0, 2, 2, 0, 3, 2, 3, 2, 0, 2, 1, 2, 3, 0, 1, 3, 1, 1, 2, 2, 3, 0, 3, 3, 3, 3, 0, 3, 3, 0, 0, 0, 2, 3, 3, 1, 0, 3, 1, 2, 2, 0, 3, 2, 1, 0, 2];
552
553 let v1: Vec<u8> = vec![
554 3, 0, 3, 0, 1, 2, 3, 3, 0, 0, 0, 1, 1, 0, 3, 3, 1, 2, 0, 1, 1, 3, 2, 1, 1, 1, 2, 3, 3,
555 2, 1, 2, 2, 1, 2, 3, 2, 1, 0, 2, 1, 1, 2, 1, 0, 1, 2, 3, 0, 3, 2, 3, 0, 1, 3, 0, 1, 0,
556 2, 3, 3, 3, 2, 3, 2, 2, 0, 2, 1, 2, 2, 2, 0, 0, 1, 0, 1, 2, 1, 0, 3, 2, 3, 1, 0, 3, 2,
557 2, 2, 2, 1, 3, 0, 2, 1, 2, 1, 3, 3, 0, 1, 1, 2, 3, 0, 2, 3, 1, 3, 2, 3, 1, 1, 0, 2, 2,
558 1, 1, 1, 2, 0, 0, 1, 2, 1, 0, 3, 3, 3, 0, 1, 1, 0, 3, 0, 2, 2, 1, 2, 1, 3, 2, 3, 1, 0,
559 2, 3, 0, 2, 2, 0, 3, 3, 2, 3, 3, 0, 3, 0, 1, 1, 3, 0, 1, 0, 1, 0, 1, 2, 2, 0, 3, 3, 3,
560 1, 3, 1, 1, 0, 1, 2, 0, 2, 1, 0, 3, 1, 1, 1, 2, 0, 3, 0, 1, 2, 1, 1, 1, 3, 1, 2, 3, 0,
561 3, 2, 2, 0, 3, 2, 3, 3, 0, 2, 3, 2, 0, 0, 1, 2, 0, 3, 2, 2, 1, 2, 2, 3, 3, 3, 2, 3, 0,
562 3, 1, 0, 2, 3, 3, 0, 3, 1, 0, 3, 2, 1, 3, 2, 1, 1, 2,
563 ];
564 let v2: Vec<u8> = vec![
565 3, 0, 3, 0, 1, 3, 3, 3, 0, 0, 0, 1, 1, 0, 3, 3, 1, 2, 0, 1, 1, 3, 2, 1, 1, 1, 2, 1, 3,
566 2, 1, 2, 2, 1, 2, 3, 2, 1, 0, 2, 1, 1, 2, 1, 1, 1, 2, 3, 0, 3, 2, 3, 0, 1, 3, 0, 1, 0,
567 2, 3, 3, 3, 2, 3, 2, 2, 0, 2, 1, 2, 2, 2, 0, 0, 1, 0, 1, 2, 1, 0, 3, 2, 3, 1, 0, 3, 2,
568 2, 2, 2, 1, 3, 0, 2, 1, 2, 1, 3, 3, 0, 1, 1, 2, 3, 0, 2, 3, 1, 3, 2, 3, 1, 1, 0, 2, 2,
569 1, 1, 0, 2, 0, 0, 1, 2, 1, 0, 3, 3, 3, 0, 1, 1, 0, 3, 0, 2, 2, 1, 2, 1, 3, 2, 3, 1, 0,
570 2, 3, 0, 2, 2, 0, 3, 3, 2, 3, 3, 0, 3, 0, 1, 1, 3, 0, 1, 0, 1, 0, 1, 2, 2, 0, 3, 3, 3,
571 1, 3, 1, 1, 0, 1, 2, 0, 2, 1, 0, 3, 1, 1, 1, 2, 0, 3, 0, 1, 2, 1, 3, 1, 3, 1, 2, 3, 0,
572 3, 2, 2, 0, 3, 2, 3, 3, 0, 2, 3, 2, 0, 0, 1, 2, 0, 3, 2, 2, 1, 2, 2, 3, 3, 3, 2, 3, 1,
573 3, 1, 0, 2, 2, 3, 0, 3, 1, 0, 3, 3, 1, 3, 2, 1, 1, 2,
574 ];
575
576 let p = 5;
577 let permutation: Vec<usize> = (0..(1 << (2 * p))).collect();
578
579 #[allow(deprecated)]
580 let s1 = super::simple_scan::<_, Kmer5>(35, &DnaSlice(&v1), &permutation, true);
581 #[allow(deprecated)]
582 let s2 = super::simple_scan::<_, Kmer5>(35, &DnaSlice(&v2), &permutation, true);
583
584 println!("{:?}", s1);
585 println!("{:?}", s2);
586 }
587}