Initial commit for this soft

git-svn-id: https://www.grenoble.prabi.fr/svn/LECASofts/ecoPrimers/trunk@175 60f365c0-8329-0410-b2a4-ec073aeeaa1d

src/libecoPCR/Makefile

@@ -0,0 +1,31 @@
+
+SOURCES = ecoapat.c \
+         ecodna.c \
+         ecoError.c \
+         ecoIOUtils.c \
+         ecoMalloc.c \
+         ecorank.c \
+         ecoseq.c \
+         ecotax.c \
+         ecofilter.c \
+         econame.c     
+
+SRCS=$(SOURCES)
+         
+OBJECTS= $(patsubst %.c,%.o,$(SOURCES))
+
+LIBFILE= libecoPCR.a
+RANLIB= ranlib
+
+
+include ../global.mk
+
+
+all: $(LIBFILE)
+
+clean:
+	rm -rf $(OBJECTS) $(LIBFILE)
+
+$(LIBFILE): $(OBJECTS)
+	ar -cr $@ $?
+	$(RANLIB) $@

src/libecoPCR/ecoError.c

@@ -0,0 +1,26 @@
+#include "ecoPCR.h"
+#include <stdio.h>
+#include <stdlib.h>
+
+/*
+ * print the message given as argument and exit the program
+ * @param error		error number	
+ * @param message 	the text explaining what's going on
+ * @param filename	the file source where the program failed
+ * @param linenumber	the line where it has failed
+ * filename and linenumber are written at pre-processing 
+ * time by a macro
+ */
+void ecoError(int32_t error,
+              const char* message,
+              const char * filename,
+              int linenumber)
+{
+	fprintf(stderr,"Error %d in file %s line %d : %s\n",
+	               error,
+	               filename,
+	               linenumber,
+	               message);
+	
+	abort();
+}

src/libecoPCR/ecoIOUtils.c

@@ -0,0 +1,122 @@
+#include "ecoPCR.h"
+#include <stdio.h>
+#include <stdlib.h>
+
+#define SWAPINT32(x)     ((((x) << 24) & 0xFF000000) | (((x) <<  8) & 0xFF0000) | \
+                         (((x) >>  8) & 0xFF00)     | (((x) >> 24) & 0xFF))
+
+
+int32_t is_big_endian()
+{
+	int32_t i=1;
+	
+	return (int32_t)((char*)&i)[0];
+}
+
+
+
+
+int32_t swap_int32_t(int32_t i)
+{
+	return SWAPINT32(i);
+}
+
+
+/**
+ * Read part of the file
+ * @param	*f	the database
+ * @param	recordSize the size to be read
+ * 
+ * @return	buffer
+ */
+void *read_ecorecord(FILE *f,int32_t *recordSize)
+{
+	static void *buffer    =NULL;
+	int32_t      buffersize=0;
+	int32_t      read;
+	
+	if (!recordSize)
+		ECOERROR(ECO_ASSERT_ERROR,
+		         "recordSize cannot be NULL");
+		
+	read = fread(recordSize,
+	      		 1,
+	      		 sizeof(int32_t),
+	             f);
+	             
+	if (feof(f))
+		return NULL;
+	             
+	if (read != sizeof(int32_t))
+		ECOERROR(ECO_IO_ERROR,"Reading record size error");
+		
+	if (is_big_endian())
+		*recordSize=swap_int32_t(*recordSize);
+		
+	if (buffersize < *recordSize)
+	{
+		if (buffer)
+			buffer = ECOREALLOC(buffer,*recordSize,
+			                    "Increase size of record buffer");
+		else
+			buffer = ECOMALLOC(*recordSize,
+			                    "Allocate record buffer");
+	}
+	
+	read = fread(buffer,
+	             1,
+				 *recordSize,
+				 f);
+				 
+	if (read != *recordSize)
+		ECOERROR(ECO_IO_ERROR,"Reading record data error");
+		
+	return buffer;	 
+};
+
+
+
+
+
+/**
+ * Open the database and check it's readable
+ * @param 	filename 		name of the database (.sdx, .rdx, .tbx)
+ * @param 	sequencecount	buffer - pointer to variable storing the number of occurence    
+ * @param 	abort_on_open_error		 	boolean to define the behaviour in case of error 
+ * 										while opening the database
+ * @return 	FILE type
+ **/
+FILE *open_ecorecorddb(const char *filename,
+                       int32_t    *sequencecount,
+                       int32_t    abort_on_open_error)
+{
+    FILE        *f;
+	int32_t      read;
+	
+	f = fopen(filename,"rb");
+	
+	if (!f)
+		{
+			if (abort_on_open_error)
+		 		ECOERROR(ECO_IO_ERROR,"Cannot open file");
+		 	else
+		 	{
+		 		*sequencecount=0;
+		 		return NULL;
+		 	}
+		}
+		
+	read = fread(sequencecount,
+	      		 1,
+	      		 sizeof(int32_t),
+	      		 f);
+	             
+	if (read != sizeof(int32_t))
+		ECOERROR(ECO_IO_ERROR,"Reading record size error");
+
+	if (is_big_endian())
+		*sequencecount=swap_int32_t(*sequencecount);
+		
+	return f;                  
+}
+

src/libecoPCR/ecoMalloc.c

