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

@@ -5,45 +5,26 @@
 	<projects>
 	</projects>
 	<buildSpec>
+		<buildCommand>
+			<name>org.python.pydev.PyDevBuilder</name>
+			<arguments>
+			</arguments>
+		</buildCommand>
 		<buildCommand>
 			<name>org.eclipse.cdt.managedbuilder.core.genmakebuilder</name>
 			<triggers>clean,full,incremental,</triggers>
 			<arguments>
-				<dictionary>
-					<key>org.eclipse.cdt.make.core.fullBuildTarget</key>
-					<value>all</value>
-				</dictionary>
 				<dictionary>
 					<key>?name?</key>
 					<value></value>
 				</dictionary>
-				<dictionary>
-					<key>org.eclipse.cdt.make.core.enableAutoBuild</key>
-					<value>false</value>
-				</dictionary>
-				<dictionary>
-					<key>org.eclipse.cdt.make.core.enableFullBuild</key>
-					<value>true</value>
-				</dictionary>
-				<dictionary>
-					<key>org.eclipse.cdt.make.core.enableCleanBuild</key>
-					<value>true</value>
-				</dictionary>
-				<dictionary>
-					<key>org.eclipse.cdt.make.core.cleanBuildTarget</key>
-					<value>clean</value>
-				</dictionary>
 				<dictionary>
 					<key>org.eclipse.cdt.make.core.append_environment</key>
 					<value>true</value>
 				</dictionary>
 				<dictionary>
-					<key>org.eclipse.cdt.make.core.contents</key>
+					<key>org.eclipse.cdt.make.core.autoBuildTarget</key>
-					<value>org.eclipse.cdt.make.core.activeConfigSettings</value>
+					<value>all</value>
-				</dictionary>
-				<dictionary>
-					<key>org.eclipse.cdt.make.core.useDefaultBuildCmd</key>
-					<value>true</value>
 				</dictionary>
 				<dictionary>
 					<key>org.eclipse.cdt.make.core.buildArguments</key>
@@ -54,13 +35,37 @@
 					<value>make</value>
 				</dictionary>
 				<dictionary>
-					<key>org.eclipse.cdt.make.core.autoBuildTarget</key>
+					<key>org.eclipse.cdt.make.core.cleanBuildTarget</key>
+					<value>clean</value>
+				</dictionary>
+				<dictionary>
+					<key>org.eclipse.cdt.make.core.contents</key>
+					<value>org.eclipse.cdt.make.core.activeConfigSettings</value>
+				</dictionary>
+				<dictionary>
+					<key>org.eclipse.cdt.make.core.enableAutoBuild</key>
+					<value>false</value>
+				</dictionary>
+				<dictionary>
+					<key>org.eclipse.cdt.make.core.enableCleanBuild</key>
+					<value>true</value>
+				</dictionary>
+				<dictionary>
+					<key>org.eclipse.cdt.make.core.enableFullBuild</key>
+					<value>true</value>
+				</dictionary>
+				<dictionary>
+					<key>org.eclipse.cdt.make.core.fullBuildTarget</key>
 					<value>all</value>
 				</dictionary>
 				<dictionary>
 					<key>org.eclipse.cdt.make.core.stopOnError</key>
 					<value>true</value>
 				</dictionary>
+				<dictionary>
+					<key>org.eclipse.cdt.make.core.useDefaultBuildCmd</key>
+					<value>true</value>
+				</dictionary>
 			</arguments>
 		</buildCommand>
 		<buildCommand>
@@ -73,5 +78,6 @@
 		<nature>org.eclipse.cdt.core.cnature</nature>
 		<nature>org.eclipse.cdt.managedbuilder.core.ScannerConfigNature</nature>
 		<nature>org.eclipse.cdt.managedbuilder.core.managedBuildNature</nature>
+		<nature>org.python.pydev.pythonNature</nature>
 	</natures>
 </projectDescription>

.pydevproject

@@ -0,0 +1,7 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<?eclipse-pydev version="1.0"?>
+
+<pydev_project>
+<pydev_property name="org.python.pydev.PYTHON_PROJECT_INTERPRETER">Default</pydev_property>
+<pydev_property name="org.python.pydev.PYTHON_PROJECT_VERSION">python 2.6</pydev_property>
+</pydev_project>

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);
 

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;

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].primer = primers->primers+i;
-				matches[j].strand=TRUE;
+			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].primer = primers->primers+i;
-				matches[j].strand=FALSE;
+			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
+		    	   //else
-		    		   pcurrent->outexample++;
+		    		 //  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;