Parasail traceback info (start of query/sequence, no. of mismatches/indels)

traceback_info.c

/* Find begin of query/sequence, together with number of mismatches/indels
 * for parasail SG trace-enabled alignment.
 * 
 * Note: defines (T, STRIPED) must match the kind of alignment used.
 *
 * Compile with gcc -Wall traceback_info.c -o traceback_info -lparasail
 */

#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "parasail.h"

const char * seqA = "ACTCCTACGGGAGGCAGCAACTCCTACGGGAGGCAGCA";
const char * seqB = "GGGGGGGGGGGGGGGGACTCCTACGGAGGCACACTCCTACGGGAGGCACCGCAAGGGGGGGGGGGGGGG";

#define ALIGN parasail_sg_trace_striped_16
#define STRIPED


#define T 16
#define CONCAT3_(X, Y, Z) X##Y##Z
#define CONCAT3(X, Y, Z) CONCAT3_(X, Y, Z)
#define D CONCAT3(int, T, _t)

#if defined(STRIPED)
/* #define NAME parasail_traceback_striped_ */
#define LOC LOC_STRIPED
#else
/* #define NAME parasail_traceback_ */
#define LOC LOC_NOVEC
#endif

#define LOC_NOVEC int64_t loc = i*lenb + j;
#define LOC_STRIPED int64_t loc = j*segLen*segWidth + (i%segLen)*segWidth + (i/segLen);


/* Adapted from parasail/src/template_traceback.c */
void compute_stats(
    const char *seqA,
    int lena,
    const char *seqB,
    int lenb,
    const parasail_result_t * result,
    int64_t *start_query,
    int64_t *start_ref,
    int64_t *n_insertions,
    int64_t *n_deletions,
    int64_t *n_mismatches)
{
    int64_t i = result->end_query;
    int64_t j = result->end_ref;
    int where = PARASAIL_DIAG;
    int64_t c_ins = 0;
    int64_t c_del = 0;
    int64_t c_mm = 0;
    D *HT = (D*)result->trace->trace_table;

#if defined(STRIPED)
    int64_t segWidth = 0;
    int64_t segLen = 0;
    if (result->flag & PARASAIL_FLAG_LANES_1) {
        segWidth = 1;
    }
    if (result->flag & PARASAIL_FLAG_LANES_2) {
        segWidth = 2;
    }
    if (result->flag & PARASAIL_FLAG_LANES_4) {
        segWidth = 4;
    }
    if (result->flag & PARASAIL_FLAG_LANES_8) {
        segWidth = 8;
    }
    if (result->flag & PARASAIL_FLAG_LANES_16) {
        segWidth = 16;
    }
    if (result->flag & PARASAIL_FLAG_LANES_32) {
        segWidth = 32;
    }
    if (result->flag & PARASAIL_FLAG_LANES_64) {
        segWidth = 64;
    }
    segLen = (lena + segWidth - 1) / segWidth;
#endif
    
    while (i >= 0 || j >= 0) {
        LOC
        if (i < 0) {
            /* Added: early stopping for sg alignment */
            if (result->flag & PARASAIL_FLAG_SG) {
                break;
            }
            
            if (!(result->flag & PARASAIL_FLAG_SW)) {
                while (j >= 0) {
                    ++c_ins;
                    --j;
                }
            }
            break;
        }
        if (j < 0) {
            if (!(result->flag & PARASAIL_FLAG_SW)) {
                while (i >= 0) {
                    ++c_del;
                    --i;
                }
            }
            break;
        }
        if (PARASAIL_DIAG == where) {
            if (HT[loc] & PARASAIL_DIAG) {
                if (seqA[i] != seqB[j])
                    ++c_mm;
                
                --i;
                --j;
            }
            else if (HT[loc] & PARASAIL_INS) {
                where = PARASAIL_INS;
            }
            else if (HT[loc] & PARASAIL_DEL) {
                where = PARASAIL_DEL;
            }
            else {
                break;
            }
        }
        else if (PARASAIL_INS == where) {
            ++c_ins;
            --j;
            if (HT[loc] & PARASAIL_DIAG_E) {
                where = PARASAIL_DIAG;
            }
            else if (HT[loc] & PARASAIL_INS_E) {
                where = PARASAIL_INS;
            }
            else {
                assert(0);
            }
        }
        else if (PARASAIL_DEL == where) {
            ++c_del;
            --i;
            if (HT[loc] & PARASAIL_DIAG_F) {
                where = PARASAIL_DIAG;
            }
            else if (HT[loc] & PARASAIL_DEL_F) {
                where = PARASAIL_DEL;
            }
            else {
                assert(0);
            }
        }
        else if (PARASAIL_ZERO == where) {
            break;
        }
        else {
            assert(0);
        }        
    }

    *start_query = i + 1;
    *start_ref = j + 1;
    *n_insertions = c_ins;
    *n_deletions = c_del;
    *n_mismatches = c_mm;
}

char *substr(const char *s, int64_t start, int64_t end)
{
    char *res = malloc(end - start + 1);
    strncpy(res, s + start, end - start);
    res[end - start] = '\0';
    return res;
}



int main()
{
    parasail_result_t *result;
    const parasail_matrix_t *matrix;
    int64_t start_query, start_ref, end_query, end_ref,
        n_insertions, n_deletions, n_mismatches;
    char *partA, *partB;
    
    matrix = parasail_matrix_lookup("nuc44");
    result = ALIGN(seqA, strlen(seqA), seqB, strlen(seqB), 12, 4, matrix);

    parasail_traceback_generic(
        seqA, strlen(seqA), seqB, strlen(seqB), "Query:", "Target:",
        matrix, result, '|', '*', '*', 80, 7, 1);

    compute_stats(
        seqA, strlen(seqA), seqB, strlen(seqB), result,
        &start_query, &start_ref, &n_insertions, &n_deletions, &n_mismatches);

    end_ref = parasail_result_get_end_ref(result) + 1;
    end_query = parasail_result_get_end_query(result) + 1;
    
    printf("query runs from %lld to %lld\n", start_query, end_query);
    printf("Ref runs from %lld to %lld\n", start_ref, end_ref);
    printf("# ins/del/mm: %lld/%lld/%lld\n", n_insertions, n_deletions, n_mismatches);

    partA = substr(seqA, start_query, end_query);
    partB = substr(seqB, start_ref, end_ref);
    
    printf("Matching part of sequence: %s\n", partB);
    printf("Matching part of query: %s\n", partA);

    free(partA);
    free(partB);
    
    return 0;
}