@@ -0,0 +1,79 @@
+#include "ecoPCR.h"
+#include <stdlib.h>
+
+static int eco_log_malloc = 0;
+
+void    eco_trace_memory_allocation()
+{
+	eco_log_malloc=1;
+}
+
+void    eco_untrace_memory_allocation()
+{
+	eco_log_malloc=0;
+}
+
+
+void   *eco_malloc(int32_t chunksize,
+                   const char *error_message,
+                   const char *filename,
+                   int32_t    line)
+{
+	void * chunk;
+	
+	chunk = calloc(1,chunksize);
+	
+	if (!chunk)
+		ecoError(ECO_MEM_ERROR,error_message,filename,line);
+		
+	if (eco_log_malloc)
+		fprintf(stderr,
+			    "Memory segment located at %p of size %d is allocated (file : %s [%d])",
+			    chunk,
+			    chunksize,
+			    filename,
+			    line);
+		
+	return chunk;
+}
+
+void   *eco_realloc(void *chunk,
+                    int32_t newsize,
+                    const char *error_message,
+                    const char *filename,
+                    int32_t    line)
+{
+	void *newchunk;
+	
+	newchunk = realloc(chunk,newsize);
+	
+	if (!newchunk)
+		ecoError(ECO_MEM_ERROR,error_message,filename,line);
+
+	if (eco_log_malloc)
+		fprintf(stderr,
+			    "Old memory segment %p is reallocated at %p with a size of %d (file : %s [%d])",
+			    chunk,
+			    newchunk,
+			    newsize,
+			    filename,
+			    line);
+		
+	return newchunk;	
+}
+
+void    eco_free(void *chunk,
+                 const char *error_message,
+                 const char *filename,
+                 int32_t    line)
+{
+	free(chunk);
+	
+	if (eco_log_malloc)
+		fprintf(stderr,
+			    "Memory segment %p is released => %s (file : %s [%d])",
+			    chunk,
+			    error_message,
+			    filename,
+			    line);
+}

src/libecoPCR/ecoPCR.h

@@ -0,0 +1,269 @@
+#ifndef ECOPCR_H_
+#define ECOPCR_H_
+
+#include <stdio.h>
+#include <inttypes.h>
+
+#ifndef H_apat
+#include "../libapat/apat.h"
+#endif
+
+/*****************************************************
+ * 
+ *  Data type declarations
+ * 
+ *****************************************************/
+
+/*
+ * 
+ *  Sequence types
+ * 
+ */
+
+typedef struct {
+	
+	int32_t  taxid;
+	char     AC[20];
+	int32_t  DE_length;
+	int32_t  SQ_length;
+	int32_t  CSQ_length;
+	
+	char     data[1];
+	
+} ecoseqformat_t;
+
+typedef struct {
+	int32_t taxid;
+	int32_t SQ_length;
+	char    *AC;
+	char    *DE;
+	char    *SQ;
+} ecoseq_t;
+
+/*
+ * 
+ * Taxonomy taxon types
+ * 
+ */
+
+
+typedef struct {
+	int32_t  taxid;
+	int32_t  rank;
+	int32_t	 parent;
+	int32_t  namelength;
+	char     name[1];
+	
+} ecotxformat_t;
+
+typedef struct ecotxnode {
+	int32_t           taxid;
+	int32_t           rank;
+	struct ecotxnode  *parent;
+	char              *name;
+} ecotx_t;
+
+typedef struct {
+	int32_t count;
+	ecotx_t taxon[1];
+} ecotxidx_t;
+ 
+	
+/*
+ * 
+ * Taxonomy rank types
+ * 
+ */
+	
+typedef struct {
+	int32_t count;
+	char*   label[1];
+} ecorankidx_t;
+
+/*
+ * 
+ * Taxonomy name types
+ * 
+ */
+
+typedef struct {
+ 	int32_t is_scientificname;
+	int32_t  namelength;
+	int32_t  classlength;
+	int32_t  taxid;
+	char     names[1];	
+} econameformat_t;
+ 
+ 
+ typedef struct {
+ 	char 	*name;
+ 	char 	*classname;
+ 	int32_t is_scientificname;
+ 	struct ecotxnode  *taxon;
+} econame_t;
+
+ 
+typedef struct {
+	int32_t count;
+	econame_t   names[1];
+} econameidx_t;
+
+
+ typedef struct {
+	ecorankidx_t *ranks;
+	econameidx_t *names;
+	ecotxidx_t   *taxons;
+} ecotaxonomy_t;
+
+ 
+/*****************************************************
+ * 
+ *  Function declarations
+ * 
+ *****************************************************/
+
+/*
+ * 
+ * Low level system functions
+ * 
+ */
+
+int32_t is_big_endian();
+int32_t swap_int32_t(int32_t);
+
+void   *eco_malloc(int32_t chunksize,
+                   const char *error_message,
+                   const char *filename,
+                   int32_t    line);
+                   
+                   
+void   *eco_realloc(void *chunk,
+                    int32_t chunksize,
+                    const char *error_message,
+                    const char *filename,
+                    int32_t    line);
+                    
+void    eco_free(void *chunk,
+                 const char *error_message,
+                 const char *filename,
+                 int32_t    line);
+                 
+void    eco_trace_memory_allocation();
+void    eco_untrace_memory_allocation();
+
+#define ECOMALLOC(size,error_message) \
+	    eco_malloc((size),(error_message),__FILE__,__LINE__)
+	   
+#define ECOREALLOC(chunk,size,error_message) \
+        eco_realloc((chunk),(size),(error_message),__FILE__,__LINE__)
+        
+#define ECOFREE(chunk,error_message) \
+        eco_free((chunk),(error_message),__FILE__,__LINE__)
+        
+        
+
+
+/*
+ * 
+ * Error managment
+ * 
+ */
+ 
+  
+void ecoError(int32_t,const char*,const char *,int);
+
+#define ECOERROR(code,message) ecoError((code),(message),__FILE__,__LINE__)
+
+#define ECO_IO_ERROR       (1)
+#define ECO_MEM_ERROR      (2)
+#define ECO_ASSERT_ERROR   (3)
+#define ECO_NOTFOUND_ERROR (4)
+
+
+/*
+ * 
+ * Low level Disk access functions
+ * 
+ */
+
+FILE *open_ecorecorddb(const char *filename,
+                       int32_t    *sequencecount,
+                       int32_t    abort_on_open_error);
+                       
+void *read_ecorecord(FILE *,int32_t *recordSize);
+
+
+
+/* 
+ *   Read function in internal binary format
+ */
+
+FILE             *open_ecoseqdb(const char *filename,
+                                int32_t    *sequencecount);
+                                                                
+ecoseq_t         *readnext_ecoseq(FILE *);
+
+ecorankidx_t     *read_rankidx(const char *filename);
+
+econameidx_t     *read_nameidx(const char *filename,ecotaxonomy_t *taxonomy);
+
+
+
+	/**
+	 * Read taxonomy data as formated by the ecoPCRFormat.py script.
+	 * 
+	 * This function is normaly uses internaly by the read_taxonomy
+	 * function and should not be called directly.
+	 * 
+	 * @arg filename  path to the *.tdx file of the reformated db
+	 * 
+	 * @return pointer to a taxonomy index structure
+	 */
+ 
+ecotxidx_t       *read_taxonomyidx(const char *filename);
+
+ecotaxonomy_t    *read_taxonomy(const char *prefix,int32_t readAlternativeName);
+
+ecotx_t *eco_findtaxonbytaxid(ecotaxonomy_t *taxonomy, int32_t taxid);
+
+int eco_isundertaxon(ecotx_t *taxon, int other_taxid);
+
+ecoseq_t *ecoseq_iterator(const char *prefix);
+
+
+
+ecoseq_t *new_ecoseq();
+int32_t   delete_ecoseq(ecoseq_t *);
+ecoseq_t *new_ecoseq_with_data( char *AC,
+								char *DE,
+								char *SQ,
+								int32_t   taxid
+								);
+
+
+int32_t delete_taxon(ecotx_t *taxon);
+int32_t delete_taxonomy(ecotxidx_t *index);
+
+
+int32_t rank_index(const char* label,ecorankidx_t* ranks);
+
+int32_t  delete_apatseq(SeqPtr pseq);
+PatternPtr buildPattern(const char *pat, int32_t error_max);
+PatternPtr complementPattern(PatternPtr pat);
+
+SeqPtr ecoseq2apatseq(ecoseq_t *in,SeqPtr out,int32_t circular);
+
+char *ecoComplementPattern(char *nucAcSeq);
+char *ecoComplementSequence(char *nucAcSeq);
+char *getSubSequence(char* nucAcSeq,int32_t begin,int32_t end);
+
+ecotx_t *eco_getspecies(ecotx_t *taxon,ecotaxonomy_t *taxonomy);
+ecotx_t *eco_getgenus(ecotx_t *taxon,ecotaxonomy_t *taxonomy);
+ecotx_t *eco_getfamily(ecotx_t *taxon,ecotaxonomy_t *taxonomy);
+ecotx_t *eco_getkingdom(ecotx_t *taxon,ecotaxonomy_t *taxonomy);
+ecotx_t *eco_getsuperkingdom(ecotx_t *taxon,ecotaxonomy_t *taxonomy);
+
+int eco_is_taxid_ignored(int32_t *ignored_taxid, int32_t tab_len, int32_t taxid);
+int eco_is_taxid_included(ecotaxonomy_t *taxonomy, int32_t *included_taxid, int32_t tab_len, int32_t taxid);
+
+#endif /*ECOPCR_H_*/

