Skip to content

Commit

Permalink
feat: add method to retrieve all event posteriors from Bayesian model.
Browse files Browse the repository at this point in the history
  • Loading branch information
johanneskoester committed Jul 9, 2021
1 parent 5853bda commit e5feda8
Show file tree
Hide file tree
Showing 3 changed files with 36 additions and 16 deletions.
12 changes: 10 additions & 2 deletions src/data_structures/fmindex.rs
Expand Up @@ -382,13 +382,21 @@ impl<DBWT: Borrow<BWT>, DLess: Borrow<Less>, DOcc: Borrow<Occ>> FMDIndex<DBWT, D
/// let intervals = fmdindex.all_smems(pattern, 0);
/// assert_eq!(intervals.len(), 2);
///
/// let solutions = vec![[0,14,0,3], [4,9,3,4]];
/// let solutions = vec![[0, 14, 0, 3], [4, 9, 3, 4]];
/// for (i, interval) in intervals.iter().enumerate() {
/// let forward_positions = interval.0.forward().occ(&sa);
/// let revcomp_positions = interval.0.revcomp().occ(&sa);
/// let pattern_position = interval.1;
/// let smem_len = interval.2;
/// assert_eq!([forward_positions[0], revcomp_positions[0], pattern_position, smem_len], solutions[i]);
/// assert_eq!(
/// [
/// forward_positions[0],
/// revcomp_positions[0],
/// pattern_position,
/// smem_len
/// ],
/// solutions[i]
/// );
/// }
/// ```
pub fn all_smems(&self, pattern: &[u8], l: usize) -> Vec<(BiInterval, usize, usize)> {
Expand Down
32 changes: 18 additions & 14 deletions src/io/fasta.rs
Expand Up @@ -116,12 +116,14 @@
//! const FAI_FILE: &[u8] = b"chr1\t16\t6\t12\t13";
//!
//! let seq_name = "chr1";
//! let start: u64 = 0; // start is 0-based, inclusive
//! let stop: u64 = 10; // stop is 0-based, exclusive
//! // load the index
//! let start: u64 = 0; // start is 0-based, inclusive
//! let stop: u64 = 10; // stop is 0-based, exclusive
//! // load the index
//! let mut faidx = IndexedReader::new(std::io::Cursor::new(FASTA_FILE), FAI_FILE).unwrap();
//! // move the pointer in the index to the desired sequence and interval
//! faidx.fetch(seq_name, start, stop).expect("Couldn't fetch interval");
//! faidx
//! .fetch(seq_name, start, stop)
//! .expect("Couldn't fetch interval");
//! // read the subsequence defined by the interval into a vector
//! let mut seq = Vec::new();
//! faidx.read(&mut seq).expect("Couldn't read the interval");
Expand Down Expand Up @@ -458,12 +460,14 @@ impl<R: io::Read + io::Seek> IndexedReader<R> {
/// const FAI_FILE: &[u8] = b"chr1\t16\t6\t12\t13";
///
/// let seq_name = "chr1";
/// let start: u64 = 0; // start is 0-based, inclusive
/// let stop: u64 = 10; // stop is 0-based, exclusive
/// // load the index
/// let start: u64 = 0; // start is 0-based, inclusive
/// let stop: u64 = 10; // stop is 0-based, exclusive
/// // load the index
/// let mut faidx = IndexedReader::new(std::io::Cursor::new(FASTA_FILE), FAI_FILE).unwrap();
/// // move the pointer in the index to the desired sequence and interval
/// faidx.fetch(seq_name, start, stop).expect("Couldn't fetch interval");
/// faidx
/// .fetch(seq_name, start, stop)
/// .expect("Couldn't fetch interval");
/// // read the subsequence defined by the interval into a vector
/// let mut seq = Vec::new();
/// faidx.read(&mut seq).expect("Couldn't read the interval");
Expand All @@ -472,7 +476,6 @@ impl<R: io::Read + io::Seek> IndexedReader<R> {
///
/// # Errors
/// If the `seq_name` does not exist within the index.
///
pub fn fetch(&mut self, seq_name: &str, start: u64, stop: u64) -> io::Result<()> {
let idx = self.idx(seq_name)?;
self.start = Some(start);
Expand All @@ -494,12 +497,14 @@ impl<R: io::Read + io::Seek> IndexedReader<R> {
/// const FAI_FILE: &[u8] = b"chr1\t16\t6\t12\t13";
///
/// let rid: usize = 0;
/// let start: u64 = 0; // start is 0-based, inclusive
/// let stop: u64 = 10; // stop is 0-based, exclusive
/// // load the index
/// let start: u64 = 0; // start is 0-based, inclusive
/// let stop: u64 = 10; // stop is 0-based, exclusive
/// // load the index
/// let mut faidx = IndexedReader::new(std::io::Cursor::new(FASTA_FILE), FAI_FILE).unwrap();
/// // move the pointer in the index to the desired sequence and interval
/// faidx.fetch_by_rid(rid, start, stop).expect("Couldn't fetch interval");
/// faidx
/// .fetch_by_rid(rid, start, stop)
/// .expect("Couldn't fetch interval");
/// // read the subsequence defined by the interval into a vector
/// let mut seq = Vec::new();
/// faidx.read(&mut seq).expect("Couldn't read the interval");
Expand All @@ -508,7 +513,6 @@ impl<R: io::Read + io::Seek> IndexedReader<R> {
///
/// # Errors
/// If `rid` does not exist within the index.
///
pub fn fetch_by_rid(&mut self, rid: usize, start: u64, stop: u64) -> io::Result<()> {
let idx = self.idx_by_rid(rid)?;
self.start = Some(start);
Expand Down
8 changes: 8 additions & 0 deletions src/stats/bayesian/model.rs
Expand Up @@ -166,6 +166,14 @@ where
.max_by_key(|(_, prob)| NotNan::new(***prob).unwrap())
.map(|(event, _)| event)
}

/// Event posteriors sorted in descending order.
pub fn event_posteriors(&self) -> impl Iterator<Item = (&Event, LogProb)> {
self.joint_probs
.iter()
.map(|(event, prob)| (event, prob - self.marginal))
.sorted_by_key(|(_, prob)| -NotNan::new(**prob).unwrap())
}
}

impl<PosteriorEvent> ModelInstance<NotNan<f64>, PosteriorEvent>
Expand Down

0 comments on commit e5feda8

Please sign in to comment.