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

rust-bio / rust-htslib / 19626566167

24 Nov 2025 07:26AM UTC coverage: 81.906% (-0.03%) from 81.935%
19626566167

Pull #494

github

web-flow
Merge d2fab2115 into 8f1cdd75c
Pull Request #494: build(deps): bump actions/checkout from 5 to 6

2784 of 3399 relevant lines covered (81.91%)

27217.84 hits per line

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

84.04
/src/bam/mod.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.
5

6
//! Module for working with SAM, BAM, and CRAM files.
7

8
pub mod buffer;
9
pub mod ext;
10
pub mod header;
11
pub mod index;
12
pub mod pileup;
13
pub mod record;
14

15
#[cfg(feature = "serde_feature")]
16
pub mod record_serde;
17

18
use std::ffi;
19
use std::os::raw::c_char;
20
use std::path::Path;
21
use std::rc::Rc;
22
use std::slice;
23
use std::str;
24

25
use url::Url;
26

27
use crate::errors::{Error, Result};
28
use crate::htslib;
29
use crate::tpool::ThreadPool;
30
use crate::utils::path_as_bytes;
31

32
pub use crate::bam::buffer::RecordBuffer;
33
pub use crate::bam::header::Header;
34
pub use crate::bam::record::Record;
35
use hts_sys::{hts_fmt_option, sam_fields};
36
use std::convert::{TryFrom, TryInto};
37
use std::mem::MaybeUninit;
38

39
/// # Safety
40
///
41
/// Implementation for `Read::set_threads` and `Writer::set_threads`.
42
unsafe fn set_threads(htsfile: *mut htslib::htsFile, n_threads: usize) -> Result<()> {
2✔
43
    assert!(n_threads != 0, "n_threads must be > 0");
3✔
44

45
    if htslib::hts_set_threads(htsfile, n_threads as i32) != 0 {
4✔
46
        Err(Error::SetThreads)
×
47
    } else {
48
        Ok(())
2✔
49
    }
50
}
51

52
unsafe fn set_thread_pool(htsfile: *mut htslib::htsFile, tpool: &ThreadPool) -> Result<()> {
5✔
53
    let mut b = tpool.handle.borrow_mut();
9✔
54

55
    if htslib::hts_set_thread_pool(htsfile, &mut b.inner as *mut _) != 0 {
15✔
56
        Err(Error::ThreadPool)
×
57
    } else {
58
        Ok(())
5✔
59
    }
60
}
61

62
/// # Safety
63
///
64
/// Set the reference FAI index path in a `htslib::htsFile` struct for reading CRAM format.
65
pub unsafe fn set_fai_filename<P: AsRef<Path>>(
6✔
66
    htsfile: *mut htslib::htsFile,
67
    fasta_path: P,
68
) -> Result<()> {
69
    let path = if let Some(ext) = fasta_path.as_ref().extension() {
17✔
70
        fasta_path
6✔
71
            .as_ref()
72
            .with_extension(format!("{}.fai", ext.to_str().unwrap()))
26✔
73
    } else {
74
        fasta_path.as_ref().with_extension(".fai")
×
75
    };
76
    let p: &Path = path.as_ref();
17✔
77
    let c_str = ffi::CString::new(p.to_str().unwrap().as_bytes()).unwrap();
31✔
78
    if htslib::hts_set_fai_filename(htsfile, c_str.as_ptr()) == 0 {
18✔
79
        Ok(())
6✔
80
    } else {
81
        Err(Error::BamInvalidReferencePath { path: p.to_owned() })
×
82
    }
83
}
84

85
/// A trait for a BAM reader with a read method.
86
pub trait Read: Sized {
87
    /// Read next BAM record into given record.
88
    /// Use this method in combination with a single allocated record to avoid the reallocations
89
    /// occurring with the iterator.
90
    ///
91
    /// # Arguments
92
    ///
93
    /// * `record` - the record to be filled
94
    ///
95
    /// # Returns
96
    ///
97
    /// Some(Ok(())) if the record was read and None if no more records to read
98
    ///
99
    /// Example:
100
    /// ```
101
    /// use rust_htslib::errors::Error;
102
    /// use rust_htslib::bam::{Read, IndexedReader, Record};
103
    ///
104
    /// let mut bam = IndexedReader::from_path("test/test.bam").unwrap();
105
    /// bam.fetch((0, 1000, 2000)); // reads on tid 0, from 1000bp to 2000bp
106
    /// let mut record = Record::new();
107
    /// while let Some(result) = bam.read(&mut record) {
108
    ///     match result {
109
    ///         Ok(_) => {
110
    ///             println!("Read sequence: {:?}", record.seq().as_bytes());
111
    ///         }
112
    ///         Err(_) => panic!("BAM parsing failed...")
113
    ///     }
114
    /// }
115
    /// ```
116
    ///
117
    /// Consider using [`rc_records`](#tymethod.rc_records) instead.
118
    fn read(&mut self, record: &mut record::Record) -> Option<Result<()>>;
119

120
    /// Iterator over the records of the seeked region.
121
    /// Note that, while being convenient, this is less efficient than pre-allocating a
122
    /// `Record` and reading into it with the `read` method, since every iteration involves
123
    /// the allocation of a new `Record`.
124
    ///
125
    /// This is about 10% slower than record in micro benchmarks.
126
    ///
127
    /// Consider using [`rc_records`](#tymethod.rc_records) instead.
128
    fn records(&mut self) -> Records<'_, Self>;
129

130
    /// Records iterator using an Rc to avoid allocating a Record each turn.
131
    /// This is about 1% slower than the [`read`](#tymethod.read) based API in micro benchmarks,
132
    /// but has nicer ergonomics (and might not actually be slower in your applications).
133
    ///
134
    /// Example:
135
    /// ```
136
    /// use rust_htslib::errors::Error;
137
    /// use rust_htslib::bam::{Read, Reader, Record};
138
    /// use rust_htslib::htslib; // for BAM_F*
139
    /// let mut bam = Reader::from_path("test/test.bam").unwrap();
140
    ///
141
    /// for read in
142
    ///     bam.rc_records()
143
    ///     .map(|x| x.expect("Failure parsing Bam file"))
144
    ///     .filter(|read|
145
    ///         read.flags()
146
    ///          & (htslib::BAM_FUNMAP
147
    ///              | htslib::BAM_FSECONDARY
148
    ///              | htslib::BAM_FQCFAIL
149
    ///              | htslib::BAM_FDUP) as u16
150
    ///          == 0
151
    ///     )
152
    ///     .filter(|read| !read.is_reverse()) {
153
    ///     println!("Found a forward read: {:?}", read.qname());
154
    /// }
155
    ///
156
    /// //or to add the read qnames into a Vec
157
    /// let collected: Vec<_> = bam.rc_records().map(|read| read.unwrap().qname().to_vec()).collect();
158
    ///
159
    ///
160
    /// ```
161
    fn rc_records(&mut self) -> RcRecords<'_, Self>;
162

163
    /// Iterator over pileups.
164
    fn pileup(&mut self) -> pileup::Pileups<'_, Self>;
165

166
    /// Return the htsFile struct
167
    fn htsfile(&self) -> *mut htslib::htsFile;
168

169
    /// Return the header.
170
    fn header(&self) -> &HeaderView;
171

172
    /// Seek to the given virtual offset in the file
173
    fn seek(&mut self, offset: i64) -> Result<()> {
7✔
174
        let htsfile = unsafe { self.htsfile().as_ref() }.expect("bug: null pointer to htsFile");
37✔
175
        let ret = match htsfile.format.format {
13✔
176
            htslib::htsExactFormat_cram => unsafe {
177
                i64::from(htslib::cram_seek(
×
178
                    htsfile.fp.cram,
×
179
                    offset as libc::off_t,
×
180
                    libc::SEEK_SET,
×
181
                ))
182
            },
183
            _ => unsafe { htslib::bgzf_seek(htsfile.fp.bgzf, offset, libc::SEEK_SET) },
19✔
184
        };
185

186
        if ret == 0 {
8✔
187
            Ok(())
7✔
188
        } else {
189
            Err(Error::FileSeek)
×
190
        }
191
    }
192

193
    /// Report the current virtual offset
194
    fn tell(&self) -> i64 {
8✔
195
        // this reimplements the bgzf_tell macro
196
        let htsfile = unsafe { self.htsfile().as_ref() }.expect("bug: null pointer to htsFile");
43✔
197
        let bgzf = unsafe { *htsfile.fp.bgzf };
15✔
198
        (bgzf.block_address << 16) | (i64::from(bgzf.block_offset) & 0xFFFF)
15✔
199
    }
200

201
    /// Activate multi-threaded BAM read support in htslib. This should permit faster
202
    /// reading of large BAM files.
203
    ///
204
    /// Setting `nthreads` to `0` does not change the current state.  Note that it is not
205
    /// possible to set the number of background threads below `1` once it has been set.
206
    ///
207
    /// # Arguments
208
    ///
209
    /// * `n_threads` - number of extra background writer threads to use, must be `> 0`.
210
    fn set_threads(&mut self, n_threads: usize) -> Result<()> {
×
211
        unsafe { set_threads(self.htsfile(), n_threads) }
×
212
    }
213

214
    /// Use a shared thread-pool for writing. This permits controlling the total
215
    /// thread count when multiple readers and writers are working simultaneously.
216
    /// A thread pool can be created with `crate::tpool::ThreadPool::new(n_threads)`
217
    ///
218
    /// # Arguments
219
    ///
220
    /// * `tpool` - thread pool to use for compression work.
221
    fn set_thread_pool(&mut self, tpool: &ThreadPool) -> Result<()>;
222

223
    /// If the underlying file is in CRAM format, allows modifying CRAM options.
224
    /// Note that this method does *not* check that the underlying file actually is in CRAM format.
225
    ///
226
    /// # Examples
227
    ///
228
    /// Set the required fields to RNAME and FLAG,
229
    /// potentially allowing htslib to skip over the rest,
230
    /// resulting in faster iteration:
231
    /// ```
232
    /// use rust_htslib::bam::{Read, Reader};
233
    /// use hts_sys;
234
    /// let mut cram = Reader::from_path("test/test_cram.cram").unwrap();
235
    /// cram.set_cram_options(hts_sys::hts_fmt_option_CRAM_OPT_REQUIRED_FIELDS,
236
    ///             hts_sys::sam_fields_SAM_RNAME | hts_sys::sam_fields_SAM_FLAG).unwrap();
237
    /// ```
238
    fn set_cram_options(&mut self, fmt_opt: hts_fmt_option, fields: sam_fields) -> Result<()> {
8✔
239
        unsafe {
240
            if hts_sys::hts_set_opt(self.htsfile(), fmt_opt, fields) != 0 {
16✔
241
                Err(Error::HtsSetOpt)
×
242
            } else {
243
                Ok(())
5✔
244
            }
245
        }
246
    }
247
}
248

249
/// A BAM reader.
250
#[derive(Debug)]
251
pub struct Reader {
252
    htsfile: *mut htslib::htsFile,
253
    header: Rc<HeaderView>,
254
    tpool: Option<ThreadPool>,
255
}
256

257
unsafe impl Send for Reader {}
258

259
impl Reader {
260
    /// Create a new Reader from path.
261
    ///
262
    /// # Arguments
263
    ///
264
    /// * `path` - the path to open.
265
    pub fn from_path<P: AsRef<Path>>(path: P) -> Result<Self> {
66✔
266
        Self::new(&path_as_bytes(path, true)?)
188✔
267
    }
268

269
    /// Create a new Reader from STDIN.
270
    pub fn from_stdin() -> Result<Self> {
×
271
        Self::new(b"-")
×
272
    }
273

274
    /// Create a new Reader from URL.
275
    pub fn from_url(url: &Url) -> Result<Self> {
×
276
        Self::new(url.as_str().as_bytes())
×
277
    }
278

279
    /// Create a new Reader.
280
    ///
281
    /// # Arguments
282
    ///
283
    /// * `path` - the path to open. Use "-" for stdin.
284
    fn new(path: &[u8]) -> Result<Self> {
64✔
285
        let htsfile = hts_open(path, b"r")?;
232✔
286

287
        let header = unsafe { htslib::sam_hdr_read(htsfile) };
170✔
288
        if header.is_null() {
116✔
289
            return Err(Error::BamOpen {
×
290
                target: String::from_utf8_lossy(path).to_string(),
×
291
            });
292
        }
293

294
        // Invalidate the `text` representation of the header
295
        unsafe {
296
            let _ = htslib::sam_hdr_line_name(header, b"SQ".as_ptr().cast::<c_char>(), 0);
224✔
297
        }
298

299
        Ok(Reader {
62✔
300
            htsfile,
108✔
301
            header: Rc::new(HeaderView::new(header)),
170✔
302
            tpool: None,
54✔
303
        })
304
    }
305

306
    extern "C" fn pileup_read(
16✔
307
        data: *mut ::std::os::raw::c_void,
308
        record: *mut htslib::bam1_t,
309
    ) -> i32 {
310
        let mut _self = unsafe { (data as *mut Self).as_mut().unwrap() };
46✔
311
        unsafe {
312
            htslib::sam_read1(
313
                _self.htsfile(),
30✔
314
                _self.header().inner_ptr() as *mut hts_sys::sam_hdr_t,
16✔
315
                record,
14✔
316
            )
317
        }
318
    }
319

320
    /// Iterator over the records between the (optional) virtual offsets `start` and `end`
321
    ///
322
    /// # Arguments
323
    ///
324
    /// * `start` - Optional starting virtual offset to seek to. Throws an error if it is not
325
    ///   a valid virtual offset.
326
    ///
327
    /// * `end` - Read until the virtual offset is less than `end`
328
    pub fn iter_chunk(&mut self, start: Option<i64>, end: Option<i64>) -> ChunkIterator<'_, Self> {
3✔
329
        if let Some(pos) = start {
3✔
330
            self.seek(pos)
×
331
                .expect("Failed to seek to the starting position");
332
        };
333

334
        ChunkIterator { reader: self, end }
335
    }
336

337
    /// Set the reference path for reading CRAM files.
338
    ///
339
    /// # Arguments
340
    ///
341
    /// * `path` - path to the FASTA reference
342
    pub fn set_reference<P: AsRef<Path>>(&mut self, path: P) -> Result<()> {
3✔
343
        unsafe { set_fai_filename(self.htsfile, path) }
7✔
344
    }
345
}
346

347
impl Read for Reader {
348
    /// Read the next BAM record into the given `Record`.
349
    /// Returns `None` if there are no more records.
350
    ///
351
    /// This method is useful if you want to read records as fast as possible as the
352
    /// `Record` can be reused. A more ergonomic approach is to use the [records](Reader::records)
353
    /// iterator.
354
    ///
355
    /// # Errors
356
    /// If there are any issues with reading the next record an error will be returned.
357
    ///
358
    /// # Examples
359
    ///
360
    /// ```
361
    /// use rust_htslib::errors::Error;
362
    /// use rust_htslib::bam::{Read, Reader, Record};
363
    ///
364
    /// let mut bam = Reader::from_path("test/test.bam")?;
365
    /// let mut record = Record::new();
366
    ///
367
    /// // Print the TID of each record
368
    /// while let Some(r) = bam.read(&mut record) {
369
    ///    r.expect("Failed to parse record");
370
    ///    println!("TID: {}", record.tid())
371
    /// }
372
    /// # Ok::<(), Error>(())
373
    /// ```
374
    fn read(&mut self, record: &mut record::Record) -> Option<Result<()>> {
30,793✔
375
        match unsafe {
30,789✔
376
            htslib::sam_read1(
61,573✔
377
                self.htsfile,
61,573✔
378
                self.header().inner_ptr() as *mut hts_sys::sam_hdr_t,
61,573✔
379
                record.inner_ptr_mut(),
30,789✔
380
            )
381
        } {
382
            -1 => None,
42✔
383
            -2 => Some(Err(Error::BamTruncatedRecord)),
×
384
            -4 => Some(Err(Error::BamInvalidRecord)),
×
385
            _ => {
386
                record.set_header(Rc::clone(&self.header));
122,993✔
387

388
                Some(Ok(()))
30,752✔
389
            }
390
        }
391
    }
392

393
    /// Iterator over the records of the fetched region.
394
    /// Note that, while being convenient, this is less efficient than pre-allocating a
395
    /// `Record` and reading into it with the `read` method, since every iteration involves
396
    /// the allocation of a new `Record`.
397
    fn records(&mut self) -> Records<'_, Self> {
41✔
398
        Records { reader: self }
399
    }
400

401
    fn rc_records(&mut self) -> RcRecords<'_, Self> {
5✔
402
        RcRecords {
403
            reader: self,
404
            record: Rc::new(record::Record::new()),
5✔
405
        }
406
    }
407

408
    fn pileup(&mut self) -> pileup::Pileups<'_, Self> {
6✔
409
        let _self = self as *const Self;
10✔
410
        let itr = unsafe {
411
            htslib::bam_plp_init(
412
                Some(Reader::pileup_read),
4✔
413
                _self as *mut ::std::os::raw::c_void,
4✔
414
            )
415
        };
416
        pileup::Pileups::new(self, itr)
14✔
417
    }
418

419
    fn htsfile(&self) -> *mut htslib::htsFile {
33✔
420
        self.htsfile
33✔
421
    }
