/// Sparse DP routine generalizing LCSk++ method above to penalize alignment gaps.
/// A gap is an unknown combination of mismatch, insertion and deletions, and incurs
/// a penalty of gap_open + d * gap_extend, where d is the distance along the diagonal of the gap.
/// # Arguments
///
/// * `matches` - a vector of tuples indicating the (string1 position, string2 position) kmer
///   matches between the strings
/// * `k` - the kmer length used for matching
/// * `match_score` - reward for each matched base
/// * `gap_open` - score of opening a gap, including a mismatch gap. Must be negative.
/// * `gap_extend` - score for extending a gap along the diagonal. Must be negative.
///
/// # Return value
///
/// The method returns a `SparseAlignmentResult` struct with the following fields:
/// * `path` is the SDP path, represented as vector of indices into the input matches vector.
/// * `score` is the score of the path, which is the number of bases covered by the matched kmers.
/// * `dp_vector` is the full DP vector, which can generally be ignored. (It may be useful for
///   testing purposes).
Sparse DP routine generalizing LCSk++ method above to penalize alignment gaps.
A gap is an unknown combination of mismatch, insertion and deletions, and incurs
a penalty of gap_open + d * gap_extend, where d is the distance along the diagonal of the gap.
Sparse DP routine generalizing LCSk ++ method above to penalize alignment gaps . A gap is an unknown combination of mismatch insertion and deletions and incurs a penalty of gap_open + d * gap_extend where d is the distance along the diagonal of the gap .
pub fn sdpkpp(
    matches: &[(u32, u32)],
    k: usize,
    match_score: u32,
    gap_open: i32,
    gap_extend: i32,
) -> SparseAlignmentResult {
    if matches.is_empty() {
        return SparseAlignmentResult {
            path: Vec::new(),
            score: 0,
            dp_vector: Vec::new(),
        };
    }
    let k = k as u32;
    if gap_open > 0 || gap_extend > 0 {
        panic!("gap parameters cannot be positive")
    }
    let _gap_open = (-gap_open) as u32;
    let _gap_extend = (-gap_extend) as u32;
    // incoming matches must be sorted to let us find the predecessor kmers by binary search.
    for i in 1..matches.len() {
        assert!(matches[i - 1] < matches[i]);
    }
    let mut events: Vec<(u32, u32, u32)> = Vec::new();
    let mut n = 0;
    for (idx, &(x, y)) in matches.iter().enumerate() {
        events.push((x, y, (idx + matches.len()) as u32));
        events.push((x + k, y + k, idx as u32));
        n = max(n, x + k);
        n = max(n, y + k);
    }
    events.sort_unstable();
    let mut max_col_dp: MaxBitTree<PrevPtr> = MaxBitTree::new(n as usize);
    let mut dp: Vec<(u32, i32)> = Vec::with_capacity(events.len());
    let mut best_dp = (k, 0);
    dp.resize(events.len(), (0, 0));
    for ev in events {
        let p = (ev.2 % matches.len() as u32) as usize;
        let j = ev.1;
        let is_start = ev.2 >= (matches.len() as u32);
        if is_start {
            // Default case -- chain starts at this node
            dp[p] = (k * match_score, -1);
            // Find best previous chain, and extend.
            let best_prev = max_col_dp.get(j as usize);
            if best_prev.score > 0 {
                let prev_x = best_prev.x;
                let prev_y = best_prev.y;
                let cur_x = ev.0;
                let cur_y = ev.1;
                let gap = max(cur_x - prev_x, cur_y - prev_y);
                let gap_penalty = if gap > 0 {
                    _gap_open + gap * _gap_extend
                } else {
                    0
                };
                let reward = k * match_score;
                let new_score = (best_prev.score + reward).saturating_sub(gap_penalty);
                dp[p] = max(dp[p], (new_score, best_prev.id as i32));
                best_dp = max(best_dp, (dp[p].0, p as i32));
            }
        } else {
            // See if this kmer continues a different kmer
            if ev.0 > k && ev.1 > k {
                if let Ok(cont_idx) = matches.binary_search(&(ev.0 - k - 1, ev.1 - k - 1)) {
                    let prev_score = dp[cont_idx].0;
                    let candidate = (prev_score + match_score, cont_idx as i32);
                    dp[p] = max(dp[p], candidate);
                    best_dp = max(best_dp, (dp[p].0, p as i32));
                }
            }
            let prev_frag = PrevPtr::new(dp[p].0, ev.0, ev.1, p, _gap_extend);
            max_col_dp.set(ev.1 as usize, prev_frag);
        }
    }
    let mut traceback = Vec::new();
    let (best_score, mut prev_match) = best_dp;
    while prev_match >= 0 {
        traceback.push(prev_match as usize);
        prev_match = dp[prev_match as usize].1;
    }
    traceback.reverse();
    SparseAlignmentResult {
        path: traceback,
        score: best_score,
        dp_vector: dp,
    }
}
pub fn sdpkpp ( matches : & [ ( u32 , u32 ) ] , k : usize , match_score : u32 , gap_open : i32 , gap_extend : i32 , ) -> SparseAlignmentResult { if matches . is_empty ( ) { return SparseAlignmentResult { path : Vec :: new ( ) , score : 0 , dp_vector : Vec :: new ( ) , } ; } let k = k as u32 ; if gap_open > 0 || gap_extend > 0 { panic ! ( " " ) } let _gap_open = ( - gap_open ) as u32 ; let _gap_extend = ( - gap_extend ) as u32 ; // incoming matches must be sorted to let us find the predecessor kmers by binary search. for i in 1 .. matches . len ( ) { assert ! ( matches [ i 1 ] matches [ i ] ) ; } let mut events : Vec < ( u32 , u32 , u32 ) > = Vec :: new ( ) ; let mut n = 0 ; for ( idx , & ( x , y ) ) in matches . iter ( ) . enumerate ( ) { events . push ( ( x , y , ( idx + matches . len ( ) ) as u32 ) ) ; events . push ( ( x + k , y + k , idx as u32 ) ) ; n = max ( n , x + k ) ; n = max ( n , y + k ) ; } events . sort_unstable ( ) ; let mut max_col_dp : MaxBitTree < PrevPtr > = MaxBitTree :: new ( n as usize ) ; let mut dp : Vec < ( u32 , i32 ) > = Vec :: with_capacity ( events . len ( ) ) ; let mut best_dp = ( k , 0 ) ; dp . resize ( events . len ( ) , ( 0 , 0 ) ) ; for ev in events { let p = ( ev . 2 % matches . len ( ) as u32 ) as usize ; let j = ev . 1 ; let is_start = ev . 2 >= ( matches . len ( ) as u32 ) ; if is_start { // Default case -- chain starts at this node dp [ p ] = ( k * match_score , - 1 ) ; // Find best previous chain, and extend. let best_prev = max_col_dp . get ( j as usize ) ; if best_prev . score > 0 { let prev_x = best_prev . x ; let prev_y = best_prev . y ; let cur_x = ev . 0 ; let cur_y = ev . 1 ; let gap = max ( cur_x - prev_x , cur_y - prev_y ) ; let gap_penalty = if gap > 0 { _gap_open + gap * _gap_extend } else { 0 } ; let reward = k * match_score ; let new_score = ( best_prev . score + reward ) . saturating_sub ( gap_penalty ) ; dp [ p ] = max ( dp [ p ] , ( new_score , best_prev . id as i32 ) ) ; best_dp = max ( best_dp , ( dp [ p ] . 0 , p as i32 ) ) ; } } else { // See if this kmer continues a different kmer if ev . 0 > k && ev . 1 > k { if let Ok ( cont_idx ) = matches . binary_search ( & ( ev . 0 - k - 1 , ev . 1 - k - 1 ) ) { let prev_score = dp [ cont_idx ] . 0 ; let candidate = ( prev_score + match_score , cont_idx as i32 ) ; dp [ p ] = max ( dp [ p ] , candidate ) ; best_dp = max ( best_dp , ( dp [ p ] . 0 , p as i32 ) ) ; } } let prev_frag = PrevPtr :: new ( dp [ p ] . 0 , ev . 0 , ev . 1 , p , _gap_extend ) ; max_col_dp . set ( ev . 1 as usize , prev_frag ) ; } } let mut traceback = Vec :: new ( ) ; let ( best_score , mut prev_match ) = best_dp ; while prev_match >= 0 { traceback . push ( prev_match as usize ) ; prev_match = dp [ prev_match as usize ] . 1 ; } traceback . reverse ( ) ; SparseAlignmentResult { path : traceback , score : best_score , dp_vector : dp , } }
Sparse DP routine generalizing LCSk++ method above to penalize alignment gaps.
Sparse DP routine generalizing LCSk ++ method above to penalize alignment gaps .
matches
&[(u32, u32)]
k
usize
match_score
u32
gap_open
i32
gap_extend
i32
return value
The method returns a `SparseAlignmentResult` struct with the following fields:
`path` is the SDP path, represented as vector of indices into the input matches vector.
`score` is the score of the path, which is the number of bases covered by the matched kmers.
`dp_vector` is the full DP vector, which can generally be ignored. (It may be useful for
testing purposes).
The method returns a ` SparseAlignmentResult ` struct with the following fields : ` path ` is the SDP path represented as vector of indices into the input matches vector . ` score ` is the score of the path which is the number of bases covered by the matched kmers . ` dp_vector ` is the full DP vector which can generally be ignored . ( It may be useful for testing purposes ) .