Struct CurrRead 

pub struct CurrRead<S: CurrReadState> { /* private fields */ }

Our main struct that receives and stores from one BAM record. Also has methods for processing this information. Also see CurrReadBuilder for a way to build this without BAM records.

We call this CurrRead as in ‘current read’. Read is used within rust-htslib, so we don’t want to create another Read struct.

The information within the struct is hard to access without the methods defined here. This is to ensure the struct doesn’t fall into an invalid state, which could cause mistakes in calculations associated with the struct. For example: if I want to measure mean modification density along windows of the raw modification data, I need a guarantee that the modification data is sorted by position. We can guarantee this when the modification data is parsed, but we cannot guarantee this if we allow free access to the struct. NOTE: we could have implemented these as a trait extension to the rust_htslib Record struct, but we have chosen not to, as we may want to persist data like modifications and do multiple operations on them. And Record has inconvenient return types like i64 instead of u64 for positions along the genome.

Implementations

impl CurrRead<NoData>

pub fn set_read_state_and_id( self, record: &Record, ) -> Result<CurrRead<OnlyAlignData>, Error>

sets the alignment of the read and read ID using BAM record

Errors

While we support normal BAM reads from ONT, PacBio etc. that contain modifications, we do not support some BAM flags like paired, duplicate, quality check failed etc. This is because of our design choices e.g. if mods are called on paired reads, then we’ll have to include both records as one read in our statistics and we do not have functionality in place to do this. So, we return errors if such flags or an invalid combination of flags (e.g. secondary and supplementary bits are set) are encountered. We also return error upon invalid read id but this is expected to be in violation of the BAM format (UTF-8 error).

use nanalogue_core::{CurrRead, Error, nanalogue_bam_reader, ReadState};
use rust_htslib::bam::Read;
let mut reader = nanalogue_bam_reader(&"examples/example_1.bam")?;
let mut count = 0;
for record in reader.records(){
    let r = record?;
    let curr_read = CurrRead::default().set_read_state_and_id(&r)?;
    match count {
        0 => assert_eq!(curr_read.read_state(), ReadState::PrimaryFwd),
        1 => assert_eq!(curr_read.read_state(), ReadState::PrimaryFwd),
        2 => assert_eq!(curr_read.read_state(), ReadState::PrimaryRev),
        3 => assert_eq!(curr_read.read_state(), ReadState::Unmapped),
        _ => unreachable!(),
    }
    count = count + 1;
}

pub fn try_from_only_alignment( self, record: &Record, ) -> Result<CurrRead<OnlyAlignDataComplete>, Error>

Runs CurrRead<NoData>::try_from_only_alignment_seq_len_optional, forcing sequence length retrieval. See comments there and the comments below.

Some BAM records have zero-length sequence fields i.e. marked by a ‘*’. This may be intentional e.g. a secondary alignment has the same sequence as a corresponding primary alignment and repeating a sequence is not space-efficient. Or it may be intentional due to other reasons. For modification processing, we cannot deal with these records as we have to match sequences across different records. So, our easy-to-use-function here forces sequence length retrieval and will fail if zero length sequences are found.

Another function below CurrRead<NoData>::try_from_only_alignment_zero_seq_len, allows zero length sequences through and can be used if zero length sequences are really needed. In such a scenario, the user has to carefully watch for errors. So we discourage its use unless really necessary.

Errors

pub fn try_from_only_alignment_zero_seq_len( self, record: &Record, ) -> Result<CurrRead<OnlyAlignDataComplete>, Error>

Runs CurrRead<NoData>::try_from_only_alignment_seq_len_optional, avoiding sequence length retrieval and setting it to zero. See comments there and the comments below.

See notes on CurrRead<NoData>::try_from_only_alignment. Use of this function is discouraged unless really necessary as we cannot parse modification information from zero-length sequences without errors.

Errors

pub fn try_from_only_alignment_seq_len_optional( self, record: &Record, is_seq_len_non_zero: bool, ) -> Result<CurrRead<OnlyAlignDataComplete>, Error>

Uses only alignment information and no modification information to create the struct. Use this if you want to perform operations that do not involve reading or manipulating the modification data. If is_seq_len_non_zero is set to false, then sequence length is not retrieved and is set to zero.

Errors

Upon failure in retrieving record information.

