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

rust-bio / rust-htslib / 19161687300

07 Nov 2025 07:43AM UTC coverage: 81.912% (-0.02%) from 81.935%
19161687300

Pull #488

github

web-flow
Merge 907dabbe9 into 8f1cdd75c
Pull Request #488: fix: Reason about Send/Sync-ness of types and change Rcs to Arcs

34 of 36 new or added lines in 7 files covered. (94.44%)

20 existing lines in 6 files now uncovered.

2785 of 3400 relevant lines covered (81.91%)

27209.26 hits per line

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

83.98
/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
use std::sync::Arc;
25

26
use url::Url;
27

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

258
unsafe impl Send for Reader {}
259

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

601
unsafe impl Send for IndexedReader {}
602

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

914
unsafe impl Send for IndexView {}
915
unsafe impl Sync for IndexView {}
916

917
impl IndexView {
918
    fn new(hts_idx: *mut hts_sys::hts_idx_t) -> Self {
16✔
919
        Self { inner: hts_idx }
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
        unsafe {
964
            htslib::hts_idx_destroy(self.inner);
16✔
965
        }
966
    }
967
}
968

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

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

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

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

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

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

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

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

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

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

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

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

1065
unsafe impl Send for Writer {}
1066

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1362
#[derive(Debug)]
1363
pub struct HeaderView {
1364
    inner: *mut htslib::bam_hdr_t,
1365
}
1366

1367
unsafe impl Send for HeaderView {}
1368
unsafe impl Sync for HeaderView {}
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
    fn new(inner: *mut htslib::bam_hdr_t) -> Self {
98✔
1403
        HeaderView { inner }
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
        }
1484
    }
1485
}
1486

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

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

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

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

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

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

1595
            let cigar = rec.cigar();
1596
            assert_eq!(*cigar, cigars[i]);
1597

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

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

1625
        let mut names_by_voffset = HashMap::new();
1626

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

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

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

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

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

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

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

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

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

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

1703
        // repeat with byte-string based fetch
1704

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

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

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

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

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

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

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

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

1773
        //all reads
1774
        bam.fetch(".").expect("Expected successful fetch.");
1775
        assert!(bam.records().count() == 6);
1776

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

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

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

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

1802
    #[test]
1803
    fn test_set_record() {
1804
        let (names, _, seqs, quals, cigars) = gold();
1805

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

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

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

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

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

1859
    #[test]
1860
    fn test_set_qname() {
1861
        let (names, _, seqs, quals, cigars) = gold();
1862

1863
        assert!(names[0] != names[1]);
1864

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

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

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

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

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

1893
            // Shorter qname
1894
            let shorter_name = b"42";
1895
            rec.set_qname(shorter_name);
1896

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

1903
            // Zero-length qname
1904
            rec.set_qname(b"");
1905

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

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

1924
        let line =
1925
            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";
1926

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

1933
    #[test]
1934
    fn test_set_cigar() {
1935
        let (names, _, seqs, quals, cigars) = gold();
1936

1937
        assert!(names[0] != names[1]);
1938

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

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

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

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

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

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

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

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

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

1998
        for record in bam.records() {
1999
            let mut rec = record.expect("Expected valid record");
2000

2001
            if rec.aux(b"XS").is_ok() {
2002
                rec.remove_aux(b"XS").unwrap();
2003
            }
2004

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

2009
            assert!(rec.remove_aux(b"ab").is_err());
2010

2011
            assert_eq!(rec.aux(b"XS"), Err(Error::BamAuxTagNotFound));
2012
            assert_eq!(rec.aux(b"YT"), Err(Error::BamAuxTagNotFound));
2013
        }
2014
    }
2015

2016
    #[test]
2017
    fn test_write() {
2018
        let (names, _, seqs, quals, cigars) = gold();
2019

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

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

2043
                bam.write(&rec).expect("Failed to write record.");
2044
            }
2045
        }
2046

2047
        {
2048
            let mut bam = Reader::from_path(bampath).expect("Error opening file.");
2049

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

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

2064
        tmp.close().expect("Failed to delete temp dir");
2065
    }
2066

2067
    #[test]
2068
    fn test_write_threaded() {
2069
        let (names, _, seqs, quals, cigars) = gold();
2070

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

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

2097
                bam.write(&rec).expect("Failed to write record.");
2098
            }
2099
        }
