diff --git a/src/libthermo/dinkelbach.c b/src/libthermo/dinkelbach.c deleted file mode 100644 index 9bacffc..0000000 --- a/src/libthermo/dinkelbach.c +++ /dev/null @@ -1,76 +0,0 @@ -/* - * dinkelbach.cpp - * PHunterLib - * - * Created by Tiayyba Riaz on 6/7/09. - * Copyright 2009 __MyCompanyName__. All rights reserved. - * - */ - -//#include -//#include -//#include -#include "nnparams.h" -#include "galign.h" -#include "dinkelbach.h" - -void dnkl_initdinkelbach(pdinkelbach dnkl, PNNParams myDiParams, GGAlign mGAlign) -{ - dnkl->myDinkParams = myDiParams; - dnkl->myGAlign = mGAlign; -} - - -int dnkl_iteration(pdinkelbach dnkl, float TempK) -{ - int iteration = 0; - nparam_AlterTM(dnkl->myDinkParams, TempK); - float precursor_oldTM = TempK; - float newTM = TempK; - float oldTM = TempK, lambda, oldlambda; - lambda = 999999999999999999.9; - int positive=0; - - do - { - iteration++; - galn_InitBorder(dnkl->myGAlign); - galn_CalculateTable(dnkl->myGAlign, 1); - oldlambda=lambda; - lambda =galn_GetFreeEnergyK(dnkl->myGAlign, dnkl->myGAlign->maxloci,dnkl->myGAlign->maxlocj,newTM); - if(lambda < 0.0) - { - if(!positive) - { - //Initial melting temperature greater than expected - newTM = galn_GetEnthalpy(dnkl->myGAlign, dnkl->myGAlign->maxloci,dnkl->myGAlign->maxlocj)/galn_GetEntropy(dnkl->myGAlign, dnkl->myGAlign->maxloci,dnkl->myGAlign->maxlocj)-1; - lambda = oldlambda; - nparam_AlterTM(dnkl->myDinkParams, newTM); - } - else - { - precursor_oldTM = oldTM; - } - } - else - { - positive = 1; - oldTM = newTM; - newTM = galn_GetMeltingTempK(dnkl->myGAlign, dnkl->myGAlign->maxloci,dnkl->myGAlign->maxlocj); - nparam_AlterTM(dnkl->myDinkParams, newTM); - } - } - while(!positive || (newTM - oldTM > 0.001 && (lambda>0.0)&&(lambda 0.001 && ((lambda<0.0)||(lambda>oldlambda))) - { - iteration++; - nparam_AlterTM(dnkl->myDinkParams, precursor_oldTM); - galn_InitBorder(dnkl->myGAlign); - galn_CalculateTable(dnkl->myGAlign, 1); - lambda = galn_GetFreeEnergyK(dnkl->myGAlign, dnkl->myGAlign->maxloci,dnkl->myGAlign->maxlocj,precursor_oldTM); - newTM = galn_GetMeltingTempK(dnkl->myGAlign, dnkl->myGAlign->maxloci,dnkl->myGAlign->maxlocj); - nparam_AlterTM(dnkl->myDinkParams, newTM); - } - return iteration; -} diff --git a/src/libthermo/dinkelbach.h b/src/libthermo/dinkelbach.h deleted file mode 100644 index c90c8c1..0000000 --- a/src/libthermo/dinkelbach.h +++ /dev/null @@ -1,36 +0,0 @@ -/* - * dinkelbach.h - * PHunterLib - * - * Created by Tiayyba Riaz on 7/2/09. - * Copyright 2009 __MyCompanyName__. All rights reserved. - * - */ - -//Module: Achievement of the Dinkelbach algorithm - -#ifndef __dinkelbach_h -#define __dinkelbach_h - -//#include -//#include -//#include -#include "nnparams.h" -#include "galign.h" - -//using namespace std; - -//typedef class dinkelbach* DINKEL; -//#pragma GCC visibility push(hidden) -typedef struct dinkelbach_st -{ - PNNParams myDinkParams; - GGAlign myGAlign; - -}dinkelbach, * pdinkelbach; -//#pragma GCC visibility pop - -void dnkl_initdinkelbach(pdinkelbach dnkl, PNNParams, GGAlign); -int dnkl_iteration(pdinkelbach dnkl, float TempK); - -#endif /* __dinkelbach_h */ diff --git a/src/libthermo/galign.c b/src/libthermo/galign.c deleted file mode 100644 index 0f96688..0000000 --- a/src/libthermo/galign.c +++ /dev/null @@ -1,685 +0,0 @@ -/* - * galign.cpp - * PHunterLib - * - * Created by Tiayyba Riaz on 6/7/09. - * Copyright 2009 __MyCompanyName__. All rights reserved. - * - */ - -//============================================================================= -// Module: galign.cpp -// Project: Cubic Project - Calculation of melting temperature and free -// energy of two DNA strands -// Type: implementation - Thermodynamic Alignment. -// 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 GAlign - Thermodynamic Alignment Algorithm -// Author: leber -// Date: 01/2002 - 02/2002 -// Copyright: (c) L. Kaderali & M. Leber, 01/2002 - 02/2002 -// -// Revision History -// $ 00sep04 LK : created -// $ 00dec30 LK : changed to do local alignment of probe against -// one entire sequence -// $ 01feb07 LK : optimized! -// $ 01aug06 LK : corrected, included salt and concentration input; -// $ made true local alignment (problem with initial mismatches!) -// #$ -//============================================================================= - -#include -#include -#include -#include - -#include "galign.h" - - -#define nndH(a,b,c) (*(galn->dH+((b)*(galn->maxseq1len+1)*3+(a)*3+(c)))) -#define nndS(a,b,c) (*(galn->dS+((b)*(galn->maxseq1len+1)*3+(a)*3+(c)))) - - - -/////////////////////////////////////////////////////////////////////////////// -// Construction/Destruction -/////////////////////////////////////////////////////////////////////////////// - -//------------------------ Constructor ---------------------------------------- -// Initialize Thermalign Object. This will allocate memory for the dynamic -// programming matrices for both entropy and enthalpy of the present -// alignment. The max length of the two sequences to align has to be passed -// to the routine. Also, an initialized CNNParams Object must be passed to -// CThermAlign. Object will be initialized only once, sequences can be -// passed to it later. - -//int s1lmax,s2lmax; // global requiered for DEBUGGING - -void galn_initGAlign(GGAlign galn, int ALength, int BLength, PNNParams myNNParams) -{ - // Create dynamic programming matrices for entropy and enthalpy. - // The first two dimensions are for the two strings respectively, the - // third dimensions differentiates the alignment according to its end: - // type 0 Alignments will be those that align x(i+1) with y(j+1), - // type 1 Alignments align x(i+1) with a gap - // type 2 Alignments align y(i+1) with a gap - // All three of these types must separately be considered to select the - // optimum alignment in the next iteration. - // The alignment algorithm proceeds analogous to the well known - // Smith-Waterman algorithm for global alignments. - - galn->dH = (float *) calloc((ALength+1)*(BLength+1)*3, sizeof (float)); - galn->dS = (float *) calloc((ALength+1)*(BLength+1)*3, sizeof (float)); -// s1lmax=ALength; -// s2lmax=BLength; - - if ((galn->dH==NULL)||(galn->dS==NULL)) - { - printf("Fatal Error: Out of memory!\n"); - exit(-1); - } - - galn->GNNParams = myNNParams; - galn->seq1len = ALength; - galn->maxseq1len = galn->seq1len; - galn->seq2len = BLength; - - memset(galn->dH,0,(ALength+1)*(BLength+1)*3*sizeof(float)); - memset(galn->dS,0,(ALength+1)*(BLength+1)*3*sizeof(float)); - - // set dH / dS entries all to zero! - #ifdef _pack - // this is FAST, but works only if #pragma pack(1) has been set... - #else - // the following is slower, but works also when padding bytes have - // been inserted in the array for alignment by the compiler... - // time is not really critical here, as this initialization is done - // only once... -/* int a, b; - for ( a=0;a<=ALength;a++) - for ( b=0;b<=BLength;b++) - { - nndH(a,b,0)=0.0f; - nndH(a,b,1)=0.0f; - nndH(a,b,2)=0.0f; - nndS(a,b,0)=0.0f; - nndS(a,b,1)=0.0f; - nndS(a,b,2)=0.0f; - }*/ - #endif - galn_InitBorder(galn); // also need to do only once... -// forbidden_entropy = (-(GNNParams->R)*(float)log((GNNParams->Ct)/4.0f)); - return; -} - -//------------------------ Destructor ----------------------------------------- -// Free memory - -void galn_finiGAlign(GGAlign galn) -{ - free (galn->dH); - free (galn->dS); - return; -} - - - -/////////////////////////////////////////////////////////////////////////////// -// FUNCTION: void galn_InitStrings(char* s1, -// char* s2, -// int s1len, -// int s2len); -// -// PURPOSE: Initialize strings for alignment table! -// -// PARAMETERS: -// s1: Target sequence -// s2: Probe sequence -// s1len: Target sequence length -// s2len: Probe sequence length -// -// RETURN VALUE: -// void -// -// REVISION HISTORY -// $ 01jan03 LK : created -// $ 01feb07 LK : removed maximum local stuff -// #$ - -void galn_InitStrings(GGAlign galn, char* s1,char* s2,int s1len ,int s2len,int minxpos) -{ - galn->seq1=s1; galn->seq1len=s1len; - galn->seq2=s2; galn->seq2len=s2len; - nparam_UpdateParams(galn->GNNParams, s1,s2); -/* lk01feb07: removed maxloc-stuff - // if maximum local is not in reused subtable, need to reset! - if ((maxlocj>=minxpos)||(maxloci!=seq1len)) - { - maxloci=0; maxlocj=0; maxlocg=0; - } -*/ -// InitBorder(); //lk01feb07: removed, as it suffices if this is done once in constructor - return; -} - - -/////////////////////////////////////////////////////////////////////////////// -// FUNCTION: void galn_InitBorder(); -// -// PURPOSE: Initialize the lower and left border of the dynamic -// programming table. -// -// PARAMETERS: -// none -// -// RETURN VALUE: -// void -// -// REVISION HISTORY -// $ 00sep06 LK : created -// $ 00dec30 LK : complete revision - simply init table now, no -// more saving etc... -// $ 01feb07 LK : removed local alignment stuff and optimized -// #$ - -void galn_InitBorder(GGAlign galn) -{ - // Initialization for new alignment! - galn->maxlocg = -99999; // maximum melting temperature found (kelvin) - galn->maxloci = 0; - galn->maxlocj = 0; - galn->maxloct = 0; - int i, j; - - for ( i=0; i<=galn->seq1len; i++) - { - nndH(i,0,0) = 0.0f; nndH(i,0,1) = 0.0f; nndH(i,0,2) = 0.0f; - nndS(i,0,0) = nparam_GetInitialEntropy(galn->GNNParams); - nndS(i,0,1) = nparam_GetInitialEntropy(galn->GNNParams); - nndS(i,0,2) = nparam_GetInitialEntropy(galn->GNNParams); - } - for ( j=1; j<=galn->seq2len; j++) { - nndH(0,j,0) = 0.0f; nndH(0,j,1) = 0.0f; nndH(0,j,2) = 0.0f; - nndS(0,j,0) = nparam_GetInitialEntropy(galn->GNNParams); - nndS(0,j,1) = nparam_GetInitialEntropy(galn->GNNParams); - nndS(0,j,2) = nparam_GetInitialEntropy(galn->GNNParams); - } - - return; -} - - -/////////////////////////////////////////////////////////////////////////////// -// FUNCTION: galn_CalculateTable(); -// -// PURPOSE: Compute the entire dynamic programming table cell by cell. -// Requieres that InitBorder has been called. -// Computations will be done in the usual rowwise manner. -// -// PARAMETERS: -// int startRow: Row that computation shall commence at! -// -// RETURN VALUE: -// void - Data will be written directly to dH and dS in this object -// -// REVISION HISTORY -// $ 00sep06 : created LK -// 00dec31 : modified to begin at a specified row! -// #$ - -void galn_CalculateTable(GGAlign galn, int startRow) -{ -// #define forbidden_enthalpy 1000000000000000000.0f - register int i; - register int iminusone; - register int j; - register int jminusone; - // Determine previous character! - register char prevchari; - register char prevcharj; - register char seq1char; - register char seq2char; - register float entropy1; - register float entropy2; - register float entropy3; - register float enthalpy1; - register float enthalpy2; - register float enthalpy3; - register float g1; - register float g2; - register float g3; - float maxg; -// ofstream dboutstream("debug.out"); -// assert(startRow>0); - for (j=startRow; j<=galn->seq2len; j++) - { - if (j>1) - prevcharj = galn->seq2[j-2]; - else - prevcharj = 0; - jminusone = j - 1; - seq2char = galn->seq2[jminusone]; - for (i=1; i<=galn->seq1len; i++) - { - // local variables - if (i>1) - prevchari = galn->seq1[i-2]; - else - prevchari = 0; - iminusone = i - 1; - seq1char=galn->seq1[iminusone]; - - // --------------- Upper cell: Alignment of seq1[i] with seg2[j] - entropy1 = nndS(iminusone,jminusone,0); - enthalpy1 = nndH(iminusone,jminusone,0) ; - entropy2 = nndS(iminusone,jminusone,1); - enthalpy2 = nndH(iminusone,jminusone,1); - entropy3 = nndS(iminusone,jminusone,2); - enthalpy3 = nndH(iminusone,jminusone,2); - if (enthalpy1GNNParams, prevchari,seq1char, prevcharj, seq2char); - enthalpy1 +=nparam_GetEnthalpy(galn->GNNParams, prevchari, seq1char, prevcharj, seq2char); - - } - if (enthalpy2GNNParams, prevchari,seq1char, 0, seq2char); // 0 is gap! - entropy2 += nparam_GetEntropy(galn->GNNParams, prevchari,seq1char, 0, seq2char); - } - if(enthalpy3GNNParams, 0, seq1char,prevcharj, seq2char); - entropy3 +=nparam_GetEntropy(galn->GNNParams, 0, seq1char, prevcharj, seq2char); - - } - // choose optimum combination - g1 = nparam_CalcG(galn->GNNParams, entropy1, enthalpy1); - g2 = nparam_CalcG(galn->GNNParams, entropy2, enthalpy2); - g3 = nparam_CalcG(galn->GNNParams, entropy3, enthalpy3); - if ((g1>=g2)&&(g1>=g3)) - { - nndS(i,j,0) = entropy1; - nndH(i,j,0) = enthalpy1; - maxg=g1; - } - else if (g2>=g3) - { - nndS(i,j,0) = entropy2; - nndH(i,j,0) = enthalpy2; - maxg=g2; - } - else - { - nndS(i,j,0) = entropy3; - nndH(i,j,0) = enthalpy3; - maxg=g3; - } - // set to zero if alignment so far is mismatches only (tm will be zero -> no need to reset maxg)! - if(nndH(i,j,0)==0) - nndS(i,j,0)=nparam_GetInitialEntropy(galn->GNNParams); - // check if local maximum found! - if (((i==galn->seq1len)||(j==galn->seq2len))&&(maxg>=galn->maxlocg)) - { - galn->maxlocg=maxg; - galn->maxloci=i; - galn->maxlocj=j; - galn->maxloct=0; - } - // set to zero if alignment so far is mismatches only! - if (nndH(i,j,0)==0) - nndS(i,j,0)=nparam_GetInitialEntropy(galn->GNNParams); - - // --------------- Middle cell: Alignment of seq1[i] with '-' - // following changes lk 01jan08 - // we do not allow $- in the beginning, therefore set to forbidden if - // j=1 - if (j==1) - { - enthalpy1=forbidden_enthalpy; enthalpy2=forbidden_enthalpy; - entropy1=forbidden_entropy; entropy2=forbidden_entropy; - } - // also, disallow -$ in the end, therefore forbid if j=seq2len-1 - else if (j==galn->seq2len-1) - { - enthalpy1=forbidden_enthalpy; enthalpy2=forbidden_enthalpy; - entropy1=forbidden_entropy; entropy2=forbidden_entropy; - } - else - { - // end lk01jan08 - // -- entropy - entropy1 = nndS(iminusone,j,0); - entropy2 = nndS(iminusone,j,1); - enthalpy1 = nndH(iminusone,j,0); - enthalpy2 = nndH(iminusone,j,1); - if (enthalpy1GNNParams, prevchari, seq1char, seq2char, 0); - enthalpy1 +=nparam_GetEnthalpy(galn->GNNParams, prevchari,seq1char, seq2char, 0); - } - if (enthalpy2GNNParams, prevchari, seq1char, 0, 0); - enthalpy2 +=nparam_GetEnthalpy(galn->GNNParams, prevchari, seq1char, 0, 0); - } - } - // -- choose optimum combination - g1 = nparam_CalcG(galn->GNNParams, entropy1, enthalpy1); - g2 = nparam_CalcG(galn->GNNParams, entropy2, enthalpy2); - if (g1>=g2) - { - nndS(i,j,1)=entropy1; - nndH(i,j,1)=enthalpy1; - maxg=g1; - } - else - { - nndS(i,j,1)=entropy2; - nndH(i,j,1)=enthalpy2; - maxg=g2; - } - // check if local maximum found! - if (((i==galn->seq1len)||(j==galn->seq2len))&&(maxg>=galn->maxlocg)) - { - galn->maxlocg=maxg; - galn->maxloci=i; - galn->maxlocj=j; - galn->maxloct=1; - } - // set to zero if alignment so far is mismatches only! - if (nndH(i,j,1)==0) - nndS(i,j,1)=nparam_GetInitialEntropy(galn->GNNParams); - - // --------------- Lower cell: Alignment of '-' with seq2[j] - // following changes lk 01jan08 - // we do not allow $- in the beginning, therefore set to forbidden if - // i=1 - if (i==1) - { - enthalpy1=forbidden_enthalpy; enthalpy2=forbidden_enthalpy; - entropy1=forbidden_entropy; entropy2=forbidden_entropy; - } - // also, disallow -$ in the end, therefore forbid if i=seq1len-1 - else if (i==galn->seq1len-1) - { - enthalpy1=forbidden_enthalpy; enthalpy2=forbidden_enthalpy; - entropy1=forbidden_entropy; entropy2=forbidden_entropy; - } - else - { - // end lk01jan08 - entropy1 = nndS(i,jminusone,0); - entropy2 = nndS(i,jminusone,2); - enthalpy1 = nndH(i,jminusone,0); - enthalpy2 = nndH(i,jminusone,2); - if (enthalpy1GNNParams, seq1char, 0, prevcharj, seq2char); - enthalpy1 +=nparam_GetEnthalpy(galn->GNNParams, seq1char, 0, prevcharj, seq2char); - } - if(enthalpy2GNNParams, 0,0, prevcharj, seq2char); - enthalpy2 +=nparam_GetEnthalpy(galn->GNNParams, 0, 0, prevcharj, seq2char); - } - } - // -- choose optimum combination - g1 = nparam_CalcG(galn->GNNParams, entropy1, enthalpy1); - g2 = nparam_CalcG(galn->GNNParams, entropy2, enthalpy2); - if (g1>=g2) - { - nndS(i,j,2)=entropy1; - nndH(i,j,2)=enthalpy1; - maxg=g1; - } - else - { - nndS(i,j,2)=entropy2; - nndH(i,j,2)=enthalpy2; - maxg=g2; - } - // check if local maximum found! - if (((i==galn->seq1len)||(j==galn->seq2len))&&(maxg>=galn->maxlocg)) - { - galn->maxlocg=maxg; - galn->maxloci=i; - galn->maxlocj=j; - galn->maxloct=2; - } - // set to zero if alignment so far is mismatches only! - if (nndH(i,j,2)==0) - nndS(i,j,2)=nparam_GetInitialEntropy(galn->GNNParams); - - - /* - dboutstream<<"After cell ("<seq1len;i++) - { - for( j=0;j<=galn->seq2len;j++) - { - printf ("%f ", nndH(i,j,level)); - } - printf ("\n"); - } - printf ("\n"); -} - -void galn_printEntropyTable(GGAlign galn, int level) -{ - int i, j; - for( i=0;i<=galn->seq1len;i++) - { - for( j=0;j<=galn->seq2len;j++) - { - printf ("%f ", nndH(i,j,level)); - } - printf ("\n"); - } - printf ("\n"); -} - - -/////////////////////////////////////////////////////////////////////////////// -// FUNCTION: float GetEntropy(int i, int j); -// -// PURPOSE: Return Entropy for given cell -// -// PARAMETERS: -// i, j - Cell coordinates -// -// RETURN VALUE: -// float - Entropy (dS) -// -// REVISION HISTORY -// $ 00sep23 : created LK -// #$ - -/* inline */ float galn_GetEntropy(GGAlign galn, int i, int j) -{ - float entropy; - - // first, need to determine optimum values for dG and dH - float tm0, tm1, tm2; - tm0 = nparam_CalcTM(nndS(i,j,0),nndH(i,j,0)); - tm1 = nparam_CalcTM(nndS(i,j,1),nndH(i,j,1)); - tm2 = nparam_CalcTM(nndS(i,j,2),nndH(i,j,2)); - if ((tm0>=tm1)&&(tm0>=tm2)) - entropy = nndS(i,j,0); - else if ((tm1>=tm0)&&(tm1>=tm2)) - entropy = nndS(i,j,1); - else - entropy = nndS(i,j,2); - - return entropy; -} - - -/////////////////////////////////////////////////////////////////////////////// -// FUNCTION: float GetEnthalpy(int i, int j); -// -// PURPOSE: Return Enthalpy for given cell -// -// PARAMETERS: -// i, j - Cell coordinates -// -// RETURN VALUE: -// float - Enthalpy (dS) -// -// REVISION HISTORY -// $ 00sep23 : created LK -// #$ - -/* inline */ float galn_GetEnthalpy(GGAlign galn, int i, int j) -{ - float enthalpy; - - // first, need to determine optimum values for dG and dH - float tm0, tm1, tm2; - tm0 = nparam_CalcTM(nndS(i,j,0),nndH(i,j,0)); - tm1 = nparam_CalcTM(nndS(i,j,1),nndH(i,j,1)); - tm2 = nparam_CalcTM(nndS(i,j,2),nndH(i,j,2)); - if ((tm0>=tm1)&&(tm0>=tm2)) - enthalpy = nndH(i,j,0); - else if ((tm1>=tm0)&&(tm1>=tm2)) - enthalpy = nndH(i,j,1); - else - enthalpy = nndH(i,j,2); - - return enthalpy; -} - - - -/////////////////////////////////////////////////////////////////////////////// -// FUNCTION: float GetFreeEnergyK(int i, int j, float t); -// -// PURPOSE: Return free Energy for given cell at given temperature -// -// PARAMETERS: -// i, j - Cell coordinates -// t - Temperature in Kelvin -// -// RETURN VALUE: -// float - free Energy value (dG) at given temperature (Kelvin) -// -// REVISION HISTORY -// $ 00sep06 : created LK -// #$ - -float galn_GetFreeEnergyK(GGAlign galn, int i, int j, float t) -{ - // dG = dH - T * dS - - float entropy; - float enthalpy; - - // first, need to determine optimum values for dG and dH - entropy = -galn_GetEntropy(galn, i,j); - enthalpy = -galn_GetEnthalpy(galn, i,j); - - return enthalpy - t * entropy; -} - - -/////////////////////////////////////////////////////////////////////////////// -// FUNCTION: float GetFreeEnergyC(int i, int j, int t); -// -// PURPOSE: Return free Energy for given cell at given temperature -// -// PARAMETERS: -// i, j - Cell coordinates -// t - Temperature in Celsius -// -// RETURN VALUE: -// float - free Energy value (dG) at melting temperature -// -// REVISION HISTORY -// $ 00sep06 : created LK -// #$ - -float galn_GetFreeEnergyC(GGAlign galn, int i, int j, float t) -{ - // dG = dH - T * dS - - float entropy; - float enthalpy; - - // first, need to determine optimum values for dG and dH - entropy = galn_GetEntropy(galn, i,j); - enthalpy = galn_GetEnthalpy(galn, i,j); - - return enthalpy - (t + 273.15f) * entropy; -} - -/////////////////////////////////////////////////////////////////////////////// -// FUNCTION: float galn_GetMeltingTempC(int i, int j) -// -// PURPOSE: Calculate TM for given cell. Return in Celsius. -// -// PARAMETERS: -// i, j - Cell coordinates -// -// RETURN VALUE: -// float - Melting Temperature TM -// -// REVISION HISTORY -// $ 00sep06 : created LK -// #$ - -float galn_GetMeltingTempC(GGAlign galn, int i, int j) -{ - return galn_GetMeltingTempK(galn, i,j) - 273.15f; -} - - -/////////////////////////////////////////////////////////////////////////////// -// FUNCTION: float galn_GetMeltingTempK(int i, int j) -// -// PURPOSE: Calculate TM for given cell. Return in Kelvin. -// -// PARAMETERS: -// i, j - Cell coordinates -// -// RETURN VALUE: -// float - Melting Temperature TM -// -// REVISION HISTORY -// $ 00sep06 : created LK -// #$ - -float galn_GetMeltingTempK(GGAlign galn, int i, int j) -{ - float tm0, tm1, tm2; - tm0 = nparam_CalcTM(nndS(i,j,0),nndH(i,j,0)); - tm1 = nparam_CalcTM(nndS(i,j,1),nndH(i,j,1)); - tm2 = nparam_CalcTM(nndS(i,j,2),nndH(i,j,2)); - float maxg; - if ((tm0>=tm1)&&(tm0>=tm2)) - maxg=tm0; - else if ((tm1>=tm0)&&(tm1>=tm2)) - maxg=tm1; - else - maxg=tm2; - if (maxg<0) - maxg=0; - return maxg; -} - - diff --git a/src/libthermo/galign.h b/src/libthermo/galign.h deleted file mode 100644 index aad5541..0000000 --- a/src/libthermo/galign.h +++ /dev/null @@ -1,97 +0,0 @@ -/* - * galign.h - * PHunterLib - * - * Created by Tiayyba Riaz on 7/2/09. - * Copyright 2009 __MyCompanyName__. All rights reserved. - * - */ - -//============================================================================= -// Module: galign.h -// Project: Cubic Project - Calculation of melting temperature and free -// energy of two DNA strands -// Type: header file - Thermodynamic Alignment. -// 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 GAlign - Thermodynamic Alignment Algorithm -// Author: leber -// Date: 01/2002 - 02/2002 -// Copyright: (c) L. Kaderali & M. Leber, 01/2002 - 02/2002 -// -// Revision History -// $ 00sep04 LK : created -// $ 00dec30 LK : changed to do local alignment of probe against -// one entire sequence -// $ 01feb07 LK : optimized -// #$ -//============================================================================= - -#if !defined(AFX_GAlign_H__1B9227F7_82AB_11D4_9FFF_000000000000__INCLUDED_) -#define AFX_GAlign_H__1B9227F7_82AB_11D4_9FFF_000000000000__INCLUDED_ - -#if _MSC_VER > 1000 -#pragma once -#endif // _MSC_VER > 1000 - -//#include -//#include -//#include -#include "nnparams.h" // Nearest Neighbor Parameters -//#include "../gsfxtree/gsfxtree.h" // Suffix Tree Stuff - -//using namespace std; - -#ifdef _pack -#pragma pack(1) -#endif - -//----------------------------------------------------------------------------- -// class GAlign -//#pragma GCC visibility push(hidden) -//typedef class GAlign* GGAlign; - -typedef struct GAlign_st -{ -//public: -// lk01feb07: removed maxlocstuff as not requiered by thermtreealign... - float maxlocg; // maximum local dG value found - int maxloci; // i position thereof - int maxlocj; // j position thereof - int maxloct; // and type of maximum alignment! - int targetNumber; // id number of target sequence -// PGSfxLeaf probeNode; // identifier for probe sequence -//private: - char *seq1; // Sequence 1 - char *seq2; // Sequence 2 - int seq1len; // Length of Sequence 1 - int seq2len; // Length of Sequence 2 - int maxseq1len; // length of longest target... - float *dH; // Dynamic Programming Table for Entropy - float *dS; // Dynamic Programming Table for Enthalpy - PNNParams GNNParams; // Nearest Neighbor parameters -// float forbidden_entropy; -}GAlign, * GGAlign; - -void galn_initGAlign(GGAlign galn, int, int, PNNParams); -void galn_finiGAlign(GGAlign galn); -//void galn_InitStrings(GGAlign galn, char*,char*,int,int,int=0); -void galn_InitStrings(GGAlign galn, char*,char*,int,int,int); -void galn_InitBorder(GGAlign galn); -//void galn_CalculateTable(GGAlign galn, int=1); -void galn_CalculateTable(GGAlign galn, int); -float galn_GetEntropy(GGAlign galn, int, int); -float galn_GetEnthalpy(GGAlign galn, int, int); -float galn_GetFreeEnergyK(GGAlign galn, int, int, float); -float galn_GetFreeEnergyC(GGAlign galn, int, int, float); -float galn_GetMeltingTempC(GGAlign galn, int, int); -float galn_GetMeltingTempK(GGAlign galn, int, int); - -void galn_printEnthalpyTable(GGAlign galn, int level); -void galn_printEntropyTable(GGAlign galn, int level); - -//#pragma GCC visibility pop -#endif //!defined(AFX_GAlign_H__1B9227F7_82AB_11D4_9FFF_000000000000__INCLUDED_) - diff --git a/src/libthermo/libthermo.a b/src/libthermo/libthermo.a index cbedd26..be29e70 100644 Binary files a/src/libthermo/libthermo.a and b/src/libthermo/libthermo.a differ diff --git a/src/libthermo/libthermo.c b/src/libthermo/libthermo.c deleted file mode 100644 index 2fea2a4..0000000 --- a/src/libthermo/libthermo.c +++ /dev/null @@ -1,190 +0,0 @@ -/* - * PHunterLib.cp - * PHunterLib - * - * Created by Tiayyba Riaz on 6/7/09. - * Copyright 2009 __MyCompanyName__. All rights reserved. - * - */ - -#include -#include "libthermo.h" -#include "nnparams.h" -#include "thermalign.h" -#include "galign.h" -#include "dinkelbach.h" - - -//#define EXPORT __attribute__((visibility("default"))) - -//char targetsFile [1000]=""; -//char nonTargetsFile [1000]=""; -//char primersFile [1000]=""; -//char primersOutputFile [1000]=""; -//int minLength = DEF_MIN_PRIMER_LENGTH; -//int maxLength = DEF_MAX_PRIMER_LENGTH; -//int productMinLength = DEF_MIN_PROD_LENGTH; -//int productMaxLength = DEF_MAX_PROD_LENGTH; -//char forwardPrimer [PRIMERS_MAX_SIZE] = DEF_FORWARD_PRIMER; -//char reversePrimer [PRIMERS_MAX_SIZE] = DEF_REVERSE_PRIMER; -//int beginPosForward = DEF_BEGINPOS_FORWARD; -//int endPosForward = DEF_ENDPOS_FORWARD; -//int beginPosReverse = DEF_BEGINPOS_REVERSE; -//int endPosReverse = DEF_ENDPOS_REVERSE; -//char maskTargets [PRIMERS_MAX_SIZE] = DEF_MASK_TARGETS; -//char maskNonTargets [PRIMERS_MAX_SIZE] = DEF_MASK_NONTARGETS; -//char degeneracy[PRIMERS_MAX_SIZE]=DEF_DEGENERACY; -//int nseqForCandidates = DEF_SEQ_FOR_CANDS; -//float minCoverageTargets = DEF_MIN_COVERAGE_TARGETS; -//float maxCoverageNonTargets = DEF_MAX_COVERAGE_NONTARGETS; -//float maxSelfScore = DEF_MAX_SELF_SCORE; -//float maxEndScore = DEF_MAX_END_SCORE; -//float minGCContent = DEF_MIN_GCCONTENT; -//float maxGCContent = DEF_MAX_GCCONTENT; -//int gcClamp = DEF_GCCLAMP; -//int maxPolyX = DEF_MAX_POLY_X; -float concPrimers = DEF_CONC_PRIMERS; -float concSequences = DEF_CONC_SEQUENCES; -float salt = DEF_SALT; -int salMethod = SALT_METHOD_SANTALUCIA; -//float minTempTargets=DEF_MIN_TEMP_TARGETS; -//float maxTempTargets=DEF_MAX_TEMP_TARGETS; -//float maxTempNonTargets=DEF_MAX_TEMP_NONTARGETS; -//float deltaTempNonTargets=DEF_DELTA_TEMP_NONTARGETS; -//float maxPairTempDiff=DEF_MAX_PAIR_TEMP_DIFF; -//char primersLabel [1000] = DEF_PRIMERS_LABEL; -//int fullStats = DEF_FULL_STATS; -//char seqRank [10000] = DEF_SEQ_RANK; - - -int countGCContent(char * seq ) { - int lseq = strlen(seq); - int k; - int count = 0; - for( k=0;k0;i--) { - seq[i]=seq[i-1]; - } - seq[0] = '$'; -} - -void removeDollarSigns(char * seq) { - int n = strlen(seq); - int i; - for( i=0;i -#include -//#include -#include -#include - -#include "thermalign.h" - - - #define dH(a,b,c) (*(thrma->dH+((b)*(thrma->maxseq1len+1)*3+(a)*3+(c)))) - #define dS(a,b,c) (*(thrma->dS+((b)*(thrma->maxseq1len+1)*3+(a)*3+(c)))) - - -/////////////////////////////////////////////////////////////////////////////// -// Construction/Destruction -/////////////////////////////////////////////////////////////////////////////// - -//------------------------ Constructor ---------------------------------------- -// Initialize Thermalign Object. This will allocate memory for the dynamic -// programming matrices for both entropy and enthalpy of the present -// alignment. The max length of the two sequences to align has to be passed -// to the routine. Also, an initialized CNNParams Object must be passed to -// CThermAlign. Object will be initialized only once, sequences can be -// passed to it later. - -//int s1lmax,s2lmax; // global requiered for DEBUGGING - -void therm_initCThermAlign(PThermAlign thrma, int ALength, int BLength, PNNParams myNNParams) { - // Create dynamic programming matrices for entropy and enthalpy. - // The first two dimensions are for the two strings respectively, the - // third dimensions differentiates the alignment according to its end: - // type 0 Alignments will be those that align x(i+1) with y(j+1), - // type 1 Alignments align x(i+1) with a gap - // type 2 Alignments align y(i+1) with a gap - // All three of these types must separately be considered to select the - // optimum alignment in the next iteration. - // The alignment algorithm proceeds analogous to the well known - // Smith-Waterman algorithm for global alignments. - - thrma->dH = (float*)calloc ((ALength+1)*(BLength+1)*3, sizeof (float)); - thrma->dS = (float*)calloc ((ALength+1)*(BLength+1)*3, sizeof (float)); -// s1lmax=ALength; -// s2lmax=BLength; - - if ((thrma->dH==NULL)||(thrma->dS==NULL)) - { - printf("Fatal Error: Out of memory!\n"); - exit(-1); - } - - thrma->NNParams = myNNParams; - thrma->seq1len = ALength; - thrma->maxseq1len = thrma->seq1len; - thrma->seq2len = BLength; - - // set dH / dS entries all to zero! - #ifdef _pack - // this is FAST, but works only if #pragma pack(1) has been set... - memset(dH,0,(ALength+1)*(BLength+1)*3*sizeof(float)); - memset(dS,0,(ALength+1)*(BLength+1)*3*sizeof(float)); - #else - // the following is slower, but works also when padding bytes have - // been inserted in the array for alignment by the compiler... - // time is not really critical here, as this initialization is done - // only once... - int a, b; - for (a=0;a<=ALength;a++) - for (b=0;b<=BLength;b++) - { - dH(a,b,0)=0.0f; - dH(a,b,1)=0.0f; - dH(a,b,2)=0.0f; - dS(a,b,0)=0.0f; - dS(a,b,1)=0.0f; - dS(a,b,2)=0.0f; - } - #endif - therm_InitBorder(thrma); // also need to do only once... -// forbidden_entropy = (-(NNParams->R)*(float)log((NNParams->Ct)/4.0f)); - return; -} - -//------------------------ Destructor ----------------------------------------- -// Free memory - -void therm_finiCThermAlign(PThermAlign thrma) -{ - free (thrma->dH); - free (thrma->dS); - return; -} - - - -/////////////////////////////////////////////////////////////////////////////// -// FUNCTION: void CThermAlign::InitStrings(char* s1, -// char* s2, -// int s1len, -// int s2len); -// -// PURPOSE: Initialize strings for alignment table! -// -// PARAMETERS: -// s1: Target sequence -// s2: Probe sequence -// s1len: Target sequence length -// s2len: Probe sequence length -// -// RETURN VALUE: -// void -// -// REVISION HISTORY -// $ 01jan03 LK : created -// $ 01feb07 LK : removed maximum local stuff -// #$ - -void therm_InitStrings(PThermAlign thrma, char* s1,char* s2,int s1len ,int s2len,int minxpos) -{ - thrma->seq1=s1; thrma->seq1len=s1len; - thrma->seq2=s2; thrma->seq2len=s2len; - nparam_UpdateParams(thrma->NNParams, s1,s2); -/* lk01feb07: removed maxloc-stuff - // if maximum local is not in reused subtable, need to reset! - if ((maxlocj>=minxpos)||(maxloci!=seq1len)) - { - maxloci=0; maxlocj=0; maxloctm=0; - } -*/ -// InitBorder(); //lk01feb07: removed, as it suffices if this is done once in constructor - return; -} - - -/////////////////////////////////////////////////////////////////////////////// -// FUNCTION: void CThermAlign::InitBorder(); -// -// PURPOSE: Initialize the lower and left border of the dynamic -// programming table. -// -// PARAMETERS: -// none -// -// RETURN VALUE: -// void -// -// REVISION HISTORY -// $ 00sep06 LK : created -// $ 00dec30 LK : complete revision - simply init table now, no -// more saving etc... -// $ 01feb07 LK : removed local alignment stuff and optimized -// #$ - -void therm_InitBorder(PThermAlign thrma) -{ - // Initialization for new alignment! - thrma->maxloctm = 0; // maximum melting temperature found (kelvin) - thrma->maxloci = 0; - thrma->maxlocj = 0; - thrma->maxloct = 0; - int i, j; - - for ( i=0; i<=thrma->seq1len; i++) - { - dH(i,0,0) = 0.0f; dH(i,0,1) = 0.0f; dH(i,0,2) = 0.0f; - dS(i,0,0) = nparam_GetInitialEntropy(thrma->NNParams); dS(i,0,1) = nparam_GetInitialEntropy(thrma->NNParams); - dS(i,0,2) = nparam_GetInitialEntropy(thrma->NNParams); - } - for ( j=1; j<=thrma->seq2len; j++) - { - dH(0,j,0) = 0.0f; dH(0,j,1) = 0.0f; dH(0,j,2) = 0.0f; - dS(0,j,0) = nparam_GetInitialEntropy(thrma->NNParams); dS(0,j,1) = nparam_GetInitialEntropy(thrma->NNParams); - dS(0,j,2) = nparam_GetInitialEntropy(thrma->NNParams); - } - - return; -} - - -/////////////////////////////////////////////////////////////////////////////// -// FUNCTION: void CThermAlign::CalculateCell(int i, int j); -// -// PURPOSE: Calculate entropy and enthalpy for a given cell using a variant -// of the dynamic programming recursion, choosing greedy for -// tm->max -// x(i,j,0) = min/max {x(i-1,j-1,t) | for t in {0,1,2}} + newnn -// x(i,j,1) = min/max {x(i-1,j-1,0) + newnn, x(i-1,j,1)} -// x(i,j,2) = min/max {x(i-1,j-1,0) + newnn, x(i, j-1,2)} -// where newnn is the contribution of the newly added nearest -// neighbor term. The underlying assumption is that two -// consecutive gaps do not contribute to x. (Can however be -// easily changed, as long as an affine relationship exists -// between the length of the gap and the associated "cost"). -// To estimate the melting temperature TM from Enthalpy and -// Entropy, the following formula is being used: -// TM = dH / (dS + R ln CT) -// where -// TM is the temperature at which 50% of the duplex -// are dissociated -// dH is the Enthalpy -// dS is the Entropy -// R is the gas constant -// CT is the concentration (assumed constant) -// and ln means the natural logarithm. -// One drawback is that the two separate computations of -// Entropy and Enthalpy may result in different alignments. -// We do not want that, and thus pick the one with the higher -// melting temperature TM. However, it must be warned that -// depending on the underlying nearest neighbor parameters this -// may only approximate the "true" optimum result! -// -// PARAMETERS: -// int i, -// int j - coordinates of the dp cell to calculate (=position in -// string). Requieres that (i-1,j), (i,j-1), (i-1,j-1) -// are known. -// -// RETURN VALUE: -// - will be written directly to dH/dS array in this object -// -// REVISION HISTORY -// $ 00sep06 LK : created -// $ 00dec30 LK : modified to align probe with full target (local -// alignment) -// $ 01feb07 LK : optimized and removed local alignment stuff -// $ 01feb08 LK : optimized further and merged into CalculateTable -// #$ -/////////////////////////////////////////////////////////////////////////////// -// FUNCTION: CThermAlign::CalculateTable(); -// -// PURPOSE: Compute the entire dynamic programming table cell by cell. -// Requieres that InitBorder has been called. -// Computations will be done in the usual rowwise manner. -// -// PARAMETERS: -// int startRow: Row that computation shall commence at! -// -// RETURN VALUE: -// void - Data will be written directly to dH and dS in this object -// -// REVISION HISTORY -// $ 00sep06 : created LK -// 00dec31 : modified to begin at a specified row! -// #$ - -void therm_CalculateTable(PThermAlign thrma, int startRow) -{ -// #define forbidden_enthalpy 1000000000000000000.0f - register int i; - register int iminusone; - register int j; - register int jminusone; - // Determine previous character! - register char prevchari; - register char prevcharj; - register char seq1char; - register char seq2char; - register float entropy1; - register float entropy2; - register float entropy3; - register float enthalpy1; - register float enthalpy2; - register float enthalpy3; - register float tm1; - register float tm2; - register float tm3; - float maxtm; -// ofstream dboutstream("debug.out"); -// assert(startRow>0); - for (j=startRow; j<=thrma->seq2len; j++) - { - if (j>1) - prevcharj = thrma->seq2[j-2]; - else - prevcharj = 0; - jminusone = j - 1; - seq2char = thrma->seq2[jminusone]; - for (i=1; i<=thrma->seq1len; i++) - { - // local variables - if (i>1) - prevchari = thrma->seq1[i-2]; - else - prevchari = 0; - iminusone = i - 1; - seq1char=thrma->seq1[iminusone]; - - // --------------- Upper cell: Alignment of seq1[i] with seg2[j] - entropy1 = dS(iminusone,jminusone,0); - enthalpy1 = dH(iminusone,jminusone,0) ; - entropy2 = dS(iminusone,jminusone,1); - enthalpy2 = dH(iminusone,jminusone,1); - entropy3 = dS(iminusone,jminusone,2); - enthalpy3 = dH(iminusone,jminusone,2); - if (enthalpy1NNParams, prevchari, seq1char, prevcharj, seq2char); - enthalpy1 += nparam_GetEnthalpy(thrma->NNParams, prevchari, seq1char, prevcharj, seq2char); - } - if (enthalpy2NNParams, prevchari, seq1char, 0, seq2char); // 0 is gap! - entropy2 += nparam_GetEntropy(thrma->NNParams, prevchari, seq1char, 0, seq2char); - } - if (enthalpy3NNParams, 0, seq1char, prevcharj, seq2char); - entropy3 += nparam_GetEntropy(thrma->NNParams, 0, seq1char, prevcharj, seq2char); - } - // choose optimum combination - tm1 = nparam_CalcTM(entropy1, enthalpy1); - tm2 = nparam_CalcTM(entropy2, enthalpy2); - tm3 = nparam_CalcTM(entropy3, enthalpy3); - if ((tm1>=tm2)&&(tm1>=tm3)) - { - dS(i,j,0) = entropy1; - dH(i,j,0) = enthalpy1; - maxtm=tm1; - } - else if (tm2>=tm3) - { - dS(i,j,0) = entropy2; - dH(i,j,0) = enthalpy2; - maxtm=tm2; - } - else - { - dS(i,j,0) = entropy3; - dH(i,j,0) = enthalpy3; - maxtm=tm3; - } - // set to zero if alignment so far is mismatches only (tm will be zero -> no need to reset maxtm)! - if (dH(i,j,0)==0) - dS(i,j,0)=nparam_GetInitialEntropy(thrma->NNParams); - // check if local maximum found! - if (((i==thrma->seq1len)||(j==thrma->seq2len))&&(maxtm>=thrma->maxloctm)) - { - thrma->maxloctm=maxtm; - thrma->maxloci=i; - thrma->maxlocj=j; - thrma->maxloct=0; - } - // set to zero if alignment so far is mismatches only! - if (dH(i,j,0)==0) - dS(i,j,0)=nparam_GetInitialEntropy(thrma->NNParams); - - // --------------- Middle cell: Alignment of seq1[i] with '-' - // following changes lk 01jan08 - // we do not allow $- in the beginning, therefore set to forbidden if - // j=1 - if (j==1) - { - enthalpy1=forbidden_enthalpy; enthalpy2=forbidden_enthalpy; - entropy1=forbidden_entropy; entropy2=forbidden_entropy; - } - // also, disallow -$ in the end, therefore forbid if j=seq2len-1 - else if (j==thrma->seq2len-1) - { - enthalpy1=forbidden_enthalpy; enthalpy2=forbidden_enthalpy; - entropy1=forbidden_entropy; entropy2=forbidden_entropy; - } - else - { - // end lk01jan08 - // -- entropy - entropy1 = dS(iminusone,j,0); - entropy2 = dS(iminusone,j,1); - enthalpy1 = dH(iminusone,j,0); - enthalpy2 = dH(iminusone,j,1); - if (enthalpy1NNParams, prevchari, seq1char, seq2char, 0); - enthalpy1 += nparam_GetEnthalpy(thrma->NNParams, prevchari, seq1char, seq2char, 0); - } - if (enthalpy2NNParams, prevchari, seq1char, 0, 0); - enthalpy2 += nparam_GetEnthalpy(thrma->NNParams, prevchari, seq1char, 0, 0); - } - } - // -- choose optimum combination - tm1 = nparam_CalcTM(entropy1, enthalpy1); - tm2 = nparam_CalcTM(entropy2, enthalpy2); - if (tm1>=tm2) - { - dS(i,j,1)=entropy1; - dH(i,j,1)=enthalpy1; - maxtm=tm1; - } - else - { - dS(i,j,1)=entropy2; - dH(i,j,1)=enthalpy2; - maxtm=tm2; - } - // check if local maximum found! - if (((i==thrma->seq1len)||(j==thrma->seq2len))&&(maxtm>=thrma->maxloctm)) - { - thrma->maxloctm=maxtm; - thrma->maxloci=i; - thrma->maxlocj=j; - thrma->maxloct=1; - } - // set to zero if alignment so far is mismatches only! - if (dH(i,j,1)==0) - dS(i,j,1)=nparam_GetInitialEntropy(thrma->NNParams); - - // --------------- Lower cell: Alignment of '-' with seq2[j] - // following changes lk 01jan08 - // we do not allow $- in the beginning, therefore set to forbidden if - // i=1 - if (i==1) - { - enthalpy1=forbidden_enthalpy; enthalpy2=forbidden_enthalpy; - entropy1=forbidden_entropy; entropy2=forbidden_entropy; - } - // also, disallow -$ in the end, therefore forbid if i=seq1len-1 - else if (i==thrma->seq1len-1) - { - enthalpy1=forbidden_enthalpy; enthalpy2=forbidden_enthalpy; - entropy1=forbidden_entropy; entropy2=forbidden_entropy; - } - else - { - // end lk01jan08 - entropy1 = dS(i,jminusone,0); - entropy2 = dS(i,jminusone,2); - enthalpy1 = dH(i,jminusone,0); - enthalpy2 = dH(i,jminusone,2); - if (enthalpy1NNParams, seq1char, 0, prevcharj, seq2char); - enthalpy1 += nparam_GetEnthalpy(thrma->NNParams, seq1char, 0, prevcharj, seq2char); - } - if (enthalpy2NNParams, 0, 0, prevcharj, seq2char); - enthalpy2 += nparam_GetEnthalpy(thrma->NNParams, 0, 0, prevcharj, seq2char); - } - } - // -- choose optimum combination - tm1 = nparam_CalcTM(entropy1, enthalpy1); - tm2 = nparam_CalcTM(entropy2, enthalpy2); - if (tm1>=tm2) - { - dS(i,j,2)=entropy1; - dH(i,j,2)=enthalpy1; - maxtm=tm1; - } - else - { - dS(i,j,2)=entropy2; - dH(i,j,2)=enthalpy2; - maxtm=tm2; - } - // check if local maximum found! - if (((i==thrma->seq1len)||(j==thrma->seq2len))&&(maxtm>=thrma->maxloctm)) - { - thrma->maxloctm=maxtm; - thrma->maxloci=i; - thrma->maxlocj=j; - thrma->maxloct=2; - } - // set to zero if alignment so far is mismatches only! - if (dH(i,j,2)==0) - dS(i,j,2)=nparam_GetInitialEntropy(thrma->NNParams); - - - /* - dboutstream<<"After cell ("<=tm1)&&(tm0>=tm2)) - entropy = dS(i,j,0); - else if ((tm1>=tm0)&&(tm1>=tm2)) - entropy = dS(i,j,1); - else - entropy = dS(i,j,2); - - return entropy; -} - - -/////////////////////////////////////////////////////////////////////////////// -// FUNCTION: float GetEnthalpy(int i, int j); -// -// PURPOSE: Return Enthalpy for given cell -// -// PARAMETERS: -// i, j - Cell coordinates -// -// RETURN VALUE: -// float - Enthalpy (dS) -// -// REVISION HISTORY -// $ 00sep23 : created LK -// #$ - -float therm_GetEnthalpy(PThermAlign thrma, int i, int j) -{ - float enthalpy; - - // first, need to determine optimum values for dG and dH - float tm0, tm1, tm2; - tm0 = nparam_CalcTM(dS(i,j,0),dH(i,j,0)); - tm1 = nparam_CalcTM(dS(i,j,1),dH(i,j,1)); - tm2 = nparam_CalcTM(dS(i,j,2),dH(i,j,2)); - if ((tm0>=tm1)&&(tm0>=tm2)) - enthalpy = dH(i,j,0); - else if ((tm1>=tm0)&&(tm1>=tm2)) - enthalpy = dH(i,j,1); - else - enthalpy = dH(i,j,2); - - return enthalpy; -} - -/////////////////////////////////////////////////////////////////////////////// -// FUNCTION: float GetFreeEnergy(int i, int j); -// -// PURPOSE: Return free Energy for given cell at melting temp -// -// PARAMETERS: -// i, j - Cell coordinates -// -// RETURN VALUE: -// float - free Energy value (dG) at melting temperature -// -// REVISION HISTORY -// $ 00sep06 : created LK -// #$ - -float therm_GetFreeEnergy(PThermAlign thrma, int i, int j) -{ - return therm_GetFreeEnergyK(thrma, i,j,therm_GetMeltingTempK(thrma, i,j)); -} - - -/////////////////////////////////////////////////////////////////////////////// -// FUNCTION: float GetFreeEnergyK(int i, int j, float t); -// -// PURPOSE: Return free Energy for given cell at given temperature -// -// PARAMETERS: -// i, j - Cell coordinates -// t - Temperature in Kelvin -// -// RETURN VALUE: -// float - free Energy value (dG) at given temperature (Kelvin) -// -// REVISION HISTORY -// $ 00sep06 : created LK -// #$ - -float therm_GetFreeEnergyK(PThermAlign thrma, int i, int j, float t) -{ - // dG = dH - T * dS - - float entropy; - float enthalpy; - - // first, need to determine optimum values for dG and dH - entropy = therm_GetEntropy(thrma, i,j); - enthalpy = therm_GetEnthalpy(thrma, i,j); - - return enthalpy - t * entropy; -} - - -/////////////////////////////////////////////////////////////////////////////// -// FUNCTION: float GetFreeEnergyC(int i, int j, int t); -// -// PURPOSE: Return free Energy for given cell at given temperature -// -// PARAMETERS: -// i, j - Cell coordinates -// t - Temperature in Celsius -// -// RETURN VALUE: -// float - free Energy value (dG) at melting temperature -// -// REVISION HISTORY -// $ 00sep06 : created LK -// #$ - -float therm_GetFreeEnergyC(PThermAlign thrma, int i, int j, float t) -{ - // dG = dH - T * dS - - float entropy; - float enthalpy; - - // first, need to determine optimum values for dG and dH - entropy = therm_GetEntropy(thrma, i,j); - enthalpy = therm_GetEnthalpy(thrma, i,j); - - return enthalpy - (t + 273.15f) * entropy; -} - -/////////////////////////////////////////////////////////////////////////////// -// FUNCTION: float CThermAlign::GetMeltingTempC(int i, int j) -// -// PURPOSE: Calculate TM for given cell. Return in Celsius. -// -// PARAMETERS: -// i, j - Cell coordinates -// -// RETURN VALUE: -// float - Melting Temperature TM -// -// REVISION HISTORY -// $ 00sep06 : created LK -// #$ - -float therm_GetMeltingTempC(PThermAlign thrma, int i, int j) -{ - return therm_GetMeltingTempK(thrma, i,j) - 273.15f; -} - - -/////////////////////////////////////////////////////////////////////////////// -// FUNCTION: float CThermAlign::GetMeltingTempK(int i, int j) -// -// PURPOSE: Calculate TM for given cell. Return in Kelvin. -// -// PARAMETERS: -// i, j - Cell coordinates -// -// RETURN VALUE: -// float - Melting Temperature TM -// -// REVISION HISTORY -// $ 00sep06 : created LK -// #$ - -float therm_GetMeltingTempK(PThermAlign thrma, int i, int j) -{ - float tm0, tm1, tm2; - tm0 = nparam_CalcTM(dS(i,j,0), dH(i,j,0)); - tm1 = nparam_CalcTM(dS(i,j,1), dH(i,j,1)); - tm2 = nparam_CalcTM(dS(i,j,2), dH(i,j,2)); - float maxtm; - if ((tm0>=tm1)&&(tm0>=tm2)) - maxtm=tm0; - else if ((tm1>=tm0)&&(tm1>=tm2)) - maxtm=tm1; - else - maxtm=tm2; - if (maxtm<0) - maxtm=0; - return maxtm; -} - diff --git a/src/libthermo/thermalign.h b/src/libthermo/thermalign.h deleted file mode 100644 index 65f368a..0000000 --- a/src/libthermo/thermalign.h +++ /dev/null @@ -1,105 +0,0 @@ -//============================================================================= -// Module: thermalign.h -// Project: Diploma Thesis - Probe Selection for DNA Microarrays -// Type: header file - Thermodynamic Alignment. -// 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 CThermAlign - Thermodynamic Alignment Algorithm -// Author: kaderali -// Date: 9/2000 - 12/2000 -// Copyright: (c) L. Kaderali, 9/2000 - 12/2000 -// -// Revision History -// $ 00sep04 LK : created -// $ 00dec30 LK : changed to do local alignment of probe against -// one entire sequence -// $ 01feb07 LK : optimized -// #$ -//============================================================================= - -#if !defined(AFX_THERMALIGN_H__1B9227F7_82AB_11D4_9FFF_000000000000__INCLUDED_) -#define AFX_THERMALIGN_H__1B9227F7_82AB_11D4_9FFF_000000000000__INCLUDED_ - -//#if _MSC_VER > 1000 -//#pragma once -//#endif // _MSC_VER > 1000 - -//#include -//#include -//#include -#include "nnparams.h" // Nearest Neighbor Parameters -//#include "../gsfxtree/gsfxtree.h" // Suffix Tree Stuff - -//using namespace std; - -//#ifdef _pack -//#pragma pack(1) -//#endif - -//----------------------------------------------------------------------------- -// class CThermAlign -//#pragma GCC visibility push(hidden) -//typedef class CThermAlign* PThermAlign; - -typedef struct CThermAlign_t -{ -//public: - #ifdef _output_alignment -/* void printAlignment(ostream &outputStream); - void OutputAlignment(ostream &outputStream); - bool OutputAlignment(ostream &outputStream, int, int, int, bool local=false); - void OutputLocalAlignment(ostream &outputStream); - void PrintDPTable(ostream&);*/ - #endif -// lk01feb07: removed maxlocstuff as not requiered by thermtreealign... - float maxloctm; // maximum local temperature found - int maxloci; // i position thereof - int maxlocj; // j position thereof - int maxloct; // and type of maximum alignment! - int targetNumber; // id number of target sequence -// PGSfxLeaf probeNode; // identifier for probe sequence - char *seq1; // Sequence 1 - char *seq2; // Sequence 2 - int seq1len; // Length of Sequence 1 - int seq2len; // Length of Sequence 2 - int maxseq1len; // length of longest target... -//private: - float *dH; // Dynamic Programming Table for Entropy - float *dS; // Dynamic Programming Table for Enthalpy - PNNParams NNParams; // Nearest Neighbor parameters -// float forbidden_entropy; -#ifdef _output_alignment - //ostrstream *s1aptr, *s2aptr, *atypptr; -// Used to buffer aligned sequences (on output) - //ostrstream s1align; - //ostrstream s2align; - //ostrstream aligntype; // insert, deletion, match, unmatch (for output) - // same for local alignment: - /* removed lk00jan08: use global instead! - ostrstream ls1align; - ostrstream ls2align; - ostrstream laligntype; // insert, deletion, match, unmatch (for output) - */ -#endif -}CThermAlign, * PThermAlign; -//#pragma GCC visibility pop - -float therm_GetEntropy(PThermAlign thrma, int, int); -float therm_GetEnthalpy(PThermAlign thrma, int, int); -float therm_GetFreeEnergy(PThermAlign thrma, int, int); -float therm_GetFreeEnergyK(PThermAlign thrma, int, int, float); -float therm_GetFreeEnergyC(PThermAlign thrma, int, int, float); -float therm_GetMeltingTempC(PThermAlign thrma, int, int); -float therm_GetMeltingTempK(PThermAlign thrma, int, int); -void therm_initCThermAlign(PThermAlign thrma, int, int, PNNParams); -void therm_finiCThermAlign(PThermAlign thrma); -//void therm_InitStrings(PThermAlign thrma, char*,char*,int,int,int=0); -void therm_InitStrings(PThermAlign thrma, char*,char*,int,int,int); -void therm_InitBorder(PThermAlign thrma); -//void therm_CalculateTable(PThermAlign thrma, int=1); -void therm_CalculateTable(PThermAlign thrma, int); - - -#endif // !defined(AFX_THERMALIGN_H__1B9227F7_82AB_11D4_9FFF_000000000000__INCLUDED_)