obitools/ecoprimers
2
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 2 Wiki Activity

Deleted some supurious files

Browse Source 
git-svn-id: https://www.grenoble.prabi.fr/svn/LECASofts/ecoPrimers/trunk@227 60f365c0-8329-0410-b2a4-ec073aeeaa1d
This commit is contained in:
Eric Coissac
2009-07-17 22:29:19 +00:00
parent b625941d72
commit 455bf63949
9 changed files with 0 additions and 1980 deletions
Download Patch File Download Diff File Expand all files Collapse all files

76
src/libthermo/dinkelbach.c
View File

@ -1,76 +0,0 @@
/*
* dinkelbach.cpp
* PHunterLib
*
* Created by Tiayyba Riaz on 6/7/09.
* Copyright 2009 __MyCompanyName__. All rights reserved.
*
*/
//#include <iostream>
//#include <strstream>
//#include <cstring>
#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<oldlambda)));
if(newTM - oldTM > 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;
}

36
src/libthermo/dinkelbach.h
View File

@ -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 <iostream>
//#include <strstream>
//#include <cstring>
#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 */

685
src/libthermo/galign.c
View File

@ -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 <math.h>
#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#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 (enthalpy1<forbidden_enthalpy)
{
entropy1 += nparam_GetEntropy(galn->GNNParams, prevchari,seq1char, prevcharj, seq2char);
enthalpy1 +=nparam_GetEnthalpy(galn->GNNParams, prevchari, seq1char, prevcharj, seq2char);
}
if (enthalpy2<forbidden_enthalpy)
{
enthalpy2 +=nparam_GetEnthalpy(galn->GNNParams, prevchari,seq1char, 0, seq2char); // 0 is gap!
entropy2 += nparam_GetEntropy(galn->GNNParams, prevchari,seq1char, 0, seq2char);
}
if(enthalpy3<forbidden_enthalpy)
{
enthalpy3 += nparam_GetEnthalpy(galn->GNNParams, 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 (enthalpy1<forbidden_enthalpy)
{
entropy1+=nparam_GetEntropy(galn->GNNParams, prevchari, seq1char, seq2char, 0);
enthalpy1 +=nparam_GetEnthalpy(galn->GNNParams, prevchari,seq1char, seq2char, 0);
}
if (enthalpy2<forbidden_enthalpy)
{
entropy2 +=nparam_GetEntropy(galn->GNNParams, 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 (enthalpy1<forbidden_enthalpy)
{
entropy1+=nparam_GetEntropy(galn->GNNParams, seq1char, 0, prevcharj, seq2char);
enthalpy1 +=nparam_GetEnthalpy(galn->GNNParams, seq1char, 0, prevcharj, seq2char);
}
if(enthalpy2<forbidden_enthalpy)
{
entropy2 += nparam_GetEntropy(galn->GNNParams, 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 ("<<i<<","<<j<<"):"<<endl;
PrintDPTable(dboutstream);
dboutstream<<endl<<flush;
*/
}
}
//printEnthalpyTable(0);
//printEntropyTable(0);
return;
}
void galn_printEnthalpyTable(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");
}
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;
}

97
src/libthermo/galign.h
View File

@ -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 <iostream>
//#include <strstream>
//#include <ostream>
#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_)

BIN
src/libthermo/libthermo.a
View File

Binary file not shown.

190
src/libthermo/libthermo.c
View File

@ -1,190 +0,0 @@
/*
* PHunterLib.cp
* PHunterLib
*
* Created by Tiayyba Riaz on 6/7/09.
* Copyright 2009 __MyCompanyName__. All rights reserved.
*
*/
#include <stdlib.h>
#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;k<lseq;k++) {
if (seq[k] == 'G' || seq[k] == 'C' ) {
count+=1;
}
}
return count;
}
char getComplement(char base) {
if(base == 'A') return 'T';
else if(base == 'G') return 'C';
else if(base == 'C') return 'G';
else if(base == 'T') return 'A';
else if(base == 'N') return 'N';
return 'A';
}
void addDollarSigns(char * seq) {
int n = strlen(seq);
int i;
seq[n+2]=0;
seq[n+1]='$';
for( i=n;i>0;i--) {
seq[i]=seq[i-1];
}
seq[0] = '$';
}
void removeDollarSigns(char * seq) {
int n = strlen(seq);
int i;
for( i=0;i<n-2;i++) {
seq[i]=seq[i+1];
}
seq[n-2]=0;
}
float calculateMeltingTemperatureInt (char * seq1, char * seq2, PNNParams params) {
int length = strlen(seq2);
int perfectComplements = 1;
int k;
if(strlen(seq1)!= length) {
perfectComplements = 0;
}
for( k=0;k<length && perfectComplements;k++) {
if(seq2[k] != getComplement(seq1[k])) {
perfectComplements = 0;
}
}
addDollarSigns(seq1);
addDollarSigns(seq2);
CThermAlign myAlignment;//(strlen(seq1),strlen(seq2),params);
therm_initCThermAlign(&myAlignment, strlen(seq1),strlen(seq2),params);
therm_InitStrings (&myAlignment, seq1,seq2,strlen(seq1),strlen(seq2), 0);
therm_InitBorder(&myAlignment);
therm_CalculateTable(&myAlignment, 1);
float answer;
if(perfectComplements) {
answer = therm_GetMeltingTempC(&myAlignment, myAlignment.maxloci,myAlignment.maxlocj);
} else {
float tempK = therm_GetMeltingTempK (&myAlignment, myAlignment.maxloci,myAlignment.maxlocj);
GAlign mydGAlign;//(strlen(seq1),strlen(seq2),params);
galn_initGAlign(&mydGAlign, strlen(seq1),strlen(seq2),params);
galn_InitStrings(&mydGAlign, seq1,seq2,strlen(seq1),strlen(seq2), 0);
dinkelbach myDinkel;//(params, &mydGAlign);
dnkl_initdinkelbach(&myDinkel, params, &mydGAlign);
dnkl_iteration(&myDinkel, tempK);
answer = galn_GetMeltingTempC(&mydGAlign, mydGAlign.maxloci,mydGAlign.maxlocj);
galn_finiGAlign(&mydGAlign);
}
removeDollarSigns(seq1);
removeDollarSigns(seq2);
therm_finiCThermAlign(&myAlignment);
return answer;
}
float calculateMeltingTemperature (char * seq1, char * seq2) {
long len = strlen (seq1);
char *seq1new = (char *) calloc(len + 3, sizeof(char));//new char[len + 3];
len = strlen (seq2);
char *seq2new = (char *) calloc(len + 3, sizeof(char));//new char[len + 3];
strcpy (seq1new, seq1);
strcpy (seq2new, seq2);
PNNParams paramsPrimerSeq = (PNNParams) calloc (1, sizeof(CNNParams));
nparam_InitParams(paramsPrimerSeq, concPrimers,concSequences,salt,salMethod);
float temp = calculateMeltingTemperatureInt(seq1new, seq2new, paramsPrimerSeq);
free (seq1new);
free (seq2new);
return temp;
}
float calculateMeltingTemperatureBasic (char * seq) {
int gccount;
float temp;
int seqlen;
seqlen = strlen (seq);
gccount = countGCContent (seq);
temp = 64.9 + 41*(gccount - 16.4)/seqlen;
return temp;
}

