diff --git a/src/ecopcr.c b/src/ecopcr.c index 671d4d5..8892e1b 100644 --- a/src/ecopcr.c +++ b/src/ecopcr.c @@ -114,7 +114,8 @@ void printRepeat(ecoseq_t *seq, char kingdom, int32_t pos1, int32_t pos2, int32_t err1, int32_t err2, - ecotaxonomy_t *taxonomy) + ecotaxonomy_t *taxonomy, + int32_t delta) { char *AC; int32_t seqlength; @@ -137,12 +138,15 @@ void printRepeat(ecoseq_t *seq, int32_t error1; int32_t error2; + int32_t ldeta,rdelta; char *amplifia = NULL; int32_t amplength; double tm1,tm2; double tm=0; + int32_t i; + AC = seq->AC; seqlength = seq->SQ_length; @@ -200,41 +204,65 @@ void printRepeat(ecoseq_t *seq, superkingdom_name = "###"; } - - amplifia = getSubSequence(seq->SQ,pos1,pos2); - amplength= strlen(amplifia); + ldeta=(pos1 <= delta)?pos1:delta; + rdelta(pos2+delta>=seqlength)?seqlength-pos2-1:delta; + + amplifia = getSubSequence(seq->SQ,pos1-ldelta,pos2+rdelta); + amplength= strlen(amplifia)-rdelta-ldelta; if (strand=='R') { ecoComplementSequence(amplifia); - strncpy(oligo1,amplifia,o2->patlen); + strncpy(oligo1,amplifia+ rdelta ,o2->patlen); oligo1[o2->patlen]=0; error1=err2; - strncpy(oligo2,amplifia + amplength - o1->patlen,o1->patlen); + strncpy(oligo2,amplifia + ldelta + amplength - o1->patlen,o1->patlen); oligo2[o1->patlen]=0; error2=err1; - amplifia+=o2->patlen; + if (delta==0) + amplifia+=o2->patlen; + else + { + delta=ldelta; + ldelta=rdelta+o2->patlen; + rdelta=delta+o1->patlen; + } } else { - strncpy(oligo1,amplifia,o1->patlen); + strncpy(oligo1,amplifia+ldelta,o1->patlen); oligo1[o1->patlen]=0; error1=err1; - strncpy(oligo2,amplifia + amplength - o2->patlen,o2->patlen); + strncpy(oligo2,amplifia + rdelta + amplength - o2->patlen,o2->patlen); oligo2[o2->patlen]=0; error2=err2; - amplifia+=o1->patlen; + if (delta==0) + amplifia+=o1->patlen; + else + { + ldelta+=o1->patlen; + rdelta+=o2->patlen; + } } ecoComplementSequence(oligo2); - amplifia[amplength - o2->patlen - o1->patlen]=0; + if(delta==0) + amplifia[amplength - o2->patlen - o1->patlen]=0; + else + { + delta=ldelta+rdelta+amplength-1; + for (i=0;i= lmin)) && (!lmax || (length <= lmax))) - printRepeat(seq,oligo1,oligo2,compute_tm,&tparm,o1,o2c,'D',kingdom_mode,posi,posj,erri,errj,taxonomy); + printRepeat(seq,oligo1,oligo2,compute_tm,&tparm,o1,o2c,'D',kingdom_mode,posi,posj,erri,errj,taxonomy,delta); //printf("%s\tD\t%s...%s (%d)\t%d\t%d\t%d\t%d\t%s\n",seq->AC,head,tail,seq->SQ_length,o1Hits,o2cHits,posi,posj,scname); } @@ -661,7 +696,7 @@ int main(int argc, char **argv) if (length && (!lmin || (length >= lmin)) && (!lmax || (length <= lmax))) - printRepeat(seq,oligo1,oligo2,compute_tm,&tparm,o2,o1c,'R',kingdom_mode,posi,posj,erri,errj,taxonomy); + printRepeat(seq,oligo1,oligo2,compute_tm,&tparm,o2,o1c,'R',kingdom_mode,posi,posj,erri,errj,taxonomy,delta); //printf("%s\tR\t%s...%s (%d)\t%d\t%d\t%d\t%d\t%s\n",seq->AC,head,tail,seq->SQ_length,o2Hits,o1cHits,posi,posj,scname); } }