impl<S: CurrReadStateWithAlign + CurrReadState> CurrRead<S>

pub fn read_state(&self) -> ReadState

gets the read state

pub fn set_seq_len(self, record: &Record) -> Result<Self, Error>

set length of sequence from BAM record

Errors

Errors are returned if sequence length is already set or sequence length is not non-zero.

use nanalogue_core::{CurrRead, Error, nanalogue_bam_reader};
use rust_htslib::bam::Read;
let mut reader = nanalogue_bam_reader(&"examples/example_1.bam")?;
let mut count = 0;
for record in reader.records(){
    let r = record?;
    let curr_read = CurrRead::default().set_read_state_and_id(&r)?.set_seq_len(&r)?;
    let Ok(len) = curr_read.seq_len() else { unreachable!() };
    match count {
        0 => assert_eq!(len, 8),
        1 => assert_eq!(len, 48),
        2 => assert_eq!(len, 33),
        3 => assert_eq!(len, 48),
        _ => unreachable!(),
    }
    count = count + 1;
}

If we call the method twice, we should hit a panic

let mut reader = nanalogue_bam_reader(&"examples/example_1.bam")?;
for record in reader.records(){
    let r = record?;
    let curr_read = CurrRead::default().set_read_state_and_id(&r)?
        .set_seq_len(&r)?.set_seq_len(&r)?;
    break;
}

pub fn seq_len(&self) -> Result<u64, Error>

gets length of sequence

Errors

Error if sequence length is not set

pub fn set_align_len(self, record: &Record) -> Result<Self, Error>

set alignment length from BAM record if available

Errors

Returns errors if alignment len is already set, instance is unmapped, or if alignment coordinates are malformed (e.g. end < start).

pub fn align_len(&self) -> Result<u64, Error>

gets alignment length

Errors

if instance is unmapped or alignment length is not set

pub fn set_contig_id_and_start(self, record: &Record) -> Result<Self, Error>

sets contig ID and start from BAM record if available

Errors

if instance is unmapped, if these data are already set and the user is trying to set them again, or if coordinates are malformed (start position is negative)

use nanalogue_core::{CurrRead, Error, nanalogue_bam_reader};
use rust_htslib::bam::Read;
let mut reader = nanalogue_bam_reader(&"examples/example_1.bam")?;
let mut count = 0;
for record in reader.records(){
    let r = record?;
    let curr_read =
        CurrRead::default().set_read_state_and_id(&r)?.set_contig_id_and_start(&r)?;
    match (count, curr_read.contig_id_and_start()) {
        (0, Ok((0, 9))) |
        (1, Ok((2, 23))) |
        (2, Ok((1, 3))) => {},
        _ => unreachable!(),
    }
    count = count + 1;
    if count == 3 { break; } // the fourth entry is unmapped, and will lead to an error.
}

If we call the method on an unmapped read, we should see an error.

let mut reader = nanalogue_bam_reader(&"examples/example_1.bam")?;
let mut count = 0;
for record in reader.records(){
    let r = record?;
    if count < 3 {
        count = count + 1;
        continue;
    }
    // the fourth read is unmapped
    let curr_read =
        CurrRead::default().set_read_state_and_id(&r)?.set_contig_id_and_start(&r)?;
}

If we call the method twice, we should hit a panic

let mut reader = nanalogue_bam_reader(&"examples/example_1.bam")?;
for record in reader.records(){
    let r = record?;
    let mut curr_read = CurrRead::default().set_read_state_and_id(&r)?
        .set_contig_id_and_start(&r)?.set_contig_id_and_start(&r)?;
    break;
}

pub fn contig_id_and_start(&self) -> Result<(i32, u64), Error>

gets contig ID and start

Errors

If instance is unmapped or if data (contig id and start) are not set

pub fn set_contig_name(self, record: &Record) -> Result<Self, Error>

sets contig name

Errors

Returns error if instance is unmapped or contig name has already been set

use nanalogue_core::{CurrRead, Error, nanalogue_bam_reader};
use rust_htslib::bam::Read;
let mut reader = nanalogue_bam_reader(&"examples/example_1.bam")?;
let mut count = 0;
for record in reader.records(){
    let r = record?;
    let mut curr_read = CurrRead::default().set_read_state_and_id(&r)?.set_contig_name(&r)?;
    let Ok(contig_name) = curr_read.contig_name() else {unreachable!()};
    match (count, contig_name) {
        (0, "dummyI") |
        (1, "dummyIII") |
        (2, "dummyII") => {},
        _ => unreachable!(),
    }
    count = count + 1;
    if count == 3 { break; } // the fourth entry is unmapped, and will lead to an error.
}

