From 5d8103f4b4a3c0a993e3ab8c84696cc3fdf26950 Mon Sep 17 00:00:00 2001
From: Eric Coissac <eric.coissac@metabarcoding.org>
Date: Mon, 13 Dec 2010 17:09:54 +0000
Subject: [PATCH] added check for filtering pairs having specificity below a
 given threshold given using -T option

git-svn-id: https://www.grenoble.prabi.fr/svn/LECASofts/ecoPrimers/trunk@288 60f365c0-8329-0410-b2a4-ec073aeeaa1d
---
 src/ecoprimer.c              | 20 ++++++++++++++++----
 src/libecoprimer/ecoprimer.h |  1 +
 2 files changed, 17 insertions(+), 4 deletions(-)

diff --git a/src/ecoprimer.c b/src/ecoprimer.c
index 88662bf..b159c5e 100644
--- a/src/ecoprimer.c
+++ b/src/ecoprimer.c
@@ -76,6 +76,7 @@ static void PrintHelp()
               PP      "-f    : Remove data mining step during  strict primer identification\n\n");
               PP      "-v    : Store statistic file about memory usage during strict primer identification\n\n");
               PP      "-p    : Print sets of primers\n\n");
+			  PP      "-T    : Ignore pairs having specificity below this Threshold\n\n");
               PP      "\n");
               PP      "------------------------------------------\n");
               PP      "Table result description : \n");
@@ -149,6 +150,7 @@ void initoptions(poptions_t options)
 	options->salt = DEF_SALT;
 	options->printAC=FALSE;
 	options->print_sets_of_primers = FALSE;
+	options->specificity_threshold = 0.6;
 }
 
 void printapair(int32_t index,ppair_t pair, poptions_t options)
@@ -345,7 +347,11 @@ uint32_t filterandsortpairs(ppair_t* sortedpairs,uint32_t count, poptions_t opti
 			sortedpairs[j]->wellIdentifiedSeqs = ECOMALLOC(options->dbsize * sizeof(int),"Cannot allocate well_identified_array");
 			(void)taxonomycoverage(sortedpairs[j],options);
 			taxonomyspecificity(sortedpairs[j], seqdb, options->dbsize);
-			j++;
+			//j++;
+			//if specificity less than user provieded threshold (default 60%) then ignore this pair
+			if (sortedpairs[j]->bs >= options->specificity_threshold)
+				j++;
+			
 		}
 
 	}
@@ -443,6 +449,7 @@ void printpairs (ppairtree_t pairs, poptions_t options,ecotaxonomy_t *taxonomy,
 		printf("# DB sequences are considered as circular\n");
 	else
 		printf("# DB sequences are considered as linear\n");
+	printf("# Pairs having specificity less than %0.2f will be ignored\n", options->specificity_threshold);
 	printf("#\n");
 
    for (i=0;i < count;i++)
@@ -538,7 +545,7 @@ int main(int argc, char **argv)
 
 	initoptions(&options);
 
-    while ((carg = getopt(argc, argv, "hAfvcUDSpE:d:l:L:e:i:r:R:q:3:s:x:t:O:m:a:")) != -1) {
+    while ((carg = getopt(argc, argv, "hAfvcUDSpE:d:l:L:e:i:r:R:q:3:s:x:t:O:m:a:T:")) != -1) {
 
      switch (carg) {
 								 /* ---------------------------- */
@@ -706,7 +713,13 @@ int main(int argc, char **argv)
 								/* --------------------------------- */
         	options.print_sets_of_primers = TRUE;
         break;
-
+			 
+								/* --------------------------------- */
+		case 'T':               /* Ignore pairs having specificity below this Threshold     */
+								/* --------------------------------- */
+			 sscanf(optarg,"%f",&(options.specificity_threshold));
+		break;
+			 
         case '?':               /* bad option           */
                                 /* -------------------- */
             errflag++;
@@ -717,7 +730,6 @@ int main(int argc, char **argv)
     if (options.saltmethod != 2) //if not SALT_METHOD_OWCZARZY
     	options.saltmethod = SALT_METHOD_SANTALUCIA; //then force SALT_METHOD_SANTALUCIA
 
-
     if (options.salt < 0.01 || options.salt > 0.3) //if salt value out of literature values
     	options.salt = DEF_SALT; //set to default
 
diff --git a/src/libecoprimer/ecoprimer.h b/src/libecoprimer/ecoprimer.h
index 2ad99bf..b2f61fd 100644
--- a/src/libecoprimer/ecoprimer.h
+++ b/src/libecoprimer/ecoprimer.h
@@ -290,6 +290,7 @@ typedef struct {
 	float		  salt;
 	PNNParams	  pnparm;
 	bool_t		  print_sets_of_primers;
+	float		  specificity_threshold;
 } options_t, *poptions_t;
 
 typedef ecoseq_t  **pecodnadb_t;