From e79738e17078301492801509039de40ab068b69f Mon Sep 17 00:00:00 2001 From: Eric Coissac Date: Mon, 13 Jul 2009 09:26:19 +0000 Subject: [PATCH] commit -m "Cleaned code for thermodynamics properties and added Melting Temperature for approximate version of strict repeats. Also cleaned printing code. git-svn-id: https://www.grenoble.prabi.fr/svn/LECASofts/ecoPrimers/trunk@223 60f365c0-8329-0410-b2a4-ec073aeeaa1d --- src/ecoprimer.c | 134 ++++----- src/libecoprimer/ecoprimer.h | 8 +- src/libecoprimer/pairs.c | 37 ++- src/libecoprimer/readdnadb.c | 15 + src/libthermo/Makefile | 7 +- src/libthermo/libthermo.a | Bin 37528 -> 13184 bytes src/libthermo/nnparams.c | 564 ++++++++++------------------------- src/libthermo/nnparams.h | 104 +++---- 8 files changed, 305 insertions(+), 564 deletions(-) diff --git a/src/ecoprimer.c b/src/ecoprimer.c index 3f481b0..75ed821 100644 --- a/src/ecoprimer.c +++ b/src/ecoprimer.c @@ -6,7 +6,6 @@ */ #include "libecoprimer/ecoprimer.h" -#include "libthermo/libthermo.h" #include #include #include @@ -15,6 +14,9 @@ #include #include #include +#include"libthermo/nnparams.h" +#include"libthermo/thermostats.h" + #define VERSION "0.3" /* TR: by default, statistics are made on species level*/ @@ -153,37 +155,10 @@ void initoptions(poptions_t options) options->r=0; options->g=0; options->no_multi_match=FALSE; + options->pnparm = NULL; strcpy(options->taxonrank, DEFAULTTAXONRANK); /*taxon level for results, species by default*/ } -void printcurrenttime () -{ - time_t now; - struct tm *ts; - char buf[80]; - - /* Get the current time */ - now = time(NULL); - - /* Format and print the time, "ddd yyyy-mm-dd hh:mm:ss zzz" */ - ts = localtime(&now); - strftime(buf, sizeof(buf), "%a %Y-%m-%d %H:%M:%S %Z", ts); - fprintf(stderr,"#%d#, %s\n",(int)now, buf); - } - -void printcurrenttimeinmilli() -{ - struct timeval tv; - struct timezone tz; - struct tm *tm; - gettimeofday(&tv, &tz); - tm=localtime(&tv.tv_sec); - fprintf(stderr, " %d:%02d:%02d %d %d \n", tm->tm_hour, tm->tm_min, - tm->tm_sec, tv.tv_usec, tv.tv_usec/1000); -} - -/*TR: Added*/ - void printapair(int32_t index,ppair_t pair, poptions_t options) { bool_t asdirect1=pair->asdirect1; @@ -198,6 +173,9 @@ void printapair(int32_t index,ppair_t pair, poptions_t options) bool_t strand; uint32_t i; float temp; + CNNParams nnparams; + + //nparam_InitParams(&nnparams, DEF_CONC_PRIMERS,DEF_CONC_SEQUENCES,DEF_SALT,SALT_METHOD_SANTALUCIA); char *c; char p1[32]; @@ -225,6 +203,7 @@ void printapair(int32_t index,ppair_t pair, poptions_t options) good2=goodtmp; } + //print serial number printf("%6d\t",index); c = ecoUnhashWord(w1,options->primer_length); @@ -232,57 +211,67 @@ void printapair(int32_t index,ppair_t pair, poptions_t options) c = ecoUnhashWord(w2,options->primer_length); strcpy (p2, c); + //print primer1 printf("%s\t", p1); + + //print primer2 printf("%s", p2); - /*For thermo: first word should be 5'-3'*/ - /*if (!asdirect1) - { - wtmp = ecoComplementWord(w1,options->primer_length); - c = ecoUnhashWord(wtmp,options->primer_length); - strcpy (p1, c); - }*/ - /*For thermo: first word should be 5'-3'*/ - /*if (asdirect2) - { - wtmp = ecoComplementWord(w2,options->primer_length); - c = ecoUnhashWord(wtmp,options->primer_length); - strcpy (p2, c); - }*/ + //print primer1 melting temperature + printf ("\t%4.3f", pair->p1temp); - /*temp = calculateMeltingTemperature (p1, p2);*/ - //float temp = calculateMeltingTemperature ("CTGTTTACCAAAAACATC", "GGTCTGAACTCAGATCAC"); - - temp = calculateMeltingTemperatureBasic(p1); - printf ("\t%4.3f", temp); - - temp = calculateMeltingTemperatureBasic(p2); - printf ("\t%4.3f", temp); - - printf ("\t%d",countGCContent(p1)); - - printf ("\t%d",countGCContent(p2)); + //print minimum melting temperature of approximate versions of primer1 + printf ("\t%4.3f", pair->p1mintemp); + + //print primer2 melting temperature + printf ("\t%4.3f", pair->p2temp); + + //print minimum melting temperature of approximate versions of primer2 + printf ("\t%4.3f", pair->p2mintemp); + + //print gc contents of primer1 + printf ("\t%d",nparam_CountGCContent(p1)); + + //print gc contents of primer2 + printf ("\t%d",nparam_CountGCContent(p2)); + //print good/bad pair indicator printf("\t%c%c", "bG"[(int)good1],"bG"[(int)good2]); - + //print inexample count printf("\t%d", pair->inexample); + + //print out example count printf("\t%d", pair->outexample); + + //print yule printf("\t%4.3f", pair->yule); + //print in taxa count printf("\t%d", pair->intaxa); + + //print out taxa count printf("\t%d", pair->outtaxa); + + //print coverage printf("\t%4.3f", (float)pair->bc); + //print well identified taxa count printf("\t%d", pair->intaxa - pair->notwellidentifiedtaxa); + //print specificity printf("\t%4.3f", pair->bs); - + //print min amplifia lenght printf("\t%d", pair->mind); + + //print max amplifia lenght printf("\t%d", pair->maxd); + + //print average amplifia lenght printf("\t%3.2f", (float)pair->sumd/pair->inexample); + //print amplifia information about reference sequence if specified if (options->refseq && pair->refsequence >=0) { printf("\t%s:",options->reference); @@ -351,31 +340,21 @@ uint32_t filterandsortpairs(ppair_t* sortedpairs,uint32_t count, poptions_t opti qfp = (float)sortedpairs[i]->outexample/options->outsamples; else qfp=0.0; - sortedpairs[i]->quorumin = q; sortedpairs[i]->quorumout = qfp; sortedpairs[i]->yule = q - qfp; - - sortedpairs[j]=sortedpairs[i]; - - if (q > options->sensitivity_quorum && qfp < options->false_positive_quorum) { (void)taxonomycoverage(sortedpairs[j],options); taxonomyspecificity(sortedpairs[j]); -// fprintf(stderr,"%4.3f %4.3f %4.3f\n",sortedpairs[j]->yule , sortedpairs[j]->bs,sortedpairs[j]->bc); - j++; } } - - qsort(sortedpairs,j,sizeof(ppair_t),cmpprintedpairs); - return j; } @@ -410,7 +389,8 @@ void printpairs (ppairtree_t pairs, poptions_t options,ecotaxonomy_t *taxonomy) sortedpairs[j]=pl->pairs+i; count=filterandsortpairs(sortedpairs,pairs->count,options); - + getThermoProperties(sortedpairs, count, options); + fprintf(stderr,"Total good pair count : %u\n",(uint32_t)count); printf("#\n"); @@ -477,7 +457,6 @@ void printpairs (ppairtree_t pairs, poptions_t options,ecotaxonomy_t *taxonomy) } - /*updateseqparams: This function counts the insample and outsample sequences * and with each sequences adds a tag of the taxon to which the sequence beongs*/ @@ -556,6 +535,12 @@ int main(int argc, char **argv) //float temp = calculateMeltingTemperatureBasic ("CTGTTTACCAAAAACATC"); //printf ("temp = %f\n", temp); //return 0; + //char *t = "ACGT"; + //printf ("\nGETNUMCODE: A=%d, C=%d, G=%d, T=%d\n", GETNUMCODE(t[0]), GETNUMCODE('C'), GETNUMCODE('G'), GETNUMCODE('T')); + //printf ("\nGETREVCODE: A=%d, C=%d, G=%d, T=%d\n", GETREVCODE(t[0]), GETREVCODE('C'), GETREVCODE('G'), GETREVCODE('T')); + //return 0; + CNNParams nnparams; + nparam_InitParams(&nnparams, DEF_CONC_PRIMERS,DEF_CONC_SEQUENCES,DEF_SALT,SALT_METHOD_SANTALUCIA); initoptions(&options); @@ -701,7 +686,7 @@ int main(int argc, char **argv) } } - + options.pnparm = &nnparams; fprintf(stderr,"Reading taxonomy database ..."); taxonomy = read_taxonomy(options.prefix,0); fprintf(stderr,"Ok\n"); @@ -741,7 +726,7 @@ int main(int argc, char **argv) words = lookforStrictPrimer(seqdb,seqdbsize,insamples,&options); fprintf(stderr,"\n Strict primer count : %d\n",words->size); - + // options.filtering=FALSE; // words2= lookforStrictPrimer(seqdb,seqdbsize,insamples,&options); // fprintf(stderr,"\n Strict primer count : %d\n",words2->size); @@ -789,17 +774,10 @@ int main(int argc, char **argv) fprintf(stderr,"\n"); - /*TR: Added*/ pairs = buildPrimerPairs(seqdb, seqdbsize, primers, &options); - //fprintf(stderr,"buildPrimerPairs finished\n"); - // setoktaxforspecificity (&pairs); - printpairs (pairs, &options,taxonomy); // closlibman (); - - //ECOFREE(pairs.pairs,"Free pairs table"); - return 0; } diff --git a/src/libecoprimer/ecoprimer.h b/src/libecoprimer/ecoprimer.h index 0cf3198..f226a01 100644 --- a/src/libecoprimer/ecoprimer.h +++ b/src/libecoprimer/ecoprimer.h @@ -13,6 +13,7 @@ #include #include "ecotype.h" #include "../libecoPCR/ecoPCR.h" +#include "../libthermo/nnparams.h" #include "apat.h" #define DEBUG @@ -190,7 +191,10 @@ typedef struct { // // uint32_t oktaxoncount; uint32_t curseqid; - + float p1temp; //strict primer1 melting temperature + float p1mintemp; //approx primer1 minimum melting temperature + float p2temp; //strict primer2 melting temperature + float p2mintemp; //approx primer2 minimum melting temperature } pair_t, *ppair_t; /*TR: Added*/ @@ -274,6 +278,7 @@ typedef struct { int32_t outsamples; int32_t intaxa; int32_t outtaxa; + PNNParams pnparm; } options_t, *poptions_t; typedef ecoseq_t **pecodnadb_t; @@ -338,5 +343,6 @@ void taxonomyspecificity (ppair_t pair); int32_t *filteringSeq(pecodnadb_t database, uint32_t seqdbsize, uint32_t exampleCount,poptions_t options,uint32_t *size,int32_t sequenceQuorum); +void printSeqTest(pecodnadb_t seqdb,uint32_t seqdbsize); #endif /* EPSORT_H_ */ diff --git a/src/libecoprimer/pairs.c b/src/libecoprimer/pairs.c index 8f2df4c..24b6dde 100644 --- a/src/libecoprimer/pairs.c +++ b/src/libecoprimer/pairs.c @@ -156,16 +156,14 @@ ppairtree_t buildPrimerPairs(pecodnadb_t seqdb,uint32_t seqdbsize,pprimercount_t { uint32_t i; ppairtree_t primerpairs; - + primerpairs = initpairtree(NULL); for (i=0; i < seqdbsize; i++) { buildPrimerPairsForOneSeq(i, seqdb, primers, primerpairs, options); } - return primerpairs; - } #define DMAX (2000000000) @@ -187,7 +185,11 @@ static void buildPrimerPairsForOneSeq(uint32_t seqid, bool_t bswp; size_t distance; bool_t strand; - + //char prmr[50]; + //float mtemp; + + //prmr[options->primer_length] = '\0'; + for (i=0;i < primers->size; i++) { matchcount+=primers->primers[i].directCount[seqid]; @@ -272,10 +274,12 @@ static void buildPrimerPairsForOneSeq(uint32_t seqid, current.outexample=0; current.curseqid = 0; current.refsequence=-1; - + //current.p1temp = 100; + //current.p1mintemp = 100; + //current.p2temp = 100; + //current.p2mintemp = 100; // Standardize the pair - strand = current.p2->word > current.p1->word; if (!strand) { @@ -309,6 +313,10 @@ static void buildPrimerPairsForOneSeq(uint32_t seqid, //else // pcurrent->outexample++; + //for each pair we save current sequence id in the pair + //when we see this pair for the first time in currnet sequence + //because we want to increment inexample & outexample count + //only once for one sequence if (pcurrent->curseqid != (seqid+1)) { if (seqdb[seqid]->isexample) @@ -359,6 +367,23 @@ static void buildPrimerPairsForOneSeq(uint32_t seqid, else pcurrent->pcr.amplifias[pcurrent->pcr.ampcount].amplifia= seqdb[seqid]->SQ + matches[j].position - 1 ; + + /*strncpy (prmr, seqdb[seqid]->SQ + matches[i].position, options->primer_length); + mtemp = nparam_CalcSelfTM (options->pnparm, prmr, options->primer_length) - 273.0; + if (mtemp < pcurrent->p1mintemp) + pcurrent->p1mintemp = mtemp; + //fprintf (stderr, "prmr1: %s\n", seqdb[seqid]->SQ); + strncpy (prmr, seqdb[seqid]->SQ + matches[j].position, options->primer_length); + mtemp = nparam_CalcSelfTM (options->pnparm, prmr, options->primer_length) - 273.0; + if (mtemp < pcurrent->p2mintemp) + pcurrent->p2mintemp = mtemp; + //fprintf (stderr, "prmr2: %s\n", prmr); + + if (pcurrent->p1temp == 100) + pcurrent->p1temp = nparam_CalcSelfTM (options->pnparm, ecoUnhashWord(pcurrent->p1->word, options->primer_length), 0) - 273.0; + if (pcurrent->p2temp == 100) + pcurrent->p2temp = nparam_CalcSelfTM (options->pnparm, ecoUnhashWord(pcurrent->p2->word, options->primer_length), 0) - 273.0; + */ pcurrent->pcr.ampcount++; // fprintf(stderr,"%c%c W1 : %s direct : %c", // "bG"[(int)pcurrent->p1->good], diff --git a/src/libecoprimer/readdnadb.c b/src/libecoprimer/readdnadb.c index a148510..98867dd 100644 --- a/src/libecoprimer/readdnadb.c +++ b/src/libecoprimer/readdnadb.c @@ -33,3 +33,18 @@ pecodnadb_t readdnadb(const char *name, uint32_t *size) return db; } + + +void printSeqTest(pecodnadb_t seqdb,uint32_t seqdbsize) +{ + uint32_t i; + char ch[11]; + ch [10] = '\0'; + + for (i=0; i < seqdbsize; i++) + { + strncpy (ch, seqdb[i]->SQ, 10); + fprintf (stderr, "seq %d = %s\n", i, ch); + } + exit (0); +} \ No newline at end of file diff --git a/src/libthermo/Makefile b/src/libthermo/Makefile index fefcbc7..60ae686 100644 --- a/src/libthermo/Makefile +++ b/src/libthermo/Makefile @@ -1,9 +1,6 @@ -SOURCES = dinkelbach.c \ - galign.c \ - libthermo.c \ - nnparams.c \ - thermalign.c +SOURCES = nnparams.c \ + thermostats.c SRCS=$(SOURCES) diff --git a/src/libthermo/libthermo.a b/src/libthermo/libthermo.a index 00900a23c85af0ed54c94872c3961baae6edd17e..cbedd2698d87b4f143a4666f7802869f42a7fabd 100644 GIT binary patch delta 6208 zcmc&&3s6+o8NQ1Pq5@|zfEXn#37RFtYZVoY8TZO=y)@8b1W|$@!9@~ztbhs5nC=?b z$j!1!=Om5O*r^?-9W$LIHrB*+8nR#lY#uG4NsLWQOdnH;c^HUUXpD!o6U=|@n^PXXD`jkF`E}J z$}#_Z#r(2jNpGD7NeWOmLz1os5WPg8U0|uedVwB+T;O8@-;^ZYYfB1WRo!xD4UaOW zPp@}YH?DPVnTrz0=tOTj>UCmOe|GDVlXrj=qRR zS6+nGn!PlaFUZfD-5myVGx?CDxwtH0JSFnllww*6+F=r-P14iRrz?V6W3uAcC6i3o zGE!%OB&G3Fix$Uk5|V3)i)Mm9BRtRJ=56NKokD8XWRjH5)9fo_4+w6X;1=^8_7wg= z%ysdJLNQfH2H|ZQe>NtG-^8vsu(8!-Vz{de$qElb9&W62&UNeb&3#}&jitoU>8R{T)~ zia&0+w1D*lz!+ByL0F;p3E00-41Mruf7m?frWA(#vKmzUMw_gT!s;B<+S>!2Ab z?a!N-rT2EZB6yr7+a%fZML2t?&(;xQzYcT^wI()g41u*UMii74_W*NfZ(K*`S8{+MR(X1jmU1tFK;~N)-v&`sV5aS)@^%Uz1ym}f8&q1Lo|j7s z)UVd6f5Q!l)JFEcC945Da~!t@9Nqxia|8BVq527Ibwsr%3d7%@cY2i&&nu{PI+FAxdSeC*oTVQW~B*4qVl^ zu4V^Jz3D&LIS`X@fA2j%IeUr%^(hjQJVJ`fRQlWdI0vL~HDUa}Okgq((gc*6fzSZX zQFlqLc98n>USi%c?Hmmh`Zsoy2>axu&!DvT8vWg&ix)3`Hk8u0>mMGp)-&o5EBmnb zjOnvJj|UHO`;g+7!uw{>v15}Q>K@m!|s-P=&${Z ztn6b?#YK|X0LspKUbIP%Gf-Me*>WRv`G}y3yA(3k`j;dJg5o_?XPr*FR`R_k?bMYQiv#=^H?5m4xD- zZ-thrYd_Tx%zJ$(6vD$^d-}UOhV{pPED8uC+mlLR(VvH+K$xs3SjYK`p%7ISZ0}k0 z;IZ7EYSIIPM|F^bLz!bUVm4VUXtC2H<{g5O{-hMo%gl`T4A_6o)0yRsv#hq8Wp3|6 zm+jt+Nhm&wC4ygG)TQ{Twx<3Ih3elGU+yssLAm6iO<0PSPLOCcLfQ2^c4~eR3O40{@9L{_;o~7Qx(>*T-08 zX&1{JWxMa8X=lE{yw6a_w-f&(kK+&CL^DVfd8)H7HKZH8z4G#(ZChtmZUnUw*K?BO z0eV(W;S0^vE7L0W_MNxH+J4EOan3P;?FqaXT3NpN=68SUKe3Vz}+hs`;Q=6oZG-{nW?$to|zWjy~WJetEStQi%}#az1qZ>?hiA@pc9RU znAURePm`X&EudD6plOYwsq4q16CGiscb;}**EIfq?NXjtXW|o1#jt_LE>DhK0nt<^ zPx5VbNql_f?AY|nB)j;l%d^7X2X$Xk<-f{$p8x!9NB)gsT?gROgFrEkfLqInv+<^`$tP}VX zf$>0+Ur5vqew3i&d3Y!ODIhsoizGtyN`b$_mcScd1{b&;Nd6muaFtO4gsWSM!cg^* zz##;bT*=W9EQ63f*-XB2xpaqGa#=JazgNM6XMMyCSAx0!PgElf>KEb z`Wndp1*8cOf**CP0-Z=G90TI^$NW(!5P~1&0pc`N5;Dkl0?7d(_)+xv0eiWU5Ko<$ zLLn!N;YC~LxJjFCJl=Gp-odnDz?gh3A0ZlnNi;(N?PH=57@Jumr9S}}c;|4}s2xO~ zLH=fnPQMKxxd%6a+Mtj}Ua!*%HkhVx7y@cVLgAcQBc(?H}+m!NmaG3 z<~ti(Yk5}pwY;inHs9S8%eOWq@I>?vHfDeAF|f)kK0tLf&AaOkm$R{^2|lF4eCx7h zl0&Jg#VAfIrYp@v``(|PmS*EC8nsO`TeVMFEv?n9EtyT+8ZmG70ScL!eCwa1at+$& zs0r*{Y5@_^NtzNTUDLGBQIAZ~qWf?Re=+=2#NycfLNa-GDtOBD6z${D0U=oreQKMa zf`lMha5{|I$0G$j8siwSP2v~r2{Xusl#KijEQ`A$jQlh+x@Ga^^SWHJ7B5|Gig%oK zG(lJKJGu+~TOv_+C<6}4ZTQS|vnZS6q)tpVfVy?~lEE2xCjHg|pFMSA zz~;gsUM2c%dcQDr!Y5nsmBafS##4L_kCNx5Y21CIPj2n^X#AGo%9@%tUol%k iBfqJ^Rqt#-eW`0TI2&qQJESXl%yiai*>J|CbogIwaHvoK literal 37528 zcmeHw4SZD9nfDEuK-3{OwBX_{=%CRCB}fD`Xw`N~j=Zc2x%O+etY2xLRu9)a@ zO&pgEf+S~;8+ZAH2`*Rm1pa53b!J)F*qdh-dAyTdH{WRW{L%|vaS1b5a1-B{qf`Z)&zNE@~Tfz4f z`~dWhg_?FovV<#R2hzoL4!f8N3cTG^ems}_{y&8{l1shaJtD4u;sRfR(4SNOdP{BvhlExI=b zQ@zlRm|R^^;ayNsedoP-(Nt074N+uKMU{Wvf;)>V<}b1(wPCHa1rYhVh1GK^s#Rtp zZuVCr^J0amM78pbs(E)>HDEPk5Q~sJBEP*_7J=1J_}O#jcotSw&91%~3{gb6r4E~L zR9H`ULm423;g+cc?%c# zE134q3V+_h`HQM5=8I@`#r%bLSHuXrcy?7$g@5kCIa(Pmk!8W6+10b>ml;)P@8 zoNb&R6T7*hsxmHtJ-RSf62&&gDq?WbE}|&5+Pv9Su%Q?V(NifQYO4}YTQmo@7l(vm z?yjizPpO$7Nf)OKm3!XeqIrww&-T;0{c|g-=O>u6#KuT&M@jJ_n5gku{E#z>gzff` z6t)||+N{Q&*oKazrZOy3m-&f_*t$f~*z<~6GiC7-4Xode?r(g*d{RH-d(NcGtnocM zLX5?zB*U*xzI$qrrX^u0(J*XjgvWT{!Vlpb{8lpFPY7X*V@Oxfii4-dpMyk%vNC_g z62H8b@nY|c;%MBu8foLlya(`aD2Wa&KwyQkvZ{F;->N1^Oj%hhKbPX;R2jyBk3tCy z>uC=WemN2w&u^T{oT@Up!OzYE(U?$HRx!7%5*?sIqT~7HD1IwdE}!DJ2=gV&I7`#! z+8HwWaQvDJ2$G~RdkFp?{&_BSJ3e$Y%6+Ux%s?U6@bB}(WhWJ~? zG;j91MzJ^h$jv_U6JPB+IUndlydA#AYmbum2YR}<1OEnhqi|BX`0p_P9l*b$LzVh@ z9!-<~eU0BwHG21ojJ)1Ofg>@pI@~xtwecQ@!0BUCRx742>nj3CD)K+iO!j1%5aueD z)oe*L4dtg7H3K*8fS%d*fnGvkZVgp}r4QNfuGFU=WN|O+V|VKh{2Le|Htmo;_JIDt zCKY>`ELo9hhkEjQdwZb|{lV|y!dw5Ye$8AHpU6rEJOl7U6q5=_+XOl5ARr3N>I7uEO$t(0WEplpAjYk% zy8t0m*7blWaM%&>Cu{>`zNZ1n{|A7~_aGq4xf_t>%mrjQtPkmz0x}=_Dr7As1iv({ zfe}gwfp@w9k-n4=MW>AbL?|Hy-uVh(8WLLL2hwOd|EchEKppg{NX&3Hex#?1n5OM! z9Q1&q=Ob>*t~JEP*T=<6;^O}h7r#6%eoFK(OXzD-PTKO7rqds)Lgs%txE`x0+VO0utus^2K5iPZ~n} zyl3%o396?3kp^lD$}A(8$h_RoJpZARi*q7_<8MSikrDjoU||oRH!q8uRb-f|GNywc zFFViyDFXU1yLKc{UaETTxtewfuSa1_1QO{Y^D^`E2bF0f@k_(+%LdNFYW7X_HQt|d zlngLvyD`C}-qzW}sC}MA`g?|(H)dikcX`bGG02W=KAF&(*KaYV6_}qH#Wy!j&kPR= z+_)zk&S&Fl;e3o`1VXL#2k&(mjovZ!A$|20pv7cnlyZIGu}bx3AD6}J=LMLxp!|)3 zMsJSM+wKdFxUDdl!JB#Qy0F7VNYZs-k4qFL(uF+^QDi~meLhh)VBMG<9qVE-hCamO zaP)yo^-ikEylN&*3(cBihWVjkzFKJBcOoz7nUm}dE=wvj=O6W%w;$TYdSQj<#zN1~ zH-~*@QAc6<(L(btqERH-%*(M5zS{5SXg)K~?O))l5BcZ$8oeDJvph5VIR1}=ePZ#m z7{-_9D<~IQh9w9s7nFw!_016HO~II!jO*;W6*k^? z!fxfh;MHtKk-Qa^?G;zi_7&^a=|;rDQKuzW-7`jCJqshCVIDVXPbKLOO+hzB1G);z z+x4x515O#i$-PEnO{Q7D1T1~!0@NSZw>k#s51oUY`mT_^tJ7#@cMjM4k%(qxdLJ{3 zHihdKAx^1!xTB*wf(6EQ=qqMHyD(7N3~W2y8wh&b$xw5<5p)dlt)n%pP;@M@(~DLm z(48xYUuQH<1Cyd0ECPd`n*<9_(h1)>eXL`-Ubg_nGz#I2tGtzCKNe9~*J3H#tTzlo1k@s?_I5-K9%BNT zJV_d}ds{?nGx{XCtoTVpNU7KkOHxroF$yE6glev77`` zUiB~%SZaPK6AN_oQ7RohYPi1DGvL(9`elg9ECr1`2 zs4Y9|En80WCA~<43qVWG6~3ewsPR>ZOr`~U>u_~FuZo(8V|#d4;wt3^+}*&Eneuwk zN^A5?*f_gAZqS<8LJfOF#JfgsI0rVt)AyGrT7!J}!0p7Y8{08FJ;GcJ8 zr2WwGi#F)BYd`7jC3S;n&h=4!yP3SsYOX#;8w?1{*p!gOi2+HDHG*;CXjA_~6skYC z3jdYl!C6qP5X~|si>uY!94n$`Vq-5d>h<3Rt0n_A()1NK;##Kt05gA)<~3_lIe3M= zVr62OUPm7yji*R*9b(P*`h({x3Hy_bqn{ds_4-dvTFVE;$JdC~68&m$1ttwPK1QM; zVyI%wi?I&lUr)X$+HLPgMLXl`kA7nBM_c-!M-`;kKmTcOFZAJsK6c{LGVWo+00TBbd?2G^0np^ug_&;x}&kW==_#Qqpu1SVn?h=rc$ zWuwqm0YMSbel4);DH!ZtjC_*T; z3dNEG4~iOWrqZ@&cvo?}B54%|tp{$fcc=jC8&_U-r&VHXaZN;Lcd|G^Aj_tRFc?VX zXqqwv<-)?(=nwq6P>{u6WIrCMNhIf*m}+DnC!;kHDh6YW$UPecE2$q&5ba`a(0bnZ zII5uRgHVxLcCz9-MTrMhiITB3PtNdi+GLXs5m#u=pTOFM{VRRMTb= z23d4PsDdN@rKg;Mxw~f~nZ~xKH#;O)cgCu>g-Fn{C7xWZrO^y8S4te8hKMVSP5*Md zUwcvf8;fK!MI;l8Fc2*QxphRSf;s-Br+gPvanBXt&((BxXU{VqMg|%=4y*OJ`H^u8 zZWYyza)J5j%~+VHdK&XHz4f7c&;N6^WzVhHJkv0!7`6M7MS0$<{27_|-()l{%@~vY zdVTMb-^JW-8jYR|tWk?IaI+@v4A^IIC)c~jlXO(gJGVMkv~WZa{EVwgGR7?H`12#Y zffl7jM=WtS9hpD^$RhW_ryKprKOOQP{37yu8w)AFx4v`9;h6i={~Pj8^90A-Z4?{k z%W5lOcw;_p(ry$Bdwuh6U*q%~sEZ>V);;*=SgzW7i1$!q9+KXi<7>RZx)0cE1U)Gl zZUfr72Ox3b00;hiQsve?E#>eR$b)qf@}#PA(u=SIUv#)S1uK}vcY#!N81X}k$4O42 zWDo8eapQU+N@m^_sc6CBqVcf7ZrKm=SO4-fnH;uE9FO}=>R4e*8FzpAKVF|X>?L1T zpM7TPTc3R-J&iuG0fC){ieu(W*p4`XT?&j1nC#Em4{^0oYRrw#n#4va~n<)eQpH%%-_)3v7$tO`amyh<3lSkzloZPt~gAj$n3a{>+onI znxuNB-L6F(l!KDOQA$xE3Lr{#%Im`#kC$vSEm($88p$ zN0DpL@ldh@Q=wKnaFtjOY91!V&noY2T#Jy9KNaYQ9IFP3=V7ap9o#Lh9Tb#ov|PmO zd}?tMB0>Tm79UJ`$N?6dLu&`EBBl^Od;VV_E*KcyadZv7M?#D{mp8I{%`Ab+yh6MevomM$&QyhGh8D!%##Ft+w(O)_Lpf6ef zc)R{~f8{Iuesn0vK5QSLj^~oo051R>sbCKnFuoj+>HL5&iL5&SVE|b}ReXSoXQS6J z{bxfY>{Rd%*ytnvUlqJp!Lfj&fgcQ*3HUCAK$S;64@mq11&aVl&rf68m#BD}ivJ0%#B{F!GTk;nrdy`svsC;-1qUho@3FtYeER{JelsBR->L9r zD(+VFe?fm_`kw)wr)gK;1$YtC--bR&*m#Dd*8-A$10eBF!<-1u20R~-&Nzf`EBH7d z`Q9>6;;&I~7$EWKXcVSDfzD0%D?sAU$3Hxh${G%MG2piW5k?LKMAKyT4A8VufUhg~ zBS4IWSt}Hz&m`sQMxzk!Q1IIda-W;{cNHuE{0`zSK>V^$hJcwU4^2@@2-T)#0D^T1 zA@I(lh``-(<;GH4p3C|Bp3Gt{gZ6_ehCj{OZ07Ta+C4^4XmI9Ja z2)uJDAiAvZCqP#e{sg2KO1#qrh)1lYgz&=^{sg3tka#EUAPqn3x$L_t{ueN1c#m-i zYCRWSi%a?wj6-D@bkONfARG+<=o?IW6Ab{NTG3so9K+>`?ndJ>lq!0+dJcS_MMq^A zHY)lim4B|HccbzQ4T?TX(VtZGkdps^RX!imBK%O%n-zch^iuwAmH%3#q5RY z92Y-3j^E(8cw=1pK3sz<;hfx;4{@J24BJ}6~dj($l=AY2G|42!7_AQ6qsRb<5<$8=6_YKp|w7Z1EoxBQHbG8;C3 z>qBPyvI~9Y@3T9b^sR$(^Yt~mbMy2yp)+;|P9_<(hkAS4e7)gGFV_44I{>L3Gjtnw z4Fb$Q%^m19T8eQ6#{LTs%Wux}{~cgnbH4xEfW`h){GZZ1#h=o`1{9lk9?P#NK_0W6 zeg}@^`s3}GL)+a+8OzSg?u_g+n8!wUK$>-v-qW`hC-<%#nq2cvdvXTlf>4!B3H*YK zPMP8QNF{h=177%{Ba^bMtrP0qb|==IzMm#oy!>`_IUovSSJ( zXYOD_iG3k=&Tm?9t?Ay5-uC)VLw|Vhqz_Pnr>?u^X9eXUZ+*yPwij@tIMQC{_t{=3 zyf~%TSXrE+)$|zcUM-{0{76)QOh8CLtp53W$jkh-;yC>z8d`w}bJ{w}7&gg%M z>3fZK&p?$P`*fl--}1d35GeXgV^5rHg~7*NkV@!=B90IuCTJl>SASBPOh^SqL0wT@ z$%MY(?IB;|v(zwPl1RKT2;L>(N?i4HUvEY*AsIxoNSzzuHTzoMEpZP|Cxz}=JD>RQ6nzHa+n=&=6 zl>@2B)i@DX7`^uR%rLYo*6Dy4e}y}M7{?`h^S6cO+fjzLU22Op^3mNN;x|@cvl_BD z)|}|v15PY%e^;S+boFrSEfXt680Q8Wi+zv#(e((Ern=Q_&FAie@_t+d-e^icdrnu_ zr(`^FA9R99jDOMexP6Z9gN`8yKWiV8%EixKuU8NkA%RaS(5q5v2ReA?^-yw*&wf~?B`v4k0UF(9V)k>>+qULH9V5Vy!9(*QA7XZ3;|=KZV>6?_|z_(Oo)m)Wo4 ztt$RQKNw)&jZB1M3z&*%M zC*TEus{om=RKY2L=OLaA$bBX5-;4s}JFVvfa$kw**$#w!_ca5sTg6{dut~vR!3>#h zpMqQmfyJ;o1@BQXOVQtju9^PF3bOk#J_B@Ap_CB$(~1ERN(g~>rURnoN(s>}X--8a z1m4+=#M~buL{p?421FWe9;r z`X=OI*r4dg6nzPBq`UBAC|2|(ivB2YBEO>FtLRTE`Y6(YtAnx_xW;5WpyFd0S9wVv zsPeZc`ea3)N+QBJiq7}M87lB&{%%GPu2FQKDxa)bp0&P*jzoN^%KtkguobWQQIn@pt3mZ^gwq4gwVk$K&F(OFR9_xcHH{^uLaaUl5l*e!oY& zKOFlial*^VF)uF1z8{QtJE|)9I(F9&l->B~Ug;@1Q zu5JIodP$e4*w=1j_X(_YxZStO6nUAqpBHjtQ(NyBMRz6Cen^ZweY{DO_{vis+bD?} zCZ~PxToA<$6`HA>@9(P(C<Yr|vm*1AnFiiJw zMDhjJyv~8^hNEZ&(o92;#_-x#8y;>hC#5mj2-hCLib*CJ)0ph5Z5wWc^Oh1(P`k%r zm@7yJs|Q|ItPJzIrH#o2wYweR$s(Pvw%uW_AesD@pRo8bnGeYzs3^C9yv+dVNNBpp z6x4=>8@*k6eKzhgjc|PcHp_y^W^xX%46{Cjg)|ZMYZTQmA6NoakVXB=AmZ(&mT+Et^?A7EHNyEWnONM=8sP@sG)Qq9*h(zFVIA@+oLdwQ+}ubN z8V;!{G(3qaWcM28iS0BDZ_t+%^;Zfy()fs1t38pT*Du6MH{$+S2P~xAU3)U6=HuSB z+K*Fe4q~^@SKI5Td8W6$_EbvE&)r*-u#qHycI@U_HhPm=Y|>f$>XQpAxLf$(Ry+EN2eN zzCO4irx^kv$AIgDbveg?D>7e0=eUxcL-dv+^FYz)-xZC1xp4IBvWpasJ|NvWmEPdG zniqCl!m0wEUs?HL2roZ!l;TT}sQ8N^8hv>3aH7K34VM^LdIJqbAeZw>;H8R`nV=I4 zhj|V6Ee|sR^e0S<4Po{|8TV~L6aEp_)9Y=ny{{ROFd{8l=!1R6cFmpKzDv`dY%|P; zbr=|posX`6mXRvlrM5h}{^_>X%F3U&87*4i-F;qjwUcRfj356D(ri0~J5J*{P5Z;X z=kQ^UcCYzFD(+da#5}og`}px2ndB*qYTz~%BiAm?r6JAgEm)%im2m-V&U~U7z}Wui z`hPTxohZYLGUpjyvmt=yH_X-$^x-wv7P2JB`H07BZp^^i-!M0JU@kSD#}pAV%(ca- z#I0U~@=8UW6eH&oQz=}yVGd@L@#BBiX6(>Hx)+FUNTF@IupO5! zRQ%EPueE{kV!WTyuz_vuI)|jgP{#Q2d#F^ffkv88Nod5wN}3H5z_n$6SCtzgwg}2v z2eoc)dr~A=`=N-n@{z+1sLX>r+$%tWRRGk|UZ9TQ8XT5Fjuw;_LKUW>awu*DRsI}y zm##tkYnOV^!tEg_d5KW@#t|q5)jou&5;_b~%jW6~q2IO507#h&<3KadFhEHay59IQ zDg{&c9vY~DlQy(;n7TySCc1}+P!R*Cfq-3ja((R+~1Y@G|0p~^YZ*qT}8{asJFAzhV}`zbI(fpw4@$@e0` zE-i2snsO~4pdf2E__prvdWU&hgn_O-;i9Q%XYYl%e>;Svp9ZKyKb|+#Gl$%CS+X<>I@ErBJ2o>ig zT3`ZOz-$OfFoL=QsTEEAvOujZ6v|xd5F)QWfn6$;ngjhHSD^%Z0iabrI-)PyoK zs5EHmD&h_c-0EUf5E|0jz+eq~nYpcsovB50kJyKHc#3u1zQ4=CL8kIrSnj2$^=9rJ z02~%gx7sNhXmuR`q@M%MXrc`$4P)07^z4nBkQ~gqDgATk+Nfm17+4ncSB0uTwoFtK zW%J!q=wgkqF=1{_BvY2R7tN+F)4^O_9;FIh0SLQxi(t~g9U8>e0ibFfyP=d%agAZ? zHRy2lIAJ^6W(cr$B^0o0f7b};q@hD7VD(Ta0K?sD*Mi9sF|stYKn>tEo0eRK)@`Y* zd>{mkOu^W&(JkuzL;xB`y^eq{>Xk7HLcri^kz#GPu%C?%8p(F-@#<9&$<3LYc1l;N^!hE(lpu6X^?!Bz%RH4p>CaQDi{;ri^x9u=-;f5~A zq6N%-7>Ql$_tC6S^UBJL*hhC#X{3c4iXo7~ZbA)({=yA8C|B$LA&1GmymU`QX{$$| zg;3~DXn%WUWtu^Cih*aF2Q{Wv(8R5)pgZVsRj5_lw95K0M-CtB72V2Tqn)y5KGzLbWv6tZ4zR}#4~#Wr?4y0U zQOE`;0{TCUlA)F(qO7$Ul6$7$*g6*+QHw2-|A#U*M%csZDwk;8fz&BI1)&SE2B0o5 z%%Nzdp;Td`BEeLepzOsR(aU0XXav+G;1Y@oajxo;_Td!vp@lLr!|(4JgL-ZA3tMa5 zB8;VBoa`Tqp#4_Oy;H=E&DlFUtl^POj}z-(W_5@Y3um$AB1531rTYik9@-on^F>V1d~H%>J&f=TURN#B2cu(->rKD8VK?KpiX?YH+a{8@@x- zw#Niz^$`(k{ZQKVUi7aWn)^TN>mGlK}W>8fkkg)+Tp|oehQu&XHblA1Ko6$No@f}&8Mq~#i+)z zWx6c`fst0&!b=W$m_SeoxD$}|PLZINtUm=7+6n;!g1Dnb!{TxsyR~e%L#Q`qO*gQ7 zvY|vMM%9G~=&hpfHsqcwI>a^9V>1#$j}yojlP!#ES7qf`D&37#08?RTJG5CFX@IkU zZ(S;syxK1W+RW!E=vS!#n8~|18t)Zz*yeO$s2e?lg~X3hAX)n(vMq8&Qq8x4eOhY%FtAGVONo_@x)WFo?scM*u3ia%rL2i*V{IX;W;U!6l(k31*wc{V zqP{l$hWgrsA!qYcsQ@2<7^=-}O_J3NIdrY#V1*V#O!bapU{vpTjA_CRM}SeiW0C0b z+nP}fG&ezxXMV&)UkNvS3XJ9!Hj~S8Fz^#2jou+KtD>E~3&Vl3nI?f1%03{lP&O9# zz@pV(1_o7LCt?kUWUK{)2o!LFqffY@L<}inir)YXR&VKi8eh_G+{W^*v`*ZWid$?0_vwBe_neR&HXF}$Q5xmRXE*ra$8MI_uW-yk-Pc)vV{H}V_A-jKRqt`qmmQ-Z4! z?w9er`+Dv!S@+9$)?PIFt%A{qibnsT!2D%_dC)idP+pi@CGLwb@;HW{S6|*%P`lTG z>q*3g9RrE-@gCV2E-=@LM7Zm1qGH0eyz7vhZI6}}r{%BvLr`Sr*B7bBhFzYG8K{8jj$@1KGH z#r`|-KcnUf@eED9j|Gf$tMdN`F6m44HXwbg-T*9aHvF#wW;AE`U+iqT7=>64S!#cx zvtXCS&?7x>r^E-L=Z10UzLirixYALok?6Zy{%d|XUN0nI-=@UhT-u2rtpen#Mg zVTc_+{{&p;89M}zOkC$1JA?=$a9wPOLB*erYbnG~aDu)Keb-Ack1$}Tulyi<^2$3o zYGL%k7p|Oe1l#pK;szd~ch~jG{Exsc9rriIN0WHgy4ce(exVD~Pf~nrKz!t`Kt8?66g?fA zw(RTIV&_v}SLuuWS-Bb8FNWdNOaP~$L$`BP3isKXOJ#o&194A1 zkRmR|25}|ZU*th@|EqzMDEnU_S>PlJ=d+N&Nr(U10H*!#pq2g4Vm#67RtY>xe@0vs z`_IKFVu+cb<^n^^B!B)pJAX0sN5_}k6nskLz_LtMC;%3@)(CNSZocu0WWSzUZ2SUW zP}u?+^^M~7hkk=05+dkz+>HPPYCOhez8a5lS**rmTnhWN+c$<2W#1EVsqCBcxw7x`aH;H@lcTclsJ?0GK8~e* zSUOw^3-}#Htoo$Ao3AO0_n5C@_?c1oFuw69K2l|U|9FPCOPVr=sb(0nuyTFH#Luw$ z89VqEm?qXai_gS2SGx3ue}w`tAZ6Otu41gP*$1rQ1-*|Zt1seXDls?qEWxFvPob)w z$tX_QuUKL9JkIhte+aGiToj=Th4-8ti5p^l(oI9k8LX&+Lq-fW-rtC}}Cc0DQz#aE*76nYAlyJFTdZkONX zL|ORO*fuhnex&|&_~+`GhnG!4BvPh$Q!n$dRcs&Sy(!8Qe(Lp5`Af_2*~cD!5+{A+ zIVfDJE0RAcLEe{Ow7qSUKBBxcaLY5}Wy<|Fgc2(|$bCz^+^`JE76oE4eg*znWe$5B z89|i&Ro|Le-fmlYC4;{w%W0#2CdZcJGf(<{eYEz?B;OeaR3pvk?71E0iT!VrZV1=n*NU<2T0U>(FY~SiFOsmoYp0t(;z$R|4Tzr>2))pI?X7dJ`u&Z;-w=<` zeeLLy)w};QbBOfoU9Y>}oBPL=&z?AP&Gx3Bj(%k8AFsi(NQC(BX+uGp13$|T_jlad zVTk&xVwFKZUkuEAkDt2wb2vI&5Bl`Cj-l>w?%H`HbA9ECZ~F0(^sJ@$t~fV&`c2O*4-umBNWCd^y!=-wye70L z_)Gotliq3;O?;?8q-T*4c_tuY*N>mg#SoQ;EvMRXRkCb6N5c2i$DkvM5Y2~eg=jhP z_A-HmGDOp>b|%5@C(ZQH4q~J3OmMb_A9hBv~i-u;o<=D1m{6+=0qOK4%>t_mb zVaIar1BA#~a{%!G<7x~R0zXaRCj*ktI{ZTsrGzLVZ4Dqo2_f*#sesUDDIr=atx(Yk zfp@w9(P2snVR&gH6rB)w@!mUj!o+*;(4Tnko&1l$o`839zQqPtDIq#)n((*)oe+5E zNEXlZi* z5lRSwce()4l%<4t&=BL74ZM?QM$s=jzyt_C!B6;2Ls*7ipdR$gJ}^qg{|$ijW<}qp z=yaZENM;evo~|1y^Qz2HW3xSLr~Boalr~;l=V8f_g38#x(AfsXotoih3anXJ0FK zj56$$7(LhtQB~NnQEk|8O-e4zjUqIIpPdKhxR~&He$CB_ z-!a7`1b%`C1Re)Yn3Fm6L;{m>Xbxaek^vaPq9Oz&7mZW#V_+I99p*UBQ(rppzSer~ z@!xcg1Z4{U?>Q3x|DGd3m-t`7CEh$PHf*J9{L6X4;=7aaMbL~S+Lylh0Ke>7SdMqG z<$KQ}Elz?+JUc?KOTko&Jrlj|Q`k;+FI?U0@jCAJeKxxU zXYFf?m;d?Zzjdy~CX^I8*J2rl4StFD3!Q5*6FI-&xfTT&!mJG+wOz&I7d+6a^#yV6 zcLmrS-~722H=!8Ycfjj5Auf*=cm!8{pKH-Y(Q0|F#Q@wUK}hRdiw{T;AKCsU&$YN3 z_ey*{q~WnMa3M})z_}JzAOgLpb1g)ZGKvp0%R?<9A8oeIvv@2b+&6fh#Xo|nRM*-l z`#w)%z;o%nNMe2IyzV2J-gchFTM-W5d#W|LvwM%{#SL5#XVRB z^>K~``y0Ot&z5Q6xiI33?F}<=#rEN`7tv42%g@KJBTk-Q!6Z1gVu*JGaS3NvJcCcn zGl@LA;w=m`_9m=5&4_|MVrKj*XIK0izia>ClladMsYAH0dv?Wqu!`9i$2kQ&$D#~) zIdpM5`iq`pkw@{Y#{zYn6Yy@?SS76(Ml2m!>cN-JW6rTSD_$PkIR!t!`Vm#`^BjxM zAERnZz}hOF`fstxAkIp_nL3KjuMQV|-l-0=DJswX;5iqih`$`E{F?E#HuB7E=V~JtH7G)bmxM6)pIb~zxi`8YQ&tAa1I8$ zQ=~%%M1PDur{U|KgYg-r0x_yS@bho(9E`Gl&cXP}H+K$3C7wr#_Wp*>!T3+iDLtRD zU1RmD&Z$u6U8r6+!~5mV!TA1v-#HjpQ|>|0b1-WAJO|@{{J8<*R_Du}8_;K_zRwNl zBdOzm&tVYbP4rxgJ`PvF{V2v*>oAPyxfXpS{VUXGe5?PpoLlf!&yBzd6mjQLNQYCo zK$QP^KmS$vAA)nj;9tQizT~+SYS=7ES&O`JewL?`N93~~!Ew4nIUw9jo?0F~ltWhj zbf-Q*9QsdKXFky0)i&zqe(x-Ui5VPxYGZpoFXD;cJHs}s4AJy>tQUQ5o{NYBKg-YJ z9qVCGG(ER-@qC`CKL6Z@ntti0!VC#=VCB%qxexLA z|6iZ0@MX_6h?l=uRKBnJrvmxSpUc4KfOb1jk5?1*ujGoCcecVOwpZWs&W6FGysR4m z(V|&b1LApm)?h#w+^`fC|1|_*d@mr={{oQd=cxEiD!va^O+SR7f-@ET2Rzg!{%Qr! zP>=_+68|nXvIw^;SO&=QiU8^Fa4Db*@XxSQ#$N+U3cpa{zoYQysW_cInC}vFlu>|H=ns^WZ!{5J1xP;Q0Y|{U;uGBDGF>hp z>6r?i1;}*CfK2yG+}ttUg$h21=3)A~6}&^i3jvw$RVbVE-3tB~ka&Jmfp~r&jj#-m z_^0s?{vKIR0A2w2Q$U1~{|69FlEv@XJe`s;#42_f*#m4IlU^NJy{4Kj6jr z3b+>&=PN+P;(P_tJ@^N_I9~xO6z40TNO8Ub=??q@-Wh@(@QAUL5YN}d`3j`(lz3+V z5M8B|5bYqo-$(jVi5K7RgG+((qU7?JD;Xd z_*u^@pA|ylxhI%0{1`v_J%pq^Xn(;q{T3`gku=8fI{@V|j6{CY$EoL?p8`jE0FA}K zb|BsI)A+v1?^5NhSM*uxx$YE2Z-NpSDy{sAe~zM$Q}iN5=krm9If`yg)b}g;IK_XI zRX&u3Y$pLrk>rh?_|aN;#NF_I76!J0?Dn?IV8pEMFw0?F;>VVj~~?Lpe4U zc65}5#UZXagl~r|(bg7mSgCm%GET)dN~}h0gs3hOkWtr@NY!oCk%tXLDA5}AOOWRl VL~_`8#ZjWXC{tlxw#}VI`!9O-I #include +#include #include #include"nnparams.h" -#ifdef _pack -#pragma pack(1) -#endif - - float forbidden_entropy; -/////////////////////////////////////////////////////////////////////////////// -// Construction/Destruction -/////////////////////////////////////////////////////////////////////////////// + +float nparam_GetInitialEntropy(PNNParams nparm) +{ + return -5.9f+nparm->rlogc; +} -/////////////////////////////////////////////////////////////////////////////// -// FUNCTION: void CNNParams::InitParams() -// -// PURPOSE: Initialize nearest neighbor parameters. For now, simply set -// params as requiered. Extend to read from file sometime. -// -// PARAMETERS: -// none -// -// RETURN VALUE: -// void -// -// REVISION HISTORY -// $ 00sep06 : created LK -// 00dec17 : modified to use new parameters LK -// 00dec29 : included dangling end data LK -// #$ +//Retrieve Enthalpy for given NN-Pair from parameter table +float nparam_GetEnthalpy(PNNParams nparm, char x0, char x1, char y0, char y1) +{ + return ndH(x0,x1,y0,y1); //xx, yx are already numbers +} + + +//Retrieve Entropy for given NN-Pair from parameter table +float nparam_GetEntropy(PNNParams nparm, char x0, char x1, char y0, char y1) +{ + //xx and yx are already numbers + char nx0=x0;//nparam_convertNum(x0); + char nx1=x1;//nparam_convertNum(x1); + char ny0=y0;//nparam_convertNum(y0); + char ny1=y1;//nparam_convertNum(y1); + float answer = ndS(nx0,nx1,ny0,ny1); + /*Salt correction Santalucia*/ + if (nparm->saltMethod == SALT_METHOD_SANTALUCIA) { + if(nx0!=5 && 1<= nx1 && nx1<=4) { + answer += 0.5*nparm->kfac; + } + if(ny1!=5 && 1<= ny0 && ny0<=4) { + answer += 0.5*nparm->kfac; + } + } + /*Salt correction Owczarzy*/ + if (nparm->saltMethod == SALT_METHOD_OWCZARZY) { + float logk = log(nparm->kplus); + answer += ndH(nx0,nx1,ny0,ny1)*((4.29 * nparm->gcContent-3.95)*0.00001*logk+ 0.0000094*logk*logk); + } + return answer; +} + +/* PURPOSE: Return melting temperature TM for given entropy and enthalpy +* Assuming a one-state transition and using the formula +* TM = dH / (dS + R ln(Ct/4)) +* entropy = dS + R ln Ct/4 (must already be included!) +* enthaklpy = dH +* where +* dH = enthalpy +* dS = entropy +* R = Boltzmann factor +* Ct = Strand Concentration +* +* PARAMETERS: +* entrypy and enthalpy +* +* RETURN VALUE: +* temperature +*/ + +float nparam_CalcTM(float entropy,float enthalpy) +{ + float tm = 0; // absolute zero - return if model fails! + if (enthalpy>=forbidden_enthalpy) //||(entropy==-cfact)) + return 0; + if (entropy<0) // avoid division by zero and model errors! + { + tm = enthalpy/entropy;// - kfac; //LKFEB + if (tm<0) + return 0; + } + return tm; +} + +/* PURPOSE: Initialize nearest neighbor parameters. +* +* PARAMETERS: +* none +* +* RETURN VALUE: +* void +*/ void nparam_InitParams(PNNParams nparm, float c1, float c2, float kp, int sm) { @@ -433,23 +466,7 @@ void nparam_InitParams(PNNParams nparm, float c1, float c2, float kp, int sm) return; } -void nparam_UpdateParams(PNNParams nparm, char * s1, char * s2) -{ - float l1 = strlen(s1); - float l2 = strlen(s2); - if(l1gcContent = nparam_CountGCContent(s1); - nparm->gcContent /= (l1-2); - } else if (l1>l2) { - nparm->gcContent = nparam_CountGCContent(s2); - nparm->gcContent /= (l2-2); - } else { - nparm->gcContent = nparam_CountGCContent(s1)+nparam_CountGCContent(s2); - nparm->gcContent /= (l1+l2-4); - } -} - -float nparam_CountGCContent(char * seq ) { +int nparam_CountGCContent(char * seq ) { int lseq = strlen(seq); int k; float count = 0; @@ -461,352 +478,89 @@ float nparam_CountGCContent(char * seq ) { return count; } -/////////////////////////////////////////////////////////////////////////////// -// FUNCTION: static char CNNParams::getComplement(char mychar) -// -// PURPOSE: return Watson-Crick complement to mychar -// -// PARAMETERS: -// mychar - Character to complement -// -// RETURN VALUE: -// char - complement to mychar -// -// REVISION HISTORY -// $ 00sep28 : created LK -// #$ - -char nparam_getComplement(char mychar, int asnum) +void nparam_CleanSeq (char* inseq, char* outseq, int len) { - if (mychar==1) // A -> T - return 4; - if (mychar==2) // C -> G - return 3; - if (mychar==3) // G -> T - return 2; - if (mychar==4) // T -> A - return 1; - if (mychar==5) // $ -> $ - return 5; - if (!asnum) + int seqlen = strlen (inseq); + int i, j; + + if (len != 0) + seqlen = len; + + for (i = 0, j = 0; i < seqlen; i++) { - if (mychar=='A') // A -> T - return 'T'; - if (mychar=='C') // C -> G - return 'G'; - if (mychar=='G') // G -> T - return 'C'; - if (mychar=='T') // T -> A - return 'A'; - if (mychar=='$') // $ -> $ - return '$'; + switch (inseq[i]) + { + case 'a': + case '\0': + case 'A': + outseq[j++] = 'A'; break; + case 'c': + case '\1': + case 'C': + outseq[j++] = 'C'; break; + case 'g': + case '\2': + case 'G': + outseq[j++] = 'G'; break; + case 't': + case '\3': + case 'T': + outseq[j++] = 'T'; break; + } } - else - { - if (mychar=='A') // A -> T - return 4; - if (mychar=='C') // C -> G - return 3; - if (mychar=='G') // G -> T - return 2; - if (mychar=='T') // T -> A - return 1; - if (mychar=='$') // $ -> $ - return 5; - } - return '*'; + outseq[j] = '\0'; } - - -/////////////////////////////////////////////////////////////////////////////// -// FUNCTION: static bool CNNParams::isMismatch(char,char) -// -// PURPOSE: Return true if char1 and char2 are not watson-crick pair -// -// PARAMETERS: -// char1 - first character -// char2 - second character -// -// RETURN VALUE: -// bool - true if char1,char2 are Watson-Crick, false otw -// -// REVISION HISTORY -// $ 00sep28 : created LK -// #$ - -int nparam_isMismatch(char c1,char c2) +//Calculate TM for given sequence against its complement +float nparam_CalcSelfTM(PNNParams nparm, char* seq, int len) { - if (c1+c2==5) - return 0; - else - return 1; + float thedH = 0; + //float thedS = nparam_GetInitialEntropy(nparm); + float thedS = -5.9f+nparm->rlogc; + float mtemp; + char c1; + char c2; + char c3; + char c4; + unsigned int i; + char nseq[50]; + char *useq = seq; + + nparam_CleanSeq (seq, nseq, len); + useq = nseq; + + int slen = strlen(useq); + + //fprintf (stderr,"Primer : %s\n",useq); + for ( i=1;idH[c3][c4][c1][c2];//nparam_GetEnthalpy(nparm, c3,c4,c1,c2); + thedS += nparam_GetEntropy(nparm, c3,c4,c1,c2); + } + //fprintf(stderr,"------------------\n"); + mtemp = nparam_CalcTM(thedS,thedH); + //if (mtemp == 0) + //{ + // fprintf(stderr,"Enthalpy: %f, entropy: %f, seq: %s\n", thedH, thedS, useq); + //exit (0); + //} + return mtemp; } - - char nparam_convertNum(char c) - { - if (c=='A') - return 1; - if (c=='C') - return 2; - if (c=='G') - return 3; - if (c=='T') - return 4; - if (c=='$') - return 5; - return 0; - } - - /////////////////////////////////////////////////////////////////////////////// - // inline function - /////////////////////////////////////////////////////////////////////////////// - // FUNCTION: float CNNParams::GetEntropy(char,char,char,char) - // - // PURPOSE: Retrieve Entropy for given NN-Pair from parameter table - // - // PARAMETERS: x0,x1,y0,y1: Pairs to look up in form - // 5'-x0-x1-3' / 3'-y0-y1-5' - // - // RETURN VALUE: - // float: Entropy dS for given NN pair - // - // REVISION HISTORY - // $ 00sep06 : created LK - // 00dec29 : included dangling end parameters - // 01feb07 : rewritten. looks ugly now, but is FAST! inline. - // #$ - - float nparam_GetEntropy(PNNParams nparm, char x0, char x1, char y0, char y1) - { - char nx0=nparam_convertNum(x0); - char nx1=nparam_convertNum(x1); - char ny0=nparam_convertNum(y0); - char ny1=nparam_convertNum(y1); - float answer = ndS(nx0,nx1,ny0,ny1); - //Salt correction Santalucia - if (nparm->saltMethod == SALT_METHOD_SANTALUCIA) { - if(nx0!=5 && 1<= nx1 && nx1<=4) { - answer += 0.5*nparm->kfac; - } - if(ny1!=5 && 1<= ny0 && ny0<=4) { - answer += 0.5*nparm->kfac; - } - } - //Salt correction Owczarzy - if (nparm->saltMethod == SALT_METHOD_OWCZARZY) { - float logk = log(nparm->kplus); - answer += ndH(nx0,nx1,ny0,ny1)*((4.29 * nparm->gcContent-3.95)*0.00001*logk+ 0.0000094*logk*logk); - } - return answer; - } +float calculateMeltingTemperatureBasic (char * seq) { + int gccount; + float temp; + int seqlen; - - /////////////////////////////////////////////////////////////////////////////// - // inline function - /////////////////////////////////////////////////////////////////////////////// - // FUNCTION: float CNNParams::GetEnthalpy(char,char,char,char) - // - // PURPOSE: Retrieve Enthalpy for given NN-Pair from parameter table - // - // PARAMETERS: x0,x1,y0,y1: Pairs to look up in form - // 5'-x0-x1-3' / 3'-y0-y1-5' - // - // RETURN VALUE: - // float: Enthalpy dH for given NN pair - // - // REVISION HISTORY - // $ 00sep06 : created LK - // $ 00dec29 : included dangling end parameters - // $ 01feb07 : rewritten. looks ugly now, but is FAST! inline. - // #$ - - float nparam_GetEnthalpy(PNNParams nparm, char x0, char x1, char y0, char y1) - { - char nx0=nparam_convertNum(x0); - char nx1=nparam_convertNum(x1); - char ny0=nparam_convertNum(y0); - char ny1=nparam_convertNum(y1); - return ndH(nx0,nx1,ny0,ny1); - } - - - - /////////////////////////////////////////////////////////////////////////////// - // inline function - /////////////////////////////////////////////////////////////////////////////// - // FUNCTION: float CNNParams::CalcTM(float entropy,float enthalpy) - // - // PURPOSE: Return melting temperature TM for given entropy and enthalpy - // Assuming a one-state transition and using the formula - // TM = dH / (dS + R ln(Ct/4)) - // entropy = dS + R ln Ct/4 (must already be included!) - // enthaklpy = dH - // where - // dH = enthalpy - // dS = entropy - // R = Boltzmann factor - // Ct = Strand Concentration - // - // PARAMETERS: - // - // RETURN VALUE: - // - // REVISION HISTORY - // $ 00sep06 : created LK - // $ 01jan07 : modified and corrected - // $ 01feb07 : optimized!!! inline - // $ 01feb09 : changed to include r ln ct in entropy!!! - // #$ - - float nparam_CalcTM(float entropy,float enthalpy) - { - float tm = 0; // absolute zero - return if model fails! - if (enthalpy>=forbidden_enthalpy) //||(entropy==-cfact)) - return 0; - if (entropy<0) // avoid division by zero and model errors! - { - tm = enthalpy/entropy;// - kfac; //LKFEB - if (tm<0) - return 0; - } - return tm; - } - - - - /////////////////////////////////////////////////////////////////////////////// - // inline function - /////////////////////////////////////////////////////////////////////////////// - // FUNCTION: void CNNParams::AlterTM(float) - // - // PURPOSE: TM can be altered by a new assignment - // - // PARAMETERS: tm_new: new value for TM - // - // RETURN VALUE: - // - // REVISION HISTORY - // - // #$ - - void nparam_AlterTM(PNNParams nparm, float tm_new) - { - nparm->new_TM = tm_new; - return; - } - - - - /////////////////////////////////////////////////////////////////////////////// - // inline function - /////////////////////////////////////////////////////////////////////////////// - // FUNCTION: float CNNParams::CalcG(float entropy, float enthalpy) - // - // PURPOSE: return free energy G for given entropy and enthalpy - // Assuming a one-state transition and using the formula - // G = dH - new_TM * dS - // dH = enthalpy - // dS = entropy - // new_TM = value for the optimal melting temperature of - // the last iteration - // - // PARAMETERS: entropy and enthalpy - // - // RETURN VALUE: free Energy value (dG) for the - // - // REVISION HISTORY - // - // #$ - - float nparam_CalcG(PNNParams nparm, float entropy, float enthalpy) - { - float freeEnergy = -999999999; - if (enthalpy>=forbidden_enthalpy) - return -999999999; - if (entropy<0) // avoid division by zero and model errors! - { - entropy = entropy * -1; - enthalpy = enthalpy * -1; - freeEnergy = enthalpy - nparm->new_TM * entropy; - } - return freeEnergy; - } - - - - - /////////////////////////////////////////////////////////////////////////////// - // inline function - /////////////////////////////////////////////////////////////////////////////// - // FUNCTION: float CNNParams::CalcSelfTM(char*) - // - // PURPOSE: Calculate TM for given sequence against its complement - // - // PARAMETERS: - // char* - Sequence to consider - // - // RETURN VALUE: - // float - Melting temperature in degrees Kelvin - // - // REVISION HISTORY - // $ 01jan09 LK : created - // $ 01feb07 LK : inline. - // #$ - - float nparam_CalcSelfTM(PNNParams nparm, char* seq) - { - float thedH = 0; - float thedS = nparam_GetInitialEntropy(nparm); - char c1; - char c2; - char c3; - char c4; - unsigned int i; - for ( i=1;i5) - { - if (c3=='A') - c3=1; - if (c3=='C') - c3=2; - if (c3=='G') - c3=3; - if (c3=='T') - c3=4; - if (c3=='$') - c3=5; - } - if (c4>5) - { - if (c4=='A') - c4=1; - if (c4=='C') - c4=2; - if (c4=='G') - c4=3; - if (c4=='T') - c4=4; - if (c4=='$') - c4=5; - } - - thedH += nparam_GetEnthalpy(nparm, c3,c4,c1,c2); - thedS += nparam_GetEntropy(nparm, c3,c4,c1,c2); - } - - return nparam_CalcTM(thedS,thedH); - } - - float nparam_GetInitialEntropy(PNNParams nparm) - { - return -5.9f+nparm->rlogc; - } - - + seqlen = strlen (seq); + gccount = nparam_CountGCContent (seq); + temp = 64.9 + 41*(gccount - 16.4)/seqlen; + return temp; +} diff --git a/src/libthermo/nnparams.h b/src/libthermo/nnparams.h index 01948c7..5c0cb36 100644 --- a/src/libthermo/nnparams.h +++ b/src/libthermo/nnparams.h @@ -2,72 +2,49 @@ * nnparams.h * PHunterLib * - * Created by Tiayyba Riaz on 7/2/09. - * Copyright 2009 __MyCompanyName__. All rights reserved. + * Nearest Neighbor Model Parameters + * + * Created by Tiayyba Riaz on 02/07/09. * */ -//============================================================================= -// Module: nnparams.h -// Project: Diploma Thesis - Probe Selection for DNA Microarrays -// Type: header file - Nearest Neighbor Parameters / Model. -// Language: c++ -// Compiler: microsoft visual c++ 6.0, unix/linux gcc -// System/OS: Windows 32, Sun solaris, Linux, other unix systems (untested) -// Database: none -// Description: class CNNParams - Nearest Neighbor Model Parameters -// Author: kaderali -// Date: 9-12/2000 -// Copyright: (c) L. Kaderali, 9/2000 - 12/2000 -// -// Revision History -// $ 00sep07 LK : created -// 00dec29 LK : changed to include dangling end data -// 01jan09 LK : included CalcSelfTM function -// 01feb07 LK : optimized -// #$ -//============================================================================= - -#if !defined(AFX_NNPARAMS_H__05604705_84E8_11D4_A001_000000000000__INCLUDED_) -#define AFX_NNPARAMS_H__05604705_84E8_11D4_A001_000000000000__INCLUDED_ - -#if _MSC_VER > 1000 -#pragma once -#endif // _MSC_VER > 1000 +#ifndef NNPARAMS_H_ +#define NNPARAMS_H_ #include #include - -#ifdef _pack -#pragma pack(1) -#endif - +//#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 forbidden_enthalpy 1000000000000000000.0f -//#define forbidden_enthalpy_div1000 1000000000000000.0f -// forbidden entropy=-rlogc -// #define forbidden_entropy 30.205986374220235304486574573422f -// Boltzmann factor (cal/degrees C*mol) #define R 1.987f #define SALT_METHOD_SANTALUCIA 1 #define SALT_METHOD_OWCZARZY 2 -// Strand concentration (assumption!) (M) -// #define Ct 0.000001f -// r*ln(ct/4) as required by many formulas -//#define rlogc -30.205986374220235304486574573422f + +#define DEF_CONC_PRIMERS 0.0000008 +#define DEF_CONC_SEQUENCES 0 +#define DEF_SALT 0.05 + +#define GETNUMCODE(a) bpencoder[a - 'A'] +#define GETREVCODE(a) 5-bpencoder[a - 'A'] + + extern float forbidden_entropy; -//----------------------------------------------------------------------------- -// class CNNParams -//typedef class CNNParams* PNNParams; +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 + -//#pragma GCC visibility push(hidden) typedef struct CNNParams_st { -//public: float Ct1; float Ct2; float rlogc; @@ -75,30 +52,19 @@ typedef struct CNNParams_st float kfac; int saltMethod; float gcContent; - float new_TM; // ge‰ndert von ML!!! - //CNNParams(); - //virtual ~CNNParams(); -//private: + float new_TM; float dH[6][6][6][6]; // A-C-G-T + gap + initiation (dangling end, $ sign) float dS[6][6][6][6]; }CNNParams, * PNNParams; - //void nparam_InitParams(PNNParams nparm, float c1=0.000001f, float c2=0.000001f, float kp=1, int sm = SALT_METHOD_SANTALUCIA ); - void nparam_InitParams(PNNParams nparm, float c1, float c2, float kp, int sm); - void nparam_UpdateParams(PNNParams nparm, char * s1, char * s2); - float nparam_CountGCContent(char * seq ); - char nparam_convertNum(char c); - float nparam_GetEntropy(PNNParams nparm, char x0, char x1, char y0, char y1); - float nparam_GetEnthalpy(PNNParams nparm, char x0, char x1, char y0, char y1); - float nparam_CalcTM(float entropy,float enthalpy); - void nparam_AlterTM(PNNParams nparm, float tm_new); - float nparam_CalcG(PNNParams nparm, float entropy, float enthalpy); - float nparam_CalcSelfTM(PNNParams nparm, char* seq); - float nparam_GetInitialEntropy(PNNParams nparm) ; - char nparam_getComplement(char mychar, int asnum); - int nparam_isMismatch(char,char); +void nparam_InitParams(PNNParams nparm, float c1, float c2, float kp, int sm); +int nparam_CountGCContent(char * seq ); +float nparam_GetEntropy(PNNParams nparm, char x0, char x1, char y0, char y1); +float nparam_GetEnthalpy(PNNParams nparm, char x0, char x1, char y0, char y1); +float nparam_CalcTM(float entropy,float enthalpy); +float nparam_CalcSelfTM(PNNParams nparm, char* seq, int len); +float nparam_GetInitialEntropy(PNNParams nparm) ; +float calculateMeltingTemperatureBasic (char * seq); +//void getThermoProperties (ppair_t* pairs, size_t count, poptions_t options); - - -//#pragma GCC visibility pop -#endif // !defined(AFX_NNPARAMS_H__05604705_84E8_11D4_A001_000000000000__INCLUDED_) +#endif