• Home
  • Features
  • Pricing
  • Docs
  • Announcements
  • Sign In

rust-bio / rust-htslib / 5327586180

pending completion
5327586180

Pull #398

github

web-flow
Merge 5d7d33350 into 261e3a3b0
Pull Request #398: fix: use correct return value in bcf_get_format and bcf_get_info_values

37 of 37 new or added lines in 1 file covered. (100.0%)

2439 of 3058 relevant lines covered (79.76%)

2.28 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

76.31
/src/bam/record.rs
1
// Copyright 2014 Christopher Schröder, Johannes Köster.
2
// Licensed under the MIT license (http://opensource.org/licenses/MIT)
3
// This file may not be copied, modified, or distributed
4
// except according to those terms.
12✔
5

6
use std::convert::TryFrom;
7
use std::ffi;
8
use std::fmt;
9
use std::marker::PhantomData;
10
use std::mem::{size_of, MaybeUninit};
11
use std::ops;
12
use std::os::raw::c_char;
13
use std::rc::Rc;
14
use std::slice;
15
use std::str;
16
use std::u32;
17

18
use byteorder::{LittleEndian, ReadBytesExt};
19

20
use crate::bam::Error;
21
use crate::bam::HeaderView;
22
use crate::errors::Result;
23
use crate::htslib;
24
use crate::utils;
25
#[cfg(feature = "serde_feature")]
26
use serde::{self, Deserialize, Serialize};
27

28
use bio_types::alignment::{Alignment, AlignmentMode, AlignmentOperation};
29
use bio_types::genome;
30
use bio_types::sequence::SequenceRead;
31
use bio_types::sequence::SequenceReadPairOrientation;
32
use bio_types::strand::ReqStrand;
33

34
/// A macro creating methods for flag access.
35
macro_rules! flag {
36
    ($get:ident, $set:ident, $unset:ident, $bit:expr) => {
37
        pub fn $get(&self) -> bool {
12✔
38
            self.inner().core.flag & $bit != 0
11✔
39
        }
40

41
        pub fn $set(&mut self) {
4✔
42
            self.inner_mut().core.flag |= $bit;
3✔
43
        }
44

45
        pub fn $unset(&mut self) {
4✔
46
            self.inner_mut().core.flag &= !$bit;
3✔
47
        }
48
    };
49
}
50

51
/// A BAM record.
52
pub struct Record {
53
    pub inner: htslib::bam1_t,
54
    own: bool,
55
    cigar: Option<CigarStringView>,
56
    header: Option<Rc<HeaderView>>,
57
}
58

59
unsafe impl Send for Record {}
60
unsafe impl Sync for Record {}
61

62
impl Clone for Record {
63
    fn clone(&self) -> Self {
1✔
64
        let mut copy = Record::new();
1✔
65
        unsafe { htslib::bam_copy1(copy.inner_ptr_mut(), self.inner_ptr()) };
2✔
66
        copy
1✔
67
    }
68
}
69

70
impl PartialEq for Record {
71
    fn eq(&self, other: &Record) -> bool {
1✔
72
        self.tid() == other.tid()
10✔
73
            && self.pos() == other.pos()
1✔
74
            && self.bin() == other.bin()
1✔
75
            && self.mapq() == other.mapq()
1✔
76
            && self.flags() == other.flags()
1✔
77
            && self.mtid() == other.mtid()
1✔
78
            && self.mpos() == other.mpos()
1✔
79
            && self.insert_size() == other.insert_size()
1✔
80
            && self.data() == other.data()
1✔
81
            && self.inner().core.l_extranul == other.inner().core.l_extranul
1✔
82
    }
83
}
84

85
impl Eq for Record {}
86

87
impl fmt::Debug for Record {
88
    fn fmt(&self, fmt: &mut fmt::Formatter<'_>) -> Result<(), fmt::Error> {
×
89
        fmt.write_fmt(format_args!(
×
90
            "Record(tid: {}, pos: {})",
91
            self.tid(),
×
92
            self.pos()
×
93
        ))
94
    }
95
}
96

97
impl Default for Record {
98
    fn default() -> Self {
×
99
        Self::new()
×
100
    }
101
}
102

103
#[inline]
104
fn extranul_from_qname(qname: &[u8]) -> usize {
1✔
105
    let qlen = qname.len() + 1;
1✔
106
    if qlen % 4 != 0 {
3✔
107
        4 - qlen % 4
2✔
108
    } else {
109
        0
1✔
110
    }
111
}
112

113
impl Record {
114
    /// Create an empty BAM record.
115
    pub fn new() -> Self {
9✔
116
        Record {
117
            inner: unsafe { MaybeUninit::zeroed().assume_init() },
9✔
118
            own: true,
119
            cigar: None,
120
            header: None,
121
        }
122
    }
123

124
    pub fn from_inner(from: *mut htslib::bam1_t) -> Self {
2✔
125
        Record {
126
            inner: {
127
                #[allow(clippy::uninit_assumed_init)]
128
                let mut inner = unsafe { MaybeUninit::uninit().assume_init() };
129
                unsafe {
130
                    ::libc::memcpy(
131
                        &mut inner as *mut htslib::bam1_t as *mut ::libc::c_void,
132
                        from as *const ::libc::c_void,
133
                        size_of::<htslib::bam1_t>(),
134
                    );
135
                }
136
                inner
137
            },
138
            own: false,
139
            cigar: None,
140
            header: None,
141
        }
142
    }
143

144
    // Create a BAM record from a line SAM text. SAM slice need not be 0-terminated.
145
    pub fn from_sam(header_view: &HeaderView, sam: &[u8]) -> Result<Record> {
3✔
146
        let mut record = Self::new();
3✔
147

148
        let mut sam_copy = Vec::with_capacity(sam.len() + 1);
3✔
149
        sam_copy.extend(sam);
3✔
150
        sam_copy.push(0);
3✔
151

152
        let mut sam_string = htslib::kstring_t {
153
            s: sam_copy.as_ptr() as *mut c_char,
3✔
154
            l: sam_copy.len() as u64,
3✔
155
            m: sam_copy.len() as u64,
3✔
156
        };
157

158
        let succ = unsafe {
159
            htslib::sam_parse1(
160
                &mut sam_string,
161
                header_view.inner_ptr() as *mut htslib::bam_hdr_t,
3✔
162
                record.inner_ptr_mut(),
3✔
163
            )
164
        };
165

166
        if succ == 0 {
6✔
167
            Ok(record)
3✔
168
        } else {
169
            Err(Error::BamParseSAM {
×
170
                rec: str::from_utf8(&sam_copy).unwrap().to_owned(),
×
171
            })
172
        }
173
    }
174

175
    pub fn set_header(&mut self, header: Rc<HeaderView>) {
6✔
176
        self.header = Some(header);
12✔
177
    }
178

179
    pub(super) fn data(&self) -> &[u8] {
4✔
180
        unsafe { slice::from_raw_parts(self.inner().data, self.inner().l_data as usize) }
4✔
181
    }
182

183
    #[inline]
184
    pub fn inner_mut(&mut self) -> &mut htslib::bam1_t {
1✔
185
        &mut self.inner
×
186
    }
187

188
    #[inline]
189
    pub(super) fn inner_ptr_mut(&mut self) -> *mut htslib::bam1_t {
7✔
190
        &mut self.inner as *mut htslib::bam1_t
×
191
    }
192

193
    #[inline]
194
    pub fn inner(&self) -> &htslib::bam1_t {
6✔
195
        &self.inner
×
196
    }
197

198
    #[inline]
199
    pub(super) fn inner_ptr(&self) -> *const htslib::bam1_t {
3✔
200
        &self.inner as *const htslib::bam1_t
×
201
    }
202

203
    /// Get target id.
204
    pub fn tid(&self) -> i32 {
2✔
205
        self.inner().core.tid
2✔
206
    }
207

208
    /// Set target id.
209
    pub fn set_tid(&mut self, tid: i32) {
×
210
        self.inner_mut().core.tid = tid;
×
211
    }
212

213
    /// Get position (0-based).
214
    pub fn pos(&self) -> i64 {
1✔
215
        self.inner().core.pos
1✔
216
    }
217

218
    /// Set position (0-based).
219
    pub fn set_pos(&mut self, pos: i64) {
1✔
220
        self.inner_mut().core.pos = pos;
1✔
221
    }
222

223
    pub fn bin(&self) -> u16 {
1✔
224
        self.inner().core.bin
1✔
225
    }
226

227
    pub fn set_bin(&mut self, bin: u16) {
×
228
        self.inner_mut().core.bin = bin;
×
229
    }
230

231
    /// Get MAPQ.
232
    pub fn mapq(&self) -> u8 {
1✔
233
        self.inner().core.qual
1✔
234
    }
235

236
    /// Set MAPQ.
237
    pub fn set_mapq(&mut self, mapq: u8) {
×
238
        self.inner_mut().core.qual = mapq;
×
239
    }
240

241
    /// Get strand information from record flags.
242
    pub fn strand(&self) -> ReqStrand {
×
243
        let reverse = self.flags() & 0x10 != 0;
×
244
        if reverse {
×
245
            ReqStrand::Reverse
×
246
        } else {
247
            ReqStrand::Forward
×
248
        }
249
    }
250

251
    /// Get raw flags.
252
    pub fn flags(&self) -> u16 {
2✔
253
        self.inner().core.flag
2✔
254
    }
255

256
    /// Set raw flags.
257
    pub fn set_flags(&mut self, flags: u16) {
×
258
        self.inner_mut().core.flag = flags;
×
259
    }
260

261
    /// Unset all flags.
262
    pub fn unset_flags(&mut self) {
×
263
        self.inner_mut().core.flag = 0;
×
264
    }
265

266
    /// Get target id of mate.
267
    pub fn mtid(&self) -> i32 {
1✔
268
        self.inner().core.mtid
1✔
269
    }
270

271
    /// Set target id of mate.
272
    pub fn set_mtid(&mut self, mtid: i32) {
×
273
        self.inner_mut().core.mtid = mtid;
×
274
    }
275

276
    /// Get mate position.
277
    pub fn mpos(&self) -> i64 {
1✔
278
        self.inner().core.mpos
1✔
279
    }
280

281
    /// Set mate position.
282
    pub fn set_mpos(&mut self, mpos: i64) {
×
283
        self.inner_mut().core.mpos = mpos;
×
284
    }
285

286
    /// Get insert size.
287
    pub fn insert_size(&self) -> i64 {
1✔
288
        self.inner().core.isize
1✔
289
    }
290

291
    /// Set insert size.
292
    pub fn set_insert_size(&mut self, insert_size: i64) {
×
293
        self.inner_mut().core.isize = insert_size;
×
294
    }
295

296
    fn qname_capacity(&self) -> usize {
4✔
297
        self.inner().core.l_qname as usize
4✔
298
    }
299

300
    fn qname_len(&self) -> usize {
2✔
301
        // discount all trailing zeros (the default one and extra nulls)
302
        self.qname_capacity() - 1 - self.inner().core.l_extranul as usize
2✔
303
    }
304

305
    /// Get qname (read name). Complexity: O(1).
306
    pub fn qname(&self) -> &[u8] {
2✔
307
        &self.data()[..self.qname_len()]
2✔
308
    }
309

310
    /// Set the variable length data buffer
311
    pub fn set_data(&mut self, new_data: &[u8]) {
1✔
312
        self.cigar = None;
2✔
313

314
        self.inner_mut().l_data = new_data.len() as i32;
1✔
315
        if (self.inner().m_data as i32) < self.inner().l_data {
1✔
316
            // Verbosity due to lexical borrowing
317
            let l_data = self.inner().l_data;
1✔
318
            self.realloc_var_data(l_data as usize);
1✔
319
        }
320

321
        // Copy new data into buffer
322
        let data =
1✔
323
            unsafe { slice::from_raw_parts_mut(self.inner.data, self.inner().l_data as usize) };
324
        utils::copy_memory(new_data, data);
1✔
325
    }
326

327
    /// Set variable length data (qname, cigar, seq, qual).
328
    /// The aux data is left unchanged.
329
    /// `qual` is Phred-scaled quality values, without any offset.
330
    /// NOTE: seq.len() must equal qual.len() or this method
331
    /// will panic. If you don't have quality values use
332
    /// `let quals = vec![ 255 as u8; seq.len()];` as a placeholder that will
333
    /// be recognized as missing QVs by `samtools`.
334
    pub fn set(&mut self, qname: &[u8], cigar: Option<&CigarString>, seq: &[u8], qual: &[u8]) {
1✔
335
        assert!(qname.len() < 255);
1✔
336
        assert_eq!(seq.len(), qual.len(), "seq.len() must equal qual.len()");
1✔
337

338
        self.cigar = None;
1✔
339

340
        let cigar_width = if let Some(cigar_string) = cigar {
2✔
341
            cigar_string.len()
1✔
342
        } else {
343
            0
×
344
        } * 4;
345
        let q_len = qname.len() + 1;
2✔
346
        let extranul = extranul_from_qname(qname);
1✔
347

348
        let orig_aux_offset = self.qname_capacity()
4✔
349
            + 4 * self.cigar_len()
1✔
350
            + (self.seq_len() + 1) / 2
2✔
351
            + self.seq_len();
1✔
352
        let new_aux_offset = q_len + extranul + cigar_width + (seq.len() + 1) / 2 + qual.len();
2✔
353
        assert!(orig_aux_offset <= self.inner.l_data as usize);
1✔
354
        let aux_len = self.inner.l_data as usize - orig_aux_offset;
2✔
355
        self.inner_mut().l_data = (new_aux_offset + aux_len) as i32;
2✔
356
        if (self.inner().m_data as i32) < self.inner().l_data {
1✔
357
            // Verbosity due to lexical borrowing
358
            let l_data = self.inner().l_data;
1✔
359
            self.realloc_var_data(l_data as usize);
1✔
360
        }
361

362
        // Copy the aux data.
363
        if aux_len > 0 && orig_aux_offset != new_aux_offset {
2✔
364
            let data =
1✔
365
                unsafe { slice::from_raw_parts_mut(self.inner.data, self.inner().m_data as usize) };
366
            data.copy_within(orig_aux_offset..orig_aux_offset + aux_len, new_aux_offset);
1✔
367
        }
368

369
        let data =
1✔
370
            unsafe { slice::from_raw_parts_mut(self.inner.data, self.inner().l_data as usize) };
371

372
        // qname
373
        utils::copy_memory(qname, data);
1✔
374
        for i in 0..=extranul {
4✔
375
            data[qname.len() + i] = b'\0';
2✔
376
        }
377
        let mut i = q_len + extranul;
2✔
378
        self.inner_mut().core.l_qname = i as u16;
1✔
379
        self.inner_mut().core.l_extranul = extranul as u8;
1✔
380

381
        // cigar
382
        if let Some(cigar_string) = cigar {
3✔
383
            let cigar_data = unsafe {
384
                //cigar is always aligned to 4 bytes (see extranul above) - so this is safe
385
                #[allow(clippy::cast_ptr_alignment)]
386
                slice::from_raw_parts_mut(data[i..].as_ptr() as *mut u32, cigar_string.len())
1✔
387
            };
388
            for (i, c) in cigar_string.iter().enumerate() {
4✔
389
                cigar_data[i] = c.encode();
2✔
390
            }
391
            self.inner_mut().core.n_cigar = cigar_string.len() as u32;
1✔
392
            i += cigar_string.len() * 4;
2✔
393
        } else {
394
            self.inner_mut().core.n_cigar = 0;
×
395
        };
396

397
        // seq
398
        {
399
            for j in (0..seq.len()).step_by(2) {
4✔
400
                data[i + j / 2] = ENCODE_BASE[seq[j] as usize] << 4
4✔
401
                    | (if j + 1 < seq.len() {
4✔
402
                        ENCODE_BASE[seq[j + 1] as usize]
2✔
403
                    } else {
404
                        0
1✔
405
                    });
406
            }
407
            self.inner_mut().core.l_qseq = seq.len() as i32;
1✔
408
            i += (seq.len() + 1) / 2;
2✔
409
        }
410

411
        // qual
412
        utils::copy_memory(qual, &mut data[i..]);
1✔
413
    }
414

415
    /// Replace current qname with a new one.
416
    pub fn set_qname(&mut self, new_qname: &[u8]) {
1✔
417
        // 251 + 1NUL is the max 32-bit aligned value that fits in u8
418
        assert!(new_qname.len() < 252);
1✔
419

420
        let old_q_len = self.qname_capacity();
1✔
421
        // We're going to add a terminal NUL
422
        let extranul = extranul_from_qname(new_qname);
1✔
423
        let new_q_len = new_qname.len() + 1 + extranul;
2✔
424

425
        // Length of data after qname
426
        let other_len = self.inner_mut().l_data - old_q_len as i32;
2✔
427

428
        if new_q_len < old_q_len && self.inner().l_data > (old_q_len as i32) {
3✔
429
            self.inner_mut().l_data -= (old_q_len - new_q_len) as i32;
2✔
430
        } else if new_q_len > old_q_len {
1✔
431
            self.inner_mut().l_data += (new_q_len - old_q_len) as i32;
1✔
432

433
            // Reallocate if necessary
434
            if (self.inner().m_data as i32) < self.inner().l_data {
1✔
435
                // Verbosity due to lexical borrowing
436
                let l_data = self.inner().l_data;
1✔
437
                self.realloc_var_data(l_data as usize);
1✔
438
            }
439
        }
440

441
        if new_q_len != old_q_len {
1✔
442
            // Move other data to new location
443
            unsafe {
444
                let data = slice::from_raw_parts_mut(self.inner.data, self.inner().l_data as usize);
1✔
445

446
                ::libc::memmove(
447
                    data.as_mut_ptr().add(new_q_len) as *mut ::libc::c_void,
448
                    data.as_mut_ptr().add(old_q_len) as *mut ::libc::c_void,
449
                    other_len as usize,
1✔
450
                );
451
            }
452
        }
453

454
        // Copy qname data
455
        let data =
1✔
456
            unsafe { slice::from_raw_parts_mut(self.inner.data, self.inner().l_data as usize) };
457
        utils::copy_memory(new_qname, data);
1✔
458
        for i in 0..=extranul {
4✔
459
            data[new_q_len - i - 1] = b'\0';
1✔
460
        }
461
        self.inner_mut().core.l_qname = new_q_len as u16;
1✔
462
        self.inner_mut().core.l_extranul = extranul as u8;
1✔
463
    }
464

465
    fn realloc_var_data(&mut self, new_len: usize) {
1✔
466
        // pad request
467
        let new_len = new_len as u32;
1✔
468
        let new_request = new_len + 32 - (new_len % 32);
1✔
469

470
        let ptr = unsafe {
471
            ::libc::realloc(
1✔
472
                self.inner().data as *mut ::libc::c_void,
1✔
473
                new_request as usize,
1✔
474
            ) as *mut u8
475
        };
476

477
        if ptr.is_null() {
1✔
478
            panic!("ran out of memory in rust_htslib trying to realloc");
×
479
        }
480

481
        // don't update m_data until we know we have
482
        // a successful allocation.
483
        self.inner_mut().m_data = new_request;
1✔
484
        self.inner_mut().data = ptr;
1✔
485

486
        // we now own inner.data
487
        self.own = true;
1✔
488
    }
489

490
    pub fn cigar_len(&self) -> usize {
3✔
491
        self.inner().core.n_cigar as usize
3✔
492
    }
493

494
    /// Get reference to raw cigar string representation (as stored in BAM file).
495
    /// Usually, the method `Record::cigar` should be used instead.
496
    pub fn raw_cigar(&self) -> &[u32] {
1✔
497
        //cigar is always aligned to 4 bytes - so this is safe
498
        #[allow(clippy::cast_ptr_alignment)]
499
        unsafe {
500
            slice::from_raw_parts(
501
                self.data()[self.qname_capacity()..].as_ptr() as *const u32,
1✔
502
                self.cigar_len(),
1✔
503
            )
504
        }
505
    }
506

507
    /// Return unpacked cigar string. This will create a fresh copy the Cigar data.
508
    pub fn cigar(&self) -> CigarStringView {
1✔
509
        match self.cigar {
1✔
510
            Some(ref c) => c.clone(),
×
511
            None => self.unpack_cigar(),
1✔
512
        }
513
    }
514

515
    // Return unpacked cigar string. This returns None unless you have first called `bam::Record::cache_cigar`.
516
    pub fn cigar_cached(&self) -> Option<&CigarStringView> {
×
517
        self.cigar.as_ref()
×
518
    }
519

520
    /// Decode the cigar string and cache it inside the `Record`
521
    pub fn cache_cigar(&mut self) {
×
522
        self.cigar = Some(self.unpack_cigar())
×
523
    }
524

525
    /// Unpack cigar string. Complexity: O(k) with k being the length of the cigar string.
526
    fn unpack_cigar(&self) -> CigarStringView {
1✔
527
        CigarString(
528
            self.raw_cigar()
1✔
529
                .iter()
530
                .map(|&c| {
2✔
531
                    let len = c >> 4;
1✔
532
                    match c & 0b1111 {
1✔
533
                        0 => Cigar::Match(len),
1✔
534
                        1 => Cigar::Ins(len),
1✔
535
                        2 => Cigar::Del(len),
1✔
536
                        3 => Cigar::RefSkip(len),
1✔
537
                        4 => Cigar::SoftClip(len),
1✔
538
                        5 => Cigar::HardClip(len),
1✔
539
                        6 => Cigar::Pad(len),
×
540
                        7 => Cigar::Equal(len),
1✔
541
                        8 => Cigar::Diff(len),
1✔
542
                        _ => panic!("Unexpected cigar operation"),
×
543
                    }
544
                })
545
                .collect(),
546
        )
547
        .into_view(self.pos())
2✔
548
    }
549

550
    pub fn seq_len(&self) -> usize {
3✔
551
        self.inner().core.l_qseq as usize
3✔
552
    }
553

554
    fn seq_data(&self) -> &[u8] {
3✔
555
        let offset = self.qname_capacity() + self.cigar_len() * 4;
3✔
556
        &self.data()[offset..][..(self.seq_len() + 1) / 2]
6✔
557
    }
558

559
    /// Get read sequence. Complexity: O(1).
560
    pub fn seq(&self) -> Seq<'_> {
3✔
561
        Seq {
562
            encoded: self.seq_data(),
3✔
563
            len: self.seq_len(),
3✔
564
        }
565
    }