2100

2101
        {
2102
            let mut bam = Reader::from_path(bampath).expect("Error opening file.");
2103

2104
            for (i, _rec) in bam.records().enumerate() {
2105
                let idx = i % names.len();
2106

2107
                let rec = _rec.expect("Failed to read record.");
2108

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

2118
        tmp.close().expect("Failed to delete temp dir");
2119
    }
2120

2121
    #[test]
2122
    fn test_write_shared_tpool() {
2123
        let (names, _, seqs, quals, cigars) = gold();
2124

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

2132
        {
2133
            let (mut bam1, mut bam2) = {
2134
                let pool = crate::tpool::ThreadPool::new(4).unwrap();
2135

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

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

2158
                bam1.set_thread_pool(&pool).unwrap();
2159
                bam2.set_thread_pool(&pool).unwrap();
2160
                (bam1, bam2)
2161
            };
2162

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

2170
                bam1.write(&rec).expect("Failed to write record.");
2171
                bam2.write(&rec).expect("Failed to write record.");
2172
            }
2173
        }
2174

2175
        {
2176
            let pool = crate::tpool::ThreadPool::new(2).unwrap();
2177

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

2182
                for (i, _rec) in bam.iter_chunk(None, None).enumerate() {
2183
                    let idx = i % names.len();
2184

2185
                    let rec = _rec.expect("Failed to read record.");
2186

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

2197
        tmp.close().expect("Failed to delete temp dir");
2198
    }
2199

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

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

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

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

2222
            for rec in input_bam.records() {
2223
                bam.write(&rec.unwrap()).expect("Failed to write record.");
2224
            }
2225
        }
2226

2227
        {
2228
            let copy_bam = Reader::from_path(bampath).expect("Error opening file.");
2229

2230
            // Verify that the header came across correctly
2231
            assert_eq!(input_bam.header().as_bytes(), copy_bam.header().as_bytes());
2232
        }
2233

2234
        tmp.close().expect("Failed to delete temp dir");
2235
    }
2236

2237
    #[test]
2238
    fn test_pileup() {
2239
        let (_, _, seqs, quals, _) = gold();
2240

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

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

2273
    #[test]
2274
    fn parse_from_sam() {
2275
        use std::fs::File;
2276
        use std::io::Read;
2277

2278
        let bamfile = "./test/bam2sam_test.bam";
2279
        let samfile = "./test/bam2sam_expected.sam";
2280

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

2285
        let mut sam = Vec::new();
2286
        assert!(File::open(samfile).unwrap().read_to_end(&mut sam).is_ok());
2287

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

2294
        for (b1, s1) in bam_recs.iter().zip(sam_recs.iter()) {
2295
            assert!(b1 == s1);
2296
        }
2297
    }
2298

2299
    #[test]
2300
    fn test_cigar_modes() {
2301
        // test the cached and uncached ways of getting the cigar string.
2302

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

2306
        for (i, record) in bam.records().enumerate() {
2307
            let rec = record.expect("Expected valid record");
2308

2309
            let cigar = rec.cigar();
2310
            assert_eq!(*cigar, cigars[i]);
2311
        }
2312

2313
        for (i, record) in bam.records().enumerate() {
2314
            let mut rec = record.expect("Expected valid record");
2315
            rec.cache_cigar();
2316

2317
            let cigar = rec.cigar_cached().unwrap();
2318
            assert_eq!(**cigar, cigars[i]);
2319

2320
            let cigar = rec.cigar();
2321
            assert_eq!(*cigar, cigars[i]);
2322
        }
2323
    }
2324

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

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

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

2340
        compare_inner_bam_cram_records(&cram_records, &bam_records);
2341
    }
2342

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

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

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

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

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

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

2407
            // Compare CRAM records to BAM records
