From f59e9ee22dddf2f63e7d49585acb842c455c9c1d Mon Sep 17 00:00:00 2001 From: Till Hartmann Date: Tue, 2 May 2023 12:59:38 +0200 Subject: [PATCH 01/20] Add an IndexView for accessing a sam/bam/cram's index --- src/bam/mod.rs | 91 ++++++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 85 insertions(+), 6 deletions(-) diff --git a/src/bam/mod.rs b/src/bam/mod.rs index 67b60431c..09e3d87b2 100644 --- a/src/bam/mod.rs +++ b/src/bam/mod.rs @@ -563,7 +563,7 @@ impl<'a, T: AsRef<[u8]>, X: Into, Y: Into> Fro pub struct IndexedReader { htsfile: *mut htslib::htsFile, header: Rc, - idx: *mut htslib::hts_idx_t, + idx: Rc, itr: Option<*mut htslib::hts_itr_t>, tpool: Option, } @@ -609,7 +609,7 @@ impl IndexedReader { Ok(IndexedReader { htsfile, header: Rc::new(HeaderView::new(header)), - idx, + idx: Rc::new(IndexView::new(idx)), itr: None, tpool: None, }) @@ -637,7 +637,7 @@ impl IndexedReader { Ok(IndexedReader { htsfile, header: Rc::new(HeaderView::new(header)), - idx, + idx: Rc::new(IndexView::new(idx)), itr: None, tpool: None, }) @@ -726,7 +726,9 @@ impl IndexedReader { if let Some(itr) = self.itr { unsafe { htslib::hts_itr_destroy(itr) } } - let itr = unsafe { htslib::sam_itr_queryi(self.idx, tid as i32, beg as i64, end as i64) }; + let itr = unsafe { + htslib::sam_itr_queryi(self.index().inner_ptr(), tid as i32, beg as i64, end as i64) + }; if itr.is_null() { self.itr = None; Err(Error::Fetch) @@ -744,7 +746,7 @@ impl IndexedReader { let rptr = rstr.as_ptr(); let itr = unsafe { htslib::sam_itr_querys( - self.idx, + self.index().inner_ptr(), self.header().inner_ptr() as *mut hts_sys::sam_hdr_t, rptr, ) @@ -783,6 +785,84 @@ impl IndexedReader { pub fn set_reference>(&mut self, path: P) -> Result<()> { unsafe { set_fai_filename(self.htsfile, path) } } + + pub fn index(&mut self) -> &IndexView { + &self.idx + } +} + +#[derive(Debug)] +pub struct IndexView { + inner: *mut hts_sys::hts_idx_t, + owned: bool, +} + +impl IndexView { + fn new(hts_idx: *mut hts_sys::hts_idx_t) -> Self { + Self { + inner: hts_idx, + owned: true, + } + } + + #[inline] + pub fn inner(&self) -> &hts_sys::hts_idx_t { + unsafe { self.inner.as_ref().unwrap() } + } + + #[inline] + // Pointer to inner bam_hdr_t struct + pub fn inner_ptr(&self) -> *const hts_sys::hts_idx_t { + self.inner + } + + #[inline] + pub fn inner_mut(&mut self) -> &mut hts_sys::hts_idx_t { + unsafe { self.inner.as_mut().unwrap() } + } + + #[inline] + // Mutable pointer to bam_hdr_t struct + pub fn inner_ptr_mut(&mut self) -> *mut hts_sys::hts_idx_t { + self.inner + } + + /// Get the number of mapped and unmapped reads for a given target id + pub fn number_mapped_unmapped(&self, tid: u32) -> (u64, u64) { + let (mut mapped, mut unmapped) = (0, 0); + unsafe { + hts_sys::hts_idx_get_stat(self.inner, tid as i32, &mut mapped, &mut unmapped); + } + (mapped, unmapped) + } + + /// Similar to samtools idxstats, this returns a vector of tuples + /// containing the target name, length, number of mapped reads, and number of unmapped reads. + pub fn stats(&self, header: HeaderView) -> Vec<(String, (u64, u64, u64))> { + header + .target_names() + .iter() + .map(|tname| { + let tid = header.tid(tname).unwrap(); + let (mapped, unmapped) = self.number_mapped_unmapped(tid); + let tname = str::from_utf8(tname).unwrap(); + let tlen = header.target_len(tid as u32).unwrap(); + (tname.to_string(), (tlen, mapped, unmapped)) + }) + .collect() + } +} + +impl Drop for IndexView { + fn drop(&mut self) { + unsafe { + if self.owned { + unsafe { + htslib::hts_idx_destroy(self.inner); + } + } + } + } } impl Read for IndexedReader { @@ -851,7 +931,6 @@ impl Drop for IndexedReader { if self.itr.is_some() { htslib::hts_itr_destroy(self.itr.unwrap()); } - htslib::hts_idx_destroy(self.idx); htslib::hts_close(self.htsfile); } } From b8727dd8eff2ae2dda8f7f33f6b9e27c94cf1463 Mon Sep 17 00:00:00 2001 From: Till Hartmann Date: Tue, 2 May 2023 13:00:21 +0200 Subject: [PATCH 02/20] fix copy and paste error --- src/bam/mod.rs | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/bam/mod.rs b/src/bam/mod.rs index 09e3d87b2..6b1c6a371 100644 --- a/src/bam/mod.rs +++ b/src/bam/mod.rs @@ -811,7 +811,7 @@ impl IndexView { } #[inline] - // Pointer to inner bam_hdr_t struct + // Pointer to inner hts_idx_t struct pub fn inner_ptr(&self) -> *const hts_sys::hts_idx_t { self.inner } @@ -822,7 +822,7 @@ impl IndexView { } #[inline] - // Mutable pointer to bam_hdr_t struct + // Mutable pointer to hts_idx_t struct pub fn inner_ptr_mut(&mut self) -> *mut hts_sys::hts_idx_t { self.inner } From 02e9d561b8ed6b3651da0dcc2c9d11b0065a3794 Mon Sep 17 00:00:00 2001 From: Till Hartmann Date: Tue, 2 May 2023 13:08:55 +0200 Subject: [PATCH 03/20] remove unnecessary cast to u32 --- src/bam/mod.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/bam/mod.rs b/src/bam/mod.rs index 6b1c6a371..c04f5a165 100644 --- a/src/bam/mod.rs +++ b/src/bam/mod.rs @@ -846,7 +846,7 @@ impl IndexView { let tid = header.tid(tname).unwrap(); let (mapped, unmapped) = self.number_mapped_unmapped(tid); let tname = str::from_utf8(tname).unwrap(); - let tlen = header.target_len(tid as u32).unwrap(); + let tlen = header.target_len(tid).unwrap(); (tname.to_string(), (tlen, mapped, unmapped)) }) .collect() From a924119cfcdcd8a5ed4093c7a5f6f5d9b7379b94 Mon Sep 17 00:00:00 2001 From: Till Hartmann Date: Tue, 2 May 2023 13:21:51 +0200 Subject: [PATCH 04/20] index.stats: return Option<_> --- src/bam/mod.rs | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/bam/mod.rs b/src/bam/mod.rs index c04f5a165..857812624 100644 --- a/src/bam/mod.rs +++ b/src/bam/mod.rs @@ -837,19 +837,19 @@ impl IndexView { } /// Similar to samtools idxstats, this returns a vector of tuples - /// containing the target name, length, number of mapped reads, and number of unmapped reads. - pub fn stats(&self, header: HeaderView) -> Vec<(String, (u64, u64, u64))> { + /// containing the target id, length, number of mapped reads, and number of unmapped reads. + pub fn stats(&self, header: HeaderView) -> Option> { header .target_names() .iter() .map(|tname| { - let tid = header.tid(tname).unwrap(); - let (mapped, unmapped) = self.number_mapped_unmapped(tid); - let tname = str::from_utf8(tname).unwrap(); - let tlen = header.target_len(tid).unwrap(); - (tname.to_string(), (tlen, mapped, unmapped)) + header.tid(tname).and_then(|tid| { + let (mapped, unmapped) = self.number_mapped_unmapped(tid); + let tlen = header.target_len(tid); + tlen.map(|tlen| (tid, (tlen, mapped, unmapped))) + }) }) - .collect() + .collect::>() } } From 5d16863daeaaa3d025dc3937472c9689f40e9327 Mon Sep 17 00:00:00 2001 From: Till Hartmann Date: Tue, 2 May 2023 15:11:02 +0200 Subject: [PATCH 05/20] clippy lints --- src/bam/mod.rs | 19 +++++-------------- 1 file changed, 5 insertions(+), 14 deletions(-) diff --git a/src/bam/mod.rs b/src/bam/mod.rs index 857812624..6474b5589 100644 --- a/src/bam/mod.rs +++ b/src/bam/mod.rs @@ -706,12 +706,7 @@ impl IndexedReader { match tid { Some(tid) => { //'large position' spec says target len must will fit into an i64. - let len: i64 = self - .header - .target_len(tid as u32) - .unwrap() - .try_into() - .unwrap(); + let len: i64 = self.header.target_len(tid).unwrap().try_into().unwrap(); self._fetch_by_coord_tuple(tid as i32, 0, len) } None => self._fetch_by_str(s), @@ -726,9 +721,7 @@ impl IndexedReader { if let Some(itr) = self.itr { unsafe { htslib::hts_itr_destroy(itr) } } - let itr = unsafe { - htslib::sam_itr_queryi(self.index().inner_ptr(), tid as i32, beg as i64, end as i64) - }; + let itr = unsafe { htslib::sam_itr_queryi(self.index().inner_ptr(), tid, beg, end) }; if itr.is_null() { self.itr = None; Err(Error::Fetch) @@ -855,11 +848,9 @@ impl IndexView { impl Drop for IndexView { fn drop(&mut self) { - unsafe { - if self.owned { - unsafe { - htslib::hts_idx_destroy(self.inner); - } + if self.owned { + unsafe { + htslib::hts_idx_destroy(self.inner); } } } From 74e6d1a81e3989b677c5ab8caf6849bedcbfde35 Mon Sep 17 00:00:00 2001 From: Till Hartmann Date: Tue, 2 May 2023 15:11:25 +0200 Subject: [PATCH 06/20] fmt --- src/bam/mod.rs | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/bam/mod.rs b/src/bam/mod.rs index 6474b5589..32a3babcb 100644 --- a/src/bam/mod.rs +++ b/src/bam/mod.rs @@ -106,7 +106,7 @@ pub trait Read: Sized { /// match result { /// Ok(_) => { /// println!("Read sequence: {:?}", record.seq().as_bytes()); - /// }, + /// } /// Err(_) => panic!("BAM parsing failed...") /// } /// } @@ -471,6 +471,7 @@ pub enum FetchDefinition<'a> { /// Only reads with the BAM flag BAM_FUNMAP (which might not be all reads with reference = -1) Unmapped, } + impl<'a, X: Into, Y: Into> From<(i32, X, Y)> for FetchDefinition<'a> { @@ -480,6 +481,7 @@ impl<'a, X: Into, Y: Into> From<(i32, X, Y)> FetchDefinition::Region(tup.0, start.0, stop.0) } } + impl<'a, X: Into, Y: Into> From<(u32, X, Y)> for FetchDefinition<'a> { @@ -1685,6 +1687,7 @@ CCCCCCCCCCCCCCCCCCC"[..], let bam = IndexedReader::from_path(&"test/test.bam").expect("Expected valid index."); _test_read_indexed_common(bam); } + #[test] fn test_read_indexed_different_index_name() { let bam = IndexedReader::from_path_and_index( From d96be2873223044abc8e66995f2484a0242108b3 Mon Sep 17 00:00:00 2001 From: Till Hartmann Date: Tue, 2 May 2023 15:22:51 +0200 Subject: [PATCH 07/20] add test for idxstats --- src/bam/mod.rs | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/src/bam/mod.rs b/src/bam/mod.rs index 32a3babcb..0eeccf651 100644 --- a/src/bam/mod.rs +++ b/src/bam/mod.rs @@ -2891,4 +2891,19 @@ CCCCCCCCCCCCCCCCCCC"[..], assert_eq!(header_refseqs[0].get("SN").unwrap(), "ref_1",); assert_eq!(header_refseqs[0].get("LN").unwrap(), "10000000",); } + + #[test] + fn test_idxstats() { + let mut reader = IndexedReader::from_path("test/test.bam").unwrap(); + let header = reader.header().clone(); + let expected = vec![ + (0, (15072423, 6, 0)), + (1, (15279345, 0, 0)), + (2, (13783700, 0, 0)), + (3, (17493793, 0, 0)), + (4, (20924149, 0, 0)), + ]; + let actual = reader.index().stats(header).unwrap(); + assert_eq!(expected, actual); + } } From 7ceb1ea1e4dc115e77d125e35cbdaf837b5cd5d0 Mon Sep 17 00:00:00 2001 From: Till Hartmann Date: Tue, 2 May 2023 15:24:07 +0200 Subject: [PATCH 08/20] add test for number_mapped_unmapped --- src/bam/mod.rs | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/bam/mod.rs b/src/bam/mod.rs index 0eeccf651..c679321ad 100644 --- a/src/bam/mod.rs +++ b/src/bam/mod.rs @@ -2906,4 +2906,12 @@ CCCCCCCCCCCCCCCCCCC"[..], let actual = reader.index().stats(header).unwrap(); assert_eq!(expected, actual); } + + #[test] + fn test_number_mapped_and_unmapped() { + let mut reader = IndexedReader::from_path("test/test.bam").unwrap(); + let expected = (6, 0); + let actual = reader.index().number_mapped_unmapped(0); + assert_eq!(expected, actual); + } } From a24b78329f9e5fbd307fc53274ed4426bf6ced6c Mon Sep 17 00:00:00 2001 From: Till Hartmann Date: Tue, 2 May 2023 15:50:21 +0200 Subject: [PATCH 09/20] add tests for cram index as well --- src/bam/mod.rs | 21 +++++++++++++++++++-- test/test_cram.cram.crai | Bin 0 -> 66 bytes 2 files changed, 19 insertions(+), 2 deletions(-) create mode 100644 test/test_cram.cram.crai diff --git a/src/bam/mod.rs b/src/bam/mod.rs index c679321ad..dc9d00ae4 100644 --- a/src/bam/mod.rs +++ b/src/bam/mod.rs @@ -2893,7 +2893,7 @@ CCCCCCCCCCCCCCCCCCC"[..], } #[test] - fn test_idxstats() { + fn test_idxstats_bam() { let mut reader = IndexedReader::from_path("test/test.bam").unwrap(); let header = reader.header().clone(); let expected = vec![ @@ -2908,10 +2908,27 @@ CCCCCCCCCCCCCCCCCCC"[..], } #[test] - fn test_number_mapped_and_unmapped() { + fn test_number_mapped_and_unmapped_bam() { let mut reader = IndexedReader::from_path("test/test.bam").unwrap(); let expected = (6, 0); let actual = reader.index().number_mapped_unmapped(0); assert_eq!(expected, actual); } + + #[test] + fn test_idxstats_cram() { + let mut reader = IndexedReader::from_path("test/test_cram.cram").unwrap(); + let header = reader.header().clone(); + let expected = vec![(0, (120, 2, 0)), (1, (120, 2, 0)), (2, (120, 2, 0))]; + let actual = reader.index().stats(header).unwrap(); + assert_eq!(expected, actual); + } + + #[test] + fn test_number_mapped_and_unmapped_cram() { + let mut reader = IndexedReader::from_path("test/test_cram.cram").unwrap(); + let expected = (2, 0); + let actual = reader.index().number_mapped_unmapped(0); + assert_eq!(expected, actual); + } } diff --git a/test/test_cram.cram.crai b/test/test_cram.cram.crai new file mode 100644 index 0000000000000000000000000000000000000000..2aa9a39531109f69c9243f222540324078266254 GIT binary patch literal 66 zcmb2|=3oE=W}cH984nn6ux$J_HO09`-A5KQNqUA$(KZr`Pzou Tnx|s=k=uab|NLiWHbAWaTfZ9r literal 0 HcmV?d00001 From 24a1a046683031237eb98fb6a4b74c6a211b45e5 Mon Sep 17 00:00:00 2001 From: Till Hartmann Date: Tue, 2 May 2023 15:51:31 +0200 Subject: [PATCH 10/20] set reference for indexed cram reader --- src/bam/mod.rs | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/bam/mod.rs b/src/bam/mod.rs index dc9d00ae4..76cfebdbf 100644 --- a/src/bam/mod.rs +++ b/src/bam/mod.rs @@ -2918,6 +2918,7 @@ CCCCCCCCCCCCCCCCCCC"[..], #[test] fn test_idxstats_cram() { let mut reader = IndexedReader::from_path("test/test_cram.cram").unwrap(); + reader.set_reference("test/test_cram.fa").unwrap(); let header = reader.header().clone(); let expected = vec![(0, (120, 2, 0)), (1, (120, 2, 0)), (2, (120, 2, 0))]; let actual = reader.index().stats(header).unwrap(); @@ -2927,6 +2928,7 @@ CCCCCCCCCCCCCCCCCCC"[..], #[test] fn test_number_mapped_and_unmapped_cram() { let mut reader = IndexedReader::from_path("test/test_cram.cram").unwrap(); + reader.set_reference("test/test_cram.fa").unwrap(); let expected = (2, 0); let actual = reader.index().number_mapped_unmapped(0); assert_eq!(expected, actual); From 5804bfde4b020bd83fe0cc3ebbf38ad8becea7d4 Mon Sep 17 00:00:00 2001 From: Till Hartmann Date: Fri, 5 May 2023 11:48:32 +0200 Subject: [PATCH 11/20] add slow_idxstats from samtools bam_index.c --- src/bam/mod.rs | 116 ++++++++++++++++++++++++++++++++++++---- test/test_unmapped.bam | Bin 0 -> 962 bytes test/test_unmapped.cram | Bin 0 -> 1165 bytes 3 files changed, 106 insertions(+), 10 deletions(-) create mode 100644 test/test_unmapped.bam create mode 100644 test/test_unmapped.cram diff --git a/src/bam/mod.rs b/src/bam/mod.rs index 76cfebdbf..3f030d946 100644 --- a/src/bam/mod.rs +++ b/src/bam/mod.rs @@ -25,14 +25,15 @@ use std::str; use url::Url; use crate::errors::{Error, Result}; -use crate::htslib; use crate::tpool::ThreadPool; use crate::utils::path_as_bytes; +use crate::{bam, htslib}; pub use crate::bam::buffer::RecordBuffer; pub use crate::bam::header::Header; pub use crate::bam::record::Record; use std::convert::TryInto; +use std::mem::MaybeUninit; /// # Safety /// @@ -831,21 +832,88 @@ impl IndexView { (mapped, unmapped) } + // Analogous to slow_idxstats in samtools, see + // /~https://github.com/samtools/samtools/blob/556c60fdff977c0e6cadc4c2581661f187098b4d/bam_index.c#L140-L199 + // TODO proper error handling, actually return values + unsafe fn slow_idxstats(&self, bam: &mut R) -> Option<(u32, (u64, u64, u64))> { + let header = bam.header(); + let h = header.inner; + let (mut ret, mut last_tid) = (-2, -2); + let fp = bam.htsfile(); + if hts_sys::hts_set_opt( + fp, + hts_sys::hts_fmt_option_CRAM_OPT_REQUIRED_FIELDS, + hts_sys::sam_fields_SAM_RNAME | hts_sys::sam_fields_SAM_FLAG, + ) > 0 + { + panic!("hts_set_opt failed"); + } + + let nref = hts_sys::sam_hdr_nref(h) as usize; + let mut counts = vec![vec![0; 2]; nref + 1]; + if counts.is_empty() { + panic!("idxstats: no sequences in the reference") + } + let mut bb: hts_sys::bam1_t = MaybeUninit::zeroed().assume_init(); + let b = &mut bb as *mut hts_sys::bam1_t; + loop { + ret = hts_sys::sam_read1(fp, h, b); + if ret < 0 { + break; + } + let tid = (*b).core.tid; + if tid >= nref as i32 || tid < -1 { + panic!("idxstats: {} is not a valid tid", tid); + } + + if tid != last_tid { + if (last_tid >= -1) && (counts[tid as usize][0] + counts[tid as usize][1]) > 0 { + panic!("idxstats: file is not position sorted"); + } + last_tid = tid; + } + + let idx = if ((*b).core.flag as u32 & hts_sys::BAM_FUNMAP) > 0 { + 1 + } else { + 0 + }; + counts[(*b).core.tid as usize][idx] += 1; + } + + if ret == -1 { + for i in 0..nref { + println!( + "{:?}\t{}\t{}\t{}", + header.tid2name(i as u32), + header.target_len(i as u32).unwrap(), + counts[i][0], + counts[i][1] + ); + } + println!("*\t0\t{}\t{}", counts[nref][0], counts[nref][1]); + } + None + } + /// Similar to samtools idxstats, this returns a vector of tuples /// containing the target id, length, number of mapped reads, and number of unmapped reads. pub fn stats(&self, header: HeaderView) -> Option> { - header - .target_names() - .iter() - .map(|tname| { - header.tid(tname).and_then(|tid| { - let (mapped, unmapped) = self.number_mapped_unmapped(tid); - let tlen = header.target_len(tid); - tlen.map(|tlen| (tid, (tlen, mapped, unmapped))) - }) + if self.inner.is_null() { + panic!("Index is null"); + } + (0..header.target_count()) + .map(|tid| { + let (mapped, unmapped) = self.number_mapped_unmapped(tid); + let tlen = header.target_len(tid); + tlen.map(|tlen| (tid, (tlen, mapped, unmapped))) }) .collect::>() } + + pub fn number_unmapped(&self) -> u64 { + unsafe { hts_sys::hts_idx_get_n_no_coor(self.inner) } + } } impl Drop for IndexView { @@ -2915,6 +2983,14 @@ CCCCCCCCCCCCCCCCCCC"[..], assert_eq!(expected, actual); } + #[test] + fn test_number_unmapped_global_bam() { + let mut reader = IndexedReader::from_path("test/test_unmapped.bam").unwrap(); + let expected = 8; + let actual = reader.index().number_unmapped(); + assert_eq!(expected, actual); + } + #[test] fn test_idxstats_cram() { let mut reader = IndexedReader::from_path("test/test_cram.cram").unwrap(); @@ -2925,6 +3001,18 @@ CCCCCCCCCCCCCCCCCCC"[..], assert_eq!(expected, actual); } + #[test] + fn test_slow_idxstats_cram() { + let mut reader = IndexedReader::from_path("test/test_cram.cram").unwrap(); + reader.set_reference("test/test_cram.fa").unwrap(); + // let header = reader.header().clone(); + let expected = vec![(0, (120, 2, 0)), (1, (120, 2, 0)), (2, (120, 2, 0))]; + let mut b = IndexedReader::from_path("test/test_cram.cram").unwrap(); + let actual = unsafe { reader.index().slow_idxstats(&mut b).unwrap() }; + // assert_eq!(expected, actual); + panic!() + } + #[test] fn test_number_mapped_and_unmapped_cram() { let mut reader = IndexedReader::from_path("test/test_cram.cram").unwrap(); @@ -2933,4 +3021,12 @@ CCCCCCCCCCCCCCCCCCC"[..], let actual = reader.index().number_mapped_unmapped(0); assert_eq!(expected, actual); } + + #[test] + fn test_number_unmapped_global_cram() { + let mut reader = IndexedReader::from_path("test/test_unmapped.cram").unwrap(); + let expected = 8; + let actual = reader.index().number_unmapped(); + assert_eq!(expected, actual); + } } diff --git a/test/test_unmapped.bam b/test/test_unmapped.bam new file mode 100644 index 0000000000000000000000000000000000000000..aa1ed45d8b36cf40db65c46c0b5fe0db4c054534 GIT binary patch literal 962 zcmV;z13mm7iwFb&00000{{{d;LjnM*07cEs3W7ishT+rlcopXc92gfaW@RZ@P$~v( zW6Bc=`8P8vT~W6ggw^ih*>+RzaU8l@xr6 z|MGE2Zm`-Df*cVO;901J;Q9U-Spw1y@D7CcE1r_S0n0l)LO>Y3*k9hH#fODzV-Pn> zC}09trmSn_3Z{xG0kHfG8}``ZBv89@l6*DWM1B5BV zuu~FM1gwioma=emONM8p>M-Fm!Di;qh+ z=4jbirCTv)HNBrq(spC07-gFF&15uF6lv9z>S=u4%icQy?s1bjne_?J^@}XNrp2eF z@kPMF+9+kOs;|OSE;d@VR(t*Hu!YerU0EHrbq$-Qo8&U>gfSK4%QK)1&ym~X0J#PD zm|_52@^G(A{RX!IH-cmkc@x6{D@aWEHZ49a#H|?^P0Cu^FS0BVMTN2|CYMuD)hK80 z?I5Ern%}6|w;1De7ZQ)z77||g4VDhzJh=VDXp}BP?hnX?1t_JKa3+51XvH k+b{+I03VA81ONa4009360763o02=@U00000000000E!d0N&o-= literal 0 HcmV?d00001 diff --git a/test/test_unmapped.cram b/test/test_unmapped.cram new file mode 100644 index 0000000000000000000000000000000000000000..e0e0525e7ccc2d4da03fafda8cb55e50963b9cce GIT binary patch literal 1165 zcmZ<`a`a_pC`m0Yi7(B|O)Mx#P0>p(O3Yv7Ru9R0{L-qQ9#nc!-X@< z&&p8Gj5FBZsx+@Szo;ZNh11>Ns<1S*sFKUU*`G5gwJf!$I8~uMvn0bxp*S(OBtJi= zSRpO3xTH`)*H9tIP!GyBQV24Fvbnx5o@WBIz}OgOHwkpS!R_;n%nS_wfq*g;+vd7#V_nr9mLbkI^xJv0;LclN%!=zq1=N3nMclLy!-U zPGMwl3;@zVZ9(qL9E@NS85w**YM2-q{DKi;eo%EG!4UN!ApJmno_?%M47^T`jEr(l zPAqJUj9QG0s=k3xLqU?DTRNy{gW#Q)GmxFN8 zLWFGwn=At>P*9)!&qqcUZk+`AZVrf=|Jvuy=$`e|@$|iDoT72;l-Z$!X2%YiIVLhQ z9CqXBN&sov5AN0SKeW^+bjR?XY{)rDd+R=L`6wm zpAj)T^iuq*MPK+e_LhGRIaohoC&R1!e7E*V^Sq9QREGB2m>Dp!vQJ@|zD6?C>{U;% z1MABN_X-2fy-~Z)vL(8}%r5Vg@SZAeKe5UEfmVu}li!|Ds59KfacH`>%oh;!_M$6?%?_>HT>Px`a(|?4(|G7aNhA= zaIXsk3qKE!Z$mE}RWTc9v#`iDyaNXHQ(!2Ag4(e~&bg({y}8SeHKh3rgKCVcj^ppQ z3>m!wdiS<87_Ezq@-JE|_+ww+=~D+^i{uJ_UoPC|pLXebR{zNp%Y7F8tC}i!a<|p5FeGpuG>chIDmnL5J;mOln yblqQnqSH3bMPW?_?@m7F2YMKsd>^>_y8+q29KyQyrazEkV*_OuMuvCEjNSmtGL1w4 literal 0 HcmV?d00001 From 236c18e365aef3d3cc776b7259cee9f9c8074a4e Mon Sep 17 00:00:00 2001 From: Till Hartmann Date: Thu, 11 May 2023 10:34:10 +0200 Subject: [PATCH 12/20] fix slow_idxstats, refactor as methods of IndexedReader --- src/bam/mod.rs | 251 ++++++++++++++++++++++++++----------------------- src/errors.rs | 9 ++ 2 files changed, 144 insertions(+), 116 deletions(-) diff --git a/src/bam/mod.rs b/src/bam/mod.rs index 3f030d946..defda671c 100644 --- a/src/bam/mod.rs +++ b/src/bam/mod.rs @@ -32,7 +32,8 @@ use crate::{bam, htslib}; pub use crate::bam::buffer::RecordBuffer; pub use crate::bam::header::Header; pub use crate::bam::record::Record; -use std::convert::TryInto; +use hts_sys::{hts_fmt_option, sam_fields}; +use std::convert::{TryFrom, TryInto}; use std::mem::MaybeUninit; /// # Safety @@ -782,78 +783,29 @@ impl IndexedReader { unsafe { set_fai_filename(self.htsfile, path) } } - pub fn index(&mut self) -> &IndexView { + pub fn index(&self) -> &IndexView { &self.idx } -} - -#[derive(Debug)] -pub struct IndexView { - inner: *mut hts_sys::hts_idx_t, - owned: bool, -} - -impl IndexView { - fn new(hts_idx: *mut hts_sys::hts_idx_t) -> Self { - Self { - inner: hts_idx, - owned: true, - } - } - - #[inline] - pub fn inner(&self) -> &hts_sys::hts_idx_t { - unsafe { self.inner.as_ref().unwrap() } - } - - #[inline] - // Pointer to inner hts_idx_t struct - pub fn inner_ptr(&self) -> *const hts_sys::hts_idx_t { - self.inner - } - - #[inline] - pub fn inner_mut(&mut self) -> &mut hts_sys::hts_idx_t { - unsafe { self.inner.as_mut().unwrap() } - } - - #[inline] - // Mutable pointer to hts_idx_t struct - pub fn inner_ptr_mut(&mut self) -> *mut hts_sys::hts_idx_t { - self.inner - } - - /// Get the number of mapped and unmapped reads for a given target id - pub fn number_mapped_unmapped(&self, tid: u32) -> (u64, u64) { - let (mut mapped, mut unmapped) = (0, 0); - unsafe { - hts_sys::hts_idx_get_stat(self.inner, tid as i32, &mut mapped, &mut unmapped); - } - (mapped, unmapped) - } // Analogous to slow_idxstats in samtools, see // /~https://github.com/samtools/samtools/blob/556c60fdff977c0e6cadc4c2581661f187098b4d/bam_index.c#L140-L199 // TODO proper error handling, actually return values - unsafe fn slow_idxstats(&self, bam: &mut R) -> Option<(u32, (u64, u64, u64))> { - let header = bam.header(); - let h = header.inner; - let (mut ret, mut last_tid) = (-2, -2); - let fp = bam.htsfile(); - if hts_sys::hts_set_opt( - fp, + unsafe fn slow_idxstats(&mut self) -> Result> { + self.set_read_options( hts_sys::hts_fmt_option_CRAM_OPT_REQUIRED_FIELDS, hts_sys::sam_fields_SAM_RNAME | hts_sys::sam_fields_SAM_FLAG, - ) > 0 - { - panic!("hts_set_opt failed"); - } + )?; + let header = self.header(); + let h = header.inner; + let (mut ret, mut last_tid) = (-2, -2); + let fp = self.htsfile(); - let nref = hts_sys::sam_hdr_nref(h) as usize; - let mut counts = vec![vec![0; 2]; nref + 1]; - if counts.is_empty() { - panic!("idxstats: no sequences in the reference") + let nref = + usize::try_from(hts_sys::sam_hdr_nref(h)).map_err(|_| Error::NoSequencesInReference)?; + if nref == 0 { + return Ok(vec![]); } + let mut counts = vec![vec![0; 2]; nref + 1]; let mut bb: hts_sys::bam1_t = MaybeUninit::zeroed().assume_init(); let b = &mut bb as *mut hts_sys::bam1_t; loop { @@ -863,7 +815,7 @@ impl IndexView { } let tid = (*b).core.tid; if tid >= nref as i32 || tid < -1 { - panic!("idxstats: {} is not a valid tid", tid); + return Err(Error::InvalidTid { tid }); } if tid != last_tid { @@ -882,36 +834,97 @@ impl IndexView { } if ret == -1 { - for i in 0..nref { - println!( - "{:?}\t{}\t{}\t{}", - header.tid2name(i as u32), - header.target_len(i as u32).unwrap(), - counts[i][0], - counts[i][1] - ); - } - println!("*\t0\t{}\t{}", counts[nref][0], counts[nref][1]); + let res = (0..nref) + .map(|i| { + ( + i as i64, + header.target_len(i as u32).unwrap(), + counts[i][0], + counts[i][1], + ) + }) + .chain([(-1, 0, counts[nref][0], counts[nref][1])]) + .collect(); + Ok(res) + } else { + Err(Error::SlowIdxStats) } - None } /// Similar to samtools idxstats, this returns a vector of tuples /// containing the target id, length, number of mapped reads, and number of unmapped reads. - pub fn stats(&self, header: HeaderView) -> Option> { - if self.inner.is_null() { + pub fn index_stats(&mut self) -> Result> { + let header = self.header(); + let index = self.index(); + if index.inner_ptr().is_null() { panic!("Index is null"); } - (0..header.target_count()) + // the quick index stats method only works for BAM files, not SAM or CRAM + unsafe { + if (*self.htsfile()).format.format != htslib::htsExactFormat_bam { + return self.slow_idxstats(); + } + } + Ok((0..header.target_count()) .map(|tid| { - let (mapped, unmapped) = self.number_mapped_unmapped(tid); - let tlen = header.target_len(tid); - tlen.map(|tlen| (tid, (tlen, mapped, unmapped))) + let (mapped, unmapped) = index.number_mapped_unmapped(tid); + let tlen = header.target_len(tid).unwrap(); + (tid as i64, tlen, mapped, unmapped) }) - .collect::>() + .chain([(-1, 0, 0, index.number_unmapped())]) + .collect::<_>()) } +} - pub fn number_unmapped(&self) -> u64 { +#[derive(Debug)] +pub struct IndexView { + inner: *mut hts_sys::hts_idx_t, + owned: bool, +} + +impl IndexView { + fn new(hts_idx: *mut hts_sys::hts_idx_t) -> Self { + Self { + inner: hts_idx, + owned: true, + } + } + + #[inline] + pub fn inner(&self) -> &hts_sys::hts_idx_t { + unsafe { self.inner.as_ref().unwrap() } + } + + #[inline] + // Pointer to inner hts_idx_t struct + pub fn inner_ptr(&self) -> *const hts_sys::hts_idx_t { + self.inner + } + + #[inline] + pub fn inner_mut(&mut self) -> &mut hts_sys::hts_idx_t { + unsafe { self.inner.as_mut().unwrap() } + } + + #[inline] + // Mutable pointer to hts_idx_t struct + pub fn inner_ptr_mut(&mut self) -> *mut hts_sys::hts_idx_t { + self.inner + } + + /// Get the number of mapped and unmapped reads for a given target id + /// FIXME only valid for BAM, not SAM/CRAM + fn number_mapped_unmapped(&self, tid: u32) -> (u64, u64) { + let (mut mapped, mut unmapped) = (0, 0); + unsafe { + hts_sys::hts_idx_get_stat(self.inner, tid as i32, &mut mapped, &mut unmapped); + } + (mapped, unmapped) + } + + /// Get the total number of unmapped reads in the file + /// FIXME only valid for BAM, not SAM/CRAM + fn number_unmapped(&self) -> u64 { unsafe { hts_sys::hts_idx_get_n_no_coor(self.inner) } } } @@ -2963,21 +2976,21 @@ CCCCCCCCCCCCCCCCCCC"[..], #[test] fn test_idxstats_bam() { let mut reader = IndexedReader::from_path("test/test.bam").unwrap(); - let header = reader.header().clone(); let expected = vec![ - (0, (15072423, 6, 0)), - (1, (15279345, 0, 0)), - (2, (13783700, 0, 0)), - (3, (17493793, 0, 0)), - (4, (20924149, 0, 0)), + (0, 15072423, 6, 0), + (1, 15279345, 0, 0), + (2, 13783700, 0, 0), + (3, 17493793, 0, 0), + (4, 20924149, 0, 0), + (-1, 0, 0, 0), ]; - let actual = reader.index().stats(header).unwrap(); + let actual = reader.index_stats().unwrap(); assert_eq!(expected, actual); } #[test] fn test_number_mapped_and_unmapped_bam() { - let mut reader = IndexedReader::from_path("test/test.bam").unwrap(); + let reader = IndexedReader::from_path("test/test.bam").unwrap(); let expected = (6, 0); let actual = reader.index().number_mapped_unmapped(0); assert_eq!(expected, actual); @@ -2985,7 +2998,7 @@ CCCCCCCCCCCCCCCCCCC"[..], #[test] fn test_number_unmapped_global_bam() { - let mut reader = IndexedReader::from_path("test/test_unmapped.bam").unwrap(); + let reader = IndexedReader::from_path("test/test_unmapped.bam").unwrap(); let expected = 8; let actual = reader.index().number_unmapped(); assert_eq!(expected, actual); @@ -2995,9 +3008,13 @@ CCCCCCCCCCCCCCCCCCC"[..], fn test_idxstats_cram() { let mut reader = IndexedReader::from_path("test/test_cram.cram").unwrap(); reader.set_reference("test/test_cram.fa").unwrap(); - let header = reader.header().clone(); - let expected = vec![(0, (120, 2, 0)), (1, (120, 2, 0)), (2, (120, 2, 0))]; - let actual = reader.index().stats(header).unwrap(); + let expected = vec![ + (0, 120, 2, 0), + (1, 120, 2, 0), + (2, 120, 2, 0), + (-1, 0, 0, 0), + ]; + let actual = reader.index_stats().unwrap(); assert_eq!(expected, actual); } @@ -3005,28 +3022,30 @@ CCCCCCCCCCCCCCCCCCC"[..], fn test_slow_idxstats_cram() { let mut reader = IndexedReader::from_path("test/test_cram.cram").unwrap(); reader.set_reference("test/test_cram.fa").unwrap(); - // let header = reader.header().clone(); - let expected = vec![(0, (120, 2, 0)), (1, (120, 2, 0)), (2, (120, 2, 0))]; - let mut b = IndexedReader::from_path("test/test_cram.cram").unwrap(); - let actual = unsafe { reader.index().slow_idxstats(&mut b).unwrap() }; - // assert_eq!(expected, actual); - panic!() - } - - #[test] - fn test_number_mapped_and_unmapped_cram() { - let mut reader = IndexedReader::from_path("test/test_cram.cram").unwrap(); - reader.set_reference("test/test_cram.fa").unwrap(); - let expected = (2, 0); - let actual = reader.index().number_mapped_unmapped(0); + let expected = vec![ + (0, 120, 2, 0), + (1, 120, 2, 0), + (2, 120, 2, 0), + (-1, 0, 0, 0), + ]; + let actual = reader.index_stats().unwrap(); assert_eq!(expected, actual); } - #[test] - fn test_number_unmapped_global_cram() { - let mut reader = IndexedReader::from_path("test/test_unmapped.cram").unwrap(); - let expected = 8; - let actual = reader.index().number_unmapped(); - assert_eq!(expected, actual); - } + // #[test] + // fn test_number_mapped_and_unmapped_cram() { + // let mut reader = IndexedReader::from_path("test/test_cram.cram").unwrap(); + // reader.set_reference("test/test_cram.fa").unwrap(); + // let expected = (2, 0); + // let actual = reader.index().number_mapped_unmapped(0); + // assert_eq!(expected, actual); + // } + // + // #[test] + // fn test_number_unmapped_global_cram() { + // let mut reader = IndexedReader::from_path("test/test_unmapped.cram").unwrap(); + // let expected = 8; + // let actual = reader.index().number_unmapped(); + // assert_eq!(expected, actual); + // } } diff --git a/src/errors.rs b/src/errors.rs index 680bd5f33..9782ae276 100644 --- a/src/errors.rs +++ b/src/errors.rs @@ -118,4 +118,13 @@ pub enum Error { #[error("invalid compression level {level}")] BgzfInvalidCompressionLevel { level: i8 }, + #[error("failed setting hts reading options")] + HtsSetOpt, + #[error("failed calculating slow index statistics")] + SlowIdxStats, + #[error("invalid tid {tid}")] + InvalidTid { tid : i32}, + #[error("No sequences in the reference")] + NoSequencesInReference, } + From 5275afe1cf78359235fa73337d6fbb4a940659b4 Mon Sep 17 00:00:00 2001 From: Till Hartmann Date: Thu, 11 May 2023 10:34:36 +0200 Subject: [PATCH 13/20] Add set_read_options 'convenience' method to bam::Read --- src/bam/mod.rs | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/bam/mod.rs b/src/bam/mod.rs index defda671c..1b635fb06 100644 --- a/src/bam/mod.rs +++ b/src/bam/mod.rs @@ -219,6 +219,16 @@ pub trait Read: Sized { /// /// * `tpool` - thread pool to use for compression work. fn set_thread_pool(&mut self, tpool: &ThreadPool) -> Result<()>; + + fn set_read_options(&mut self, fmt_opt: hts_fmt_option, fields: sam_fields) -> Result<()> { + unsafe { + if hts_sys::hts_set_opt(self.htsfile(), fmt_opt, fields) != 0 { + Err(Error::HtsSetOpt) + } else { + Ok(()) + } + } + } } /// A BAM reader. From cdab234dccd942d62d7c21730d0dc526fe8d047b Mon Sep 17 00:00:00 2001 From: Till Hartmann Date: Thu, 11 May 2023 10:34:57 +0200 Subject: [PATCH 14/20] cargo fmt --- src/errors.rs | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/errors.rs b/src/errors.rs index 9782ae276..fbe50285c 100644 --- a/src/errors.rs +++ b/src/errors.rs @@ -123,8 +123,7 @@ pub enum Error { #[error("failed calculating slow index statistics")] SlowIdxStats, #[error("invalid tid {tid}")] - InvalidTid { tid : i32}, + InvalidTid { tid: i32 }, #[error("No sequences in the reference")] NoSequencesInReference, } - From 4f3d298c9f62ee2c4fbec9c831fb37d2b6821d20 Mon Sep 17 00:00:00 2001 From: Till Hartmann Date: Thu, 11 May 2023 12:06:58 +0200 Subject: [PATCH 15/20] add missing test bam/cram indices --- test/test_cram.bam.bai | Bin 0 -> 256 bytes test/test_unmapped.bam.bai | Bin 0 -> 16 bytes test/test_unmapped.cram.crai | Bin 0 -> 41 bytes 3 files changed, 0 insertions(+), 0 deletions(-) create mode 100644 test/test_cram.bam.bai create mode 100644 test/test_unmapped.bam.bai create mode 100644 test/test_unmapped.cram.crai diff --git a/test/test_cram.bam.bai b/test/test_cram.bam.bai new file mode 100644 index 0000000000000000000000000000000000000000..437840c4fbc2894a15229500c42f7d2a62c6f6a3 GIT binary patch literal 256 zcmZ>A^kigaU|?VZVoxCk21X#wz`zNnwm@mGNf3E-agZKlfUFnUJeWDfj1Y5S=AnyY cGY_T)**usz3=BMAEim)Y#j%-(Py=x}06NkN^Z)<= literal 0 HcmV?d00001 diff --git a/test/test_unmapped.bam.bai b/test/test_unmapped.bam.bai new file mode 100644 index 0000000000000000000000000000000000000000..3bbdf839ebe336ced8d2966743fa3246d0d57ae8 GIT binary patch literal 16 RcmZ>A^kifJ0uB%X1ON<50M-Bi literal 0 HcmV?d00001 diff --git a/test/test_unmapped.cram.crai b/test/test_unmapped.cram.crai new file mode 100644 index 0000000000000000000000000000000000000000..ac8c91848f8425ccdbaf7c3c22de1dfe95a00289 GIT binary patch literal 41 tcmb2|=3oE==1ZneOdgm#F)=cDZep_ifr Date: Thu, 11 May 2023 12:20:59 +0200 Subject: [PATCH 16/20] rename to set_cram_options, add doc --- src/bam/mod.rs | 24 ++++++++++++++++++++---- 1 file changed, 20 insertions(+), 4 deletions(-) diff --git a/src/bam/mod.rs b/src/bam/mod.rs index 1b635fb06..c8a05efe1 100644 --- a/src/bam/mod.rs +++ b/src/bam/mod.rs @@ -25,9 +25,9 @@ use std::str; use url::Url; use crate::errors::{Error, Result}; +use crate::htslib; use crate::tpool::ThreadPool; use crate::utils::path_as_bytes; -use crate::{bam, htslib}; pub use crate::bam::buffer::RecordBuffer; pub use crate::bam::header::Header; @@ -220,7 +220,22 @@ pub trait Read: Sized { /// * `tpool` - thread pool to use for compression work. fn set_thread_pool(&mut self, tpool: &ThreadPool) -> Result<()>; - fn set_read_options(&mut self, fmt_opt: hts_fmt_option, fields: sam_fields) -> Result<()> { + /// If the underlying file is in CRAM format, allows modifying CRAM options. + /// Note that this method does *not* check that the underlying file actually is in CRAM format. + /// + /// # Examples + /// + /// Set the required fields to RNAME and FLAG, + /// potentially allowing htslib to skip over the rest, + /// resulting in faster iteration: + /// ``` + /// use rust_htslib::bam::{Read, Reader}; + /// use hts_sys; + /// let mut cram = Reader::from_path(&"test/test.cram").unwrap(); + /// cram.set_cram_options(hts_sys::hts_fmt_option_CRAM_OPT_REQUIRED_FIELDS, + /// hts_sys::sam_fields_SAM_RNAME | hts_sys::sam_fields_SAM_FLAG).unwrap(); + /// ``` + fn set_cram_options(&mut self, fmt_opt: hts_fmt_option, fields: sam_fields) -> Result<()> { unsafe { if hts_sys::hts_set_opt(self.htsfile(), fmt_opt, fields) != 0 { Err(Error::HtsSetOpt) @@ -801,13 +816,14 @@ impl IndexedReader { // /~https://github.com/samtools/samtools/blob/556c60fdff977c0e6cadc4c2581661f187098b4d/bam_index.c#L140-L199 // TODO proper error handling, actually return values unsafe fn slow_idxstats(&mut self) -> Result> { - self.set_read_options( + self.set_cram_options( hts_sys::hts_fmt_option_CRAM_OPT_REQUIRED_FIELDS, hts_sys::sam_fields_SAM_RNAME | hts_sys::sam_fields_SAM_FLAG, )?; let header = self.header(); let h = header.inner; - let (mut ret, mut last_tid) = (-2, -2); + let mut ret; + let mut last_tid = -2; let fp = self.htsfile(); let nref = From 9475c999a919958ecc822a024710dea4fe56ded7 Mon Sep 17 00:00:00 2001 From: Till Hartmann Date: Thu, 11 May 2023 12:25:09 +0200 Subject: [PATCH 17/20] replace panic with BamUnsorted error --- src/bam/mod.rs | 2 +- src/errors.rs | 2 ++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/src/bam/mod.rs b/src/bam/mod.rs index c8a05efe1..65d8f174e 100644 --- a/src/bam/mod.rs +++ b/src/bam/mod.rs @@ -846,7 +846,7 @@ impl IndexedReader { if tid != last_tid { if (last_tid >= -1) && (counts[tid as usize][0] + counts[tid as usize][1]) > 0 { - panic!("idxstats: file is not position sorted"); + return Err(Error::BamUnsorted); } last_tid = tid; } diff --git a/src/errors.rs b/src/errors.rs index fbe50285c..18e3d6b5f 100644 --- a/src/errors.rs +++ b/src/errors.rs @@ -71,6 +71,8 @@ pub enum Error { BamBuildIndex, #[error("failed to create SAM/BAM/CRAM pileup")] BamPileup, + #[error("file is not sorted by position")] + BamUnsorted, // Errors for BAM auxiliary fields #[error("failed to add aux field (out of memory?)")] From bf42d1ea4a6139bb166e4df5788193e7b04f16bf Mon Sep 17 00:00:00 2001 From: Till Hartmann Date: Thu, 11 May 2023 12:26:59 +0200 Subject: [PATCH 18/20] Add tid == -1 comment --- src/bam/mod.rs | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/bam/mod.rs b/src/bam/mod.rs index 65d8f174e..d5bc63a9b 100644 --- a/src/bam/mod.rs +++ b/src/bam/mod.rs @@ -879,6 +879,8 @@ impl IndexedReader { /// Similar to samtools idxstats, this returns a vector of tuples /// containing the target id, length, number of mapped reads, and number of unmapped reads. + /// The last entry in the vector corresponds to the unmapped reads for the entire file, with + /// the tid set to -1. pub fn index_stats(&mut self) -> Result> { let header = self.header(); let index = self.index(); From babb711fbc2e8e3cc8cab68c8c9c4d8b17f98ed7 Mon Sep 17 00:00:00 2001 From: Till Hartmann Date: Thu, 11 May 2023 12:27:13 +0200 Subject: [PATCH 19/20] remove todo comment --- src/bam/mod.rs | 1 - 1 file changed, 1 deletion(-) diff --git a/src/bam/mod.rs b/src/bam/mod.rs index d5bc63a9b..afdf67837 100644 --- a/src/bam/mod.rs +++ b/src/bam/mod.rs @@ -814,7 +814,6 @@ impl IndexedReader { // Analogous to slow_idxstats in samtools, see // /~https://github.com/samtools/samtools/blob/556c60fdff977c0e6cadc4c2581661f187098b4d/bam_index.c#L140-L199 - // TODO proper error handling, actually return values unsafe fn slow_idxstats(&mut self) -> Result> { self.set_cram_options( hts_sys::hts_fmt_option_CRAM_OPT_REQUIRED_FIELDS, From 36c5cf03e70292812cdc00a78b1d898f9a7c97ca Mon Sep 17 00:00:00 2001 From: Till Hartmann Date: Thu, 11 May 2023 12:32:04 +0200 Subject: [PATCH 20/20] fix doctest path to file --- src/bam/mod.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/bam/mod.rs b/src/bam/mod.rs index afdf67837..dae37c511 100644 --- a/src/bam/mod.rs +++ b/src/bam/mod.rs @@ -231,7 +231,7 @@ pub trait Read: Sized { /// ``` /// use rust_htslib::bam::{Read, Reader}; /// use hts_sys; - /// let mut cram = Reader::from_path(&"test/test.cram").unwrap(); + /// let mut cram = Reader::from_path(&"test/test_cram.cram").unwrap(); /// cram.set_cram_options(hts_sys::hts_fmt_option_CRAM_OPT_REQUIRED_FIELDS, /// hts_sys::sam_fields_SAM_RNAME | hts_sys::sam_fields_SAM_FLAG).unwrap(); /// ```