src/libecoPCR/ecoapat.c

@@ -0,0 +1,199 @@
+#include "../libapat/libstki.h"
+#include "../libapat/apat.h"
+
+#include "ecoPCR.h"
+
+#include <string.h>
+
+static void EncodeSequence(SeqPtr seq);
+static void UpperSequence(char *seq);
+
+/* -------------------------------------------- */
+/* uppercase sequence                           */
+/* -------------------------------------------- */
+
+#define IS_LOWER(c) (((c) >= 'a') && ((c) <= 'z'))
+#define TO_UPPER(c) ((c) - 'a' + 'A')
+
+void UpperSequence(char *seq)
+{
+        char *cseq;
+
+        for (cseq = seq ; *cseq ; cseq++) 
+            if (IS_LOWER(*cseq))
+                *cseq = TO_UPPER(*cseq);
+}
+ 
+#undef IS_LOWER
+#undef TO_UPPER
+
+
+
+
+/* -------------------------------------------- */
+/* encode sequence                              */
+/* IS_UPPER is slightly faster than isupper     */
+/* -------------------------------------------- */
+
+#define IS_UPPER(c) (((c) >= 'A') && ((c) <= 'Z'))
+
+
+
+void EncodeSequence(SeqPtr seq)
+{
+        int   i;
+        UInt8 *data;
+        char  *cseq;
+
+        data = seq->data;
+        cseq = seq->cseq;
+
+        while (*cseq) {
+
+            *data = (IS_UPPER(*cseq) ? *cseq - 'A' : 0x0);
+            data++;
+            cseq++;
+        }
+        
+        for (i=0,cseq=seq->cseq;i < seq->circular; i++,cseq++,data++)
+            *data = (IS_UPPER(*cseq) ? *cseq - 'A' : 0x0);
+
+        for (i = 0 ; i < MAX_PATTERN ; i++)
+            seq->hitpos[i]->top = seq->hiterr[i]->top = 0;
+
+}
+
+#undef IS_UPPER
+
+
+SeqPtr ecoseq2apatseq(ecoseq_t *in,SeqPtr out,int32_t circular)
+{
+        int    i;
+
+		if (!out)
+		{
+			out = ECOMALLOC(sizeof(Seq),
+			                "Error in Allocation of a new Seq structure");
+			          
+	        for (i  = 0 ; i < MAX_PATTERN ; i++) 
+	        {
+			   
+	           if (! (out->hitpos[i] = NewStacki(kMinStackiSize)))
+	             	ECOERROR(ECO_MEM_ERROR,"Error in hit stack Allocation");
+	
+	           if (! (out->hiterr[i] = NewStacki(kMinStackiSize)))
+	            	ECOERROR(ECO_MEM_ERROR,"Error in error stack Allocation");
+	        }
+		}
+
+		
+		out->name    = in->AC;
+		out->seqsiz = out->seqlen = in->SQ_length;
+		out->circular = circular;
+		
+		if (!out->data)
+		{
+			out->data = ECOMALLOC((out->seqlen+circular) *sizeof(UInt8),
+		    	     			  "Error in Allocation of a new Seq data member");    
+		   	out->datsiz=  out->seqlen+circular;
+		}
+		else if ((out->seqlen +circular) >= out->datsiz)
+		{
+			out->data = ECOREALLOC(out->data,(out->seqlen+circular),
+			                      "Error during Seq data buffer realloc");
+		   	out->datsiz=  out->seqlen+circular;			                      
+		}
+
+		out->cseq = in->SQ;
+		
+		EncodeSequence(out);
+
+        return out;
+}
+
+int32_t delete_apatseq(SeqPtr pseq)
+{
+         int i;
+
+        if (pseq) {
+
+            if (pseq->data) 
+            	ECOFREE(pseq->data,"Freeing sequence data buffer");
+
+            for (i = 0 ; i < MAX_PATTERN ; i++) {
+                if (pseq->hitpos[i]) FreeStacki(pseq->hitpos[i]);
+                if (pseq->hiterr[i]) FreeStacki(pseq->hiterr[i]);
+            }
+
+            ECOFREE(pseq,"Freeing apat sequence structure");
+            
+            return 0;
+        }
+        
+        return 1;
+}
+
+PatternPtr buildPattern(const char *pat, int32_t error_max)
+{
+	PatternPtr pattern;
+	int32_t    patlen;
+	
+	pattern = ECOMALLOC(sizeof(Pattern),
+						"Error in pattern allocation");
+						
+	pattern->ok      = Vrai;
+	pattern->hasIndel= Faux;
+	pattern->maxerr  = error_max;
+	patlen  = strlen(pat);
+	
+	pattern->cpat    = ECOMALLOC(sizeof(char)*patlen+1,
+	                             "Error in sequence pattern allocation");
+	                             
+	strncpy(pattern->cpat,pat,patlen);
+	pattern->cpat[patlen]=0;
+	UpperSequence(pattern->cpat);
+	
+	if (!CheckPattern(pattern))
+		ECOERROR(ECO_ASSERT_ERROR,"Error in pattern checking");
+		
+	if (! EncodePattern(pattern, dna))
+		ECOERROR(ECO_ASSERT_ERROR,"Error in pattern encoding");
+
+   	if (! CreateS(pattern, ALPHA_LEN))
+		ECOERROR(ECO_ASSERT_ERROR,"Error in pattern compiling");
+	
+	return pattern;
+		
+}
+
+PatternPtr complementPattern(PatternPtr pat)
+{
+	PatternPtr pattern;
+	
+	pattern = ECOMALLOC(sizeof(Pattern),
+						"Error in pattern allocation");
+						
+	pattern->ok      = Vrai;
+	pattern->hasIndel= pat->hasIndel;
+	pattern->maxerr  = pat->maxerr;
+	pattern->patlen  = pat->patlen;
+	
+	pattern->cpat    = ECOMALLOC(sizeof(char)*(strlen(pat->cpat)+1),
+	                             "Error in sequence pattern allocation");
+	                             
+	strcpy(pattern->cpat,pat->cpat);
+	
+	ecoComplementPattern(pattern->cpat);
+	
+	if (!CheckPattern(pattern))
+		ECOERROR(ECO_ASSERT_ERROR,"Error in pattern checking");
+		
+	if (! EncodePattern(pattern, dna))
+		ECOERROR(ECO_ASSERT_ERROR,"Error in pattern encoding");
+
+   	if (! CreateS(pattern, ALPHA_LEN))
+		ECOERROR(ECO_ASSERT_ERROR,"Error in pattern compiling");
+	
+	return pattern;
+		
+}