566

567
    /// Get base qualities (PHRED-scaled probability that base is wrong).
568
    /// This does not entail any offsets, hence the qualities can be used directly without
569
    /// e.g. subtracting 33. Complexity: O(1).
570
    pub fn qual(&self) -> &[u8] {
1✔
571
        &self.data()[self.qname_capacity() + self.cigar_len() * 4 + (self.seq_len() + 1) / 2..]
2✔
572
            [..self.seq_len()]
1✔
573
    }
574

575
    /// Look up an auxiliary field by its tag.
576
    ///
577
    /// Only the first two bytes of a given tag are used for the look-up of a field.
578
    /// See [`Aux`] for more details.
579
    pub fn aux(&self, tag: &[u8]) -> Result<Aux<'_>> {
3✔
580
        let c_str = ffi::CString::new(tag).map_err(|_| Error::BamAuxStringError)?;
3✔
581
        let aux = unsafe {
582
            htslib::bam_aux_get(
583
                &self.inner as *const htslib::bam1_t,
584
                c_str.as_ptr() as *mut c_char,
6✔
585
            )
586
        };
587
        unsafe { Self::read_aux_field(aux).map(|(aux_field, _length)| aux_field) }
9✔
588
    }
589

590
    unsafe fn read_aux_field<'a>(aux: *const u8) -> Result<(Aux<'a>, usize)> {
3✔
591
        const TAG_LEN: isize = 2;
592
        // Used for skipping type identifier
593
        const TYPE_ID_LEN: isize = 1;
594

595
        if aux.is_null() {
3✔
596
            return Err(Error::BamAuxTagNotFound);
3✔
597
        }
598

599
        let (data, type_size) = match *aux {
6✔
600
            b'A' => {
×
601
                let type_size = size_of::<u8>();
×
602
                (Aux::Char(*aux.offset(TYPE_ID_LEN)), type_size)
×
603
            }
604
            b'c' => {
×
605
                let type_size = size_of::<i8>();
1✔
606
                (Aux::I8(*aux.offset(TYPE_ID_LEN).cast::<i8>()), type_size)
1✔
607
            }
608
            b'C' => {
×
609
                let type_size = size_of::<u8>();
1✔
610
                (Aux::U8(*aux.offset(TYPE_ID_LEN)), type_size)
1✔
611
            }
612
            b's' => {
×
613
                let type_size = size_of::<i16>();
1✔
614
                (
615
                    Aux::I16(
1✔
616
                        slice::from_raw_parts(aux.offset(TYPE_ID_LEN), type_size)
2✔
617
                            .read_i16::<LittleEndian>()
×
618
                            .map_err(|_| Error::BamAuxParsingError)?,
×
619
                    ),
620
                    type_size,
×
621
                )
622
            }
623
            b'S' => {
×
624
                let type_size = size_of::<u16>();
1✔
625
                (
626
                    Aux::U16(
1✔
627
                        slice::from_raw_parts(aux.offset(TYPE_ID_LEN), type_size)
2✔
628
                            .read_u16::<LittleEndian>()
×
629
                            .map_err(|_| Error::BamAuxParsingError)?,
×
630
                    ),
631
                    type_size,
×
632
                )
633
            }
634
            b'i' => {
×
635
                let type_size = size_of::<i32>();
2✔
636
                (
637
                    Aux::I32(
2✔
638
                        slice::from_raw_parts(aux.offset(TYPE_ID_LEN), type_size)
4✔
639
                            .read_i32::<LittleEndian>()
×
640
                            .map_err(|_| Error::BamAuxParsingError)?,
×
641
                    ),
642
                    type_size,
×
643
                )
644
            }
645
            b'I' => {
×
646
                let type_size = size_of::<u32>();
1✔
647
                (
648
                    Aux::U32(
1✔
649
                        slice::from_raw_parts(aux.offset(TYPE_ID_LEN), type_size)
2✔
650
                            .read_u32::<LittleEndian>()
×
651
                            .map_err(|_| Error::BamAuxParsingError)?,
×
652
                    ),
653
                    type_size,
×
654
                )
655
            }
656
            b'f' => {
×
657
                let type_size = size_of::<f32>();
1✔
658
                (
659
                    Aux::Float(
1✔
660
                        slice::from_raw_parts(aux.offset(TYPE_ID_LEN), type_size)
2✔
661
                            .read_f32::<LittleEndian>()
×
662
                            .map_err(|_| Error::BamAuxParsingError)?,
×
663
                    ),
664
                    type_size,
×
665
                )
666
            }
667
            b'd' => {
×
668
                let type_size = size_of::<f64>();
1✔
669
                (
670
                    Aux::Double(
1✔
671
                        slice::from_raw_parts(aux.offset(TYPE_ID_LEN), type_size)
2✔
672
                            .read_f64::<LittleEndian>()
×
673
                            .map_err(|_| Error::BamAuxParsingError)?,
×
674
                    ),
675
                    type_size,
×
676
                )
677
            }
678
            b'Z' | b'H' => {
×
679
                let c_str = ffi::CStr::from_ptr(aux.offset(TYPE_ID_LEN).cast::<c_char>());
1✔
680
                let rust_str = c_str.to_str().map_err(|_| Error::BamAuxParsingError)?;
2✔
681
                (Aux::String(rust_str), c_str.to_bytes_with_nul().len())
1✔
682
            }
683
            b'B' => {
×
684
                const ARRAY_INNER_TYPE_LEN: isize = 1;
685
                const ARRAY_COUNT_LEN: isize = 4;
686

687
                // Used for skipping metadata
688
                let array_data_offset = TYPE_ID_LEN + ARRAY_INNER_TYPE_LEN + ARRAY_COUNT_LEN;
6✔
689

690
                let length =
6✔
691
                    slice::from_raw_parts(aux.offset(TYPE_ID_LEN + ARRAY_INNER_TYPE_LEN), 4)
×
692
                        .read_u32::<LittleEndian>()
×
693
                        .map_err(|_| Error::BamAuxParsingError)? as usize;
×
694

695
                // Return tuples of an `Aux` enum and the length of data + metadata in bytes
696
                let (array_data, array_size) = match *aux.offset(TYPE_ID_LEN) {
6✔
697
                    b'c' => (
1✔
698
                        Aux::ArrayI8(AuxArray::<'a, i8>::from_bytes(slice::from_raw_parts(
1✔
699
                            aux.offset(array_data_offset),
×
700
                            length,
×
701
                        ))),
702
                        length,
×
703
                    ),
704
                    b'C' => (
1✔
705
                        Aux::ArrayU8(AuxArray::<'a, u8>::from_bytes(slice::from_raw_parts(
1✔
706
                            aux.offset(array_data_offset),
×
707
                            length,
×
708
                        ))),
709
                        length,
×
710
                    ),
711
                    b's' => (
1✔
712
                        Aux::ArrayI16(AuxArray::<'a, i16>::from_bytes(slice::from_raw_parts(
1✔
713
                            aux.offset(array_data_offset),
×
714
                            length * size_of::<i16>(),
1✔
715
                        ))),
716
                        length * std::mem::size_of::<i16>(),
1✔
717
                    ),
718
                    b'S' => (
1✔
719
                        Aux::ArrayU16(AuxArray::<'a, u16>::from_bytes(slice::from_raw_parts(
1✔
720
                            aux.offset(array_data_offset),
×
721
                            length * size_of::<u16>(),
1✔
722
                        ))),
723
                        length * std::mem::size_of::<u16>(),
1✔
724
                    ),
725
                    b'i' => (
1✔
726
                        Aux::ArrayI32(AuxArray::<'a, i32>::from_bytes(slice::from_raw_parts(
1✔
727
                            aux.offset(array_data_offset),
×
728
                            length * size_of::<i32>(),
1✔
729
                        ))),
730
                        length * std::mem::size_of::<i32>(),
1✔
731
                    ),
732
                    b'I' => (
1✔
733
                        Aux::ArrayU32(AuxArray::<'a, u32>::from_bytes(slice::from_raw_parts(
1✔
734
                            aux.offset(array_data_offset),
×
735
                            length * size_of::<u32>(),
1✔
736
                        ))),
737
                        length * std::mem::size_of::<u32>(),
1✔
738
                    ),
739
                    b'f' => (
3✔
740
                        Aux::ArrayFloat(AuxArray::<f32>::from_bytes(slice::from_raw_parts(
3✔
741
                            aux.offset(array_data_offset),
×
742
                            length * size_of::<f32>(),
3✔
743
                        ))),
744
                        length * std::mem::size_of::<f32>(),
3✔
745
                    ),
746
                    _ => {
×
747
                        return Err(Error::BamAuxUnknownType);
×
748
                    }
749
                };
750
                (
751
                    array_data,
×
752
                    // Offset: array-specific metadata + array size
753
                    ARRAY_INNER_TYPE_LEN as usize + ARRAY_COUNT_LEN as usize + array_size,
6✔
754
                )
755
            }
756
            _ => {
×
757
                return Err(Error::BamAuxUnknownType);
×
758
            }
759
        };
760

761
        // Offset: metadata + type size
762
        Ok((data, TAG_LEN as usize + TYPE_ID_LEN as usize + type_size))
6✔
763
    }
764

765
    /// Returns an iterator over the auxiliary fields of the record.
