diff --git a/src/libecoprimer/strictprimers.c b/src/libecoprimer/strictprimers.c index ad78b08..e6e7f4e 100644 --- a/src/libecoprimer/strictprimers.c +++ b/src/libecoprimer/strictprimers.c @@ -132,7 +132,7 @@ pwordcount_t lookforStrictPrimer(pecodnadb_t database, uint32_t seqdbsize, double seconde; char *logfilename; FILE *logfile; - uint32_t i; + uint32_t i, j; bool_t first=TRUE; pwordcount_t strictprimers=NULL; uint64_t totallength=0; @@ -232,6 +232,36 @@ pwordcount_t lookforStrictPrimer(pecodnadb_t database, uint32_t seqdbsize, if (neededWords) ECOFREE(neededWords,"Clean needed word table"); + //TR: Somehow for some primers strictcount value is extremely large hence invalid + //we need to remove these primers from the list + j = strictprimers->size+1; + for (i=0; isize; i++) + { + if (strictprimers->strictcount[i] > seqdbsize) + { + if (j == (strictprimers->size+1)) + j = i; + } + + if (j < i && strictprimers->strictcount[i] <= seqdbsize) + { + strictprimers->words[j] = strictprimers->words[i]; + strictprimers->strictcount[j] = strictprimers->strictcount[i]; + j++; + } + } + + if (j < strictprimers->size) + { + strictprimers->size = j; + strictprimers->strictcount = ECOREALLOC(strictprimers->strictcount, + sizeof(uint32_t)*strictprimers->size, + "Cannot reallocate strict primer count table"); + strictprimers->words = ECOREALLOC(strictprimers->words, + sizeof(word_t)*strictprimers->size, + "Cannot reallocate strict primer table"); + } + return strictprimers; } diff --git a/src/libecoprimer/taxstats.c b/src/libecoprimer/taxstats.c index 5b2d64e..cdc97a5 100644 --- a/src/libecoprimer/taxstats.c +++ b/src/libecoprimer/taxstats.c @@ -126,7 +126,7 @@ float taxonomycoverage(ppair_t pair, poptions_t options) return pair->bc; } - +/* static int cmpamp(const void *ampf1, const void* ampf2) { int i; @@ -170,6 +170,55 @@ static int cmpamp(const void *ampf1, const void* ampf2) if (pampf1->length > pampf2->length) return chd ? -1: 1; if (pampf2->length > pampf1->length) return chd ? 1: -1; + return 0; +}*/ + +static int cmpamp(const void *ampf1, const void* ampf2) +{ + int i; + char cd1; + char cd2; + int len = 0; + char *ch1; + char *ch2; + int incr1; + int incr2; + + pamptotaxon_t pampf1 = (pamptotaxon_t) ampf1; + pamptotaxon_t pampf2 = (pamptotaxon_t) ampf2; + + ch1 = pampf1->amplifia; + ch2 = pampf2->amplifia; + + incr1 = 1; + incr2 = 1; + + if (!pampf1->strand) + incr1 = -1; + if (!pampf2->strand) + incr2 = -1; + + len = (pampf1->length <= pampf2->length)? pampf1->length: pampf2->length; + for (i = 0; i < len; i++) + { + cd1 = *ch1; + if (incr1 == -1) + cd1 = ecoComplementChar(*ch1); + + cd2 = *ch2; + if (incr2 == -1) + cd2 = ecoComplementChar(*ch2); + + if (cd1 < cd2) return -1; + if (cd2 < cd1) return 1; + + ch1 += incr1; + ch2 += incr2; + } + + if (pampf1->length > pampf2->length) return 1; + if (pampf2->length > pampf1->length) return -1; + return 0; }