obitools/ecopcr
2
0
Fork 0
Code Issues 4 Pull Requests Actions Packages Projects Releases 3 Wiki Activity

add an hidden -D option

Browse Source 
git-svn-id: https://www.grenoble.prabi.fr/svn/LECASofts/ecoPCR/trunk@296 60f365c0-8329-0410-b2a4-ec073aeeaa1d
This commit is contained in:
Eric Coissac
2011-05-25 15:41:22 +00:00
parent 87c2496447
commit 01173d22cf
1 changed files with 49 additions and 14 deletions
Download Patch File Download Diff File Expand all files Collapse all files

63
src/ecopcr.c
View File

@ -114,7 +114,8 @@ void printRepeat(ecoseq_t *seq,
char kingdom, char kingdom,
int32_t pos1, int32_t pos2, int32_t pos1, int32_t pos2,
int32_t err1, int32_t err2, int32_t err1, int32_t err2,
ecotaxonomy_t *taxonomy) ecotaxonomy_t *taxonomy,
int32_t delta)
{ {
char *AC; char *AC;
int32_t seqlength; int32_t seqlength;
@ -137,12 +138,15 @@ void printRepeat(ecoseq_t *seq,
int32_t error1; int32_t error1;
int32_t error2; int32_t error2;
int32_t ldeta,rdelta;
char *amplifia = NULL; char *amplifia = NULL;
int32_t amplength; int32_t amplength;
double tm1,tm2; double tm1,tm2;
double tm=0; double tm=0;
int32_t i;
AC = seq->AC; AC = seq->AC;
seqlength = seq->SQ_length; seqlength = seq->SQ_length;
@ -200,41 +204,65 @@ void printRepeat(ecoseq_t *seq,
superkingdom_name = "###"; superkingdom_name = "###";
} }
ldeta=(pos1 <= delta)?pos1:delta;
amplifia = getSubSequence(seq->SQ,pos1,pos2); rdelta(pos2+delta>=seqlength)?seqlength-pos2-1:delta;
amplength= strlen(amplifia);
amplifia = getSubSequence(seq->SQ,pos1-ldelta,pos2+rdelta);
amplength= strlen(amplifia)-rdelta-ldelta;
if (strand=='R') if (strand=='R')
{ {
ecoComplementSequence(amplifia); ecoComplementSequence(amplifia);
strncpy(oligo1,amplifia,o2->patlen); strncpy(oligo1,amplifia+ rdelta ,o2->patlen);
oligo1[o2->patlen]=0; oligo1[o2->patlen]=0;
error1=err2; error1=err2;
strncpy(oligo2,amplifia + amplength - o1->patlen,o1->patlen); strncpy(oligo2,amplifia + ldelta + amplength - o1->patlen,o1->patlen);
oligo2[o1->patlen]=0; oligo2[o1->patlen]=0;
error2=err1; error2=err1;
amplifia+=o2->patlen; if (delta==0)
amplifia+=o2->patlen;
else
{
delta=ldelta;
ldelta=rdelta+o2->patlen;
rdelta=delta+o1->patlen;
}
} }
else else
{ {
strncpy(oligo1,amplifia,o1->patlen); strncpy(oligo1,amplifia+ldelta,o1->patlen);
oligo1[o1->patlen]=0; oligo1[o1->patlen]=0;
error1=err1; error1=err1;
strncpy(oligo2,amplifia + amplength - o2->patlen,o2->patlen); strncpy(oligo2,amplifia + rdelta + amplength - o2->patlen,o2->patlen);
oligo2[o2->patlen]=0; oligo2[o2->patlen]=0;
error2=err2; error2=err2;
amplifia+=o1->patlen; if (delta==0)
amplifia+=o1->patlen;
else
{
ldelta+=o1->patlen;
rdelta+=o2->patlen;
}
} }
ecoComplementSequence(oligo2); 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<ldelta;i++)
amplifia[i]|=32;
for (i=0;i<rdelta;i++)
amplifia[delta-i]|=32;
}
if (compute_tm) if (compute_tm)
{ {
@ -309,6 +337,7 @@ int main(int argc, char **argv)
PatternPtr o1c; PatternPtr o1c;
PatternPtr o2c; PatternPtr o2c;
int32_t delta=0;
int32_t lmin=0; int32_t lmin=0;
int32_t lmax=0; int32_t lmax=0;
int32_t error_max=0; int32_t error_max=0;
@ -350,7 +379,7 @@ int main(int argc, char **argv)
double salt=0.05; double salt=0.05;
CNNParams tparm; CNNParams tparm;
while ((carg = getopt(argc, argv, "hcd:l:L:e:i:r:km:a:t")) != -1) { while ((carg = getopt(argc, argv, "hcd:l:L:e:i:r:km:a:tD:")) != -1) {
switch (carg) { switch (carg) {
/* -------------------- */ /* -------------------- */
@ -368,6 +397,12 @@ int main(int argc, char **argv)
exit(0); exit(0);
break; break;
/* ------------------------- */
case 'D': /* min amplification lenght */
/* ------------------------- */
sscanf(optarg,"%d",&delta);
break;
/* ------------------------- */ /* ------------------------- */
case 'l': /* min amplification lenght */ case 'l': /* min amplification lenght */
/* ------------------------- */ /* ------------------------- */
@ -607,7 +642,7 @@ int main(int argc, char **argv)
if (length && if (length &&
(!lmin || (length >= lmin)) && (!lmin || (length >= lmin)) &&
(!lmax || (length <= lmax))) (!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); //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 && if (length &&
(!lmin || (length >= lmin)) && (!lmin || (length >= lmin)) &&
(!lmax || (length <= lmax))) (!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); //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);
} }
} }