Skip to content

Concatenate the reference contigs into one "flat" reference #607

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
marcelm wants to merge 1 commit into
base: main
Choose a base branch
from
Open

Concatenate the reference contigs into one "flat" reference #607

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
40 changes: 18 additions & 22 deletions src/chainer.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
}
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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!(
Expand Down Expand Up @@ -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(),
});
Expand All @@ -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,
});
Expand Down Expand Up @@ -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,
Expand All @@ -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);
Expand All @@ -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;
Expand All @@ -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.
Expand Down
47 changes: 19 additions & 28 deletions src/index.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -35,28 +34,22 @@ 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,
}
}

pub fn hash(&self) -> RandstrobeHash {
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 {
Expand Down Expand Up @@ -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,
}

Expand Down Expand Up @@ -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
Expand All @@ -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"
Expand Down Expand Up @@ -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);
}

Expand Down Expand Up @@ -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);

Expand Down
32 changes: 16 additions & 16 deletions src/indexer.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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};

Expand Down Expand Up @@ -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(),
);
}
});
}
Expand All @@ -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;
}
Expand All @@ -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());
}
Expand All @@ -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;
}
});
Expand All @@ -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,
Expand Down Expand Up @@ -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);
Expand All @@ -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), &parameters), 1090);
assert_eq!(count_randstrobes(&refseq.contig(0), &parameters), 1090);
}
}
Loading
Loading