If we try to set contig name on an unmapped read, we will get an error

let mut reader = nanalogue_bam_reader(&"examples/example_1.bam")?;
let mut count = 0;
for record in reader.records(){
    if count < 3 {
        count = count + 1;
        continue;
    }
    let r = record?;
    let mut curr_read = CurrRead::default().set_read_state_and_id(&r)?;
    curr_read.set_contig_name(&r)?;
}

If we call the method twice, we should hit a panic

let mut reader = nanalogue_bam_reader(&"examples/example_1.bam")?;
for record in reader.records(){
    let r = record?;
    let mut curr_read = CurrRead::default().set_read_state_and_id(&r)?
        .set_contig_name(&r)?.set_contig_name(&r)?;
    break;
}

pub fn contig_name(&self) -> Result<&str, Error>

gets contig name

Errors

If instance is unmapped or contig name has not been set

pub fn read_id(&self) -> &str

gets read id

pub fn strand(&self) -> char

Returns the character corresponding to the strand

use nanalogue_core::{CurrRead, Error, nanalogue_bam_reader};
use rust_htslib::bam::Read;
let mut reader = nanalogue_bam_reader(&"examples/example_1.bam")?;
let mut count = 0;
for record in reader.records(){
    let r = record?;
    let mut curr_read = CurrRead::default().set_read_state_and_id(&r)?;
    let strand = curr_read.strand() else { unreachable!() };
    match (count, strand) {
        (0, '+') | (1, '+') | (2, '-') | (3, '.') => {},
        _ => unreachable!(),
    }
    count = count + 1;
}

pub fn seq_on_ref_coords( &self, record: &Record, region: &Bed3<i32, u64>, ) -> Result<Vec<u8>, Error>

Returns read sequence overlapping with a genomic region

Errors

If getting sequence coordinates from reference coordinates fails, see CurrRead::seq_and_qual_on_ref_coords

Example

use bedrs::Bed3;
use nanalogue_core::{CurrRead, Error, nanalogue_bam_reader};
use rust_htslib::bam::Read;

let mut reader = nanalogue_bam_reader(&"examples/example_1.bam")?;
let mut count = 0;
for record in reader.records() {
    let r = record?;
    let curr_read = CurrRead::default().try_from_only_alignment(&r)?;

    // Skip unmapped reads
    if !curr_read.read_state().is_unmapped() {
        let (contig_id, start) = curr_read.contig_id_and_start()?;
        let align_len = curr_read.align_len()?;

        // Create a region that overlaps with the read but is short of one bp.
        // Note that this BAM file has reads with all bases matching perfectly
        // with the reference.
        let region = Bed3::new(contig_id, start, start + align_len - 1);
        let seq_subset = curr_read.seq_on_ref_coords(&r, &region)?;

        // Check for sequence length match
        assert_eq!(curr_read.seq_len()? - 1, u64::try_from(seq_subset.len())?);

        // Create a region with no overlap at all and check we get no data
        let region = Bed3::new(contig_id, start + align_len, start + align_len + 2);
        match curr_read.seq_on_ref_coords(&r, &region){
            Err(Error::UnavailableData(_)) => (),
            _ => unreachable!(),
        };

    }
    count += 1;
}

pub fn seq_and_qual_on_ref_coords( &self, record: &Record, region: &Bed3<i32, u64>, ) -> Result<Vec<Option<(bool, u8, u8)>>, Error>

Returns match-or-mismatch, read sequence, base-quality values overlapping with genomic region.

Returns a vector of Options:

• is a None at a deletion i.e. a base on the reference and not on the read.
• is a Some(bool, u8, u8) at a match/mismatch/insertion. The first u8 is the base itself and the bool is true if a match/mismatch and false if an insertion, and the second u8 is the base quality (0-93), which the BAM format sets to 255 if no quality is available for the entire read.

Because sequences are encoded using 4-bit values into a [u8], we need to use rust-htslib functions to convert them into 8-bit values and then use the Index<usize> trait on sequences from rust-htslib.

