bsalign 0.7.0

Rust bindings for the bsalign C library
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
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
use std::{
    ops::{Deref, DerefMut},
    os::raw::c_char,
};

use crate::{AlignMode, AlignScore};
use bsalign_sys::bindings;

#[derive(Debug)]
pub struct BsPoaParam(bindings::BSPOAPar);

impl Deref for BsPoaParam {
    type Target = bindings::BSPOAPar;
    fn deref(&self) -> &Self::Target {
        &self.0
    }
}

impl DerefMut for BsPoaParam {
    fn deref_mut(&mut self) -> &mut Self::Target {
        &mut self.0
    }
}

#[derive(Debug)]
pub struct BsPoaAligner {
    pub params: BsPoaParam,
    pub poa: *mut bindings::BSPOA,
    metainfo: Option<String>,
    aligned: bool,
    pub nseq: usize,
}

impl Drop for BsPoaAligner {
    fn drop(&mut self) {
        unsafe {
            bindings::bspoa_free(self.poa);
        }
    }
}

impl BsPoaAligner {
    pub fn new(params: BsPoaParam) -> Self {
        let poa = unsafe {
            let poa = bindings::bspoa_init(params.0.clone());
            bindings::bspoa_begin(poa);
            poa
        };
        BsPoaAligner {
            params,
            poa,
            metainfo: None,
            aligned: false,
            nseq: 0,
        }
    }
}

#[inline]
fn base_from_2bit(base: u8) -> u8 {
    match base {
        0 => b'A',
        1 => b'C',
        2 => b'G',
        3 => b'T',
        _ => b'N',
    }
}

#[derive(Debug)]
pub struct Consensus<'a> {
    inner: &'a [u8],
    convert: bool,
}

impl<'a> Consensus<'a> {
    pub fn as_string(&self) -> String {
        if self.convert {
            let mut s = String::with_capacity(self.len());
            for i in 0..self.len() {
                let base = self.get_base(i).unwrap();
                s.push(base as char);
            }
            return s;
        }
        unsafe {
            let s = str::from_utf8_unchecked(self.inner);
            s.into()
        }
    }

    pub fn get_bit(&self, index: usize) -> Option<u8> {
        if index >= self.len() {
            return None;
        }
        let bit = self.inner[index];
        Some(bit)
    }

    pub fn get_base(&self, index: usize) -> Option<u8> {
        if index >= self.len() {
            return None;
        }
        let bit = self.inner[index];
        Some(base_from_2bit(bit))
    }

    pub fn len(&self) -> usize {
        self.inner.len()
    }
}

#[derive(Debug)]
pub struct Quality<'a> {
    inner: &'a [u8],
}

impl<'a> Quality<'a> {
    pub fn len(&self) -> usize {
        self.inner.len()
    }

    pub fn as_string(&self) -> String {
        let mut s = String::with_capacity(self.len());
        for i in 0..self.len() {
            let base = self.get(i).unwrap();
            s.push((b'!' + base) as char);
        }
        s
    }

    pub fn get(&self, index: usize) -> Option<u8> {
        if index >= self.len() {
            return None;
        }
        Some(self.inner[index])
    }
}

#[derive(Debug)]
pub struct Alt<'a> {
    inner: &'a [u8],
}

impl<'a> Alt<'a> {
    pub fn len(&self) -> usize {
        self.inner.len()
    }

    pub fn as_string(&self) -> String {
        let mut s = String::with_capacity(self.len());
        for i in 0..self.len() {
            let base = self.get(i).unwrap();
            s.push((b'!' + base) as char);
        }
        s
    }

    pub fn get(&self, index: usize) -> Option<u8> {
        if index >= self.len() {
            return None;
        }
        Some(self.inner[index])
    }
}

pub struct PoaAlignResult<'a> {
    idx: usize,
    poa: &'a BsPoaAligner,
}

impl<'a> Iterator for PoaAlignResult<'a> {
    type Item = AlignmentString<'a>;

    fn next(&mut self) -> Option<Self::Item> {
        if self.idx >= self.poa.nseq {
            return None;
        }
        let a = self.poa.get_alignment(self.idx);
        self.idx += 1;
        a
    }
}

#[derive(Debug)]
pub struct AlignmentString<'a> {
    inner: &'a [u8],
}

impl<'a> AlignmentString<'a> {
    pub fn len(&self) -> usize {
        self.inner.len()
    }