766
    ///
767
    /// When an error occurs, the `Err` variant will be returned
768
    /// and the iterator will not be able to advance anymore.
769
    pub fn aux_iter(&self) -> AuxIter {
1✔
770
        AuxIter {
771
            // In order to get to the aux data section of a `bam::Record`
772
            // we need to skip fields in front of it
773
            aux: &self.data()[
2✔
774
                // NUL terminated read name:
775
                self.qname_capacity()
776
                // CIGAR (uint32_t):
777
                + self.cigar_len() * std::mem::size_of::<u32>()
778
                // Read sequence (4-bit encoded):
779
                + (self.seq_len() + 1) / 2
780
                // Base qualities (char):
781
                + self.seq_len()..],
782
        }
783
    }
784

785
    /// Add auxiliary data.
786
    pub fn push_aux(&mut self, tag: &[u8], value: Aux<'_>) -> Result<()> {
3✔
787
        // Don't allow pushing aux data when the given tag is already present in the record.
788
        // `htslib` seems to allow this (for non-array values), which can lead to problems
789
        // since retrieving aux fields consumes &[u8; 2] and yields one field only.
790
        if self.aux(tag).is_ok() {
3✔
791
            return Err(Error::BamAuxTagAlreadyPresent);
1✔
792
        }
793

794
        let ctag = tag.as_ptr() as *mut c_char;
3✔
795
        let ret = unsafe {
796
            match value {
3✔
797
                Aux::Char(v) => htslib::bam_aux_append(
798
                    self.inner_ptr_mut(),
×
799
                    ctag,
800
                    b'A' as c_char,
801
                    size_of::<u8>() as i32,
×
802
                    [v].as_mut_ptr() as *mut u8,
×
803
                ),
804
                Aux::I8(v) => htslib::bam_aux_append(
805
                    self.inner_ptr_mut(),
1✔
806
                    ctag,
807
                    b'c' as c_char,
808
                    size_of::<i8>() as i32,
1✔
809
                    [v].as_mut_ptr() as *mut u8,
1✔
810
                ),
811
                Aux::U8(v) => htslib::bam_aux_append(
812
                    self.inner_ptr_mut(),
1✔
813
                    ctag,
814
                    b'C' as c_char,
815
                    size_of::<u8>() as i32,
1✔
816
                    [v].as_mut_ptr() as *mut u8,
1✔
817
                ),
818
                Aux::I16(v) => htslib::bam_aux_append(
819
                    self.inner_ptr_mut(),
1✔
820
                    ctag,
821
                    b's' as c_char,
822
                    size_of::<i16>() as i32,
1✔
823
                    [v].as_mut_ptr() as *mut u8,
1✔
824
                ),
825
                Aux::U16(v) => htslib::bam_aux_append(
826
                    self.inner_ptr_mut(),
1✔
827
                    ctag,
828
                    b'S' as c_char,
829
                    size_of::<u16>() as i32,
1✔
830
                    [v].as_mut_ptr() as *mut u8,
1✔
831
                ),
832
                Aux::I32(v) => htslib::bam_aux_append(
833
                    self.inner_ptr_mut(),
2✔
834
                    ctag,
835
                    b'i' as c_char,
836
                    size_of::<i32>() as i32,
2✔
837
                    [v].as_mut_ptr() as *mut u8,
2✔
838
                ),
839
                Aux::U32(v) => htslib::bam_aux_append(
840
                    self.inner_ptr_mut(),
1✔
841
                    ctag,
842
                    b'I' as c_char,
843
                    size_of::<u32>() as i32,
1✔
844
                    [v].as_mut_ptr() as *mut u8,
1✔
845
                ),
846
                Aux::Float(v) => htslib::bam_aux_append(
847
                    self.inner_ptr_mut(),
1✔
848
                    ctag,
849
                    b'f' as c_char,
850
                    size_of::<f32>() as i32,
1✔
851
                    [v].as_mut_ptr() as *mut u8,
1✔
852
                ),
853
                // Not part of specs but implemented in `htslib`:
854
                Aux::Double(v) => htslib::bam_aux_append(
855
                    self.inner_ptr_mut(),
1✔
856
                    ctag,
857
                    b'd' as c_char,
858
                    size_of::<f64>() as i32,
1✔
859
                    [v].as_mut_ptr() as *mut u8,
1✔
860
                ),
861
                Aux::String(v) => {
1✔
862
                    let c_str = ffi::CString::new(v).map_err(|_| Error::BamAuxStringError)?;
2✔
863
                    htslib::bam_aux_append(
864
                        self.inner_ptr_mut(),
1✔
865
                        ctag,
866
                        b'Z' as c_char,
867
                        (v.len() + 1) as i32,
1✔
868
                        c_str.as_ptr() as *mut u8,
2✔
869
                    )
870
                }
871
                Aux::HexByteArray(v) => {
×
872
                    let c_str = ffi::CString::new(v).map_err(|_| Error::BamAuxStringError)?;
×
873
                    htslib::bam_aux_append(
874
                        self.inner_ptr_mut(),
×
875
                        ctag,
876
                        b'H' as c_char,
877
                        (v.len() + 1) as i32,
×
878
                        c_str.as_ptr() as *mut u8,
×
879
                    )
880
                }
881
                // Not sure it's safe to cast an immutable slice to a mutable pointer in the following branches
882
                Aux::ArrayI8(aux_array) => match aux_array {
1✔
883
                    AuxArray::TargetType(inner) => htslib::bam_aux_update_array(
884
                        self.inner_ptr_mut(),
1✔
885
                        ctag,
886
                        b'c',
887
                        inner.len() as u32,
1✔
888
                        inner.slice.as_ptr() as *mut ::libc::c_void,
1✔
889
                    ),
890
                    AuxArray::RawLeBytes(inner) => htslib::bam_aux_update_array(
891
                        self.inner_ptr_mut(),
1✔
892
                        ctag,
893
                        b'c',
894
                        inner.len() as u32,
1✔
895
                        inner.slice.as_ptr() as *mut ::libc::c_void,
1✔
896
                    ),
897
                },
898
                Aux::ArrayU8(aux_array) => match aux_array {
1✔
899
                    AuxArray::TargetType(inner) => htslib::bam_aux_update_array(
900
                        self.inner_ptr_mut(),
1✔
901
                        ctag,
902
                        b'C',
903
                        inner.len() as u32,
1✔
904
                        inner.slice.as_ptr() as *mut ::libc::c_void,
1✔
905
                    ),
906
                    AuxArray::RawLeBytes(inner) => htslib::bam_aux_update_array(
907
                        self.inner_ptr_mut(),
1✔
908
                        ctag,
909
                        b'C',
910
                        inner.len() as u32,
1✔
911
                        inner.slice.as_ptr() as *mut ::libc::c_void,
1✔
912
                    ),
913
                },
914
                Aux::ArrayI16(aux_array) => match aux_array {
1✔
915
                    AuxArray::TargetType(inner) => htslib::bam_aux_update_array(
916
                        self.inner_ptr_mut(),
1✔
917
                        ctag,
918
                        b's',
919
                        inner.len() as u32,
1✔
920
                        inner.slice.as_ptr() as *mut ::libc::c_void,
1✔
921
                    ),
922
                    AuxArray::RawLeBytes(inner) => htslib::bam_aux_update_array(
923
                        self.inner_ptr_mut(),
1✔
924
                        ctag,
925
                        b's',
926
                        inner.len() as u32,
1✔
927
                        inner.slice.as_ptr() as *mut ::libc::c_void,
1✔
928
                    ),
929
                },
930
                Aux::ArrayU16(aux_array) => match aux_array {
1✔
931
                    AuxArray::TargetType(inner) => htslib::bam_aux_update_array(
932
                        self.inner_ptr_mut(),
1✔
933
                        ctag,
934
                        b'S',
935
                        inner.len() as u32,
1✔
936
                        inner.slice.as_ptr() as *mut ::libc::c_void,
1✔
937
                    ),
938
                    AuxArray::RawLeBytes(inner) => htslib::bam_aux_update_array(
939
                        self.inner_ptr_mut(),
1✔
940
                        ctag,
941
                        b'S',
942
                        inner.len() as u32,
1✔
943
                        inner.slice.as_ptr() as *mut ::libc::c_void,
1✔
944
                    ),
945
                },
946
                Aux::ArrayI32(aux_array) => match aux_array {
1✔
947
                    AuxArray::TargetType(inner) => htslib::bam_aux_update_array(
948
                        self.inner_ptr_mut(),
1✔
949
                        ctag,
950
                        b'i',
951
                        inner.len() as u32,
1✔
952
                        inner.slice.as_ptr() as *mut ::libc::c_void,
1✔
953
                    ),
954
                    AuxArray::RawLeBytes(inner) => htslib::bam_aux_update_array(
955
                        self.inner_ptr_mut(),
1✔
956
                        ctag,
957
                        b'i',
958
                        inner.len() as u32,
1✔
959
                        inner.slice.as_ptr() as *mut ::libc::c_void,
1✔
960
                    ),
961
                },
962
                Aux::ArrayU32(aux_array) => match aux_array {
1✔
963
                    AuxArray::TargetType(inner) => htslib::bam_aux_update_array(
964
                        self.inner_ptr_mut(),
1✔
965
                        ctag,
966
                        b'I',
967
                        inner.len() as u32,
1✔
968
                        inner.slice.as_ptr() as *mut ::libc::c_void,
1✔
969
                    ),
970
                    AuxArray::RawLeBytes(inner) => htslib::bam_aux_update_array(
971
                        self.inner_ptr_mut(),
1✔
972
                        ctag,
973
                        b'I',
974
                        inner.len() as u32,
1✔
975
                        inner.slice.as_ptr() as *mut ::libc::c_void,
1✔
976
                    ),
977
                },
978
                Aux::ArrayFloat(aux_array) => match aux_array {
3✔
979
                    AuxArray::TargetType(inner) => htslib::bam_aux_update_array(
980
                        self.inner_ptr_mut(),
3✔
981
                        ctag,
982
                        b'f',
983
                        inner.len() as u32,
3✔
984
                        inner.slice.as_ptr() as *mut ::libc::c_void,
3✔
985
                    ),
986
                    AuxArray::RawLeBytes(inner) => htslib::bam_aux_update_array(
987
                        self.inner_ptr_mut(),
1✔
988
                        ctag,
989
                        b'f',
990
                        inner.len() as u32,
1✔
991
                        inner.slice.as_ptr() as *mut ::libc::c_void,
1✔
992
                    ),
993
                },
994
            }
995
        };
996

997
        if ret < 0 {
6✔
998
            Err(Error::BamAux)
×
999
        } else {
1000
            Ok(())
3✔
1001
        }
1002
    }
1003

1004
    // Delete auxiliary tag.
1005
    pub fn remove_aux(&mut self, tag: &[u8]) -> Result<()> {
1✔
1006
        let c_str = ffi::CString::new(tag).map_err(|_| Error::BamAuxStringError)?;
1✔
1007
        let aux = unsafe {
1008
            htslib::bam_aux_get(
1009
                &self.inner as *const htslib::bam1_t,
1010
                c_str.as_ptr() as *mut c_char,
2✔
1011
            )
1012
        };
1013
        unsafe {
1014
            if aux.is_null() {
2✔
1015
                Err(Error::BamAuxTagNotFound)
1✔
1016
            } else {
1017
                htslib::bam_aux_del(self.inner_ptr_mut(), aux);
2✔
1018
                Ok(())
1✔
1019
            }
1020
        }
1021
    }
1022

1023
    /// Access the base modifications associated with this Record through the MM tag.
1024
    /// Example:
1025
    /// ```
1026
    ///    use rust_htslib::bam::{Read, Reader, Record};
1027
    ///    let mut bam = Reader::from_path(&"test/base_mods/MM-orient.sam").unwrap();
1028
    ///    let mut mod_count = 0;
1029
    ///    for r in bam.records() {
1030
    ///        let record = r.unwrap();
1031
    ///        if let Ok(mods) = record.basemods_iter() {
1032
    ///            // print metadata for the modifications present in this record
1033
    ///            for mod_code in mods.recorded() {
1034
    ///                if let Ok(mod_metadata) = mods.query_type(*mod_code) {
1035
    ///                    println!("mod found with code {}/{} flags: [{} {} {}]",
1036
    ///                              mod_code, *mod_code as u8 as char,
1037
    ///                              mod_metadata.strand, mod_metadata.implicit, mod_metadata.canonical as u8 as char);
1038
    ///                }
1039
    ///            }
1040
    ///
1041
    ///            // iterate over the modifications in this record
1042
    ///            // the modifications are returned as a tuple with the
1043
    ///            // position within SEQ and an hts_base_mod struct
1044
    ///            for res in mods {
1045
    ///                if let Ok( (position, m) ) = res {
1046
    ///                    println!("{} {},{}", position, m.modified_base as u8 as char, m.qual);
1047
    ///                    mod_count += 1;
1048
    ///                }
1049
    ///            }
1050
    ///        };
1051
    ///    }
1052
    ///    assert_eq!(mod_count, 14);
1053
    /// ```
1054
    pub fn basemods_iter(&self) -> Result<BaseModificationsIter> {
2✔
1055
        BaseModificationsIter::new(self)
2✔
1056
    }
1057

1058
    /// An iterator that returns all of the modifications for each position as a vector.
1059
    /// This is useful for the case where multiple possible modifications can be annotated
1060
    /// at a single position (for example a C could be 5-mC or 5-hmC)
1061
    pub fn basemods_position_iter(&self) -> Result<BaseModificationsPositionIter> {
1✔
1062
        BaseModificationsPositionIter::new(self)
1✔
1063
    }
1064

1065
    /// Infer read pair orientation from record. Returns `SequenceReadPairOrientation::None` if record
1066
    /// is not paired, mates are not mapping to the same contig, or mates start at the
1067
    /// same position.
1068
    pub fn read_pair_orientation(&self) -> SequenceReadPairOrientation {
1✔
1069
        if self.is_paired()
4✔
1070
            && !self.is_unmapped()
1✔
1071
            && !self.is_mate_unmapped()
1✔
1072
            && self.tid() == self.mtid()
1✔
1073
        {
1074
            if self.pos() == self.mpos() {
1✔
1075
                // both reads start at the same position, we cannot decide on the orientation.
1076
                return SequenceReadPairOrientation::None;
×
1077
            }
1078

1079
            let (is_reverse, is_first_in_template, is_mate_reverse) = if self.pos() < self.mpos() {
3✔
1080
                // given record is the left one
1081
                (
1082
                    self.is_reverse(),
1✔
1083
                    self.is_first_in_template(),
1✔
1084
                    self.is_mate_reverse(),
1✔
1085
                )
1086
            } else {
1087
                // given record is the right one
1088
                (
1089
                    self.is_mate_reverse(),
1✔
1090
                    self.is_last_in_template(),
1✔
1091
                    self.is_reverse(),
1✔
1092
                )
1093
            };
1094
            match (is_reverse, is_first_in_template, is_mate_reverse) {
1✔
1095
                (false, false, false) => SequenceReadPairOrientation::F2F1,
×
1096
                (false, false, true) => SequenceReadPairOrientation::F2R1,
×
1097
                (false, true, false) => SequenceReadPairOrientation::F1F2,
×
1098
                (true, false, false) => SequenceReadPairOrientation::R2F1,
×
1099
                (false, true, true) => SequenceReadPairOrientation::F1R2,
1✔
1100
                (true, false, true) => SequenceReadPairOrientation::R2R1,
×
1101
                (true, true, false) => SequenceReadPairOrientation::R1F2,
×
1102
                (true, true, true) => SequenceReadPairOrientation::R1R2,
×
1103
            }
1104
        } else {
1105
            SequenceReadPairOrientation::None
×
1106
        }
1107
    }
1108

1109
    flag!(is_paired, set_paired, unset_paired, 1u16);