src/libecoPCR/ecodna.c

@@ -0,0 +1,153 @@
+#include <string.h>
+#include "ecoPCR.h"
+
+/*
+ * @doc: DNA alphabet (IUPAC)
+ */
+#define LX_BIO_DNA_ALPHA   "ABCDEFGHIJKLMNOPQRSTUVWXYZ#![]"
+
+/*
+ * @doc: complementary DNA alphabet (IUPAC)
+ */
+#define LX_BIO_CDNA_ALPHA  "TVGHEFCDIJMLKNOPQYSAABWXRZ#!]["
+
+
+static char sNuc[]     = LX_BIO_DNA_ALPHA;
+static char sAnuc[]    = LX_BIO_CDNA_ALPHA;
+
+static char LXBioBaseComplement(char nucAc);
+static char *LXBioSeqComplement(char *nucAcSeq);
+static char *reverseSequence(char *str,char isPattern);
+
+ 
+/* ---------------------------- */
+
+char LXBioBaseComplement(char nucAc)
+{
+    char *c;
+
+    if ((c = strchr(sNuc, nucAc)))
+        return sAnuc[(c - sNuc)];
+    else
+        return nucAc;
+}
+
+/* ---------------------------- */
+
+char *LXBioSeqComplement(char *nucAcSeq)
+{
+    char *s;
+
+    for (s = nucAcSeq ; *s ; s++)
+        *s = LXBioBaseComplement(*s);
+
+    return nucAcSeq;
+}
+
+
+char *reverseSequence(char *str,char isPattern)
+{
+        char *sb, *se, c;
+
+        if (! str)
+            return str;
+            
+        sb = str;
+        se = str + strlen(str) - 1;
+
+        while(sb <= se) {
+           c    = *sb;
+          *sb++ = *se;
+          *se-- = c;
+        }
+
+		sb = str;
+		se = str + strlen(str) - 1;
+		
+		if (isPattern)
+			for (;sb < se; sb++)
+			{
+				if (*sb=='#')
+				{
+					if (((se - sb) > 2) && (*(sb+2)=='!'))
+					{
+						*sb='!';
+						sb+=2;
+						*sb='#';
+					}
+					else
+					{
+						*sb=*(sb+1);
+						sb++;
+						*sb='#';
+					}
+				}
+				else if (*sb=='!')
+					{
+						*sb=*(sb-1);
+						*(sb-1)='!';
+					}
+			}
+
+        return str;
+}
+
+char *ecoComplementPattern(char *nucAcSeq)
+{
+    return reverseSequence(LXBioSeqComplement(nucAcSeq),1);
+}
+
+char *ecoComplementSequence(char *nucAcSeq)
+{
+    return reverseSequence(LXBioSeqComplement(nucAcSeq),0);
+}
+
+
+char *getSubSequence(char* nucAcSeq,int32_t begin,int32_t end)
+{
+	static char *buffer  = NULL;
+	static int32_t buffSize= 0;
+	int32_t length;
+	
+	if (begin < end)
+	{
+		length = end - begin;
+		
+		if (length >= buffSize)
+		{
+			buffSize = length+1;
+			if (buffer)
+				buffer=ECOREALLOC(buffer,buffSize,
+						   	      "Error in reallocating sub sequence buffer");
+			else
+				buffer=ECOMALLOC(buffSize,
+				          		 "Error in allocating sub sequence buffer");
+				
+		}
+		
+		strncpy(buffer,nucAcSeq + begin,length);
+		buffer[length]=0;
+	}
+	else
+	{
+		length = end + strlen(nucAcSeq) - begin;
+		
+		if (length >= buffSize)
+		{
+			buffSize = length+1;
+			if (buffer)
+				buffer=ECOREALLOC(buffer,buffSize,
+						   	      "Error in reallocating sub sequence buffer");
+			else
+				buffer=ECOMALLOC(buffSize,
+				          		 "Error in allocating sub sequence buffer");
+				
+		}
+		strncpy(buffer,nucAcSeq+begin,length - end);
+		strncpy(buffer+(length-end),nucAcSeq ,end);
+		buffer[length]=0;
+	}
+	
+	return buffer;
+}
+

