From 3a617890cadf01eaff845738b57e28ad27590bec Mon Sep 17 00:00:00 2001 From: Eric Coissac Date: Tue, 23 Jun 2009 14:11:39 +0000 Subject: [PATCH] New option for reference sequence and bug correction for insequence count git-svn-id: https://www.grenoble.prabi.fr/svn/LECASofts/ecoPrimers/trunk@218 60f365c0-8329-0410-b2a4-ec073aeeaa1d --- .project | 60 ++++++++++++++++++++---------------- .pydevproject | 7 +++++ src/ecoprimer.c | 45 +++++++++++++++++++++++++-- src/libecoprimer/ecoprimer.h | 7 +++++ src/libecoprimer/pairs.c | 42 +++++++++++++++++++------ 5 files changed, 123 insertions(+), 38 deletions(-) create mode 100644 .pydevproject diff --git a/.project b/.project index ffa37e5..50a83f9 100644 --- a/.project +++ b/.project @@ -5,45 +5,26 @@ + + org.python.pydev.PyDevBuilder + + + org.eclipse.cdt.managedbuilder.core.genmakebuilder clean,full,incremental, - - org.eclipse.cdt.make.core.fullBuildTarget - all - ?name? - - org.eclipse.cdt.make.core.enableAutoBuild - false - - - org.eclipse.cdt.make.core.enableFullBuild - true - - - org.eclipse.cdt.make.core.enableCleanBuild - true - - - org.eclipse.cdt.make.core.cleanBuildTarget - clean - org.eclipse.cdt.make.core.append_environment true - org.eclipse.cdt.make.core.contents - org.eclipse.cdt.make.core.activeConfigSettings - - - org.eclipse.cdt.make.core.useDefaultBuildCmd - true + org.eclipse.cdt.make.core.autoBuildTarget + all org.eclipse.cdt.make.core.buildArguments @@ -54,13 +35,37 @@ make - org.eclipse.cdt.make.core.autoBuildTarget + org.eclipse.cdt.make.core.cleanBuildTarget + clean + + + org.eclipse.cdt.make.core.contents + org.eclipse.cdt.make.core.activeConfigSettings + + + org.eclipse.cdt.make.core.enableAutoBuild + false + + + org.eclipse.cdt.make.core.enableCleanBuild + true + + + org.eclipse.cdt.make.core.enableFullBuild + true + + + org.eclipse.cdt.make.core.fullBuildTarget all org.eclipse.cdt.make.core.stopOnError true + + org.eclipse.cdt.make.core.useDefaultBuildCmd + true + @@ -73,5 +78,6 @@ org.eclipse.cdt.core.cnature org.eclipse.cdt.managedbuilder.core.ScannerConfigNature org.eclipse.cdt.managedbuilder.core.managedBuildNature + org.python.pydev.pythonNature diff --git a/.pydevproject b/.pydevproject new file mode 100644 index 0000000..6725cde --- /dev/null +++ b/.pydevproject @@ -0,0 +1,7 @@ + + + + +Default +python 2.6 + diff --git a/src/ecoprimer.c b/src/ecoprimer.c index 3483031..98e4f76 100644 --- a/src/ecoprimer.c +++ b/src/ecoprimer.c @@ -106,6 +106,8 @@ void initoptions(poptions_t options) options->restricted_taxid=NULL; //**< limit amplification below these taxid options->ignored_taxid=NULL; //**< no amplification below these taxid options->prefix=NULL; + options->reference=NULL; + options->refseq=NULL; options->circular=0; options->doublestrand=1; options->strict_quorum=0.7; @@ -204,7 +206,29 @@ void printapair(int32_t index,ppair_t pair, poptions_t options) printf("\t%d", pair->mind); printf("\t%d", pair->maxd); - printf("\t%3.2f\n", (float)pair->sumd/pair->inexample); + printf("\t%3.2f", (float)pair->sumd/pair->inexample); + + if (options->refseq && pair->refsequence >=0) + { + printf("\t%s:",options->reference); + + + if (pair->pcr.amplifias[pair->refsequence].strand) + printf("join("); + else + printf("complement("); + + printf("%d..%d,%d..%d",pair->pcr.amplifias[pair->refsequence].begin - options->primer_length + 1, + pair->pcr.amplifias[pair->refsequence].begin, + pair->pcr.amplifias[pair->refsequence].end + 2, + pair->pcr.amplifias[pair->refsequence].end + options->primer_length + 1 + ); + printf(")"); + printf("\t"); + + } + + printf("\n"); } @@ -444,7 +468,7 @@ int main(int argc, char **argv) initoptions(&options); - while ((carg = getopt(argc, argv, "hfvcUDSd:l:L:e:i:r:q:3:s:x:t:O:")) != -1) { + while ((carg = getopt(argc, argv, "hfvcUDSd:l:L:e:i:r:R:q:3:s:x:t:O:")) != -1) { switch (carg) { /* ---------------------------- */ @@ -549,6 +573,14 @@ int main(int argc, char **argv) "Error on restricted_taxid reallocation"); sscanf(optarg,"%d",&(options.restricted_taxid[options.r])); options.r++; + break; + + /* -------------------- */ + case 'R': /* reference sequence */ + /* -------------------- */ + options.reference = ECOMALLOC(strlen(optarg)+1, + "Error on prefix allocation"); + strcpy(options.reference,optarg); break; /* --------------------------------- */ @@ -589,6 +621,15 @@ int main(int argc, char **argv) seqdb = readdnadb(options.prefix,&seqdbsize); + if (options.reference) + for (i=0; i < seqdbsize;i++) + if (strcmp(seqdb[i]->AC,options.reference)==0) + { + options.refseq=seqdb[i]; + options.refseqid=i; + fprintf(stderr,"Reference sequence %s identified\n",options.reference); + } + fprintf(stderr,"Ok\n"); fprintf(stderr,"Sequence read : %d\n",(int32_t)seqdbsize); diff --git a/src/libecoprimer/ecoprimer.h b/src/libecoprimer/ecoprimer.h index b984e76..0cf3198 100644 --- a/src/libecoprimer/ecoprimer.h +++ b/src/libecoprimer/ecoprimer.h @@ -133,6 +133,8 @@ typedef struct { bool_t strand; const char *amplifia; int32_t length; + uint32_t begin; + uint32_t end; } amplifia_t, *pamplifia_t; typedef struct { @@ -180,12 +182,14 @@ typedef struct { float quorumout; float bs; float bc; + int32_t refsequence; // // uint32_t taxsetcount; // uint32_t taxsetindex; // ptaxampset_t taxset; // // uint32_t oktaxoncount; + uint32_t curseqid; } pair_t, *ppair_t; @@ -247,6 +251,9 @@ typedef struct { int32_t *restricted_taxid; //**< limit amplification below these taxid int32_t *ignored_taxid; //**< no amplification below these taxid char *prefix; + char *reference; + pecoseq_t refseq; + uint32_t refseqid; uint32_t circular; uint32_t doublestrand; float strict_quorum; diff --git a/src/libecoprimer/pairs.c b/src/libecoprimer/pairs.c index 88e6c21..8f2df4c 100644 --- a/src/libecoprimer/pairs.c +++ b/src/libecoprimer/pairs.c @@ -205,8 +205,8 @@ static void buildPrimerPairsForOneSeq(uint32_t seqid, { if (primers->primers[i].directCount[seqid]==1) { - matches[j].primer = primers->primers+i; - matches[j].strand=TRUE; + matches[j].primer = primers->primers+i; + matches[j].strand=TRUE; matches[j].position=primers->primers[i].directPos[seqid].value; j++; } @@ -222,8 +222,8 @@ static void buildPrimerPairsForOneSeq(uint32_t seqid, { if (primers->primers[i].reverseCount[seqid]==1) { - matches[j].primer = primers->primers+i; - matches[j].strand=FALSE; + matches[j].primer = primers->primers+i; + matches[j].strand=FALSE; matches[j].position=primers->primers[i].reversePos[seqid].value; j++; } @@ -270,6 +270,8 @@ static void buildPrimerPairsForOneSeq(uint32_t seqid, current.sumd=0; current.inexample=0; current.outexample=0; + current.curseqid = 0; + current.refsequence=-1; // Standardize the pair @@ -285,7 +287,7 @@ static void buildPrimerPairsForOneSeq(uint32_t seqid, current.asdirect1=current.asdirect2; current.asdirect2=bswp; } - + // Look for the new pair in already seen pairs @@ -295,7 +297,7 @@ static void buildPrimerPairsForOneSeq(uint32_t seqid, if (seqdb[seqid]->isexample) { - pcurrent->inexample++; + //pcurrent->inexample++; pcurrent->sumd+=distance; if ((pcurrent->maxd==DMAX) || (distance > pcurrent->maxd)) @@ -304,11 +306,29 @@ static void buildPrimerPairsForOneSeq(uint32_t seqid, if (distance < pcurrent->mind) pcurrent->mind = distance; } - else - pcurrent->outexample++; + //else + // pcurrent->outexample++; - if ((pcurrent->outexample+pcurrent->inexample)==1) + if (pcurrent->curseqid != (seqid+1)) + { + if (seqdb[seqid]->isexample) + pcurrent->inexample++; + else + pcurrent->outexample++; + + if (pcurrent->curseqid != 0) + pcurrent->curseqid = seqid+1; + } + + /*if ((pcurrent->outexample+pcurrent->inexample)==0) + { + fprintf(stderr,"pcurrent->outexample+pcurrent->inexample=0!\n"); + exit(0); + }*/ + + if (pcurrent->curseqid == 0)//((pcurrent->outexample+pcurrent->inexample)==1) { + pcurrent->curseqid = seqid+1; paircount++; pcurrent->pcr.ampslot=200; pcurrent->pcr.ampcount=0; @@ -326,9 +346,13 @@ static void buildPrimerPairsForOneSeq(uint32_t seqid, } } + if (seqid==options->refseqid) + pcurrent->refsequence=seqid; pcurrent->pcr.amplifias[pcurrent->pcr.ampcount].length=distance; pcurrent->pcr.amplifias[pcurrent->pcr.ampcount].sequence=seqdb[seqid]; pcurrent->pcr.amplifias[pcurrent->pcr.ampcount].strand=strand; + pcurrent->pcr.amplifias[pcurrent->pcr.ampcount].begin=matches[i].position + options->primer_length; + pcurrent->pcr.amplifias[pcurrent->pcr.ampcount].end= matches[j].position - 1; if (strand) pcurrent->pcr.amplifias[pcurrent->pcr.ampcount].amplifia= seqdb[seqid]->SQ + matches[i].position + options->primer_length;