1110
    flag!(is_proper_pair, set_proper_pair, unset_proper_pair, 2u16);
1111
    flag!(is_unmapped, set_unmapped, unset_unmapped, 4u16);
1112
    flag!(
1113
        is_mate_unmapped,
1114
        set_mate_unmapped,
1115
        unset_mate_unmapped,
1116
        8u16
1117
    );
1118
    flag!(is_reverse, set_reverse, unset_reverse, 16u16);
1119
    flag!(is_mate_reverse, set_mate_reverse, unset_mate_reverse, 32u16);
1120
    flag!(
1121
        is_first_in_template,
1122
        set_first_in_template,
1123
        unset_first_in_template,
1124
        64u16
1125
    );
1126
    flag!(
1127
        is_last_in_template,
1128
        set_last_in_template,
1129
        unset_last_in_template,
1130
        128u16
1131
    );
1132
    flag!(is_secondary, set_secondary, unset_secondary, 256u16);
1133
    flag!(
1134
        is_quality_check_failed,
1135
        set_quality_check_failed,
1136
        unset_quality_check_failed,
1137
        512u16
1138
    );
1139
    flag!(is_duplicate, set_duplicate, unset_duplicate, 1024u16);
1140
    flag!(
1141
        is_supplementary,
1142
        set_supplementary,
1143
        unset_supplementary,
1144
        2048u16
1145
    );
1146
}
1147

1148
impl Drop for Record {
1149
    fn drop(&mut self) {
10✔
1150
        if self.own {
10✔
1151
            unsafe { ::libc::free(self.inner.data as *mut ::libc::c_void) }
9✔
1152
        }
1153
    }
1154
}
1155

1156
impl SequenceRead for Record {
1157
    fn name(&self) -> &[u8] {
×
1158
        self.qname()
×
1159
    }
1160

1161
    fn base(&self, i: usize) -> u8 {
×
1162
        *decode_base_unchecked(encoded_base(self.seq_data(), i))
×
1163
    }
1164

1165
    fn base_qual(&self, i: usize) -> u8 {
×
1166
        self.qual()[i]
×
1167
    }
1168

1169
    fn len(&self) -> usize {
×
1170
        self.seq_len()
×
1171
    }
1172

1173
    fn is_empty(&self) -> bool {
×
1174
        self.len() == 0
×
1175
    }
1176
}
1177

1178
impl genome::AbstractInterval for Record {
1179
    /// Return contig name. Panics if record does not know its header (which happens if it has not been read from a file).
1180
    fn contig(&self) -> &str {
×
1181
        let tid = self.tid();
×
1182
        if tid < 0 {
×
1183
            panic!("invalid tid, must be at least zero");
×
1184
        }
1185
        str::from_utf8(
1186
            self.header
×
1187
                .as_ref()
1188
                .expect(
1189
                    "header must be set (this is the case if the record has been read from a file)",
1190
                )
1191
                .tid2name(tid as u32),
1192
        )
1193
        .expect("unable to interpret contig name as UTF-8")
1194
    }
1195

1196
    /// Return genomic range covered by alignment. Panics if `Record::cache_cigar()` has not been called first or `Record::pos()` is less than zero.
1197
    fn range(&self) -> ops::Range<genome::Position> {
×
1198
        let end_pos = self
×
1199
            .cigar_cached()
1200
            .expect("cigar has not been cached yet, call cache_cigar() first")
1201
            .end_pos() as u64;
1202

1203
        if self.pos() < 0 {
×
1204
            panic!("invalid position, must be positive")
×
1205
        }
1206

1207
        self.pos() as u64..end_pos
×
1208
    }
1209
}
1210

1211
/// Auxiliary record data
1212
///
1213
/// The specification allows a wide range of types to be stored as an auxiliary data field of a BAM record.
1214
///
1215
/// Please note that the [`Aux::Double`] variant is _not_ part of the specification, but it is supported by `htslib`.
1216
///
1217
/// # Examples
1218
///
1219
/// ```
1220
/// use rust_htslib::{
1221
///     bam,
1222
///     bam::record::{Aux, AuxArray},
1223
///     errors::Error,
1224
/// };
1225
///
1226
/// //Set up BAM record
1227
/// let bam_header = bam::Header::new();
1228
/// let mut record = bam::Record::from_sam(
1229
///     &mut bam::HeaderView::from_header(&bam_header),
1230
///     "ali1\t4\t*\t0\t0\t*\t*\t0\t0\tACGT\tFFFF".as_bytes(),
1231
/// )
1232
/// .unwrap();
1233
///
1234
/// // Add an integer field
1235
/// let aux_integer_field = Aux::I32(1234);
1236
/// record.push_aux(b"XI", aux_integer_field).unwrap();
1237
///
1238
/// match record.aux(b"XI") {
1239
///     Ok(value) => {
1240
///         // Typically, callers expect an aux field to be of a certain type.
1241
///         // If that's not the case, the value can be `match`ed exhaustively.
1242
///         if let Aux::I32(v) = value {
1243
///             assert_eq!(v, 1234);
1244
///         }
1245
///     }
1246
///     Err(e) => {
1247
///         panic!("Error reading aux field: {}", e);
1248
///     }
1249
/// }
1250
///
1251
/// // Add an array field
1252
/// let array_like_data = vec![0.4, 0.3, 0.2, 0.1];
1253
/// let slice_of_data = &array_like_data;
1254
/// let aux_array: AuxArray<f32> = slice_of_data.into();
1255
/// let aux_array_field = Aux::ArrayFloat(aux_array);
1256
/// record.push_aux(b"XA", aux_array_field).unwrap();
1257
///
1258
/// if let Ok(Aux::ArrayFloat(array)) = record.aux(b"XA") {
1259
///     let read_array = array.iter().collect::<Vec<_>>();
1260
///     assert_eq!(read_array, array_like_data);
1261
/// } else {
1262
///     panic!("Could not read array data");
1263
/// }
1264
/// ```
1265
#[derive(Debug, PartialEq)]
1266
pub enum Aux<'a> {
1267
    Char(u8),
1268
    I8(i8),
1269
    U8(u8),
1270
    I16(i16),
1271
    U16(u16),
1272
    I32(i32),
1273
    U32(u32),
1274
    Float(f32),
1275
    Double(f64), // Not part of specs but implemented in `htslib`