src/libecoPCR/ecofilter.c

@@ -0,0 +1,19 @@
+#include "ecoPCR.h"
+
+int eco_is_taxid_included(	ecotaxonomy_t *taxonomy, 
+							int32_t *restricted_taxid, 
+							int32_t tab_len, 
+							int32_t taxid)
+{
+	int i;
+	ecotx_t *taxon;
+	
+	taxon = eco_findtaxonbytaxid(taxonomy, taxid);
+	
+	for (i=0; i < tab_len; i++)
+		if ( (taxon->taxid == restricted_taxid[i]) ||
+			 (eco_isundertaxon(taxon, restricted_taxid[i])) )
+			return 1;
+	
+	return 0;
+}

src/libecoPCR/econame.c

@@ -0,0 +1,61 @@
+#include "ecoPCR.h"
+#include <string.h>
+#include <stdlib.h>
+
+static econame_t *readnext_econame(FILE *f,econame_t *name,ecotaxonomy_t *taxonomy);
+
+econameidx_t *read_nameidx(const char *filename,ecotaxonomy_t *taxonomy)
+{
+
+	int32_t      		count;
+	FILE         		*f;
+	econameidx_t		*indexname;
+	int32_t      		i;
+	
+	f = open_ecorecorddb(filename,&count,1);
+	
+	indexname = (econameidx_t*) ECOMALLOC(sizeof(econameidx_t) + sizeof(econame_t) * (count-1),"Allocate names");
+	
+	indexname->count=count;
+	                                    
+	for (i=0; i < count; i++){
+		readnext_econame(f,(indexname->names)+i,taxonomy);
+	}
+
+	return indexname;
+}
+
+econame_t *readnext_econame(FILE *f,econame_t *name,ecotaxonomy_t *taxonomy)
+{
+	
+	econameformat_t *raw;
+	int32_t  rs;
+	
+	raw = read_ecorecord(f,&rs);
+	
+	if (!raw)
+		return NULL;
+
+	if (is_big_endian())
+	{
+		raw->is_scientificname 	= swap_int32_t(raw->is_scientificname);
+		raw->namelength 	    = swap_int32_t(raw->namelength);
+		raw->classlength        = swap_int32_t(raw->classlength);
+		raw->taxid  	        = swap_int32_t(raw->taxid); 
+	}
+	
+	name->is_scientificname=raw->is_scientificname;
+	
+	name->name   	= ECOMALLOC((raw->namelength+1) * sizeof(char),"Allocate name");
+	strncpy(name->name,raw->names,raw->namelength);
+	name->name[raw->namelength]=0;
+	
+	name->classname = ECOMALLOC((raw->classlength+1) * sizeof(char),"Allocate classname");
+	strncpy(name->classname,(raw->names+raw->namelength),raw->classlength);
+	name->classname[raw->classlength]=0;
+	
+	name->taxon = taxonomy->taxons->taxon + raw->taxid;
+
+	return name;
+}
+

src/libecoPCR/ecorank.c