    pub fn as_string(&self) -> String {
        unsafe {
            let s = str::from_utf8_unchecked(self.inner);
            s.into()
        }
    }

    pub fn as_bytes(&self) -> &'a [u8] {
        self.inner
    }
}

impl BsPoaAligner {
    pub fn add_sequence<T>(&mut self, seq: &T)
    where
        T: AsRef<[u8]> + ?Sized,
    {
        let s = seq.as_ref();
        unsafe {
            bindings::bspoa_add_sequence(self.poa, s.as_ptr(), s.len() as u32);
            self.aligned = false;
        }
        self.nseq += 1;
    }

    /// Reset poa so that start next loop align.
    /// !!Note: Call this method before `add_sequence`
    pub fn reset(&mut self) {
        unsafe {
            self.aligned = false;
            self.metainfo = None;
            self.nseq = 0;
            bindings::bspoa_begin(self.poa);
        }
    }

    pub fn align(&mut self) {
        unsafe {
            bindings::bspoa_end(self.poa);
            bindings::bspoa_tidy_msa(self.poa);
            self.aligned = true;
        }
    }

    pub fn get_alignment_result(&self) -> PoaAlignResult<'_> {
        if !self.aligned {
            panic!("Align sequences before calling dump_msa, call `align` first");
        }
        return PoaAlignResult { idx: 0, poa: self };
    }

    pub fn call_snvs(&mut self) {
        if !self.aligned {
            panic!("Align sequences before calling call_snvs, call `align` first");
        }
        unsafe {
            bindings::bspoa_call_snvs(self.poa);
        }
    }

    pub fn metainfo(&self) -> Option<&String> {
        self.metainfo.as_ref()
    }

    pub fn set_metainfo(&mut self, metainfo: String) {
        self.metainfo = Some(metainfo);
    }

    pub fn dump_msa(&self, filename: &str) {
        if !self.aligned {
            panic!("Align sequences before calling dump_msa, call `align` first");
        }
        unsafe {
            let cfilename = std::ffi::CString::new(filename).unwrap();
            if let Some(metainfo) = self.metainfo() {
                bindings::bspoa_dump_binary_msa(
                    self.poa,
                    metainfo.as_ptr() as *mut i8,
                    metainfo.len() as u32,
                    cfilename.as_ptr() as *mut i8,
                );
            } else {
                bindings::bspoa_dump_binary_msa(
                    self.poa,
                    std::ptr::null_mut(),
                    0,
                    cfilename.as_ptr() as *mut i8,
                );
            }
        }
    }

    pub fn get_cns_with_indel(&self) -> Consensus<'_> {
        if !self.aligned {
            panic!("POA must be aligned before getting consensus");
        }
        let alignment = self.get_alignment(self.nseq);
        Consensus {
            inner: alignment.unwrap().inner,
            convert: false,
        }
    }

    pub fn get_cns(&self) -> Consensus<'_> {
        if !self.aligned {
            panic!("Align sequences before calling get_consensus, call `align` first");
        }
        let c = unsafe {
            let mut len = 0;
            let cns = bindings::bspoa_get_cns(self.poa, &mut len);
            std::slice::from_raw_parts(cns, len as usize)
        };
        Consensus {
            inner: c,
            convert: true,
        }
    }

    pub fn get_qlt(&self) -> Quality<'_> {
        if !self.aligned {
            panic!("Align sequences before calling get_qlt, call `align` first");
        }
        let q = unsafe {
            let mut len = 0;
            let qlt = bindings::bspoa_get_qlt(self.poa, &mut len);
            std::slice::from_raw_parts(qlt, len as usize)
        };
        Quality { inner: q }
    }

    pub fn get_alignment(&self, idx: usize) -> Option<AlignmentString<'_>> {
        if !self.aligned {
            panic!("Align sequences before calling get_qlt, call `align` first");
        }
        if idx > self.nseq {
            return None;
        }
        let idx = idx + 1;
        let mut len: usize = 0;
        let inner = unsafe {
            let buffer = bindings::bspoa_get_rid_alignment(self.poa, idx as i32, &mut len);
            std::slice::from_raw_parts(buffer, len)
        };
        Some(AlignmentString { inner })
    }

    pub fn get_alt(&self) -> Alt<'_> {
        if !self.aligned {
            panic!("Align sequences before calling get_alt, call `align` first");
        }
        let a = unsafe {
            let mut len = 0;
            let alt = bindings::bspoa_get_alt(self.poa, &mut len);
            std::slice::from_raw_parts(alt, len as usize)
        };
        Alt { inner: a }
    }
    /// Print SNVs to a file
    pub fn print_snvs(&self, label: Option<&str>, filename: &str) {
        if !self.aligned {
            panic!("Align sequences before calling print_snvs, call `align` first");
        }
        unsafe {
            let cfilename = std::ffi::CString::new(filename).unwrap();
            let clabel = match label {
                None => std::ffi::CString::new("SNV").unwrap(),
                Some(s) => std::ffi::CString::new(s).unwrap(),
            };

            let cfile = bindings::c_file_open(cfilename.as_ptr() as *mut c_char);

            bindings::bspoa_print_snvs(self.poa, clabel.as_ptr() as *mut i8, cfile);
        }
    }
}