1276
    String(&'a str),
1277
    HexByteArray(&'a str),
1278
    ArrayI8(AuxArray<'a, i8>),
1279
    ArrayU8(AuxArray<'a, u8>),
1280
    ArrayI16(AuxArray<'a, i16>),
1281
    ArrayU16(AuxArray<'a, u16>),
1282
    ArrayI32(AuxArray<'a, i32>),
1283
    ArrayU32(AuxArray<'a, u32>),
1284
    ArrayFloat(AuxArray<'a, f32>),
1285
}
1286

1287
unsafe impl<'a> Send for Aux<'a> {}
1288
unsafe impl<'a> Sync for Aux<'a> {}
1289

1290
/// Types that can be used in aux arrays.
1291
pub trait AuxArrayElement: Copy {
1292
    fn from_le_bytes(bytes: &[u8]) -> Option<Self>;
1293
}
1294

1295
impl AuxArrayElement for i8 {
1296
    fn from_le_bytes(bytes: &[u8]) -> Option<Self> {
1✔
1297
        std::io::Cursor::new(bytes).read_i8().ok()
1✔
1298
    }
1299
}
1300
impl AuxArrayElement for u8 {
1301
    fn from_le_bytes(bytes: &[u8]) -> Option<Self> {
1✔
1302
        std::io::Cursor::new(bytes).read_u8().ok()
1✔
1303
    }
1304
}
1305
impl AuxArrayElement for i16 {
1306
    fn from_le_bytes(bytes: &[u8]) -> Option<Self> {
1✔
1307
        std::io::Cursor::new(bytes).read_i16::<LittleEndian>().ok()
1✔
1308
    }
1309
}
1310
impl AuxArrayElement for u16 {
1311
    fn from_le_bytes(bytes: &[u8]) -> Option<Self> {
1✔
1312
        std::io::Cursor::new(bytes).read_u16::<LittleEndian>().ok()
1✔
1313
    }
1314
}
1315
impl AuxArrayElement for i32 {
1316
    fn from_le_bytes(bytes: &[u8]) -> Option<Self> {
1✔
1317
        std::io::Cursor::new(bytes).read_i32::<LittleEndian>().ok()
1✔
1318
    }
1319
}
1320
impl AuxArrayElement for u32 {
1321
    fn from_le_bytes(bytes: &[u8]) -> Option<Self> {
1✔
1322
        std::io::Cursor::new(bytes).read_u32::<LittleEndian>().ok()
1✔
1323
    }
1324
}
1325
impl AuxArrayElement for f32 {
1326
    fn from_le_bytes(bytes: &[u8]) -> Option<Self> {
3✔
1327
        std::io::Cursor::new(bytes).read_f32::<LittleEndian>().ok()
3✔
1328
    }
1329
}
1330

1331
/// Provides access to aux arrays.
1332
///
1333
/// Provides methods to either retrieve single elements or an iterator over the
1334
/// array.
1335
///
1336
/// This type is used for wrapping both, array data that was read from a
1337
/// BAM record and slices of data that are going to be stored in one.
1338
///
1339
/// In order to be able to add an `AuxArray` field to a BAM record, `AuxArray`s
1340
/// can be constructed via the `From` trait which is implemented for all
1341
/// supported types (see [`AuxArrayElement`] for a list).
1342
///
1343
/// # Examples
1344
///
1345
/// ```
1346
/// use rust_htslib::{
1347
///     bam,
1348
///     bam::record::{Aux, AuxArray},
1349
/// };
1350
///
1351
/// //Set up BAM record
1352
/// let bam_header = bam::Header::new();
1353
/// let mut record = bam::Record::from_sam(
1354
///     &mut bam::HeaderView::from_header(&bam_header),
1355
///     "ali1\t4\t*\t0\t0\t*\t*\t0\t0\tACGT\tFFFF".as_bytes(),
1356
/// ).unwrap();
1357
///
1358
/// let data = vec![0.4, 0.3, 0.2, 0.1];
1359
/// let slice_of_data = &data;
1360
/// let aux_array: AuxArray<f32> = slice_of_data.into();
1361
/// let aux_field = Aux::ArrayFloat(aux_array);
1362
/// record.push_aux(b"XA", aux_field);
1363
///
1364
/// if let Ok(Aux::ArrayFloat(array)) = record.aux(b"XA") {
1365
///     // Retrieve the second element from the array
1366
///     assert_eq!(array.get(1).unwrap(), 0.3);
1367
///     // Iterate over the array and collect it into a `Vec`
1368
///     let read_array = array.iter().collect::<Vec<_>>();
1369
///     assert_eq!(read_array, data);
1370
/// } else {
1371
///     panic!("Could not read array data");
1372
/// }
1373
/// ```
1374
#[derive(Debug)]
1375
pub enum AuxArray<'a, T> {
1376
    TargetType(AuxArrayTargetType<'a, T>),
1377
    RawLeBytes(AuxArrayRawLeBytes<'a, T>),
1378
}
1379

1380
impl<T> PartialEq<AuxArray<'_, T>> for AuxArray<'_, T>
1381
where
1382
    T: AuxArrayElement + PartialEq,
1383
{
1384
    fn eq(&self, other: &AuxArray<'_, T>) -> bool {
7✔
1385
        use AuxArray::*;
1386
        match (self, other) {
7✔
1387
            (TargetType(v), TargetType(v_other)) => v == v_other,
1✔
1388
            (RawLeBytes(v), RawLeBytes(v_other)) => v == v_other,
1✔
1389
            (TargetType(_), RawLeBytes(_)) => self.iter().eq(other.iter()),
×
1390
            (RawLeBytes(_), TargetType(_)) => self.iter().eq(other.iter()),
7✔
1391
        }
1392
    }
1393
}
1394

1395
/// Create AuxArrays from slices of allowed target types.
1396
impl<'a, I, T> From<&'a T> for AuxArray<'a, I>
1397
where
1398
    I: AuxArrayElement,
1399
    T: AsRef<[I]> + ?Sized,
1400
{
1401
    fn from(src: &'a T) -> Self {
9✔
1402
        AuxArray::TargetType(AuxArrayTargetType {
9✔
1403
            slice: src.as_ref(),
9✔
1404
        })
1405
    }
1406
}
1407

1408
impl<'a, T> AuxArray<'a, T>
1409
where
1410
    T: AuxArrayElement,
1411
{
1412
    /// Returns the element at a position or None if out of bounds.
1413
    pub fn get(&self, index: usize) -> Option<T> {
9✔
1414
        match self {
9✔
1415
            AuxArray::TargetType(v) => v.get(index),
7✔
1416
            AuxArray::RawLeBytes(v) => v.get(index),
9✔
1417
        }
1418
    }
1419

1420
    /// Returns the number of elements in the array.
1421
    pub fn len(&self) -> usize {
9✔
1422
        match self {
9✔
1423
            AuxArray::TargetType(a) => a.len(),
×
1424
            AuxArray::RawLeBytes(a) => a.len(),
9✔
1425
        }
1426
    }
1427

1428
    /// Returns true if the array contains no elements.
1429
    pub fn is_empty(&self) -> bool {
×
1430
        self.len() == 0
×
1431
    }
1432

1433
    /// Returns an iterator over the array.
1434
    pub fn iter(&self) -> AuxArrayIter<T> {
9✔
1435
        AuxArrayIter {
1436
            index: 0,
1437
            array: self,
1438
        }
1439
    }
1440

1441
    /// Create AuxArrays from raw byte slices borrowed from `bam::Record`.
1442
    fn from_bytes(bytes: &'a [u8]) -> Self {
9✔
1443
        Self::RawLeBytes(AuxArrayRawLeBytes {
9✔
1444
            slice: bytes,
×
1445
            phantom_data: PhantomData,
×
1446
        })
1447
    }
1448
}
1449

1450
/// Encapsulates slice of target type.
1451
#[doc(hidden)]
1452
#[derive(Debug, PartialEq)]
1453
pub struct AuxArrayTargetType<'a, T> {
1454
    slice: &'a [T],
1455
}
1456

1457
impl<'a, T> AuxArrayTargetType<'a, T>
1458
where
1459
    T: AuxArrayElement,
1460
{
1461
    fn get(&self, index: usize) -> Option<T> {
7✔
1462
        self.slice.get(index).copied()
7✔
1463
    }
1464

1465
    fn len(&self) -> usize {
9✔
1466
        self.slice.len()
9✔
1467
    }
1468
}
1469

1470
/// Encapsulates slice of raw bytes to prevent it from being accidentally accessed.
1471
#[doc(hidden)]
1472
#[derive(Debug, PartialEq)]
1473
pub struct AuxArrayRawLeBytes<'a, T> {
1474
    slice: &'a [u8],
1475
    phantom_data: PhantomData<T>,
1476
}
1477

1478
impl<'a, T> AuxArrayRawLeBytes<'a, T>
1479
where
1480
    T: AuxArrayElement,
1481
{
1482
    fn get(&self, index: usize) -> Option<T> {
9✔
1483
        let type_size = std::mem::size_of::<T>();
9✔
1484
        if index * type_size + type_size > self.slice.len() {
9✔
1485
            return None;
9✔
1486
        }
1487
        T::from_le_bytes(&self.slice[index * type_size..][..type_size])
18✔
1488
    }
1489

1490
    fn len(&self) -> usize {
9✔
1491
        self.slice.len() / std::mem::size_of::<T>()
9✔
1492
    }
1493
}
1494

1495
/// Aux array iterator
1496
///
1497
/// This struct is created by the [`AuxArray::iter`] method.
1498
pub struct AuxArrayIter<'a, T> {
1499
    index: usize,
1500
    array: &'a AuxArray<'a, T>,
1501
}
1502

1503
impl<T> Iterator for AuxArrayIter<'_, T>
1504
where
1505
    T: AuxArrayElement,
1506
{
1507
    type Item = T;
1508
    fn next(&mut self) -> Option<Self::Item> {
9✔
1509
        let value = self.array.get(self.index);
9✔
1510
        self.index += 1;
9✔
1511
        value
×
1512
    }
1513

1514
    fn size_hint(&self) -> (usize, Option<usize>) {
9✔
1515
        let array_length = self.array.len() - self.index;
9✔
1516
        (array_length, Some(array_length))
9✔
1517
    }
1518
}
1519

1520
/// Auxiliary data iterator
1521
///
1522
/// This struct is created by the [`Record::aux_iter`] method.
1523
///
1524
/// This iterator returns `Result`s that wrap tuples containing
1525
/// a slice which represents the two-byte tag (`&[u8; 2]`) as
1526
/// well as an `Aux` enum that wraps the associated value.
1527
///
1528
/// When an error occurs, the `Err` variant will be returned
1529
/// and the iterator will not be able to advance anymore.
1530
pub struct AuxIter<'a> {
1531
    aux: &'a [u8],
1532
}
1533

1534
impl<'a> Iterator for AuxIter<'a> {
1535
    type Item = Result<(&'a [u8], Aux<'a>)>;
1536

1537
    fn next(&mut self) -> Option<Self::Item> {
1✔
1538
        // We're finished
1539
        if self.aux.is_empty() {
1✔
1540
            return None;
1✔
1541
        }
1542
        // Incomplete aux data
1543
        if (1..=3).contains(&self.aux.len()) {
1✔
1544
            // In the case of an error, we can not safely advance in the aux data, so we terminate the Iteration
1545
            self.aux = &[];
×
1546
            return Some(Err(Error::BamAuxParsingError));
×
1547
        }
1548
        let tag = &self.aux[..2];
1✔
1549
        Some(unsafe {
1✔
1550
            let data_ptr = self.aux[2..].as_ptr();
1✔
1551
            Record::read_aux_field(data_ptr)
3✔
1552
                .map(|(aux, offset)| {
3✔
1553
                    self.aux = &self.aux[offset..];
1✔
1554
                    (tag, aux)
1✔
1555
                })
1556
                .map_err(|e| {
1✔
1557
                    // In the case of an error, we can not safely advance in the aux data, so we terminate the Iteration
1558
                    self.aux = &[];
×
1559
                    e
×
1560
                })
1561
        })
1562
    }
1563
}
1564

1565
static DECODE_BASE: &[u8] = b"=ACMGRSVTWYHKDBN";
1566
static ENCODE_BASE: [u8; 256] = [
1567
    15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15,
1568
    15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15,
1569
    1, 2, 4, 8, 15, 15, 15, 15, 15, 15, 15, 15, 15, 0, 15, 15, 15, 1, 14, 2, 13, 15, 15, 4, 11, 15,
1570
    15, 12, 15, 3, 15, 15, 15, 15, 5, 6, 8, 15, 7, 9, 15, 10, 15, 15, 15, 15, 15, 15, 15, 1, 14, 2,
1571
    13, 15, 15, 4, 11, 15, 15, 12, 15, 3, 15, 15, 15, 15, 5, 6, 8, 15, 7, 9, 15, 10, 15, 15, 15,
1572
    15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15,
1573
    15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15,
1574
    15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15,
1575
    15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15,
1576
    15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15,
1577
    15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15,
1578
];
1579

1580
#[inline]
1581
fn encoded_base(encoded_seq: &[u8], i: usize) -> u8 {
3✔
1582
    (encoded_seq[i / 2] >> ((!i & 1) << 2)) & 0b1111
3✔
1583
}
1584

1585
#[inline]
1586
unsafe fn encoded_base_unchecked(encoded_seq: &[u8], i: usize) -> u8 {
×
1587
    (encoded_seq.get_unchecked(i / 2) >> ((!i & 1) << 2)) & 0b1111
×
1588
}
1589

1590
#[inline]
1591
fn decode_base_unchecked(base: u8) -> &'static u8 {
3✔
1592
    unsafe { DECODE_BASE.get_unchecked(base as usize) }
3✔
1593
}
1594

1595
/// The sequence of a record.
1596
#[derive(Debug, Copy, Clone)]
1597
pub struct Seq<'a> {
1598
    pub encoded: &'a [u8],
1599
    len: usize,
1600
}
1601

1602
impl<'a> Seq<'a> {
1603
    /// Return encoded base. Complexity: O(1).
1604
    #[inline]
1605
    pub fn encoded_base(&self, i: usize) -> u8 {
3✔
1606
        encoded_base(self.encoded, i)
3✔
1607
    }
1608

1609
    /// Return encoded base. Complexity: O(1).
1610
    #[inline]
1611
    pub unsafe fn encoded_base_unchecked(&self, i: usize) -> u8 {
×
1612
        encoded_base_unchecked(self.encoded, i)
×
1613
    }
1614

1615
    /// Obtain decoded base without performing bounds checking.
1616
    /// Use index based access seq()[i], for checked, safe access.
1617
    /// Complexity: O(1).
1618
    #[inline]
1619
    pub unsafe fn decoded_base_unchecked(&self, i: usize) -> u8 {
×
1620
        *decode_base_unchecked(self.encoded_base_unchecked(i))
×
1621
    }
1622

1623
    /// Return decoded sequence. Complexity: O(m) with m being the read length.
1624
    pub fn as_bytes(&self) -> Vec<u8> {
2✔
1625
        (0..self.len()).map(|i| self[i]).collect()
6✔
1626
    }
1627

1628
    /// Return length (in bases) of the sequence.
1629
    pub fn len(&self) -> usize {
2✔
1630
        self.len
2✔
1631
    }
1632

1633
    pub fn is_empty(&self) -> bool {
×
1634
        self.len() == 0
×
1635
    }
1636
}
1637

1638
impl<'a> ops::Index<usize> for Seq<'a> {
1639
    type Output = u8;
1640

1641
    /// Return decoded base at given position within read. Complexity: O(1).
1642
    fn index(&self, index: usize) -> &u8 {
3✔
1643
        decode_base_unchecked(self.encoded_base(index))
3✔
1644
    }
1645
}
1646

1647
unsafe impl<'a> Send for Seq<'a> {}
1648
unsafe impl<'a> Sync for Seq<'a> {}
1649

1650
#[cfg_attr(feature = "serde_feature", derive(Serialize, Deserialize))]
1651
#[derive(PartialEq, PartialOrd, Eq, Debug, Clone, Copy, Hash)]
1652
pub enum Cigar {
1653
    Match(u32),    // M
1654
    Ins(u32),      // I
1655
    Del(u32),      // D
1656
    RefSkip(u32),  // N
1657
    SoftClip(u32), // S
1658
    HardClip(u32), // H
1659
    Pad(u32),      // P
1660
    Equal(u32),    // =
1661
    Diff(u32),     // X
1662
}
1663

1664
impl Cigar {
1665
    fn encode(self) -> u32 {
1✔
1666
        match self {
1✔
1667
            Cigar::Match(len) => len << 4, // | 0,
1✔
1668
            Cigar::Ins(len) => len << 4 | 1,
1✔
1669
            Cigar::Del(len) => len << 4 | 2,
1✔
1670
            Cigar::RefSkip(len) => len << 4 | 3,
×
1671
            Cigar::SoftClip(len) => len << 4 | 4,
1✔
1672
            Cigar::HardClip(len) => len << 4 | 5,
1✔
1673
            Cigar::Pad(len) => len << 4 | 6,
×
1674
            Cigar::Equal(len) => len << 4 | 7,
1✔
1675
            Cigar::Diff(len) => len << 4 | 8,
1✔
1676
        }
1677
    }
1678

1679
    /// Return the length of the CIGAR.
1680
    pub fn len(self) -> u32 {
2✔
1681
        match self {
2✔
1682
            Cigar::Match(len) => len,
2✔
1683
            Cigar::Ins(len) => len,
1✔
1684
            Cigar::Del(len) => len,
1✔
1685
            Cigar::RefSkip(len) => len,
×
1686
            Cigar::SoftClip(len) => len,
2✔
1687
            Cigar::HardClip(len) => len,
1✔
1688
            Cigar::Pad(len) => len,
×
1689
            Cigar::Equal(len) => len,
1✔
1690
            Cigar::Diff(len) => len,
1✔
1691
        }
1692
    }
1693

1694
    pub fn is_empty(self) -> bool {
×
1695
        self.len() == 0
×
1696
    }
1697

1698
    /// Return the character representing the CIGAR.
1699
    pub fn char(self) -> char {
2✔
1700
        match self {
2✔
1701
            Cigar::Match(_) => 'M',
2✔
1702
            Cigar::Ins(_) => 'I',
1✔
1703
            Cigar::Del(_) => 'D',
1✔
1704
            Cigar::RefSkip(_) => 'N',
×
1705
            Cigar::SoftClip(_) => 'S',
2✔
1706
            Cigar::HardClip(_) => 'H',
1✔
1707
            Cigar::Pad(_) => 'P',
×
1708
            Cigar::Equal(_) => '=',
1✔
1709
            Cigar::Diff(_) => 'X',
1✔
1710
        }
1711
    }
1712
}
1713

1714
impl fmt::Display for Cigar {
1715
    fn fmt(&self, fmt: &mut fmt::Formatter<'_>) -> Result<(), fmt::Error> {
2✔
1716
        fmt.write_fmt(format_args!("{}{}", self.len(), self.char()))
6✔
1717
    }
1718
}
1719

1720
unsafe impl Send for Cigar {}
1721
unsafe impl Sync for Cigar {}
1722

1723
custom_derive! {
1724
    /// A CIGAR string. This type wraps around a `Vec<Cigar>`.
1725
    ///
1726
    /// # Example
1727
    ///
1728
    /// ```
1729
    /// use rust_htslib::bam::record::{Cigar, CigarString};
1730
    ///
1731
    /// let cigar = CigarString(vec![Cigar::Match(100), Cigar::SoftClip(10)]);
1732
    ///
1733
    /// // access by index
1734
    /// assert_eq!(cigar[0], Cigar::Match(100));
1735
    /// // format into classical string representation
1736
    /// assert_eq!(format!("{}", cigar), "100M10S");
1737
    /// // iterate
1738
    /// for op in &cigar {
1739
    ///    println!("{}", op);
1740
    /// }
1741
    /// ```
1742
    #[cfg_attr(feature = "serde_feature", derive(Serialize, Deserialize))]
×
1743
    #[derive(NewtypeDeref,
1744
            NewtypeDerefMut,
1745
             NewtypeIndex(usize),
1746
             NewtypeIndexMut(usize),
1747
             NewtypeFrom,
1748
             PartialEq,
1749
             PartialOrd,
1750
             Eq,
1751
             NewtypeDebug,
1752
             Clone,
1753
             Hash
1754
    )]
1755
    pub struct CigarString(pub Vec<Cigar>);
3✔
1756
}
1757

1758
impl CigarString {
1759
    /// Create a `CigarStringView` from this CigarString at position `pos`
1760
    pub fn into_view(self, pos: i64) -> CigarStringView {
1✔
1761
        CigarStringView::new(self, pos)
1✔
1762
    }
1763

1764
    /// Calculate the bam cigar from the alignment struct. x is the target string
1765
    /// and y is the reference. `hard_clip` controls how unaligned read bases are encoded in the
1766
    /// cigar string. Set to true to use the hard clip (`H`) code, or false to use soft clip
1767
    /// (`S`) code. See the [SAM spec](https://samtools.github.io/hts-specs/SAMv1.pdf) for more details.
1768
    pub fn from_alignment(alignment: &Alignment, hard_clip: bool) -> Self {
1✔
1769
        match alignment.mode {
1✔
1770
            AlignmentMode::Global => {
1771
                panic!(" Bam cigar fn not supported for Global Alignment mode")
×
1772
            }
1773
            AlignmentMode::Local => panic!(" Bam cigar fn not supported for Local Alignment mode"),
×
1774
            _ => {}
1775
        }
1776

1777
        let mut cigar = Vec::new();
1✔
1778
        if alignment.operations.is_empty() {
2✔
1779
            return CigarString(cigar);
×
1780
        }
1781

1782
        let add_op = |op: AlignmentOperation, length: u32, cigar: &mut Vec<Cigar>| match op {
3✔
1783
            AlignmentOperation::Del => cigar.push(Cigar::Del(length)),
1✔
1784
            AlignmentOperation::Ins => cigar.push(Cigar::Ins(length)),
1✔
1785
            AlignmentOperation::Subst => cigar.push(Cigar::Diff(length)),
1✔
1786
            AlignmentOperation::Match => cigar.push(Cigar::Equal(length)),
1✔
1787
            _ => {}
1788
        };
1789

1790
        if alignment.xstart > 0 {
1✔
1791
            cigar.push(if hard_clip {
3✔
1792
                Cigar::HardClip(alignment.xstart as u32)
×
1793
            } else {
1794
                Cigar::SoftClip(alignment.xstart as u32)
1✔
1795
            });
1796
        }
1797

1798
        let mut last = alignment.operations[0];
2✔
1799
        let mut k = 1u32;
1✔
1800
        for &op in alignment.operations[1..].iter() {
3✔
1801
            if op == last {
3✔
1802
                k += 1;
2✔
1803
            } else {
1804
                add_op(last, k, &mut cigar);
1✔
1805
                k = 1;
1✔
1806
            }
1807
            last = op;
1✔
1808
        }
1809
        add_op(last, k, &mut cigar);
1✔
1810
        if alignment.xlen > alignment.xend {
1✔
1811
            cigar.push(if hard_clip {
5✔
1812
                Cigar::HardClip((alignment.xlen - alignment.xend) as u32)
2✔
1813
            } else {
1814
                Cigar::SoftClip((alignment.xlen - alignment.xend) as u32)
2✔
1815
            });
1816
        }
1817

1818
        CigarString(cigar)
1✔
1819
    }
1820
}
1821

1822
impl TryFrom<&[u8]> for CigarString {
1823
    type Error = Error;
1824

1825
    /// Create a CigarString from given &[u8].
1826
    /// # Example
1827
    /// ```
1828
    /// use rust_htslib::bam::record::*;
1829
    /// use rust_htslib::bam::record::CigarString;
1830
    /// use rust_htslib::bam::record::Cigar::*;
1831
    /// use std::convert::TryFrom;
1832
    ///
1833
    /// let cigar_str = "2H10M5X3=2H".as_bytes();
1834
    /// let cigar = CigarString::try_from(cigar_str)
1835
    ///     .expect("Unable to parse cigar string.");
1836
    /// let expected_cigar = CigarString(vec![
1837
    ///     HardClip(2),
1838
    ///     Match(10),
1839
    ///     Diff(5),
1840
    ///     Equal(3),
1841
    ///     HardClip(2),
1842
    /// ]);
1843
    /// assert_eq!(cigar, expected_cigar);
1844
    /// ```
1845
    fn try_from(bytes: &[u8]) -> Result<Self> {
3✔
1846
        let mut inner = Vec::new();
3✔
1847
        let mut i = 0;
3✔
1848
        let text_len = bytes.len();
3✔
1849
        while i < text_len {
6✔
1850
            let mut j = i;
3✔
1851
            while j < text_len && bytes[j].is_ascii_digit() {
12✔
1852
                j += 1;
6✔
1853
            }
1854
            // check that length is provided
1855
            if i == j {
3✔
1856
                return Err(Error::BamParseCigar {
×
1857
                    msg: "Expected length before cigar operation [0-9]+[MIDNSHP=X]".to_owned(),
×
1858
                });
1859
            }
1860
            // get the length of the operation
1861
            let s = str::from_utf8(&bytes[i..j]).map_err(|_| Error::BamParseCigar {
6✔
1862
                msg: format!("Invalid utf-8 bytes '{:?}'.", &bytes[i..j]),
×
1863
            })?;
1864
            let n = s.parse().map_err(|_| Error::BamParseCigar {
6✔
1865
                msg: format!("Unable to parse &str '{:?}' to u32.", s),
×
1866
            })?;
1867
            // get the operation
1868
            let op = &bytes[j];
6✔
1869
            inner.push(match op {
6✔
1870
                b'M' => Cigar::Match(n),
3✔
1871
                b'I' => Cigar::Ins(n),
1✔
1872
                b'D' => Cigar::Del(n),
1✔
1873
                b'N' => Cigar::RefSkip(n),
1✔
1874
                b'H' => {
1875
                    if i == 0 || j + 1 == text_len {
6✔
1876
                        Cigar::HardClip(n)
3✔
1877
                    } else {
1878
                        return Err(Error::BamParseCigar {
×
1879
                            msg: "Hard clipping ('H') is only valid at the start or end of a cigar."
×
1880
                                .to_owned(),
1881
                        });
1882
                    }
1883
                }
1884
                b'S' => {
1885
                    if i == 0
9✔
1886
                        || j + 1 == text_len
2✔
1887
                        || bytes[i-1] == b'H'
2✔
1888
                        || bytes[j+1..].iter().all(|c| c.is_ascii_digit() || *c == b'H') {
4✔
1889
                        Cigar::SoftClip(n)
1✔
1890
                    } else {
1891
                        return Err(Error::BamParseCigar {
×
1892
                        msg: "Soft clips ('S') can only have hard clips ('H') between them and the end of the CIGAR string."
×
1893
                            .to_owned(),
1894
                        });
1895
                    }
1896
                },
1897
                b'P' => Cigar::Pad(n),
1✔
1898
                b'=' => Cigar::Equal(n),
3✔
1899
                b'X' => Cigar::Diff(n),
3✔
1900
                op => {
×
1901
                    return Err(Error::BamParseCigar {
×
1902
                        msg: format!("Expected cigar operation [MIDNSHP=X] but got [{}]", op),
×
1903
                    })
1904
                }
1905
            });
1906
            i = j + 1;
3✔
1907
        }
1908
        Ok(CigarString(inner))
3✔
1909
    }
1910
}
1911

1912
impl TryFrom<&str> for CigarString {
1913
    type Error = Error;
1914

1915
    /// Create a CigarString from given &str.
1916
    /// # Example
1917
    /// ```
1918
    /// use rust_htslib::bam::record::*;
1919
    /// use rust_htslib::bam::record::CigarString;
1920
    /// use rust_htslib::bam::record::Cigar::*;
1921
    /// use std::convert::TryFrom;
1922
    ///
1923
    /// let cigar_str = "2H10M5X3=2H";
1924
    /// let cigar = CigarString::try_from(cigar_str)
1925
    ///     .expect("Unable to parse cigar string.");
1926
    /// let expected_cigar = CigarString(vec![
1927
    ///     HardClip(2),
1928
    ///     Match(10),
1929
    ///     Diff(5),
1930
    ///     Equal(3),
1931
    ///     HardClip(2),
1932
    /// ]);
1933
    /// assert_eq!(cigar, expected_cigar);
1934
    /// ```
1935
    fn try_from(text: &str) -> Result<Self> {
2✔
1936
        let bytes = text.as_bytes();
2✔
1937
        if text.chars().count() != bytes.len() {
2✔
1938
            return Err(Error::BamParseCigar {
1✔
1939
                msg: "CIGAR string contained non-ASCII characters, which are not valid. Valid are [0-9MIDNSHP=X].".to_owned(),
1✔
1940
            });
1941
        }
1942
        CigarString::try_from(bytes)
2✔
1943
    }
1944
}
1945

1946
impl<'a> CigarString {
1947
    pub fn iter(&'a self) -> ::std::slice::Iter<'a, Cigar> {
1✔
1948
        self.into_iter()
1✔
1949
    }
1950
}
1951

1952
impl<'a> IntoIterator for &'a CigarString {
1953
    type Item = &'a Cigar;
1954
    type IntoIter = ::std::slice::Iter<'a, Cigar>;
1955

1956
    fn into_iter(self) -> Self::IntoIter {
2✔
1957
        self.0.iter()
2✔
1958
    }
1959
}
1960

1961
impl fmt::Display for CigarString {
1962
    fn fmt(&self, fmt: &mut fmt::Formatter<'_>) -> Result<(), fmt::Error> {
2✔
1963
        for op in self {
4✔
1964
            fmt.write_fmt(format_args!("{}{}", op.len(), op.char()))?;
8✔
1965
        }
1966
        Ok(())
2✔
1967
    }
1968
}
1969

1970
// Get number of leading/trailing softclips if a CigarString taking hardclips into account
1971
fn calc_softclips<'a>(mut cigar: impl DoubleEndedIterator<Item = &'a Cigar>) -> i64 {
2✔
1972
    match (cigar.next(), cigar.next()) {
8✔
1973
        (Some(Cigar::HardClip(_)), Some(Cigar::SoftClip(s))) | (Some(Cigar::SoftClip(s)), _) => {
4✔
1974
            *s as i64
2✔
1975
        }
1976
        _ => 0,
×
1977
    }
1978
}
1979