75
src/libthermo/libthermo.h
View File

@ -1,75 +0,0 @@
/*
* PHunterLib.h
* PHunterLib
*
* Created by Tiayyba Riaz on 7/2/09.
* Copyright 2009 __MyCompanyName__. All rights reserved.
*
*/
#ifndef PHunterLib_
#define PHunterLib_
#define TARGETS_MAX_SIZE 2000
#define NON_TARGETS_MAX_SIZE 10000
#define SEQUENCES_MAX_SIZE 100000
#define PRIMERS_MAX_SIZE 100
#define CANDIDATES_MAX_SIZE 100000
#define SEEDS_MAX_SIZE 9
#define SEEDSTABLE_TAM 1000000
#define MAX_DEGENERACY 256
#define STAT_COUNTS 12
#define STAT_TOTAL_CANDIDATES 0
#define STAT_TEST_DEGENERACY 1
#define STAT_TEST_GCCONTENT 2
#define STAT_TEST_GCCLAMP 3
#define STAT_TEST_MAX_POLY_X 4
#define STAT_TEST_SELFCOMP 5
#define STAT_TEST_TMSELF 6
#define STAT_TEST_TARGETS 7
#define STAT_TEST_NONTARGETS 8
#define STAT_ALLTESTS 9
#define STAT_FAILED_SELFCOMP_LOW 10
#define STAT_FAILED_SELFCOMP_HIGH 11
#define DEF_MIN_PRIMER_LENGTH 20
#define DEF_MAX_PRIMER_LENGTH 25
#define DEF_MIN_PROD_LENGTH 75
#define DEF_MAX_PROD_LENGTH 200
#define DEF_FORWARD_PRIMER "NONE"
#define DEF_REVERSE_PRIMER "NONE"
#define DEF_BEGINPOS_FORWARD 0
#define DEF_ENDPOS_FORWARD SEQUENCES_MAX_SIZE
#define DEF_BEGINPOS_REVERSE 0
#define DEF_ENDPOS_REVERSE SEQUENCES_MAX_SIZE
#define DEF_MASK_TARGETS "11"
#define DEF_MASK_NONTARGETS "NONE"
#define DEF_DEGENERACY "1"
#define DEF_SEQ_FOR_CANDS 1
#define DEF_MIN_COVERAGE_TARGETS 100
#define DEF_MAX_COVERAGE_NONTARGETS 0
#define DEF_MAX_SELF_SCORE 800
#define DEF_MAX_END_SCORE 300
#define DEF_MIN_GCCONTENT 25
#define DEF_MAX_GCCONTENT 75
#define DEF_GCCLAMP 0
#define DEF_MAX_POLY_X 5
#define DEF_CONC_PRIMERS 0.0000008
#define DEF_CONC_SEQUENCES 0
#define DEF_SALT 0.05
#define DEF_MIN_TEMP_TARGETS 40
#define DEF_MAX_TEMP_TARGETS 70
#define DEF_MAX_TEMP_NONTARGETS 40
#define DEF_DELTA_TEMP_NONTARGETS 0
#define DEF_MAX_PAIR_TEMP_DIFF 40
#define DEF_PRIMERS_LABEL "P"
#define DEF_FULL_STATS 0
#define DEF_SEQ_RANK "NONE"
float calculateMeltingTemperature (char * seq1, char * seq2);
float calculateMeltingTemperatureBasic (char * seq);
int countGCContent(char * seq );
#endif

