Fixed a bug when primers would touch or overlap, producing amplicons of

size <= 0

VERSION

@@ -1 +1 @@
-0.7.0
+0.8.0

src/ecopcr.c

@@ -6,7 +6,7 @@
 #include <getopt.h>
 
 
-#define VERSION "0.7.0"
+#define VERSION "1.0.0"
 
 
 /* ----------------------------------------------- */
@@ -80,7 +80,7 @@ static void PrintHelp()
         PP      "column 21 : sequence\n");
         PP      "column 22 : definition\n");
         PP      "------------------------------------------\n");
-        PP		" http://www.grenoble.prabi.fr/trac/ecoPCR/\n");
+        PP		" https://git.metabarcoding.org/obitools/ecopcr/wikis/home\n");
         PP      "------------------------------------------\n\n");        
         PP      "\n");
 
@@ -566,11 +566,11 @@ int main(int argc, char **argv)
 		        )
 		     {
 		
-				scname = taxonomy->taxons->taxon[seq->taxid].name;
+				//scname = taxonomy->taxons->taxon[seq->taxid].name;
-				strncpy(head,seq->SQ,10);
+				//strncpy(head,seq->SQ,10);
-				head[10]=0;
+				//head[10]=0;
-				strncpy(tail,seq->SQ+seq->SQ_length-10,10);
+				//strncpy(tail,seq->SQ+seq->SQ_length-10,10);
-				tail[10]=0;
+				//tail[10]=0;
 		
 				apatseq=ecoseq2apatseq(seq,apatseq,circular);
 				
@@ -613,15 +613,17 @@ int main(int argc, char **argv)
 										errj  =apatseq->hiterr[1]->val[j];
 										length = 0;
 										if (posj > posi)
-											length=posj - posi  - o1->patlen - o2->patlen; 
+											length = posj - posi - o1->patlen - o2->patlen;
 										if (posj < posi)
-											length= posj + apatseq->seqlen - posi - o1->patlen - o2->patlen;
+											//length = posj + apatseq->seqlen - posi - o1->patlen - o2->patlen;   // TODO
-										if (length &&
+											length = posi - posj - o1->patlen - o2->patlen;
+										if ((length>0) &&	// For when primers touch or overlap
 											(!lmin || (length >= lmin)) &&
-										    (!lmax || (length <= lmax)))
+											(!lmax || (length <= lmax)))
+										{
 										    printRepeat(seq,oligo1,oligo2,&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);
-	
+										}
 									}
 								}
 							}
@@ -666,15 +668,18 @@ int main(int argc, char **argv)
 
 										length = 0;
 										if (posj > posi)
-											length=posj - posi + 1  - o2->patlen - o1->patlen; /* - o1->patlen : suppress by <EC> */
+											//length = posj - posi + 1  - o2->patlen - o1->patlen; /* - o1->patlen : deleted by <EC> (prior to the OBITools3) */ TODO ????
+											length = posj - posi - o2->patlen - o1->patlen;
 										if (posj < posi)
-											length= posj + apatseq->seqlen - posi - o1->patlen - o2->patlen;
+											//length = posj + apatseq->seqlen - posi - o1->patlen - o2->patlen; TODO ????
-
+											length = posi - posj - o2->patlen - o1->patlen;
-										if (length &&
+										if ((length>0) &&	// For when primers touch or overlap
 											(!lmin || (length >= lmin)) &&
-										    (!lmax || (length <= lmax)))
+											(!lmax || (length <= lmax)))
+										{
 										    printRepeat(seq,oligo1,oligo2,&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);
+										}
 									}
 								}
 							}

tools/ecoPCRFormat.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/env python2.7
 
 import re
 import gzip
@@ -7,11 +7,8 @@ import sys
 import time
 import getopt
 
-try:
+_dbenable=False
-    import psycopg2
+
-    _dbenable=True
-except ImportError:
-    _dbenable=False
 
 #####
 #
@@ -221,56 +218,7 @@ def readTaxonomyDump(taxdir):
     
     return taxonomy,ranks,alternativeName,index
 
-def readTaxonomyDB(dbname):
-    connection = psycopg2.connect(database=dbname)
-    
-    cursor = connection.cursor()
-    cursor.execute("select numid,rank,parent from ncbi_taxonomy.taxon")
-    taxonomy=[list(x) for x in cursor]
 
-    cursor.execute("select rank_class from ncbi_taxonomy.taxon_rank_class order by rank_class")
-    ranks=cursor.fetchall()
-    ranks = dict(map(None,(x[0] for x in ranks),xrange(len(ranks))))
-    
-    print >>sys.stderr,"Sorting taxons..."
-    taxonomy.sort(taxonCmp)
-
-    print >>sys.stderr,"Indexing taxonomy..."
-    index = {}
-    for t in taxonomy:
-        index[t[0]]=bsearchTaxon(taxonomy, t[0])
-    
-    print >>sys.stderr,"Indexing parent and rank..."
-    for t in taxonomy:
-        t[1]=ranks[t[1]]
-        try:
-            t[2]=index[t[2]]
-        except KeyError,e:
-            if t[2] is None and t[0]==1:
-                t[2]=index[t[0]]
-            else:
-                raise e
-        
-    cursor.execute("select taxid,name,category from ncbi_taxonomy.name")
-    
-    alternativeName=[]
-    for taxid,name,classname in cursor:
-        alternativeName.append((name,classname,index[taxid]))
-        if classname == 'scientific name':
-            taxonomy[index[taxid]].append(name)
-    
-    cursor.execute("select old_numid,current_numid from ncbi_taxonomy.taxon_id_alias")
- 
-    print >>sys.stderr,"Adding taxid alias..."
-    for taxid,current in cursor:
-        if current is not None:
-            index[taxid]=index[current]
-        else:
-            index[taxid]=None
-            
-       
-    return taxonomy,ranks,alternativeName,index
-    
 #####
 #
 #
@@ -585,10 +533,9 @@ def ecoParseOptions(arguments):
            }
     
     o,filenames = getopt.getopt(arguments,
-                                'ht:T:n:gfe',
+                                'ht:n:gfe',
                                 ['help',
                                  'taxonomy=',
-                                 'taxonomy_db=',
                                  'name=',
                                  'genbank',
                                  'fasta',
@@ -601,9 +548,6 @@ def ecoParseOptions(arguments):
         elif name in ('-t','--taxonomy'):
             opt['taxmod']='dump'
             opt['taxdir']=value
-        elif name in ('-T','--taxonomy_db'):
-            opt['taxmod']='db'
-            opt['taxdb']=value
         elif name in ('-n','--name'):
             opt['prefix']=value
         elif name in ('-g','--genbank'):
@@ -622,6 +566,7 @@ def ecoParseOptions(arguments):
 
     return opt,filenames
 
+
 def printHelp():
     print "-----------------------------------"
     print " ecoPCRFormat.py"
@@ -641,11 +586,7 @@ if __name__ == '__main__':
     
     opt,filenames = ecoParseOptions(sys.argv[1:])
     
-    if opt['taxmod']=='dump':
+    taxonomy = readTaxonomyDump(opt['taxdir'])
-        taxonomy = readTaxonomyDump(opt['taxdir'])
+    
-    elif opt['taxmod']=='db':
-        taxonomy = readTaxonomyDB(opt['taxdb'])
-        
-   
     ecoDBWriter(opt['prefix'], taxonomy, filenames, opt['parser'])