1980
#[derive(Eq, PartialEq, Clone, Debug)]
1981
pub struct CigarStringView {
1982
    inner: CigarString,
1983
    pos: i64,
1984
}
1985

1986
impl CigarStringView {
1987
    /// Construct a new CigarStringView from a CigarString at a position
1988
    pub fn new(c: CigarString, pos: i64) -> CigarStringView {
1✔
1989
        CigarStringView { inner: c, pos }
1990
    }
1991

1992
    /// Get (exclusive) end position of alignment.
1993
    pub fn end_pos(&self) -> i64 {
1✔
1994
        let mut pos = self.pos;
1✔
1995
        for c in self {
2✔
1996
            match c {
3✔
1997
                Cigar::Match(l)
1✔
1998
                | Cigar::RefSkip(l)
1999
                | Cigar::Del(l)
2000
                | Cigar::Equal(l)
2001
                | Cigar::Diff(l) => pos += *l as i64,
2002
                // these don't add to end_pos on reference
2003
                Cigar::Ins(_) | Cigar::SoftClip(_) | Cigar::HardClip(_) | Cigar::Pad(_) => (),
2004
            }
2005
        }
2006
        pos
1✔
2007
    }
2008

2009
    /// Get the start position of the alignment (0-based).
2010
    pub fn pos(&self) -> i64 {
1✔
2011
        self.pos
1✔
2012
    }
2013

2014
    /// Get number of bases softclipped at the beginning of the alignment.
2015
    pub fn leading_softclips(&self) -> i64 {
1✔
2016
        calc_softclips(self.iter())
1✔
2017
    }
2018

2019
    /// Get number of bases softclipped at the end of the alignment.
2020
    pub fn trailing_softclips(&self) -> i64 {
1✔
2021
        calc_softclips(self.iter().rev())
1✔
2022
    }
2023

2024
    /// Get number of bases hardclipped at the beginning of the alignment.
2025
    pub fn leading_hardclips(&self) -> i64 {
×
2026
        self.first().map_or(0, |cigar| {
×
2027
            if let Cigar::HardClip(s) = cigar {
×
2028
                *s as i64
×
2029
            } else {
2030
                0
×
2031
            }
2032
        })
2033
    }
2034

2035
    /// Get number of bases hardclipped at the end of the alignment.
2036
    pub fn trailing_hardclips(&self) -> i64 {
×
2037
        self.last().map_or(0, |cigar| {
×
2038
            if let Cigar::HardClip(s) = cigar {
×
2039
                *s as i64
×
2040
            } else {
2041
                0
×
2042
            }
2043
        })
2044
    }
2045

2046
    /// For a given position in the reference, get corresponding position within read.
2047
    /// If reference position is outside of the read alignment, return None.
2048
    ///
2049
    /// # Arguments
2050
    ///
2051
    /// * `ref_pos` - the reference position
2052
    /// * `include_softclips` - if true, softclips will be considered as matches or mismatches
2053
    /// * `include_dels` - if true, positions within deletions will be considered (first reference matching read position after deletion will be returned)
2054
    ///
2055
    pub fn read_pos(
1✔
2056
        &self,
2057
        ref_pos: u32,
2058
        include_softclips: bool,
2059
        include_dels: bool,
2060
    ) -> Result<Option<u32>> {
2061
        let mut rpos = self.pos as u32; // reference position
1✔
2062
        let mut qpos = 0u32; // position within read
1✔
2063
        let mut j = 0; // index into cigar operation vector
1✔
2064

2065
        // find first cigar operation referring to qpos = 0 (and thus bases in record.seq()),
2066
        // because all augmentations of qpos and rpos before that are invalid
2067
        for (i, c) in self.iter().enumerate() {
2✔
2068
            match c {
1✔
2069
                Cigar::Match(_) |
2070
                Cigar::Diff(_)  |
2071
                Cigar::Equal(_) |
2072
                // this is unexpected, but bwa + GATK indel realignment can produce insertions
2073
                // before matching positions
2074
                Cigar::Ins(_) => {
2075
                    j = i;
1✔
2076
                    break;
1✔
2077
                },
2078
                Cigar::SoftClip(l) => {
1✔
2079
                    j = i;
1✔
2080
                    if include_softclips {
2✔
2081
                        // Alignment starts with softclip and we want to include it in the
2082
                        // projection of the reference position. However, the POS field does not
2083
                        // include the softclip. Hence we have to subtract its length.
2084
                        rpos = rpos.saturating_sub(*l);
2✔
2085
                    }
2086
                    break;
2087
                },
2088
                Cigar::Del(_) => {
2089
                    return Err(Error::BamUnexpectedCigarOperation {
×
2090
                        msg: "'deletion' (D) found before any operation describing read sequence".to_owned()
×
2091
                    });
2092
                },
2093
                Cigar::RefSkip(_) => {
2094
                    return Err(Error::BamUnexpectedCigarOperation {
×
2095
                        msg: "'reference skip' (N) found before any operation describing read sequence".to_owned()
×
2096
                    });
2097
                },
2098
                Cigar::HardClip(_) if i > 0 && i < self.len()-1 => {
2✔
2099
                    return Err(Error::BamUnexpectedCigarOperation{
1✔
2100
                        msg: "'hard clip' (H) found in between operations, contradicting SAMv1 spec that hard clips can only be at the ends of reads".to_owned()
1✔
2101
                    });
2102
                },
2103
                // if we have reached the end of the CigarString with only pads and hard clips, we have no read position matching the variant
2104
                Cigar::Pad(_) | Cigar::HardClip(_) if i == self.len()-1 => return Ok(None),
3✔
2105
                // skip leading HardClips and Pads, as they consume neither read sequence nor reference sequence
2106
                Cigar::Pad(_) | Cigar::HardClip(_) => ()
2107
            }
2108
        }
2109

2110
        let contains_ref_pos = |cigar_op_start: u32, cigar_op_length: u32| {
2✔
2111
            cigar_op_start <= ref_pos && cigar_op_start + cigar_op_length > ref_pos
2✔
2112
        };
2113

2114
        while rpos <= ref_pos && j < self.len() {
2✔
2115
            match self[j] {
9✔
2116
                // potential SNV evidence
2117
                Cigar::Match(l) | Cigar::Diff(l) | Cigar::Equal(l) if contains_ref_pos(rpos, l) => {
5✔
2118
                    // difference between desired position and first position of current cigar
2119
                    // operation
2120
                    qpos += ref_pos - rpos;
2✔
2121
                    return Ok(Some(qpos));
1✔
2122
                }
2123
                Cigar::SoftClip(l) if include_softclips && contains_ref_pos(rpos, l) => {
2✔
2124
                    qpos += ref_pos - rpos;
2✔
2125
                    return Ok(Some(qpos));
1✔
2126
                }
2127
                Cigar::Del(l) if include_dels && contains_ref_pos(rpos, l) => {
2✔
2128
                    // qpos shall resemble the start of the deletion
2129
                    return Ok(Some(qpos));
1✔
2130
                }
2131
                // for others, just increase pos and qpos as needed
2132
                Cigar::Match(l) | Cigar::Diff(l) | Cigar::Equal(l) => {
3✔
2133
                    rpos += l;
2✔
2134
                    qpos += l;
2✔
2135
                    j += 1;
2✔
2136
                }
2137
                Cigar::SoftClip(l) => {
1✔
2138
                    qpos += l;
2✔
2139
                    j += 1;
2✔
2140
                    if include_softclips {
2✔
2141
                        rpos += l;
1✔
2142
                    }
2143
                }
2144
                Cigar::Ins(l) => {
1✔
2145
                    qpos += l;
2✔
2146
                    j += 1;
2✔
2147
                }
2148
                Cigar::RefSkip(l) | Cigar::Del(l) => {
2✔
2149
                    rpos += l;
1✔
2150
                    j += 1;
2✔
2151
                }
2152
                Cigar::Pad(_) => {
1✔
2153
                    j += 1;
2✔
2154
                }
2155
                Cigar::HardClip(_) if j < self.len() - 1 => {
2✔
2156
                    return Err(Error::BamUnexpectedCigarOperation{
1✔
2157
                        msg: "'hard clip' (H) found in between operations, contradicting SAMv1 spec that hard clips can only be at the ends of reads".to_owned()
1✔
2158
                    });
2159
                }
2160
                Cigar::HardClip(_) => return Ok(None),
×
2161
            }
2162
        }
2163

2164
        Ok(None)
1✔
2165
    }
2166

2167
    /// transfer ownership of the Cigar out of the CigarView
2168
    pub fn take(self) -> CigarString {
1✔
2169
        self.inner
1✔
2170
    }
2171
}
2172

2173
impl ops::Deref for CigarStringView {
2174
    type Target = CigarString;
2175

2176
    fn deref(&self) -> &CigarString {
1✔
2177
        &self.inner
2178
    }
2179
}
2180

2181
impl ops::Index<usize> for CigarStringView {
2182
    type Output = Cigar;
2183

2184
    fn index(&self, index: usize) -> &Cigar {
1✔
2185
        self.inner.index(index)
1✔
2186
    }
2187
}
2188

2189
impl ops::IndexMut<usize> for CigarStringView {
2190
    fn index_mut(&mut self, index: usize) -> &mut Cigar {
×
2191
        self.inner.index_mut(index)
×
2192
    }
2193
}
2194

2195
impl<'a> CigarStringView {
2196
    pub fn iter(&'a self) -> ::std::slice::Iter<'a, Cigar> {
1✔
2197
        self.inner.into_iter()
1✔
2198
    }
2199
}
2200

2201
impl<'a> IntoIterator for &'a CigarStringView {
2202
    type Item = &'a Cigar;
2203
    type IntoIter = ::std::slice::Iter<'a, Cigar>;
2204

2205
    fn into_iter(self) -> Self::IntoIter {
1✔
2206
        self.inner.into_iter()
1✔
2207
    }
2208
}
2209

2210
impl fmt::Display for CigarStringView {
2211
    fn fmt(&self, fmt: &mut fmt::Formatter<'_>) -> Result<(), fmt::Error> {
×
2212
        self.inner.fmt(fmt)
×
2213
    }
2214
}
2215

2216
pub struct BaseModificationMetadata {
2217
    pub strand: i32,
2218
    pub implicit: i32,
2219
    pub canonical: i8,
2220
}
2221

2222
/// struct containing the internal state required to access
2223
/// the base modifications for a bam::Record
2224
pub struct BaseModificationState<'a> {
2225
    record: &'a Record,
2226
    state: *mut htslib::hts_base_mod_state,
2227
    buffer: Vec<htslib::hts_base_mod>,
2228
    buffer_pos: i32,
2229
}
2230

