diff --git a/CHANGES.md b/CHANGES.md index 3197565a4..0fb087e43 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -17,6 +17,11 @@ * #686: Avoid unbounded memory usage when data cannot be written to the output as fast as reads are mapped (that is, when the output file/stdout blocks). +* #682: Report a supplementary alignment for a single-end read that crosses a + large insertion, deletion or tandem duplication, instead of soft-clipping the + part of the read that does not fit the primary alignment (options `--sm`, + `--ss`, `--sd`, `--si` and `--so`, see the README). Chains are now kept per read + interval (`--vs`, `--vn`), which also changes mapping for paired-end reads. ## v0.20.0 (2026-11-06) diff --git a/README.md b/README.md index 099047170..679cd5f45 100644 --- a/README.md +++ b/README.md @@ -281,6 +281,40 @@ controls the chaining look-back window heuristic, `--gd` sets the diagonal gap c `--gl` sets the gap length cost, `--vp` determines the best-chain score threshold, and `--sg` controls the maximum skip distance on the reference. +Chains are grouped by the read interval they cover, and only the best scoring +chains of each interval are kept: `--vp` and a cap of `--vn` chains (1000 by +default) apply per interval. A chain is also kept if it scores at least `--vs` +(40 by default), even when it scores below `--vp` of the best chain found, so a +chain that covers a part of the read the best chain misses is kept as well. + +### Supplementary alignments + +Strobealign searches for a single supplementary alignment when the primary +alignment does not explain the whole read, which happens when a read crosses a +large variation. Once the primary is extended, candidate chains are checked +against two variation patterns: for an insertion, the reference carries on +within `--si` bases (200 by default) of the primary while the read jumps the +inserted bases, and for a deletion or a tandem duplication, the read carries on +while the reference jumps up to `--sd` bases (100,000 by default). A candidate +also has to be on the same contig and the same strand as the primary, so +translocations and inversions are not detected. Candidate chains may overlap +the primary on the read by a fraction `--so` (0.5 by default) of what the +primary covers. + +Candidates are extended like the primary, but only over the read interval the +primary does not cover, so that the end bonus `-L` pulls the extension towards +the primary. Strobealign reports the best extension as a supplementary record +if it scores at least `--sm` (80 by default) and meets the primary within +`--ss` bases (50 by default), either on the read or on the reference. The +supplementary gets its own mapping quality, computed from the second best +candidate. + +Duplication detection is run on the supplementary alignment the same way it is +run on the primary. + +The supplementary alignment search runs for single-end reads only. Use +`--no-supplementary` to turn it off. + ### Duplication detection Strobealign attempts to reinterpret insertions inside an alignment as tandem diff --git a/src/chainer.rs b/src/chainer.rs index 770021096..05fc4d787 100644 --- a/src/chainer.rs +++ b/src/chainer.rs @@ -54,6 +54,8 @@ pub struct ChainingParameters { pub diag_diff_penalty: f32, pub gap_length_penalty: f32, pub valid_score_threshold: f32, + pub min_chain_score: f32, + pub max_chains_per_interval: usize, pub max_ref_gap: Option, pub matches_weight: f32, pub max_diagonal_ratio: f32, @@ -66,6 +68,8 @@ impl Default for ChainingParameters { diag_diff_penalty: 0.1, gap_length_penalty: 0.05, valid_score_threshold: 0.7, + min_chain_score: 40.0, + max_chains_per_interval: 1000, max_ref_gap: None, matches_weight: 0.01, max_diagonal_ratio: 10.0, @@ -452,6 +456,70 @@ fn hits_to_anchors(hits: &Vec, index: &StrobemerIndex) -> Vec { anchors } +/// A part of the query covered by chains, with the best chain score and the +/// number of chains admitted for it +struct QueryInterval { + start: usize, + end: usize, + best_score: f32, + n_chains: usize, +} + +/// Decides whether a chain covering `start..end` on the query is kept, and +/// records it in `intervals`. +/// +/// The chain belongs to the first interval that overlaps it by at least half +/// of the shorter of the two spans. If there is none, a new interval is added +/// for the chain and the chain is kept. +/// +/// Within its interval, a chain scoring higher than the interval's best is +/// kept, and the interval takes over the chain's span and score. Otherwise the +/// chain is rejected if it scores below `valid_score_threshold` times the +/// interval's best, or if the interval already holds `max_chains_per_interval` +/// chains. Each kept chain increments the interval's chain count. +/// +/// Returns whether the chain is kept. +fn admit_chain( + intervals: &mut Vec, + start: usize, + end: usize, + score: f32, + parameters: &ChainingParameters, +) -> bool { + for interval in intervals.iter_mut() { + let overlap = interval + .end + .min(end) + .saturating_sub(interval.start.max(start)); + if overlap * 2 < (end - start).min(interval.end - interval.start) { + continue; + } + + if score > interval.best_score { + interval.start = start; + interval.end = end; + interval.best_score = score; + } else if score < interval.best_score * parameters.valid_score_threshold + || interval.n_chains >= parameters.max_chains_per_interval + { + return false; + } + + interval.n_chains += 1; + + return true; + } + + intervals.push(QueryInterval { + start, + end, + best_score: score, + n_chains: 1, + }); + + true +} + impl ChainingResult { fn extract_chains( &self, @@ -461,7 +529,8 @@ impl ChainingResult { chains: &mut Vec, ) { let n = self.anchors.len(); - let valid_score = self.best_score * self.parameters.valid_score_threshold; + let valid_score = (self.best_score * self.parameters.valid_score_threshold) + .min(self.parameters.min_chain_score); let mut candidates = vec![]; for i in 0..n { @@ -473,12 +542,12 @@ impl ChainingResult { candidates.sort_by(|a, b| b.1.total_cmp(&a.1)); let mut used = vec![false; n]; - // Each candidate yields at most one chain - chains.reserve(candidates.len()); - for (i, score) in candidates { + let mut intervals = Vec::new(); + for (i, mut score) in candidates { if used[i] { continue; } + used[i] = true; let mut j = i; let mut overlaps = false; @@ -488,11 +557,12 @@ impl ChainingResult { let mut ref_coverage = self.anchors[i].ref_start(); while self.predecessors[j] != usize::MAX { - j = self.predecessors[j]; - if used[j] { + let predecessor = self.predecessors[j]; + if used[predecessor] { overlaps = true; break; } + j = predecessor; n_chain_anchors += 1; used[j] = true; @@ -503,6 +573,23 @@ impl ChainingResult { } if overlaps { + score = self.dp[i] - self.dp[j] + k as f32; + if score < valid_score { + continue; + } + } + + let first = &self.anchors[j]; + let last = &self.anchors[i]; + let query_start = first.query_start(); + let query_end = last.query_start() + k; + if !admit_chain( + &mut intervals, + query_start, + query_end, + score, + &self.parameters, + ) { continue; } @@ -515,13 +602,10 @@ impl ChainingResult { current = self.predecessors[current]; } - let first = &self.anchors[j]; - let last = &self.anchors[i]; - chains.push(Chain { id: chains.len(), - query_start: first.query_start(), - query_end: last.query_start() + k, + query_start, + query_end, ref_contig_start: contig_starts.ref_contig_start(first.ref_start()), ref_start: first.ref_start(), ref_end: last.ref_start() + k, @@ -600,6 +684,49 @@ pub mod test { assert!(chaining_result.best_score > 42.0); } + #[test] + fn recover_lower_scoring_sub_chain_for_supplementary() { + let chainer = Chainer::new(20, ChainingParameters::default()); + let starts = ContigStarts::new(vec![0], 100_000); + let extract = |anchors| { + let mut chains = vec![]; + chainer + .collinear_chaining(anchors, &starts, 2000) + .extract_chains(20, false, &starts, &mut chains); + chains + }; + + #[rustfmt::skip] + let chains = extract(anchors![ + (0, 0), (20, 20), (40, 40), (60, 60), (80, 80), (100, 100), + (1090, 200), (1110, 220), (1130, 240), + ]); + let alone = extract(anchors![(1090, 200), (1110, 220), (1130, 240)]); + + assert_eq!(chains.len(), 2); + assert_eq!(alone.len(), 1); + let (cut, whole) = (&chains[1], &alone[0]); + assert_eq!(cut.anchors, whole.anchors); + assert_eq!( + ( + cut.query_start, + cut.query_end, + cut.ref_start, + cut.ref_end, + cut.matching_bases + ), + ( + whole.query_start, + whole.query_end, + whole.ref_start, + whole.ref_end, + whole.matching_bases + ), + ); + + assert!((cut.score - whole.score).abs() < 1e-4); + } + #[test] fn linear_score_adjacent_anchors() { let chainer = Chainer::new(20, ChainingParameters::default()); diff --git a/src/details.rs b/src/details.rs index 74a7f24f7..c46df3bff 100644 --- a/src/details.rs +++ b/src/details.rs @@ -71,7 +71,20 @@ pub struct AlignmentDetails { /// No. of insertions re-interpreted as tandem duplications pub duplications: usize, + /// No. of chains beside a primary that fit a variant pattern and were kept as candidates + pub supplementary_candidates: usize, + + /// No. of candidates whose extension produced an alignment + pub supplementary_extended: usize, + + /// No. of supplementary alignments reported + pub supplementary: usize, + pub time_extend: f64, + + pub time_extend_primary: f64, + + pub time_extend_supplementary: f64, } impl ops::AddAssign for AlignmentDetails { @@ -84,7 +97,12 @@ impl ops::AddAssign for AlignmentDetails { self.best_alignments += rhs.best_alignments; self.duplication_candidates += rhs.duplication_candidates; self.duplications += rhs.duplications; + self.supplementary_candidates += rhs.supplementary_candidates; + self.supplementary_extended += rhs.supplementary_extended; + self.supplementary += rhs.supplementary; self.time_extend += rhs.time_extend; + self.time_extend_primary += rhs.time_extend_primary; + self.time_extend_supplementary += rhs.time_extend_supplementary; } } diff --git a/src/duplications.rs b/src/duplications.rs index 0b82589f1..62cd31151 100644 --- a/src/duplications.rs +++ b/src/duplications.rs @@ -1,13 +1,56 @@ +use std::cmp::Reverse; + use crate::aligner::{Aligner, Scores}; use crate::cigar::{Cigar, CigarOperation, OpLen}; use crate::details::AlignmentDetails; -use crate::modes::align::Alignment; +use crate::modes::align::{Alignment, AlignmentWithMapq, MappingParameters}; use crate::packed_seq::PackedSeqSlice; use crate::read::Read; use crate::refseq::RefSequence; use crate::simdaligner::SplitReferenceAlignment; -/// Cuts the alignment into one segment per copy wherever an insertion of at least `min_length` +/// Splits each alignment that contains a tandem duplication into one alignment per copy. Each +/// copy keeps the mapping quality of the alignment it came from. +/// +/// The returned alignments keep the order of the input alignments: the copies of the first +/// input alignment come first, then those of the second, and so on. The copies of one input +/// alignment are sorted by score, best first. With the primary alignment passed first, the +/// first returned alignment is therefore the best copy of the primary. +pub fn split_duplications( + aligner: &Aligner, + refseq: &RefSequence, + read: &Read, + alignments: Vec, + parameters: &MappingParameters, + details: &mut AlignmentDetails, +) -> Vec { + let mut resolved = Vec::new(); + for AlignmentWithMapq { alignment, mapq } in alignments { + let mut split = reinterpret_insertions( + aligner, + alignment, + refseq, + read, + parameters.min_duplication_length, + parameters.min_duplication_matches, + parameters.duplication_band, + details, + ); + + // Sort only within one split, so that the copies of the primary alignment stay ahead of + // those of the supplementary alignment + split.sort_by_key(|alignment| Reverse(alignment.score)); + resolved.extend( + split + .into_iter() + .map(|alignment| AlignmentWithMapq { alignment, mapq }), + ); + } + + resolved +} + +/// Cuts the alignment into one alignment per copy wherever an insertion of at least `min_length` /// bases turns out to be a tandem duplication. A `min_length` of zero turns this off. Always /// returns at least the alignment it was given. pub fn reinterpret_insertions( @@ -30,10 +73,10 @@ pub fn reinterpret_insertions( }; let contig = refseq.contig(alignment.contig_id); - let mut segments = Vec::new(); - let mut segment_cigar = Cigar::new(); - let mut segment_start = alignment.contig_start; - let mut segment_clip_left = alignment.soft_clip_left; + let mut alignments = Vec::new(); + let mut current_cigar = Cigar::new(); + let mut current_start = alignment.contig_start; + let mut current_clip_left = alignment.soft_clip_left; let mut ref_pos = alignment.contig_start; let mut query_pos = alignment.soft_clip_left; @@ -51,52 +94,52 @@ pub fn reinterpret_insertions( details, ) { - // The head arm carries the segment being built on past the insertion point, and + // The head arm carries the alignment being built on past the insertion point, and // the tail arm opens the next one at the second copy. - let mut upstream = segment_cigar.clone(); + let mut upstream = current_cigar.clone(); upstream.extend(&split.head.cigar); - if let Some(segment) = build_segment( + if let Some(head) = build_alignment( &alignment, upstream, - segment_start, - segment_clip_left, + current_start, + current_clip_left, query.len(), &aligner.scores, ) { - segments.push(segment); - segment_cigar = split.tail.cigar; - segment_start = copy_start; - segment_clip_left = query_pos + split.tail.query_start; + alignments.push(head); + current_cigar = split.tail.cigar; + current_start = copy_start; + current_clip_left = query_pos + split.tail.query_start; query_pos += len; continue; } } - segment_cigar.push(op, len); + current_cigar.push(op, len); let (query_len, ref_len) = query_ref_len(op, len); query_pos += query_len; ref_pos += ref_len; } - if segments.is_empty() { + if alignments.is_empty() { return vec![alignment]; } - let Some(last) = build_segment( + let Some(last) = build_alignment( &alignment, - segment_cigar, - segment_start, - segment_clip_left, + current_cigar, + current_start, + current_clip_left, query.len(), &aligner.scores, ) else { return vec![alignment]; }; - segments.push(last); - details.duplications += segments.len() - 1; + alignments.push(last); + details.duplications += alignments.len() - 1; - segments + alignments } /// Aligns the bases inserted at `ref_pos` against both references surrounding the @@ -136,10 +179,10 @@ fn find_duplication( Some((left_start + split.tail.ref_start, split)) } -/// Makes an alignment for a segment of a duplication split. It needs its own +/// Makes an alignment for one copy of a duplication split. It needs its own /// cigar, reference start, soft clips, score and edit distance, the rest comes /// from the original alignment. Returns None if it scores zero or less. -fn build_segment( +fn build_alignment( source: &Alignment, cigar: Cigar, ref_start: usize, diff --git a/src/lib.rs b/src/lib.rs index d35cb535c..d81f7f58a 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -22,3 +22,4 @@ pub mod seeding; pub mod shuffle; pub mod simdaligner; pub mod ssw; +pub mod supplementary; diff --git a/src/main.rs b/src/main.rs index 2999c3ff1..dd5d03b66 100644 --- a/src/main.rs +++ b/src/main.rs @@ -249,6 +249,14 @@ struct Args { #[arg(long = "vp", default_value_t = ChainingParameters::default().valid_score_threshold, help_heading = "Collinear chaining")] valid_score_threshold: f32, + /// Chains scoring at least this are kept even when below the --vp fraction of the best chain + #[arg(long = "vs", default_value_t = ChainingParameters::default().min_chain_score, help_heading = "Collinear chaining")] + min_chain_score: f32, + + /// Maximum number of chains kept over one query interval + #[arg(long = "vn", default_value_t = ChainingParameters::default().max_chains_per_interval, help_heading = "Collinear chaining")] + max_chains_per_interval: usize, + /// Maximum allowed distance between anchors on the reference [default: same as read length] #[arg(long = "sg", help_heading = "Collinear chaining")] max_ref_gap: Option, @@ -299,6 +307,32 @@ struct Args { #[arg(long = "bd", default_value_t = 1036, value_name = "N", help_heading = "Alignment")] score_drop: u32, + // Supplementary parameters + + /// Do not search for supplementary alignments + #[arg(long = "no-supplementary", help_heading = "Supplementary")] + no_supplementary: bool, + + /// Minimum score for a supplementary alignment to be reported + #[arg(long = "sm", default_value_t = MappingParameters::default().min_supplementary_score, value_name = "N", help_heading = "Supplementary")] + min_supplementary_score: u32, + + /// Maximum gap between the primary and a supplementary, on the read or the reference + #[arg(long = "ss", default_value_t = MappingParameters::default().max_supplementary_gap, value_name = "N", help_heading = "Supplementary")] + max_supplementary_gap: usize, + + /// Longest deletion or tandem duplication a supplementary alignment may describe + #[arg(long = "sd", default_value_t = MappingParameters::default().max_gap_distance, value_name = "N", help_heading = "Supplementary")] + max_gap_distance: usize, + + /// Maximum reference gap between the primary and an insertion candidate + #[arg(long = "si", default_value_t = MappingParameters::default().max_insertion_ref_gap, value_name = "N", help_heading = "Supplementary")] + max_insertion_ref_gap: usize, + + /// Fraction of the primary's read span a supplementary candidate may overlap + #[arg(long = "so", default_value_t = MappingParameters::default().max_supplementary_overlap, value_name = "FRACTION", help_heading = "Supplementary")] + max_supplementary_overlap: f32, + /// Minimum insertion length that triggers tandem duplication detection (single-end only); 0 to disable #[arg(long = "dl", default_value_t = MappingParameters::default().min_duplication_length, value_name = "N", help_heading = "Duplications")] min_duplication_length: usize, @@ -584,6 +618,12 @@ fn run() -> Result<(), CliError> { args.max_secondary as usize }, secondary_threshold: args.secondary_threshold, + max_supplementary_gap: args.max_supplementary_gap, + max_gap_distance: args.max_gap_distance, + max_insertion_ref_gap: args.max_insertion_ref_gap, + min_supplementary_score: args.min_supplementary_score, + max_supplementary_overlap: args.max_supplementary_overlap, + supplementary: !args.no_supplementary, max_tries: args.max_tries, dropoff_threshold: args.dropoff_threshold, rescue_distance: args.rescue_distance, @@ -600,6 +640,8 @@ fn run() -> Result<(), CliError> { diag_diff_penalty: args.diag_diff_penalty, gap_length_penalty: args.gap_length_penalty, valid_score_threshold: args.valid_score_threshold, + min_chain_score: args.min_chain_score, + max_chains_per_interval: args.max_chains_per_interval, max_ref_gap: args.max_ref_gap, matches_weight: args.matches_weight, max_diagonal_ratio: args.max_diagonal_ratio, @@ -958,6 +1000,21 @@ fn print_statistics( "Insertions reinterpreted as duplications: {:12}", details.align.duplications ); + debug!( + "Supplementary candidates: {:12} Per read: {:7.1}", + details.align.supplementary_candidates, + details.align.supplementary_candidates as f64 / details.chain.n_reads as f64 + ); + debug!( + "Supplementary candidates extended: {:12} Per read: {:7.1}", + details.align.supplementary_extended, + details.align.supplementary_extended as f64 / details.chain.n_reads as f64 + ); + debug!( + "Supplementary alignments reported: {:12} Per read: {:7.1}", + details.align.supplementary, + details.align.supplementary as f64 / details.chain.n_reads as f64 + ); debug!(""); debug!("## Mapping runtime breakdown"); @@ -1002,6 +1059,16 @@ fn print_statistics( details.align.time_extend, details.align.time_extend * 100. / total ); + debug!( + " primary extension: {:8.2} s {:5.1} %", + details.align.time_extend_primary, + details.align.time_extend_primary * 100. / total + ); + debug!( + " supplementary extension: {:8.2} s {:5.1} %", + details.align.time_extend_supplementary, + details.align.time_extend_supplementary * 100. / total + ); debug!( "Total: {:8.2} s {:5.1} %", total, 100. diff --git a/src/modes/align.rs b/src/modes/align.rs index b88b6573b..e659b5480 100644 --- a/src/modes/align.rs +++ b/src/modes/align.rs @@ -6,6 +6,7 @@ use std::cmp::{Reverse, min}; use std::collections::hash_map::Entry; use std::collections::{HashMap, HashSet}; use std::mem; +use std::ops::Range; use std::time::Instant; use fastrand::Rng; @@ -18,7 +19,7 @@ use crate::chain::Chain; use crate::chainer::Anchor; use crate::cigar::{Cigar, CigarOperation}; use crate::details::{AlignmentDetails, ChainDetails, Details}; -use crate::duplications::reinterpret_insertions; +use crate::duplications::split_duplications; use crate::hit::HitsDetails; use crate::index::StrobemerIndex; use crate::insertsize::InsertSizeDistribution; @@ -34,6 +35,7 @@ use crate::read::Read; use crate::refseq::{ContigPosition, RefSequence}; use crate::revcomp::reverse_complement; use crate::seeding::SeedingParameters; +use crate::supplementary::get_supplementary_alignment; const MAX_PAIR_CHAINS: usize = 1000; @@ -41,6 +43,12 @@ const MAX_PAIR_CHAINS: usize = 1000; pub struct MappingParameters { pub max_secondary: usize, pub secondary_threshold: f32, + pub max_supplementary_gap: usize, + pub max_gap_distance: usize, + pub max_insertion_ref_gap: usize, + pub min_supplementary_score: u32, + pub max_supplementary_overlap: f32, + pub supplementary: bool, pub dropoff_threshold: f32, pub rescue_distance: usize, pub max_tries: usize, @@ -57,6 +65,12 @@ impl Default for MappingParameters { MappingParameters { max_secondary: 0, secondary_threshold: 0.95, + max_supplementary_gap: 50, + max_gap_distance: 100_000, + max_insertion_ref_gap: 200, + min_supplementary_score: 80, + max_supplementary_overlap: 0.5, + supplementary: true, dropoff_threshold: 0.5, rescue_distance: 100, max_tries: 20, @@ -126,22 +140,33 @@ impl Alignment { } } +/// An alignment with its mapping quality. +pub struct AlignmentWithMapq { + pub alignment: Alignment, + pub mapq: u8, +} + +impl AlignmentWithMapq { + fn sa_entry(&self, refseq: &RefSequence, cigar_eqx: bool) -> String { + self.alignment.sa_entry(refseq, self.mapq, cigar_eqx) + } +} + fn build_sa_tags( - segments: &[Alignment], - mapq: u8, + alignments: &[AlignmentWithMapq], refseq: &RefSequence, cigar_eqx: bool, ) -> Vec> { - if segments.len() < 2 { - return vec![None; segments.len()]; + if alignments.len() < 2 { + return vec![None; alignments.len()]; } - let mut entries = Vec::with_capacity(segments.len()); - for segment in segments { - entries.push(segment.sa_entry(refseq, mapq, cigar_eqx)); + let mut entries = Vec::with_capacity(alignments.len()); + for alignment in alignments { + entries.push(alignment.sa_entry(refseq, cigar_eqx)); } - let mut tags = Vec::with_capacity(segments.len()); + let mut tags = Vec::with_capacity(alignments.len()); for i in 0..entries.len() { let mut tag = String::new(); for (j, entry) in entries.iter().enumerate() { @@ -156,7 +181,7 @@ fn build_sa_tags( } impl AlignmentDetails { - fn add_alignment(&mut self, alignment: &Alignment) { + pub fn add_alignment(&mut self, alignment: &Alignment) { self.tried_alignment += 1; self.gapped += alignment.gapped as usize; self.collisions += alignment.collisions; @@ -474,6 +499,7 @@ pub fn align_single_end_read( let mut best_alignment = None; let read = Read::new(&record.sequence); + let primary_timer = Instant::now(); for (tries, chain) in chains .iter_mut() .take(mapping_parameters.max_tries) @@ -533,6 +559,8 @@ pub fn align_single_end_read( alignments.push(alignment); } } + details.align.time_extend_primary += primary_timer.elapsed().as_secs_f64(); + let Some(best_alignment) = best_alignment else { return ( vec![sam_output.make_unmapped_record(record, details.clone())], @@ -541,32 +569,50 @@ pub fn align_single_end_read( }; let mapq = (60 * (best_score - second_best_score)).div_ceil(best_score) as u8; - let mut segments = reinterpret_insertions( + let supplementary = if mapping_parameters.supplementary { + let supplementary_timer = Instant::now(); + let supplementary = get_supplementary_alignment( + aligner, + refseq, + &read, + index.k(), + &best_alignment, + chains, + mapping_parameters, + &mut details.align, + ); + details.align.time_extend_supplementary += supplementary_timer.elapsed().as_secs_f64(); + supplementary + } else { + None + }; + + let mut reported = vec![AlignmentWithMapq { + alignment: best_alignment, + mapq, + }]; + reported.extend(supplementary); + let reported = split_duplications( aligner, - best_alignment, refseq, &read, - mapping_parameters.min_duplication_length, - mapping_parameters.min_duplication_matches, - mapping_parameters.duplication_band, + reported, + mapping_parameters, &mut details.align, ); + let mut sa_tags = build_sa_tags(&reported, refseq, sam_output.cigar_eqx); - // The best scoring segment becomes the primary and the rest supplementary - segments.sort_by_key(|segment| Reverse(segment.score)); - let mut sa_tags = build_sa_tags(&segments, mapq, refseq, sam_output.cigar_eqx); - - for i in 0..segments.len() { + for i in 0..reported.len() { let status = if i == 0 { AlignmentStatus::Primary } else { AlignmentStatus::Supplementary }; sam_records.push(sam_output.make_mapped_record( - &segments[i], + &reported[i].alignment, refseq, record, - mapq, + reported[i].mapq, status, details.clone(), sa_tags[i].take(), @@ -649,7 +695,7 @@ fn trace_log( /// Extend a chain so that it covers the entire read and return the resulting /// alignment. -fn extend_seed( +pub fn extend_seed( aligner: &Aligner, chain: &mut Chain, refseq: &RefSequence, @@ -657,19 +703,38 @@ fn extend_seed( use_ssw: bool, k: usize, ) -> Option { - let query = if chain.is_revcomp { + extend_seed_in_interval(aligner, chain, refseq, read, use_ssw, k, 0..read.len()) +} + +/// Extend a chain over the given interval of the read (in the chain's +/// orientation) and return the resulting alignment. Read bases outside the +/// interval are soft clipped. +pub fn extend_seed_in_interval( + aligner: &Aligner, + chain: &mut Chain, + refseq: &RefSequence, + read: &Read, + use_ssw: bool, + k: usize, + read_interval: Range, +) -> Option { + let full = if chain.is_revcomp { read.rc() } else { read.seq() }; + let start = read_interval.start; + let query = &full[read_interval]; + let chain_start = chain.query_start.saturating_sub(start); + let chain_end = chain.query_end.saturating_sub(start).min(query.len()); let (contig_id, contig_start) = refseq.unflatten(chain.ref_start); let contig_start = contig_start.0 as usize; let contig_end = contig_start + chain.ref_span(); let contig = refseq.contig(contig_id); assert!(contig_end <= contig.len()); - let projected_start = contig_start.saturating_sub(chain.query_start); - let projected_end = min(contig_end + query.len() - chain.query_end, contig.len()); + let projected_start = contig_start.saturating_sub(chain_start); + let projected_end = min(contig_end + (query.len() - chain_end), contig.len()); // TODO ugly let mut info = AlignmentInfo { @@ -722,11 +787,14 @@ fn extend_seed( .iter() .filter_map(|a| { let ref_start = a.ref_start() - decode_start - refseq.starts.0[contig_id]; + // Remove anchors outside the interval + let query_start = a.query_start().checked_sub(start)?; + if query_start + k > query.len() { + return None; + } // Skip anchors resulting from hash collisions - if query[a.query_start()..a.query_start() + k] - == segment[ref_start..ref_start + k] - { - Some(Anchor::new(ref_start, a.query_start())) + if query[query_start..query_start + k] == segment[ref_start..ref_start + k] { + Some(Anchor::new(ref_start, query_start)) } else { collisions += 1; @@ -742,8 +810,8 @@ fn extend_seed( Some(Alignment { cigar: mem::take(&mut info.cigar), edit_distance: info.edit_distance, - soft_clip_left: info.query_start, - soft_clip_right: query.len() - info.query_end, + soft_clip_left: start + info.query_start, + soft_clip_right: full.len() - start - info.query_end, score: info.score, contig_id, contig_start: result_start, @@ -1631,11 +1699,25 @@ mod tests { collisions: 0, }; - let tags = build_sa_tags(&[first.clone(), second], 60, &refseq, true); + let alignments = [ + AlignmentWithMapq { + alignment: first.clone(), + mapq: 60, + }, + AlignmentWithMapq { + alignment: second, + mapq: 60, + }, + ]; + let tags = build_sa_tags(&alignments, &refseq, true); assert_eq!(tags[0].as_deref(), Some("chr1,11,-,40S70=,60,0;")); assert_eq!(tags[1].as_deref(), Some("chr1,101,+,40=70S,60,0;")); - assert_eq!(build_sa_tags(&[first], 60, &refseq, true), vec![None]); + let only = [AlignmentWithMapq { + alignment: first, + mapq: 60, + }]; + assert_eq!(build_sa_tags(&only, &refseq, true), vec![None]); } #[test] diff --git a/src/supplementary.rs b/src/supplementary.rs new file mode 100644 index 000000000..de96ff313 --- /dev/null +++ b/src/supplementary.rs @@ -0,0 +1,160 @@ +use crate::aligner::Aligner; +use crate::chain::Chain; +use crate::details::AlignmentDetails; +use crate::modes::align::{ + Alignment, AlignmentWithMapq, MappingParameters, extend_seed_in_interval, +}; +use crate::read::Read; +use crate::refseq::RefSequence; + +/// Extend the chains that fit the pattern of a large variation with the +/// primary alignment and return the best one. Returns `None` if no chain fits, +/// or if the best extension scores below `min_supplementary_score`. +pub fn get_supplementary_alignment( + aligner: &Aligner, + refseq: &RefSequence, + read: &Read, + k: usize, + primary: &Alignment, + chains: &mut [Chain], + parameters: &MappingParameters, + details: &mut AlignmentDetails, +) -> Option { + let mut candidates: Vec<&mut Chain> = chains + .iter_mut() + .filter(|chain| { + fits_variant_pattern( + primary, + refseq, + chain, + read.len(), + parameters.max_gap_distance as i64, + parameters.max_insertion_ref_gap as i64, + parameters.max_supplementary_overlap, + ) + }) + .collect(); + details.supplementary_candidates += candidates.len(); + + let dropoff = candidates.first()?.score * parameters.dropoff_threshold; + let max_gap = parameters.max_supplementary_gap as i64; + let mut best_alignment = None; + let mut best_score = 0; + let mut second_best_score = 0; + + let primary_start = primary.soft_clip_left; + let primary_end = read.len() - primary.soft_clip_right; + for chain in candidates.iter_mut().take(parameters.max_tries) { + if chain.score < dropoff { + break; + } + + let interval = if chain.query_start >= primary_start { + primary_end..read.len() + } else { + 0..primary_start + }; + let Some(alignment) = extend_seed_in_interval( + aligner, + chain, + refseq, + read, + parameters.use_ssw, + k, + interval, + ) else { + continue; + }; + details.supplementary_extended += 1; + details.add_alignment(&alignment); + + let (read_gap, ref_gap) = alignment_gaps(primary, &alignment, read.len()); + if read_gap.abs() > max_gap && ref_gap.abs() > max_gap { + continue; + } + + if alignment.score > best_score { + second_best_score = best_score; + best_score = alignment.score; + best_alignment = Some(alignment); + } else if alignment.score > second_best_score { + second_best_score = alignment.score; + } + } + + if best_score < parameters.min_supplementary_score { + return None; + } + + let alignment = best_alignment?; + let mapq = (60 * best_score.saturating_sub(second_best_score)).div_ceil(best_score) as u8; + details.supplementary += 1; + + Some(AlignmentWithMapq { alignment, mapq }) +} + +/// Whether a chain and the primary alignment are the two sides of a deletion, insertion +/// or duplication. +fn fits_variant_pattern( + primary: &Alignment, + refseq: &RefSequence, + chain: &Chain, + read_len: usize, + max_gap_distance: i64, + max_insertion_ref_gap: i64, + overlap_fraction: f32, +) -> bool { + let (contig_id, contig_start) = refseq.unflatten(chain.ref_start); + let contig_start = contig_start.0 as usize; + if primary.is_revcomp != chain.is_revcomp || primary.contig_id != contig_id { + return false; + } + + let (read_gap, ref_gap) = if chain.query_start >= primary.soft_clip_left { + ( + chain.query_start as i64 - (read_len - primary.soft_clip_right) as i64, + contig_start as i64 - (primary.contig_start + primary.length) as i64, + ) + } else { + ( + primary.soft_clip_left as i64 - chain.query_end as i64, + primary.contig_start as i64 - (contig_start + chain.ref_span()) as i64, + ) + }; + + let primary_span = (read_len - primary.soft_clip_left - primary.soft_clip_right) as i64; + let overlap_tolerance = (primary_span as f32 * overlap_fraction) as i64; + if read_gap < -overlap_tolerance { + return false; + } + + // An insertion: the contig carries straight on while the read jumps the + // inserted bases + if (0..=max_insertion_ref_gap).contains(&ref_gap) { + return true; + } + + // A deletion/duplication: the read carries straight on while the contig + // jumps bases + if read_gap <= overlap_tolerance && (-max_gap_distance..max_gap_distance).contains(&ref_gap) { + return true; + } + + false +} + +/// Returns the gaps between two alignments on the read and on the reference, +/// as (read_gap, ref_gap). A negative gap means the alignments overlap. +fn alignment_gaps(primary: &Alignment, candidate: &Alignment, read_len: usize) -> (i64, i64) { + if candidate.soft_clip_left >= primary.soft_clip_left { + ( + candidate.soft_clip_left as i64 - (read_len - primary.soft_clip_right) as i64, + candidate.contig_start as i64 - (primary.contig_start + primary.length) as i64, + ) + } else { + ( + primary.soft_clip_left as i64 - (read_len - candidate.soft_clip_right) as i64, + primary.contig_start as i64 - (candidate.contig_start + candidate.length) as i64, + ) + } +} diff --git a/tests/phix.pe.sam b/tests/phix.pe.sam index b4b2e4458..e8476211c 100644 --- a/tests/phix.pe.sam +++ b/tests/phix.pe.sam @@ -23,7 +23,7 @@ SRR1377138.10 99 NC_001422.1 1 60 26S275= = 58 358 NAAATGATTGGCGTATCCAACCTGCAGAG SRR1377138.10 147 NC_001422.1 58 60 300=1X = 1 -358 GAAAAATTATCTTGATAAAGCAGGAATTACTACTGCTTGTTTACGAATTAAATCGAAGTGGACTGCTGGCGGAAAATGAGAAAATTCGACCTATCCTTGCGCAGCTCGAGAAGCTCTTACTTTGCGACCTTTCGCCATCAACTAACGATTCTGTCAAAAACTGACGCGTTGGATGAGGAGAAGTGGCTTAATATGCTTGGCACGTTCGTCAAGGACTGGTTTAGATATGAGTCACATTTTGTTCATGGTAGAGATTCTCTTGTTGACATTTTAAAAGAGCGTGGATTACTATCTGAGTCCN ###FBDA<5)2?>FFDFFECC49FFD?3<@FFFFEDEC?3E?C?0EFFFFFFF>@5FECFFGGGFGD:FGGCDGFGGGGGG>GGGGGGGDGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGEGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGCCB8# NM:i:1 AS:i:612 na:i:1 al:i:1 ga:i:0 X0:i:1 du:i:0 mr:i:0 RG:Z:1 SRR1377138.11 99 NC_001422.1 5117 60 1X268=32S = 5220 270 NAAGCTGTTCAGAATCAGAATGAGCCGCAACTTCGGGATGAAAATGCTCACAATGACAAATCTGTCCACGGAGTGCTTAATCCAACTTACCAAGCTGGGTTACGACGCGACGCCGTTCAACCAGATATTGAAGCAGAACGCAAAAAGAGAGATGAGATTGAGGCTGGGAAAAGTTACTGTAGCCGACGTTTTGGCGGCGCAACCTGTGACGACAAATCTGCTCAAATTTATGCGCGCTTCGATAAAAATGATTGGCGTATCCAACCTGCGGAGTTTTATCGCTTCCATGACGCAGAAGTTA #8ACCGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFFGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGDGGGGF8FFGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGFGGGCFGGGGGGGGGGEE5EGGGGGGGGGFGGGGFDGEGDGGDGGGE7FFGGG>*?*:<>9BDECDFGF<*9@F?*7B?;F;F;B:AFFFFFF7*(03:0:AAAFFF:18?BA:82><9B>4<2BAA2A<<:?:4,.(B?C?.37;)EE<;)C?6(((;FFFFFEE?7)?9/);FDFB>A=FDAFFFFFFFDEDDCCFAGFFGGGGGGF:GGGFGGGGGFGGGGFGCGGEECDEFCCGGFGGGGFGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGFFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGCCA8# NM:i:0 AS:i:344 na:i:2 al:i:1 ga:i:1 X0:i:1 du:i:0 mr:i:0 RG:Z:1 -SRR1377138.12 99 NC_001422.1 1 60 57S238=1X5= = 28 328 NAATCTGCTCAAATTTATGCGCGCTTCGATAAAAATGATTGGCGTATCCAACCTGCAGAGTTTTATCGCTTCCATGACGCAGAAGTTAACACTTTCGGATATTTCTGATGAGTCGAAAAATTATCTTGATAAAGCAGGAATTACTACTGCTTGTTTACGAATTAAATCGAAGTGGACTGCTGGCGGAAAATGAGAAAATTCGACCTATCCTTGCGCAGCTCGAGAAGCTCTTACTTTGCGACCTTTCGCCATCAACTAACGATTCTGTCAAAAACTGACGCGTTGGATGAGGAGATGTGGC #8BCCGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGDCGGGGGGGGGGFGGGGGGGFGFGGFGFGFFFFFFFFFFFAFFFFFFFFFFFA3.1.48?A<:?FFFFFFBDFFFF?1GGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGCCB8# NM:i:2 AS:i:602 na:i:1 al:i:1 ga:i:0 X0:i:1 du:i:0 mr:i:0 RG:Z:1 SRR1377138.13 99 NC_001422.1 5006 60 1X300= = 5049 344 NTATGTGGCTAAATACGTTAACAAAAAGTCAGATATGGACCTTGCTGCTAAAGGTCTAGGAGCTAAAGAATGGAACAACTCACTAAAAACCAAGCTGTCGCTACTTCCCAAGAAGCTGTTCAGAATCAGAATGAGCCGCAACTTCGGGATGAAAATGCTCACAATGACAAATCTGTCCACGGAGTGCTTAATCCAACTTACCAAGCTGGGTTACGACGCGACGCCGTTCAACCAGATATTGAAGCAGAACGCAAAAAGAGAGATGAGATTGAGGCTGGGAAAAGTTACTGTAGCCGACGTT #8BCCGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFCFGGFGFGFFFFFFFFFFFDFFFFB:>FFFFF2.1.4::?FFFFF?BAAFFFFFFDFF>>9>9? NM:i:1 AS:i:612 na:i:1 al:i:1 ga:i:0 X0:i:1 du:i:0 mr:i:0 RG:Z:1 SRR1377138.13 147 NC_001422.1 5049 60 300=1X = 5006 -344 GCTGCTAAAGGTCTAGGAGCTAAAGAATGGAACAACTCACTAAAAACCAAGCTGTCGCTACTTCCCAAGAAGCTGTTCAGAATCAGAATGAGCCGCAACTTCGGGATGAAAATGCTCACAATGACAAATCTGTCCACGGAGTGCTTAATCCAACTTACCAAGCTGGGTTACGACGCGACGCCGTTCAACCAGATATTGAAGCAGAACGCAAAAAGAGAGATGAGATTGAGGCTGGGAAAAGTTACTGTAGCCGACGTTTTGGCGGCGCAACCTGTGACGACAAATCTGCTCAAATTTATGN #FFFFBFFFBD>A<::64>AAFFFB><?FFBABBAFFFFFFBAFFFFFFFFFFFFFE;@F@FFFFFFFD@FE?FFFFFFFFFFFFFEBFFFFFFFFGFFGGGDGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGDGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGCGGCGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGCCC8# NM:i:1 AS:i:612 na:i:1 al:i:1 ga:i:0 X0:i:1 du:i:0 mr:i:0 RG:Z:1 @@ -59,8 +59,8 @@ SRR1377138.28 99 NC_001422.1 3016 60 1X116=1X183= = 3273 558 NCTCCGGTGGCATTCAAGG SRR1377138.28 147 NC_001422.1 3273 60 52=1X247=1X = 3016 -558 AAGGATACTCGTGATTATCTTGCTGCTGCATTTCCTGAGCTTAATGCTTGGGGGCGTGCTGGTGCTGATGCTTCCTCTGCTGGTATGGTTGACGCCGGATTTGAGAATCAAAAAGAGCTTACTAAAATGCAACTGGACAATCAGAAAGAGATTGCCGAGATGCAAAATGAGACTCAAAAAGAGATTGCTGGCATTCAGTCGGCGACTTCACGCCAGAATACGAAAGACCAGGTATATGCACAAAATGAGATGCTTGCTTATCAACAGAAGGAGTCTACTGCTCGCGTTGCGTCTATTATGN #####F?9FFFD<4<6,22:(FFFDA((FFFFEAFFFEAFDDFFFFEEFAFFFFFCDGGGGGGGFDDFGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGEFFCGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGFGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGCCA8# NM:i:2 AS:i:602 na:i:1 al:i:1 ga:i:0 X0:i:1 du:i:0 mr:i:0 RG:Z:1 SRR1377138.29 99 NC_001422.1 155 60 1X300= = 250 396 NGCGCAGCTCGAGAAGCTCTTACTTTGCGACCTTTCGCCATCAACTAACGATTCTGTCAAAAACTGACGCGTTGGATGAGGAGAAGTGGCTTAATATGCTTGGCACGTTCGTCAAGGACTGGTTTAGATATGAGTCACATTTTGTTCATGGTAGAGATTCTCTTGTTGACATTTTAAAAGAGCGTGGATTACTATCTGAGTCCGATGCTGTTCAACCACTAATAGGTAAGAAATCATGAGTCAAGTTACTGAACAATCCGTACGTTTCCAGACCGCTTTGGCCTCTATTAAGCTCATTCAG #8BCCGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGFGGGGGGGGGGGGFGGGGGGGDGGGGGGGGGGFFGFGGGGGGGGGGGGGGGGGEGGGGGFGFGGFGGFGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFFGFACFGGFGEFGGGGFEGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGFGGGGGGG6C<>=4A:AEDFFFFFFFFFFFFFGFFGGGGGGGGGGGGGGGGGGGGGGGGGGCFGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGFGGGGFCGEGGGGGGGGGGD9FGGGGGGGGGFGFFGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGEGGFCCB8# NM:i:3 AS:i:592 na:i:1 al:i:1 ga:i:0 X0:i:1 du:i:0 mr:i:0 RG:Z:1 -SRR1377138.30 81 NC_001422.1 1 60 97S203=1X = 5176 5386 TTACTGTAGCCGACGTTTTGGCGGCGCAACCTGTGACGACAAATCTGCTCAAATTTATGCGCGCTTCGATAAAAATGATTGGCGTATCCAACCTGCAGAGTTTTATCGCTTCCATGACGCAGAAGTTAACACTTTCGGATATTTCTGATGAGTCGAAAAATTATCTTGATAAAGCAGGAATTACTACTGCTTGTTTACGAATTAAATCGAAGTGGACTGCTGGCGGAAAATGAGAAAATTCGACCTATCCTTGCGCAGCTCGAGAAGCTCTTACTTTGCGACCTTTCGCCATCAACTAACN FFFBFF?4(>7,.6FF?;96FFFB>??450.3:FFFFFBB,>EGFGGGGFGCEGCFEEEGDEGGGFGGGCF4CFGGGGC3*85CFF4:3?@FF6*.9*96<4:>?F24><?FB9AAFBFF2796;??BFFF99:A96-)))46?F>(:0<>4>5<2))-5>BF7,.6FF?;96FFFB>??450.3:FFFFFBB,>EGFGGGGFGCEGCFEEEGDEGGGFGGGCF4CFGGGGC3*85CFF4:3?@FF6*.9*96<4:>?F24><?FB9AAFBFF2796;??BFFF99:A96-)))46?F>(:0<>4>5<2))-5>BFBFBFFFAF?FF9B? NM:i:1 AS:i:612 na:i:1 al:i:1 ga:i:0 X0:i:1 du:i:0 mr:i:0 RG:Z:1 SRR1377138.31 147 NC_001422.1 4401 60 300=1X = 4320 -382 GCTGAACGCCCTCTTAAGGATATTCGCGATGAGTATAATTACCCCAAAAAGAAAGGTATTAAGGATGAGTGTTCAAGATTGCTGGAGGCCTCCACTATGAAATCGCGTAGAGGCTTTGCTATTCAGCGTTTGATGAATGCAATGCGACAGGCTCATGCTGATGGTTGGTTTATCGTTTTTGACACTCTCACGTTGGCTGACGACCGATTAGAGGCGTTTTATGATAATCCCAATGCTTTGCGTGACTATTTTCGTGATATTGGTCGTATGGTTCTTGCTGCCGAGGGTCGCAAGGCTAATN ################AA4)19>B913<96.))6<6<976;6B>9EGGGGGGGGGFGFEF9GEGGGGGGGGGDGGGFGGGGGGGGGGGGGGGFGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGCCB8# NM:i:1 AS:i:612 na:i:1 al:i:1 ga:i:0 X0:i:1 du:i:0 mr:i:0 RG:Z:1 SRR1377138.32 99 NC_001422.1 1433 60 1X268=1X31= = 1487 355 NCCGAGGCTAACCCTAATGAGCTTAATCAAGATGATGCTCGTTATGGTTTCCGTTGCTGCCATCTCAAAAACATTTGGACTGCTCCGCTTCCTCCTGAGACTGAGCTTTCTCGCCAAATGACGACTTCTACCACATCTATTGACATTATGGGTCTGCAAGCTGCTTATGCTAATTTGCATACTGACCAAGAACGTGATTACTTCATGCAGCGTTACCATGATGTTATTTCTTCATTTGGAGGTAAAACCTCTTATGACGCTGACAACCGCCCTTTACTTGTCATGCGCTCTAATCTCTGGG #8BCCGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGFGGFGFGFFFFFFFFFFFFFFFFFFFFFFFF*.:6=;CFFFFFEFFEF9;(3:AFFB)66>?: NM:i:2 AS:i:602 na:i:1 al:i:1 ga:i:0 X0:i:1 du:i:0 mr:i:0 RG:Z:1 diff --git a/tests/phix.se.sam b/tests/phix.se.sam index f091aab66..40aba8fcc 100644 --- a/tests/phix.se.sam +++ b/tests/phix.se.sam @@ -12,7 +12,8 @@ SRR1377138.8 16 NC_001422.1 1254 60 124=1X132=1X42=1X * 0 0 TCATGAAGGATGGTGTTAAT SRR1377138.9 16 NC_001422.1 4731 60 300=1X * 0 0 TGTGTGCCTGAGTATGGTACAGCTAATGGCCGTCTTCATTTCCATGCGGTGCACTTTATGCGGACACTTCCTACAGGTAGCGTTGACCCTAATTTTGGTCGTCGGGTACGCAATCGCCGCCAGTTAAATAGCTTGCAAAATACGTGGCCTTATGGTTACAGTATGCCCATCGCAGTTCGCTACACGCAGGACGCTTTTTCACGTTCTGGTTGGTTGTGGCCTGTTGATGCTAAAGGTGAGCCGCTTAAAGCTACCAGTTATATGGCTGTTGGTTTCTATGTGGCTAAATACGTTAACAAAN ?BFF?<:?*?*:<>9BDECDFGF<*9@F?*7B?;F;F;B:AFFFFFF7*(03:0:AAAFFF:18?BA:FFFFF2.1.4::?FFFFF?BAAFFFFFFDFF>>9>9? NM:i:1 AS:i:612 na:i:1 al:i:1 ga:i:0 X0:i:0 du:i:0 RG:Z:1 SRR1377138.14 0 NC_001422.1 66 60 1X300= * 0 0 NTCTTGATAAAGCAGGAATTACTACTGCTTGTTTACGAATTAAATCGAAGTGGACTGCTGGCGGAAAATGAGAAAATTCGACCTATCCTTGCGCAGCTCGAGAAGCTCTTACTTTGCGACCTTTCGCCATCAACTAACGATTCTGTCAAAAACTGACGCGTTGGATGAGGAGAAGTGGCTTAATATGCTTGGCACGTTCGTCAAGGACTGGTTTAGATATGAGTCACATTTTGTTCATGGTAGAGATTCTCTTGTTGACATTTTAAAAGAGCGTGGATTACTATCTGAGTCCGATGCTGTT #8@CCGGGGGGFFGGGGGGGGGGGGGFGGGGGGGGGGGGGFGGEEFGGGGGGGGGGGGGFGDCFEFGGGGGFGGGGGGGGGGGGGGGGGGGG7FGGGGGGGGGGGGGGGGGGGDFGGGGGGGFGGEGGGGGFGGGGGGGGGGGGGGGGGFGGGGGGGGGGDGGGGGGGGDGGGGGGGFEGGGGGGGFGGFFGFGGGGGGGGGFGGGFG7FGGFGDGGGGGFFCGGGGGGGGGGGGGFGFGGC>FGCF=A@FFFFFFFFFFEFFFFFFE<=/04(34:501:AB>4 NM:i:1 AS:i:612 na:i:1 al:i:1 ga:i:0 X0:i:0 du:i:0 RG:Z:1 SRR1377138.15 0 NC_001422.1 1907 60 1X185=115S * 0 0 NTTAACGCTAAAGGTGCTTTGACTTATACCGATATTGCTGGCGACCCTGTTTTGTATGGCAACTTGCCGCCGCGTGAAATTTCTATGAAGGATGTTTTCCGTTCTGGTGATTCGTCTAAGAAGTTTAAGATTGCTGAGGGTCAGTGGTATCGTTATGCGCCTTCGTATGTTTCTCCTGCTTATCACAAGAACGGAAAACATCCTTCATAGAAATTTCACGCAGATCGGAAGAGCGGTTCAGCAGGAATGCCGAGACCGATCTCGTATGCCGTCTTCTGCTTGAAAAAAAAAAAAAAAGATA #8BCCGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGEGFGGFGFGFFFFFFFFFFFFF>73?FFFFFFF3-048:BDFFAFFFFFF06>BFF######### NM:i:1 AS:i:372 na:i:1 al:i:1 ga:i:1 X0:i:0 du:i:0 RG:Z:1 @@ -30,7 +31,8 @@ SRR1377138.26 16 NC_001422.1 1642 60 300=1X * 0 0 GCGTTACCATGATGTTATTTCTTCATTTGG SRR1377138.27 16 NC_001422.1 4289 60 31=1X268=1X * 0 0 TGTGCTCCCCCAACTTGATATTAATAACACTTTAGACCACCGCCCCGAAGGGGACGAAAAATGGTTTTTAGAGAACGAGAAGACGGTTACGCAGTTTTGCCGCAAGCTGGCTGCTGAACGCCCTCTTAAGGATATTCGCGATGAGTATAATTACCCCAAAAAGAAAGGTATTAAGGATGAGTGTTCAAGATTGCTGGAGGCCTCCACTATGAAATCGCGTAGAGGCTTTGCTATTCAGCGTTTGATGAATGCAATGCGACAGGCTCATGCTGATGGTTGGTTTATCGTTTTTGACACTCTN 0FFF>>96BBAE9FFFF@FFFDFFFFGFGFGGDGFGGDGGGGGGGGGFGGGGEGGGGGGGGGGGGGGGGEGDFDGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFFFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFFGFGGGFGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGEGGGGGGGGGGFGGGGGGGGGGGGGGGGF?GGGGGGGGGGGGCCA8# NM:i:2 AS:i:602 na:i:1 al:i:1 ga:i:0 X0:i:0 du:i:0 RG:Z:1 SRR1377138.28 0 NC_001422.1 3016 60 1X116=1X183= * 0 0 NCTCCGGTGGCATTCAAGGTGATGTGCTTGCTACCGATAACAATACTGTAGGCATGGGTGATGCTGGTATTAAATCTGCCATTCAAGGCTCTAATGTTCCTAACCCTGATGAGGCCGTCCCTAGTTTTGTTTCTGGTGCTATGGCTAAAGCTGGTAAAGGACTTCTTGAAGGTACGTTGCAGGCTGGCACTTCTGCCGTTTCTGATAAGTTGCTTGATTTGGTTGGACTTGGTGGCAAGTCTGCCGCTGATAAAGGAAAGGATACTCGTGATTATCTTGCTGCTGCATTTCCTGAGCTTAA #8ACCGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGCFGGGGGGGGGGGGGGFFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGEGGGGGGGGGGFGGGGGGGFGGFGGGGGDFGGFGGGGGFGFGGFGFGFFFFFFFFFFF:FEF@@F@=FEF@8/;?9))9CCEFFFFEFFEFF<7,.6FF?;96FFFB>??450.3:FFFFFB7,.6FF?;96FFFB>??450.3:FFFFFB7,.6FF?;96FFFB>??450.3:FFFFFBBFBFFFAF?FF9B? NM:i:1 AS:i:612 na:i:1 al:i:1 ga:i:0 X0:i:0 du:i:0 RG:Z:1 SRR1377138.32 0 NC_001422.1 1433 60 1X268=1X31= * 0 0 NCCGAGGCTAACCCTAATGAGCTTAATCAAGATGATGCTCGTTATGGTTTCCGTTGCTGCCATCTCAAAAACATTTGGACTGCTCCGCTTCCTCCTGAGACTGAGCTTTCTCGCCAAATGACGACTTCTACCACATCTATTGACATTATGGGTCTGCAAGCTGCTTATGCTAATTTGCATACTGACCAAGAACGTGATTACTTCATGCAGCGTTACCATGATGTTATTTCTTCATTTGGAGGTAAAACCTCTTATGACGCTGACAACCGCCCTTTACTTGTCATGCGCTCTAATCTCTGGG #8BCCGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGFGGFGFGFFFFFFFFFFFFFFFFFFFFFFFF*.:6=;CFFFFFEFFEF9;(3:AFFB)66>?: NM:i:2 AS:i:602 na:i:1 al:i:1 ga:i:0 X0:i:0 du:i:0 RG:Z:1 SRR1377138.33 0 NC_001422.1 3817 60 1X268=1X31= * 0 0 NTATTTTTCATGGTATTGATAAAGCTGTTGCCGATACTTGGAACAATTTCTGGAAAGACGGTAAAGCTGATGGTATTGGCTCTAATTTGTCTAGGAAATAACCGTCAGGATTGACACCCTCCCAATTGTATGTTTTCATGCCTCCAAATCTTGGAGGCTTTTTTATGGTTCGTTCTTATTACCCTTCTGAATGTCACGCTGATTATTTTGACTTTGAGCGTATCGAGGCTCTTAAACCTGCTATTGAGGCTTGTGGCATTTCTACTCTTCCTCAATCCCCAATGCTTGGCTTCCATAAGCA #8BCCGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGFEFGGGGGGGFGGFGGGGGGFFFGGFGFGFFFAFFFFFFBDFEFFEFFFFFFC,)217@EFFFDF<@EFFFF?FFFDE####### NM:i:2 AS:i:602 na:i:1 al:i:1 ga:i:0 X0:i:0 du:i:0 RG:Z:1