422

423
    fn header(&self) -> &HeaderView {
30,822✔
424
        &self.header
30,822✔
425
    }
426

427
    fn set_thread_pool(&mut self, tpool: &ThreadPool) -> Result<()> {
3✔
428
        unsafe { set_thread_pool(self.htsfile(), tpool)? }
9✔
429
        self.tpool = Some(tpool.clone());
6✔
430
        Ok(())
3✔
431
    }
432
}
433

434
impl Drop for Reader {
435
    fn drop(&mut self) {
62✔
436
        unsafe {
437
            htslib::hts_close(self.htsfile);
62✔
438
        }
439
    }
440
}
441

442
/// Conversion type for start/stop coordinates
443
/// only public because it's leaked by the conversions
444
#[doc(hidden)]
445
pub struct FetchCoordinate(i64);
446

447
//the old sam spec
448
impl From<i32> for FetchCoordinate {
449
    fn from(coord: i32) -> FetchCoordinate {
32✔
450
        FetchCoordinate(coord as i64)
32✔
451
    }
452
}
453

454
// to support un-annotated literals (type interference fails on those)
455
impl From<u32> for FetchCoordinate {
456
    fn from(coord: u32) -> FetchCoordinate {
7✔
457
        FetchCoordinate(coord as i64)
7✔
458
    }
459
}
460

461
//the new sam spec
462
impl From<i64> for FetchCoordinate {
463
    fn from(coord: i64) -> FetchCoordinate {
5✔
464
        FetchCoordinate(coord)
4✔
465
    }
466
}
467

468
//what some of our header methods return
469
impl From<u64> for FetchCoordinate {
470
    fn from(coord: u64) -> FetchCoordinate {
5✔
471
        FetchCoordinate(coord.try_into().expect("Coordinate exceeded 2^^63-1"))
13✔
472
    }
473
}
474

475
/// Enum for [IndexdReader.fetch()](struct.IndexedReader.html#method.fetch) arguments.
476
///
477
/// tids may be converted From<>:
478
/// * i32 (correct as per spec)
479
/// * u32 (because of header.tid. Will panic if above 2^31-1).
480
///
481
///Coordinates may be (via FetchCoordinate)
482
/// * i32 (as of the sam v1 spec)
483
/// * i64 (as of the htslib 'large coordinate' extension (even though they are not supported in BAM)
484
/// * u32 (because that's what rust literals will default to)
485
/// * u64 (because of header.target_len(). Will panic if above 2^^63-1).
486
#[derive(Debug)]
487
pub enum FetchDefinition<'a> {
488
    /// tid, start, stop,
489
    Region(i32, i64, i64),
490
    /// 'named-reference', start, stop tuple.
491
    RegionString(&'a [u8], i64, i64),
492
    ///complete reference. May be i32 or u32 (which panics if above 2^31-')
493
    CompleteTid(i32),
494
    ///complete reference by name (&[u8] or &str)
495
    String(&'a [u8]),
496
    /// Every read
497
    All,
498
    /// Only reads with the BAM flag BAM_FUNMAP (which might not be all reads with reference = -1)
499
    Unmapped,
500
}
501

502
impl<'a, X: Into<FetchCoordinate>, Y: Into<FetchCoordinate>> From<(i32, X, Y)>
503
    for FetchDefinition<'a>
504
{
505
    fn from(tup: (i32, X, Y)) -> FetchDefinition<'a> {
5✔
506
        let start: FetchCoordinate = tup.1.into();
16✔
507
        let stop: FetchCoordinate = tup.2.into();
14✔
508
        FetchDefinition::Region(tup.0, start.0, stop.0)
8✔
509
    }
510
}
511

512
impl<'a, X: Into<FetchCoordinate>, Y: Into<FetchCoordinate>> From<(u32, X, Y)>
513
    for FetchDefinition<'a>
514
{
515
    fn from(tup: (u32, X, Y)) -> FetchDefinition<'a> {
13✔
516
        let start: FetchCoordinate = tup.1.into();
46✔
517
        let stop: FetchCoordinate = tup.2.into();
43✔
518
        FetchDefinition::Region(
519
            tup.0.try_into().expect("Tid exceeded 2^31-1"),
53✔
520
            start.0,
10✔
521
            stop.0,
10✔
522
        )
523
    }
524
}
525

526
//non tuple impls
527
impl<'a> From<i32> for FetchDefinition<'a> {
528
    fn from(tid: i32) -> FetchDefinition<'a> {
3✔
529
        FetchDefinition::CompleteTid(tid)
3✔
530
    }
531
}
532

533
impl<'a> From<u32> for FetchDefinition<'a> {
534
    fn from(tid: u32) -> FetchDefinition<'a> {
3✔
535
        let tid: i32 = tid.try_into().expect("tid exceeded 2^31-1");
13✔
536
        FetchDefinition::CompleteTid(tid)
3✔
537
    }
538
}
539

540
impl<'a> From<&'a str> for FetchDefinition<'a> {
541
    fn from(s: &'a str) -> FetchDefinition<'a> {
7✔
542
        FetchDefinition::String(s.as_bytes())
7✔
543
    }
544
}
545

546
//also accept &[u8;n] literals
547
impl<'a> From<&'a [u8]> for FetchDefinition<'a> {
548
    fn from(s: &'a [u8]) -> FetchDefinition<'a> {
7✔
549
        FetchDefinition::String(s)
7✔
550
    }
551
}
552

553
//also accept &[u8;n] literals
554
impl<'a, T: AsRef<[u8]>> From<&'a T> for FetchDefinition<'a> {
555
    fn from(s: &'a T) -> FetchDefinition<'a> {
3✔
556
        FetchDefinition::String(s.as_ref())
3✔
557
    }
558
}
559

560
impl<'a, X: Into<FetchCoordinate>, Y: Into<FetchCoordinate>> From<(&'a str, X, Y)>
561
    for FetchDefinition<'a>
562
{
563
    fn from(tup: (&'a str, X, Y)) -> FetchDefinition<'a> {
7✔
564
        let start: FetchCoordinate = tup.1.into();
22✔
565
        let stop: FetchCoordinate = tup.2.into();
19✔
566
        FetchDefinition::RegionString(tup.0.as_bytes(), start.0, stop.0)
11✔
567
    }
568
}
569

570
impl<'a, X: Into<FetchCoordinate>, Y: Into<FetchCoordinate>> From<(&'a [u8], X, Y)>
571
    for FetchDefinition<'a>
572
{
573
    fn from(tup: (&'a [u8], X, Y)) -> FetchDefinition<'a> {
3✔
574
        let start: FetchCoordinate = tup.1.into();
10✔
575
        let stop: FetchCoordinate = tup.2.into();
9✔
576
        FetchDefinition::RegionString(tup.0, start.0, stop.0)
5✔
577
    }
578
}
579

580
//also accept &[u8;n] literals
581
impl<'a, T: AsRef<[u8]>, X: Into<FetchCoordinate>, Y: Into<FetchCoordinate>> From<(&'a T, X, Y)>
582
    for FetchDefinition<'a>
583
{
584
    fn from(tup: (&'a T, X, Y)) -> FetchDefinition<'a> {
3✔
585
        let start: FetchCoordinate = tup.1.into();
10✔
586
        let stop: FetchCoordinate = tup.2.into();
9✔
587
        FetchDefinition::RegionString(tup.0.as_ref(), start.0, stop.0)
5✔
588
    }
589
}
590

591
#[derive(Debug)]
592
pub struct IndexedReader {
593
    htsfile: *mut htslib::htsFile,
594
    header: Rc<HeaderView>,
595
    idx: Rc<IndexView>,
596
    itr: Option<*mut htslib::hts_itr_t>,
597
    tpool: Option<ThreadPool>,
598
}
599

600
unsafe impl Send for IndexedReader {}
601

602
impl IndexedReader {
603
    /// Create a new Reader from path.
604
    ///
605
    /// # Arguments
606
    ///
607
    /// * `path` - the path to open.
608
    pub fn from_path<P: AsRef<Path>>(path: P) -> Result<Self> {
16✔
609
        Self::new(&path_as_bytes(path, true)?)
44✔
610
    }
611

612
    pub fn from_path_and_index<P: AsRef<Path>>(path: P, index_path: P) -> Result<Self> {
2✔
613
        Self::new_with_index_path(
614
            &path_as_bytes(path, true)?,
4✔
615
            &path_as_bytes(index_path, true)?,
3✔
616
        )
617
    }
618

619
    pub fn from_url(url: &Url) -> Result<Self> {
×
620
        Self::new(url.as_str().as_bytes())
×
621
    }
622

623
    /// Create a new Reader.
624
    ///
625
    /// # Arguments
626
    ///
627
    /// * `path` - the path. Use "-" for stdin.
628
    fn new(path: &[u8]) -> Result<Self> {
16✔
629
        let htsfile = hts_open(path, b"r")?;
52✔
630
        let header = unsafe { htslib::sam_hdr_read(htsfile) };
37✔
631
        let c_str = ffi::CString::new(path).unwrap();
48✔
632
        let idx = unsafe { htslib::sam_index_load(htsfile, c_str.as_ptr()) };
52✔
633
        if idx.is_null() {
26✔
634
            Err(Error::BamInvalidIndex {
×
635
                target: str::from_utf8(path).unwrap().to_owned(),
×
636
            })
637
        } else {
638
            Ok(IndexedReader {
15✔
639
                htsfile,
22✔
640
                header: Rc::new(HeaderView::new(header)),
52✔
641
                idx: Rc::new(IndexView::new(idx)),
52✔
642
                itr: None,
11✔
643
                tpool: None,
11✔
644
            })
645
        }
646
    }
647
    /// Create a new Reader.
648
    ///
649
    /// # Arguments
650
    ///
651
    /// * `path` - the path. Use "-" for stdin.
652
    /// * `index_path` - the index path to use
653
    fn new_with_index_path(path: &[u8], index_path: &[u8]) -> Result<Self> {
2✔
654
        let htsfile = hts_open(path, b"r")?;
5✔
655
        let header = unsafe { htslib::sam_hdr_read(htsfile) };
4✔
656
        let c_str_path = ffi::CString::new(path).unwrap();
5✔
657
        let c_str_index_path = ffi::CString::new(index_path).unwrap();
6✔
658
        let idx = unsafe {
659
            htslib::sam_index_load2(htsfile, c_str_path.as_ptr(), c_str_index_path.as_ptr())
6✔
660
        };
661
        if idx.is_null() {
3✔
662
            Err(Error::BamInvalidIndex {
×
663
                target: str::from_utf8(path).unwrap().to_owned(),
×
664
            })
665
        } else {
666
            Ok(IndexedReader {
2✔
667
                htsfile,
2✔
668
                header: Rc::new(HeaderView::new(header)),
6✔
669
                idx: Rc::new(IndexView::new(idx)),
6✔
670
                itr: None,
1✔
671
                tpool: None,
1✔
672
            })
673
        }
674
    }
675

676
    /// Define the region from which .read() or .records will retrieve reads.
677
    ///
678
    /// Both iterating (with [.records()](trait.Read.html#tymethod.records)) and looping without allocation (with [.read()](trait.Read.html#tymethod.read) are a two stage process:
679
    /// 1. 'fetch' the region of interest
680
    /// 2. iter/loop through the reads.
681
    ///
682
    /// Example:
683
    /// ```
684
    /// use rust_htslib::bam::{IndexedReader, Read};
685
    /// let mut bam = IndexedReader::from_path("test/test.bam").unwrap();
686
    /// bam.fetch(("chrX", 10000, 20000)); // coordinates 10000..20000 on reference named "chrX"
687
    /// for read in bam.records() {
688
    ///     println!("read name: {:?}", read.unwrap().qname());
689
    /// }
690
    /// ```
691
    ///
692
    /// The arguments may be anything that can be converted into a FetchDefinition
693
    /// such as
694
    ///
695
    /// * fetch(tid: u32) -> fetch everything on this reference
696
    /// * fetch(reference_name: &[u8] | &str) -> fetch everything on this reference
697
    /// * fetch((tid: i32, start: i64, stop: i64)): -> fetch in this region on this tid
698
    /// * fetch((reference_name: &[u8] | &str, start: i64, stop: i64) -> fetch in this region on this tid
699
    /// * fetch(FetchDefinition::All) or fetch(".") -> Fetch overything
700
    /// * fetch(FetchDefinition::Unmapped) or fetch("*") -> Fetch unmapped (as signified by the 'unmapped' flag in the BAM - might be unreliable with some aligners.
701
    ///
702
    /// The start / stop coordinates will take i64 (the correct type as of htslib's 'large
703
    /// coordinates' expansion), i32, u32, and u64 (with a possible panic! if the coordinate
704
    /// won't fit an i64).
705
    ///
706
    /// `start` and `stop` are zero-based. `start` is inclusive, `stop` is exclusive.
707
    ///
708
    /// This replaces the old fetch and fetch_str implementations.
709
    pub fn fetch<'a, T: Into<FetchDefinition<'a>>>(&mut self, fetch_definition: T) -> Result<()> {
59✔
710
        //this 'compile time redirect' safes us
711
        //from monomorphing the 'meat' of the fetch function
712
        self._inner_fetch(fetch_definition.into())
188✔
713
    }
714

715
    fn _inner_fetch(&mut self, fetch_definition: FetchDefinition) -> Result<()> {
47✔
716
        match fetch_definition {
47✔
717
            FetchDefinition::Region(tid, start, stop) => {
41✔
718
                self._fetch_by_coord_tuple(tid, start, stop)
67✔
719
            }
720
            FetchDefinition::RegionString(s, start, stop) => {
27✔
721
                let tid = self.header().tid(s);
35✔
722
                match tid {
11✔
723
                    Some(tid) => self._fetch_by_coord_tuple(tid as i32, start, stop),
44✔
724
                    None => Err(Error::Fetch),
2✔
725
                }
726
            }
727
            FetchDefinition::CompleteTid(tid) => {
5✔
728
                let len = self.header().target_len(tid as u32);
17✔
729
                match len {
5✔
730
                    Some(len) => self._fetch_by_coord_tuple(tid, 0, len as i64),
21✔
731
                    None => Err(Error::Fetch),
×
732
                }
733
            }
734
            FetchDefinition::String(s) => {
15✔
735
                // either a target-name or a samtools style definition
736
                let tid = self.header().tid(s);
57✔
737
                match tid {
15✔
738
                    Some(tid) => {
5✔
739
                        //'large position' spec says target len must will fit into an i64.
740
                        let len: i64 = self.header.target_len(tid).unwrap().try_into().unwrap();
30✔
741
                        self._fetch_by_coord_tuple(tid as i32, 0, len)
17✔
742
                    }
743
                    None => self._fetch_by_str(s),
31✔
744
                }
745
            }
746
            FetchDefinition::All => self._fetch_by_str(b"."),
7✔
747
            FetchDefinition::Unmapped => self._fetch_by_str(b"*"),
7✔
748
        }
749
    }
750

751
    fn _fetch_by_coord_tuple(&mut self, tid: i32, beg: i64, end: i64) -> Result<()> {
31✔
752
        if let Some(itr) = self.itr {
53✔
753
            unsafe { htslib::hts_itr_destroy(itr) }
45✔
754
        }
755
        let itr = unsafe { htslib::sam_itr_queryi(self.index().inner_ptr(), tid, beg, end) };
199✔
756
        if itr.is_null() {
62✔
757
            self.itr = None;
×
758
            Err(Error::Fetch)
×
759
        } else {
760
            self.itr = Some(itr);
31✔
761
            Ok(())
31✔
762
        }
763
    }
764

765
    fn _fetch_by_str(&mut self, region: &[u8]) -> Result<()> {
15✔
766
        if let Some(itr) = self.itr {
29✔
767
            unsafe { htslib::hts_itr_destroy(itr) }
29✔
768
        }
769
        let rstr = ffi::CString::new(region).unwrap();
57✔
770
        let rptr = rstr.as_ptr();
30✔
771
        let itr = unsafe {
772
            htslib::sam_itr_querys(
773
                self.index().inner_ptr(),
29✔
774
                self.header().inner_ptr() as *mut hts_sys::sam_hdr_t,
15✔
775
                rptr,
14✔
776
            )
777
        };
778
        if itr.is_null() {
30✔
779
            self.itr = None;
×
780
            Err(Error::Fetch)
×
781
        } else {
782
            self.itr = Some(itr);
15✔
783
            Ok(())
15✔
784
        }
785
    }
786

787
    extern "C" fn pileup_read(
15✔
788
        data: *mut ::std::os::raw::c_void,
789
        record: *mut htslib::bam1_t,
790
    ) -> i32 {
791
        let _self = unsafe { (data as *mut Self).as_mut().unwrap() };
44✔
792
        match _self.itr {
15✔
793
            Some(itr) => itr_next(_self.htsfile, itr, record), // read fetched region
37✔
794
            None => unsafe {
795
                htslib::sam_read1(
796
                    _self.htsfile,
8✔
797
                    _self.header().inner_ptr() as *mut hts_sys::sam_hdr_t,
8✔
798
                    record,
7✔
799
                )
800
            }, // ordinary reading
801
        }
802
    }
803

804
    /// Set the reference path for reading CRAM files.
805
    ///
806
    /// # Arguments
807
    ///
808
    /// * `path` - path to the FASTA reference
809
    pub fn set_reference<P: AsRef<Path>>(&mut self, path: P) -> Result<()> {
3✔
810
        unsafe { set_fai_filename(self.htsfile, path) }
7✔
811
    }
812

813
    pub fn index(&self) -> &IndexView {
50✔
814
        &self.idx
50✔
815
    }
816

817
    // Analogous to slow_idxstats in samtools, see
818
    // https://github.com/samtools/samtools/blob/556c60fdff977c0e6cadc4c2581661f187098b4d/bam_index.c#L140-L199
819
    unsafe fn slow_idxstats(&mut self) -> Result<Vec<(i64, u64, u64, u64)>> {
3✔
820
        self.set_cram_options(
5✔
821
            hts_sys::hts_fmt_option_CRAM_OPT_REQUIRED_FIELDS,
822
            hts_sys::sam_fields_SAM_RNAME | hts_sys::sam_fields_SAM_FLAG,
2✔
823
        )?;
824
        let header = self.header();
7✔
825
        let h = header.inner;
5✔
826
        let mut ret;
827
        let mut last_tid = -2;
5✔
828
        let fp = self.htsfile();
7✔
829

830
        let nref =
3✔
831
            usize::try_from(hts_sys::sam_hdr_nref(h)).map_err(|_| Error::NoSequencesInReference)?;
8✔
832
        if nref == 0 {
3✔
833
            return Ok(vec![]);
×
834
        }
835
        let mut counts = vec![vec![0; 2]; nref + 1];
10✔
836
        let mut bb: hts_sys::bam1_t = MaybeUninit::zeroed().assume_init();
9✔
837
        let b = &mut bb as *mut hts_sys::bam1_t;
5✔
838
        loop {
1✔
839
            ret = hts_sys::sam_read1(fp, h, b);
43✔
840
            if ret < 0 {
15✔
841
                break;
2✔
842
            }
843
            let tid = (*b).core.tid;
26✔
844
            if tid >= nref as i32 || tid < -1 {
26✔
845
                return Err(Error::InvalidTid { tid });
×
846
            }
847

848
            if tid != last_tid {
14✔
849
                if (last_tid >= -1) && (counts[tid as usize][0] + counts[tid as usize][1]) > 0 {
16✔
850
                    return Err(Error::BamUnsorted);
×
851
                }
852
                last_tid = tid;
7✔
853
            }
854

855
            let idx = if ((*b).core.flag as u32 & hts_sys::BAM_FUNMAP) > 0 {
27✔
856
                1
×
857
            } else {
858
                0
13✔
859
            };
860
            counts[(*b).core.tid as usize][idx] += 1;
13✔
861
        }
862

863
        if ret == -1 {
3✔
864
            let res = (0..nref)
4✔
865
                .map(|i| {
10✔
866
                    (
867
                        i as i64,
12✔
868
                        header.target_len(i as u32).unwrap(),
31✔
869
                        counts[i][0],
13✔
870
                        counts[i][1],
7✔
871
                    )
872
                })
873
                .chain([(-1, 0, counts[nref][0], counts[nref][1])])
7✔
874
                .collect();
875
            Ok(res)
3✔
876
        } else {
877
            Err(Error::SlowIdxStats)
×
878
        }
879
    }
880

881
    /// Similar to samtools idxstats, this returns a vector of tuples
882
    /// containing the target id, length, number of mapped reads, and number of unmapped reads.
883
    /// The last entry in the vector corresponds to the unmapped reads for the entire file, with
884
    /// the tid set to -1.
885
    pub fn index_stats(&mut self) -> Result<Vec<(i64, u64, u64, u64)>> {
4✔
886
        let header = self.header();
10✔
887
        let index = self.index();
10✔
888
        if index.inner_ptr().is_null() {
10✔
889
            panic!("Index is null");
×
890
        }
891
        // the quick index stats method only works for BAM files, not SAM or CRAM
892
        unsafe {
893
            if (*self.htsfile()).format.format != htslib::htsExactFormat_bam {
5✔
894
                return self.slow_idxstats();
5✔
895
            }
896
        }
897
        Ok((0..header.target_count())
3✔
898
            .map(|tid| {
8✔
899
                let (mapped, unmapped) = index.number_mapped_unmapped(tid);
21✔
900
                let tlen = header.target_len(tid).unwrap();
27✔
901
                (tid as i64, tlen, mapped, unmapped)
16✔
902
            })
903
            .chain([(-1, 0, 0, index.number_unmapped())])
3✔
904
            .collect::<_>())
2✔
905
    }
906
}
907

908
#[derive(Debug)]
909
pub struct IndexView {
910
    inner: *mut hts_sys::hts_idx_t,
911
    owned: bool,
912
}
913

914
impl IndexView {
915
    fn new(hts_idx: *mut hts_sys::hts_idx_t) -> Self {
16✔
916
        Self {
917
            inner: hts_idx,
918
            owned: true,
919
        }
920
    }
921

922
    #[inline]
923
    pub fn inner(&self) -> &hts_sys::hts_idx_t {
×
924
        unsafe { self.inner.as_ref().unwrap() }
×
925
    }
926

927
    #[inline]
928
    // Pointer to inner hts_idx_t struct
929
    pub fn inner_ptr(&self) -> *const hts_sys::hts_idx_t {
48✔
930
        self.inner
48✔
931
    }
932

933
    #[inline]
934
    pub fn inner_mut(&mut self) -> &mut hts_sys::hts_idx_t {
×
935
        unsafe { self.inner.as_mut().unwrap() }
×
936
    }
937

938
    #[inline]
939
    // Mutable pointer to hts_idx_t struct
940
    pub fn inner_ptr_mut(&mut self) -> *mut hts_sys::hts_idx_t {
×
941
        self.inner
×
942
    }
943

944
    /// Get the number of mapped and unmapped reads for a given target id
945
    /// FIXME only valid for BAM, not SAM/CRAM
946
    fn number_mapped_unmapped(&self, tid: u32) -> (u64, u64) {
7✔
947
        let (mut mapped, mut unmapped) = (0, 0);
19✔
948
        unsafe {
949
            hts_sys::hts_idx_get_stat(self.inner, tid as i32, &mut mapped, &mut unmapped);
25✔
950
        }
951
        (mapped, unmapped)
7✔
952
    }
953

954
    /// Get the total number of unmapped reads in the file
955
    /// FIXME only valid for BAM, not SAM/CRAM
956
    fn number_unmapped(&self) -> u64 {
3✔
957
        unsafe { hts_sys::hts_idx_get_n_no_coor(self.inner) }
5✔
958
    }
959
}
960

961
impl Drop for IndexView {
962
    fn drop(&mut self) {
16✔
963
        if self.owned {
16✔
964
            unsafe {
965
                htslib::hts_idx_destroy(self.inner);
16✔
966
            }
967
        }
968
    }
969
}
970

971
impl Read for IndexedReader {
972
    fn read(&mut self, record: &mut record::Record) -> Option<Result<()>> {
235✔
973
        match self.itr {
235✔
974
            Some(itr) => {
233✔
975
                match itr_next(self.htsfile, itr, &mut record.inner as *mut htslib::bam1_t) {
695✔
976
                    -1 => None,
42✔
977
                    -2 => Some(Err(Error::BamTruncatedRecord)),
5✔
978
                    -4 => Some(Err(Error::BamInvalidRecord)),
×
979
                    _ => {
980
                        record.set_header(Rc::clone(&self.header));
750✔
981

982
                        Some(Ok(()))
189✔
983
                    }
984
                }
985
            }
986
            None => None,
2✔
987
        }
988
    }
989

990
    /// Iterator over the records of the fetched region.
991
    /// Note that, while being convenient, this is less efficient than pre-allocating a
992
    /// `Record` and reading into it with the `read` method, since every iteration involves
993
    /// the allocation of a new `Record`.
994
    fn records(&mut self) -> Records<'_, Self> {
41✔
995
        Records { reader: self }
996
    }
997

998
    fn rc_records(&mut self) -> RcRecords<'_, Self> {
×
999
        RcRecords {
1000
            reader: self,
1001
            record: Rc::new(record::Record::new()),
×
1002
        }
1003
    }
1004

1005
    fn pileup(&mut self) -> pileup::Pileups<'_, Self> {
3✔
1006
        let _self = self as *const Self;
5✔
1007
        let itr = unsafe {
1008
            htslib::bam_plp_init(
1009
                Some(IndexedReader::pileup_read),
2✔
1010
                _self as *mut ::std::os::raw::c_void,
2✔
1011
            )
1012
        };
1013
        pileup::Pileups::new(self, itr)
7✔
1014
    }
1015

1016
    fn htsfile(&self) -> *mut htslib::htsFile {
8✔
1017
        self.htsfile
8✔
1018
    }
1019

1020
    fn header(&self) -> &HeaderView {
56✔
1021
        &self.header
56✔
1022
    }
1023

1024
    fn set_thread_pool(&mut self, tpool: &ThreadPool) -> Result<()> {
×
1025
        unsafe { set_thread_pool(self.htsfile(), tpool)? }
×
1026
        self.tpool = Some(tpool.clone());
×
1027
        Ok(())
×
1028
    }
1029
}
1030

1031
impl Drop for IndexedReader {
1032
    fn drop(&mut self) {
16✔
1033
        unsafe {
1034
            if self.itr.is_some() {
34✔
1035
                htslib::hts_itr_destroy(self.itr.unwrap());
16✔
1036
            }
1037
            htslib::hts_close(self.htsfile);
28✔
1038
        }
1039
    }
1040
}
1041

1042
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
1043
pub enum Format {
1044
    Sam,
1045
    Bam,
1046
    Cram,
1047
}
1048

1049
impl Format {
1050
    fn write_mode(self) -> &'static [u8] {
15✔
1051
        match self {
17✔
1052
            Format::Sam => b"w",
2✔
1053
            Format::Bam => b"wb",
13✔
1054
            Format::Cram => b"wc",
2✔
1055
        }
1056
    }
1057
}
1058