2231
impl BaseModificationState<'_> {
2232
    /// Initialize a new BaseModification struct from a bam::Record
2233
    /// This function allocates memory for the state structure
2234
    /// and initializes the iterator to the start of the modification
2235
    /// records.
2236
    fn new<'a>(r: &'a Record) -> Result<BaseModificationState<'a>> {
2✔
2237
        let mut bm = unsafe {
2238
            BaseModificationState {
2239
                record: r,
2240
                state: hts_sys::hts_base_mod_state_alloc(),
2✔
2241
                buffer: Vec::new(),
2✔
2242
                buffer_pos: -1,
×
2243
            }
2244
        };
2245

2246
        if bm.state.is_null() {
4✔
2247
            panic!("Unable to allocate memory for hts_base_mod_state");
×
2248
        }
2249

2250
        // parse the MM tag to initialize the state
2251
        unsafe {
2252
            let ret = hts_sys::bam_parse_basemod(bm.record.inner_ptr(), bm.state);
4✔
2253
            if ret != 0 {
2✔
2254
                return Err(Error::BamBaseModificationTagNotFound);
×
2255
            }
2256
        }
2257

2258
        let types = bm.recorded();
4✔
2259
        bm.buffer.reserve(types.len());
2✔
2260
        return Ok(bm);
2✔
2261
    }
2262

2263
    pub fn buffer_next_mods(&mut self) -> Result<usize> {
2✔
2264
        unsafe {
2265
            let ret = hts_sys::bam_next_basemod(
2266
                self.record.inner_ptr(),
2✔
2267
                self.state,
2✔
2268
                self.buffer.as_mut_ptr(),
2✔
2269
                self.buffer.capacity() as i32,
2✔
2270
                &mut self.buffer_pos,
2✔
2271
            );
2272

2273
            if ret < 0 {
2✔
2274
                return Err(Error::BamBaseModificationIterationFailed);
×
2275
            }
2276

2277
            // the htslib API won't write more than buffer.capacity() mods to the output array but it will
2278
            // return the actual number of modifications found. We return an error to the caller
2279
            // in the case where there was insufficient storage to return all mods.
2280
            if ret as usize > self.buffer.capacity() {
2✔
2281
                return Err(Error::BamBaseModificationTooManyMods);
×
2282
            }
2283

2284
            // we read the modifications directly into the vector, which does
2285
            // not update the length so needs to be manually set
2286
            self.buffer.set_len(ret as usize);
2✔
2287

2288
            return Ok(ret as usize);
2✔
2289
        }
2290
    }
2291

2292
    /// Return an array containing the modification codes listed for this record.
2293
    /// Positive values are ascii character codes (eg m), negative values are chEBI codes.
2294
    pub fn recorded<'a>(&self) -> &'a [i32] {
2✔
2295
        unsafe {
2296
            let mut n: i32 = 0;
2✔
2297
            let data_ptr: *const i32 = hts_sys::bam_mods_recorded(self.state, &mut n);
2✔
2298

2299
            // htslib should not return a null pointer, even when there are no base mods
2300
            if data_ptr.is_null() {
2✔
2301
                panic!("Unable to obtain pointer to base modifications");
×
2302
            }
2303
            assert!(n >= 0);
2✔
2304
            return slice::from_raw_parts(data_ptr, n as usize);
2✔
2305
        }
2306
    }
2307

2308
    /// Return metadata for the specified character code indicating the strand
2309
    /// the base modification was called on, whether the tag uses implicit mode
2310
    /// and the ascii code for the canonical base.
2311
    /// If there are multiple modifications with the same code this will return the data
2312
    /// for the first mod.  See https://github.com/samtools/htslib/issues/1635
2313
    pub fn query_type<'a>(&self, code: i32) -> Result<BaseModificationMetadata> {
2✔
2314
        unsafe {
2315
            let mut strand: i32 = 0;
2✔
2316
            let mut implicit: i32 = 0;
2✔
2317
            let mut canonical: i8 = 0;
2✔
2318

2319
            let ret = hts_sys::bam_mods_query_type(
2320
                self.state,
2✔
2321
                code,
×
2322
                &mut strand,
×
2323
                &mut implicit,
×
2324
                &mut canonical,
×
2325
            );
2326
            if ret == -1 {
2✔
2327
                return Err(Error::BamBaseModificationTypeNotFound);
×
2328
            } else {
2329
                return Ok(BaseModificationMetadata {
2✔
2330
                    strand,
2✔
2331
                    implicit,
2✔
2332
                    canonical,
2✔
2333
                });
2334
            }
2335
        }
2336
    }
2337
}
2338

2339
impl Drop for BaseModificationState<'_> {
2340
    fn drop<'a>(&mut self) {
2✔
2341
        unsafe {
2342
            hts_sys::hts_base_mod_state_free(self.state);
2✔
2343
        }
2344
    }
2345
}
2346

2347
/// Iterator over the base modifications that returns
2348
/// a vector for all of the mods at each position
2349
pub struct BaseModificationsPositionIter<'a> {
2350
    mod_state: BaseModificationState<'a>,
2351
}
2352

2353
impl BaseModificationsPositionIter<'_> {
2354
    fn new<'a>(r: &'a Record) -> Result<BaseModificationsPositionIter<'a>> {
1✔
2355
        let state = BaseModificationState::new(r)?;
1✔
2356
        Ok(BaseModificationsPositionIter { mod_state: state })
1✔
2357
    }
2358

2359
    pub fn recorded<'a>(&self) -> &'a [i32] {
×
2360
        return self.mod_state.recorded();
×
2361
    }
2362

2363
    pub fn query_type<'a>(&self, code: i32) -> Result<BaseModificationMetadata> {
×
2364
        return self.mod_state.query_type(code);
×
2365
    }
2366
}
2367

2368
impl<'a> Iterator for BaseModificationsPositionIter<'a> {
2369
    type Item = Result<(i32, Vec<hts_sys::hts_base_mod>)>;
2370

2371
    fn next(&mut self) -> Option<Self::Item> {
1✔
2372
        let ret = self.mod_state.buffer_next_mods();
1✔
2373

2374
        // Three possible things happened in buffer_next_mods:
2375
        // 1. the htslib API call was successful but there are no more mods
2376
        // 2. ths htslib API call was successful and we read some mods
2377
        // 3. the htslib API call failed, we propogate the error wrapped in an option
2378
        match ret {
1✔
2379
            Ok(num_mods) => {
1✔
2380
                if num_mods == 0 {
1✔
2381
                    return None;
1✔
2382
                } else {
2383
                    let data = (self.mod_state.buffer_pos, self.mod_state.buffer.clone());
1✔
2384
                    return Some(Ok(data));
1✔
2385
                }
2386
            }
2387
            Err(e) => return Some(Err(e)),
×
2388
        }
2389
    }
2390
}
2391

2392
/// Iterator over the base modifications that returns
2393
/// the next modification found, one by one
2394
pub struct BaseModificationsIter<'a> {
2395
    mod_state: BaseModificationState<'a>,
2396
    buffer_idx: usize,
2397
}
2398

2399
impl BaseModificationsIter<'_> {
2400
    fn new<'a>(r: &'a Record) -> Result<BaseModificationsIter<'a>> {
2✔
2401
        let state = BaseModificationState::new(r)?;
2✔
2402
        Ok(BaseModificationsIter {
2✔
2403
            mod_state: state,
×
2404
            buffer_idx: 0,
×
2405
        })
2406
    }
2407

2408
    pub fn recorded<'a>(&self) -> &'a [i32] {
2✔
2409
        return self.mod_state.recorded();
2✔
2410
    }
2411

2412
    pub fn query_type<'a>(&self, code: i32) -> Result<BaseModificationMetadata> {
2✔
2413
        return self.mod_state.query_type(code);
2✔
2414
    }
2415
}
2416

2417
impl<'a> Iterator for BaseModificationsIter<'a> {
2418
    type Item = Result<(i32, hts_sys::hts_base_mod)>;
2419

2420
    fn next(&mut self) -> Option<Self::Item> {
2✔
2421
        if self.buffer_idx == self.mod_state.buffer.len() {
4✔
2422
            // need to use the internal state to read the next
2423
            // set of modifications into the buffer
2424
            let ret = self.mod_state.buffer_next_mods();
2✔
2425

2426
            match ret {
2✔
2427
                Ok(num_mods) => {
2✔
2428
                    if num_mods == 0 {
2✔
2429
                        // done iterating
2430
                        return None;
2✔
2431
                    } else {
2432
                        // we read some mods, reset the position in the buffer then fall through
2433
                        self.buffer_idx = 0;
2✔
2434
                    }
2435
                }
2436
                Err(e) => return Some(Err(e)),
×
2437
            }
2438
        }
2439

2440
        // if we got here when there are mods buffered that we haven't emitted yet
2441
        assert!(self.buffer_idx < self.mod_state.buffer.len());
2✔
2442
        let data = (
2✔
2443
            self.mod_state.buffer_pos,
2✔
2444
            self.mod_state.buffer[self.buffer_idx],
2✔
2445
        );
2446
        self.buffer_idx += 1;
4✔
2447
        return Some(Ok(data));
2✔
2448
    }
2449
}
2450

2451
#[cfg(test)]
2452
mod tests {
2453
    use super::*;
2454

2455
    #[test]
2456
    fn test_cigar_string() {
2457
        let cigar = CigarString(vec![Cigar::Match(100), Cigar::SoftClip(10)]);
2458

2459
        assert_eq!(cigar[0], Cigar::Match(100));
2460
        assert_eq!(format!("{}", cigar), "100M10S");
2461
        for op in &cigar {
2462
            println!("{}", op);
2463
        }
2464
    }
2465

2466
    #[test]
2467
    fn test_cigar_string_view_pos() {
2468
        let cigar = CigarString(vec![Cigar::Match(100), Cigar::SoftClip(10)]).into_view(5);
2469
        assert_eq!(cigar.pos(), 5);
2470
    }
2471

2472
    #[test]
2473
    fn test_cigar_string_leading_softclips() {
2474
        let cigar = CigarString(vec![Cigar::SoftClip(10), Cigar::Match(100)]).into_view(0);
2475
        assert_eq!(cigar.leading_softclips(), 10);
2476
        let cigar2 = CigarString(vec![
2477
            Cigar::HardClip(5),
2478
            Cigar::SoftClip(10),
2479
            Cigar::Match(100),
2480
        ])
2481
        .into_view(0);
2482
        assert_eq!(cigar2.leading_softclips(), 10);
2483
    }
2484

2485
    #[test]
2486
    fn test_cigar_string_trailing_softclips() {
2487
        let cigar = CigarString(vec![Cigar::Match(100), Cigar::SoftClip(10)]).into_view(0);
2488
        assert_eq!(cigar.trailing_softclips(), 10);
2489
        let cigar2 = CigarString(vec![
2490
            Cigar::Match(100),
2491
            Cigar::SoftClip(10),
2492
            Cigar::HardClip(5),
2493
        ])
2494
        .into_view(0);
2495
        assert_eq!(cigar2.trailing_softclips(), 10);
2496
    }
2497

2498
    #[test]
2499
    fn test_cigar_read_pos() {
2500
        let vpos = 5; // variant position
2501

2502
        // Ignore leading HardClip
2503
        // ref:       00 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15
2504
        // var:                       V
2505
        // c01: 7H                 M  M
2506
        // qpos:                  00 01
2507
        let c01 = CigarString(vec![Cigar::HardClip(7), Cigar::Match(2)]).into_view(4);
2508
        assert_eq!(c01.read_pos(vpos, false, false).unwrap(), Some(1));
2509

2510
        // Skip leading SoftClip or use as pre-POS matches
2511
        // ref:       00 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15
2512
        // var:                       V
2513
        // c02: 5H2S         M  M  M  M  M  M
2514
        // qpos:  00        02 03 04 05 06 07
2515
        // c02: 5H     S  S  M  M  M  M  M  M
2516
        // qpos:      00 01 02 03 04 05 06 07
2517
        let c02 = CigarString(vec![Cigar::SoftClip(2), Cigar::Match(6)]).into_view(2);
2518
        assert_eq!(c02.read_pos(vpos, false, false).unwrap(), Some(5));
2519
        assert_eq!(c02.read_pos(vpos, true, false).unwrap(), Some(5));
2520

2521
        // Skip leading SoftClip returning None for unmatched reference positiong or use as
2522
        // pre-POS matches
2523
        // ref:       00 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15
2524
        // var:                       V
2525
        // c03:  3S                      M  M
2526
        // qpos: 00                     03 04
2527
        // c03:                 S  S  S  M  M
2528
        // qpos:               00 01 02 03 04
2529
        let c03 = CigarString(vec![Cigar::SoftClip(3), Cigar::Match(6)]).into_view(6);
2530
        assert_eq!(c03.read_pos(vpos, false, false).unwrap(), None);
2531
        assert_eq!(c03.read_pos(vpos, true, false).unwrap(), Some(2));
2532

2533
        // Skip leading Insertion before variant position
2534
        // ref:       00 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15
2535
        // var:                       V
2536
        // c04:  3I                X  X  X
2537
        // qpos: 00               03 04 05
2538
        let c04 = CigarString(vec![Cigar::Ins(3), Cigar::Diff(3)]).into_view(4);
2539
        assert_eq!(c04.read_pos(vpos, true, false).unwrap(), Some(4));
2540

2541
        // Matches and deletion before variant position
2542
        // ref:       00 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15
2543
        // var:                       V
2544
        // c05:        =  =  D  D  X  =  =
2545
        // qpos:      00 01       02 03 04 05
2546
        let c05 = CigarString(vec![
2547
            Cigar::Equal(2),
2548
            Cigar::Del(2),
2549
            Cigar::Diff(1),
2550
            Cigar::Equal(2),
2551
        ])
2552
        .into_view(0);
2553
        assert_eq!(c05.read_pos(vpos, true, false).unwrap(), Some(3));
2554

2555
        // single nucleotide Deletion covering variant position
2556
        // ref:       00 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15
2557
        // var:                       V
2558
        // c06:                 =  =  D  X  X
2559
        // qpos:               00 01    02 03
2560
        let c06 = CigarString(vec![Cigar::Equal(2), Cigar::Del(1), Cigar::Diff(2)]).into_view(3);
2561
        assert_eq!(c06.read_pos(vpos, false, true).unwrap(), Some(2));
2562
        assert_eq!(c06.read_pos(vpos, false, false).unwrap(), None);
2563

2564
        // three nucleotide Deletion covering variant position
2565
        // ref:       00 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15
2566
        // var:                       V
2567
        // c07:              =  =  D  D  D  M  M
2568
        // qpos:            00 01          02 03
2569
        let c07 = CigarString(vec![Cigar::Equal(2), Cigar::Del(3), Cigar::Match(2)]).into_view(2);
2570
        assert_eq!(c07.read_pos(vpos, false, true).unwrap(), Some(2));
2571
        assert_eq!(c07.read_pos(vpos, false, false).unwrap(), None);
2572

2573
        // three nucleotide RefSkip covering variant position
2574
        // ref:       00 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15
2575
        // var:                       V
2576
        // c08:              =  X  N  N  N  M  M
2577
        // qpos:            00 01          02 03
2578
        let c08 = CigarString(vec![
2579
            Cigar::Equal(1),
2580
            Cigar::Diff(1),
2581
            Cigar::RefSkip(3),
2582
            Cigar::Match(2),
2583
        ])
2584
        .into_view(2);
2585
        assert_eq!(c08.read_pos(vpos, false, true).unwrap(), None);
2586
        assert_eq!(c08.read_pos(vpos, false, false).unwrap(), None);
2587

2588
        // internal hard clip before variant pos
2589
        // ref:       00 01 02 03    04 05 06 07 08 09 10 11 12 13 14 15
2590
        // var:                          V
2591
        // c09: 3H           =  = 3H  =  =
2592
        // qpos:            00 01    02 03
2593
        let c09 = CigarString(vec![
2594
            Cigar::HardClip(3),
2595
            Cigar::Equal(2),
2596
            Cigar::HardClip(3),
2597
            Cigar::Equal(2),
2598
        ])
2599
        .into_view(2);
2600
        assert_eq!(c09.read_pos(vpos, false, true).is_err(), true);
2601

2602
        // Deletion right before variant position
2603
        // ref:       00 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15
2604
        // var:                       V
2605
        // c10:           M  M  D  D  M  M
2606
        // qpos:         00 01       02 03
2607
        let c10 = CigarString(vec![Cigar::Match(2), Cigar::Del(2), Cigar::Match(2)]).into_view(1);
2608
        assert_eq!(c10.read_pos(vpos, false, false).unwrap(), Some(2));
2609

2610
        // Insertion right before variant position
2611
        // ref:       00 01 02 03 04    05 06 07 08 09 10 11 12 13 14 15
2612
        // var:                          V
2613
        // c11:                 M  M 3I  M
2614
        // qpos:               00 01 02 05 06
2615
        let c11 = CigarString(vec![Cigar::Match(2), Cigar::Ins(3), Cigar::Match(2)]).into_view(3);
2616
        assert_eq!(c11.read_pos(vpos, false, false).unwrap(), Some(5));
2617

2618
        // Insertion right after variant position
2619
        // ref:       00 01 02 03 04 05    06 07 08 09 10 11 12 13 14 15
2620
        // var:                       V
2621
        // c12:                 M  M  M 2I  =
2622
        // qpos:               00 01 02 03 05
2623
        let c12 = CigarString(vec![Cigar::Match(3), Cigar::Ins(2), Cigar::Equal(1)]).into_view(3);
2624
        assert_eq!(c12.read_pos(vpos, false, false).unwrap(), Some(2));
2625

2626
        // Deletion right after variant position
2627
        // ref:       00 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15
2628
        // var:                       V
2629
        // c13:                 M  M  M  D  =
2630
        // qpos:               00 01 02    03
2631
        let c13 = CigarString(vec![Cigar::Match(3), Cigar::Del(1), Cigar::Equal(1)]).into_view(3);
2632
        assert_eq!(c13.read_pos(vpos, false, false).unwrap(), Some(2));
2633

2634
        // A messy and complicated example, including a Pad operation
2635
        let vpos2 = 15;
2636
        // ref:       00    01 02    03 04 05    06 07 08 09 10 11 12 13 14 15
2637
        // var:                                                           V
2638
        // c14: 5H3S   = 2P  M  X 3I  M  M  D 2I  =  =  N  N  N  M  M  M  =  =  5S2H
2639
        // qpos:  00  03    04 05 06 09 10    11 13 14          15 16 17 18 19
2640
        let c14 = CigarString(vec![
2641
            Cigar::HardClip(5),
2642
            Cigar::SoftClip(3),
2643
            Cigar::Equal(1),
2644
            Cigar::Pad(2),
2645
            Cigar::Match(1),
2646
            Cigar::Diff(1),
2647
            Cigar::Ins(3),
2648
            Cigar::Match(2),
2649
            Cigar::Del(1),
2650
            Cigar::Ins(2),
2651
            Cigar::Equal(2),
2652
            Cigar::RefSkip(3),
2653
            Cigar::Match(3),
2654
            Cigar::Equal(2),
2655
            Cigar::SoftClip(5),
2656
            Cigar::HardClip(2),
2657
        ])
2658
        .into_view(0);
2659
        assert_eq!(c14.read_pos(vpos2, false, false).unwrap(), Some(19));
2660

2661
        // HardClip after Pad
2662
        // ref:       00 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15
2663
        // var:                       V
2664
        // c15: 5P1H            =  =  =
2665
        // qpos:               00 01 02
2666
        let c15 =
2667
            CigarString(vec![Cigar::Pad(5), Cigar::HardClip(1), Cigar::Equal(3)]).into_view(3);
2668
        assert_eq!(c15.read_pos(vpos, false, false).is_err(), true);
2669

2670
        // only HardClip and Pad operations
2671
        // c16: 7H5P2H
2672
        let c16 =
2673
            CigarString(vec![Cigar::HardClip(7), Cigar::Pad(5), Cigar::HardClip(2)]).into_view(3);
2674
        assert_eq!(c16.read_pos(vpos, false, false).unwrap(), None);
2675
    }
2676

2677
    #[test]
2678
    fn test_clone() {
2679
        let mut rec = Record::new();
2680
        rec.set_pos(300);
2681
        rec.set_qname(b"read1");
2682
        let clone = rec.clone();
2683
        assert_eq!(rec, clone);
2684
    }
2685

2686
    #[test]
2687
    fn test_flags() {
2688
        let mut rec = Record::new();
2689

2690
        rec.set_paired();
2691
        assert_eq!(rec.is_paired(), true);
2692

2693
        rec.set_supplementary();
2694
        assert_eq!(rec.is_supplementary(), true);
2695
        assert_eq!(rec.is_supplementary(), true);
2696

2697
        rec.unset_paired();
2698
        assert_eq!(rec.is_paired(), false);
2699
        assert_eq!(rec.is_supplementary(), true);
2700

2701
        rec.unset_supplementary();
2702
        assert_eq!(rec.is_paired(), false);
2703
        assert_eq!(rec.is_supplementary(), false);
2704
    }
2705

2706
    #[test]
2707
    fn test_cigar_parse() {
2708
        let cigar = "1S20M1D2I3X1=2H";
2709
        let parsed = CigarString::try_from(cigar).unwrap();
2710
        assert_eq!(parsed.to_string(), cigar);
2711
    }
2712
}
2713

