vcmp.c

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <htslib/hts.h>
#include "vcmp.h"

struct _vcmp_t
{
    char *dref;
    int ndref, mdref;   // ndref: positive when ref1 longer, negative when ref2 is longer
    int nmatch;
};

vcmp_t *vcmp_init()
{
    return (vcmp_t*)calloc(1,sizeof(vcmp_t));
}

void vcmp_destroy(vcmp_t *vcmp)
{
    free(vcmp->dref);
    free(vcmp);
}

int vcmp_set_ref(vcmp_t *vcmp, char *ref1, char *ref2)
{
    vcmp->ndref = 0;

    char *a = ref1, *b = ref2;
    while ( *a && *b && *a==*b ) { a++; b++; }
    if ( !*a && !*b ) return 0;
    if ( *a && *b ) return -1;  // refs not compatible

    if ( *a )   // ref1 is longer
    {
        vcmp->nmatch = b-ref2;
        while ( *a ) a++;
        vcmp->ndref = (a-ref1) - vcmp->nmatch;
        hts_expand(char,vcmp->ndref+1,vcmp->mdref,vcmp->dref);
        memcpy(vcmp->dref,ref1+vcmp->nmatch,vcmp->ndref);
        vcmp->dref[vcmp->ndref] = 0;
        return 0;
    }

    // ref2 is longer
    vcmp->nmatch = a-ref1;
    while ( *b ) b++;
    vcmp->ndref = (b-ref2) - vcmp->nmatch;
    hts_expand(char,vcmp->ndref+1,vcmp->mdref,vcmp->dref);
    memcpy(vcmp->dref,ref2+vcmp->nmatch,vcmp->ndref);
    vcmp->dref[vcmp->ndref] = 0;
    vcmp->ndref *= -1;
    return 0;
}

int vcmp_find_allele(vcmp_t *vcmp, char **als1, int nals1, char *al2)
{
    int i, j;
    for (i=0; i<nals1; i++)
    {
        char *a = als1[i], *b = al2;
        while ( *a && *b && *a==*b ) { a++; b++; }
        if ( *a && *b ) continue;   // mismatch
        if ( !vcmp->ndref )
        {
            if ( !*a && !*b ) break;    // found
            continue;
        }

        // the prefixes match
        if ( *a )
        {
            if ( vcmp->ndref<0 ) continue;
            for (j=0; j<vcmp->ndref; j++)
                if ( !a[j] || a[j]!=vcmp->dref[j] ) break;
            if ( j!=vcmp->ndref || a[j] ) continue;
            break;  // found
        }
        
        if ( vcmp->ndref>0 ) continue;
        for (j=0; j<-vcmp->ndref; j++)
            if ( !b[j] || b[j]!=vcmp->dref[j] ) break;
        if ( j!=-vcmp->ndref || b[j] ) continue;
        break;  // found
    }
    if (i==nals1) return -1;
    return i;
}