1059
/// A BAM writer.
1060
#[derive(Debug)]
1061
pub struct Writer {
1062
    f: *mut htslib::htsFile,
1063
    header: Rc<HeaderView>,
1064
    tpool: Option<ThreadPool>,
1065
}
1066

1067
unsafe impl Send for Writer {}
1068

1069
impl Writer {
1070
    /// Create a new SAM/BAM/CRAM file.
1071
    ///
1072
    /// # Arguments
1073
    ///
1074
    /// * `path` - the path.
1075
    /// * `header` - header definition to use
1076
    /// * `format` - the format to use (SAM/BAM/CRAM)
1077
    pub fn from_path<P: AsRef<Path>>(
16✔
1078
        path: P,
1079
        header: &header::Header,
1080
        format: Format,
1081
    ) -> Result<Self> {
1082
        Self::new(&path_as_bytes(path, false)?, format.write_mode(), header)
71✔
1083
    }
1084

1085
    /// Create a new SAM/BAM/CRAM file at STDOUT.
1086
    ///
1087
    /// # Arguments
1088
    ///
1089
    /// * `header` - header definition to use
1090
    /// * `format` - the format to use (SAM/BAM/CRAM)
1091
    pub fn from_stdout(header: &header::Header, format: Format) -> Result<Self> {
×
1092
        Self::new(b"-", format.write_mode(), header)
×
1093
    }
1094

1095
    /// Create a new SAM/BAM/CRAM file.
1096
    ///
1097
    /// # Arguments
1098
    ///
1099
    /// * `path` - the path. Use "-" for stdout.
1100
    /// * `mode` - write mode, refer to htslib::hts_open()
1101
    /// * `header` - header definition to use
1102
    fn new(path: &[u8], mode: &[u8], header: &header::Header) -> Result<Self> {
15✔
1103
        let f = hts_open(path, mode)?;
54✔
1104

1105
        // sam_hdr_parse does not populate the text and l_text fields of the header_record.
1106
        // This causes non-SQ headers to be dropped in the output BAM file.
1107
        // To avoid this, we copy the All header to a new C-string that is allocated with malloc,
1108
        // and set this into header_record manually.
1109
        let header_record = unsafe {
1110
            let mut header_string = header.to_bytes();
41✔
1111
            if !header_string.is_empty() && header_string[header_string.len() - 1] != b'\n' {
58✔
1112
                header_string.push(b'\n');
15✔
1113
            }
1114
            let l_text = header_string.len();
43✔
1115
            let text = ::libc::malloc(l_text + 1);
41✔
1116
            libc::memset(text, 0, l_text + 1);
43✔
1117
            libc::memcpy(
1118
                text,
13✔
1119
                header_string.as_ptr() as *const ::libc::c_void,
15✔
1120
                header_string.len(),
28✔
1121
            );
1122

1123
            //println!("{}", str::from_utf8(&header_string).unwrap());
1124
            let rec = htslib::sam_hdr_parse(l_text + 1, text as *const c_char);
54✔
1125

1126
            (*rec).text = text as *mut c_char;
17✔
1127
            (*rec).l_text = l_text;
17✔
1128
            rec
13✔
1129
        };
1130

1131
        unsafe {
1132
            htslib::sam_hdr_write(f, header_record);
28✔
1133
        }
1134

1135
        Ok(Writer {
15✔
1136
            f,
26✔
1137
            header: Rc::new(HeaderView::new(header_record)),
41✔
1138
            tpool: None,
13✔
1139
        })
1140
    }
1141

1142
    /// Activate multi-threaded BAM write support in htslib. This should permit faster
1143
    /// writing of large BAM files.
1144
    ///
1145
    /// # Arguments
1146
    ///
1147
    /// * `n_threads` - number of extra background writer threads to use, must be `> 0`.
1148
    pub fn set_threads(&mut self, n_threads: usize) -> Result<()> {
2✔
1149
        unsafe { set_threads(self.f, n_threads) }
4✔
1150
    }
1151

1152
    /// Use a shared thread-pool for writing. This permits controlling the total
1153
    /// thread count when multiple readers and writers are working simultaneously.
1154
    /// A thread pool can be created with `crate::tpool::ThreadPool::new(n_threads)`
1155
    ///
1156
    /// # Arguments
1157
    ///
1158
    /// * `tpool` - thread pool to use for compression work.
1159
    pub fn set_thread_pool(&mut self, tpool: &ThreadPool) -> Result<()> {
3✔
1160
        unsafe { set_thread_pool(self.f, tpool)? }
7✔
1161
        self.tpool = Some(tpool.clone());
6✔
1162
        Ok(())
3✔
1163
    }
1164

1165
    /// Write record to BAM.
1166
    ///
1167
    /// # Arguments
1168
    ///
1169
    /// * `record` - the record to write
1170
    pub fn write(&mut self, record: &record::Record) -> Result<()> {
30,056✔
1171
        if unsafe { htslib::sam_write1(self.f, self.header.inner(), record.inner_ptr()) } == -1 {
120,218✔
1172
            Err(Error::WriteRecord)
×
1173
        } else {
1174
            Ok(())
30,056✔
1175
        }
1176
    }
1177

1178
    /// Return the header.
1179
    pub fn header(&self) -> &HeaderView {
×
1180
        &self.header
×
1181
    }
1182

1183
    /// Set the reference path for reading CRAM files.
1184
    ///
1185
    /// # Arguments
1186
    ///
1187
    /// * `path` - path to the FASTA reference
1188
    pub fn set_reference<P: AsRef<Path>>(&mut self, path: P) -> Result<()> {
2✔
1189
        unsafe { set_fai_filename(self.f, path) }
4✔
1190
    }
1191

1192
    /// Set the compression level for writing BAM/CRAM files.
1193
    ///
1194
    /// # Arguments
1195
    ///
1196
    /// * `compression_level` - `CompressionLevel` enum variant
1197
    pub fn set_compression_level(&mut self, compression_level: CompressionLevel) -> Result<()> {
5✔
1198
        let level = compression_level.convert()?;
13✔
1199
        match unsafe {
5✔
1200
            htslib::hts_set_opt(
9✔
1201
                self.f,
5✔
1202
                htslib::hts_fmt_option_HTS_OPT_COMPRESSION_LEVEL,
4✔
1203
                level,
4✔
1204
            )
1205
        } {
1206
            0 => Ok(()),
5✔
1207
            _ => Err(Error::BamInvalidCompressionLevel { level }),
×
1208
        }
1209
    }
1210
}
1211

1212
/// Compression levels in BAM/CRAM files
1213
///
1214
/// * Uncompressed: No compression, zlib level 0
1215
/// * Fastest: Lowest compression level, zlib level 1
1216
/// * Maximum: Highest compression level, zlib level 9
1217
/// * Level(i): Custom compression level in the range [0, 9]
1218
#[derive(Debug, Clone, Copy)]
1219
pub enum CompressionLevel {
1220
    Uncompressed,
1221
    Fastest,
1222
    Maximum,
1223
    Level(u32),
1224
}
1225

1226
impl CompressionLevel {
1227
    // Convert and check the variants of the `CompressionLevel` enum to a numeric level
1228
    fn convert(self) -> Result<u32> {
19✔
1229
        match self {
19✔
1230
            CompressionLevel::Uncompressed => Ok(0),
3✔
1231
            CompressionLevel::Fastest => Ok(1),
3✔
1232
            CompressionLevel::Maximum => Ok(9),
3✔
1233
            CompressionLevel::Level(i @ 0..=9) => Ok(i),
36✔
1234
            CompressionLevel::Level(i) => Err(Error::BamInvalidCompressionLevel { level: i }),
3✔
1235
        }
1236
    }
1237
}
1238