716
src/libthermo/thermalign.c
View File

@ -1,716 +0,0 @@
/*
* thermalign.cpp
* PHunterLib
*
* Created by Tiayyba Riaz on 7/2/09.
* Copyright 2009 __MyCompanyName__. All rights reserved.
*
*/
//=============================================================================
// Module: thermalign.cpp
// Project: Diploma Thesis - Probe Selection for DNA Microarrays
// 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 CThermAlign - Thermodynamic Alignment Algorithm
// Author: kaderali
// Date: 9/2000 - 01/2001
// Copyright: (c) L. Kaderali, 9/2000 - 01/2001
//
// 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 <math.h>
#include <assert.h>
//#include <iomanip>
#include <stdio.h>
#include <stdlib.h>
#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 (enthalpy1<forbidden_enthalpy)
{
entropy1 += nparam_GetEntropy(thrma->NNParams, prevchari, seq1char, prevcharj, seq2char);
enthalpy1 += nparam_GetEnthalpy(thrma->NNParams, prevchari, seq1char, prevcharj, seq2char);
}
if (enthalpy2<forbidden_enthalpy)
{
enthalpy2 += nparam_GetEnthalpy(thrma->NNParams, prevchari, seq1char, 0, seq2char); // 0 is gap!
entropy2 += nparam_GetEntropy(thrma->NNParams, prevchari, seq1char, 0, seq2char);
}
if (enthalpy3<forbidden_enthalpy)
{
enthalpy3 += nparam_GetEnthalpy(thrma->NNParams, 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 (enthalpy1<forbidden_enthalpy)
{
entropy1 += nparam_GetEntropy(thrma->NNParams, prevchari, seq1char, seq2char, 0);
enthalpy1 += nparam_GetEnthalpy(thrma->NNParams, prevchari, seq1char, seq2char, 0);
}
if (enthalpy2<forbidden_enthalpy)
{
entropy2 += nparam_GetEntropy(thrma->NNParams, 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 (enthalpy1<forbidden_enthalpy)
{
entropy1 += nparam_GetEntropy(thrma->NNParams, seq1char, 0, prevcharj, seq2char);
enthalpy1 += nparam_GetEnthalpy(thrma->NNParams, seq1char, 0, prevcharj, seq2char);
}
if (enthalpy2<forbidden_enthalpy)
{
entropy2 += nparam_GetEntropy(thrma->NNParams, 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 ("<<i<<","<<j<<"):"<<endl;
PrintDPTable(dboutstream);
dboutstream<<endl<<flush;
*/
}
}
return;
}
///////////////////////////////////////////////////////////////////////////////
// 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
// #$
float therm_GetEntropy(PThermAlign thrma, int i, int j)
{
float entropy;
// 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))
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;
}

105
src/libthermo/thermalign.h
View File

@ -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 <iostream>
//#include <strstream>
//#include <ostream>
#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_)