2408
            compare_inner_bam_cram_records(&cram_records, &bam_records);
2409
        }
2410

2411
        tmp.close().expect("Failed to delete temp dir");
2412
    }
2413

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

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

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

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

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

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

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

2475
        tmp.close().expect("Failed to delete temp dir");
2476
    }
2477

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

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

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

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

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

2556
    //     assert_eq!(r.header().target_names()[0], b"chr1");
2557
    // }
2558

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

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

2573
            let cigar = rec.cigar();
2574
            assert_eq!(*cigar, cigars[i]);
2575

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

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

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

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

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

2645
                // Copy the stored aux array to another record
2646
                {
2647
                    let mut copy_test_record = test_record.clone();
2648

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

2652
                    // Remove aux array from target record
2653
                    copy_test_record.remove_aux(tag).unwrap();
2654
                    assert!(copy_test_record.aux(tag).is_err());
2655

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

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

2678
                // Copy the stored aux array to another record
2679
                {
2680
                    let mut copy_test_record = test_record.clone();
2681

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

2685
                    // Remove aux array from target record
2686
                    copy_test_record.remove_aux(tag).unwrap();
2687
                    assert!(copy_test_record.aux(tag).is_err());
2688

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

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

2711
                // Copy the stored aux array to another record
2712
                {
2713
                    let mut copy_test_record = test_record.clone();
2714

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

2718
                    // Remove aux array from target record
2719
                    copy_test_record.remove_aux(tag).unwrap();
2720
                    assert!(copy_test_record.aux(tag).is_err());
2721

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

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

2744
                // Copy the stored aux array to another record
2745
                {
2746
                    let mut copy_test_record = test_record.clone();
2747

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

2751
                    // Remove aux array from target record
2752
                    copy_test_record.remove_aux(tag).unwrap();
2753
                    assert!(copy_test_record.aux(tag).is_err());
2754

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

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

2777
                // Copy the stored aux array to another record
2778
                {
2779
                    let mut copy_test_record = test_record.clone();
2780

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

2784
                    // Remove aux array from target record
2785
                    copy_test_record.remove_aux(tag).unwrap();
2786
                    assert!(copy_test_record.aux(tag).is_err());
2787

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

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

2810
                // Copy the stored aux array to another record
2811
                {
2812
                    let mut copy_test_record = test_record.clone();
2813

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

2817
                    // Remove aux array from target record
2818
                    copy_test_record.remove_aux(tag).unwrap();
2819
                    assert!(copy_test_record.aux(tag).is_err());
2820

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

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

2843
                // Copy the stored aux array to another record
2844
                {
2845
                    let mut copy_test_record = test_record.clone();
2846

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

2850
                    // Remove aux array from target record
2851
                    copy_test_record.remove_aux(tag).unwrap();
2852
                    assert!(copy_test_record.aux(tag).is_err());
2853

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

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

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

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

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

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

2986
    #[test]
2987
    fn test_aux_array_partial_eq() {
2988
        use record::AuxArray;
2989

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

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

2997
        assert_ne!(&one_data, &two_data);
2998
        assert_ne!(&one_aux_array, &two_aux_array);
2999

3000
        let one_aux = Aux::ArrayI8(one_aux_array);
3001
        let two_aux = Aux::ArrayI8(two_aux_array);
3002
        assert_ne!(&one_aux, &two_aux);
3003

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

3012
        test_record.push_aux(b"XA", one_aux).unwrap();
3013
        test_record.push_aux(b"XB", two_aux).unwrap();
3014

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

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

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

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

3065
        // write an unmapped BAM record (uBAM)
3066
        {
3067
            // Build the header
3068
            let mut header = Header::new();
3069

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

3077
            // Build the writer
3078
            let mut writer = Writer::from_path(&bampath, &header, Format::Bam).unwrap();
3079

3080
            // Build an empty record
3081
            let record = Record::new();
3082

3083
            // Write the record (this previously seg-faulted)
3084
            assert!(writer.write(&record).is_ok());
3085
        }
3086

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

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

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

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

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

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

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

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

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