Skip to content
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

feat: skip too short reads #32

Merged
merged 3 commits into from
Jul 28, 2023
Merged

feat: skip too short reads #32

Changes from 1 commit
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
6fe2ce9
feat: add an option to skip too short reads
nh13 Jun 26, 2023
884c449
nit
nh13 Jul 13, 2023
92279c6
fixups
nh13 Jul 28, 2023
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
Prev Previous commit
fixups
  • Loading branch information
@nh13
nh13 committed Jul 28, 2023
commit 92279c6ab0eeb31d54da06eebf13bb7b40cd582d
73 changes: 59 additions & 14 deletions src/bin/commands/demux.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -273,7 +273,7 @@ struct ReadSetIterator {
/// The file containing FASTQ records.
source: FastqReader<Box<dyn BufRead + Send>>,
/// Valid reasons for skipping reads, otherwise panic!
skip_reasons: HashSet<SkipReason>,
skip_reasons: Vec<SkipReason>,
}

impl Iterator for ReadSetIterator {
Expand All @@ -286,18 +286,28 @@ impl Iterator for ReadSetIterator {
let read_name = next_fq_rec.head();
let bases = next_fq_rec.seq();
let quals = next_fq_rec.qual();
for (read_segment_index, read_segment) in self.read_structure.iter().enumerate() {
// Check that we have enough bases for this segment. For variable length segments,
// we must have at least one base.
if self.skip_reasons.contains(&SkipReason::TooFewBases)
&& bases.len() < read_segment.length().unwrap_or(1)
{

// Check that we have enough bases for this segment. For variable length segments,
// we must have at least one base.
let min_len = self.read_structure.iter().map(|s| s.length().unwrap_or(1)).sum();
if bases.len() < min_len {
if self.skip_reasons.contains(&SkipReason::TooFewBases) {
return Some(ReadSet {
header: read_name.to_vec(),
segments: vec![],
skip_reason: Some(SkipReason::TooFewBases),
});
}
panic!(
"Read {} had too few bases to demux {} vs. {} needed in read structure {}.",
String::from_utf8(read_name.to_vec()).unwrap(),
bases.len(),
min_len,
self.read_structure
);
}

for (read_segment_index, read_segment) in self.read_structure.iter().enumerate() {
// Extract the bases and qualities for this segment
let (seq, quals) =
read_segment.extract_bases_and_quals(bases, quals).unwrap_or_else(|e| {
Expand Down Expand Up @@ -332,7 +342,7 @@ impl ReadSetIterator {
pub fn new(
read_structure: ReadStructure,
source: FastqReader<Box<dyn BufRead + Send>>,
skip_reasons: HashSet<SkipReason>,
skip_reasons: Vec<SkipReason>,
) -> Self {
Self { read_structure, source, skip_reasons }
}
Expand Down Expand Up @@ -878,12 +888,8 @@ impl Command for Demux {
let mut fq_iterators = fq_sources
.zip(self.read_structures.clone().into_iter())
.map(|(source, read_structure)| {
ReadSetIterator::new(
read_structure,
source,
self.skip_reasons.iter().copied().collect::<HashSet<_>>(),
)
.read_ahead(1000, 1000)
ReadSetIterator::new(read_structure, source, self.skip_reasons.clone())
.read_ahead(1000, 1000)
})
.collect::<Vec<_>>();

Expand Down Expand Up @@ -1735,6 +1741,45 @@ mod tests {
demux_inputs.execute().unwrap();
}

#[test]
#[should_panic(
expected = "Read ex_2 had too few bases to demux 0 vs. 1 needed in read structure +T."
)]
fn test_fails_if_reads_too_short() {
let tmp = TempDir::new().unwrap();
let read_structures =
vec![ReadStructure::from_str("+T").unwrap(), ReadStructure::from_str("7B").unwrap()];

let records = vec![
vec!["AAAAAAA", &SAMPLE1_BARCODE[0..7]], // barcode too short
vec!["CCCCCCC", SAMPLE1_BARCODE], // barcode the correct length
vec!["", SAMPLE1_BARCODE], // template basese too short
vec!["G", SAMPLE1_BARCODE], // barcode the correct length
];

let input_files = vec![
fastq_file(&tmp, "read1", "ex", &[records[0][0], records[1][0], records[2][0]]),
fastq_file(&tmp, "index1", "ex", &[records[0][1], records[1][1], records[2][1]]),
];
let sample_metadata = metadata(&tmp);
let output_dir: PathBuf = tmp.path().to_path_buf().join("output");

let demux_inputs = Demux {
inputs: input_files,
read_structures,
sample_metadata,
output_types: vec!['T', 'B'],
output: output_dir.clone(),
unmatched_prefix: "unmatched".to_owned(),
max_mismatches: 1,
min_mismatch_delta: 2,
threads: 5,
compression_level: 5,
skip_reasons: vec![],
};
demux_inputs.execute().unwrap();
}

#[test]
fn test_skip_reads_too_short() {
let tmp = TempDir::new().unwrap();
Expand Down