diff --git a/src/ecoprimer.c b/src/ecoprimer.c index f3c6fd8..c6bc79c 100644 --- a/src/ecoprimer.c +++ b/src/ecoprimer.c @@ -121,13 +121,13 @@ void printapair(int32_t index,ppair_t pair, poptions_t options) wellidentifiedtaxa = (pair->intaxa + pair->outtaxa) - pair->notwellidentifiedtaxa; - printf("\t%d", pair->notwellidentifiedtaxa); + //printf("\t%d", pair->notwellidentifiedtaxa); - printf("\t%d", (pair->intaxa + pair->outtaxa)); + //printf("\t%d", (pair->intaxa + pair->outtaxa)); - //printf("\t%4.3f", (float)wellidentifiedtaxa/(options->intaxa + options->outtaxa)); + printf("\t%4.3f", (float)wellidentifiedtaxa/(options->intaxa + options->outtaxa)); printf("\t%d", pair->mind); @@ -164,6 +164,7 @@ uint32_t filterandsortpairs(ppair_t* sortedpairs,uint32_t count, poptions_t opti qfp < options->false_positive_quorum) { (void)taxonomycoverage(sortedpairs[j],options); + taxonomyspecificity(sortedpairs[j]); j++; } } diff --git a/src/libecoprimer/ecoprimer.h b/src/libecoprimer/ecoprimer.h index 5ab2c28..c65a7fd 100644 --- a/src/libecoprimer/ecoprimer.h +++ b/src/libecoprimer/ecoprimer.h @@ -218,8 +218,8 @@ typedef struct { const char *amplifia; bool_t strand; int32_t length; - void *taxontree; int32_t taxoncount; + void *taxontree; }amptotaxon_t, *pamptotaxon_t; typedef struct { diff --git a/src/libecoprimer/taxstats.c b/src/libecoprimer/taxstats.c index 1939592..83c28af 100644 --- a/src/libecoprimer/taxstats.c +++ b/src/libecoprimer/taxstats.c @@ -25,18 +25,6 @@ static int cmptaxon(const void *t1, const void* t2) return 0; } -static int cmptaxamp (const void *t1, const void* t2) -{ - const ptaxontoamp_t taxid1=(ptaxontoamp_t)t1; - const ptaxontoamp_t taxid2=(ptaxontoamp_t)t2; - - if (taxid1->taxid < taxid2->taxid) - return -1; - if (taxid1->taxid > taxid2->taxid) - return +1; - return 0; -} - int32_t counttaxon(int32_t taxid) { static void* taxontree=NULL; @@ -144,6 +132,7 @@ static int cmpamp(const void *ampf1, const void* ampf2) char cd1; char cd2; int chd = 0; + int len = 0; pamptotaxon_t pampf1 = (pamptotaxon_t) ampf1; pamptotaxon_t pampf2 = (pamptotaxon_t) ampf2; @@ -161,7 +150,9 @@ static int cmpamp(const void *ampf1, const void* ampf2) } } - for (i = 0; i < pampf1->length; i++, j += incr) + len = (pampf1->length <= pampf2->length)? pampf1->length: pampf2->length; + + for (i = 0; i < len; i++, j += incr) { cd1 = pampf1->amplifia[i]; if (incr == -1) @@ -172,6 +163,10 @@ static int cmpamp(const void *ampf1, const void* ampf2) if (cd1 < cd2) return chd ? 1: -1; if (cd2 < cd1) return chd ? -1: 1; } + + if (pampf1->length > pampf2->length) return chd ? -1: 1; + if (pampf2->length > pampf1->length) return chd ? 1: -1; + return 0; } @@ -186,11 +181,11 @@ void twalkaction (const void *node, VISIT order, int level) void taxonomyspecificity (ppair_t pair) { uint32_t i; - int32_t j; - int32_t ampfindex = 0; + uint32_t ampfindex = 0; int32_t taxid; void *ampftree = NULL; pamptotaxon_t pcurrentampf; + pamptotaxon_t *ptmp; pamptotaxon_t ampfwithtaxtree = ECOMALLOC(sizeof(amptotaxon_t) * pair->pcr.ampcount,"Cannot allocate amplifia tree"); @@ -203,23 +198,28 @@ void taxonomyspecificity (ppair_t pair) ampfwithtaxtree[ampfindex].length = pair->pcr.amplifias[i].length; pcurrentampf = &fwithtaxtree[ampfindex]; taxid = pair->pcr.amplifias[i].sequence->ranktaxonid; - - if ((pcurrentampf = tfind((void*)pcurrentampf, &ftree, cmpamp)) == NULL) + ptmp = tfind((const void*)pcurrentampf, &ftree, cmpamp); + if (ptmp == NULL) { pcurrentampf = &fwithtaxtree[ampfindex]; tsearch((void*)pcurrentampf,&ftree,cmpamp); ampfindex++; } - if (tfind((void*)((size_t)taxid), &pcurrentampf->taxontree, cmptaxon) == NULL) + else + pcurrentampf = *ptmp; + + if (tfind((void*)((size_t)taxid), &(pcurrentampf->taxontree), cmptaxon) == NULL) + { pcurrentampf->taxoncount++; - tsearch((void*)((size_t)taxid),&pcurrentampf->taxontree,cmptaxon); + tsearch((void*)((size_t)taxid),&(pcurrentampf->taxontree),cmptaxon); + } } counttaxon(-1); - for (j = 0; j < ampfindex; j++) + for (i = 0; i < ampfindex; i++) { - if (ampfwithtaxtree[j].taxoncount > 1) - twalk(ampfwithtaxtree[j].taxontree, twalkaction); + if (ampfwithtaxtree[i].taxoncount > 1) + twalk(ampfwithtaxtree[i].taxontree, twalkaction); } pair->notwellidentifiedtaxa = counttaxon(-2);