2714
#[cfg(test)]
2715
mod alignment_cigar_tests {
2716
    use super::*;
2717
    use crate::bam::{Read, Reader};
2718
    use bio_types::alignment::AlignmentOperation::{Del, Ins, Match, Subst, Xclip, Yclip};
2719
    use bio_types::alignment::{Alignment, AlignmentMode};
2720

2721
    #[test]
2722
    fn test_cigar() {
2723
        let alignment = Alignment {
2724
            score: 5,
2725
            xstart: 3,
2726
            ystart: 0,
2727
            xend: 9,
2728
            yend: 10,
2729
            ylen: 10,
2730
            xlen: 10,
2731
            operations: vec![Match, Match, Match, Subst, Ins, Ins, Del, Del],
2732
            mode: AlignmentMode::Semiglobal,
2733
        };
2734
        assert_eq!(alignment.cigar(false), "3S3=1X2I2D1S");
2735
        assert_eq!(
2736
            CigarString::from_alignment(&alignment, false).0,
2737
            vec![
2738
                Cigar::SoftClip(3),
2739
                Cigar::Equal(3),
2740
                Cigar::Diff(1),
2741
                Cigar::Ins(2),
2742
                Cigar::Del(2),
2743
                Cigar::SoftClip(1),
2744
            ]
2745
        );
2746

2747
        let alignment = Alignment {
2748
            score: 5,
2749
            xstart: 0,
2750
            ystart: 5,
2751
            xend: 4,
2752
            yend: 10,
2753
            ylen: 10,
2754
            xlen: 5,
2755
            operations: vec![Yclip(5), Match, Subst, Subst, Ins, Del, Del, Xclip(1)],
2756
            mode: AlignmentMode::Custom,
2757
        };
2758
        assert_eq!(alignment.cigar(false), "1=2X1I2D1S");
2759
        assert_eq!(alignment.cigar(true), "1=2X1I2D1H");
2760
        assert_eq!(
2761
            CigarString::from_alignment(&alignment, false).0,
2762
            vec![
2763
                Cigar::Equal(1),
2764
                Cigar::Diff(2),
2765
                Cigar::Ins(1),
2766
                Cigar::Del(2),
2767
                Cigar::SoftClip(1),
2768
            ]
2769
        );
2770
        assert_eq!(
2771
            CigarString::from_alignment(&alignment, true).0,
2772
            vec![
2773
                Cigar::Equal(1),
2774
                Cigar::Diff(2),
2775
                Cigar::Ins(1),
2776
                Cigar::Del(2),
2777
                Cigar::HardClip(1),
2778
            ]
2779
        );
2780

2781
        let alignment = Alignment {
2782
            score: 5,
2783
            xstart: 0,
2784
            ystart: 5,
2785
            xend: 3,
2786
            yend: 8,
2787
            ylen: 10,
2788
            xlen: 3,
2789
            operations: vec![Yclip(5), Subst, Match, Subst, Yclip(2)],
2790
            mode: AlignmentMode::Custom,
2791
        };
2792
        assert_eq!(alignment.cigar(false), "1X1=1X");
2793
        assert_eq!(
2794
            CigarString::from_alignment(&alignment, false).0,
2795
            vec![Cigar::Diff(1), Cigar::Equal(1), Cigar::Diff(1)]
2796
        );
2797

2798
        let alignment = Alignment {
2799
            score: 5,
2800
            xstart: 0,
2801
            ystart: 5,
2802
            xend: 3,
2803
            yend: 8,
2804
            ylen: 10,
2805
            xlen: 3,
2806
            operations: vec![Subst, Match, Subst],
2807
            mode: AlignmentMode::Semiglobal,
2808
        };
2809
        assert_eq!(alignment.cigar(false), "1X1=1X");
2810
        assert_eq!(
2811
            CigarString::from_alignment(&alignment, false).0,
2812
            vec![Cigar::Diff(1), Cigar::Equal(1), Cigar::Diff(1)]
2813
        );
2814
    }
2815

2816
    #[test]
2817
    fn test_read_orientation_f1r2() {
2818
        let mut bam = Reader::from_path(&"test/test_paired.sam").unwrap();
2819

2820
        for res in bam.records() {
2821
            let record = res.unwrap();
2822
            assert_eq!(
2823
                record.read_pair_orientation(),
2824
                SequenceReadPairOrientation::F1R2
2825
            );
2826
        }
2827
    }
2828

2829
    #[test]
2830
    pub fn test_cigar_parsing_non_ascii_error() {
2831
        let cigar_str = "43ጷ";
2832
        let expected_error = Err(Error::BamParseCigar {
2833
                msg: "CIGAR string contained non-ASCII characters, which are not valid. Valid are [0-9MIDNSHP=X].".to_owned(),
2834
            });
2835

2836
        let result = CigarString::try_from(cigar_str);
2837
        assert_eq!(expected_error, result);
2838
    }
2839

2840
    #[test]
2841
    pub fn test_cigar_parsing() {
2842
        // parsing test cases
2843
        let cigar_strs = vec![
2844
            "1H10M4D100I300N1102=10P25X11S", // test every cigar opt
2845
            "100M",                          // test a single op
2846
            "",                              // test empty input
2847
            "1H1=1H",                        // test simple hardclip
2848
            "1S1=1S",                        // test simple softclip
2849
            "11H11S11=11S11H",               // test complex softclip
2850
            "10H",
2851
            "10S",
2852
        ];
2853
        // expected results
2854
        let cigars = vec![
2855
            CigarString(vec![
2856
                Cigar::HardClip(1),
2857
                Cigar::Match(10),
2858
                Cigar::Del(4),
2859
                Cigar::Ins(100),
2860
                Cigar::RefSkip(300),
2861
                Cigar::Equal(1102),
2862
                Cigar::Pad(10),
2863
                Cigar::Diff(25),
2864
                Cigar::SoftClip(11),
2865
            ]),
2866
            CigarString(vec![Cigar::Match(100)]),
2867
            CigarString(vec![]),
2868
            CigarString(vec![
2869
                Cigar::HardClip(1),
2870
                Cigar::Equal(1),
2871
                Cigar::HardClip(1),
2872
            ]),
2873
            CigarString(vec![
2874
                Cigar::SoftClip(1),
2875
                Cigar::Equal(1),
2876
                Cigar::SoftClip(1),
2877
            ]),
2878
            CigarString(vec![
2879
                Cigar::HardClip(11),
2880
                Cigar::SoftClip(11),
2881
                Cigar::Equal(11),
2882
                Cigar::SoftClip(11),
2883
                Cigar::HardClip(11),
2884
            ]),
2885
            CigarString(vec![Cigar::HardClip(10)]),
2886
            CigarString(vec![Cigar::SoftClip(10)]),
2887
        ];
2888
        // compare
2889
        for (&cigar_str, truth) in cigar_strs.iter().zip(cigars.iter()) {
2890
            let cigar_parse = CigarString::try_from(cigar_str)
2891
                .expect(&format!("Unable to parse cigar: {}", cigar_str));
2892
            assert_eq!(&cigar_parse, truth);
2893
        }
2894
    }
2895
}
2896

2897
#[cfg(test)]
2898
mod basemod_tests {
2899
    use crate::bam::{Read, Reader};
2900

2901
    #[test]
2902
    pub fn test_count_recorded() {
2903
        let mut bam = Reader::from_path(&"test/base_mods/MM-double.sam").unwrap();
2904

2905
        for r in bam.records() {
2906
            let record = r.unwrap();
2907
            if let Ok(mods) = record.basemods_iter() {
2908
                let n = mods.recorded().len();
2909
                assert_eq!(n, 3);
2910
            };
2911
        }
2912
    }
2913

2914
    #[test]
2915
    pub fn test_query_type() {
2916
        let mut bam = Reader::from_path(&"test/base_mods/MM-orient.sam").unwrap();
2917

2918
        let mut n_fwd = 0;
2919
        let mut n_rev = 0;
2920

2921
        for r in bam.records() {
2922
            let record = r.unwrap();
2923
            if let Ok(mods) = record.basemods_iter() {
2924
                for mod_code in mods.recorded() {
2925
                    if let Ok(mod_metadata) = mods.query_type(*mod_code) {
2926
                        if mod_metadata.strand == 0 {
2927
                            n_fwd += 1;
2928
                        }
2929
                        if mod_metadata.strand == 1 {
2930
                            n_rev += 1;
2931
                        }
2932
                    }
2933
                }
2934
            };
2935
        }
2936
        assert_eq!(n_fwd, 2);
2937
        assert_eq!(n_rev, 2);
2938
    }
2939

2940
    #[test]
2941
    pub fn test_mod_iter() {
2942
        let mut bam = Reader::from_path(&"test/base_mods/MM-double.sam").unwrap();
2943
        let expected_positions = [1, 7, 12, 13, 13, 22, 30, 31];
2944
        let mut i = 0;
2945

2946
        for r in bam.records() {
2947
            let record = r.unwrap();
2948
            for res in record.basemods_iter().unwrap() {
2949
                if let Ok((position, _m)) = res {
2950
                    assert_eq!(position, expected_positions[i]);
2951
                    i += 1;
2952
                }
2953
            }
2954
        }
2955
    }
2956

2957
    #[test]
2958
    pub fn test_position_iter() {
2959
        let mut bam = Reader::from_path(&"test/base_mods/MM-double.sam").unwrap();
2960
        let expected_positions = [1, 7, 12, 13, 22, 30, 31];
2961
        let expected_counts = [1, 1, 1, 2, 1, 1, 1];
2962
        let mut i = 0;
2963

2964
        for r in bam.records() {
2965
            let record = r.unwrap();
2966
            for res in record.basemods_position_iter().unwrap() {
2967
                if let Ok((position, elements)) = res {
2968
                    assert_eq!(position, expected_positions[i]);
2969
                    assert_eq!(elements.len(), expected_counts[i]);
2970
                    i += 1;
2971
                }
2972
            }
2973
        }
2974
    }
2975
}
STATUS · Troubleshooting · Open an Issue · Sales · Support · CAREERS · ENTERPRISE · START FREE · SCHEDULE DEMO
ANNOUNCEMENTS · TWITTER · TOS & SLA · Supported CI Services · What's a CI service? · Automated Testing

© 2025 Coveralls, Inc