biohazard-2.1: cbits/trim.c

#include <stdint.h>
#include <limits.h>

#define MIN(a,b) ((a)<(b)?(a):(b))
#define MAX(a,b) ((a)>(b)?(a):(b))

/* This version receives the second sequence reverse-complemented and
 * the second quality vector reversed.  Allows for a forward loop,
 * doesn't need the lookup table. */
int prim_match_reads( int i1       // matches from i2
                    , int i2       // ???
                    , int r        // length to match
                    , const uint8_t *rd1
                    , const uint8_t *qs1
                    , const uint8_t *rc_rd2
                    , const uint8_t *rv_qs2 )
{
    int acc = 0 ;
    while( r != 0 )
    {
        uint8_t n1 = rd1[ i1 ] ;
        uint8_t n2 = rc_rd2[ i2 ] ;
        uint8_t q1 = qs1[ i1 ] ;
        uint8_t q2 = rv_qs2[ i2 ] ;

        acc += (n1 & 0xF) == (n2 & 0xF) ? 0 : 5 + (q1 < q2 ? q1 : q2) ;

        ++i1 ;
        ++i2 ;
        --r ;
    }
    return acc ;
}

int prim_match_ad( int off
                 , int i
                 , const uint8_t *rd
                 , const uint8_t *qs
                 , const uint8_t *ad )
{
    int acc = 0 ;
    while( i > 0 )
    {
        --i;
        acc += rd[ i+off ] == ad[ i ] ? 0 : 5 +
               (qs[ i+off ] < 25 ? qs[ i+off ] : 25) ;
    }
    return acc ;
}

// return the length that minimizes the score
int prim_merge_score (
        uint8_t **p_fwd_ads,
        int *p_fwd_lns,
        int n_fwd_adapters,

        uint8_t **p_rev_ads,
        int *p_rev_lns,
        int n_rev_adapters,

        uint8_t *p_rd1,
        uint8_t *p_qs1,
        int l1,
        uint8_t *p_rd2,
        uint8_t *p_qs2,
        int l2,
        uint8_t *p_rrd2,
        uint8_t *p_rqs2,

        int *scores )
{
    int minl = -1 ;
    int min0 = 6 * (l1+l2) ;
    int min1= INT_MAX ;

    for( int l = 0 ; l <= l1+l2 ; ++l )
    {
        int ff = INT_MAX ;
        for( int i = 0 ; i != n_fwd_adapters ; ++i )
        {
            int efflength = MIN(l1 - l, p_fwd_lns[i]) ;
            int x = prim_match_ad( l, efflength, p_rd1, p_qs1, p_fwd_ads[i] ) ;
            ff = MIN(ff, x + 6 * MAX( 0, l1 - p_fwd_lns[i] - l)) ;
        }

        int rr = INT_MAX ;
        for( int i = 0 ; i != n_rev_adapters ; ++i )
        {
            int efflength = MIN(l2 - l,p_rev_lns[i]) ;
            int x = prim_match_ad( l, efflength, p_rd2, p_qs2, p_rev_ads[i] ) ;
            rr = MIN(rr, x + 6 * MAX( 0, l1 - p_rev_lns[i] - l)) ;
        }

        int minidx1 = MAX(l - l2, 0) ; // vec1, forward
        int minidx2 = MAX(l2 - l, 0) ; // vec2, forward-on-reverse
        int efflength = MIN(l1 + l2 - l, l) ; // effective length

        int mm = efflength <= 0 ? 0 : prim_match_reads( minidx1, minidx2, efflength, p_rd1, p_qs1, p_rrd2, p_rqs2 ) ;
        int z = 6 * l + ff + rr + mm ;

        if( z < min0 ) {
            min1 = min0 ;
            min0 = z ;
            minl = l ;
        } else if( z < min1 ) {
            min1 = z ;
        }
    }
    scores[0] = min1 - min0 ;
    scores[1] = 6*(l1+l2) - min0 ;
    return minl ;
}