impl BsPoaAligner {
    pub fn load_msa(filename: &str) -> Option<BsPoaAligner> {
        let p = BsPoaParam::default();
        let mut poa = BsPoaAligner::new(p);
        let cfilename = std::ffi::CString::new(filename).unwrap();
        unsafe {
            let s = bindings::string_init(32);
            let ret = bindings::bspoa_load_binary_msa(poa.poa, cfilename.as_ptr() as *mut i8, s);
            if ret != 0 {
                bindings::string_free(s);
                return None;
            }
            let olds = s;
            let s = s.as_ref().unwrap();
            if s.size > 0 {
                let slice = std::slice::from_raw_parts(s.string as *const u8, s.size as usize);
                let metainfo = String::from_utf8_lossy(slice).to_string();
                poa.set_metainfo(metainfo);
            }
            bindings::string_free(olds);
            Some(poa)
        }
    }
}

impl Default for BsPoaParam {
    fn default() -> Self {
        let p = bindings::BSPOAPar {
            refmode: 0,
            shuffle: 1,
            alnmode: bindings::SEQALIGN_MODE_OVERLAP as i32,
            realn: 3,
            seqcore: 40,
            nrec: 20,
            ksz: 15,
            bwtrigger: 1,
            bandwidth: 128,
            M: 2,
            X: -6,
            O: -3,
            E: -2,
            Q: -8,
            P: -1,
            T: 20,
            refbonus: 1,
            editbw: 64,
            althi: 5,
            qlthi: 70,
            psub: 0.10,
            pins: 0.10,
            pdel: 0.15,
            piex: 0.15,
            pdex: 0.20,
            hins: 0.20,
            hdel: 0.40,
            min_varcnt: 3,
            min_covfrq: 0.5,
            min_snvqlt: 5,
        };
        BsPoaParam(p)
    }
}

#[derive(Debug)]
pub enum PoaParamRefMode {
    NoReference = 0,
    Reference = 1,
}

impl Default for PoaParamRefMode {
    fn default() -> Self {
        PoaParamRefMode::NoReference
    }
}

#[derive(Debug)]
pub struct PoaParamCns {
    pub psub: f32,
    pub pins: f32,
    pub pdel: f32,
    pub piex: f32,
    pub pdex: f32,
    pub hins: f32,
    pub hdel: f32,
}

impl Default for PoaParamCns {
    fn default() -> Self {
        PoaParamCns {
            psub: 0.10,
            pins: 0.10,
            pdel: 0.15,
            piex: 0.15,
            pdex: 0.20,
            hins: 0.20,
            hdel: 0.40,
        }
    }
}

impl BsPoaParam {
    /// Set reference mode
    /// 1: The first sequence is the reference
    /// 0: No reference
    pub fn set_refmode(self, refmode: PoaParamRefMode) -> Self {
        let mut p = self;
        p.refmode = refmode as i32;
        p
    }

    /// whether to shuffle/sort the reads according to most kmers matches first
    pub fn set_shuffle(self, enable: bool) -> Self {
        let mut p = self;
        if enable {
            p.shuffle = 1;
        } else {
            p.shuffle = 0;
        }
        p
    }

    /// Set alignment mode
    ///
    /// global/extend/overlap
    pub fn set_alnmode(self, alnmode: AlignMode) -> Self {
        let mut p = self;
        p.alnmode = alnmode as i32;
        p
    }

    /// Set rounds of realignment
    pub fn set_realn(self, realn: u32) -> Self {
        let mut p = self;
        p.realn = realn as i32;
        p
    }