Errors

If the read does not intersect with the specified region, see CurrRead::seq_coords_from_ref_coords

Example

Example 1

use bedrs::Bed3;
use nanalogue_core::{CurrRead, Error, nanalogue_bam_reader};
use rust_htslib::bam::Read;

let mut reader = nanalogue_bam_reader(&"examples/example_5_valid_basequal.sam")?;
for record in reader.records() {
    let r = record?;
    let curr_read = CurrRead::default().try_from_only_alignment(&r)?;

    let region = Bed3::new(0, 0, 20);
    let seq_subset = curr_read.seq_and_qual_on_ref_coords(&r, &region)?;
    assert_eq!(seq_subset, [Some((true, b'T', 32)), Some((true, b'C', 0)),
        Some((true, b'G', 69)), Some((true, b'T', 80)), Some((true, b'T', 79)),
        Some((true, b'T', 81)), Some((true, b'C', 29)), Some((true, b'T', 30))]);

    // Create a region with no overlap at all and check we get no data
    let region = Bed3::new(0, 20, 22);
    match curr_read.seq_and_qual_on_ref_coords(&r, &region){
        Err(Error::UnavailableData(_)) => (),
        _ => unreachable!(),
    };

}

Example 2

use bedrs::Bed3;
use nanalogue_core::{CurrRead, Error, nanalogue_bam_reader};
use rust_htslib::bam::Read;

let mut reader = nanalogue_bam_reader(&"examples/example_7.sam")?;
for record in reader.records() {
    let r = record?;
    let curr_read = CurrRead::default().try_from_only_alignment(&r)?;

    let region = Bed3::new(0, 0, 20);
    let seq_subset = curr_read.seq_and_qual_on_ref_coords(&r, &region)?;
    assert_eq!(seq_subset, [Some((true, b'T', 32)), None, None,
        Some((false, b'A', 0)), Some((true, b'T', 0)), Some((true, b'T', 79)),
        Some((true, b'T', 81)), Some((true, b'G', 29)), Some((true, b'T', 30))]);

}

pub fn seq_coords_from_ref_coords( &self, record: &Record, region: &Bed3<i32, u64>, ) -> Result<Vec<Option<(bool, usize)>>, Error>

Extract sequence coordinates corresponding to a region on the reference genome.

The vector we return contains Some((bool, usize)) entries where both the reference and the read have bases, and None where bases from the reference are missing on the read:

• matches or mismatches, we record the coordinate. so SNPs for example (i.e. a 1 bp difference from the ref) will show up as Some((true, _)).
• a deletion or a ref skip (“N” in cigar) will show up as None.
• insertions are preserved i.e. bases in the middle of a read present on the read but not on the reference are Some((false, _))
• clipped bases at the end of the read are not preserved. These are bases on the read but not on the reference and are denoted as soft or hard clips on the CIGAR string e.g. barcodes from sequencing
• some CIGAR combinations are illogical and we are assuming they do not happen e.g. a read’s CIGAR can end with, say 10D20S, this means last 10 bp are in a deletion and the next 20 are a softclip. This is illogical, as they must be combined into a 30-bp softclip i.e. 30S. So if the aligner produces such illogical states, then the sequence coordinates reported here may be erroneous.

If the read does not intersect with the region, we return an Error (see below). If the read does intersect with the region but we cannot retrieve any bases, we return an empty vector (I am not sure if we will run into this scenario).

Errors

If the read does not intersect with the specified region, or if usize conversions fail.

Example

use bedrs::Bed3;
use nanalogue_core::{CurrRead, Error, nanalogue_bam_reader};
use rust_htslib::bam::Read;

let mut reader = nanalogue_bam_reader(&"examples/example_7.sam")?;
for record in reader.records() {
    let r = record?;
    let curr_read = CurrRead::default().try_from_only_alignment(&r)?;

    let region = Bed3::new(0, 9, 13);
    let seq_subset = curr_read.seq_coords_from_ref_coords(&r, &region)?;
    // there are deletions on the read above
    assert_eq!(seq_subset, vec![Some((true, 0)), None, None, Some((false, 1)), Some((true, 2))]);

    // Create a region with no overlap at all and check we get no data
    let region = Bed3::new(0, 20, 22);
    match curr_read.seq_coords_from_ref_coords(&r, &region){
        Err(Error::UnavailableData(_)) => (),
        _ => unreachable!(),
    };

}