1239
impl Drop for Writer {
1240
    fn drop(&mut self) {
15✔
1241
        unsafe {
1242
            htslib::hts_close(self.f);
15✔
1243
        }
1244
    }
1245
}
1246

1247
/// Iterator over the records of a BAM.
1248
#[derive(Debug)]
1249
pub struct Records<'a, R: Read> {
1250
    reader: &'a mut R,
1251
}
1252

1253
impl<R: Read> Iterator for Records<'_, R> {
1254
    type Item = Result<record::Record>;
1255

1256
    fn next(&mut self) -> Option<Result<record::Record>> {
10,968✔
1257
        let mut record = record::Record::new();
21,931✔
1258
        match self.reader.read(&mut record) {
32,899✔
1259
            None => None,
74✔
1260
            Some(Ok(_)) => Some(Ok(record)),
10,894✔
1261
            Some(Err(err)) => Some(Err(err)),
9✔
1262
        }
1263
    }
1264
}
1265

1266
/// Iterator over the records of a BAM, using an Rc.
1267
///
1268
/// See [rc_records](trait.Read.html#tymethod.rc_records).
1269
#[derive(Debug)]
1270
pub struct RcRecords<'a, R: Read> {
1271
    reader: &'a mut R,
1272
    record: Rc<record::Record>,
1273
}
1274

1275
impl<R: Read> Iterator for RcRecords<'_, R> {
1276
    type Item = Result<Rc<record::Record>>;
1277

1278
    fn next(&mut self) -> Option<Self::Item> {
17✔
1279
        let record = match Rc::get_mut(&mut self.record) {
32✔
1280
            //not make_mut, we don't need a clone
1281
            Some(x) => x,
32✔
1282
            None => {
×
1283
                self.record = Rc::new(record::Record::new());
×
1284
                Rc::get_mut(&mut self.record).unwrap()
×
1285
            }
1286
        };
1287

1288
        match self.reader.read(record) {
47✔
1289
            None => None,
5✔
1290
            Some(Ok(_)) => Some(Ok(Rc::clone(&self.record))),
14✔
1291
            Some(Err(err)) => Some(Err(err)),
×
1292
        }
1293
    }
1294
}
1295

1296
/// Iterator over the records of a BAM until the virtual offset is less than `end`
1297
pub struct ChunkIterator<'a, R: Read> {
1298
    reader: &'a mut R,
1299
    end: Option<i64>,
1300
}
1301

1302
impl<R: Read> Iterator for ChunkIterator<'_, R> {
1303
    type Item = Result<record::Record>;
1304
    fn next(&mut self) -> Option<Result<record::Record>> {
20,003✔
1305
        if let Some(pos) = self.end {
20,003✔
1306
            if self.reader.tell() >= pos {
×
1307
                return None;
×
1308
            }
1309
        }
1310
        let mut record = record::Record::new();
40,005✔
1311
        match self.reader.read(&mut record) {
60,008✔
1312
            None => None,
3✔
1313
            Some(Ok(_)) => Some(Ok(record)),
20,001✔
1314
            Some(Err(err)) => Some(Err(err)),
×
1315
        }
1316
    }
1317
}
1318

1319
/// Wrapper for opening a BAM file.
1320
fn hts_open(path: &[u8], mode: &[u8]) -> Result<*mut htslib::htsFile> {
93✔
1321
    let cpath = ffi::CString::new(path).unwrap();
339✔
1322
    let path = str::from_utf8(path).unwrap();
350✔
1323
    let c_str = ffi::CString::new(mode).unwrap();
339✔
1324
    let ret = unsafe { htslib::hts_open(cpath.as_ptr(), c_str.as_ptr()) };
350✔
1325
    if ret.is_null() {
186✔
1326
        Err(Error::BamOpen {
×
1327
            target: path.to_owned(),
×
1328
        })
1329
    } else {
1330
        if !mode.contains(&b'w') {
186✔
1331
            unsafe {
1332
                // Comparison against 'htsFormatCategory_sequence_data' doesn't handle text files correctly
1333
                // hence the explicit checks against all supported exact formats
1334
                if (*ret).format.format != htslib::htsExactFormat_sam
91✔
1335
                    && (*ret).format.format != htslib::htsExactFormat_bam
70✔
1336
                    && (*ret).format.format != htslib::htsExactFormat_cram
10✔
1337
                {
1338
                    return Err(Error::BamOpen {
4✔
1339
                        target: path.to_owned(),
4✔
1340
                    });
1341
                }
1342
            }
1343
        }
1344
        Ok(ret)
90✔
1345
    }
1346
}
1347

1348
/// Wrapper for iterating an indexed BAM file.
1349
fn itr_next(
240✔
1350
    htsfile: *mut htslib::htsFile,
1351
    itr: *mut htslib::hts_itr_t,
1352
    record: *mut htslib::bam1_t,
1353
) -> i32 {
1354
    unsafe {
1355
        htslib::hts_itr_next(
1356
            (*htsfile).fp.bgzf,
240✔
1357
            itr,
238✔
1358
            record as *mut ::std::os::raw::c_void,
238✔
1359
            htsfile as *mut ::std::os::raw::c_void,
238✔
1360
        )
1361
    }
1362
}
1363

1364
#[derive(Debug)]
1365
pub struct HeaderView {
1366
    inner: *mut htslib::bam_hdr_t,
1367
    owned: bool,
1368
}
1369

1370
impl HeaderView {
1371
    /// Create a new HeaderView from a pre-populated Header object
1372
    pub fn from_header(header: &Header) -> Self {
9✔
1373
        let mut header_string = header.to_bytes();
21✔
1374
        if !header_string.is_empty() && header_string[header_string.len() - 1] != b'\n' {
16✔
1375
            header_string.push(b'\n');
2✔
1376
        }
1377
        Self::from_bytes(&header_string)
18✔
1378
    }
1379

1380
    /// Create a new HeaderView from bytes
1381
    pub fn from_bytes(header_string: &[u8]) -> Self {
9✔
1382
        let header_record = unsafe {
1383
            let l_text = header_string.len();
21✔
1384
            let text = ::libc::malloc(l_text + 1);
21✔
1385
            ::libc::memset(text, 0, l_text + 1);
24✔
1386
            ::libc::memcpy(
1387
                text,
6✔
1388
                header_string.as_ptr() as *const ::libc::c_void,
9✔
1389
                header_string.len(),
12✔
1390
            );
1391

1392
            let rec = htslib::sam_hdr_parse(l_text + 1, text as *const c_char);
30✔
1393
            (*rec).text = text as *mut c_char;
12✔
1394
            (*rec).l_text = l_text;
12✔
1395
            rec
6✔
1396
        };
1397

1398
        HeaderView::new(header_record)
15✔
1399
    }
1400

1401
    /// Create a new HeaderView from the underlying Htslib type, and own it.
1402
    pub fn new(inner: *mut htslib::bam_hdr_t) -> Self {
98✔
1403
        HeaderView { inner, owned: true }
1404
    }
1405

1406
    #[inline]
1407
    pub fn inner(&self) -> &htslib::bam_hdr_t {
30,057✔
1408
        unsafe { self.inner.as_ref().unwrap() }
60,112✔
1409
    }
1410

1411
    #[inline]
1412
    // Pointer to inner bam_hdr_t struct
1413
    pub fn inner_ptr(&self) -> *const htslib::bam_hdr_t {
30,839✔
1414
        self.inner
30,839✔
1415
    }
1416

1417
    #[inline]
1418
    pub fn inner_mut(&mut self) -> &mut htslib::bam_hdr_t {
×
1419
        unsafe { self.inner.as_mut().unwrap() }
×
1420
    }
1421

1422
    #[inline]
1423
    // Mutable pointer to bam_hdr_t struct
1424
    pub fn inner_ptr_mut(&mut self) -> *mut htslib::bam_hdr_t {
×
1425
        self.inner
×
1426
    }
1427

1428
    pub fn tid(&self, name: &[u8]) -> Option<u32> {
31✔
1429
        let c_str = ffi::CString::new(name).expect("Expected valid name.");
143✔
1430
        let tid = unsafe { htslib::sam_hdr_name2tid(self.inner, c_str.as_ptr()) };
118✔
1431
        if tid < 0 {
33✔
1432
            None
13✔
1433
        } else {
1434
            Some(tid as u32)
19✔
1435
        }
1436
    }
1437

1438
    pub fn tid2name(&self, tid: u32) -> &[u8] {
×
1439
        unsafe { ffi::CStr::from_ptr(htslib::sam_hdr_tid2name(self.inner, tid as i32)).to_bytes() }
×
1440
    }
1441

1442
    pub fn target_count(&self) -> u32 {
2✔
1443
        self.inner().n_targets as u32
2✔
1444
    }
1445

1446
    pub fn target_names(&self) -> Vec<&[u8]> {
×
1447
        let names = unsafe {
1448
            slice::from_raw_parts(self.inner().target_name, self.target_count() as usize)
×
1449
        };
1450
        names
1451
            .iter()
1452
            .map(|name| unsafe { ffi::CStr::from_ptr(*name).to_bytes() })
×
1453
            .collect()
1454
    }
1455

1456
    pub fn target_len(&self, tid: u32) -> Option<u64> {
23✔
1457
        let inner = unsafe { *self.inner };
45✔
1458
        if (tid as i32) < inner.n_targets {
24✔
1459
            let l: &[u32] =
45✔
1460
                unsafe { slice::from_raw_parts(inner.target_len, inner.n_targets as usize) };
44✔
1461
            Some(l[tid as usize] as u64)
24✔
1462
        } else {
1463
            None
×
1464
        }
1465
    }
1466

1467
    /// Retrieve the textual SAM header as bytes
1468
    pub fn as_bytes(&self) -> &[u8] {
15✔
1469
        unsafe {
1470
            let rebuilt_hdr = htslib::sam_hdr_str(self.inner);
39✔
1471
            if rebuilt_hdr.is_null() {
27✔
1472
                return b"";
×
1473
            }
1474
            ffi::CStr::from_ptr(rebuilt_hdr).to_bytes()
15✔
1475
        }
1476
    }
1477
}
1478

1479
impl Clone for HeaderView {
1480
    fn clone(&self) -> Self {
×
1481
        HeaderView {
1482
            inner: unsafe { htslib::sam_hdr_dup(self.inner) },
×
1483
            owned: true,
1484
        }
1485
    }
1486
}
1487

1488
impl Drop for HeaderView {
1489
    fn drop(&mut self) {
98✔
1490
        if self.owned {
98✔
1491
            unsafe {
1492
                htslib::sam_hdr_destroy(self.inner);
98✔
1493
            }
1494
        }
1495
    }
1496
}
1497