    /// number of seqs in core MSA, addtional reads will be realigned (realn > 0) against core MSA profile to build a full MSA
    pub fn set_seqcore(self, seqcore: u32) -> Self {
        let mut p = self;
        p.seqcore = seqcore;
        p
    }

    /// every query read is aligning against previous <nrec> reads on graph, 0 to all the previous
    pub fn set_nrec(self, nrec: u32) -> Self {
        let mut p = self;
        p.nrec = nrec as i32;
        p
    }

    /// kmer size for graph construction
    pub fn set_ksz(self, ksz: u32) -> Self {
        let mut p = self;
        p.ksz = ksz as i32;
        p
    }

    /// when <trigger> > 0 and <-W> < query length, genrates CNS per after <trigger> reads, and trigger banded alignment
    pub fn set_bwtrigger(self, bwtrigger: i32) -> Self {
        let mut p = self;
        p.bwtrigger = bwtrigger;
        p
    }

    /// bandwidth for realignment
    pub fn set_bandwidth(self, bandwidth: u32) -> Self {
        let mut p = self;
        p.bandwidth = bandwidth as i32;
        p
    }

    /// Set score parameters
    pub fn set_score(self, score: AlignScore) -> Self {
        let mut p = self;
        p.M = score.M;
        p.X = score.X;
        p.O = score.O;
        p.E = score.E;
        p.Q = score.Q;
        p.P = score.P;
        p
    }

    /// base match score on reference will be M + refbonus
    pub fn set_refbonus(self, refbonus: u32) -> Self {
        let mut p = self;
        p.refbonus = refbonus as i32;
        p
    }

    /// Set edit bandwidth
    pub fn set_editbw(self, editbw: u32) -> Self {
        let mut p = self;
        p.editbw = editbw as i32;
        p
    }

    /// cutoff of high alt base quality, used in tidying out possible SNV sites
    pub fn set_althi(self, althi: u32) -> Self {
        let mut p = self;
        p.althi = althi as i32;
        p
    }

    /// cutoff of high cns base quality, used in print colorful MSA
    pub fn set_qlthi(self, qlthi: u32) -> Self {
        let mut p = self;
        p.qlthi = qlthi as i32;
        p
    }

    /// min count of variant base in SNV sites
    pub fn set_varcnt(self, min_varcnt: u32) -> Self {
        let mut p = self;
        p.min_varcnt = min_varcnt as i32;
        p
    }

    /// min spanned_reads / total_reads in SNV sites
    pub fn set_covfrq(self, min_covfrq: f32) -> Self {
        let mut p = self;
        p.min_covfrq = min_covfrq;
        p
    }

    /// -log10(p-value), 5 ~= 1e-5
    pub fn set_snvqlt(self, min_snvqlt: u32) -> Self {
        let mut p = self;
        p.min_snvqlt = min_snvqlt;
        p
    }

    pub fn set_consensus_param(self, cns: PoaParamCns) -> Self {
        let mut p = self;
        p.psub = cns.psub;
        p.pins = cns.pins;
        p.pdel = cns.pdel;
        p.piex = cns.piex;
        p.pdex = cns.pdex;
        p.hins = cns.hins;
        p.hdel = cns.hdel;
        p
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    #[test]
    fn test_poa() {
        let bspoa_params = BsPoaParam::default().set_varcnt(1);
        let seq1 = "TCATCTGATGCTAGTAcccccccc";
        let seq2 = "TCCTGATCTAGCTAGTA";
        let seq3 = "TCGATCTGATAGCTAGTA";
        let mut poa = BsPoaAligner::new(bspoa_params);
        poa.add_sequence(seq1);
        poa.add_sequence(seq2);
        poa.add_sequence(seq3);
        poa.align();
        let cns = poa.get_cns();
        println!("Consensus: {}", cns.as_string());
        let cns_with_indel = poa.get_cns_with_indel();
        println!("Consensus with indel: {}", cns_with_indel.as_string());
        assert_eq!(cns.as_string(), "TCATCTGATTAGCTAGTACCCCCCCC");
        let qlt = poa.get_qlt();
        assert_eq!(qlt.len(), cns.len());
        poa.set_metainfo("test".to_string());
        poa.dump_msa("test.msa");
        let poa = BsPoaAligner::load_msa("test.msa");
        assert!(poa.is_some());
        if let Some(poa) = poa {
            assert_eq!(poa.metainfo().unwrap(), "test");
        }
        std::fs::remove_file("test.msa").unwrap();
    }
}