pub fn set_mod_data( self, record: &Record, mod_thres: ThresholdState, min_qual: u8, ) -> Result<CurrRead<AlignAndModData>, Error>

sets modification data using the BAM record

Errors

If tags in the BAM record containing the modification information (MM, ML) contain mistakes.

pub fn set_mod_data_restricted<G, H>( self, record: &Record, mod_thres: ThresholdState, mod_fwd_pos_filter: G, mod_filter_base_strand_tag: H, min_qual: u8, ) -> Result<CurrRead<AlignAndModData>, Error>

where G: Fn(&usize) -> bool, H: Fn(&u8, &char, &ModChar) -> bool,

sets modification data using BAM record but restricted to the specified filters

Errors

If tags in the BAM record containing the modification information (MM, ML) contain mistakes.

impl CurrRead<OnlyAlignDataComplete>

pub fn set_mod_data_restricted_options<S: InputModOptions + InputRegionOptions>( self, record: &Record, mod_options: &S, ) -> Result<CurrRead<AlignAndModData>, Error>

sets modification data using BAM record but with restrictions applied by the crate::InputMods options for example.

Errors

If a region filter is specified and we fail to convert current instance to Bed, and if parsing the MM/ML BAM tags fails (presumably because they are malformed).

impl CurrRead<AlignAndModData>

pub fn mod_data(&self) -> &(BaseMods, ThresholdState)

gets modification data

pub fn windowed_mod_data_restricted<F>( &self, window_function: &F, win_options: InputWindowing, tag: ModChar, ) -> Result<Vec<F32Bw0and1>, Error>

where F: Fn(&[u8]) -> Result<F32Bw0and1, Error>,

window modification data with restrictions. If a read has the same modification on both the basecalled strand and its complement, then windows along both are returned.

If win_size exceeds the modification data length, no windows are produced.

Errors

Returns an error if the window function returns an error.

pub fn base_count_per_mod(&self) -> HashMap<ModChar, u32>

Performs a count of number of bases per modified type. Note that this result depends on the type of filtering done while the struct was created e.g. by modification threshold.

Panics

Panics if the number of modifications exceeds u32::MAX (approximately 4.2 billion).

use nanalogue_core::{CurrRead, Error, ModChar, nanalogue_bam_reader, ThresholdState};
use rust_htslib::bam::Read;
use std::collections::HashMap;
let mut reader = nanalogue_bam_reader(&"examples/example_1.bam")?;
let mut count = 0;
for record in reader.records(){
    let r = record?;
    let curr_read = CurrRead::default().set_read_state_and_id(&r)?
        .set_mod_data(&r, ThresholdState::GtEq(180), 0)?;
    let mod_count = curr_read.base_count_per_mod();
    let zero_count = HashMap::from([(ModChar::new('T'), 0)]);
    let a = HashMap::from([(ModChar::new('T'), 3)]);
    let b = HashMap::from([(ModChar::new('T'), 1)]);
    let c = HashMap::from([(ModChar::new('T'), 3),(ModChar::new('ᰠ'), 0)]);
    match (count, mod_count) {
        (0, v) => assert_eq!(v, zero_count),
        (1, v) => assert_eq!(v, a),
        (2, v) => assert_eq!(v, b),
        (3, v) => assert_eq!(v, c),
        _ => unreachable!(),
    }
    count = count + 1;
}

Trait Implementations

impl<S: Clone + CurrReadState> Clone for CurrRead<S>

fn clone(&self) -> CurrRead<S>

Returns a duplicate of the value.

fn clone_from(&mut self, source: &Self)

Performs copy-assignment from source.

impl<S: Debug + CurrReadState> Debug for CurrRead<S>

