From e1473a6d4108143b43560986ba74d67ea9bdaa88 Mon Sep 17 00:00:00 2001 From: Andrew Frantz Date: Sun, 11 Feb 2024 16:53:47 -0500 Subject: [PATCH] tests(derive): all derive commands have tests --- src/derive/encoding/compute.rs | 247 +++++++++++++++++ src/derive/endedness/compute.rs | 49 +++- src/derive/junction_annotation/compute.rs | 319 +++++++++++----------- src/derive/junction_annotation/results.rs | 2 +- src/utils/alignment.rs | 16 +- 5 files changed, 447 insertions(+), 186 deletions(-) diff --git a/src/derive/encoding/compute.rs b/src/derive/encoding/compute.rs index eb60962..421ce79 100644 --- a/src/derive/encoding/compute.rs +++ b/src/derive/encoding/compute.rs @@ -77,3 +77,250 @@ pub fn predict(score_set: HashSet) -> Result = HashSet::new(); + score_set.insert(40); + score_set.insert(41); + score_set.insert(42); + score_set.insert(43); + score_set.insert(44); + score_set.insert(45); + score_set.insert(46); + score_set.insert(47); + score_set.insert(48); + score_set.insert(49); + score_set.insert(50); + score_set.insert(51); + score_set.insert(52); + score_set.insert(53); + score_set.insert(54); + score_set.insert(55); + score_set.insert(56); + score_set.insert(57); + score_set.insert(58); + score_set.insert(59); + score_set.insert(60); + score_set.insert(61); + score_set.insert(62); + score_set.insert(63); + score_set.insert(64); + score_set.insert(65); + score_set.insert(66); + score_set.insert(67); + score_set.insert(68); + score_set.insert(69); + score_set.insert(70); + score_set.insert(71); + score_set.insert(72); + score_set.insert(73); + score_set.insert(74); + score_set.insert(75); + score_set.insert(76); + score_set.insert(77); + score_set.insert(78); + score_set.insert(79); + score_set.insert(80); + score_set.insert(81); + score_set.insert(82); + score_set.insert(83); + score_set.insert(84); + score_set.insert(85); + score_set.insert(86); + score_set.insert(87); + score_set.insert(88); + score_set.insert(89); + score_set.insert(90); + score_set.insert(91); + score_set.insert(92); + score_set.insert(93); + + let result = predict(score_set).unwrap(); + assert!(result.succeeded); + assert_eq!(result.encoding, Some("Illumina 1.3".to_string())); + assert_eq!(result.observed_min, 40); + assert_eq!(result.observed_max, 93); + } + + #[test] + fn test_predict_illumina_1_0() { + let mut score_set: HashSet = HashSet::new(); + score_set.insert(26); + score_set.insert(27); + score_set.insert(28); + score_set.insert(29); + score_set.insert(30); + score_set.insert(31); + score_set.insert(32); + score_set.insert(33); + score_set.insert(34); + score_set.insert(35); + score_set.insert(36); + score_set.insert(37); + score_set.insert(38); + score_set.insert(39); + score_set.insert(40); + score_set.insert(41); + score_set.insert(42); + score_set.insert(43); + score_set.insert(44); + score_set.insert(45); + score_set.insert(46); + score_set.insert(47); + score_set.insert(48); + score_set.insert(49); + score_set.insert(50); + score_set.insert(51); + score_set.insert(52); + score_set.insert(53); + score_set.insert(54); + score_set.insert(55); + score_set.insert(56); + score_set.insert(57); + score_set.insert(58); + score_set.insert(59); + score_set.insert(60); + score_set.insert(61); + score_set.insert(62); + score_set.insert(63); + score_set.insert(64); + score_set.insert(65); + score_set.insert(66); + score_set.insert(67); + score_set.insert(68); + score_set.insert(69); + score_set.insert(70); + score_set.insert(71); + score_set.insert(72); + score_set.insert(73); + score_set.insert(74); + score_set.insert(75); + score_set.insert(76); + score_set.insert(77); + score_set.insert(78); + score_set.insert(79); + score_set.insert(80); + score_set.insert(81); + score_set.insert(82); + score_set.insert(83); + score_set.insert(84); + score_set.insert(85); + score_set.insert(86); + score_set.insert(87); + score_set.insert(88); + score_set.insert(89); + score_set.insert(90); + score_set.insert(91); + score_set.insert(92); + score_set.insert(93); + + let result = predict(score_set).unwrap(); + assert!(result.succeeded); + assert_eq!(result.encoding, Some("Illumina 1.0".to_string())); + assert_eq!(result.observed_min, 26); + assert_eq!(result.observed_max, 93); + } + + #[test] + fn test_predict_sanger() { + let mut score_set: HashSet = HashSet::new(); + score_set.insert(0); + score_set.insert(1); + score_set.insert(2); + score_set.insert(3); + score_set.insert(4); + score_set.insert(5); + score_set.insert(6); + score_set.insert(7); + score_set.insert(8); + score_set.insert(9); + score_set.insert(10); + score_set.insert(11); + score_set.insert(12); + score_set.insert(13); + score_set.insert(14); + score_set.insert(15); + score_set.insert(16); + score_set.insert(17); + score_set.insert(18); + score_set.insert(19); + score_set.insert(20); + score_set.insert(21); + score_set.insert(22); + score_set.insert(23); + score_set.insert(24); + score_set.insert(25); + score_set.insert(26); + score_set.insert(27); + score_set.insert(28); + score_set.insert(29); + score_set.insert(30); + score_set.insert(31); + score_set.insert(32); + score_set.insert(33); + score_set.insert(34); + score_set.insert(35); + score_set.insert(36); + score_set.insert(37); + score_set.insert(38); + score_set.insert(39); + score_set.insert(40); + score_set.insert(41); + score_set.insert(42); + score_set.insert(43); + score_set.insert(44); + score_set.insert(45); + score_set.insert(46); + score_set.insert(47); + score_set.insert(48); + score_set.insert(49); + score_set.insert(50); + score_set.insert(51); + score_set.insert(52); + score_set.insert(53); + score_set.insert(54); + score_set.insert(55); + score_set.insert(56); + score_set.insert(57); + score_set.insert(58); + score_set.insert(59); + score_set.insert(60); + score_set.insert(61); + score_set.insert(62); + score_set.insert(63); + score_set.insert(64); + score_set.insert(65); + score_set.insert(66); + score_set.insert(67); + score_set.insert(68); + + let result = predict(score_set).unwrap(); + assert!(result.succeeded); + assert_eq!(result.encoding, Some("Sanger/Illumina 1.8".to_string())); + assert_eq!(result.observed_min, 0); + assert_eq!(result.observed_max, 68); + } + + #[test] + fn test_predict_fail() { + let score_set: HashSet = HashSet::new(); + let result = predict(score_set); + assert!(result.is_err()); + } + + #[test] + fn test_predict_too_high_max_score() { + let mut score_set: HashSet = HashSet::new(); + score_set.insert(94); + let result = predict(score_set).unwrap(); + assert!(!result.succeeded); + assert_eq!(result.encoding, None); + assert_eq!(result.observed_min, 94); + assert_eq!(result.observed_max, 94); + } +} diff --git a/src/derive/endedness/compute.rs b/src/derive/endedness/compute.rs index 265fc2d..4179d2e 100644 --- a/src/derive/endedness/compute.rs +++ b/src/derive/endedness/compute.rs @@ -307,9 +307,8 @@ pub fn predict( mod tests { use super::*; - // TODO add tests for unsegmented reads #[test] - fn test_predict_endedness() { + fn test_predict_endedness_from_first_and_last() { let mut ordering_flags: HashMap = HashMap::new(); ordering_flags.insert( Arc::new("overall".to_string()), @@ -340,7 +339,38 @@ mod tests { } #[test] - fn test_derive_endedness_from_all_zero_counts() { + fn test_predict_endedness_from_unsegmented() { + let mut ordering_flags: HashMap = HashMap::new(); + ordering_flags.insert( + Arc::new("overall".to_string()), + OrderingFlagsCounts { + unsegmented: 1, + first: 0, + last: 0, + both: 0, + neither: 0, + }, + ); + let result = predict_endedness( + "overall".to_string(), + ordering_flags + .get(&Arc::new("overall".to_string())) + .unwrap(), + 0.0, + None, + false, + ); + assert!(result.succeeded); + assert_eq!(result.endedness, "Single-End"); + assert_eq!(result.first, 0); + assert_eq!(result.last, 0); + assert_eq!(result.both, 0); + assert_eq!(result.neither, 0); + assert_eq!(result.rpt, None); + } + + #[test] + fn test_predict_endedness_from_all_zero_counts() { let mut ordering_flags: HashMap = HashMap::new(); ordering_flags.insert(Arc::new(String::from("rg1")), OrderingFlagsCounts::new()); let result = predict_endedness( @@ -360,7 +390,7 @@ mod tests { } #[test] - fn test_derive_endedness_from_only_first() { + fn test_predict_from_only_first() { let mut ordering_flags: HashMap = HashMap::new(); ordering_flags.insert( Arc::new("overall".to_string()), @@ -384,7 +414,7 @@ mod tests { } #[test] - fn test_derive_endedness_from_only_last() { + fn test_predict_from_only_last() { let mut ordering_flags: HashMap = HashMap::new(); ordering_flags.insert( Arc::new("overall".to_string()), @@ -408,7 +438,7 @@ mod tests { } #[test] - fn test_derive_endedness_from_only_both() { + fn test_predict_from_only_both() { let mut ordering_flags: HashMap = HashMap::new(); ordering_flags.insert( Arc::new("overall".to_string()), @@ -432,7 +462,7 @@ mod tests { } #[test] - fn test_derive_endedness_from_only_neither() { + fn test_predict_from_only_neither() { let mut ordering_flags: HashMap = HashMap::new(); ordering_flags.insert( Arc::new("overall".to_string()), @@ -456,7 +486,7 @@ mod tests { } #[test] - fn test_derive_endedness_from_first_and_last() { + fn test_predict_from_first_and_last() { let mut ordering_flags: HashMap = HashMap::new(); ordering_flags.insert( Arc::new("overall".to_string()), @@ -518,7 +548,7 @@ mod tests { } #[test] - fn test_derive_endedness_from_first_and_last_with_rpt() { + fn test_predict_with_rpt_complex() { let mut ordering_flags: HashMap = HashMap::new(); let rg_paired = Arc::new("rg_paired".to_string()); let rg_single = Arc::new("rg_single".to_string()); @@ -580,7 +610,6 @@ mod tests { assert_eq!(result.read_groups.len(), 2); // We can't know which read group will be first in the vector. // But both should succeed. - print!("{:?}", result.read_groups); assert!(result.read_groups[0].succeeded && result.read_groups[1].succeeded); } } diff --git a/src/derive/junction_annotation/compute.rs b/src/derive/junction_annotation/compute.rs index 069a5dc..428daac 100644 --- a/src/derive/junction_annotation/compute.rs +++ b/src/derive/junction_annotation/compute.rs @@ -598,172 +598,157 @@ mod tests { assert_eq!(results.records.not_spliced, 0); assert_eq!(results.records.bad_mapq, 0); - // TODO: Below tests are not working as expected. Need to fix them. - // Test complete novel junction - // let mut record = Record::default(); - // let r6_name: ReadName = "complete1".parse().unwrap(); - // *record.read_name_mut() = Some(r6_name); - // *record.reference_sequence_id_mut() = Some(0); - // *record.alignment_start_mut() = Position::new(1); - // *record.cigar_mut() = "10M10N10M10N10M".parse().unwrap(); - // *record.mapping_quality_mut() = MappingQuality::new(60); - // record.flags_mut().set(0x4.into(), false); - // process(&record, &exons, &header, ¶ms, &mut results).unwrap(); - // assert_eq!(results.records.processed, 6); - // assert_eq!(results.records.filtered_by_flags, 2); - // assert_eq!(results.records.not_spliced, 0); - // assert_eq!(results.records.low_mapq, 0); - // assert_eq!(results.records.missing_mapq, 0); - - // // Test complete novel junction (again for more read support) - // let mut record = Record::default(); - // let r6_name: ReadName = "complete2".parse().unwrap(); - // *record.read_name_mut() = Some(r6_name); - // *record.reference_sequence_id_mut() = Some(0); - // *record.alignment_start_mut() = Position::new(1); - // *record.cigar_mut() = "10M10N10M10N10M".parse().unwrap(); - // *record.mapping_quality_mut() = MappingQuality::new(60); - // record.flags_mut().set(0x4.into(), false); - // process(&record, &exons, &header, ¶ms, &mut results).unwrap(); - // assert_eq!(results.records.processed, 7); - // assert_eq!(results.records.filtered_by_flags, 2); - // assert_eq!(results.records.not_spliced, 0); - // assert_eq!(results.records.low_mapq, 0); - // assert_eq!(results.records.missing_mapq, 0); - - // // Test fails MAPQ filter - // let mut record = Record::default(); - // let r7_name: ReadName = "low_mapq".parse().unwrap(); - // *record.read_name_mut() = Some(r7_name); - // *record.reference_sequence_id_mut() = Some(0); - // *record.alignment_start_mut() = Position::new(1); - // *record.cigar_mut() = "10M10N10M".parse().unwrap(); - // *record.mapping_quality_mut() = MappingQuality::new(20); - // record.flags_mut().set(0x4.into(), false); - // process(&record, &exons, &header, ¶ms, &mut results).unwrap(); - // assert_eq!(results.records.processed, 6); - // assert_eq!(results.records.filtered_by_flags, 2); - // assert_eq!(results.records.not_spliced, 0); - // assert_eq!(results.records.low_mapq, 1); - // assert_eq!(results.records.missing_mapq, 0); - - // // Test missing MAPQ - // let mut record = Record::default(); - // let r8_name: ReadName = "missing_mapq".parse().unwrap(); - // *record.read_name_mut() = Some(r8_name); - // *record.reference_sequence_id_mut() = Some(0); - // *record.alignment_start_mut() = Position::new(1); - // *record.cigar_mut() = "10M10N10M".parse().unwrap(); - // *record.mapping_quality_mut() = MappingQuality::new(255); - // record.flags_mut().set(0x4.into(), false); - // process(&record, &exons, &header, ¶ms, &mut results).unwrap(); - // assert_eq!(results.records.processed, 6); - // assert_eq!(results.records.filtered_by_flags, 2); - // assert_eq!(results.records.not_spliced, 0); - // assert_eq!(results.records.low_mapq, 1); - // assert_eq!(results.records.missing_mapq, 1); - - // // Test that intron is too short - // let mut record = Record::default(); - // let r9_name: ReadName = "short".parse().unwrap(); - // *record.read_name_mut() = Some(r9_name); - // *record.reference_sequence_id_mut() = Some(0); - // *record.alignment_start_mut() = Position::new(1); - // *record.cigar_mut() = "5M5N5M".parse().unwrap(); - // *record.mapping_quality_mut() = MappingQuality::new(60); - // record.flags_mut().set(0x4.into(), false); - // process(&record, &exons, &header, ¶ms, &mut results).unwrap(); - // assert_eq!(results.records.processed, 7); // Still gets processed, will be filtered later - // assert_eq!(results.records.filtered_by_flags, 2); - // assert_eq!(results.records.not_spliced, 0); - // assert_eq!(results.records.low_mapq, 1); - // assert_eq!(results.records.missing_mapq, 1); - - // // Test that that reads not spliced are ignored - // let mut record = Record::default(); - // let r10_name: ReadName = "not_spliced".parse().unwrap(); - // *record.read_name_mut() = Some(r10_name); - // *record.reference_sequence_id_mut() = Some(0); - // *record.alignment_start_mut() = Position::new(1); - // *record.cigar_mut() = "10M".parse().unwrap(); - // *record.mapping_quality_mut() = MappingQuality::new(60); - // record.flags_mut().set(0x4.into(), false); - // process(&record, &exons, &header, ¶ms, &mut results).unwrap(); - // assert_eq!(results.records.processed, 7); - // assert_eq!(results.records.filtered_by_flags, 2); - // assert_eq!(results.records.not_spliced, 1); - // assert_eq!(results.records.low_mapq, 1); - // assert_eq!(results.records.missing_mapq, 1); - - // // Test unannoted reference - // let mut record = Record::default(); - // let r11_name: ReadName = "unannotated1".parse().unwrap(); - // *record.read_name_mut() = Some(r11_name); - // *record.reference_sequence_id_mut() = Some(1); - // *record.alignment_start_mut() = Position::new(1); - // *record.cigar_mut() = "10M10N10M".parse().unwrap(); - // *record.mapping_quality_mut() = MappingQuality::new(60); - // record.flags_mut().set(0x4.into(), false); - // process(&record, &exons, &header, ¶ms, &mut results).unwrap(); - // assert_eq!(results.records.processed, 8); - // assert_eq!(results.records.filtered_by_flags, 2); - // assert_eq!(results.records.not_spliced, 1); - // assert_eq!(results.records.low_mapq, 1); - // assert_eq!(results.records.missing_mapq, 1); - - // // Test unannoted reference (again for more read support) - // let mut record = Record::default(); - // let r11_name: ReadName = "unannotated2".parse().unwrap(); - // *record.read_name_mut() = Some(r11_name); - // *record.reference_sequence_id_mut() = Some(1); - // *record.alignment_start_mut() = Position::new(1); - // *record.cigar_mut() = "10M10N10M".parse().unwrap(); - // *record.mapping_quality_mut() = MappingQuality::new(60); - // record.flags_mut().set(0x4.into(), false); - // process(&record, &exons, &header, ¶ms, &mut results).unwrap(); - // assert_eq!(results.records.processed, 9); - // assert_eq!(results.records.filtered_by_flags, 2); - // assert_eq!(results.records.not_spliced, 1); - // assert_eq!(results.records.low_mapq, 1); - // assert_eq!(results.records.missing_mapq, 1); - - // // Test summarize - // summarize(&mut results, ¶ms); - - // assert_eq!(results.summary.total_rejected_junctions, 1); - // assert_eq!(results.summary.intron_too_short, 1); - // assert_eq!(results.summary.junctions_with_not_enough_read_support, 1); - // assert_eq!(results.summary.known_junctions, 1); - // assert_eq!(results.summary.known_junctions_read_support, 2); - // assert_eq!(results.summary.partial_novel_junctions, 1); - // assert_eq!(results.summary.partial_novel_junctions_read_support, 2); - // assert_eq!(results.summary.complete_novel_junctions, 1); - // assert_eq!(results.summary.complete_novel_junctions_read_support, 2); - // assert_eq!(results.summary.unannotated_reference_junctions, 1); - // assert_eq!( - // results.summary.unannotated_reference_junctions_read_support, - // 2 - // ); - // assert_eq!(results.summary.total_junctions, 4); - // assert_eq!(results.summary.total_junctions_read_support, 8); - // assert_eq!(results.summary.known_junctions_percent, 33.33333333333333); - // assert_eq!( - // results.summary.partial_novel_junctions_percent, - // 33.33333333333333 - // ); - // assert_eq!( - // results.summary.complete_novel_junctions_percent, - // 33.33333333333333 - // ); - // assert_eq!(results.summary.average_junction_read_support, 2.0); - // assert_eq!(results.summary.average_known_junction_read_support, 2.0); - // assert_eq!( - // results.summary.average_partial_novel_junction_read_support, - // 2.0 - // ); - // assert_eq!( - // results.summary.average_complete_novel_junction_read_support, - // 2.0 - // ); + // Test complete novel junction with 2 junctions + let mut record = Record::default(); + let r6_name: ReadName = "complete_twice1".parse().unwrap(); + *record.read_name_mut() = Some(r6_name); + *record.reference_sequence_id_mut() = Some(0); + *record.alignment_start_mut() = Position::new(200); + *record.cigar_mut() = "10M10N10M10N10M".parse().unwrap(); + *record.mapping_quality_mut() = MappingQuality::new(60); + record.flags_mut().set(0x4.into(), false); + process(&record, &exons, &header, ¶ms, &mut results).unwrap(); + assert_eq!(results.records.processed, 5); + assert_eq!(results.records.filtered_by_flags, 2); + assert_eq!(results.records.not_spliced, 0); + assert_eq!(results.records.bad_mapq, 0); + + // Test complete novel junction with 2 junctions (again for more read support) + let mut record = Record::default(); + let r6_name: ReadName = "complete_twice2".parse().unwrap(); + *record.read_name_mut() = Some(r6_name); + *record.reference_sequence_id_mut() = Some(0); + *record.alignment_start_mut() = Position::new(200); + *record.cigar_mut() = "10M10N10M10N10M".parse().unwrap(); + *record.mapping_quality_mut() = MappingQuality::new(60); + record.flags_mut().set(0x4.into(), false); + process(&record, &exons, &header, ¶ms, &mut results).unwrap(); + assert_eq!(results.records.processed, 6); + assert_eq!(results.records.filtered_by_flags, 2); + assert_eq!(results.records.not_spliced, 0); + assert_eq!(results.records.bad_mapq, 0); + + // Test fails MAPQ filter + let mut record = Record::default(); + let r7_name: ReadName = "low_mapq".parse().unwrap(); + *record.read_name_mut() = Some(r7_name); + *record.reference_sequence_id_mut() = Some(0); + *record.alignment_start_mut() = Position::new(1); + *record.cigar_mut() = "10M10N10M".parse().unwrap(); + *record.mapping_quality_mut() = MappingQuality::new(20); + record.flags_mut().set(0x4.into(), false); + process(&record, &exons, &header, ¶ms, &mut results).unwrap(); + assert_eq!(results.records.processed, 6); + assert_eq!(results.records.filtered_by_flags, 2); + assert_eq!(results.records.not_spliced, 0); + assert_eq!(results.records.bad_mapq, 1); + + // Test missing MAPQ + let mut record = Record::default(); + let r8_name: ReadName = "bad_mapq".parse().unwrap(); + *record.read_name_mut() = Some(r8_name); + *record.reference_sequence_id_mut() = Some(0); + *record.alignment_start_mut() = Position::new(1); + *record.cigar_mut() = "10M10N10M".parse().unwrap(); + *record.mapping_quality_mut() = None; + record.flags_mut().set(0x4.into(), false); + process(&record, &exons, &header, ¶ms, &mut results).unwrap(); + assert_eq!(results.records.processed, 6); + assert_eq!(results.records.filtered_by_flags, 2); + assert_eq!(results.records.not_spliced, 0); + assert_eq!(results.records.bad_mapq, 2); + + // Test that intron is too short + let mut record = Record::default(); + let r9_name: ReadName = "short".parse().unwrap(); + *record.read_name_mut() = Some(r9_name); + *record.reference_sequence_id_mut() = Some(0); + *record.alignment_start_mut() = Position::new(1); + *record.cigar_mut() = "5M5N5M".parse().unwrap(); + *record.mapping_quality_mut() = MappingQuality::new(60); + record.flags_mut().set(0x4.into(), false); + process(&record, &exons, &header, ¶ms, &mut results).unwrap(); + assert_eq!(results.records.processed, 7); // Still gets processed, will be filtered later + assert_eq!(results.records.filtered_by_flags, 2); + assert_eq!(results.records.not_spliced, 0); + assert_eq!(results.records.bad_mapq, 2); + + // Test that that reads not spliced are ignored + let mut record = Record::default(); + let r10_name: ReadName = "not_spliced".parse().unwrap(); + *record.read_name_mut() = Some(r10_name); + *record.reference_sequence_id_mut() = Some(0); + *record.alignment_start_mut() = Position::new(1); + *record.cigar_mut() = "10M".parse().unwrap(); + *record.mapping_quality_mut() = MappingQuality::new(60); + record.flags_mut().set(0x4.into(), false); + process(&record, &exons, &header, ¶ms, &mut results).unwrap(); + assert_eq!(results.records.processed, 7); + assert_eq!(results.records.filtered_by_flags, 2); + assert_eq!(results.records.not_spliced, 1); + assert_eq!(results.records.bad_mapq, 2); + + // Test unannoted reference + let mut record = Record::default(); + let r11_name: ReadName = "unannotated1".parse().unwrap(); + *record.read_name_mut() = Some(r11_name); + *record.reference_sequence_id_mut() = Some(1); + *record.alignment_start_mut() = Position::new(1); + *record.cigar_mut() = "10M10N10M".parse().unwrap(); + *record.mapping_quality_mut() = MappingQuality::new(60); + record.flags_mut().set(0x4.into(), false); + process(&record, &exons, &header, ¶ms, &mut results).unwrap(); + assert_eq!(results.records.processed, 8); + assert_eq!(results.records.filtered_by_flags, 2); + assert_eq!(results.records.not_spliced, 1); + assert_eq!(results.records.bad_mapq, 2); + + // Test unannoted reference (again for more read support) + let mut record = Record::default(); + let r11_name: ReadName = "unannotated2".parse().unwrap(); + *record.read_name_mut() = Some(r11_name); + *record.reference_sequence_id_mut() = Some(1); + *record.alignment_start_mut() = Position::new(1); + *record.cigar_mut() = "10M10N10M".parse().unwrap(); + *record.mapping_quality_mut() = MappingQuality::new(60); + record.flags_mut().set(0x4.into(), false); + process(&record, &exons, &header, ¶ms, &mut results).unwrap(); + assert_eq!(results.records.processed, 9); + assert_eq!(results.records.filtered_by_flags, 2); + assert_eq!(results.records.not_spliced, 1); + assert_eq!(results.records.bad_mapq, 2); + + // Test summarize + summarize(&mut results, ¶ms); + + assert_eq!(results.summary.total_rejected_junctions, 1); + assert_eq!(results.summary.intron_too_short, 1); + assert_eq!(results.summary.junctions_with_not_enough_read_support, 1); + assert_eq!(results.summary.known_junctions, 1); + assert_eq!(results.summary.known_junctions_read_support, 2); + assert_eq!(results.summary.partial_novel_junctions, 1); + assert_eq!(results.summary.partial_novel_junctions_read_support, 2); + assert_eq!(results.summary.complete_novel_junctions, 2); + assert_eq!(results.summary.complete_novel_junctions_read_support, 4); + assert_eq!(results.summary.unannotated_reference_junctions, 1); + assert_eq!( + results.summary.unannotated_reference_junctions_read_support, + 2 + ); + assert_eq!(results.summary.total_junctions, 5); + assert_eq!(results.summary.total_junctions_read_support, 10); + assert_eq!(results.summary.known_junctions_percent, 25.0); + assert_eq!(results.summary.partial_novel_junctions_percent, 25.0); + assert_eq!(results.summary.complete_novel_junctions_percent, 50.0); + assert_eq!(results.summary.average_junction_read_support, 2.0); + assert_eq!(results.summary.average_known_junction_read_support, 2.0); + assert_eq!( + results.summary.average_partial_novel_junction_read_support, + 2.0 + ); + assert_eq!( + results.summary.average_complete_novel_junction_read_support, + 2.0 + ); } } diff --git a/src/derive/junction_annotation/results.rs b/src/derive/junction_annotation/results.rs index 29daf20..edd1dd3 100644 --- a/src/derive/junction_annotation/results.rs +++ b/src/derive/junction_annotation/results.rs @@ -113,7 +113,7 @@ pub struct RecordMetrics { } /// Summary statistics for the junction-annotation subcommand. -#[derive(Clone, Default, Serialize)] +#[derive(Clone, Default, Debug, Serialize)] pub struct SummaryResults { /// The total number of junctions observed in the file. pub total_junctions: usize, diff --git a/src/utils/alignment.rs b/src/utils/alignment.rs index c3cd230..678114f 100644 --- a/src/utils/alignment.rs +++ b/src/utils/alignment.rs @@ -14,9 +14,9 @@ pub fn filter_by_mapq( match min_mapq { Some(min_mapq) => match record.mapping_quality() { Some(mapq) => mapq.get() < min_mapq.get(), - None => false, + None => true, // Filter if no MAPQ is present }, - None => false, + None => false, // Do not filter if no min MAPQ is specified } } @@ -143,29 +143,29 @@ impl<'a> ReferenceRecordStepThrough<'a> { mod tests { use noodles::sam::record::{Cigar, MappingQuality, Sequence}; - use super::ReferenceRecordStepThrough; + use super::*; #[test] pub fn it_filters_by_mapq() -> anyhow::Result<()> { let mut record = noodles::sam::alignment::Record::default(); - assert!(super::filter_by_mapq( + assert!(filter_by_mapq( &record, Some(MappingQuality::new(0).unwrap()) )); // Get filtered because MAPQ is missing - assert!(!super::filter_by_mapq(&record, None)); // Do not get filtered because filter is disabled + assert!(!filter_by_mapq(&record, None)); // Do not get filtered because filter is disabled record .mapping_quality_mut() .replace(MappingQuality::new(10).unwrap()); - assert!(!super::filter_by_mapq( + assert!(!filter_by_mapq( &record, Some(MappingQuality::new(0).unwrap()) )); // Do not get filtered because MAPQ is present - assert!(!super::filter_by_mapq( + assert!(!filter_by_mapq( &record, Some(MappingQuality::new(1).unwrap()) )); // Do not get filtered because MAPQ is greater than 1 - assert!(super::filter_by_mapq( + assert!(filter_by_mapq( &record, Some(MappingQuality::new(11).unwrap()) )); // Do get filtered because MAPQ is less than 11