1498
#[cfg(test)]
1499
mod tests {
1500
    use super::header::HeaderRecord;
1501
    use super::record::{Aux, Cigar, CigarString};
1502
    use super::*;
1503
    use std::collections::HashMap;
1504
    use std::fs;
1505
    use std::path::Path;
1506
    use std::str;
1507

1508
    type GoldType = (
1509
        [&'static [u8]; 6],
1510
        [u16; 6],
1511
        [&'static [u8]; 6],
1512
        [&'static [u8]; 6],
1513
        [CigarString; 6],
1514
    );
1515
    fn gold() -> GoldType {
1516
        let names = [
1517
            &b"I"[..],
1518
            &b"II.14978392"[..],
1519
            &b"III"[..],
1520
            &b"IV"[..],
1521
            &b"V"[..],
1522
            &b"VI"[..],
1523
        ];
1524
        let flags = [16u16, 16u16, 16u16, 16u16, 16u16, 2048u16];
1525
        let seqs = [
1526
            &b"CCTAGCCCTAACCCTAACCCTAACCCTAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCC\
1527
TAAGCCTAAGCCTAAGCCTAA"[..],
1528
            &b"CCTAGCCCTAACCCTAACCCTAACCCTAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCC\
1529
TAAGCCTAAGCCTAAGCCTAA"[..],
1530
            &b"CCTAGCCCTAACCCTAACCCTAACCCTAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCC\
1531
TAAGCCTAAGCCTAAGCCTAA"[..],
1532
            &b"CCTAGCCCTAACCCTAACCCTAACCCTAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCC\
1533
TAAGCCTAAGCCTAAGCCTAA"[..],
1534
            &b"CCTAGCCCTAACCCTAACCCTAACCCTAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCC\
1535
TAAGCCTAAGCCTAAGCCTAA"[..],
1536
            &b"ACTAAGCCTAAGCCTAAGCCTAAGCCAATTATCGATTTCTGAAAAAATTATCGAATTTTCTAGAAATTTTGCAAATTTT\
1537
TTCATAAAATTATCGATTTTA"[..],
1538
        ];
1539
        let quals = [
1540
            &b"#############################@B?8B?BA@@DDBCDDCBC@CDCDCCCCCCCCCCCCCCCCCCCCCCCCCCCC\
1541
CCCCCCCCCCCCCCCCCCC"[..],
1542
            &b"#############################@B?8B?BA@@DDBCDDCBC@CDCDCCCCCCCCCCCCCCCCCCCCCCCCCCCC\
1543
CCCCCCCCCCCCCCCCCCC"[..],
1544
            &b"#############################@B?8B?BA@@DDBCDDCBC@CDCDCCCCCCCCCCCCCCCCCCCCCCCCCCCC\
1545
CCCCCCCCCCCCCCCCCCC"[..],
1546
            &b"#############################@B?8B?BA@@DDBCDDCBC@CDCDCCCCCCCCCCCCCCCCCCCCCCCCCCCC\
1547
CCCCCCCCCCCCCCCCCCC"[..],
1548
            &b"#############################@B?8B?BA@@DDBCDDCBC@CDCDCCCCCCCCCCCCCCCCCCCCCCCCCCCC\
1549
CCCCCCCCCCCCCCCCCCC"[..],
1550
            &b"#############################@B?8B?BA@@DDBCDDCBC@CDCDCCCCCCCCCCCCCCCCCCCCCCCCCCCC\
1551
CCCCCCCCCCCCCCCCCCC"[..],
1552
        ];
1553
        let cigars = [
1554
            CigarString(vec![Cigar::Match(27), Cigar::Del(1), Cigar::Match(73)]),
1555
            CigarString(vec![Cigar::Match(27), Cigar::Del(1), Cigar::Match(73)]),
1556
            CigarString(vec![Cigar::Match(27), Cigar::Del(1), Cigar::Match(73)]),
1557
            CigarString(vec![Cigar::Match(27), Cigar::Del(1), Cigar::Match(73)]),
1558
            CigarString(vec![Cigar::Match(27), Cigar::Del(1), Cigar::Match(73)]),
1559
            CigarString(vec![Cigar::Match(27), Cigar::Del(100000), Cigar::Match(73)]),
1560
        ];
1561
        (names, flags, seqs, quals, cigars)
1562
    }
1563

1564
    fn compare_inner_bam_cram_records(cram_records: &[Record], bam_records: &[Record]) {
1565
        // Selectively compares bam1_t struct fields from BAM and CRAM
1566
        for (c1, b1) in cram_records.iter().zip(bam_records.iter()) {
1567
            // CRAM vs BAM l_data is off by 3, see: https://github.com/rust-bio/rust-htslib/pull/184#issuecomment-590133544
1568
            // The rest of the fields should be identical:
1569
            assert_eq!(c1.cigar(), b1.cigar());
1570
            assert_eq!(c1.inner().core.pos, b1.inner().core.pos);
1571
            assert_eq!(c1.inner().core.mpos, b1.inner().core.mpos);
1572
            assert_eq!(c1.inner().core.mtid, b1.inner().core.mtid);
1573
            assert_eq!(c1.inner().core.tid, b1.inner().core.tid);
1574
            assert_eq!(c1.inner().core.bin, b1.inner().core.bin);
1575
            assert_eq!(c1.inner().core.qual, b1.inner().core.qual);
1576
            assert_eq!(c1.inner().core.l_extranul, b1.inner().core.l_extranul);
1577
            assert_eq!(c1.inner().core.flag, b1.inner().core.flag);
1578
            assert_eq!(c1.inner().core.l_qname, b1.inner().core.l_qname);
1579
            assert_eq!(c1.inner().core.n_cigar, b1.inner().core.n_cigar);
1580
            assert_eq!(c1.inner().core.l_qseq, b1.inner().core.l_qseq);
1581
            assert_eq!(c1.inner().core.isize_, b1.inner().core.isize_);
1582
            //... except m_data
1583
        }
1584
    }
1585

1586
    #[test]
1587
    fn test_read() {
1588
        let (names, flags, seqs, quals, cigars) = gold();
1589
        let mut bam = Reader::from_path(Path::new("test/test.bam")).expect("Error opening file.");
1590
        let del_len = [1, 1, 1, 1, 1, 100000];
1591

1592
        for (i, record) in bam.records().enumerate() {
1593
            let rec = record.expect("Expected valid record");
1594
            assert_eq!(rec.qname(), names[i]);
1595
            assert_eq!(rec.flags(), flags[i]);
1596
            assert_eq!(rec.seq().as_bytes(), seqs[i]);
1597

1598
            let cigar = rec.cigar();
1599
            assert_eq!(*cigar, cigars[i]);
1600

1601
            let end_pos = cigar.end_pos();
1602
            assert_eq!(end_pos, rec.pos() + 100 + del_len[i]);
1603
            assert_eq!(
1604
                cigar
1605
                    .read_pos(end_pos as u32 - 10, false, false)
1606
                    .unwrap()
1607
                    .unwrap(),
1608
                90
1609
            );
1610
            assert_eq!(
1611
                cigar
1612
                    .read_pos(rec.pos() as u32 + 20, false, false)
1613
                    .unwrap()
1614
                    .unwrap(),
1615
                20
1616
            );
1617
            assert_eq!(cigar.read_pos(4000000, false, false).unwrap(), None);
1618
            // fix qual offset
1619
            let qual: Vec<u8> = quals[i].iter().map(|&q| q - 33).collect();
1620
            assert_eq!(rec.qual(), &qual[..]);
1621
        }
1622
    }
1623

1624
    #[test]
1625
    fn test_seek() {
1626
        let mut bam = Reader::from_path(Path::new("test/test.bam")).expect("Error opening file.");
1627

1628
        let mut names_by_voffset = HashMap::new();
1629

1630
        let mut offset = bam.tell();
1631
        let mut rec = Record::new();
1632
        while let Some(r) = bam.read(&mut rec) {
1633
            r.expect("error reading bam");
1634
            let qname = str::from_utf8(rec.qname()).unwrap().to_string();
1635
            println!("{} {}", offset, qname);
1636
            names_by_voffset.insert(offset, qname);
1637
            offset = bam.tell();
1638
        }
1639

1640
        for (offset, qname) in names_by_voffset.iter() {
1641
            println!("{} {}", offset, qname);
1642
            bam.seek(*offset).unwrap();
1643
            if let Some(r) = bam.read(&mut rec) {
1644
                r.unwrap();
1645
            };
1646
            let rec_qname = str::from_utf8(rec.qname()).unwrap().to_string();
1647
            assert_eq!(qname, &rec_qname);
1648
        }
1649
    }
1650

1651
    #[test]
1652
    fn test_read_sam_header() {
1653
        let bam = Reader::from_path("test/test.bam").expect("Error opening file.");
1654

1655
        let true_header = "@SQ\tSN:CHROMOSOME_I\tLN:15072423\n@SQ\tSN:CHROMOSOME_II\tLN:15279345\
1656
             \n@SQ\tSN:CHROMOSOME_III\tLN:13783700\n@SQ\tSN:CHROMOSOME_IV\tLN:17493793\n@SQ\t\
1657
             SN:CHROMOSOME_V\tLN:20924149\n"
1658
            .to_string();
1659
        let header_text = String::from_utf8(bam.header.as_bytes().to_owned()).unwrap();
1660
        assert_eq!(header_text, true_header);
1661
    }
1662

1663
    #[test]
1664
    fn test_read_against_sam() {
1665
        let mut bam = Reader::from_path("./test/bam2sam_out.sam").unwrap();
1666
        for read in bam.records() {
1667
            let _read = read.unwrap();
1668
        }
1669
    }
1670

1671
    fn _test_read_indexed_common(mut bam: IndexedReader) {
1672
        let (names, flags, seqs, quals, cigars) = gold();
1673
        let sq_1 = b"CHROMOSOME_I";
1674
        let sq_2 = b"CHROMOSOME_II";
1675
        let tid_1 = bam.header.tid(sq_1).expect("Expected tid.");
1676
        let tid_2 = bam.header.tid(sq_2).expect("Expected tid.");
1677
        assert!(bam.header.target_len(tid_1).expect("Expected target len.") == 15072423);
1678

1679
        // fetch to position containing reads
1680
        bam.fetch((tid_1, 0, 2))
1681
            .expect("Expected successful fetch.");
1682
        assert!(bam.records().count() == 6);
1683

1684
        // compare reads
1685
        bam.fetch((tid_1, 0, 2))
1686
            .expect("Expected successful fetch.");
1687
        for (i, record) in bam.records().enumerate() {
1688
            let rec = record.expect("Expected valid record");
1689

1690
            println!("{}", str::from_utf8(rec.qname()).unwrap());
1691
            assert_eq!(rec.qname(), names[i]);
1692
            assert_eq!(rec.flags(), flags[i]);
1693
            assert_eq!(rec.seq().as_bytes(), seqs[i]);
1694
            assert_eq!(*rec.cigar(), cigars[i]);
1695
            // fix qual offset
1696
            let qual: Vec<u8> = quals[i].iter().map(|&q| q - 33).collect();
1697
            assert_eq!(rec.qual(), &qual[..]);
1698
            assert_eq!(rec.aux(b"NotAvailableAux"), Err(Error::BamAuxTagNotFound));
1699
        }
1700

1701
        // fetch to empty position
1702
        bam.fetch((tid_2, 1, 1))
1703
            .expect("Expected successful fetch.");
1704
        assert!(bam.records().count() == 0);
1705

1706
        // repeat with byte-string based fetch
1707

1708
        // fetch to position containing reads
1709
        // using coordinate-string chr:start-stop
1710
        bam.fetch(format!("{}:{}-{}", str::from_utf8(sq_1).unwrap(), 0, 2).as_bytes())
1711
            .expect("Expected successful fetch.");
1712
        assert!(bam.records().count() == 6);
1713
        // using &str and exercising some of the coordinate conversion funcs
1714
        bam.fetch((str::from_utf8(sq_1).unwrap(), 0_u32, 2_u64))
1715
            .expect("Expected successful fetch.");
1716
        assert!(bam.records().count() == 6);
1717
        // using a slice
1718
        bam.fetch((&sq_1[..], 0, 2))
1719
            .expect("Expected successful fetch.");
1720
        assert!(bam.records().count() == 6);
1721
        // using a literal
1722
        bam.fetch((sq_1, 0, 2)).expect("Expected successful fetch.");
1723
        assert!(bam.records().count() == 6);
1724

1725
        // using a tid
1726
        bam.fetch((0i32, 0u32, 2i64))
1727
            .expect("Expected successful fetch.");
1728
        assert!(bam.records().count() == 6);
1729
        // using a tid:u32
1730
        bam.fetch((0u32, 0u32, 2i64))
1731
            .expect("Expected successful fetch.");
1732
        assert!(bam.records().count() == 6);
1733

1734
        // compare reads
1735
        bam.fetch(format!("{}:{}-{}", str::from_utf8(sq_1).unwrap(), 0, 2).as_bytes())
1736
            .expect("Expected successful fetch.");
1737
        for (i, record) in bam.records().enumerate() {
1738
            let rec = record.expect("Expected valid record");
1739

1740
            println!("{}", str::from_utf8(rec.qname()).unwrap());
1741
            assert_eq!(rec.qname(), names[i]);
1742
            assert_eq!(rec.flags(), flags[i]);
1743
            assert_eq!(rec.seq().as_bytes(), seqs[i]);
1744
            assert_eq!(*rec.cigar(), cigars[i]);
1745
            // fix qual offset
1746
            let qual: Vec<u8> = quals[i].iter().map(|&q| q - 33).collect();
1747
            assert_eq!(rec.qual(), &qual[..]);
1748
            assert_eq!(rec.aux(b"NotAvailableAux"), Err(Error::BamAuxTagNotFound));
1749
        }
1750

1751
        // fetch to empty position
1752
        bam.fetch(format!("{}:{}-{}", str::from_utf8(sq_2).unwrap(), 1, 1).as_bytes())
1753
            .expect("Expected successful fetch.");
1754
        assert!(bam.records().count() == 0);
1755

1756
        //all on a tid
1757
        bam.fetch(0).expect("Expected successful fetch.");
1758
        assert!(bam.records().count() == 6);
1759
        //all on a tid:u32
1760
        bam.fetch(0u32).expect("Expected successful fetch.");
1761
        assert!(bam.records().count() == 6);
1762

1763
        //all on a tid - by &[u8]
1764
        bam.fetch(sq_1).expect("Expected successful fetch.");
1765
        assert!(bam.records().count() == 6);
1766
        //all on a tid - by str
1767
        bam.fetch(str::from_utf8(sq_1).unwrap())
1768
            .expect("Expected successful fetch.");
1769
        assert!(bam.records().count() == 6);
1770

1771
        //all reads
1772
        bam.fetch(FetchDefinition::All)
1773
            .expect("Expected successful fetch.");
1774
        assert!(bam.records().count() == 6);
1775

1776
        //all reads
1777
        bam.fetch(".").expect("Expected successful fetch.");
1778
        assert!(bam.records().count() == 6);
1779

1780
        //all unmapped
1781
        bam.fetch(FetchDefinition::Unmapped)
1782
            .expect("Expected successful fetch.");
1783
        assert_eq!(bam.records().count(), 1); // expect one 'truncade record' Record.
1784

1785
        bam.fetch("*").expect("Expected successful fetch.");
1786
        assert_eq!(bam.records().count(), 1); // expect one 'truncade record' Record.
1787
    }
1788

1789
    #[test]
1790
    fn test_read_indexed() {
1791
        let bam = IndexedReader::from_path("test/test.bam").expect("Expected valid index.");
1792
        _test_read_indexed_common(bam);
1793
    }
1794

1795
    #[test]
1796
    fn test_read_indexed_different_index_name() {
1797
        let bam = IndexedReader::from_path_and_index(
1798
            &"test/test_different_index_name.bam",
1799
            &"test/test.bam.bai",
1800
        )
1801
        .expect("Expected valid index.");
1802
        _test_read_indexed_common(bam);
1803
    }
1804

1805
    #[test]
1806
    fn test_set_record() {
1807
        let (names, _, seqs, quals, cigars) = gold();
1808

1809
        let mut rec = record::Record::new();
1810
        rec.set_reverse();
1811
        rec.set(names[0], Some(&cigars[0]), seqs[0], quals[0]);
1812
        // note: this segfaults if you push_aux() before set()
1813
        //       because set() obliterates aux
1814
        rec.push_aux(b"NM", Aux::I32(15)).unwrap();
1815

1816
        assert_eq!(rec.qname(), names[0]);
1817
        assert_eq!(*rec.cigar(), cigars[0]);
1818
        assert_eq!(rec.seq().as_bytes(), seqs[0]);
1819
        assert_eq!(rec.qual(), quals[0]);
1820
        assert!(rec.is_reverse());
1821
        assert_eq!(rec.aux(b"NM").unwrap(), Aux::I32(15));
1822
    }
1823

1824
    #[test]
1825
    fn test_set_repeated() {
1826
        let mut rec = Record::new();
1827
        rec.set(
1828
            b"123",
1829
            Some(&CigarString(vec![Cigar::Match(3)])),
1830
            b"AAA",
1831
            b"III",
1832
        );
1833
        rec.push_aux(b"AS", Aux::I32(12345)).unwrap();
1834
        assert_eq!(rec.qname(), b"123");
1835
        assert_eq!(rec.seq().as_bytes(), b"AAA");
1836
        assert_eq!(rec.qual(), b"III");
1837
        assert_eq!(rec.aux(b"AS").unwrap(), Aux::I32(12345));
1838

1839
        rec.set(
1840
            b"1234",
1841
            Some(&CigarString(vec![Cigar::SoftClip(1), Cigar::Match(3)])),
1842
            b"AAAA",
1843
            b"IIII",
1844
        );
1845
        assert_eq!(rec.qname(), b"1234");
1846
        assert_eq!(rec.seq().as_bytes(), b"AAAA");
1847
        assert_eq!(rec.qual(), b"IIII");
1848
        assert_eq!(rec.aux(b"AS").unwrap(), Aux::I32(12345));
1849

1850
        rec.set(
1851
            b"12",
1852
            Some(&CigarString(vec![Cigar::Match(2)])),
1853
            b"AA",
1854
            b"II",
1855
        );
1856
        assert_eq!(rec.qname(), b"12");
1857
        assert_eq!(rec.seq().as_bytes(), b"AA");
1858
        assert_eq!(rec.qual(), b"II");
1859
        assert_eq!(rec.aux(b"AS").unwrap(), Aux::I32(12345));
1860
    }
1861

1862
    #[test]
1863
    fn test_set_qname() {
1864
        let (names, _, seqs, quals, cigars) = gold();
1865

1866
        assert!(names[0] != names[1]);
1867

1868
        for i in 0..names.len() {
1869
            let mut rec = record::Record::new();
1870
            rec.set(names[i], Some(&cigars[i]), seqs[i], quals[i]);
1871
            rec.push_aux(b"NM", Aux::I32(15)).unwrap();
1872

1873
            assert_eq!(rec.qname(), names[i]);
1874
            assert_eq!(*rec.cigar(), cigars[i]);
1875
            assert_eq!(rec.seq().as_bytes(), seqs[i]);
1876
            assert_eq!(rec.qual(), quals[i]);
1877
            assert_eq!(rec.aux(b"NM").unwrap(), Aux::I32(15));
1878

1879
            // Equal length qname
1880
            assert!(rec.qname()[0] != b'X');
1881
            rec.set_qname(b"X");
1882
            assert_eq!(rec.qname(), b"X");
1883

1884
            // Longer qname
1885
            let mut longer_name = names[i].to_owned().clone();
1886
            let extension = b"BuffaloBUffaloBUFFaloBUFFAloBUFFALoBUFFALO";
1887
            longer_name.extend(extension.iter());
1888
            rec.set_qname(&longer_name);
1889

1890
            assert_eq!(rec.qname(), longer_name.as_slice());
1891
            assert_eq!(*rec.cigar(), cigars[i]);
1892
            assert_eq!(rec.seq().as_bytes(), seqs[i]);
1893
            assert_eq!(rec.qual(), quals[i]);
1894
            assert_eq!(rec.aux(b"NM").unwrap(), Aux::I32(15));
1895

1896
            // Shorter qname
1897
            let shorter_name = b"42";
1898
            rec.set_qname(shorter_name);
1899

1900
            assert_eq!(rec.qname(), shorter_name);
1901
            assert_eq!(*rec.cigar(), cigars[i]);
1902
            assert_eq!(rec.seq().as_bytes(), seqs[i]);
1903
            assert_eq!(rec.qual(), quals[i]);
1904
            assert_eq!(rec.aux(b"NM").unwrap(), Aux::I32(15));
1905

1906
            // Zero-length qname
1907
            rec.set_qname(b"");
1908

1909
            assert_eq!(rec.qname(), b"");
1910
            assert_eq!(*rec.cigar(), cigars[i]);
1911
            assert_eq!(rec.seq().as_bytes(), seqs[i]);
1912
            assert_eq!(rec.qual(), quals[i]);
1913
            assert_eq!(rec.aux(b"NM").unwrap(), Aux::I32(15));
1914
        }
1915
    }
1916

1917
    #[test]
1918
    fn test_set_qname2() {
1919
        let mut _header = Header::new();
1920
        _header.push_record(
1921
            HeaderRecord::new(b"SQ")
1922
                .push_tag(b"SN", "1")
1923
                .push_tag(b"LN", 10000000),
1924
        );
1925
        let header = HeaderView::from_header(&_header);
1926

1927
        let line =
1928
            b"blah1        0        1        1        255        1M        *        0        0        A        F        CB:Z:AAAA-1        UR:Z:AAAA        UB:Z:AAAA        GX:Z:G1        xf:i:1        fx:Z:G1\tli:i:0\ttf:Z:cC";
1929

1930
        let mut rec = Record::from_sam(&header, line).unwrap();
1931
        assert_eq!(rec.qname(), b"blah1");
1932
        rec.set_qname(b"r0");
1933
        assert_eq!(rec.qname(), b"r0");
1934
    }
1935

1936
    #[test]
1937
    fn test_set_cigar() {
1938
        let (names, _, seqs, quals, cigars) = gold();
1939

1940
        assert!(names[0] != names[1]);
1941

1942
        for i in 0..names.len() {
1943
            let mut rec = record::Record::new();
1944
            rec.set(names[i], Some(&cigars[i]), seqs[i], quals[i]);
1945
            rec.push_aux(b"NM", Aux::I32(15)).unwrap();
1946

1947
            assert_eq!(rec.qname(), names[i]);
1948
            assert_eq!(*rec.cigar(), cigars[i]);
1949
            assert_eq!(rec.seq().as_bytes(), seqs[i]);
1950
            assert_eq!(rec.qual(), quals[i]);
1951
            assert_eq!(rec.aux(b"NM").unwrap(), Aux::I32(15));
1952

1953
            // boring cigar
1954
            let new_cigar = CigarString(vec![Cigar::Match(rec.seq_len() as u32)]);
1955
            assert_ne!(*rec.cigar(), new_cigar);
1956
            rec.set_cigar(Some(&new_cigar));
1957
            assert_eq!(*rec.cigar(), new_cigar);
1958

1959
            assert_eq!(rec.qname(), names[i]);
1960
            assert_eq!(rec.seq().as_bytes(), seqs[i]);
1961
            assert_eq!(rec.qual(), quals[i]);
1962
            assert_eq!(rec.aux(b"NM").unwrap(), Aux::I32(15));
1963

1964
            // bizarre cigar
1965
            let new_cigar = (0..rec.seq_len())
1966
                .map(|i| {
1967
                    if i % 2 == 0 {
1968
                        Cigar::Match(1)
1969
                    } else {
1970
                        Cigar::Ins(1)
1971
                    }
1972
                })
1973
                .collect::<Vec<_>>();
1974
            let new_cigar = CigarString(new_cigar);
1975
            assert_ne!(*rec.cigar(), new_cigar);
1976
            rec.set_cigar(Some(&new_cigar));
1977
            assert_eq!(*rec.cigar(), new_cigar);
1978

1979
            assert_eq!(rec.qname(), names[i]);
1980
            assert_eq!(rec.seq().as_bytes(), seqs[i]);
1981
            assert_eq!(rec.qual(), quals[i]);
1982
            assert_eq!(rec.aux(b"NM").unwrap(), Aux::I32(15));
1983

1984
            // empty cigar
1985
            let new_cigar = CigarString(Vec::new());
1986
            assert_ne!(*rec.cigar(), new_cigar);
1987
            rec.set_cigar(None);
1988
            assert_eq!(*rec.cigar(), new_cigar);
1989

1990
            assert_eq!(rec.qname(), names[i]);
1991
            assert_eq!(rec.seq().as_bytes(), seqs[i]);
1992
            assert_eq!(rec.qual(), quals[i]);
1993
            assert_eq!(rec.aux(b"NM").unwrap(), Aux::I32(15));
1994
        }
1995
    }
1996

1997
    #[test]
1998
    fn test_remove_aux() {
1999
        let mut bam = Reader::from_path(Path::new("test/test.bam")).expect("Error opening file.");
2000

2001
        for record in bam.records() {
2002
            let mut rec = record.expect("Expected valid record");
2003

2004
            if rec.aux(b"XS").is_ok() {
2005
                rec.remove_aux(b"XS").unwrap();
2006
            }
2007

2008
            if rec.aux(b"YT").is_ok() {
2009
                rec.remove_aux(b"YT").unwrap();
2010
            }
2011

2012
            assert!(rec.remove_aux(b"ab").is_err());
2013

2014
            assert_eq!(rec.aux(b"XS"), Err(Error::BamAuxTagNotFound));
2015
            assert_eq!(rec.aux(b"YT"), Err(Error::BamAuxTagNotFound));
2016
        }
2017
    }
2018

2019
    #[test]
2020
    fn test_write() {
2021
        let (names, _, seqs, quals, cigars) = gold();
2022

2023
        let tmp = tempfile::Builder::new()
2024
            .prefix("rust-htslib")
2025
            .tempdir()
2026
            .expect("Cannot create temp dir");
2027
        let bampath = tmp.path().join("test.bam");
2028
        println!("{:?}", bampath);
2029
        {
2030
            let mut bam = Writer::from_path(
2031
                &bampath,
2032
                Header::new().push_record(
2033
                    HeaderRecord::new(b"SQ")
2034
                        .push_tag(b"SN", "chr1")
2035
                        .push_tag(b"LN", 15072423),
2036
                ),
2037
                Format::Bam,
2038
            )
2039
            .expect("Error opening file.");
2040

2041
            for i in 0..names.len() {
2042
                let mut rec = record::Record::new();
2043
                rec.set(names[i], Some(&cigars[i]), seqs[i], quals[i]);
2044
                rec.push_aux(b"NM", Aux::I32(15)).unwrap();
2045

2046
                bam.write(&rec).expect("Failed to write record.");
2047
            }
2048
        }
2049

2050
        {
2051
            let mut bam = Reader::from_path(bampath).expect("Error opening file.");
2052

2053
            for i in 0..names.len() {
2054
                let mut rec = record::Record::new();
2055
                if let Some(r) = bam.read(&mut rec) {
2056
                    r.expect("Failed to read record.");
2057
                };
2058

2059
                assert_eq!(rec.qname(), names[i]);
2060
                assert_eq!(*rec.cigar(), cigars[i]);
2061
                assert_eq!(rec.seq().as_bytes(), seqs[i]);
2062
                assert_eq!(rec.qual(), quals[i]);
2063
                assert_eq!(rec.aux(b"NM").unwrap(), Aux::I32(15));
2064
            }
2065
        }
2066

2067
        tmp.close().expect("Failed to delete temp dir");
2068
    }
2069

2070
    #[test]
2071
    fn test_write_threaded() {
2072
        let (names, _, seqs, quals, cigars) = gold();
2073

2074
        let tmp = tempfile::Builder::new()
2075
            .prefix("rust-htslib")
2076
            .tempdir()
2077
            .expect("Cannot create temp dir");
2078
        let bampath = tmp.path().join("test.bam");
2079
        println!("{:?}", bampath);
2080
        {
2081
            let mut bam = Writer::from_path(
2082
                &bampath,
2083
                Header::new().push_record(
2084
                    HeaderRecord::new(b"SQ")
2085
                        .push_tag(b"SN", "chr1")
2086
                        .push_tag(b"LN", 15072423),
2087
                ),
2088
                Format::Bam,
2089
            )
2090
            .expect("Error opening file.");
2091
            bam.set_threads(4).unwrap();
2092

2093
            for i in 0..10000 {
2094
                let mut rec = record::Record::new();
2095
                let idx = i % names.len();
2096
                rec.set(names[idx], Some(&cigars[idx]), seqs[idx], quals[idx]);
2097
                rec.push_aux(b"NM", Aux::I32(15)).unwrap();
2098
                rec.set_pos(i as i64);
2099

2100
                bam.write(&rec).expect("Failed to write record.");
2101
            }
2102
        }
2103

2104
        {
2105
            let mut bam = Reader::from_path(bampath).expect("Error opening file.");
2106

2107
            for (i, _rec) in bam.records().enumerate() {
2108
                let idx = i % names.len();
2109

2110
                let rec = _rec.expect("Failed to read record.");
2111

2112
                assert_eq!(rec.pos(), i as i64);
2113
                assert_eq!(rec.qname(), names[idx]);
2114
                assert_eq!(*rec.cigar(), cigars[idx]);
2115
                assert_eq!(rec.seq().as_bytes(), seqs[idx]);
2116
                assert_eq!(rec.qual(), quals[idx]);
2117
                assert_eq!(rec.aux(b"NM").unwrap(), Aux::I32(15));
2118
            }
2119
        }
2120

2121
        tmp.close().expect("Failed to delete temp dir");
2122
    }
2123

2124
    #[test]
2125
    fn test_write_shared_tpool() {
2126
        let (names, _, seqs, quals, cigars) = gold();
2127

2128
        let tmp = tempfile::Builder::new()
2129
            .prefix("rust-htslib")
2130
            .tempdir()
2131
            .expect("Cannot create temp dir");
2132
        let bampath1 = tmp.path().join("test1.bam");
2133
        let bampath2 = tmp.path().join("test2.bam");
2134

2135
        {
2136
            let (mut bam1, mut bam2) = {
2137
                let pool = crate::tpool::ThreadPool::new(4).unwrap();
2138

2139
                let mut bam1 = Writer::from_path(
2140
                    &bampath1,
2141
                    Header::new().push_record(
2142
                        HeaderRecord::new(b"SQ")
2143
                            .push_tag(b"SN", "chr1")
2144
                            .push_tag(b"LN", 15072423),
2145
                    ),
2146
                    Format::Bam,
2147
                )
2148
                .expect("Error opening file.");
2149

2150
                let mut bam2 = Writer::from_path(
2151
                    &bampath2,
2152
                    Header::new().push_record(
2153
                        HeaderRecord::new(b"SQ")
2154
                            .push_tag(b"SN", "chr1")
2155
                            .push_tag(b"LN", 15072423),
2156
                    ),
2157
                    Format::Bam,
2158
                )
2159
                .expect("Error opening file.");
2160

2161
                bam1.set_thread_pool(&pool).unwrap();
2162
                bam2.set_thread_pool(&pool).unwrap();
2163
                (bam1, bam2)
2164
            };
2165

2166
            for i in 0..10000 {
2167
                let mut rec = record::Record::new();
2168
                let idx = i % names.len();
2169
                rec.set(names[idx], Some(&cigars[idx]), seqs[idx], quals[idx]);
2170
                rec.push_aux(b"NM", Aux::I32(15)).unwrap();
2171
                rec.set_pos(i as i64);
2172

2173
                bam1.write(&rec).expect("Failed to write record.");
2174
                bam2.write(&rec).expect("Failed to write record.");
2175
            }
2176
        }
2177

2178
        {
2179
            let pool = crate::tpool::ThreadPool::new(2).unwrap();
2180

2181
            for p in [bampath1, bampath2] {
2182
                let mut bam = Reader::from_path(p).expect("Error opening file.");
2183
                bam.set_thread_pool(&pool).unwrap();
2184

2185
                for (i, _rec) in bam.iter_chunk(None, None).enumerate() {
2186
                    let idx = i % names.len();
2187

2188
                    let rec = _rec.expect("Failed to read record.");
2189

2190
                    assert_eq!(rec.pos(), i as i64);
2191
                    assert_eq!(rec.qname(), names[idx]);
2192
                    assert_eq!(*rec.cigar(), cigars[idx]);
2193
                    assert_eq!(rec.seq().as_bytes(), seqs[idx]);
2194
                    assert_eq!(rec.qual(), quals[idx]);
2195
                    assert_eq!(rec.aux(b"NM").unwrap(), Aux::I32(15));
2196
                }
2197
            }
2198
        }
2199

2200
        tmp.close().expect("Failed to delete temp dir");
2201
    }
2202

2203
    #[test]
2204
    fn test_copy_template() {
2205
        // Verify that BAM headers are transmitted correctly when using an existing BAM as a
2206
        // template for headers.
2207

2208
        let tmp = tempfile::Builder::new()
2209
            .prefix("rust-htslib")
2210
            .tempdir()
2211
            .expect("Cannot create temp dir");
2212
        let bampath = tmp.path().join("test.bam");
2213
        println!("{:?}", bampath);
2214

2215
        let mut input_bam = Reader::from_path("test/test.bam").expect("Error opening file.");
2216

2217
        {
2218
            let mut bam = Writer::from_path(
2219
                &bampath,
2220
                &Header::from_template(input_bam.header()),
2221
                Format::Bam,
2222
            )
2223
            .expect("Error opening file.");
2224

2225
            for rec in input_bam.records() {
2226
                bam.write(&rec.unwrap()).expect("Failed to write record.");
2227
            }
2228
        }
2229

2230
        {
2231
            let copy_bam = Reader::from_path(bampath).expect("Error opening file.");
2232

2233
            // Verify that the header came across correctly
2234
            assert_eq!(input_bam.header().as_bytes(), copy_bam.header().as_bytes());
2235
        }
2236

2237
        tmp.close().expect("Failed to delete temp dir");
2238
    }
2239

2240
    #[test]
2241
    fn test_pileup() {
2242
        let (_, _, seqs, quals, _) = gold();
2243

2244
        let mut bam = Reader::from_path("test/test.bam").expect("Error opening file.");
2245
        let pileups = bam.pileup();
2246
        for pileup in pileups.take(26) {
2247
            let _pileup = pileup.expect("Expected successful pileup.");
2248
            let pos = _pileup.pos() as usize;
2249
            assert_eq!(_pileup.depth(), 6);
2250
            assert!(_pileup.tid() == 0);
2251
            for (i, a) in _pileup.alignments().enumerate() {
2252
                assert_eq!(a.indel(), pileup::Indel::None);
2253
                let qpos = a.qpos().unwrap();
2254
                assert_eq!(qpos, pos - 1);
2255
                assert_eq!(a.record().seq()[qpos], seqs[i][qpos]);
2256
                assert_eq!(a.record().qual()[qpos], quals[i][qpos] - 33);
2257
            }
2258
        }
2259
    }
2260

2261
    #[test]
2262
    fn test_idx_pileup() {
2263
        let mut bam = IndexedReader::from_path("test/test.bam").expect("Error opening file.");
2264
        // read without fetch
2265
        for pileup in bam.pileup() {
2266
            pileup.unwrap();
2267
        }
2268
        // go back again
2269
        let tid = bam.header().tid(b"CHROMOSOME_I").unwrap();
2270
        bam.fetch((tid, 0, 5)).unwrap();
2271
        for p in bam.pileup() {
2272
            println!("{}", p.unwrap().pos())
2273
        }
2274
    }
2275

2276
    #[test]
2277
    fn parse_from_sam() {
2278
        use std::fs::File;
2279
        use std::io::Read;
2280

2281
        let bamfile = "./test/bam2sam_test.bam";
2282
        let samfile = "./test/bam2sam_expected.sam";
2283

2284
        // Load BAM file:
2285
        let mut rdr = Reader::from_path(bamfile).unwrap();
2286
        let bam_recs: Vec<Record> = rdr.records().map(|v| v.unwrap()).collect();
2287

2288
        let mut sam = Vec::new();
2289
        assert!(File::open(samfile).unwrap().read_to_end(&mut sam).is_ok());
2290

2291
        let sam_recs: Vec<Record> = sam
2292
            .split(|x| *x == b'\n')
2293
            .filter(|x| !x.is_empty() && x[0] != b'@')
2294
            .map(|line| Record::from_sam(rdr.header(), line).unwrap())
2295
            .collect();
2296

2297
        for (b1, s1) in bam_recs.iter().zip(sam_recs.iter()) {
2298
            assert!(b1 == s1);
2299
        }
2300
    }
2301

2302
    #[test]
2303
    fn test_cigar_modes() {
2304
        // test the cached and uncached ways of getting the cigar string.
2305

2306
        let (_, _, _, _, cigars) = gold();
2307
        let mut bam = Reader::from_path(Path::new("test/test.bam")).expect("Error opening file.");
2308

2309
        for (i, record) in bam.records().enumerate() {
2310
            let rec = record.expect("Expected valid record");
2311

2312
            let cigar = rec.cigar();
2313
            assert_eq!(*cigar, cigars[i]);
2314
        }
2315

2316
        for (i, record) in bam.records().enumerate() {
2317
            let mut rec = record.expect("Expected valid record");
2318
            rec.cache_cigar();
2319

2320
            let cigar = rec.cigar_cached().unwrap();
2321
            assert_eq!(**cigar, cigars[i]);
2322

2323
            let cigar = rec.cigar();
2324
            assert_eq!(*cigar, cigars[i]);
2325
        }
2326
    }
2327

2328
    #[test]
2329
    fn test_read_cram() {
2330
        let cram_path = "./test/test_cram.cram";
2331
        let bam_path = "./test/test_cram.bam";
2332
        let ref_path = "./test/test_cram.fa";
2333

2334
        // Load CRAM file, records
2335
        let mut cram_reader = Reader::from_path(cram_path).unwrap();
2336
        cram_reader.set_reference(ref_path).unwrap();
2337
        let cram_records: Vec<Record> = cram_reader.records().map(|v| v.unwrap()).collect();
2338

2339
        // Load BAM file, records
2340
        let mut bam_reader = Reader::from_path(bam_path).unwrap();
2341
        let bam_records: Vec<Record> = bam_reader.records().map(|v| v.unwrap()).collect();
2342

2343
        compare_inner_bam_cram_records(&cram_records, &bam_records);
2344
    }
2345

2346
    #[test]
2347
    fn test_write_cram() {
2348
        // BAM file, records
2349
        let bam_path = "./test/test_cram.bam";
2350
        let ref_path = "./test/test_cram.fa";
2351
        let mut bam_reader = Reader::from_path(bam_path).unwrap();
2352
        let bam_records: Vec<Record> = bam_reader.records().map(|v| v.unwrap()).collect();
2353

2354
        // New CRAM file
2355
        let tmp = tempfile::Builder::new()
2356
            .prefix("rust-htslib")
2357
            .tempdir()
2358
            .expect("Cannot create temp dir");
2359
        let cram_path = tmp.path().join("test.cram");
2360

2361
        // Write BAM records to new CRAM file
2362
        {
2363
            let mut header = Header::new();
2364
            header.push_record(
2365
                HeaderRecord::new(b"HD")
2366
                    .push_tag(b"VN", "1.5")
2367
                    .push_tag(b"SO", "coordinate"),
2368
            );
2369
            header.push_record(
2370
                HeaderRecord::new(b"SQ")
2371
                    .push_tag(b"SN", "chr1")
2372
                    .push_tag(b"LN", 120)
2373
                    .push_tag(b"M5", "20a9a0fb770814e6c5e49946750f9724")
2374
                    .push_tag(b"UR", "test/test_cram.fa"),
2375
            );
2376
            header.push_record(
2377
                HeaderRecord::new(b"SQ")
2378
                    .push_tag(b"SN", "chr2")
2379
                    .push_tag(b"LN", 120)
2380
                    .push_tag(b"M5", "7a2006ccca94ea92b6dae5997e1b0d70")
2381
                    .push_tag(b"UR", "test/test_cram.fa"),
2382
            );
2383
            header.push_record(
2384
                HeaderRecord::new(b"SQ")
2385
                    .push_tag(b"SN", "chr3")
2386
                    .push_tag(b"LN", 120)
2387
                    .push_tag(b"M5", "a66b336bfe3ee8801c744c9545c87e24")
2388
                    .push_tag(b"UR", "test/test_cram.fa"),
2389
            );
2390

2391
            let mut cram_writer = Writer::from_path(&cram_path, &header, Format::Cram)
2392
                .expect("Error opening CRAM file.");
2393
            cram_writer.set_reference(ref_path).unwrap();
2394

2395
            // Write BAM records to CRAM file
2396
            for rec in bam_records.iter() {
2397
                cram_writer
2398
                    .write(rec)
2399
                    .expect("Faied to write record to CRAM.");
2400
            }
2401
        }
2402

2403
        // Compare written CRAM records with BAM records
2404
        {
2405
            // Load written CRAM file
2406
            let mut cram_reader = Reader::from_path(cram_path).unwrap();
2407
            cram_reader.set_reference(ref_path).unwrap();
2408
            let cram_records: Vec<Record> = cram_reader.records().map(|v| v.unwrap()).collect();
2409

2410
            // Compare CRAM records to BAM records
2411
            compare_inner_bam_cram_records(&cram_records, &bam_records);
2412
        }
2413

2414
        tmp.close().expect("Failed to delete temp dir");
2415
    }
2416

2417
    #[test]
2418
    fn test_compression_level_conversion() {
2419
        // predefined compression levels
2420
        assert_eq!(CompressionLevel::Uncompressed.convert().unwrap(), 0);
2421
        assert_eq!(CompressionLevel::Fastest.convert().unwrap(), 1);
2422
        assert_eq!(CompressionLevel::Maximum.convert().unwrap(), 9);
2423

2424
        // numeric compression levels
2425
        for level in 0..=9 {
2426
            assert_eq!(CompressionLevel::Level(level).convert().unwrap(), level);
2427
        }
2428
        // invalid levels
2429
        assert!(CompressionLevel::Level(10).convert().is_err());
2430
    }
2431

2432
    #[test]
2433
    fn test_write_compression() {
2434
        let tmp = tempfile::Builder::new()
2435
            .prefix("rust-htslib")
2436
            .tempdir()
2437
            .expect("Cannot create temp dir");
2438
        let input_bam_path = "test/test.bam";
2439

2440
        // test levels with decreasing compression factor
2441
        let levels_to_test = vec![
2442
            CompressionLevel::Maximum,
2443
            CompressionLevel::Level(6),
2444
            CompressionLevel::Fastest,
2445
            CompressionLevel::Uncompressed,
2446
        ];
2447
        let file_sizes: Vec<_> = levels_to_test
2448
            .iter()
2449
            .map(|level| {
2450
                let output_bam_path = tmp.path().join("test.bam");
2451
                {
2452
                    let mut reader = Reader::from_path(input_bam_path).unwrap();
2453
                    let header = Header::from_template(reader.header());
2454
                    let mut writer =
2455
                        Writer::from_path(&output_bam_path, &header, Format::Bam).unwrap();
2456
                    writer.set_compression_level(*level).unwrap();
2457
                    for record in reader.records() {
2458
                        let r = record.unwrap();
2459
                        writer.write(&r).unwrap();
2460
                    }
2461
                }
2462
                fs::metadata(output_bam_path).unwrap().len()
2463
            })
2464
            .collect();
2465

2466
        // check that out BAM file sizes are in decreasing order, in line with the expected compression factor
2467
        println!("testing compression leves: {:?}", levels_to_test);
2468
        println!("got compressed sizes: {:?}", file_sizes);
2469

2470
        // libdeflate comes out with a slightly bigger file on Max compression
2471
        // than on Level(6), so skip that check
2472
        #[cfg(feature = "libdeflate")]
2473
        assert!(file_sizes[1..].windows(2).all(|size| size[0] <= size[1]));
2474

2475
        #[cfg(not(feature = "libdeflate"))]
2476
        assert!(file_sizes.windows(2).all(|size| size[0] <= size[1]));
2477

2478
        tmp.close().expect("Failed to delete temp dir");
2479
    }
2480

2481
    #[test]
2482
    fn test_bam_fails_on_vcf() {
2483
        let bam_path = "./test/test_left.vcf";
2484
        let bam_reader = Reader::from_path(bam_path);
2485
        assert!(bam_reader.is_err());
2486
    }
2487

2488
    #[test]
2489
    fn test_indexde_bam_fails_on_vcf() {
2490
        let bam_path = "./test/test_left.vcf";
2491
        let bam_reader = IndexedReader::from_path(bam_path);
2492
        assert!(bam_reader.is_err());
2493
    }
2494

2495
    #[test]
2496
    fn test_bam_fails_on_toml() {
2497
        let bam_path = "./Cargo.toml";
2498
        let bam_reader = Reader::from_path(bam_path);
2499
        assert!(bam_reader.is_err());
2500
    }
2501

2502
    #[test]
2503
    fn test_sam_writer_example() {
2504
        fn from_bam_with_filter<F>(bamfile: &str, samfile: &str, f: F) -> bool
2505
        where
2506
            F: Fn(&record::Record) -> Option<bool>,
2507
        {
2508
            let mut bam_reader = Reader::from_path(bamfile).unwrap(); // internal functions, just unwrap
2509
            let header = header::Header::from_template(bam_reader.header());
2510
            let mut sam_writer = Writer::from_path(samfile, &header, Format::Sam).unwrap();
2511
            for record in bam_reader.records() {
2512
                if record.is_err() {
2513
                    return false;
2514
                }
2515
                let parsed = record.unwrap();
2516
                match f(&parsed) {
2517
                    None => return true,
2518
                    Some(false) => {}
2519
                    Some(true) => {
2520
                        if sam_writer.write(&parsed).is_err() {
2521
                            return false;
2522
                        }
2523
                    }
2524
                }
2525
            }
2526
            true
2527
        }
2528
        use std::fs::File;
2529
        use std::io::Read;
2530
        let bamfile = "./test/bam2sam_test.bam";
2531
        let samfile = "./test/bam2sam_out.sam";
2532
        let expectedfile = "./test/bam2sam_expected.sam";
2533
        let result = from_bam_with_filter(bamfile, samfile, |_| Some(true));
2534
        assert!(result);
2535
        let mut expected = Vec::new();
2536
        let mut written = Vec::new();
2537
        assert!(File::open(expectedfile)
2538
            .unwrap()
2539
            .read_to_end(&mut expected)
2540
            .is_ok());
2541
        assert!(File::open(samfile)
2542
            .unwrap()
2543
            .read_to_end(&mut written)
2544
            .is_ok());
2545
        assert_eq!(expected, written);
2546
    }
2547

2548
    // #[cfg(feature = "curl")]
2549
    // #[test]
2550
    // fn test_http_connect() {
2551
    //     let url: Url = Url::parse(
2552
    //         "https://raw.githubusercontent.com/brainstorm/tiny-test-data/master/wgs/mt.bam",
2553
    //     )
2554
    //     .unwrap();
2555
    //     let r = Reader::from_url(&url);
2556
    //     println!("{:#?}", r);
2557
    //     let r = r.unwrap();
2558

2559
    //     assert_eq!(r.header().target_names()[0], b"chr1");
2560
    // }
2561

2562
    #[test]
2563
    fn test_rc_records() {
2564
        let (names, flags, seqs, quals, cigars) = gold();
2565
        let mut bam = Reader::from_path(Path::new("test/test.bam")).expect("Error opening file.");
2566
        let del_len = [1, 1, 1, 1, 1, 100000];
2567

2568
        for (i, record) in bam.rc_records().enumerate() {
2569
            //let rec = record.expect("Expected valid record");
2570
            let rec = record.unwrap();
2571
            println!("{}", str::from_utf8(rec.qname()).ok().unwrap());
2572
            assert_eq!(rec.qname(), names[i]);
2573
            assert_eq!(rec.flags(), flags[i]);
2574
            assert_eq!(rec.seq().as_bytes(), seqs[i]);
2575

2576
            let cigar = rec.cigar();
2577
            assert_eq!(*cigar, cigars[i]);
2578

2579
            let end_pos = cigar.end_pos();
2580
            assert_eq!(end_pos, rec.pos() + 100 + del_len[i]);
2581
            assert_eq!(
2582
                cigar
2583
                    .read_pos(end_pos as u32 - 10, false, false)
2584
                    .unwrap()
2585
                    .unwrap(),
2586
                90
2587
            );
2588
            assert_eq!(
2589
                cigar
2590
                    .read_pos(rec.pos() as u32 + 20, false, false)
2591
                    .unwrap()
2592
                    .unwrap(),
2593
                20
2594
            );
2595
            assert_eq!(cigar.read_pos(4000000, false, false).unwrap(), None);
2596
            // fix qual offset
2597
            let qual: Vec<u8> = quals[i].iter().map(|&q| q - 33).collect();
2598
            assert_eq!(rec.qual(), &qual[..]);
2599
        }
2600
    }
2601

2602
    #[test]
2603
    fn test_aux_arrays() {
2604
        let bam_header = Header::new();
2605
        let mut test_record = Record::from_sam(
2606
            &HeaderView::from_header(&bam_header),
2607
            "ali1\t4\t*\t0\t0\t*\t*\t0\t0\tACGT\tFFFF".as_bytes(),
2608
        )
2609
        .unwrap();
2610

2611
        let array_i8: Vec<i8> = vec![i8::MIN, -1, 0, 1, i8::MAX];
2612
        let array_u8: Vec<u8> = vec![u8::MIN, 0, 1, u8::MAX];
2613
        let array_i16: Vec<i16> = vec![i16::MIN, -1, 0, 1, i16::MAX];
2614
        let array_u16: Vec<u16> = vec![u16::MIN, 0, 1, u16::MAX];
2615
        let array_i32: Vec<i32> = vec![i32::MIN, -1, 0, 1, i32::MAX];
2616
        let array_u32: Vec<u32> = vec![u32::MIN, 0, 1, u32::MAX];
2617
        let array_f32: Vec<f32> = vec![f32::MIN, 0.0, -0.0, 0.1, 0.99, f32::MAX];
2618

2619
        test_record
2620
            .push_aux(b"XA", Aux::ArrayI8((&array_i8).into()))
2621
            .unwrap();
2622
        test_record
2623
            .push_aux(b"XB", Aux::ArrayU8((&array_u8).into()))
2624
            .unwrap();
2625
        test_record
2626
            .push_aux(b"XC", Aux::ArrayI16((&array_i16).into()))
2627
            .unwrap();
2628
        test_record
2629
            .push_aux(b"XD", Aux::ArrayU16((&array_u16).into()))
2630
            .unwrap();
2631
        test_record
2632
            .push_aux(b"XE", Aux::ArrayI32((&array_i32).into()))
2633
            .unwrap();
2634
        test_record
2635
            .push_aux(b"XF", Aux::ArrayU32((&array_u32).into()))
2636
            .unwrap();
2637
        test_record
2638
            .push_aux(b"XG", Aux::ArrayFloat((&array_f32).into()))
2639
            .unwrap();
2640

2641
        {
2642
            let tag = b"XA";
2643
            if let Ok(Aux::ArrayI8(array)) = test_record.aux(tag) {
2644
                // Retrieve aux array
2645
                let aux_array_content = array.iter().collect::<Vec<_>>();
2646
                assert_eq!(aux_array_content, array_i8);
2647

2648
                // Copy the stored aux array to another record
2649
                {
2650
                    let mut copy_test_record = test_record.clone();
2651

2652
                    // Pushing a field with an existing tag should fail
2653
                    assert!(copy_test_record.push_aux(tag, Aux::I8(3)).is_err());
2654

2655
                    // Remove aux array from target record
2656
                    copy_test_record.remove_aux(tag).unwrap();
2657
                    assert!(copy_test_record.aux(tag).is_err());
2658

2659
                    // Copy array to target record
2660
                    let src_aux = test_record.aux(tag).unwrap();
2661
                    assert!(copy_test_record.push_aux(tag, src_aux).is_ok());
2662
                    if let Ok(Aux::ArrayI8(array)) = copy_test_record.aux(tag) {
2663
                        let aux_array_content_copied = array.iter().collect::<Vec<_>>();
2664
                        assert_eq!(aux_array_content_copied, array_i8);
2665
                    } else {
2666
                        panic!("Aux tag not found");
2667
                    }
2668
                }
2669
            } else {
2670
                panic!("Aux tag not found");
2671
            }
2672
        }
2673

2674
        {
2675
            let tag = b"XB";
2676
            if let Ok(Aux::ArrayU8(array)) = test_record.aux(tag) {
2677
                // Retrieve aux array
2678
                let aux_array_content = array.iter().collect::<Vec<_>>();
2679
                assert_eq!(aux_array_content, array_u8);
2680

2681
                // Copy the stored aux array to another record
2682
                {
2683
                    let mut copy_test_record = test_record.clone();
2684

2685
                    // Pushing a field with an existing tag should fail
2686
                    assert!(copy_test_record.push_aux(tag, Aux::U8(3)).is_err());
2687

2688
                    // Remove aux array from target record
2689
                    copy_test_record.remove_aux(tag).unwrap();
2690
                    assert!(copy_test_record.aux(tag).is_err());
2691

2692
                    // Copy array to target record
2693
                    let src_aux = test_record.aux(tag).unwrap();
2694
                    assert!(copy_test_record.push_aux(tag, src_aux).is_ok());
2695
                    if let Ok(Aux::ArrayU8(array)) = copy_test_record.aux(tag) {
2696
                        let aux_array_content_copied = array.iter().collect::<Vec<_>>();
2697
                        assert_eq!(aux_array_content_copied, array_u8);
2698
                    } else {
2699
                        panic!("Aux tag not found");
2700
                    }
2701
                }
2702
            } else {
2703
                panic!("Aux tag not found");
2704
            }
2705
        }
2706

2707
        {
2708
            let tag = b"XC";
2709
            if let Ok(Aux::ArrayI16(array)) = test_record.aux(tag) {
2710
                // Retrieve aux array
2711
                let aux_array_content = array.iter().collect::<Vec<_>>();
2712
                assert_eq!(aux_array_content, array_i16);
2713

2714
                // Copy the stored aux array to another record
2715
                {
2716
                    let mut copy_test_record = test_record.clone();
2717

2718
                    // Pushing a field with an existing tag should fail
2719
                    assert!(copy_test_record.push_aux(tag, Aux::I16(3)).is_err());
2720

2721
                    // Remove aux array from target record
2722
                    copy_test_record.remove_aux(tag).unwrap();
2723
                    assert!(copy_test_record.aux(tag).is_err());
2724

2725
                    // Copy array to target record
2726
                    let src_aux = test_record.aux(tag).unwrap();
2727
                    assert!(copy_test_record.push_aux(tag, src_aux).is_ok());
2728
                    if let Ok(Aux::ArrayI16(array)) = copy_test_record.aux(tag) {
2729
                        let aux_array_content_copied = array.iter().collect::<Vec<_>>();
2730
                        assert_eq!(aux_array_content_copied, array_i16);
2731
                    } else {
2732
                        panic!("Aux tag not found");
2733
                    }
2734
                }
2735
            } else {
2736
                panic!("Aux tag not found");
2737
            }
2738
        }
2739

2740
        {
2741
            let tag = b"XD";
2742
            if let Ok(Aux::ArrayU16(array)) = test_record.aux(tag) {
2743
                // Retrieve aux array
2744
                let aux_array_content = array.iter().collect::<Vec<_>>();
2745
                assert_eq!(aux_array_content, array_u16);
2746

2747
                // Copy the stored aux array to another record
2748
                {
2749
                    let mut copy_test_record = test_record.clone();
2750

2751
                    // Pushing a field with an existing tag should fail
2752
                    assert!(copy_test_record.push_aux(tag, Aux::U16(3)).is_err());
2753

2754
                    // Remove aux array from target record
2755
                    copy_test_record.remove_aux(tag).unwrap();
2756
                    assert!(copy_test_record.aux(tag).is_err());
2757

2758
                    // Copy array to target record
2759
                    let src_aux = test_record.aux(tag).unwrap();
2760
                    assert!(copy_test_record.push_aux(tag, src_aux).is_ok());
2761
                    if let Ok(Aux::ArrayU16(array)) = copy_test_record.aux(tag) {
2762
                        let aux_array_content_copied = array.iter().collect::<Vec<_>>();
2763
                        assert_eq!(aux_array_content_copied, array_u16);
2764
                    } else {
2765
                        panic!("Aux tag not found");
2766
                    }
2767
                }
2768
            } else {
2769
                panic!("Aux tag not found");
2770
            }
2771
        }
2772

2773
        {
2774
            let tag = b"XE";
2775
            if let Ok(Aux::ArrayI32(array)) = test_record.aux(tag) {
2776
                // Retrieve aux array
2777
                let aux_array_content = array.iter().collect::<Vec<_>>();
2778
                assert_eq!(aux_array_content, array_i32);
2779

2780
                // Copy the stored aux array to another record
2781
                {
2782
                    let mut copy_test_record = test_record.clone();
2783

2784
                    // Pushing a field with an existing tag should fail
2785
                    assert!(copy_test_record.push_aux(tag, Aux::I32(3)).is_err());
2786

2787
                    // Remove aux array from target record
2788
                    copy_test_record.remove_aux(tag).unwrap();
2789
                    assert!(copy_test_record.aux(tag).is_err());
2790

2791
                    // Copy array to target record
2792
                    let src_aux = test_record.aux(tag).unwrap();
2793
                    assert!(copy_test_record.push_aux(tag, src_aux).is_ok());
2794
                    if let Ok(Aux::ArrayI32(array)) = copy_test_record.aux(tag) {
2795
                        let aux_array_content_copied = array.iter().collect::<Vec<_>>();
2796
                        assert_eq!(aux_array_content_copied, array_i32);
2797
                    } else {
2798
                        panic!("Aux tag not found");
2799
                    }
2800
                }
2801
            } else {
2802
                panic!("Aux tag not found");
2803
            }
2804
        }
2805

2806
        {
2807
            let tag = b"XF";
2808
            if let Ok(Aux::ArrayU32(array)) = test_record.aux(tag) {
2809
                // Retrieve aux array
2810
                let aux_array_content = array.iter().collect::<Vec<_>>();
2811
                assert_eq!(aux_array_content, array_u32);
2812

2813
                // Copy the stored aux array to another record
2814
                {
2815
                    let mut copy_test_record = test_record.clone();
2816

2817
                    // Pushing a field with an existing tag should fail
2818
                    assert!(copy_test_record.push_aux(tag, Aux::U32(3)).is_err());
2819

2820
                    // Remove aux array from target record
2821
                    copy_test_record.remove_aux(tag).unwrap();
2822
                    assert!(copy_test_record.aux(tag).is_err());
2823

2824
                    // Copy array to target record
2825
                    let src_aux = test_record.aux(tag).unwrap();
2826
                    assert!(copy_test_record.push_aux(tag, src_aux).is_ok());
2827
                    if let Ok(Aux::ArrayU32(array)) = copy_test_record.aux(tag) {
2828
                        let aux_array_content_copied = array.iter().collect::<Vec<_>>();
2829
                        assert_eq!(aux_array_content_copied, array_u32);
2830
                    } else {
2831
                        panic!("Aux tag not found");
2832
                    }
2833
                }
2834
            } else {
2835
                panic!("Aux tag not found");
2836
            }
2837
        }
2838

2839
        {
2840
            let tag = b"XG";
2841
            if let Ok(Aux::ArrayFloat(array)) = test_record.aux(tag) {
2842
                // Retrieve aux array
2843
                let aux_array_content = array.iter().collect::<Vec<_>>();
2844
                assert_eq!(aux_array_content, array_f32);
2845

2846
                // Copy the stored aux array to another record
2847
                {
2848
                    let mut copy_test_record = test_record.clone();
2849

2850
                    // Pushing a field with an existing tag should fail
2851
                    assert!(copy_test_record.push_aux(tag, Aux::Float(3.0)).is_err());
2852

2853
                    // Remove aux array from target record
2854
                    copy_test_record.remove_aux(tag).unwrap();
2855
                    assert!(copy_test_record.aux(tag).is_err());
2856

2857
                    // Copy array to target record
2858
                    let src_aux = test_record.aux(tag).unwrap();
2859
                    assert!(copy_test_record.push_aux(tag, src_aux).is_ok());
2860
                    if let Ok(Aux::ArrayFloat(array)) = copy_test_record.aux(tag) {
2861
                        let aux_array_content_copied = array.iter().collect::<Vec<_>>();
2862
                        assert_eq!(aux_array_content_copied, array_f32);
2863
                    } else {
2864
                        panic!("Aux tag not found");
2865
                    }
2866
                }
2867
            } else {
2868
                panic!("Aux tag not found");
2869
            }
2870
        }
2871

2872
        // Test via `Iterator` impl
2873
        for item in test_record.aux_iter() {
2874
            match item.unwrap() {
2875
                (b"XA", Aux::ArrayI8(array)) => {
2876
                    assert_eq!(&array.iter().collect::<Vec<_>>(), &array_i8);
2877
                }
2878
                (b"XB", Aux::ArrayU8(array)) => {
2879
                    assert_eq!(&array.iter().collect::<Vec<_>>(), &array_u8);
2880
                }
2881
                (b"XC", Aux::ArrayI16(array)) => {
2882
                    assert_eq!(&array.iter().collect::<Vec<_>>(), &array_i16);
2883
                }
2884
                (b"XD", Aux::ArrayU16(array)) => {
2885
                    assert_eq!(&array.iter().collect::<Vec<_>>(), &array_u16);
2886
                }
2887
                (b"XE", Aux::ArrayI32(array)) => {
2888
                    assert_eq!(&array.iter().collect::<Vec<_>>(), &array_i32);
2889
                }
2890
                (b"XF", Aux::ArrayU32(array)) => {
2891
                    assert_eq!(&array.iter().collect::<Vec<_>>(), &array_u32);
2892
                }
2893
                (b"XG", Aux::ArrayFloat(array)) => {
2894
                    assert_eq!(&array.iter().collect::<Vec<_>>(), &array_f32);
2895
                }
2896
                _ => {
2897
                    panic!();
2898
                }
2899
            }
2900
        }
2901

2902
        // Test via `PartialEq` impl
2903
        assert_eq!(
2904
            test_record.aux(b"XA").unwrap(),
2905
            Aux::ArrayI8((&array_i8).into())
2906
        );
2907
        assert_eq!(
2908
            test_record.aux(b"XB").unwrap(),
2909
            Aux::ArrayU8((&array_u8).into())
2910
        );
2911
        assert_eq!(
2912
            test_record.aux(b"XC").unwrap(),
2913
            Aux::ArrayI16((&array_i16).into())
2914
        );
2915
        assert_eq!(
2916
            test_record.aux(b"XD").unwrap(),
2917
            Aux::ArrayU16((&array_u16).into())
2918
        );
2919
        assert_eq!(
2920
            test_record.aux(b"XE").unwrap(),
2921
            Aux::ArrayI32((&array_i32).into())
2922
        );
2923
        assert_eq!(
2924
            test_record.aux(b"XF").unwrap(),
2925
            Aux::ArrayU32((&array_u32).into())
2926
        );
2927
        assert_eq!(
2928
            test_record.aux(b"XG").unwrap(),
2929
            Aux::ArrayFloat((&array_f32).into())
2930
        );
2931
    }
2932

2933
    #[test]
2934
    fn test_aux_scalars() {
2935
        let bam_header = Header::new();
2936
        let mut test_record = Record::from_sam(
2937
            &HeaderView::from_header(&bam_header),
2938
            "ali1\t4\t*\t0\t0\t*\t*\t0\t0\tACGT\tFFFF".as_bytes(),
2939
        )
2940
        .unwrap();
2941

2942
        test_record.push_aux(b"XA", Aux::I8(i8::MIN)).unwrap();
2943
        test_record.push_aux(b"XB", Aux::I8(i8::MAX)).unwrap();
2944
        test_record.push_aux(b"XC", Aux::U8(u8::MIN)).unwrap();
2945
        test_record.push_aux(b"XD", Aux::U8(u8::MAX)).unwrap();
2946
        test_record.push_aux(b"XE", Aux::I16(i16::MIN)).unwrap();
2947
        test_record.push_aux(b"XF", Aux::I16(i16::MAX)).unwrap();
2948
        test_record.push_aux(b"XG", Aux::U16(u16::MIN)).unwrap();
2949
        test_record.push_aux(b"XH", Aux::U16(u16::MAX)).unwrap();
2950
        test_record.push_aux(b"XI", Aux::I32(i32::MIN)).unwrap();
2951
        test_record.push_aux(b"XJ", Aux::I32(i32::MAX)).unwrap();
2952
        test_record.push_aux(b"XK", Aux::U32(u32::MIN)).unwrap();
2953
        test_record.push_aux(b"XL", Aux::U32(u32::MAX)).unwrap();
2954
        test_record
2955
            .push_aux(b"XM", Aux::Float(std::f32::consts::PI))
2956
            .unwrap();
2957
        test_record
2958
            .push_aux(b"XN", Aux::Double(std::f64::consts::PI))
2959
            .unwrap();
2960
        test_record
2961
            .push_aux(b"XO", Aux::String("Test str"))
2962
            .unwrap();
2963
        test_record.push_aux(b"XP", Aux::I8(0)).unwrap();
2964

2965
        let collected_aux_fields = test_record.aux_iter().collect::<Result<Vec<_>>>().unwrap();
2966
        assert_eq!(
2967
            collected_aux_fields,
2968
            vec![
2969
                (&b"XA"[..], Aux::I8(i8::MIN)),
2970
                (&b"XB"[..], Aux::I8(i8::MAX)),
2971
                (&b"XC"[..], Aux::U8(u8::MIN)),
2972
                (&b"XD"[..], Aux::U8(u8::MAX)),
2973
                (&b"XE"[..], Aux::I16(i16::MIN)),
2974
                (&b"XF"[..], Aux::I16(i16::MAX)),
2975
                (&b"XG"[..], Aux::U16(u16::MIN)),
2976
                (&b"XH"[..], Aux::U16(u16::MAX)),
2977
                (&b"XI"[..], Aux::I32(i32::MIN)),
2978
                (&b"XJ"[..], Aux::I32(i32::MAX)),
2979
                (&b"XK"[..], Aux::U32(u32::MIN)),
2980
                (&b"XL"[..], Aux::U32(u32::MAX)),
2981
                (&b"XM"[..], Aux::Float(std::f32::consts::PI)),
2982
                (&b"XN"[..], Aux::Double(std::f64::consts::PI)),
2983
                (&b"XO"[..], Aux::String("Test str")),
2984
                (&b"XP"[..], Aux::I8(0)),
2985
            ]
2986
        );
2987
    }
2988

2989
    #[test]
2990
    fn test_aux_array_partial_eq() {
2991
        use record::AuxArray;
2992

2993
        // Target types
2994
        let one_data: Vec<i8> = vec![0, 1, 2, 3, 4, 5, 6];
2995
        let one_aux_array = AuxArray::from(&one_data);
2996

2997
        let two_data: Vec<i8> = vec![0, 1, 2, 3, 4, 5];
2998
        let two_aux_array = AuxArray::from(&two_data);
2999

3000
        assert_ne!(&one_data, &two_data);
3001
        assert_ne!(&one_aux_array, &two_aux_array);
3002

3003
        let one_aux = Aux::ArrayI8(one_aux_array);
3004
        let two_aux = Aux::ArrayI8(two_aux_array);
3005
        assert_ne!(&one_aux, &two_aux);
3006

3007
        // Raw bytes
3008
        let bam_header = Header::new();
3009
        let mut test_record = Record::from_sam(
3010
            &HeaderView::from_header(&bam_header),
3011
            "ali1\t4\t*\t0\t0\t*\t*\t0\t0\tACGT\tFFFF".as_bytes(),
3012
        )
3013
        .unwrap();
3014

3015
        test_record.push_aux(b"XA", one_aux).unwrap();
3016
        test_record.push_aux(b"XB", two_aux).unwrap();
3017

3018
        // RawLeBytes == RawLeBytes
3019
        assert_eq!(
3020
            test_record.aux(b"XA").unwrap(),
3021
            test_record.aux(b"XA").unwrap()
3022
        );
3023
        // RawLeBytes != RawLeBytes
3024
        assert_ne!(
3025
            test_record.aux(b"XA").unwrap(),
3026
            test_record.aux(b"XB").unwrap()
3027
        );
3028

3029
        // RawLeBytes == TargetType
3030
        assert_eq!(
3031
            test_record.aux(b"XA").unwrap(),
3032
            Aux::ArrayI8((&one_data).into())
3033
        );
3034
        assert_eq!(
3035
            test_record.aux(b"XB").unwrap(),
3036
            Aux::ArrayI8((&two_data).into())
3037
        );
3038
        // RawLeBytes != TargetType
3039
        assert_ne!(
3040
            test_record.aux(b"XA").unwrap(),
3041
            Aux::ArrayI8((&two_data).into())
3042
        );
3043
        assert_ne!(
3044
            test_record.aux(b"XB").unwrap(),
3045
            Aux::ArrayI8((&one_data).into())
3046
        );
3047
    }
3048

3049
    /// Test if both text and binary representations of a BAM header are in sync (#156)
3050
    #[test]
3051
    fn test_bam_header_sync() {
3052
        let reader = Reader::from_path("test/test_issue_156_no_text.bam").unwrap();
3053
        let header_hashmap = Header::from_template(reader.header()).to_hashmap();
3054
        let header_refseqs = header_hashmap.get("SQ").unwrap();
3055
        assert_eq!(header_refseqs[0].get("SN").unwrap(), "ref_1",);
3056
        assert_eq!(header_refseqs[0].get("LN").unwrap(), "10000000",);
3057
    }
3058

3059
    #[test]
3060
    fn test_bam_new() {
3061
        // Create the path to write the tmp test BAM
3062
        let tmp = tempfile::Builder::new()
3063
            .prefix("rust-htslib")
3064
            .tempdir()
3065
            .expect("Cannot create temp dir");
3066
        let bampath = tmp.path().join("test.bam");
3067

3068
        // write an unmapped BAM record (uBAM)
3069
        {
3070
            // Build the header
3071
            let mut header = Header::new();
3072

3073
            // Add the version
3074
            header.push_record(
3075
                HeaderRecord::new(b"HD")
3076
                    .push_tag(b"VN", "1.6")
3077
                    .push_tag(b"SO", "unsorted"),
3078
            );
3079

3080
            // Build the writer
3081
            let mut writer = Writer::from_path(&bampath, &header, Format::Bam).unwrap();
3082

3083
            // Build an empty record
3084
            let record = Record::new();
3085

3086
            // Write the record (this previously seg-faulted)
3087
            assert!(writer.write(&record).is_ok());
3088
        }
3089

3090
        // Read the record
3091
        {
3092
            // Build th reader
3093
            let mut reader = Reader::from_path(bampath).expect("Error opening file.");
3094

3095
            // Read the record
3096
            let mut rec = Record::new();
3097
            match reader.read(&mut rec) {
3098
                Some(r) => r.expect("Failed to read record."),
3099
                None => panic!("No record read."),
3100
            };
3101

3102
            // Check a few things
3103
            assert!(rec.is_unmapped());
3104
            assert_eq!(rec.tid(), -1);
3105
            assert_eq!(rec.pos(), -1);
3106
            assert_eq!(rec.mtid(), -1);
3107
            assert_eq!(rec.mpos(), -1);
3108
        }
3109
    }
3110

3111
    #[test]
3112
    fn test_idxstats_bam() {
3113
        let mut reader = IndexedReader::from_path("test/test.bam").unwrap();
3114
        let expected = vec![
3115
            (0, 15072423, 6, 0),
3116
            (1, 15279345, 0, 0),
3117
            (2, 13783700, 0, 0),
3118
            (3, 17493793, 0, 0),
3119
            (4, 20924149, 0, 0),
3120
            (-1, 0, 0, 0),
3121
        ];
3122
        let actual = reader.index_stats().unwrap();
3123
        assert_eq!(expected, actual);
3124
    }
3125

3126
    #[test]
3127
    fn test_number_mapped_and_unmapped_bam() {
3128
        let reader = IndexedReader::from_path("test/test.bam").unwrap();
3129
        let expected = (6, 0);
3130
        let actual = reader.index().number_mapped_unmapped(0);
3131
        assert_eq!(expected, actual);
3132
    }
3133

3134
    #[test]
3135
    fn test_number_unmapped_global_bam() {
3136
        let reader = IndexedReader::from_path("test/test_unmapped.bam").unwrap();
3137
        let expected = 8;
3138
        let actual = reader.index().number_unmapped();
3139
        assert_eq!(expected, actual);
3140
    }
3141

3142
    #[test]
3143
    fn test_idxstats_cram() {
3144
        let mut reader = IndexedReader::from_path("test/test_cram.cram").unwrap();
3145
        reader.set_reference("test/test_cram.fa").unwrap();
3146
        let expected = vec![
3147
            (0, 120, 2, 0),
3148
            (1, 120, 2, 0),
3149
            (2, 120, 2, 0),
3150
            (-1, 0, 0, 0),
3151
        ];
3152
        let actual = reader.index_stats().unwrap();
3153
        assert_eq!(expected, actual);
3154
    }
3155

3156
    #[test]
3157
    fn test_slow_idxstats_cram() {
3158
        let mut reader = IndexedReader::from_path("test/test_cram.cram").unwrap();
3159
        reader.set_reference("test/test_cram.fa").unwrap();
3160
        let expected = vec![
3161
            (0, 120, 2, 0),
3162
            (1, 120, 2, 0),
3163
            (2, 120, 2, 0),
3164
            (-1, 0, 0, 0),
3165
        ];
3166
        let actual = reader.index_stats().unwrap();
3167
        assert_eq!(expected, actual);
3168
    }
3169

3170
    // #[test]
3171
    // fn test_number_mapped_and_unmapped_cram() {
3172
    //     let mut reader = IndexedReader::from_path("test/test_cram.cram").unwrap();
3173
    //     reader.set_reference("test/test_cram.fa").unwrap();
3174
    //     let expected = (2, 0);
3175
    //     let actual = reader.index().number_mapped_unmapped(0);
3176
    //     assert_eq!(expected, actual);
3177
    // }
3178
    //
3179
    // #[test]
3180
    // fn test_number_unmapped_global_cram() {
3181
    //     let mut reader = IndexedReader::from_path("test/test_unmapped.cram").unwrap();
3182
    //     let expected = 8;
3183
    //     let actual = reader.index().number_unmapped();
3184
    //     assert_eq!(expected, actual);
3185
    // }
3186
}
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