diff --git a/src/chainer.rs b/src/chainer.rs index 101aded3..a9bc719c 100644 --- a/src/chainer.rs +++ b/src/chainer.rs @@ -13,7 +13,6 @@ const N_PRECOMPUTED: usize = 1024; #[derive(Debug, Ord, PartialOrd, Eq, PartialEq, Clone, Copy)] pub struct Anchor { - pub ref_id: usize, pub ref_start: usize, pub query_start: usize, } @@ -102,9 +101,10 @@ impl Chainer { for j in (lookup_end..i).rev() { let aj = &anchors[j]; - if ai.ref_id != aj.ref_id { - break; - } + // TODO + // if ai.ref_id != aj.ref_id { + // break; + // } let Some(dq) = ai.query_start.checked_sub(aj.query_start) else { // Not collinear @@ -143,7 +143,8 @@ impl Chainer { // The above runtime heuristic sometimes skips the anchor that // represents the so-far optimal chain and can then give suboptimal // results. To mitigate the issue, we explicitly check that anchor. - if best_index != usize::MAX && ai.ref_id == anchors[best_index].ref_id { + if best_index != usize::MAX { + // TODO && ai.ref_id == anchors[best_index].ref_id { let aj = &anchors[best_index]; if ai.query_start > aj.query_start { let dq = ai.query_start - aj.query_start; @@ -225,7 +226,7 @@ impl Chainer { n_anchors += anchors.len(); let chaining_timer = Instant::now(); trace!("Chaining {} anchors", anchors.len()); - anchors.sort_unstable_by_key(|a| (a.ref_id, a.ref_start, a.query_start)); + anchors.sort_unstable_by_key(|a| (a.ref_start, a.query_start)); anchors.dedup(); // trace!( @@ -315,18 +316,15 @@ fn add_to_anchors_full( if entry.hash() != forward_hash { break; } - let ref_start = entry.position(); + let ref_start = entry.ref_start(); let ref_end = ref_start + entry.strobe2_offset() + index.k(); let length_diff = (query_end - query_start).abs_diff(ref_end - ref_start); if length_diff <= min_length_diff { - let ref_id = entry.reference_index(); anchors.push(Anchor { - ref_id, ref_start, query_start, }); anchors.push(Anchor { - ref_id, ref_start: ref_end - index.k(), query_start: query_end - index.k(), }); @@ -348,11 +346,9 @@ fn add_to_anchors_partial( if entry.get_hash_partial_forward() != forward_hash { break; } - let ref_id = entry.reference_index(); let ref_start = entry.strobe_extent_partial().0; anchors.push(Anchor { - ref_id, ref_start, query_start, }); @@ -439,7 +435,7 @@ impl ChainingResult { ref_start: first.ref_start, ref_end: last.ref_start + k, matching_bases, - ref_id: last.ref_id, + // TODO ref_id: last.ref_id, score: score + chain_anchors.len() as f32 * self.parameters.matches_weight, is_revcomp, anchors: chain_anchors, @@ -457,10 +453,10 @@ mod test { let chainer = Chainer::new(20, ChainingParameters::default()); #[rustfmt::skip] let anchors = vec![ - Anchor { ref_id: 0, ref_start: 0, query_start: 0, }, - Anchor { ref_id: 0, ref_start: 30, query_start: 20, }, - Anchor { ref_id: 0, ref_start: 60, query_start: 0, }, - Anchor { ref_id: 0, ref_start: 95, query_start: 35, }, + Anchor { ref_start: 0, query_start: 0, }, + Anchor { ref_start: 30, query_start: 20, }, + Anchor { ref_start: 60, query_start: 0, }, + Anchor { ref_start: 95, query_start: 35, }, ]; let chaining_result = chainer.collinear_chaining(anchors, 2000); @@ -477,9 +473,9 @@ mod test { let chainer = Chainer::new(20, ChainingParameters::default()); #[rustfmt::skip] let anchors = [ - Anchor { ref_id: 0, ref_start: 0, query_start: 0, }, - Anchor { ref_id: 0, ref_start: 20, query_start: 20, }, - Anchor { ref_id: 0, ref_start: 40, query_start: 40, }, + Anchor { ref_start: 0, query_start: 0, }, + Anchor { ref_start: 20, query_start: 20, }, + Anchor { ref_start: 40, query_start: 40, }, ]; let chaining_result = chainer.collinear_chaining(anchors[0..1].to_vec(), 200); let score1 = chaining_result.best_score; @@ -502,8 +498,8 @@ mod test { let chainer = Chainer::new(20, ChainingParameters::default()); #[rustfmt::skip] let anchors = vec![ - Anchor { ref_id: 0, ref_start: 0, query_start: 0, }, - Anchor { ref_id: 0, ref_start: 11, query_start: 1, }, + Anchor { ref_start: 0, query_start: 0, }, + Anchor { ref_start: 11, query_start: 1, }, ]; let chaining_result = chainer.collinear_chaining(anchors, 2000); // dr=11, dq=1 gives a diagonal ratio of 11.0, exceeding the default max of 10. diff --git a/src/index.rs b/src/index.rs index 7fa40014..327c6044 100644 --- a/src/index.rs +++ b/src/index.rs @@ -25,8 +25,7 @@ pub struct RefRandstrobe { /// [`REF_RANDSTROBE_HASH_MASK`] is the mask for strobe1 and strobe2 hashes (including their orientations) /// The [`STROBE2_OFFSET_BITS`] constant specifies the number of bits reserved for the offset hash_offset: u64, - position: u32, - ref_index: u32, + ref_start: u64, } /// Mask for the part of the randstrobe hash that includes individual strobe hashes and orientations @@ -35,15 +34,13 @@ pub const REF_RANDSTROBE_HASH_MASK: u64 = 0xFFFFFFFFFFFFFF00; pub const STROBE2_OFFSET_BITS: u32 = 8; /// Mask for the part of the randstrobe hash that includes the offset between first and second strobe pub const STROBE2_OFFSET_MASK: u64 = (1u64 << STROBE2_OFFSET_BITS) - 1; -pub const REF_RANDSTROBE_MAX_NUMBER_OF_REFERENCES: usize = u32::MAX as usize; impl RefRandstrobe { - pub fn new(hash: RandstrobeHash, ref_index: u32, position: u32, offset: u8) -> Self { + pub fn new(hash: RandstrobeHash, ref_start: usize, offset: u8) -> Self { let hash_offset = (hash & REF_RANDSTROBE_HASH_MASK) | (offset as u64); RefRandstrobe { hash_offset, - position, - ref_index, + ref_start: ref_start as u64, } } @@ -51,12 +48,8 @@ impl RefRandstrobe { self.hash_offset & REF_RANDSTROBE_HASH_MASK } - pub fn position(&self) -> usize { - self.position as usize - } - - pub fn reference_index(&self) -> usize { - self.ref_index as usize + pub fn ref_start(&self) -> usize { + self.ref_start as usize } pub fn strobe2_offset(&self) -> usize { @@ -88,6 +81,8 @@ pub struct StrobemerIndex { pub struct IndexEntry<'a> { strobemer_index: &'a StrobemerIndex, + + /// An index that points to an item in the randstrobes vector pub position: usize, } @@ -218,18 +213,14 @@ impl<'a> IndexEntry<'a> { self.randstrobe().hash() } - pub fn position(&self) -> usize { - self.randstrobe().position() + pub fn ref_start(&self) -> usize { + self.randstrobe().ref_start() } pub fn strobe2_offset(&self) -> usize { self.randstrobe().strobe2_offset() } - pub fn reference_index(&self) -> usize { - self.randstrobe().reference_index() - } - pub fn get_hash_partial(&self) -> RandstrobeHash { self.strobemer_index.randstrobes[self.position].hash() & self.strobemer_index.parameters.randstrobe.main_hash_mask @@ -245,9 +236,9 @@ impl<'a> IndexEntry<'a> { } pub fn strobe_extent_partial(&self) -> (usize, usize) { - let p = self.strobemer_index.randstrobes[self.position].position; + let p = self.strobemer_index.randstrobes[self.position].ref_start(); - (p as usize, p as usize + self.k()) + (p, p + self.k()) } /// Count number of hits for the randstrobe *and* its "reverse complement" @@ -550,14 +541,12 @@ mod tests { #[test] fn ref_randstrobe() { let hash: u64 = 0x1234567890ABCDEFu64 & REF_RANDSTROBE_HASH_MASK; - let ref_index: u32 = (REF_RANDSTROBE_MAX_NUMBER_OF_REFERENCES - 1) as u32; + let ref_start = u64::MAX as usize; let offset = 255; - let position = !0; - let rr = RefRandstrobe::new(hash, ref_index, position, offset); + let rr = RefRandstrobe::new(hash, ref_start, offset); assert_eq!(rr.hash(), hash); - assert_eq!(rr.position(), position as usize); - assert_eq!(rr.reference_index(), ref_index as usize); + assert_eq!(rr.ref_start(), ref_start); assert_eq!(rr.strobe2_offset(), offset as usize); } @@ -591,9 +580,11 @@ mod tests { let refseq = read_ref("tests/phix.fasta").unwrap(); let seq_decoded = refseq.contig(0).decode_all(); let rc_seq = reverse_complement(&seq_decoded); - let mut rc_refseq = RefSequence::new(); - rc_refseq.push("phix_rc".to_string(), PackedSeq::from_slice(&rc_seq)); - + let rc_refseq = RefSequence::new( + PackedSeq::from_slice(&rc_seq), + vec![0], + vec!["phix_rc".to_string()], + ); let parameters = SeedingParameters::new(300); let bits = parameters.syncmer.pick_bits(&refseq); diff --git a/src/indexer.rs b/src/indexer.rs index 70dc291e..24161453 100644 --- a/src/indexer.rs +++ b/src/indexer.rs @@ -1,5 +1,5 @@ use crate::index::{BucketIndex, RandstrobeHash, RefRandstrobe, StrobemerIndex}; -use crate::packed_seq::PackedSeq; +use crate::packed_seq::PackedSeqSlice; use crate::refseq::RefSequence; use crate::seeding::{RandstrobeIterator, SeedingParameters, SyncmerIterator, SyncmerParameters}; @@ -195,7 +195,13 @@ fn make_randstrobes_parallel( if j >= refseq.names.len() { break; } - assign_randstrobes(refseq, parameters, j, *slices[j].lock().unwrap()); + let start = refseq.contig_start(j); + assign_randstrobes( + &refseq.contig(j), + parameters, + start, + *slices[j].lock().unwrap(), + ); } }); } @@ -206,12 +212,11 @@ fn make_randstrobes_parallel( /// Compute randstrobes of one reference contig and assign them to the provided slice fn assign_randstrobes( - refseq: &RefSequence, + seq: &PackedSeqSlice, parameters: &SeedingParameters, - ref_index: usize, + start: usize, randstrobes: &mut [RefRandstrobe], ) { - let seq = refseq.contig(ref_index); if seq.len() < parameters.randstrobe.w_max { return; } @@ -228,12 +233,8 @@ fn assign_randstrobes( for (i, randstrobe) in randstrobe_iter.enumerate() { n += 1; let offset = randstrobe.strobe2_pos - randstrobe.strobe1_pos; - randstrobes[i] = RefRandstrobe::new( - randstrobe.hash, - ref_index as u32, - randstrobe.strobe1_pos as u32, - offset as u8, - ); + let strobe1_start = randstrobe.strobe1_pos + start; + randstrobes[i] = RefRandstrobe::new(randstrobe.hash, strobe1_start, offset as u8); } debug_assert_eq!(n, randstrobes.len()); } @@ -254,7 +255,7 @@ fn count_all_randstrobes( if j >= refseq.names.len() { break; } - let count = count_randstrobes(refseq.contig(j), parameters); + let count = count_randstrobes(&refseq.contig(j), parameters); mutex.lock().unwrap()[j] = count; } }); @@ -265,7 +266,7 @@ fn count_all_randstrobes( } /// Count randstrobes by counting syncmers (operates directly on PackedSeq) -fn count_randstrobes(seq: &PackedSeq, parameters: &SeedingParameters) -> usize { +fn count_randstrobes(seq: &PackedSeqSlice, parameters: &SeedingParameters) -> usize { SyncmerIterator::new( seq, parameters.syncmer.k, @@ -368,8 +369,7 @@ mod test { #[test] fn index_empty_reference() { - let mut refseq = RefSequence::new(); - refseq.push("name".to_string(), PackedSeq::new()); + let refseq = RefSequence::new(PackedSeq::new(), vec![0], vec!["name".to_string()]); let parameters = SeedingParameters::new(150); let bits = parameters.syncmer.pick_bits(&refseq); let (_index2, stats) = make_index(&refseq, parameters, bits, 0.1, 1); @@ -381,6 +381,6 @@ mod test { let refseq = read_ref("tests/phix.fasta").unwrap(); let parameters = SeedingParameters::new(150); - assert_eq!(count_randstrobes(refseq.contig(0), ¶meters), 1090); + assert_eq!(count_randstrobes(&refseq.contig(0), ¶meters), 1090); } } diff --git a/src/io/fasta.rs b/src/io/fasta.rs index 90e42959..657cdb2a 100644 --- a/src/io/fasta.rs +++ b/src/io/fasta.rs @@ -26,8 +26,8 @@ fn is_valid_name(name: &[u8]) -> bool { true } -fn check_duplicate_names(refseq: &RefSequence) -> Result<(), SequenceIOError> { - let mut names = refseq.names.clone(); +fn check_duplicate_names(names: &[String]) -> Result<(), SequenceIOError> { + let mut names: Vec<_> = names.iter().collect(); names.sort(); for window in names.windows(2) { if window[0] == window[1] { @@ -109,9 +109,9 @@ impl Iterator for FastaReader { } pub fn read_fasta(reader: &mut R) -> Result { - let mut refseq = RefSequence::new(); - let mut current_name: Option = None; - let mut current_seq = PackedSeq::new(); + let mut names = vec![]; + let mut starts = vec![]; + let mut seq = PackedSeq::new(); let mut line = String::new(); loop { @@ -120,17 +120,14 @@ pub fn read_fasta(reader: &mut R) -> Result break, Ok(_) => { if let Some(header) = line.strip_prefix('>') { - if let Some(name) = current_name.take() { - refseq.push(name, current_seq); - current_seq = PackedSeq::new(); - } - let (name, _comment) = split_header(header); - current_name = Some(name.to_string()); - } else if current_name.is_some() { - for c in line.trim_ascii_end().bytes() { - current_seq.push(c); - } - } else if !line.trim_ascii().is_empty() { + starts.push(seq.len()); + names.push(split_header(header).0); + } else if !names.is_empty() { + let line = line.trim_ascii_end(); + seq.extend(line.bytes()); + } else if line.trim_ascii().is_empty() { + // ignore empty lines + } else { return Err(SequenceIOError::Fasta( "FASTA file must start with '>'".to_string(), )); @@ -139,18 +136,16 @@ pub fn read_fasta(reader: &mut R) -> Result return Err(SequenceIOError::IO(e)), } } - if let Some(name) = current_name { - refseq.push(name, current_seq); - } - for name in &refseq.names { + for name in &names { if !is_valid_name(name.as_bytes()) { return Err(SequenceIOError::Name); } } - check_duplicate_names(&refseq)?; + check_duplicate_names(&names)?; + assert_eq!(names.len(), starts.len()); - Ok(refseq) + Ok(RefSequence::new(seq, starts, names)) } pub fn read_ref>(path: P) -> Result { @@ -171,7 +166,7 @@ mod tests { let refseq = read_ref("tests/phix.fasta").unwrap(); assert_eq!(refseq.names.len(), 1); assert_eq!(refseq.names[0], "NC_001422.1"); - assert_eq!(refseq.contig_len(0), 5386); + assert_eq!(refseq.contig(0).len(), 5386); assert_eq!(refseq.contig(0).decode(0, 5), b"GAGTT"); } diff --git a/src/io/paf.rs b/src/io/paf.rs index b1ae53c4..64038199 100644 --- a/src/io/paf.rs +++ b/src/io/paf.rs @@ -1,6 +1,6 @@ use std::fmt::Display; -use crate::io::record::End; +use crate::{io::record::End, refseq::ContigPosition}; /* PAF columns (see https://github.com/lh3/miniasm/blob/master/PAF.md): * 1 query name @@ -27,8 +27,8 @@ pub struct PafRecord { pub is_revcomp: bool, pub target_name: String, pub target_length: u64, - pub target_start: u64, - pub target_end: u64, + pub target_start: ContigPosition, + pub target_end: ContigPosition, pub matching_bases: u64, pub alignment_length: u64, pub mapping_quality: Option, @@ -51,8 +51,8 @@ impl Display for PafRecord { if self.is_revcomp { "-" } else { "+" }, self.target_name, self.target_length, - self.target_start, - self.target_end, + self.target_start.0, + self.target_end.0, self.matching_bases, self.alignment_length, self.mapping_quality.unwrap_or(255), diff --git a/src/io/sam.rs b/src/io/sam.rs index ba078260..3c3917a5 100644 --- a/src/io/sam.rs +++ b/src/io/sam.rs @@ -2,7 +2,7 @@ use std::fmt::{Display, Formatter}; use crate::cigar::Cigar; use crate::details::Details; -use crate::refseq::RefSequence; +use crate::refseq::{ContigPosition, RefSequence}; pub const PAIRED: u16 = 1; pub const PROPER_PAIR: u16 = 2; @@ -23,12 +23,11 @@ pub struct SamRecord { pub query_name: String, pub flags: u16, pub reference_name: Option, - /// 0-based position, converted to 1-based for output - pub pos: Option, + pub pos: Option, pub mapq: u8, pub cigar: Option, pub mate_reference_name: Option, - pub mate_pos: Option, + pub mate_pos: Option, pub template_len: Option, pub query_sequence: Option>, pub query_qualities: Option>, @@ -68,7 +67,7 @@ impl Display for SamRecord { // 12 optional fields let mate_pos = match self.mate_pos { - Some(pos) => pos + 1, + Some(pos) => pos.0 + 1, None => 0, }; let query_sequence = match &self.query_sequence { @@ -84,7 +83,7 @@ impl Display for SamRecord { _ => "*", }; let pos = match self.pos { - Some(pos) => pos + 1, + Some(pos) => pos.0 + 1, None => 0, }; let cigar = match &self.cigar { @@ -183,7 +182,7 @@ impl Display for SamHeader<'_> { f, "@SQ\tSN:{}\tLN:{}", self.refseq.names[i], - self.refseq.contig_len(i) + self.refseq.contig(i).len() )?; } if let Some(read_group) = &self.read_group { diff --git a/src/main.rs b/src/main.rs index 51fccbaf..692cc553 100644 --- a/src/main.rs +++ b/src/main.rs @@ -21,9 +21,7 @@ use thiserror::Error; use strobealign::aligner::{Aligner, Scores}; use strobealign::chainer::{Chainer, ChainingParameters}; use strobealign::details::Details; -use strobealign::index::{ - IndexReadingError, REF_RANDSTROBE_MAX_NUMBER_OF_REFERENCES, StrobemerIndex, read_index, -}; +use strobealign::index::{IndexReadingError, StrobemerIndex, read_index}; use strobealign::insertsize::InsertSizeDistribution; use strobealign::io::SequenceIOError; use strobealign::io::fasta; @@ -437,13 +435,6 @@ fn run() -> Result<(), CliError> { if refseq.names.len() != 1 { "s" } else { "" }, max_contig_size as f64 / 1E6 ); - if refseq.names.len() > REF_RANDSTROBE_MAX_NUMBER_OF_REFERENCES { - error!( - "Too many reference sequences. Current maximum is {}.", - REF_RANDSTROBE_MAX_NUMBER_OF_REFERENCES - ); - exit(1); - } // Create the index debug!("Auxiliary hash length: {}", args.aux_len); @@ -959,6 +950,7 @@ impl Mapper<'_> { &r1, &r2, self.index, + self.refseq, &mut self.abundances, self.mapping_parameters.rescue_distance, &mut isizedist, @@ -969,6 +961,7 @@ impl Mapper<'_> { } else { abundances_single_end_read( &r1, + self.refseq, self.index, &mut self.abundances, self.mapping_parameters.rescue_distance, @@ -991,7 +984,7 @@ pub fn output_abundances( ) -> io::Result<()> { #[allow(clippy::needless_range_loop)] for i in 0..refseq.names.len() { - let normalized = abundances[i] / refseq.contig_len(i) as f64; + let normalized = abundances[i] / refseq.contig(i).len() as f64; writeln!(out, "{}\t{:.6}", refseq.names[i], normalized)?; } Ok(()) diff --git a/src/maponly.rs b/src/maponly.rs index 2c1530c4..8e3c0991 100644 --- a/src/maponly.rs +++ b/src/maponly.rs @@ -38,14 +38,16 @@ pub fn map_single_end_read( } else { let mapq = mapping_quality(&nams); ( - vec![paf_record_from_nam( + paf_record_from_nam( &nams[0], &record.name, refseq, record.sequence.len(), Some(mapq), End::None, - )], + ) + .into_iter() + .collect::>(), nam_details.into(), ) } @@ -56,6 +58,7 @@ pub fn map_single_end_read( /// This implements abundance estimation mode (`--aemb`) pub fn abundances_single_end_read( record: &SequenceRecord, + refseq: &RefSequence, index: &StrobemerIndex, abundances: &mut [f64], rescue_distance: usize, @@ -78,7 +81,8 @@ pub fn abundances_single_end_read( .count(); let weight = record.sequence.len() as f64 / n_best as f64; for nam in &nams[0..n_best] { - abundances[nam.ref_id] += weight; + let reference_id = refseq.unflatten(nam.ref_start).0; + abundances[reference_id] += weight; } } @@ -90,22 +94,28 @@ fn paf_record_from_nam( query_length: usize, mapq: Option, end: End, -) -> PafRecord { - PafRecord { +) -> Option { + let (reference_id, target_start) = refseq.unflatten(nam.ref_start); + let (end_reference_id, target_end) = refseq.unflatten(nam.ref_end); + if reference_id != end_reference_id { + return None; + } + + Some(PafRecord { query_name: name.into(), end, query_length: query_length as u64, query_start: nam.query_start as u64, query_end: nam.query_end as u64, is_revcomp: nam.is_revcomp, - target_name: refseq.names[nam.ref_id].clone(), - target_length: refseq.contig_len(nam.ref_id) as u64, - target_start: nam.ref_start as u64, - target_end: nam.ref_end as u64, + target_name: refseq.names[reference_id].clone(), + target_length: refseq.contig(reference_id).len() as u64, + target_start, + target_end, matching_bases: nam.matching_bases as u64, alignment_length: (nam.ref_end - nam.ref_start) as u64, mapping_quality: mapq, - } + }) } /// Map a paired-end read pair to the reference and return PAF records @@ -148,7 +158,7 @@ pub fn map_paired_end_read( match get_best_paired_mapping_location(&nam_pairs, &nams1, &nams2, insert_size_distribution) { MappedNams::Individual(best1, best2) => { if let Some(nam1) = best1 { - records.push(paf_record_from_nam( + records.extend(paf_record_from_nam( nam1, &r1.name, refseq, @@ -158,7 +168,7 @@ pub fn map_paired_end_read( )); } if let Some(nam2) = best2 { - records.push(paf_record_from_nam( + records.extend(paf_record_from_nam( nam2, &r2.name, refseq, @@ -169,7 +179,7 @@ pub fn map_paired_end_read( } } MappedNams::Pair(nam1, nam2, _) => { - records.push(paf_record_from_nam( + records.extend(paf_record_from_nam( nam1, &r1.name, refseq, @@ -177,7 +187,7 @@ pub fn map_paired_end_read( None, End::One, )); - records.push(paf_record_from_nam( + records.extend(paf_record_from_nam( nam2, &r2.name, refseq, @@ -189,6 +199,7 @@ pub fn map_paired_end_read( } nam_details1 += nam_details2; + (records, nam_details1.into()) } @@ -199,6 +210,7 @@ pub fn abundances_paired_end_read( r1: &SequenceRecord, r2: &SequenceRecord, index: &StrobemerIndex, + refseq: &RefSequence, abundances: &mut [f64], rescue_distance: usize, insert_size_distribution: &mut InsertSizeDistribution, @@ -236,7 +248,8 @@ pub fn abundances_paired_end_read( .count(); let weight = read_len as f64 / n_best as f64; for nam in &nams[0..n_best] { - abundances[nam.ref_id] += weight; + let reference_id = refseq.unflatten(nam.ref_start).0; + abundances[reference_id] += weight; } } } @@ -253,8 +266,8 @@ pub fn abundances_paired_end_read( score: _, } in &nam_pairs[..n_best] { - abundances[nam1.ref_id] += weight_r1; - abundances[nam2.ref_id] += weight_r2; + abundances[refseq.unflatten(nam1.ref_start).0] += weight_r1; + abundances[refseq.unflatten(nam2.ref_start).0] += weight_r2; } } } @@ -330,13 +343,13 @@ fn get_nam_pairs( split_nams_by_orientation_checked(nams2, details2.both_orientations); if !fwd1.is_empty() && !rev2.is_empty() { - fwd1.sort_unstable_by_key(|nam| (nam.ref_id, nam.projected_ref_start())); - rev2.sort_unstable_by_key(|nam| (nam.ref_id, nam.projected_ref_start())); + fwd1.sort_unstable_by_key(|nam| nam.projected_ref_start()); + rev2.sort_unstable_by_key(|nam| nam.projected_ref_start()); nam_pairs.extend(find_pairs(fwd1, rev2, mu, sigma, false)); } if !fwd2.is_empty() && !rev1.is_empty() { - fwd2.sort_unstable_by_key(|nam| (nam.ref_id, nam.projected_ref_start())); - rev1.sort_unstable_by_key(|nam| (nam.ref_id, nam.projected_ref_start())); + fwd2.sort_unstable_by_key(|nam| nam.projected_ref_start()); + rev1.sort_unstable_by_key(|nam| nam.projected_ref_start()); nam_pairs.extend(find_pairs(fwd2, rev1, mu, sigma, true)); } @@ -376,7 +389,7 @@ fn split_nams_by_orientation(nams: &mut [Nam]) -> (&mut [Nam], &mut [Nam]) { } /// Find most forward/revcomp pairs using a two-pointer scan. -/// Assumes both slices are sorted by (ref_id, projected_ref_start). +/// Assumes both slices are sorted by projected_ref_start. fn find_pairs(fwd: &[Nam], rev: &[Nam], mu: f32, sigma: f32, swap_order: bool) -> Vec { let mut out = Vec::new(); let max_dist = (mu + 10.0 * sigma).ceil() as usize; // distance cutoff from insert size distribution @@ -385,25 +398,28 @@ fn find_pairs(fwd: &[Nam], rev: &[Nam], mu: f32, sigma: f32, swap_order: bool) - for f in fwd { // Advance revcomp pointer to the first possible candidate - while rev_ptr < rev.len() - && (rev[rev_ptr].ref_id < f.ref_id - || rev[rev_ptr].ref_id == f.ref_id - && rev[rev_ptr].projected_ref_start() < f.projected_ref_start()) - { + // TODO this used to be: + // while rev_ptr < rev.len() + // && (rev[rev_ptr].ref_id < f.ref_id + // || rev[rev_ptr].ref_id == f.ref_id + // && rev[rev_ptr].projected_ref_start() < f.projected_ref_start()) + + while rev_ptr < rev.len() && rev[rev_ptr].projected_ref_start() < f.projected_ref_start() { rev_ptr += 1; } if rev_ptr == rev.len() { break; } - if rev[rev_ptr].ref_id > f.ref_id { - continue; - } + // TODO this used to be uncommented: + // if rev[rev_ptr].ref_id > f.ref_id { + // continue; + // } // Scan window of revcomp nams within distance limit. let mut best = None; let mut i = rev_ptr; while i < rev.len() - && rev[i].ref_id == f.ref_id + // && rev[i].ref_id == f.ref_id TODO && (rev[i].projected_ref_start() - f.projected_ref_start()) <= max_dist { let r = &rev[i]; diff --git a/src/mapper.rs b/src/mapper.rs index 82f9e099..6e372299 100644 --- a/src/mapper.rs +++ b/src/mapper.rs @@ -56,7 +56,6 @@ impl Default for MappingParameters { #[derive(Debug, Clone)] struct Alignment { - reference_id: usize, ref_start: usize, cigar: Cigar, edit_distance: usize, @@ -173,7 +172,9 @@ impl SamOutput { cigar.push(CigarOperation::Softclip, alignment.soft_clip_left); cigar.extend(&alignment.cigar); cigar.push(CigarOperation::Softclip, alignment.soft_clip_right); - let reference_name = Some(refseq.names[alignment.reference_id].clone()); + let (reference_id, reference_start) = refseq.unflatten(alignment.ref_start); + let reference_name = Some(refseq.names[reference_id].clone()); + let pos = Some(reference_start); // TODO let details = if self.details { Some(details) } else { None }; let cigar = if self.cigar_eqx { Some(cigar) @@ -189,7 +190,7 @@ impl SamOutput { query_name: record.name.clone(), flags, reference_name, - pos: Some(alignment.ref_start as u32), + pos, mapq, cigar, query_sequence: Some(query_sequence), @@ -286,15 +287,22 @@ impl SamOutput { if mate.is_revcomp { sam_records[i].flags |= MREVERSE; } + + // TODO This second lookup could be avoided because the two + // SamRecords contain the unflattened coordinates + let (mate_reference_id, mate_ref_start) = refseq.unflatten(mate.ref_start); + // PNEXT (position of mate) - sam_records[i].mate_pos = Some(mate.ref_start as u32); + sam_records[i].mate_pos = Some(mate_ref_start); // RNEXT (reference name of mate) - sam_records[i].mate_reference_name = Some(refseq.names[mate.reference_id].clone()); + sam_records[i].mate_reference_name = Some(refseq.names[mate_reference_id].clone()); if let Some(this) = alignments[i] { // both aligned - if this.reference_id == mate.reference_id { + // TODO see comment above: this lookup could be avoided + let this_reference_id = refseq.unflatten(this.ref_start).0; + if this_reference_id == mate_reference_id { // aligned to same reference sam_records[i].mate_reference_name = Some("=".to_string()); @@ -315,8 +323,8 @@ impl SamOutput { // mate-pair read whose mate is mapped, the unmapped read should have // RNAME and POS identical to its mate." sam_records[i].mate_reference_name = Some("=".to_string()); - sam_records[i].reference_name = Some(refseq.names[mate.reference_id].clone()); - sam_records[i].pos = Some(mate.ref_start as u32); + sam_records[i].reference_name = Some(refseq.names[mate_reference_id].clone()); + sam_records[i].pos = Some(mate_ref_start); } } else { // mate unmapped @@ -325,7 +333,8 @@ impl SamOutput { // Set RNAME and POS identical to mate if let Some(this) = alignments[i] { sam_records[i].mate_reference_name = Some("=".to_string()); - sam_records[i].mate_pos = Some(this.ref_start as u32); + let this_reference_start = refseq.unflatten(this.ref_start).1; + sam_records[i].mate_pos = Some(this_reference_start); } } } @@ -540,7 +549,7 @@ fn extend_seed( } else { read.seq() }; - let seq = &refseq.contig(nam.ref_id); + let seq = &refseq.sequence(); let projected_ref_start = nam.ref_start.saturating_sub(nam.query_start); let projected_ref_end = min(nam.ref_end + query.len() - nam.query_end, seq.len()); @@ -594,7 +603,6 @@ fn extend_seed( .anchors .iter() .map(|a| Anchor { - ref_id: a.ref_id, ref_start: a.ref_start - decode_start, query_start: a.query_start, }) @@ -612,7 +620,6 @@ fn extend_seed( ref_start: result_ref_start, length: info.ref_span(), is_revcomp: nam.is_revcomp, - reference_id: nam.ref_id, gapped, }) } @@ -1114,7 +1121,7 @@ fn rescue_align( ) }; - let ref_len = refseq.contig_len(mate_nam.ref_id); + let ref_len = refseq.total_length(); let ref_start = ref_start.min(ref_len); let ref_end = ref_end.min(ref_len); @@ -1122,7 +1129,7 @@ fn rescue_align( // std::cerr << "RESCUE: Caught Bug3! ref start: " << ref_start << " ref end: " << ref_end << " ref len: " << ref_len << std::endl; return None; } - let ref_segm = refseq.contig(mate_nam.ref_id).decode(ref_start, ref_end); + let ref_segm = refseq.decode(ref_start, ref_end); if !has_shared_substring(r_tmp, &ref_segm, k) { return None; @@ -1130,7 +1137,6 @@ fn rescue_align( let info = aligner.align(r_tmp, &ref_segm); if let Some(info) = info { Some(Alignment { - reference_id: mate_nam.ref_id, ref_start: ref_start + info.ref_start, edit_distance: info.edit_distance, soft_clip_left: info.query_start, @@ -1179,7 +1185,7 @@ fn is_proper_pair_opt( /// on the same reference, in opposite orientations, within the expected insert size. fn is_proper_pair(a1: &Alignment, a2: &Alignment, mu: f32, sigma: f32) -> PairStatus { let dist = a2.ref_start as isize - a1.ref_start as isize; - let same_reference = a1.reference_id == a2.reference_id; + let same_reference = true; // TODO a1.reference_id == a2.reference_id; let r1_r2 = !a1.is_revcomp && a2.is_revcomp && dist >= 0; // r1 ---> <---- r2 let r2_r1 = !a2.is_revcomp && a1.is_revcomp && dist <= 0; // r2 ---> <---- r1 let rel_orientation_good = r1_r2 || r2_r1; @@ -1193,7 +1199,7 @@ fn is_proper_pair(a1: &Alignment, a2: &Alignment, mu: f32, sigma: f32) -> PairSt } pub fn is_proper_nam_pair(nam1: &Nam, nam2: &Nam, mu: f32, sigma: f32) -> bool { - if nam1.ref_id != nam2.ref_id || nam1.is_revcomp == nam2.is_revcomp { + if nam1.is_revcomp == nam2.is_revcomp { return false; } let r1_ref_start = nam1.projected_ref_start(); @@ -1332,9 +1338,7 @@ fn deduplicate_scored_pairs(pairs: &mut Vec) { (None, Some(_)) => false, (Some(_), None) => false, (None, None) => true, - (Some(a1), Some(a2)) => { - a1.ref_start == a2.ref_start && a1.reference_id == a2.reference_id - } + (Some(a1), Some(a2)) => a1.ref_start == a2.ref_start, } } let mut k = 0; @@ -1468,14 +1472,8 @@ mod tests { use super::*; use crate::cigar::Cigar; - pub fn make_alignment( - reference_id: usize, - ref_start: usize, - score: u32, - is_revcomp: bool, - ) -> Alignment { + pub fn make_alignment(ref_start: usize, score: u32, is_revcomp: bool) -> Alignment { Alignment { - reference_id, ref_start, score, is_revcomp, @@ -1490,7 +1488,6 @@ mod tests { fn dummy_alignment() -> Alignment { Alignment { - reference_id: 0, ref_start: 0, cigar: Cigar::default(), edit_distance: 0, @@ -1531,12 +1528,10 @@ mod tests { #[test] fn deduplicate_scored_pairs_works() { let a1 = Some(Alignment { - reference_id: 0, ref_start: 1906, ..dummy_alignment() }); let a2 = Some(Alignment { - reference_id: 0, ref_start: 123, ..dummy_alignment() }); @@ -1589,8 +1584,8 @@ mod tests { fn is_proper_pair_within_distance() { let mu = 300.0_f32; let sigma = 50.0_f32; - let a1 = make_alignment(0, 100, 40, false); - let a2 = make_alignment(0, 350, 40, true); + let a1 = make_alignment(100, 40, false); + let a2 = make_alignment(350, 40, true); assert_eq!(is_proper_pair(&a1, &a2, mu, sigma), PairStatus::Proper); } @@ -1598,8 +1593,8 @@ mod tests { fn is_proper_pair_rev_within_distance() { let mu = 300.0_f32; let sigma = 50.0_f32; - let a1 = make_alignment(0, 350, 40, true); - let a2 = make_alignment(0, 100, 40, false); + let a1 = make_alignment(350, 40, true); + let a2 = make_alignment(100, 40, false); assert_eq!(is_proper_pair(&a1, &a2, mu, sigma), PairStatus::Proper); } @@ -1607,8 +1602,8 @@ mod tests { fn is_proper_pair_incompatible() { let mu = 300.0_f32; let sigma = 50.0_f32; - let a1 = make_alignment(0, 100, 40, false); - let a2 = make_alignment(0, 350, 40, false); + let a1 = make_alignment(100, 40, false); + let a2 = make_alignment(350, 40, false); assert_eq!(is_proper_pair(&a1, &a2, mu, sigma), PairStatus::NotProper); } @@ -1616,8 +1611,8 @@ mod tests { fn is_proper_pair_too_far() { let mu = 300.0_f32; let sigma = 50.0_f32; - let a1 = make_alignment(0, 100, 40, false); - let a2 = make_alignment(0, 10_000, 40, true); + let a1 = make_alignment(100, 40, false); + let a2 = make_alignment(10_000, 40, true); assert_eq!(is_proper_pair(&a1, &a2, mu, sigma), PairStatus::NotProper); } } diff --git a/src/nam.rs b/src/nam.rs index f41d0bc4..e10357a4 100644 --- a/src/nam.rs +++ b/src/nam.rs @@ -24,7 +24,6 @@ pub struct Nam { pub query_start: usize, pub query_end: usize, pub matching_bases: usize, - pub ref_id: usize, pub score: f32, pub is_revcomp: bool, pub anchors: Vec, @@ -47,12 +46,8 @@ impl Nam { /// reference by comparing the nucleotide sequences of the first and last /// strobe (taking orientation into account). pub fn is_consistent(&self, read: &Read, refseq: &RefSequence, k: usize) -> bool { - let ref_start_kmer = refseq - .contig(self.ref_id) - .decode(self.ref_start, self.ref_start + k); - let ref_end_kmer = refseq - .contig(self.ref_id) - .decode(self.ref_end - k, self.ref_end); + let ref_start_kmer = refseq.decode(self.ref_start, self.ref_start + k); + let ref_end_kmer = refseq.decode(self.ref_end - k, self.ref_end); let seq = if self.is_revcomp { read.rc() @@ -70,8 +65,7 @@ impl Display for Nam { fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result { write!( f, - "Nam(ref_id={}, query: {}..{}, ref: {}..{}, rc={}, score={})", - self.ref_id, + "Nam(query: {}..{}, ref: {}..{}, rc={}, score={})", self.query_start, self.query_end, self.ref_start, @@ -92,12 +86,8 @@ impl Display for Nam { /// in place and return true. /// - If first and last strobe do not match consistently, return false. pub fn reverse_nam_if_needed(nam: &mut Nam, read: &Read, refseq: &RefSequence, k: usize) -> bool { - let ref_start_kmer = refseq - .contig(nam.ref_id) - .decode(nam.ref_start, nam.ref_start + k); - let ref_end_kmer = refseq - .contig(nam.ref_id) - .decode(nam.ref_end - k, nam.ref_end); + let ref_start_kmer = refseq.decode(nam.ref_start, nam.ref_start + k); + let ref_end_kmer = refseq.decode(nam.ref_end - k, nam.ref_end); let (seq, seq_rc) = if nam.is_revcomp { (read.rc(), read.seq()) diff --git a/src/packed_seq.rs b/src/packed_seq.rs index f8f45485..481ca8f1 100644 --- a/src/packed_seq.rs +++ b/src/packed_seq.rs @@ -14,6 +14,14 @@ pub struct PackedSeq { len: usize, } +#[derive(Debug)] +#[allow(clippy::len_without_is_empty)] +pub struct PackedSeqSlice<'a> { + seq: &'a PackedSeq, + start: usize, + len: usize, +} + const fn make_encode_table() -> [u8; 256] { let mut t = [0u8; 256]; t[b'C' as usize] = 1; @@ -91,6 +99,12 @@ impl PackedSeq { self.len += 1; } + pub fn extend>(&mut self, bytes: I) { + for b in bytes.into_iter() { + self.push(b); + } + } + pub fn len(&self) -> usize { self.len } @@ -124,3 +138,42 @@ impl PackedSeq { self.decode(0, self.len) } } + +impl<'a> PackedSeqSlice<'a> { + pub fn new(seq: &'a PackedSeq, start: usize, len: usize) -> Self { + PackedSeqSlice { seq, start, len } + } + + pub fn nucleotide_bits(&self, i: usize) -> u8 { + self.seq.nucleotide_bits(self.start + i) + } + + pub fn len(&self) -> usize { + self.len + } + + /// Decode bases `start..end` as ASCII bytes. + pub fn decode(&self, start: usize, end: usize) -> Vec { + (start..end) + .map(|i| self.seq.decode_at(self.start + i)) + .collect() + } + + /// Decode bases `start..end` as ASCII bytes. + pub fn decode_all(&self) -> Vec { + (self.start..self.start + self.len) + .map(|i| self.seq.decode_at(i)) + .collect() + } + + pub fn to_owned(&self) -> PackedSeq { + // This is probably inefficient because it unpacks and then + // re-packs the sequence, but it is only used for testing + // with tiny references. + let mut ps = PackedSeq::new(); + let decoded = self.decode_all(); + ps.extend(decoded); + + ps + } +} diff --git a/src/piecewisealigner.rs b/src/piecewisealigner.rs index a7b4aa59..74b3f168 100644 --- a/src/piecewisealigner.rs +++ b/src/piecewisealigner.rs @@ -1252,57 +1252,46 @@ mod tests { fn remove_spurious_anchors_control() { let expected = vec![ Anchor { - ref_id: 0, ref_start: 0, query_start: 0, }, Anchor { - ref_id: 0, ref_start: 10, query_start: 10, }, Anchor { - ref_id: 0, ref_start: 20, query_start: 20, }, Anchor { - ref_id: 0, ref_start: 30, query_start: 30, }, Anchor { - ref_id: 0, ref_start: 40, query_start: 40, }, Anchor { - ref_id: 0, ref_start: 50, query_start: 50, }, Anchor { - ref_id: 0, ref_start: 60, query_start: 60, }, Anchor { - ref_id: 0, ref_start: 70, query_start: 70, }, Anchor { - ref_id: 0, ref_start: 80, query_start: 80, }, Anchor { - ref_id: 0, ref_start: 90, query_start: 90, }, Anchor { - ref_id: 0, ref_start: 100, query_start: 100, }, @@ -1316,57 +1305,46 @@ mod tests { fn remove_spurious_anchors_inside_trim() { let mut result = vec![ Anchor { - ref_id: 0, ref_start: 0, query_start: 0, }, Anchor { - ref_id: 0, ref_start: 10, query_start: 10, }, Anchor { - ref_id: 0, ref_start: 20, query_start: 20, }, Anchor { - ref_id: 0, ref_start: 130, query_start: 30, }, Anchor { - ref_id: 0, ref_start: 140, query_start: 40, }, Anchor { - ref_id: 0, ref_start: 150, query_start: 50, }, Anchor { - ref_id: 0, ref_start: 60, query_start: 60, }, Anchor { - ref_id: 0, ref_start: 70, query_start: 50, }, Anchor { - ref_id: 0, ref_start: 80, query_start: 80, }, Anchor { - ref_id: 0, ref_start: 90, query_start: 90, }, Anchor { - ref_id: 0, ref_start: 100, query_start: 100, }, @@ -1374,37 +1352,30 @@ mod tests { remove_spurious_anchors(&mut result); let expected = vec![ Anchor { - ref_id: 0, ref_start: 0, query_start: 0, }, Anchor { - ref_id: 0, ref_start: 10, query_start: 10, }, Anchor { - ref_id: 0, ref_start: 20, query_start: 20, }, Anchor { - ref_id: 0, ref_start: 60, query_start: 60, }, Anchor { - ref_id: 0, ref_start: 80, query_start: 80, }, Anchor { - ref_id: 0, ref_start: 90, query_start: 90, }, Anchor { - ref_id: 0, ref_start: 100, query_start: 100, }, @@ -1416,57 +1387,46 @@ mod tests { fn remove_spurious_anchors_ends_trim() { let mut result = vec![ Anchor { - ref_id: 0, ref_start: 5, query_start: 0, }, Anchor { - ref_id: 0, ref_start: 15, query_start: 10, }, Anchor { - ref_id: 0, ref_start: 20, query_start: 20, }, Anchor { - ref_id: 0, ref_start: 30, query_start: 30, }, Anchor { - ref_id: 0, ref_start: 40, query_start: 40, }, Anchor { - ref_id: 0, ref_start: 50, query_start: 50, }, Anchor { - ref_id: 0, ref_start: 60, query_start: 60, }, Anchor { - ref_id: 0, ref_start: 70, query_start: 70, }, Anchor { - ref_id: 0, ref_start: 80, query_start: 80, }, Anchor { - ref_id: 0, ref_start: 90, query_start: 95, }, Anchor { - ref_id: 0, ref_start: 100, query_start: 105, }, @@ -1474,37 +1434,30 @@ mod tests { remove_spurious_anchors(&mut result); let expected = vec![ Anchor { - ref_id: 0, ref_start: 20, query_start: 20, }, Anchor { - ref_id: 0, ref_start: 30, query_start: 30, }, Anchor { - ref_id: 0, ref_start: 40, query_start: 40, }, Anchor { - ref_id: 0, ref_start: 50, query_start: 50, }, Anchor { - ref_id: 0, ref_start: 60, query_start: 60, }, Anchor { - ref_id: 0, ref_start: 70, query_start: 70, }, Anchor { - ref_id: 0, ref_start: 80, query_start: 80, }, @@ -1529,22 +1482,18 @@ mod tests { let refseq = b"AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA"; let anchors = vec![ Anchor { - ref_id: 0, ref_start: 40, query_start: 30, }, Anchor { - ref_id: 0, ref_start: 30, query_start: 20, }, Anchor { - ref_id: 0, ref_start: 20, query_start: 10, }, Anchor { - ref_id: 0, ref_start: 10, query_start: 0, }, @@ -1578,17 +1527,14 @@ mod tests { let refseq = b"TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT"; let anchors = vec![ Anchor { - ref_id: 0, ref_start: 30, query_start: 20, }, Anchor { - ref_id: 0, ref_start: 20, query_start: 20, }, Anchor { - ref_id: 0, ref_start: 10, query_start: 10, }, @@ -1616,37 +1562,30 @@ mod tests { b"AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA"; let anchors = vec![ Anchor { - ref_id: 0, ref_start: 40, query_start: 30, }, Anchor { - ref_id: 0, ref_start: 37, query_start: 27, }, Anchor { - ref_id: 0, ref_start: 35, query_start: 25, }, Anchor { - ref_id: 0, ref_start: 30, query_start: 20, }, Anchor { - ref_id: 0, ref_start: 25, query_start: 20, }, Anchor { - ref_id: 0, ref_start: 20, query_start: 15, }, Anchor { - ref_id: 0, ref_start: 10, query_start: 5, }, @@ -1680,17 +1619,14 @@ mod tests { let refseq = b"TTTTTAAAAATTTTTAAAAATTTTTAAAAATTTTTAAAAATTTTTAAAAA"; let anchors = vec![ Anchor { - ref_id: 0, ref_start: 35, query_start: 36, }, Anchor { - ref_id: 0, ref_start: 15, query_start: 14, }, Anchor { - ref_id: 0, ref_start: 5, query_start: 5, }, diff --git a/src/refseq.rs b/src/refseq.rs index 5d5d2923..09ff6fcb 100644 --- a/src/refseq.rs +++ b/src/refseq.rs @@ -1,37 +1,132 @@ -use crate::packed_seq::PackedSeq; +use crate::packed_seq::{PackedSeq, PackedSeqSlice}; + +/// A position on a contig. This separate type is here to prevent confusion +/// with "flat" reference coordinates, which are used everywhere else. +#[derive(Debug, Clone, Copy, Default, PartialEq)] +pub struct ContigPosition(pub u32); #[derive(Default, Debug, Clone)] pub struct RefSequence { + /// Contig names pub names: Vec, - sequences: Vec, + + /// Concatenated sequence of all contigs + sequence: PackedSeq, + + /// Start positions for all contigs. `starts[i]` is the start position of + /// contig with index `i` within `sequence`. This contains one more item + /// than there are contigs, which is set to the length of `sequence`. + starts: Vec, } impl RefSequence { - pub fn new() -> Self { - RefSequence { - names: vec![], - sequences: vec![], - } + pub fn new(sequence: PackedSeq, starts: Vec, names: Vec) -> Self { + assert_eq!(starts.len(), names.len()); + let total_length = sequence.len(); + let mut refseq = RefSequence { + sequence, + starts, + names, + }; + refseq.starts.push(total_length); + + refseq + } + + pub fn sequence(&self) -> &PackedSeq { + &self.sequence } - pub fn push(&mut self, name: String, seq: PackedSeq) { - self.names.push(name); - self.sequences.push(seq); + /// Returns (reference_index, reference_start) + pub fn unflatten(&self, start: usize) -> (usize, ContigPosition) { + let ref_index = self.starts[1..].partition_point(|&x| x <= start); + let start_within_contig = start - self.starts[ref_index]; + + (ref_index, ContigPosition(start_within_contig as u32)) + } + + pub fn decode(&self, start: usize, end: usize) -> Vec { + self.sequence.decode(start, end) } - pub fn contig(&self, index: usize) -> &PackedSeq { - &self.sequences[index] + pub fn contig<'a>(&'a self, index: usize) -> PackedSeqSlice<'a> { + let start = self.starts[index]; + let len = self.starts[index + 1] - start; + + PackedSeqSlice::new(&self.sequence, start, len) } - pub fn contig_len(&self, index: usize) -> usize { - self.sequences[index].len() + pub fn contig_start(&self, index: usize) -> usize { + self.starts[index] } pub fn max_contig_len(&self) -> Option { - self.sequences.iter().map(|seq| seq.len()).max() + (0..self.names.len()).map(|i| self.contig(i).len()).max() } pub fn total_length(&self) -> usize { - self.sequences.iter().map(|seq| seq.len()).sum() + *self.starts.last().unwrap() + } +} + +#[cfg(test)] +mod test { + use super::*; + + fn testref() -> RefSequence { + let mut sequence = PackedSeq::new(); + sequence.extend(b"AAAACCCGGT".to_vec()); + let names = vec![ + "n1".to_string(), + "n2".to_string(), + "n3".to_string(), + "n4".to_string(), + ]; + let starts = vec![0, 4, 7, 9]; + RefSequence::new(sequence, starts, names) + } + + #[test] + fn new() { + let rs = testref(); + + assert_eq!(rs.total_length(), 10); + assert_eq!(rs.max_contig_len(), Some(4)); + + assert_eq!(rs.contig(0).len(), 4); + assert_eq!(rs.contig(1).len(), 3); + assert_eq!(rs.contig(2).len(), 2); + assert_eq!(rs.contig(3).len(), 1); + assert_eq!(rs.contig(0).len(), 4); + } + + #[test] + fn decode() { + let rs = testref(); + + assert_eq!(rs.decode(0, 0), b""); + assert_eq!(rs.decode(0, 1), b"A"); + assert_eq!(rs.decode(0, 2), b"AA"); + assert_eq!(rs.decode(0, 3), b"AAA"); + assert_eq!(rs.decode(4, 5), b"C"); + assert_eq!(rs.decode(4, 6), b"CC"); + assert_eq!(rs.decode(4, 7), b"CCC"); + assert_eq!(rs.decode(9, 10), b"T"); + } + + #[test] + fn unflatten() { + let rs = testref(); + + assert_eq!(rs.unflatten(0), (0, ContigPosition(0))); + assert_eq!(rs.unflatten(1), (0, ContigPosition(1))); + assert_eq!(rs.unflatten(2), (0, ContigPosition(2))); + assert_eq!(rs.unflatten(3), (0, ContigPosition(3))); + assert_eq!(rs.unflatten(4), (1, ContigPosition(0))); + assert_eq!(rs.unflatten(5), (1, ContigPosition(1))); + assert_eq!(rs.unflatten(6), (1, ContigPosition(2))); + assert_eq!(rs.unflatten(7), (2, ContigPosition(0))); + assert_eq!(rs.unflatten(8), (2, ContigPosition(1))); + assert_eq!(rs.unflatten(9), (3, ContigPosition(0))); } } diff --git a/src/seeding/strobes.rs b/src/seeding/strobes.rs index 1c340338..189b8b50 100644 --- a/src/seeding/strobes.rs +++ b/src/seeding/strobes.rs @@ -179,10 +179,10 @@ impl> Iterator for RandstrobeIterator { #[cfg(test)] mod test { use super::*; + use crate::io::fasta::read_ref; use crate::packed_seq::PackedSeq; - use crate::seeding::SeedingParameters; - use crate::seeding::SyncmerIterator; + use crate::seeding::{SeedingParameters, SyncmerIterator}; fn read_phix() -> PackedSeq { read_ref("tests/phix.fasta").unwrap().contig(0).to_owned() diff --git a/src/seeding/syncmers.rs b/src/seeding/syncmers.rs index aac65afd..8a40b3fd 100644 --- a/src/seeding/syncmers.rs +++ b/src/seeding/syncmers.rs @@ -3,7 +3,7 @@ use std::collections::VecDeque; use super::InvalidSeedingParameter; use super::hash::xxh64; -use crate::packed_seq::PackedSeq; +use crate::packed_seq::{PackedSeq, PackedSeqSlice}; /// Trait for types that can serve as a sequence source for `SyncmerIterator`. /// @@ -27,7 +27,7 @@ impl SeqAccess for &[u8] { } } -impl SeqAccess for &PackedSeq { +impl<'a> SeqAccess for &'a PackedSeqSlice<'a> { /// Direct 2-bit extract - no table lookups, never returns 4. /// `**self` reaches `PackedSeq` so Rust resolves the inherent method, /// not this trait method (no recursion). @@ -40,6 +40,16 @@ impl SeqAccess for &PackedSeq { } } +impl SeqAccess for &PackedSeq { + fn nucleotide_bits(&self, i: usize) -> u8 { + (**self).nucleotide_bits(i) + } + + fn len(&self) -> usize { + (**self).len() + } +} + #[derive(Debug, Clone, Copy, PartialEq, Eq)] pub struct Syncmer { hash: u64,