use crate as libradicl;
use self::libradicl::rad_types::{MappedFragmentOrientation, RadIntId};
use self::libradicl::record::AlevinFryReadRecord;
use self::libradicl::record::AtacSeqReadRecord;
use self::libradicl::record::{
CollatableMappedRecord, CollatableRecordHeader, ConvertiblePrimitiveInteger, KnownSize,
MappedRecord, RecordHeader,
};
use self::libradicl::schema::{CollateKey, TempCellInfo};
#[allow(unused_imports)]
use ahash::{AHasher, RandomState};
use bio_types::strand::*;
use scroll::Pread;
use serde::{Deserialize, Serialize};
use std::collections::HashMap;
use std::fs::File;
use std::io::{BufReader, Cursor, Read, Seek, SeekFrom};
use std::io::{BufWriter, Write};
use std::sync::atomic::{AtomicU32, AtomicU64, Ordering};
use std::sync::{Arc, Mutex};
use std::vec::Vec;
pub mod chunk;
pub mod codec;
pub mod collation;
pub mod constants;
pub mod exit_codes;
pub mod header;
pub mod io;
pub mod rad_types;
pub mod readers;
pub mod record;
pub mod schema;
pub mod unmapped;
pub mod utils;
pub mod writers;
pub use chunk::ChunkBuf;
pub use codec::{CHUNK_CODEC_TAG, ChunkCodec};
pub use libradicl_macros::UmiTagged;
pub use writers::{ConcurrentChunkWriter, RadFileWriter};
#[macro_use]
mod macros;
static LIB_NAME: &str = "libradicl";
pub fn lib_name() -> &'static str {
LIB_NAME
}
#[derive(Serialize, Deserialize, Debug)]
#[allow(dead_code)]
pub struct BarcodeLookupMap {
pub barcodes: Vec<u64>,
offsets: Vec<usize>,
bclen: u32,
prefix_len: u32,
suffix_len: u32,
}
impl BarcodeLookupMap {
pub fn new(mut kv: Vec<u64>, bclen: u32) -> BarcodeLookupMap {
let prefix_len = bclen.div_ceil(2) as u64;
let suffix_len = bclen - prefix_len as u32;
let _prefix_bits = 2 * prefix_len;
let suffix_bits = 2 * suffix_len;
kv.sort_unstable();
let pref_mask = ((4usize.pow(prefix_len as u32) - 1) as u64) << (suffix_bits);
let mut offsets = vec![0; 4usize.pow(prefix_len as u32) + 1];
let mut prev_ind = 0xFFFF;
for (n, &v) in kv.iter().enumerate() {
let ind = ((v & pref_mask) >> (suffix_bits)) as usize;
if ind != prev_ind {
for item in offsets.iter_mut().take(ind).skip(prev_ind + 1) {
*item = n;
}
offsets[ind] = n;
prev_ind = ind;
}
}
for item in offsets.iter_mut().skip(prev_ind + 1) {
*item = kv.len();
}
BarcodeLookupMap {
barcodes: kv,
offsets,
bclen,
prefix_len: prefix_len as u32,
suffix_len,
}
}
#[allow(dead_code)]
pub fn barcode_for_idx(&self, idx: usize) -> u64 {
self.barcodes[idx]
}
pub fn find_exact(&self, query: u64) -> Option<usize> {
let mut ret: Option<usize> = None;
let suffix_bits = 2 * self.suffix_len;
let query_pref = query >> suffix_bits;
let qrange = std::ops::Range {
start: self.offsets[query_pref as usize],
end: self.offsets[(query_pref + 1) as usize],
};
let qs = qrange.start;
if let Ok(res) = self.barcodes[qrange].binary_search(&query) {
ret = Some(qs + res);
}
ret
}
pub fn find_neighbors(&self, query: u64, try_exact: bool) -> (Option<usize>, usize) {
let mut ret: Option<usize> = None;
let pref_bits = 2 * self.prefix_len;
let suffix_bits = 2 * self.suffix_len;
let mut query_pref = query >> suffix_bits;
let mut num_neighbors = 0usize;
let qrange = std::ops::Range {
start: self.offsets[query_pref as usize],
end: self.offsets[(query_pref + 1) as usize],
};
let qs = qrange.start;
if try_exact {
if let Ok(res) = self.barcodes[qrange.clone()].binary_search(&query) {
ret = Some(qs + res);
num_neighbors += 1;
return (ret, num_neighbors);
}
}
if !(std::ops::Range::<usize>::is_empty(&qrange)) {
let qs = qrange.start;
for i in (0..suffix_bits).step_by(2) {
let bit_mask = 3 << (i);
for nmod in 1..4 {
let nucl = 0x3 & ((query >> i) + nmod);
let nquery = (query & (!bit_mask)) | (nucl << i);
if let Ok(res) = self.barcodes[qrange.clone()].binary_search(&nquery) {
ret = Some(qs + res);
num_neighbors += 1;
if num_neighbors >= 2 {
return (ret, num_neighbors);
}
}
}
}
}
{
for i in (suffix_bits..(suffix_bits + pref_bits)).step_by(2) {
let bit_mask = 3 << i;
for nmod in 1..4 {
let nucl = 0x3 & ((query >> i) + nmod);
let nquery = (query & (!bit_mask)) | (nucl << i);
query_pref = nquery >> suffix_bits;
let qrange = std::ops::Range {
start: self.offsets[query_pref as usize],
end: self.offsets[(query_pref + 1) as usize],
};
let qs = qrange.start;
if let Ok(res) = self.barcodes[qrange].binary_search(&nquery) {
ret = Some(qs + res);
num_neighbors += 1;
if num_neighbors >= 2 {
return (ret, num_neighbors);
}
}
}
}
}
(ret, num_neighbors)
}
}
#[derive(Debug, Clone)]
pub struct GlobalEqCellList {
cell_ids: Vec<usize>,
count: u32,
}
impl GlobalEqCellList {
pub fn from_umi_and_count(bc_mer: usize, count: u32) -> GlobalEqCellList {
let mut cc = GlobalEqCellList {
cell_ids: Vec::new(),
count: 0,
};
cc.cell_ids.push(bc_mer);
cc.count += count;
cc
}
pub fn add_element(&mut self, bc_mer: usize, count: u32) {
self.cell_ids.push(bc_mer);
self.count += count;
}
}
pub fn collate_temporary_bucket_twopass_generic<
B: ConvertiblePrimitiveInteger,
T: Read + Seek,
U: Write,
R: MappedRecord + KnownSize + CollatableMappedRecord<B>,
>(
reader: &mut BufReader<T>,
rec_context: &<R as MappedRecord>::ParsingContext,
nrec: u32,
owriter: &Mutex<U>,
compress: bool,
cb_byte_map: &mut HashMap<u64, TempCellInfo, ahash::RandomState>,
) -> usize
where
u64: From<B>,
{
let mut tbuf = vec![0u8; 65536];
let mut total_bytes = 0usize;
let chunk_header_size = 2 * std::mem::size_of::<u32>() as u64;
let size_of_aln = R::nbytes_aln(rec_context);
let calc_record_bytes = |num_aln: usize| -> usize { R::nbytes(num_aln as u32, rec_context) };
for _ in 0..(nrec as usize) {
let tup =
<R as CollatableMappedRecord<B>>::from_bytes_collatable_header(reader, rec_context)
.expect("can read header");
let v = cb_byte_map
.entry(tup.collation_group_key(rec_context))
.or_insert(TempCellInfo {
offset: chunk_header_size,
nbytes: chunk_header_size as u32,
nrec: 0_u32,
});
let na = tup.naln() as usize;
let req_size = size_of_aln * na;
if tbuf.len() < req_size {
tbuf.resize(req_size, 0);
}
reader.read_exact(&mut tbuf[0..(size_of_aln * na)]).unwrap();
let nbytes = calc_record_bytes(na);
v.offset += nbytes as u64;
v.nbytes += nbytes as u32;
v.nrec += 1;
total_bytes += nbytes;
}
total_bytes += cb_byte_map.len() * chunk_header_size as usize;
let mut output_buffer = Cursor::new(vec![0u8; total_bytes]);
let mut next_offset = 0u64;
for (_, v) in cb_byte_map.iter_mut() {
output_buffer.set_position(next_offset);
let cell_bytes = v.nbytes;
let cell_rec = v.nrec;
output_buffer.write_all(&cell_bytes.to_le_bytes()).unwrap();
output_buffer.write_all(&cell_rec.to_le_bytes()).unwrap();
v.offset = output_buffer.position();
let nbytes = v.nbytes as u64;
next_offset += nbytes;
}
reader
.get_mut()
.seek(SeekFrom::Start(0))
.expect("could not get read pointer.");
for _ in 0..(nrec as usize) {
let tup =
<R as CollatableMappedRecord<B>>::from_bytes_collatable_header(reader, rec_context)
.expect("can read header");
if let Some(v) = cb_byte_map.get_mut(&tup.collation_group_key(rec_context)) {
output_buffer.set_position(v.offset);
let na = tup.naln() as usize;
tup.write_fields(&mut output_buffer, rec_context)
.expect("could write header");
let bytes_for_aln_rec = R::nbytes_aln(rec_context);
reader
.read_exact(&mut tbuf[0..(bytes_for_aln_rec * na)])
.unwrap();
output_buffer
.write_all(&tbuf[..(bytes_for_aln_rec * na)])
.unwrap();
v.offset = output_buffer.position();
} else {
panic!("should not have any barcodes we can't find");
}
}
output_buffer.set_position(0);
if compress {
let mut compressed_output =
snap::write::FrameEncoder::new(Cursor::new(Vec::<u8>::with_capacity(total_bytes)));
compressed_output
.write_all(output_buffer.get_ref())
.expect("could not compress the output chunk.");
output_buffer = compressed_output
.into_inner()
.expect("couldn't unwrap the FrameEncoder.");
output_buffer.set_position(0);
}
owriter
.lock()
.unwrap()
.write_all(output_buffer.get_ref())
.unwrap();
cb_byte_map.len()
}
#[deprecated(
since = "0.10.0",
note = "This function is highly-specalized and works only with the AlevinFryReadRecordT<u64> type. \
This function has been deprecated in favor of the more generic `collate_temporary_bucket_twopass_generic`. \
Please use that function instead."
)]
pub fn collate_temporary_bucket_twopass<T: Read + Seek, U: Write>(
reader: &mut BufReader<T>,
bct: &RadIntId,
umit: &RadIntId,
nrec: u32,
owriter: &Mutex<U>,
compress: bool,
cb_byte_map: &mut HashMap<u64, TempCellInfo, ahash::RandomState>,
) -> usize {
let mut tbuf = vec![0u8; 65536];
let mut total_bytes = 0usize;
let chunk_header_size = 2 * std::mem::size_of::<u32>() as u64;
let size_of_u32 = std::mem::size_of::<u32>();
let size_of_bc = bct.bytes_for_type();
let size_of_umi = umit.bytes_for_type();
let calc_record_bytes = |num_aln: usize| -> usize {
size_of_u32 + size_of_bc + size_of_umi + (size_of_u32 * num_aln)
};
for _ in 0..(nrec as usize) {
let tup = AlevinFryReadRecord::from_bytes_record_header(reader, bct, umit);
let v = cb_byte_map.entry(tup.0).or_insert(TempCellInfo {
offset: chunk_header_size,
nbytes: chunk_header_size as u32,
nrec: 0_u32,
});
let na = tup.2 as usize;
let req_size = size_of_u32 * na;
if tbuf.len() < req_size {
tbuf.resize(req_size, 0);
}
reader.read_exact(&mut tbuf[0..(size_of_u32 * na)]).unwrap();
let nbytes = calc_record_bytes(na);
v.offset += nbytes as u64;
v.nbytes += nbytes as u32;
v.nrec += 1;
total_bytes += nbytes;
}
total_bytes += cb_byte_map.len() * chunk_header_size as usize;
let mut output_buffer = Cursor::new(vec![0u8; total_bytes]);
let mut next_offset = 0u64;
for (_, v) in cb_byte_map.iter_mut() {
output_buffer.set_position(next_offset);
let cell_bytes = v.nbytes;
let cell_rec = v.nrec;
output_buffer.write_all(&cell_bytes.to_le_bytes()).unwrap();
output_buffer.write_all(&cell_rec.to_le_bytes()).unwrap();
v.offset = output_buffer.position();
let nbytes = v.nbytes as u64;
next_offset += nbytes;
}
reader
.get_mut()
.seek(SeekFrom::Start(0))
.expect("could not get read pointer.");
for _ in 0..(nrec as usize) {
let tup = AlevinFryReadRecord::from_bytes_record_header(reader, bct, umit);
if let Some(v) = cb_byte_map.get_mut(&tup.0) {
output_buffer.set_position(v.offset);
let na = tup.2 as usize;
let nau32 = na as u32;
output_buffer.write_all(&nau32.to_le_bytes()).unwrap();
bct.write_to(tup.0, &mut output_buffer).unwrap();
umit.write_to(tup.1, &mut output_buffer).unwrap();
reader.read_exact(&mut tbuf[0..(size_of_u32 * na)]).unwrap();
output_buffer
.write_all(&tbuf[..(size_of_u32 * na)])
.unwrap();
v.offset = output_buffer.position();
} else {
panic!("should not have any barcodes we can't find");
}
}
output_buffer.set_position(0);
if compress {
let mut compressed_output =
snap::write::FrameEncoder::new(Cursor::new(Vec::<u8>::with_capacity(total_bytes)));
compressed_output
.write_all(output_buffer.get_ref())
.expect("could not compress the output chunk.");
output_buffer = compressed_output
.into_inner()
.expect("couldn't unwrap the FrameEncoder.");
output_buffer.set_position(0);
}
owriter
.lock()
.unwrap()
.write_all(output_buffer.get_ref())
.unwrap();
cb_byte_map.len()
}
pub fn collate_temporary_bucket_twopass_atac<T: Read + Seek, U: Write>(
reader: &mut BufReader<T>,
bct: &RadIntId,
nrec: u32,
owriter: &Mutex<U>,
compress: bool,
cb_byte_map: &mut HashMap<u64, TempCellInfo, ahash::RandomState>,
) -> usize {
let mut tbuf = vec![0u8; 65536];
let mut total_bytes = 0usize;
let header_size = 2 * std::mem::size_of::<u32>() as u64;
let size_of_bc = bct.bytes_for_type();
let size_of_rec = std::mem::size_of::<u32>()
+ std::mem::size_of::<u8>()
+ std::mem::size_of::<u32>()
+ std::mem::size_of::<u16>();
let size_of_u32 = std::mem::size_of::<u32>();
let calc_record_bytes =
|num_aln: usize| -> usize { size_of_u32 + size_of_bc + (size_of_rec * num_aln) };
for _ in 0..(nrec as usize) {
let tup = AtacSeqReadRecord::from_bytes_record_header(reader, bct);
let v = cb_byte_map.entry(tup.0).or_insert(TempCellInfo {
offset: header_size,
nbytes: header_size as u32,
nrec: 0_u32,
});
let na = tup.1 as usize;
let req_size = size_of_rec * na;
if tbuf.len() < req_size {
tbuf.resize(req_size, 0);
}
reader.read_exact(&mut tbuf[0..(size_of_rec * na)]).unwrap();
let nbytes = calc_record_bytes(na);
v.offset += nbytes as u64;
v.nbytes += nbytes as u32;
v.nrec += 1;
total_bytes += nbytes;
}
total_bytes += cb_byte_map.len() * header_size as usize;
let mut output_buffer = Cursor::new(vec![0u8; total_bytes]);
let mut next_offset = 0u64;
for (_, v) in cb_byte_map.iter_mut() {
output_buffer.set_position(next_offset);
let cell_bytes = v.nbytes;
let cell_rec = v.nrec;
output_buffer.write_all(&cell_bytes.to_le_bytes()).unwrap();
output_buffer.write_all(&cell_rec.to_le_bytes()).unwrap();
v.offset = output_buffer.position();
let nbytes = v.nbytes as u64;
next_offset += nbytes;
}
reader
.get_mut()
.seek(SeekFrom::Start(0))
.expect("could not get read pointer.");
for _ in 0..(nrec as usize) {
let tup = AtacSeqReadRecord::from_bytes_record_header(reader, bct);
if let Some(v) = cb_byte_map.get_mut(&tup.0) {
output_buffer.set_position(v.offset);
let na = tup.1 as usize;
let nau32 = na as u32;
output_buffer.write_all(&nau32.to_le_bytes()).unwrap();
bct.write_to(tup.0, &mut output_buffer).unwrap();
reader.read_exact(&mut tbuf[0..(size_of_rec * na)]).unwrap();
output_buffer
.write_all(&tbuf[..(size_of_rec * na)])
.unwrap();
v.offset = output_buffer.position();
} else {
panic!("should not have any barcodes we can't find");
}
}
output_buffer.set_position(0);
if compress {
let mut compressed_output =
snap::write::FrameEncoder::new(Cursor::new(Vec::<u8>::with_capacity(total_bytes)));
compressed_output
.write_all(output_buffer.get_ref())
.expect("could not compress the output chunk.");
output_buffer = compressed_output
.into_inner()
.expect("couldn't unwrap the FrameEncoder.");
output_buffer.set_position(0);
}
owriter
.lock()
.unwrap()
.write_all(output_buffer.get_ref())
.unwrap();
cb_byte_map.len()
}
pub struct TempBucket {
pub bucket_id: u32,
pub bucket_writer: Arc<Mutex<BufWriter<File>>>,
pub num_chunks: u32,
pub num_records: u32,
pub num_records_written: AtomicU32,
pub num_bytes_written: AtomicU64,
}
impl TempBucket {
pub fn from_id_and_parent(bucket_id: u32, parent: &std::path::Path) -> Self {
TempBucket {
bucket_id,
bucket_writer: Arc::new(Mutex::new(BufWriter::with_capacity(
4096_usize,
File::create(parent.join(format!("bucket_{}.tmp", bucket_id))).unwrap(),
))),
num_chunks: 0u32,
num_records: 0u32,
num_records_written: AtomicU32::new(0u32),
num_bytes_written: AtomicU64::new(0u64),
}
}
}
#[allow(clippy::too_many_arguments)]
pub fn dump_corrected_cb_chunk_to_temp_file_generic<
B: ConvertiblePrimitiveInteger + std::convert::From<u64>,
T: Read,
R: MappedRecord + KnownSize + CollatableMappedRecord<B>,
>(
reader: &mut BufReader<T>,
rec_context: &<R as MappedRecord>::ParsingContext,
correct_map: &HashMap<u64, u64>,
expected_ori: &Strand,
output_cache: &HashMap<u64, Arc<TempBucket>>,
local_buffers: &mut [Cursor<&mut [u8]>],
flush_limit: usize,
) where
u64: From<B>,
<R as MappedRecord>::ParsingContext: std::fmt::Debug,
{
let mut buf = [0u8; 8];
let mut tbuf = vec![0u8; 4096];
reader.read_exact(&mut buf).unwrap();
let _nbytes = buf.pread::<u32>(0).unwrap();
let nrec = buf.pread::<u32>(4).unwrap();
let expected_ori: MappedFragmentOrientation = expected_ori.into();
for _ in 0..(nrec as usize) {
let mut tup =
<R as CollatableMappedRecord<B>>::from_bytes_collatable_header(reader, rec_context)
.expect("could read header");
if let Some(corrected_id) = correct_map.get(&tup.collate_key().into()) {
let mut rr =
R::from_bytes_with_header_retain_ori(reader, &mut tup, rec_context, &expected_ori);
if rr.is_empty() {
continue;
}
if let Some(v) = output_cache.get(corrected_id) {
let na = tup.naln() as usize;
let nb = R::nbytes(na as u32, rec_context) as u64;
let buffidx = v.bucket_id as usize;
let bcursor = &mut local_buffers[buffidx];
let len = bcursor.position() as usize;
if len + nb as usize >= flush_limit {
let mut filebuf = v.bucket_writer.lock().unwrap();
filebuf.write_all(&bcursor.get_ref()[0..len]).unwrap();
bcursor.set_position(0);
}
rr.set_collate_key((*corrected_id).into());
let blen = bcursor.position() as usize;
rr.write(bcursor, rec_context).expect("can write record");
let alen = bcursor.position() as usize;
let actual = alen - blen;
let expected = R::nbytes(na as u32, rec_context);
assert_eq!(
expected, actual,
"Expected to write {} bytes, but wrote {}.",
expected, actual
);
v.num_records_written.fetch_add(1, Ordering::SeqCst);
v.num_bytes_written.fetch_add(nb, Ordering::SeqCst);
}
} else {
let req_len = R::nbytes_aln(rec_context) * tup.naln() as usize;
let do_resize = req_len > tbuf.len();
if do_resize {
tbuf.resize(req_len, 0);
}
reader.read_exact(&mut tbuf[0..req_len]).unwrap();
if do_resize {
tbuf.resize(4096, 0);
tbuf.shrink_to_fit();
}
}
}
}
#[deprecated(
since = "0.10.0",
note = "This function is highly-specalized and works only with the AlevinFryReadRecordT<u64> type. \
This function has been deprecated in favor of the more generic `dump_corrected_cb_chunk_to_temp_file_generic`. \
Please use that function instead."
)]
#[allow(clippy::too_many_arguments)]
pub fn dump_corrected_cb_chunk_to_temp_file<T: Read>(
reader: &mut BufReader<T>,
bct: &RadIntId,
umit: &RadIntId,
correct_map: &HashMap<u64, u64>,
expected_ori: &Strand,
output_cache: &HashMap<u64, Arc<TempBucket>>,
local_buffers: &mut [Cursor<&mut [u8]>],
flush_limit: usize,
) {
let mut buf = [0u8; 8];
let mut tbuf = vec![0u8; 4096];
reader.read_exact(&mut buf).unwrap();
let _nbytes = buf.pread::<u32>(0).unwrap();
let nrec = buf.pread::<u32>(4).unwrap();
let bc_bytes = bct.bytes_for_type();
let umi_bytes = umit.bytes_for_type();
let na_bytes = std::mem::size_of::<u32>();
let target_id_bytes = std::mem::size_of::<u32>();
for _ in 0..(nrec as usize) {
let tup = AlevinFryReadRecord::from_bytes_record_header(reader, bct, umit);
if let Some(corrected_id) = correct_map.get(&tup.0) {
let rr = AlevinFryReadRecord::from_bytes_with_header_keep_ori(
reader,
tup.0,
tup.1,
tup.2,
expected_ori,
);
if rr.is_empty() {
continue;
}
if let Some(v) = output_cache.get(corrected_id) {
let nb = (rr.refs.len() * target_id_bytes + na_bytes + bc_bytes + umi_bytes) as u64;
let buffidx = v.bucket_id as usize;
let bcursor = &mut local_buffers[buffidx];
let len = bcursor.position() as usize;
if len + nb as usize >= flush_limit {
let mut filebuf = v.bucket_writer.lock().unwrap();
filebuf.write_all(&bcursor.get_ref()[0..len]).unwrap();
bcursor.set_position(0);
}
let na = rr.refs.len() as u32;
bcursor.write_all(&na.to_le_bytes()).unwrap();
bct.write_to(*corrected_id, bcursor).unwrap();
umit.write_to(rr.umi, bcursor).unwrap();
bcursor.write_all(as_u8_slice(&rr.refs[..])).unwrap();
v.num_records_written.fetch_add(1, Ordering::SeqCst);
v.num_bytes_written.fetch_add(nb, Ordering::SeqCst);
}
} else {
let req_len = target_id_bytes * (tup.2 as usize);
let do_resize = req_len > tbuf.len();
if do_resize {
tbuf.resize(req_len, 0);
}
reader
.read_exact(&mut tbuf[0..(target_id_bytes * (tup.2 as usize))])
.unwrap();
if do_resize {
tbuf.resize(4096, 0);
tbuf.shrink_to_fit();
}
}
}
}
#[allow(clippy::too_many_arguments)]
pub fn dump_corrected_cb_chunk_to_temp_file_atac<T: Read>(
reader: &mut BufReader<T>,
bct: &RadIntId,
correct_map: &HashMap<u64, u64>,
output_cache: &HashMap<u64, Arc<TempBucket>>,
local_buffers: &mut [Cursor<&mut [u8]>],
flush_limit: usize,
ck: CollateKey, )
where
{
let mut buf = [0u8; 8];
let mut tbuf = vec![0u8; 4096];
reader.read_exact(&mut buf).unwrap();
let _nbytes = buf.pread::<u32>(0).unwrap();
let nrec = buf.pread::<u32>(4).unwrap();
let bc_bytes = bct.bytes_for_type();
let na_bytes = std::mem::size_of::<u32>();
let target_id_bytes = std::mem::size_of::<u32>() + std::mem::size_of::<u8>() + std::mem::size_of::<u32>() + std::mem::size_of::<u16>();
for _ in 0..(nrec as usize) {
let tup = AtacSeqReadRecord::from_bytes_record_header(reader, bct);
if let Some(corrected_id) = correct_map.get(&tup.0) {
let rr = AtacSeqReadRecord::from_bytes_with_header(reader, tup.0, tup.1);
if rr.is_empty() || tup.1 > 1 {
continue;
}
let pos = rr.start_pos[0];
let ref_id = rr.refs[0];
let check_id = match ck {
CollateKey::Barcode => *corrected_id,
CollateKey::Pos(ref f) => f(pos, ref_id as usize) as u64,
};
if let Some(v) = output_cache.get(&check_id) {
let nb = (rr.refs.len() * target_id_bytes + na_bytes + bc_bytes) as u64;
let buffidx = v.bucket_id as usize;
let bcursor = &mut local_buffers[buffidx];
let len = bcursor.position() as usize;
if len + nb as usize >= flush_limit {
let mut filebuf = v.bucket_writer.lock().unwrap();
filebuf.write_all(&bcursor.get_ref()[0..len]).unwrap();
bcursor.set_position(0);
}
let na = rr.refs.len() as u32;
bcursor.write_all(&na.to_le_bytes()).unwrap();
bct.write_to(*corrected_id, bcursor).unwrap();
bcursor.write_all(as_u8_slice(&rr.refs[..])).unwrap();
bcursor.write_all(as_u8_slice_u8(&rr.map_type[..])).unwrap();
bcursor.write_all(as_u8_slice(&rr.start_pos[..])).unwrap();
bcursor
.write_all(as_u8_slice_u16(&rr.frag_lengths[..]))
.unwrap();
v.num_records_written.fetch_add(1, Ordering::SeqCst);
v.num_bytes_written.fetch_add(nb, Ordering::SeqCst);
}
} else {
let req_len = target_id_bytes * (tup.1 as usize);
let do_resize = req_len > tbuf.len();
if do_resize {
tbuf.resize(req_len, 0);
}
reader
.read_exact(&mut tbuf[0..(target_id_bytes * (tup.1 as usize))])
.unwrap();
if do_resize {
tbuf.resize(4096, 0);
tbuf.shrink_to_fit();
}
}
}
}
pub fn as_u8_slice(v: &[u32]) -> &[u8] {
unsafe { std::slice::from_raw_parts(v.as_ptr() as *const u8, std::mem::size_of_val(v)) }
}
pub fn as_u8_slice_u16(v: &[u16]) -> &[u8] {
unsafe { std::slice::from_raw_parts(v.as_ptr() as *const u8, std::mem::size_of_val(v)) }
}
pub fn as_u8_slice_u8(v: &[u8]) -> &[u8] {
unsafe { std::slice::from_raw_parts(v.as_ptr(), std::mem::size_of_val(v)) }
}
#[cfg(test)]
mod tests {
use crate::BarcodeLookupMap;
#[test]
fn test_barcode_lookup_map() {
let barcode_sv_even = vec![
b"AACC", b"AAGG", b"CAGT", b"CATT", b"GACC", b"GATA", b"TCAG", b"TCGT",
];
let barcode_sv_odd = vec![
b"AACCA", b"AAGGC", b"CAGTA", b"CATTG", b"GACCG", b"GATAC", b"TCAGA", b"TCGTG",
];
let mut barcode_even = Vec::with_capacity(barcode_sv_even.len());
for b in barcode_sv_even.clone() {
if let Some((_, km, _)) =
needletail::bitkmer::BitNuclKmer::new(&b[..], b.len() as u8, false).next()
{
barcode_even.push(km.0);
}
}
let mut barcode_odd = Vec::with_capacity(barcode_sv_odd.len());
for b in barcode_sv_odd.clone() {
if let Some((_, km, _)) =
needletail::bitkmer::BitNuclKmer::new(&b[..], b.len() as u8, false).next()
{
barcode_odd.push(km.0);
}
}
let me = BarcodeLookupMap::new(barcode_even, 4);
let mo = BarcodeLookupMap::new(barcode_odd, 5);
let x = b"CAGA";
if let Some((_, et, _)) =
needletail::bitkmer::BitNuclKmer::new(&x[..], x.len() as u8, false).next()
{
assert_eq!((Some(2), 1), me.find_neighbors(et.0, false));
}
let x = b"CAATG";
if let Some((_, et, _)) =
needletail::bitkmer::BitNuclKmer::new(&x[..], x.len() as u8, false).next()
{
assert_eq!((Some(3), 1), mo.find_neighbors(et.0, false));
}
}
}