fn fmt(&self, f: &mut Formatter<'_>) -> Result

Formats the value using the given formatter.

impl Default for CurrRead<NoData>

Implements defaults for CurrRead

fn default() -> Self

Returns the “default value” for a type.

impl<'de> Deserialize<'de> for CurrRead<AlignAndModData>

fn deserialize<D>(deserializer: D) -> Result<Self, D::Error>

where D: Deserializer<'de>,

Deserialize this value from the given Serde deserializer.

impl<S> Display for CurrRead<S>

where S: CurrReadState, CurrRead<S>: DisplayCondensedModData,

fn fmt(&self, f: &mut Formatter<'_>) -> Result

Formats the value using the given formatter.

impl FilterByRefCoords for CurrRead<AlignAndModData>

Implements filter by reference coordinates for our CurrRead. This only filters modification data.

fn filter_by_ref_pos(&mut self, start: i64, end: i64) -> Result<(), Error>

filters modification data by reference position i.e. all pos such that start <= pos < end are retained. does not use contig in filtering.

impl<S: PartialEq + CurrReadState> PartialEq for CurrRead<S>

fn eq(&self, other: &CurrRead<S>) -> bool

Tests for self and other values to be equal, and is used by ==.

fn ne(&self, other: &Rhs) -> bool

Tests for !=. The default implementation is almost always sufficient, and should not be overridden without very good reason.

impl Serialize for CurrRead<AlignAndModData>

fn serialize<S>(&self, serializer: S) -> Result<S::Ok, S::Error>

where S: Serializer,

Serialize this value into the given Serde serializer.

impl<S: CurrReadStateWithAlign + CurrReadState> TryFrom<&CurrRead<S>> for StrandedBed3<i32, u64>

Converts CurrRead to StrandedBed3

use bedrs::{Coordinates, Strand};
use bedrs::prelude::StrandedBed3;
use nanalogue_core::{CurrRead, Error, nanalogue_bam_reader};
use rust_htslib::bam::Read;
let mut reader = nanalogue_bam_reader(&"examples/example_1.bam")?;
let mut count = 0;
for record in reader.records(){
    let r = record?;
    let mut curr_read = CurrRead::default().try_from_only_alignment(&r)?;
    let Ok(bed3_stranded) = StrandedBed3::try_from(&curr_read) else {unreachable!()};
    let exp_bed3_stranded = match count {
        0 => StrandedBed3::new(0, 9, 17, Strand::Forward),
        1 => StrandedBed3::new(2, 23, 71, Strand::Forward),
        2 => StrandedBed3::new(1, 3, 36, Strand::Reverse),
        3 => StrandedBed3::empty(),
        _ => unreachable!(),
    };
    assert_eq!(*bed3_stranded.chr(), *exp_bed3_stranded.chr());
    assert_eq!(bed3_stranded.start(), exp_bed3_stranded.start());
    assert_eq!(bed3_stranded.end(), exp_bed3_stranded.end());
    assert_eq!(bed3_stranded.strand(), exp_bed3_stranded.strand());
    count = count + 1;
}

type Error = Error

The type returned in the event of a conversion error.

fn try_from(value: &CurrRead<S>) -> Result<Self, Self::Error>

Performs the conversion.

impl TryFrom<CurrRead<AlignAndModData>> for CurrReadBuilder

type Error = Error

The type returned in the event of a conversion error.

fn try_from(curr_read: CurrRead<AlignAndModData>) -> Result<Self, Self::Error>

Performs the conversion.

impl TryFrom<CurrReadBuilder> for CurrRead<AlignAndModData>

type Error = Error

The type returned in the event of a conversion error.

fn try_from(serialized: CurrReadBuilder) -> Result<Self, Self::Error>

Performs the conversion.

impl TryFrom<Rc<Record>> for CurrRead<AlignAndModData>

Convert a rust htslib rc record into our struct. I think the rc datatype is just like the normal record, except the record datatype is not destroyed and created every time a new record is read (or something like that). All comments I’ve made for the TryFrom<Record> function apply here as well.

type Error = Error

The type returned in the event of a conversion error.

fn try_from(record: Rc<Record>) -> Result<Self, Self::Error>

Performs the conversion.

impl TryFrom<Record> for CurrRead<AlignAndModData>

Convert a rust htslib record to our CurrRead struct. NOTE: This operation loads many types of data from the record and you may not be interested in all of them. So, unless you know for sure that you are dealing with a small number of reads, please do not use this function, and call only a subset of the individual invocations below in your program for the sake of speed and/or memory.

type Error = Error

The type returned in the event of a conversion error.

fn try_from(record: Record) -> Result<Self, Self::Error>

Performs the conversion.

impl<S: CurrReadState> StructuralPartialEq for CurrRead<S>