extended_htslib/faidx/
mod.rs

1// Copyright 2020 Manuel Landesfeind, Evotec International GmbH
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//!
7//! Module for working with faidx-indexed FASTA files.
8//!
9
10use std::ffi;
11use std::path::{Path, PathBuf};
12use std::ptr::null;
13use std::str::FromStr;
14use url::Url;
15
16use crate::{errors, htslib};
17
18use crate::errors::{Error, Result};
19use crate::utils::path_as_bytes;
20
21/// A Fasta reader.
22#[derive(Debug)]
23pub struct Reader {
24    inner: *mut htslib::faidx_t,
25}
26
27///
28/// Build a faidx for input path.
29///
30/// # Errors
31/// If indexing fails. Could be malformatted or file could not be accessible.
32///
33///```
34/// use extended_htslib::faidx::build;
35/// let path = std::path::PathBuf::from(concat!(env!("CARGO_MANIFEST_DIR"),"/test/test_cram.fa"));
36/// build(&path).expect("Failed to build fasta index");
37///```
38///
39pub fn build(
40    path: impl Into<std::path::PathBuf>,
41) -> Result<(), std::boxed::Box<dyn std::error::Error>> {
42    let path = path.into();
43    let os_path = std::ffi::CString::new(path.display().to_string())?;
44    let rc = unsafe { htslib::fai_build(os_path.as_ptr()) };
45    if rc < 0 {
46        Err(Error::FaidxBuildFailed { path })?
47    } else {
48        Ok(())
49    }
50}
51
52impl Reader {
53    /// Create a new Reader from a path.
54    ///
55    /// # Arguments
56    ///
57    /// * `path` - the path to open.
58    pub fn from_path<P: AsRef<Path>>(path: P) -> Result<Self, Error> {
59        Self::new(&path_as_bytes(path, true)?)
60    }
61
62    /// Create a new Reader from an URL.
63    ///
64    /// # Arguments
65    ///
66    /// * `url` - the url to open
67    pub fn from_url(url: &Url) -> Result<Self, Error> {
68        Self::new(url.as_str().as_bytes())
69    }
70    /// Internal function to create a Reader from some sort of path and index (could be file path but also URL).
71    /// The path or URL will be handled by the c-implementation transparently.
72    ///
73    /// # Arguments
74    ///
75    /// * `path` - the path or URL to open
76    pub fn from_path_and_index<P: AsRef<Path>>(path: P, index: P) -> Result<Self, Error> {
77        let cpath = ffi::CString::new(path_as_bytes(&path, true)?).unwrap();
78        let cindex = ffi::CString::new(path_as_bytes(index, true)?).unwrap();
79        let inner = unsafe { htslib::fai_load3(cpath.as_ptr(), cindex.as_ptr(), null(), 0) };
80        if inner.is_null() {
81            return Err(errors::Error::FaidxBuildFailed {
82                path: path.as_ref().to_path_buf(),
83            });
84        }
85        Ok(Self { inner })
86    }
87    /// Internal function to create a Reader from some sort of path (could be file path but also URL).
88    /// The path or URL will be handled by the c-implementation transparently.
89    ///
90    /// # Arguments
91    ///
92    /// * `path` - the path or URL to open
93    fn new(path: &[u8]) -> Result<Self, Error> {
94        let cpath = ffi::CString::new(path).unwrap();
95        let inner = unsafe { htslib::fai_load(cpath.as_ptr()) };
96        if inner.is_null() {
97            return Err(errors::Error::FaidxBuildFailed {
98                path: PathBuf::from_str(&String::from_utf8_lossy(path)).unwrap_or(PathBuf::new()),
99            });
100        }
101        Ok(Self { inner })
102    }
103
104    /// Fetch the sequence as a byte array.
105    ///
106    /// # Arguments
107    ///
108    /// * `name` - the name of the template sequence (e.g., "chr1")
109    /// * `begin` - the offset within the template sequence (starting with 0)
110    /// * `end` - the end position to return (if smaller than `begin`, the behavior is undefined).
111    pub fn fetch_seq<N: AsRef<str>>(&self, name: N, begin: usize, end: usize) -> Result<Vec<u8>> {
112        if begin > i64::MAX as usize {
113            return Err(Error::FaidxPositionTooLarge);
114        }
115        if end > i64::MAX as usize {
116            return Err(Error::FaidxPositionTooLarge);
117        }
118        let cname = ffi::CString::new(name.as_ref().as_bytes()).unwrap();
119        let mut len_out: htslib::hts_pos_t = 0;
120        let ptr = unsafe {
121            htslib::faidx_fetch_seq64(
122                self.inner,                 //*const faidx_t,
123                cname.as_ptr(),             // c_name
124                begin as htslib::hts_pos_t, // p_beg_i
125                end as htslib::hts_pos_t,   // p_end_i
126                &mut len_out,               //len
127            )
128        };
129        let vec =
130            unsafe { Vec::from_raw_parts(ptr as *mut u8, len_out as usize, len_out as usize) };
131        Ok(vec)
132    }
133
134    /// Fetches the sequence and returns it as string.
135    ///
136    /// # Arguments
137    ///
138    /// * `name` - the name of the template sequence (e.g., "chr1")
139    /// * `begin` - the offset within the template sequence (starting with 0)
140    /// * `end` - the end position to return (if smaller than `begin`, the behavior is undefined).
141    pub fn fetch_seq_string<N: AsRef<str>>(
142        &self,
143        name: N,
144        begin: usize,
145        end: usize,
146    ) -> Result<String> {
147        let bytes = self.fetch_seq(name, begin, end)?;
148        Ok(std::str::from_utf8(&bytes).unwrap().to_owned())
149    }
150
151    /// Fetches the number of sequences in the fai index
152    pub fn n_seqs(&self) -> u64 {
153        let n = unsafe { htslib::faidx_nseq(self.inner) };
154        n as u64
155    }
156
157    /// Fetches the i-th sequence name
158    ///
159    /// # Arguments
160    ///
161    /// * `i` - index to query
162    pub fn seq_name(&self, i: i32) -> Result<String> {
163        let cname = unsafe {
164            let ptr = htslib::faidx_iseq(self.inner, i);
165            ffi::CStr::from_ptr(ptr)
166        };
167
168        let out = match cname.to_str() {
169            Ok(s) => s.to_string(),
170            Err(_) => {
171                return Err(Error::FaidxBadSeqName);
172            }
173        };
174
175        Ok(out)
176    }
177
178    /// Fetches the length of the given sequence name.
179    ///
180    /// # Arguments
181    ///
182    /// * `name` - the name of the template sequence (e.g., "chr1")
183    pub fn fetch_seq_len<N: AsRef<str>>(&self, name: N) -> u64 {
184        let cname = ffi::CString::new(name.as_ref().as_bytes()).unwrap();
185        let seq_len = unsafe { htslib::faidx_seq_len(self.inner, cname.as_ptr()) };
186        seq_len as u64
187    }
188
189    /// Returns a Result<Vector<String>> for all seq names.
190    /// # Errors
191    ///
192    /// * `errors::Error::FaidxBadSeqName` - missing sequence name for sequence id.
193    ///
194    /// If thrown, the index is malformed, and the number of sequences in the index does not match the number of sequence names available.
195    ///```
196    /// use extended_htslib::faidx::build;
197    /// let path = std::path::PathBuf::from(concat!(env!("CARGO_MANIFEST_DIR"),"/test/test_cram.fa"));
198    /// build(&path).expect("Failed to build fasta index");
199    /// let reader = extended_htslib::faidx::Reader::from_path(path).expect("Failed to open faidx");
200    /// assert_eq!(reader.seq_names(), Ok(vec!["chr1".to_string(), "chr2".to_string(), "chr3".to_string()]));
201    ///```
202    ///
203    pub fn seq_names(&self) -> Result<Vec<String>> {
204        let num_seq = self.n_seqs();
205        let mut ret = Vec::with_capacity(num_seq as usize);
206        for seq_id in 0..num_seq {
207            ret.push(self.seq_name(seq_id as i32)?);
208        }
209        Ok(ret)
210    }
211}
212
213impl Drop for Reader {
214    fn drop(&mut self) {
215        unsafe {
216            htslib::fai_destroy(self.inner);
217        }
218    }
219}
220
221unsafe impl Send for Reader {}
222
223#[cfg(test)]
224mod tests {
225    use super::*;
226
227    fn open_reader() -> Reader {
228        Reader::from_path(format!("{}/test/test_cram.fa", env!("CARGO_MANIFEST_DIR")))
229            .ok()
230            .unwrap()
231    }
232    #[test]
233    fn faidx_open() {
234        open_reader();
235    }
236
237    #[test]
238    fn faidx_read_chr_first_base() {
239        let r = open_reader();
240
241        let bseq = r.fetch_seq("chr1", 0, 0).unwrap();
242        assert_eq!(bseq.len(), 1);
243        assert_eq!(bseq, b"G");
244
245        let seq = r.fetch_seq_string("chr1", 0, 0).unwrap();
246        assert_eq!(seq.len(), 1);
247        assert_eq!(seq, "G");
248    }
249
250    #[test]
251    fn faidx_read_chr_start() {
252        let r = open_reader();
253
254        //for _i in 0..100_000_000 { // loop to check for memory leaks
255        let bseq = r.fetch_seq("chr1", 0, 9).unwrap();
256        assert_eq!(bseq.len(), 10);
257        assert_eq!(bseq, b"GGGCACAGCC");
258
259        let seq = r.fetch_seq_string("chr1", 0, 9).unwrap();
260        assert_eq!(seq.len(), 10);
261        assert_eq!(seq, "GGGCACAGCC");
262        //}
263    }
264
265    #[test]
266    fn faidx_read_chr_between() {
267        let r = open_reader();
268
269        let bseq = r.fetch_seq("chr1", 4, 14).unwrap();
270        assert_eq!(bseq.len(), 11);
271        assert_eq!(bseq, b"ACAGCCTCACC");
272
273        let seq = r.fetch_seq_string("chr1", 4, 14).unwrap();
274        assert_eq!(seq.len(), 11);
275        assert_eq!(seq, "ACAGCCTCACC");
276    }
277
278    #[test]
279    fn faidx_read_chr_end() {
280        let r = open_reader();
281
282        let bseq = r.fetch_seq("chr1", 110, 120).unwrap();
283        assert_eq!(bseq.len(), 10);
284        assert_eq!(bseq, b"CCCCTCCGTG");
285
286        let seq = r.fetch_seq_string("chr1", 110, 120).unwrap();
287        assert_eq!(seq.len(), 10);
288        assert_eq!(seq, "CCCCTCCGTG");
289    }
290
291    #[test]
292    fn faidx_read_twice_string() {
293        let r = open_reader();
294        let seq = r.fetch_seq_string("chr1", 110, 120).unwrap();
295        assert_eq!(seq.len(), 10);
296        assert_eq!(seq, "CCCCTCCGTG");
297
298        let seq = r.fetch_seq_string("chr1", 5, 9).unwrap();
299        assert_eq!(seq.len(), 5);
300        assert_eq!(seq, "CAGCC");
301    }
302
303    #[test]
304    fn faidx_read_twice_bytes() {
305        let r = open_reader();
306        let seq = r.fetch_seq("chr1", 110, 120).unwrap();
307        assert_eq!(seq.len(), 10);
308        assert_eq!(seq, b"CCCCTCCGTG");
309
310        let seq = r.fetch_seq("chr1", 5, 9).unwrap();
311        assert_eq!(seq.len(), 5);
312        assert_eq!(seq, b"CAGCC");
313    }
314
315    #[test]
316    fn faidx_position_too_large() {
317        let r = open_reader();
318        let position_too_large = i64::MAX as usize;
319        let res = r.fetch_seq("chr1", position_too_large, position_too_large + 1);
320        assert_eq!(res, Err(Error::FaidxPositionTooLarge));
321    }
322
323    #[test]
324    fn faidx_n_seqs() {
325        let r = open_reader();
326        assert_eq!(r.n_seqs(), 3);
327    }
328
329    #[test]
330    fn faidx_seq_name() {
331        let r = open_reader();
332        let n = r.seq_name(1).unwrap();
333        assert_eq!(n, "chr2");
334    }
335
336    #[test]
337    fn faidx_get_seq_len() {
338        let r = open_reader();
339        let chr1_len = r.fetch_seq_len("chr1");
340        let chr2_len = r.fetch_seq_len("chr2");
341        assert_eq!(chr1_len, 120u64);
342        assert_eq!(chr2_len, 120u64);
343    }
344
345    #[test]
346    fn open_many_readers() {
347        for _ in 0..500_000 {
348            let reader = open_reader();
349            drop(reader);
350        }
351    }
352}
extended_htslib/faidx/mod.rs

extended_htslib/faidx/
mod.rs