@@ -0,0 +1,52 @@
+#include "ecoPCR.h"
+#include <string.h>
+#include <stdlib.h>
+
+static int compareRankLabel(const void *label1, const void *label2);
+
+ecorankidx_t     *read_rankidx(const char *filename)
+{
+	int32_t      count;
+	FILE         *f;
+	ecorankidx_t *index;
+	int32_t      i;
+	int32_t      rs;
+	char         *buffer;
+	
+	f = open_ecorecorddb(filename,&count,1);
+	
+	index = (ecorankidx_t*) ECOMALLOC(sizeof(ecorankidx_t) + sizeof(char*) * (count-1),
+	                                  "Allocate rank index");
+	 
+	index->count=count;                                 
+	          
+	for (i=0; i < count; i++)
+		{
+			buffer = read_ecorecord(f,&rs);
+			index->label[i]=(char*) ECOMALLOC(rs+1,
+			                                  "Allocate rank label");
+			strncpy(index->label[i],buffer,rs);
+		}
+		
+	return index;
+}
+
+int32_t rank_index(const char* label,ecorankidx_t* ranks)
+{
+	char **rep;
+	
+	rep = bsearch(label,ranks->label,ranks->count,sizeof(char*),compareRankLabel);
+	
+	if (rep)
+		return rep-ranks->label;
+	else
+		ECOERROR(ECO_NOTFOUND_ERROR,"Rank label not found");
+		
+	return -1;
+}
+
+
+int compareRankLabel(const void *label1, const void *label2)
+{
+	return strcmp((const char*)label1,*(const char**)label2);
+}

src/libecoPCR/ecoseq.c

@@ -0,0 +1,223 @@
+#include "ecoPCR.h"
+#include <stdlib.h>
+#include <string.h>
+#include <zlib.h>
+#include <string.h>
+#include <stdio.h>
+
+static FILE *open_seqfile(const char *prefix,int32_t index);
+
+
+ecoseq_t *new_ecoseq()
+{
+	void *tmp;
+	
+	tmp = ECOMALLOC(sizeof(ecoseq_t),"Allocate new ecoseq structure");
+	
+	return tmp;
+}
+
+int32_t delete_ecoseq(ecoseq_t * seq)
+{
+	
+	if (seq)
+	{
+		if (seq->AC)
+			ECOFREE(seq->AC,"Free sequence AC");
+			
+		if (seq->DE)
+			ECOFREE(seq->DE,"Free sequence DE");
+			
+		if (seq->SQ)
+			ECOFREE(seq->SQ,"Free sequence SQ");
+			
+		ECOFREE(seq,"Free sequence structure");
+		
+		return 0;
+			
+	}
+	
+	return 1;
+}
+
+ecoseq_t *new_ecoseq_with_data( char *AC,
+								char *DE,
+								char *SQ,
+								int32_t   taxid_idx
+								)
+{
+	ecoseq_t *tmp;
+	int32_t lstr;
+	tmp = new_ecoseq();
+	
+	tmp->taxid=taxid_idx;
+	
+	if (AC)
+		{
+			lstr =strlen(AC);
+			tmp->AC=ECOMALLOC((lstr+1) * sizeof(char),
+			                  "Allocate sequence accession");
+			strcpy(tmp->AC,AC);
+		}
+
+	if (DE)
+		{
+			lstr =strlen(DE);
+			tmp->DE=ECOMALLOC((lstr+1) * sizeof(char),
+			                  "Allocate sequence definition");
+			strcpy(tmp->DE,DE);
+		}
+
+	if (SQ)
+		{
+			lstr =strlen(SQ);
+			tmp->SQ=ECOMALLOC((lstr+1) * sizeof(char),
+			                  "Allocate sequence data");
+			strcpy(tmp->SQ,SQ);
+		}
+	return tmp;
+
+}
+
+/**
+ * ?? used ??
+ **/
+FILE *open_ecoseqdb(const char *filename,
+                    int32_t    *sequencecount)
+{
+	return open_ecorecorddb(filename,sequencecount,1);
+}
+
+ecoseq_t *readnext_ecoseq(FILE *f)
+{
+	char     *compressed=NULL;
+
+	ecoseqformat_t *raw;
+	ecoseq_t *seq;
+	int32_t  comp_status;
+	unsigned long int seqlength;
+	int32_t  rs;
+	
+	raw = read_ecorecord(f,&rs);
+	
+	if (!raw)
+		return NULL;
+		
+	if (is_big_endian())
+	{
+		raw->CSQ_length = swap_int32_t(raw->CSQ_length);
+		raw->DE_length  = swap_int32_t(raw->DE_length);
+		raw->SQ_length  = swap_int32_t(raw->SQ_length);
+		raw->taxid      = swap_int32_t(raw->taxid);
+	}
+	
+	seq = new_ecoseq();
+	
+	seq->taxid = raw->taxid;
+	
+    seq->AC    = ECOMALLOC(strlen(raw->AC) +1,
+                           "Allocate Sequence Accesion number");
+    strncpy(seq->AC,raw->AC,strlen(raw->AC));
+
+                           
+    seq->DE    = ECOMALLOC(raw->DE_length+1,
+                           "Allocate Sequence definition");
+    strncpy(seq->DE,raw->data,raw->DE_length);
+    
+	seqlength = seq->SQ_length = raw->SQ_length;
+    
+    compressed = raw->data + raw->DE_length;
+
+    seq->SQ = ECOMALLOC(seqlength+1,
+                        "Allocate sequence buffer");
+                        
+    comp_status = uncompress((unsigned char*)seq->SQ,
+                             &seqlength,
+                             (unsigned char*)compressed,
+                             raw->CSQ_length);
+    
+    if (comp_status != Z_OK)
+    	ECOERROR(ECO_IO_ERROR,"I cannot uncompress sequence data");
+    
+	return seq;
+}
+
+/**
+ * Open the sequences database (.sdx file)
+ * @param	prefix	name of the database (radical without extension) 
+ * @param	index 	integer
+ * 
+ * @return	file object 
+ */
+FILE *open_seqfile(const char *prefix,int32_t index)
+{
+	char           filename_buffer[1024];
+	int32_t        filename_length;
+	FILE           *input;
+	int32_t        seqcount;
+
+	filename_length = snprintf(filename_buffer,
+								1023,
+	                           "%s_%03d.sdx",
+	                           prefix,
+	                           index);
+	                           
+			fprintf(stderr,"# Coucou %s\n",filename_buffer);
+	                           
+
+	if (filename_length >= 1024)
+		ECOERROR(ECO_ASSERT_ERROR,"file name is too long");
+
+	filename_buffer[filename_length]=0;
+		
+	input=open_ecorecorddb(filename_buffer,&seqcount,0);
+	
+	if (input)
+		fprintf(stderr,"# Reading file %s containing %d sequences...\n",
+				filename_buffer,
+				seqcount);
+				
+	return input;
+}
+
+ecoseq_t *ecoseq_iterator(const char *prefix)
+{
+	static FILE    *current_seq_file= NULL;
+	static int32_t current_file_idx = 1;
+	static char    current_prefix[1024];
+	ecoseq_t       *seq;
+	
+	if (prefix)
+	{
+		current_file_idx = 1;
+
+		if (current_seq_file)
+			fclose(current_seq_file);
+			
+		strncpy(current_prefix,prefix,1023);
+		current_prefix[1024]=0;
+		
+		current_seq_file = open_seqfile(current_prefix,
+		 							    current_file_idx);
+		
+		if (!current_seq_file)
+			return NULL;
+			
+	}
+	
+	seq = readnext_ecoseq(current_seq_file);
+	
+	if (!seq && feof(current_seq_file))
+	{
+		current_file_idx++;
+		fclose(current_seq_file);
+		current_seq_file = open_seqfile(current_prefix,
+		 							    current_file_idx);
+		 							    
+		
+		if (current_seq_file)
+			seq = readnext_ecoseq(current_seq_file);
+	}
+	
+	return seq;
+}

