14 Commits

25 changed files with 685 additions and 613 deletions

15
.gitignore vendored
View File

@ -1,16 +1,15 @@
# /src/
*.o
/src/*.P
/src/*.a
/src/ecoPrimer
# /src/libecoPCR/
/src/libecoPCR/*.P
/src/libecoPCR/*.a
# /src/libecoprimer/
/src/libecoPCR/*.o
/src/libecoprimer/*.P
/src/libecoprimer/*.a
# /src/libthermo/
/src/libecoprimer/*.o
/src/libthermo/*.P
/src/libthermo/*.a
/src/libthermo/*.o
phylonorway.*
src/ecoPrimers

View File

@ -1 +1 @@
0.4
0.5

View File

@ -28,7 +28,7 @@ all: $(EXEC)
# executable compilation and link
ecoPrimers: $(PRIMER_OBJ) $(LIBFILE)
$(CC) -g $(LDFLAGS) -O5 -m64 -o $@ $< $(LIBPATH) $(LIB)
$(CC) $(LDFLAGS) -O5 -o $@ $< $(LIBPATH) $(LIB)
########

View File

@ -20,7 +20,7 @@
#include"libthermo/thermostats.h"
#define VERSION "0.4"
#define VERSION "0.5"
/* TR: by default, statistics are made on species level*/
#define DEFAULTTAXONRANK "species"
@ -171,13 +171,9 @@ void printapair(int32_t index,ppair_t pair, poptions_t options)
bool_t good2=pair->p2->good;
bool_t goodtmp;
bool_t strand;
uint32_t i, j;
float temp;
CNNParams nnparams;
uint32_t i;
//nparam_InitParams(&nnparams, DEF_CONC_PRIMERS,DEF_CONC_SEQUENCES,DEF_SALT,SALT_METHOD_SANTALUCIA);
char *c;
const char *c;
char p1[32];
char p2[32];
@ -376,21 +372,20 @@ uint32_t filterandsortpairs(ppair_t* sortedpairs,uint32_t count, poptions_t opti
void printpairs (ppairtree_t pairs, poptions_t options,ecotaxonomy_t *taxonomy, pecodnadb_t seqdb)
{
ppair_t* sortedpairs;
ppair_t* index;
// ppair_t* index;
ppairlist_t pl;
size_t i,j;
size_t count;
char *taxon[]={"taxon","taxa"};
ecotx_t *current_taxon;
//pairset pair_sets;
pairset *pset = NULL;
//printf("Index\tPrimer1\tPrimer2\tGB\tInexampleCount\tOutexampleCount\tYule\tIntaxaCount\tOuttaxaCount\tCoverage\tSpecificity\tMinAmplifiedLength\tMaxAmplifiedLength\tAvgAmplifiedLength\n");
fprintf(stderr,"Total pair count : %d\n",pairs->count);
sortedpairs = ECOMALLOC(pairs->count*sizeof(ppair_t),"Cannot Allocate ordered pairs");
index=sortedpairs;
// index=sortedpairs;
pl=pairs->first;
j=0;
while(pl->next)
@ -475,10 +470,10 @@ void printpairs (ppairtree_t pairs, poptions_t options,ecotaxonomy_t *taxonomy,
if (options->filter_on_links)
{
fprintf (stderr, "Old size: %d, ", count);
fprintf (stderr, "Old size: %ld, ", count);
count = primers_changeSortedArray (&sortedpairs, count, options);
//count = primers_filterWithGivenLinks (&sortedpairs, count, options);
fprintf (stderr, "New size: %d\n", count);
fprintf (stderr, "New size: %ld\n", count);
if (count == 0)
{
@ -799,13 +794,17 @@ int main(int argc, char **argv)
}
}
if (errflag) {
ExitUsage(1);
}
options.pnparm = &nnparams;
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
if (options.salt < 0.01 || options.salt > 0.3) {//if salt value out of literature values
options.salt = DEF_SALT; //set to default
}
nparam_InitParams(&nnparams, DEF_CONC_PRIMERS,DEF_CONC_SEQUENCES,options.salt,options.saltmethod);
fprintf(stderr,"Reading taxonomy database ...");
@ -851,7 +850,7 @@ int main(int argc, char **argv)
fprintf(stderr,"\nIndexing words in sequences\n");
words = lookforStrictPrimer(seqdb,seqdbsize,insamples,&options);
fprintf(stderr,"\n Strict primer count : %d\n",words->size);
fprintf(stderr,"\n Strict primer count : %lld\n",words->size);
/*/TR Testing
fprintf(stderr,"\nReducing for debugging\n");
@ -871,7 +870,7 @@ int main(int argc, char **argv)
if (options.no_multi_match)
{
(void)filterMultiStrictPrimer(words);
fprintf(stderr,"\n Strict primer with single match count : %d\n",words->size);
fprintf(stderr,"\n Strict primer with single match count : %lld\n",words->size);
}
@ -921,7 +920,7 @@ pwordcount_t reduce_words_to_debug (pwordcount_t words, poptions_t options)
{
uint32_t i, k;
pwordcount_t new_words;
char *rwrd;
const char *rwrd;
char dwrd[20];
/*char *strict_words[DEBUG_WORDS_CNT] = {"GAGTCTCTGCACCTATCC", "GCAATCCTGAGCCAAATC", "ACCCCTAACCACAACTCA",
"TCCGAACCGACTGATGTT", "GAAGCTTGGGTGAAACTA", "GGAGAACCAGCTAGCTCT", "GCTGGTTCTCCCCGAAAT",
@ -981,7 +980,7 @@ pwordcount_t reduce_words_to_debug (pwordcount_t words, poptions_t options)
void print_wordwith_positions (primer_t prm, uint32_t seqdbsize, poptions_t options)
{
char *wrd;
const char *wrd;
uint32_t i, j;
char *twrd = "GCCTGTTTACCAAAAACA";

View File

@ -1,9 +1,9 @@
MACHINE=MAC_OS_X
LIBPATH= -LlibecoPCR -Llibecoprimer -Llibthermo
MACHINE=LINUX
LIBPATH= -LlibecoPCR -Llibecoprimer -Llibthermo -L/usr/local/lib
MAKEDEPEND = gcc -D$(MACHINE) -M $(CPPFLAGS) -o $*.d $<
CC=gcc
CFLAGS= -W -Wall -m64 -g
CFLAGS= -W -Wall -O5
#CFLAGS= -W -Wall -O5 -m64 -g
#CFLAGS= -W -Wall -O0 -m64 -g
#CFLAGS= -W -Wall -O5 -fast -g

View File

@ -2,7 +2,7 @@
#include <stdlib.h>
static int eco_log_malloc = 0;
static size_t eco_amount_malloc=0;
//static size_t eco_amount_malloc=0;
static size_t eco_chunk_malloc=0;
void eco_trace_memory_allocation()
@ -37,7 +37,7 @@ void *eco_malloc(int64_t chunksize,
if (eco_log_malloc)
fprintf(stderr,
"Memory segment located at %p of size %d is allocated (file : %s [%d])",
"Memory segment located at %p of size %lld is allocated (file : %s [%d])",
chunk,
chunksize,
filename,
@ -65,7 +65,7 @@ void *eco_realloc(void *chunk,
if (!newchunk)
{
fprintf(stderr,"Requested memory : %d\n",newsize);
fprintf(stderr,"Requested memory : %lld\n",newsize);
ecoError(ECO_MEM_ERROR,error_message,filename,line);
}
if (!chunk)
@ -73,7 +73,7 @@ void *eco_realloc(void *chunk,
if (eco_log_malloc)
fprintf(stderr,
"Old memory segment %p is reallocated at %p with a size of %d (file : %s [%d])",
"Old memory segment %p is reallocated at %p with a size of %lld (file : %s [%d])",
chunk,
newchunk,
newsize,

View File

@ -21,8 +21,8 @@ typedef struct {
int32_t taxid;
char AC[20];
int32_t DE_length;
int32_t SQ_length;
int32_t CSQ_length; /*what is this CSQ_length ? */
uint32_t SQ_length;
uint32_t CSQ_length; /*what is this CSQ_length ? */
char data[1];
@ -30,7 +30,7 @@ typedef struct {
typedef struct {
int32_t taxid;
int32_t SQ_length;
uint32_t SQ_length;
int32_t isexample;
char *AC;
char *DE;
@ -219,7 +219,7 @@ econameidx_t *read_nameidx(const char *filename,ecotaxonomy_t *taxonomy);
* @return pointer to a taxonomy index structure
*/
ecotxidx_t *read_taxonomyidx(const char *filename);
ecotxidx_t *read_taxonomyidx(const char *filename,const char *filename2);
ecotaxonomy_t *read_taxonomy(const char *prefix,int32_t readAlternativeName);

View File

@ -34,14 +34,11 @@ ecorankidx_t *read_rankidx(const char *filename)
int32_t rank_index(const char* label,ecorankidx_t* ranks)
{
char **rep;
fprintf(stderr,"Looking for rank -%s-... ",label);
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;
}

View File

@ -2,6 +2,7 @@
#include <stdlib.h>
#include <string.h>
#include <zlib.h>
#include <ctype.h>
#include <string.h>
#include <stdio.h>
@ -194,7 +195,7 @@ 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];
static char current_prefix[1025];
ecoseq_t *seq;
if (prefix)

View File

@ -10,23 +10,41 @@ static ecotx_t *readnext_ecotaxon(FILE *f,ecotx_t *taxon);
* @param pointer to the database (.tdx file)
* @return a ecotxidx_t structure
*/
ecotxidx_t *read_taxonomyidx(const char *filename)
ecotxidx_t *read_taxonomyidx(const char *filename,const char *filename2)
{
int32_t count;
int32_t count2;
FILE *f;
FILE *f2;
ecotxidx_t *index;
int32_t i;
f = open_ecorecorddb(filename,&count,1);
f = open_ecorecorddb(filename,&count,1);
f2 = open_ecorecorddb(filename2,&count2,0);
index = (ecotxidx_t*) ECOMALLOC(sizeof(ecotxidx_t) + sizeof(ecotx_t) * (count-1),
index = (ecotxidx_t*) ECOMALLOC(sizeof(ecotxidx_t) + sizeof(ecotx_t) * (count+count2-1),
"Allocate taxonomy");
index->count=count;
index->count=count+count2;
fprintf(stderr,"Reading %d taxa...\n",count);
for (i=0; i < count; i++){
readnext_ecotaxon(f,&(index->taxon[i]));
index->taxon[i].parent=index->taxon + (size_t)index->taxon[i].parent;
index->taxon[i].parent=index->taxon + (int32_t)index->taxon[i].parent;
}
if (count2>0)
fprintf(stderr,"Reading %d local taxa...\n",count2);
else
fprintf(stderr,"No local taxon\n");
for (i=0; i < count2; i++){
readnext_ecotaxon(f2,&(index->taxon[count+i]));
index->taxon[count+i].parent=index->taxon + (int32_t)index->taxon[count+i].parent;
}
return index;
}
@ -94,7 +112,7 @@ ecotx_t *readnext_ecotaxon(FILE *f,ecotx_t *taxon)
raw->taxid = swap_int32_t(raw->taxid);
}
taxon->parent = (ecotx_t*)(size_t)raw->parent;
taxon->parent = (ecotx_t*)raw->parent;
taxon->taxid = raw->taxid;
taxon->rank = raw->rank;
@ -111,6 +129,7 @@ ecotaxonomy_t *read_taxonomy(const char *prefix,int32_t readAlternativeName)
{
ecotaxonomy_t *tax;
char *filename;
char *filename2;
int buffsize;
tax = ECOMALLOC(sizeof(ecotaxonomy_t),
@ -120,14 +139,17 @@ ecotaxonomy_t *read_taxonomy(const char *prefix,int32_t readAlternativeName)
filename = ECOMALLOC(buffsize,
"Allocate filename");
filename2= ECOMALLOC(buffsize,
"Allocate filename");
snprintf(filename,buffsize,"%s.rdx",prefix);
tax->ranks = read_rankidx(filename);
snprintf(filename,buffsize,"%s.tdx",prefix);
snprintf(filename2,buffsize,"%s.ldx",prefix);
tax->taxons = read_taxonomyidx(filename);
tax->taxons = read_taxonomyidx(filename,filename2);
if (readAlternativeName)
{
@ -319,6 +341,9 @@ ecotx_t *eco_getsuperkingdom(ecotx_t *taxon,
if (taxonomy && tax!=taxonomy)
{
rankindex = rank_index("superkingdom",taxonomy->ranks);
if (rankindex < 0) {
rankindex = rank_index("domain",taxonomy->ranks);
}
tax=taxonomy;
}

View File

@ -1139,7 +1139,7 @@ void sets_by_SimulatedAnealing (pairset *pair_set,
{
pair_set = extend_set_randomly (NULL, &params, 3);
printf("\nStart Random seed set for Simulated :\n");
print_set_info (&pair_set, &params);
print_set_info (pair_set, &params);
}
min_spc = max_spc = pair_set->set_specificity;
min_cov = max_cov = pair_set->set_coverage;
@ -1566,7 +1566,8 @@ int32_t *addinset (int32_t *set, int32_t i, int32_t j, int32_t* slots, int32_t *
size_t primers_changeSortedArray (ppair_t ** pairs,
size_t sorted_count, poptions_t options)
{
int32_t i, j, k, l, total_links;
int32_t i, k, total_links;
uint32_t j;
int *owi;
int *iwi;
int allowedtaxa;
@ -1579,14 +1580,14 @@ size_t primers_changeSortedArray (ppair_t ** pairs,
idx_set = ECOMALLOC(slots*sizeof (int32_t),
"Could not allocate memory for index set.");
for (i=0; i<sorted_count; i++)
for (i=0; ((uint32_t)i)<sorted_count; i++)
{
owi = sortedpairs[i]->wellIdentifiedSeqs;
passed = FALSE;
for (j=0; j<sorted_count; j++)
{
if (i == j) continue;
if ((uint32_t)i == j) continue;
iwi = sortedpairs[j]->wellIdentifiedSeqs;
total_links = 0;
@ -1610,9 +1611,10 @@ size_t primers_changeSortedArray (ppair_t ** pairs,
if (options->max_links_percent > 0)
{
allowedtaxa = options->max_links_percent;
if (total_links > allowedtaxa)
if (total_links > allowedtaxa){
passed = TRUE;
break;
}
break;
}
else
if (!(total_links > 5 && total_links <= allowedtaxa))
@ -1630,7 +1632,7 @@ size_t primers_changeSortedArray (ppair_t ** pairs,
for (j=0; j<sorted_count; j++)
{
for (k=0; k<index; k++)
if (j == idx_set[k]) break;
if (j == (uint32_t)(idx_set[k])) break;
//need to remove this element
if (k == index)
{
@ -1711,9 +1713,9 @@ int32_t *addinset_withLinks (int32_t *set, int32_t i, int32_t* slots, int32_t *i
size_t primers_filterWithGivenLinks (ppair_t ** pairs,
size_t sorted_count, poptions_t options)
{
int32_t i, j, k;
int32_t i, k;
uint32_t j;
ppair_t *sortedpairs = *pairs;
bool_t passed;
int32_t *idx_set = NULL;
int32_t slots=50, index=0;
@ -1732,7 +1734,7 @@ size_t primers_filterWithGivenLinks (ppair_t ** pairs,
for (j=0; j<sorted_count; j++)
{
for (k=0; k<index; k++)
if (j == idx_set[k]) break;
if (j == (uint32_t) (idx_set[k])) break;
//need to remove this element
if (k == index)
{
@ -1755,7 +1757,7 @@ size_t primers_filterWithGivenLinks (ppair_t ** pairs,
else i=sorted_count;
for (j=0; j<i; j++)
for (j=0; j<(uint32_t)i; j++)
for (k=0; k<options->dbsize; k++)
if ((*pairs)[j]->coveredSeqs[k] == 1)
cov[k] = 1;

View File

@ -213,7 +213,6 @@ pprimercount_t ahoc_lookforStrictPrimers (pecodnadb_t database, uint32_t seqdbsi
char *base;
int8_t code;
uint32_t goodPrimers=0;
static int iii=0;
//inSequenceQuorum = (uint32_t)floor((float)exampleCount * options->sensitivity_quorum);

View File

@ -36,10 +36,10 @@ int32_t ManberNoErr(pecoseq_t pseq,ppattern_t pat,
ppatternParam_t param,
StackiPtr stkpos)
{
int32_t pos;
uint32_t pos;
uint32_t smask, r;
uint8_t *data;
int32_t end;
uint32_t end;
end = (size_t)(pseq->SQ_length);
@ -84,11 +84,11 @@ int32_t ManberSub(pecoseq_t pseq,ppattern_t pat,
StackiPtr stkpos)
{
int e, found;
int32_t pos;
uint32_t pos;
uint32_t smask, cmask, sindx;
uint32_t *pr, r[2 * MAX_PAT_ERR + 2];
uint8_t *data;
int32_t end;
uint32_t end;
end = (size_t)(pseq->SQ_length);

View File

@ -59,7 +59,7 @@ void encodeSequence(ecoseq_t *seq)
data = (uint8_t*)(seq->SQ);
cseq = seq->SQ;
for (i=0;i<seq->SQ_length;i++,data++,cseq++)
for (i=0;i<(int)(seq->SQ_length);i++,data++,cseq++)
{
*data = encoder[(IS_UPPER(*cseq) ? *cseq : 'Z') - 'A'];
}
@ -201,7 +201,7 @@ pprimercount_t lookforAproxPrimer(pecodnadb_t database, uint32_t seqdbsize,uint3
data[w].good = data[w].inexample >= inSequenceQuorum && data[w].outexample <= outSequenceQuorum;
goodPrimers+=data[w].good? 1:0;
fprintf(stderr,"Primers %5d/%d analyzed => sequence : %s in %d example and %d counterexample sequences \r",
fprintf(stderr,"Primers %5d/%lld analyzed => sequence : %s in %d example and %d counterexample sequences \r",
i+1,words->size,ecoUnhashWord(data[w].word,options->primer_length),
data[w].inexample,data[w].outexample);
@ -224,8 +224,8 @@ pprimercount_t lookforAproxPrimer(pecodnadb_t database, uint32_t seqdbsize,uint3
ECOFREE(data[w].reversePos,"Free direct count table");
}
fprintf(stderr,"\n\nOn %d analyzed primers %d respect quorum conditions\n",words->size,goodPrimers);
fprintf(stderr,"Conserved primers for further analysis : %d/%d\n",w,words->size);
fprintf(stderr,"\n\nOn %lld analyzed primers %d respect quorum conditions\n",words->size,goodPrimers);
fprintf(stderr,"Conserved primers for further analysis : %d/%lld\n",w,words->size);
primers = ECOMALLOC(sizeof(primercount_t),"Cannot allocate memory for primer table");
primers->primers=ECOREALLOC(data,

View File

@ -8,14 +8,6 @@
#ifndef HASHENCODER_H_
#define HASHENCODER_H_
static int8_t encoder[] = {0, // A
-1, // b
1, // C
-1,-1,-1, // d, e, f
2, // G
-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, // h,i,j,k,l,m,n,o,p,q,r,s
3,3, // T,U
-1,-1,-1,-1,-1}; // v,w,x,y,z
extern int8_t encoder[];
#endif /* HASHENCODER_H_ */

View File

@ -12,6 +12,16 @@ static int cmpword(const void *x,const void *y);
#include "hashencoder.h"
int8_t encoder[] = {0, // A
-1, // b
1, // C
-1,-1,-1, // d, e, f
2, // G
-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, // h,i,j,k,l,m,n,o,p,q,r,s
3,3, // T,U
-1,-1,-1,-1,-1}; // v,w,x,y,z
uint32_t ecoWordCount(uint32_t wordsize, uint32_t circular, ecoseq_t *seq)
{
uint32_t wordcount;
@ -257,6 +267,6 @@ uint32_t ecoFindWord(pwordcount_t table,word_t word)
char ecoComplementChar(char base)
{
return (base < 4)? !base & 3: 4;
return (base < 4)? ~base & 3: 4;
}

View File

@ -250,8 +250,7 @@ void CursiToTop(StackiPtr stki)
} /* end of CursiToTop */
void CursiToBottom(stki)
StackiPtr stki;
void CursiToBottom(StackiPtr stki)
{
stki->cursor = 0;

View File

@ -16,6 +16,7 @@
#include "ecotype.h"
#include <stdint.h>
/* ==================================================== */
/* Constantes de dimensionnement */
@ -85,5 +86,6 @@ bool_t SearchDownStacki (StackiPtr stki , int32_t sval );
bool_t BinSearchStacki (StackiPtr stki , int32_t sval );
bool_t SameStacki (StackiPtr stki1 , StackiPtr stki2 );
bool_t ReverseStacki (StackiPtr stki );
void CursiToBottom (StackiPtr stki );
#endif /* _H_libstki */

View File

@ -41,7 +41,7 @@ void ecomerge(pwordcount_t data,uint32_t s1,uint32_t s2,uint32_t remainingSeq,ui
uint32_t currentcount,tmpcount;
int same;
queue_t queue;
int nsame=0;
// int nsame=0;
uint32_t maxcount=0;
bool_t writed=TRUE;
@ -103,7 +103,7 @@ void ecomerge(pwordcount_t data,uint32_t s1,uint32_t s2,uint32_t remainingSeq,ui
* record if the two words in the both subtable are the same
*/
same = (source != S2) && (WORD(currentword) == WORD(merged.words[merged.read2]));
nsame+=same;
// nsame+=same;
// DEBUG_LOG("Merging : r1 = %d s1 = %d r2 = %d size = %d word = %s source=%u same=%u",merged.read1,s1,merged.read2-s1,merged.size,ecoUnhashWord(currentword,18),source,same)

View File

@ -205,8 +205,8 @@ pwordcount_t lookforStrictPrimer(pecodnadb_t database, uint32_t seqdbsize,
}
else
{
uint32_t s;
s = strictprimers->size;
// uint32_t s;
// s = strictprimers->size;
// DEBUG_LOG("stack size : %u",s);
addSeqToWordCountTable(strictprimers,options->primer_length,
options->circular,
@ -223,7 +223,7 @@ pwordcount_t lookforStrictPrimer(pecodnadb_t database, uint32_t seqdbsize,
logfile = fopen(logfilename,"a");
seconde = timeval_subtract(&(usage.ru_utime),&(start.ru_utime)) +
timeval_subtract(&(usage.ru_stime),&(start.ru_stime));
fprintf(logfile,"%d\t%llu\t%lu\t%8.3f\t%8.3e\n",i,
fprintf(logfile,"%d\t%llu\t%llu\t%8.3f\t%8.3e\n",i,
(long long unsigned)totallength,
strictprimers->size*(sizeof(int64_t)+sizeof(int32_t)),
seconde,seconde/(double)totallength);
@ -248,9 +248,9 @@ pwordcount_t lookforStrictPrimer(pecodnadb_t database, uint32_t seqdbsize,
sizeof(word_t)*strictprimers->size,
"Cannot reallocate strict primer table");
if (neededWords)
if (neededWords){
ECOFREE(neededWords,"Clean needed word table");
}
//TR: Somehow for some primers strictcount value is extremely large hence invalid
//we need to remove these primers from the list
j = strictprimers->size+1;

View File

@ -231,8 +231,8 @@ static int cmpamp(const void *ampf1, const void* ampf2)
char cd1;
char cd2;
int len = 0;
char *ch1;
char *ch2;
const char *ch1;
const char *ch2;
int incr1;
int incr2;

View File

@ -17,6 +17,15 @@
double forbidden_entropy;
char bpencoder[] = { 1, // A
0, // b
2, // C
0,0,0, // d, e, f
3, // G
0,0,0,0,0,0,0,0,0,0,0,0, // h,i,j,k,l,m,n,o,p,q,r,s
4,0, // T,U
0,0,0,0,0}; // v,w,x,y,z
double nparam_GetInitialEntropy(PNNParams nparm)
{
@ -470,7 +479,7 @@ int nparam_CountGCContent(char * seq ) {
return count;
}
void nparam_CleanSeq (char* inseq, char* outseq, int len)
void nparam_CleanSeq (const char* inseq, char* outseq, int len)
{
int seqlen = strlen (inseq);
int i, j;
@ -508,7 +517,7 @@ void nparam_CleanSeq (char* inseq, char* outseq, int len)
}
//Calculate TM for given sequence against its complement
double nparam_CalcSelfTM(PNNParams nparm, char* seq, int len)
double nparam_CalcSelfTM(PNNParams nparm, const char* seq, size_t len)
{
double thedH = 0;
//double thedS = nparam_GetInitialEntropy(nparm);
@ -520,7 +529,7 @@ double nparam_CalcSelfTM(PNNParams nparm, char* seq, int len)
char c4;
unsigned int i;
char nseq[50];
char *useq = seq;
const char *useq = seq;
nparam_CleanSeq (seq, nseq, len);
useq = nseq;
@ -533,7 +542,7 @@ double nparam_CalcSelfTM(PNNParams nparm, char* seq, int len)
c4 = GETNUMCODE(useq[i]);
thedH += nparm->dH[c3][c4][c1][c2];//nparam_GetEnthalpy(nparm, c3,c4,c1,c2);
thedH += nparm->dH[(uint8_t) c3][(uint8_t)c4][(uint8_t)c1][(uint8_t)c2];//nparam_GetEnthalpy(nparm, c3,c4,c1,c2);
thedS += nparam_GetEntropy(nparm, c3,c4,c1,c2);
}
//printf("------------------\n");
@ -543,7 +552,7 @@ double nparam_CalcSelfTM(PNNParams nparm, char* seq, int len)
return mtemp;
}
double nparam_CalcTwoTM(PNNParams nparm, char* seq1, char* seq2, int len)
double nparam_CalcTwoTM(PNNParams nparm, const char* seq1, const char* seq2, size_t len)
{
double thedH = 0;
//double thedS = nparam_GetInitialEntropy(nparm);
@ -575,7 +584,7 @@ double nparam_CalcTwoTM(PNNParams nparm, char* seq1, char* seq2, int len)
//fprintf (stderr,"Primer : %s %f %f %d %d, %d %d %f\n",useq,thedH,thedS,(int)c3,(int)c4,(int)c1,(int)c2,nparam_GetEnthalpy(nparm, c3,c4,c1,c2));
thedH += nparm->dH[c3][c4][c1][c2];//nparam_GetEnthalpy(nparm, c3,c4,c1,c2);
thedH += nparm->dH[(uint8_t)c3][(uint8_t)c4][(uint8_t)c1][(uint8_t)c2];//nparam_GetEnthalpy(nparm, c3,c4,c1,c2);
thedS += nparam_GetEntropy(nparm, c3,c4,c1,c2);
}
//fprintf(stderr,"------------------\n");

View File

@ -13,11 +13,13 @@
#include <math.h>
#include <string.h>
#include <stdint.h>
//#include "../libecoprimer/ecoprimer.h"
// following defines to simplify coding...
#define ndH(a,b,c,d) nparm->dH[a][b][c][d]
#define ndS(a,b,c,d) nparm->dS[a][b][c][d]
#define ndH(a,b,c,d) nparm->dH[(uint8_t) a][(uint8_t) b][(uint8_t) c][(uint8_t) d]
#define ndS(a,b,c,d) nparm->dS[(uint8_t) a][(uint8_t) b][(uint8_t) c][(uint8_t) d]
#define forbidden_enthalpy 1000000000000000000.0f
#define R 1.987f
#define SALT_METHOD_SANTALUCIA 1
@ -33,14 +35,7 @@
extern double forbidden_entropy;
static char bpencoder[] = { 1, // A
0, // b
2, // C
0,0,0, // d, e, f
3, // G
0,0,0,0,0,0,0,0,0,0,0,0, // h,i,j,k,l,m,n,o,p,q,r,s
4,0, // T,U
0,0,0,0,0}; // v,w,x,y,z
extern char bpencoder[]; // v,w,x,y,z
typedef struct CNNParams_st
@ -62,8 +57,8 @@ int nparam_CountGCContent(char * seq );
double nparam_GetEntropy(PNNParams nparm, char x0, char x1, char y0, char y1);
double nparam_GetEnthalpy(PNNParams nparm, char x0, char x1, char y0, char y1);
double nparam_CalcTM(double entropy,double enthalpy);
double nparam_CalcSelfTM(PNNParams nparm, char* seq, int len);
double nparam_CalcTwoTM(PNNParams nparm, char* seq1, char* seq2, int len);
double nparam_CalcSelfTM(PNNParams nparm, const char* seq, size_t len);
double nparam_CalcTwoTM(PNNParams nparm, const char* seq1, const char* seq2, size_t len);
double nparam_GetInitialEntropy(PNNParams nparm) ;
double calculateMeltingTemperatureBasic (char * seq);

View File

@ -28,17 +28,15 @@ word_t extractSite(char* sequence, size_t begin, size_t length, bool_t strand)
void getThermoProperties (ppair_t* pairs, size_t count, poptions_t options)
{
size_t i, j,k,l;
size_t i, j;
uint32_t bp1,bp2;
uint32_t ep1,ep2;
word_t w1;
word_t w2;
bool_t strand;
char *sq,*sq1,*sq2,*c;
char *sq;
char prmrd[50];
char prmrr[50];
char sqsite[50];
double mtemp;
for (i = 0; i < count; i++)

File diff suppressed because it is too large Load Diff