src/libecoPCR/ecotax.c

@@ -0,0 +1,329 @@
+#include "ecoPCR.h"
+#include <string.h>
+#include <stdlib.h>
+#include <stdio.h>
+
+static ecotx_t *readnext_ecotaxon(FILE *f,ecotx_t *taxon);
+
+ /** 
+ * Open the taxonomy database 
+ * @param	pointer to the database (.tdx file) 
+ * @return	a ecotxidx_t structure 
+ */
+ecotxidx_t     *read_taxonomyidx(const char *filename)
+{
+	int32_t      count;
+	FILE         *f;
+	ecotxidx_t *index;
+	int32_t      i;
+	
+	f = open_ecorecorddb(filename,&count,1);
+	
+	index = (ecotxidx_t*) ECOMALLOC(sizeof(ecotxidx_t) + sizeof(ecotx_t) * (count-1),
+	                                  "Allocate taxonomy");
+	 
+	index->count=count;                                         
+	for (i=0; i < count; i++){
+		readnext_ecotaxon(f,&(index->taxon[i]));
+		index->taxon[i].parent=index->taxon + (int32_t)index->taxon[i].parent;
+	}					
+	return index;
+}
+
+
+int32_t delete_taxonomy(ecotxidx_t *index)
+{
+	int32_t i;
+	
+	if (index)
+	{
+		for (i=0; i< index->count; i++)
+			if (index->taxon[i].name)
+				ECOFREE(index->taxon[i].name,"Free scientific name");
+				
+		ECOFREE(index,"Free Taxonomy");
+		
+		return 0;
+	}
+	
+	return 1;
+}
+
+
+
+int32_t delete_taxon(ecotx_t *taxon)
+{
+	if (taxon)
+	{
+		if (taxon->name)
+			ECOFREE(taxon->name,"Free scientific name");
+			
+		ECOFREE(taxon,"Free Taxon");
+		
+		return 0;
+	}
+		
+	return 1;
+}
+
+
+/**
+ * Read the database for a given taxon a save the data 
+ * into the taxon structure(if any found)
+ * @param	*f	pointer to FILE type returned by fopen
+ * @param	*taxon	pointer to the structure
+ * 
+ * @return	a ecotx_t structure if any taxon found else NULL 
+ */
+ecotx_t *readnext_ecotaxon(FILE *f,ecotx_t *taxon)
+{
+	
+	ecotxformat_t *raw;
+	int32_t  rs;
+	
+	raw = read_ecorecord(f,&rs);
+	
+	if (!raw)
+		return NULL;
+
+	if (is_big_endian())
+	{
+		raw->namelength = swap_int32_t(raw->namelength);
+		raw->parent     = swap_int32_t(raw->parent);
+		raw->rank       = swap_int32_t(raw->rank);
+		raw->taxid      = swap_int32_t(raw->taxid); 
+	}
+	
+	taxon->parent = (ecotx_t*)raw->parent;
+	taxon->taxid  = raw->taxid;
+	taxon->rank   = raw->rank;
+	
+	taxon->name   = ECOMALLOC((raw->namelength+1) * sizeof(char),
+	                          "Allocate taxon scientific name");
+	                          
+	strncpy(taxon->name,raw->name,raw->namelength);
+		
+	return taxon;
+}
+
+
+ecotaxonomy_t    *read_taxonomy(const char *prefix,int32_t readAlternativeName)
+{
+	ecotaxonomy_t *tax;
+	char          *filename;
+	int           buffsize;
+	
+	tax = ECOMALLOC(sizeof(ecotaxonomy_t),
+	                "Allocate taxonomy structure");
+	
+	buffsize = strlen(prefix)+10;
+	
+	filename = ECOMALLOC(buffsize,
+	                     "Allocate filename");
+	
+	snprintf(filename,buffsize,"%s.rdx",prefix);
+	
+	tax->ranks = read_rankidx(filename);
+	
+	snprintf(filename,buffsize,"%s.tdx",prefix);
+	
+	tax->taxons = read_taxonomyidx(filename);
+	
+	if (readAlternativeName)
+	{
+  	   snprintf(filename,buffsize,"%s.ndx",prefix);
+	   tax->names=read_nameidx(filename,tax);
+	}
+	else
+	   tax->names=NULL;
+	return tax;
+	
+}
+
+
+
+int32_t delete_ecotaxonomy(ecotaxonomy_t *taxonomy)
+{
+	if (taxonomy)
+	{
+		if (taxonomy->ranks)
+			ECOFREE(taxonomy->ranks,"Free rank index");
+			
+		if (taxonomy->taxons)
+			ECOFREE(taxonomy->taxons,"Free taxon index");
+			
+		ECOFREE(taxonomy,"Free taxonomy structure");
+		
+		return 0;
+	}
+	
+	return 1;
+}
+
+ecotx_t *eco_findtaxonatrank(ecotx_t *taxon,
+                                int32_t rankidx)
+{
+	ecotx_t *current_taxon;
+	ecotx_t *next_taxon;
+	
+	current_taxon = taxon;
+	next_taxon    = current_taxon->parent;
+	
+	while ((current_taxon!=next_taxon) &&  // I' am the root node
+		   (current_taxon->rank!=rankidx))
+		   {
+		   	current_taxon = next_taxon;
+		   	next_taxon    = current_taxon->parent;
+		   }
+		   
+	if (current_taxon->rank==rankidx)
+		return current_taxon;
+	else
+		return NULL;
+}
+
+/**
+ * Get back information concerning a taxon from a taxonomic id
+ * @param 	*taxonomy 	the taxonomy database
+ * @param	taxid		the taxonomic id 
+ * 
+ * @result	a ecotx_t structure containing the taxonimic information  
+ **/
+ecotx_t *eco_findtaxonbytaxid(ecotaxonomy_t *taxonomy, 
+							  int32_t taxid)
+{
+	ecotx_t    *current_taxon;
+	int32_t     taxoncount;
+	int32_t     i;
+	
+	taxoncount=taxonomy->taxons->count;
+	
+	for (current_taxon=taxonomy->taxons->taxon,
+	     i=0;
+	     i < taxoncount;
+	     i++,
+	     current_taxon++){
+	     if (current_taxon->taxid==taxid){
+	     	return current_taxon;
+	     }
+	 }
+	
+	return (ecotx_t*)NULL;	
+}
+
+/**
+ * Find out if taxon is son of other taxon (identified by its taxid)
+ * @param	*taxon son 		taxon
+ * @param	parent_taxid 	taxonomic id of the other taxon
+ * 
+ * @return 	1 is the other taxid math a parent taxid, else 0
+ **/
+int eco_isundertaxon(ecotx_t *taxon, 
+						int other_taxid)
+{
+	ecotx_t *next_parent;
+	
+	next_parent = taxon->parent;	
+	
+	while ( (other_taxid != next_parent->taxid) && 
+			(strcmp(next_parent->name, "root")) )
+	{
+		next_parent = next_parent->parent;
+	}
+		
+	if (other_taxid == next_parent->taxid)
+		return 1;
+	else
+		return 0;
+}
+
+ecotx_t *eco_getspecies(ecotx_t *taxon,
+						ecotaxonomy_t *taxonomy)
+{
+	static ecotaxonomy_t *tax=NULL;
+	static int32_t		 rankindex=-1;
+	
+	if (taxonomy && tax!=taxonomy)
+	{
+		rankindex = rank_index("species",taxonomy->ranks);
+		tax=taxonomy;
+	}
+		
+	if (!tax || rankindex < 0)
+		ECOERROR(ECO_ASSERT_ERROR,"No taxonomy defined");
+		
+	return eco_findtaxonatrank(taxon,rankindex);
+}
+
+ecotx_t *eco_getgenus(ecotx_t *taxon,
+						ecotaxonomy_t *taxonomy)
+{
+	static ecotaxonomy_t *tax=NULL;
+	static int32_t		 rankindex=-1;
+	
+	if (taxonomy && tax!=taxonomy)
+	{
+		rankindex = rank_index("genus",taxonomy->ranks);
+		tax=taxonomy;
+	}
+		
+	if (!tax || rankindex < 0)
+		ECOERROR(ECO_ASSERT_ERROR,"No taxonomy defined");
+		
+	return eco_findtaxonatrank(taxon,rankindex);
+}
+
+
+ecotx_t *eco_getfamily(ecotx_t *taxon,
+						ecotaxonomy_t *taxonomy)
+{
+	static ecotaxonomy_t *tax=NULL;
+	static int32_t		 rankindex=-1;
+	
+	if (taxonomy && tax!=taxonomy)
+	{
+		rankindex = rank_index("family",taxonomy->ranks);
+		tax=taxonomy;
+	}
+		
+	if (!tax || rankindex < 0)
+		ECOERROR(ECO_ASSERT_ERROR,"No taxonomy defined");
+		
+	return eco_findtaxonatrank(taxon,rankindex);
+}
+
+ecotx_t *eco_getkingdom(ecotx_t *taxon,
+						ecotaxonomy_t *taxonomy)
+{
+	static ecotaxonomy_t *tax=NULL;
+	static int32_t		 rankindex=-1;
+	
+	if (taxonomy && tax!=taxonomy)
+	{
+		rankindex = rank_index("kingdom",taxonomy->ranks);
+		tax=taxonomy;
+	}
+		
+	if (!tax || rankindex < 0)
+		ECOERROR(ECO_ASSERT_ERROR,"No taxonomy defined");
+		
+	return eco_findtaxonatrank(taxon,rankindex);
+}
+
+ecotx_t *eco_getsuperkingdom(ecotx_t *taxon,
+						ecotaxonomy_t *taxonomy)
+{
+	static ecotaxonomy_t *tax=NULL;
+	static int32_t		 rankindex=-1;
+	
+	if (taxonomy && tax!=taxonomy)
+	{
+		rankindex = rank_index("superkingdom",taxonomy->ranks);
+		tax=taxonomy;
+	}
+		
+	if (!tax || rankindex < 0)
+		ECOERROR(ECO_ASSERT_ERROR,"No taxonomy defined");
+		
+	return eco_findtaxonatrank(taxon,rankindex);
+}