/*****************************************************************************
# Copyright (C) 1994-2008 by David Gordon.
# All rights reserved.
#
# This software is part of a beta-test version of the Consed/Autofinish
# package. It should not be redistributed or
# used for any commercial purpose, including commercially funded
# sequencing, without written permission from the author and the
# University of Washington.
#
# This software is provided ``AS IS'' and any express or implied
# warranties, including, but not limited to, the implied warranties of
# merchantability and fitness for a particular purpose, are disclaimed.
# In no event shall the authors or the University of Washington be
# liable for any direct, indirect, incidental, special, exemplary, or
# consequential damages (including, but not limited to, procurement of
# substitute goods or services; loss of use, data, or profits; or
# business interruption) however caused and on any theory of liability,
# whether in contract, strict liability, or tort (including negligence
# or otherwise) arising in any way out of the use of this software, even
# if advised of the possibility of such damage.
#
# Building Consed from source is error prone and not simple which is
# why I provide executables. Due to time limitations I cannot
# provide any assistance in building Consed. Even if you do not
# modify the source, you may introduce errors due to using a
# different version of the compiler, a different version of motif,
# different versions of other libraries than I used, etc. For this
# reason, if you discover Consed bugs, I can only offer help with
# those bugs if you first reproduce those bugs with an executable
# provided by me--not an executable you have built.
#
# Modifying Consed is also difficult. Although Consed is modular,
# some modules are used by many other modules. Thus making a change
# in one place can have unforeseen effects on many other features.
# It may takes months for you to notice these other side-effects
# which may not seen connected at all. It is not feasable for me to
# provide help with modifying Consed sources because of the
# potentially huge amount of time involved.
#
#*****************************************************************************/
#include "kimuraBranchProb.h"
#include "transitionsTransversions.h"
#include <math.h>
#include "mbt_exception.h"
#include "rwcstring.h"
#include "assert.h"
kimuraBranchProb :: kimuraBranchProb() {
setTransitionsTransversionsTableIfNecessary();
float fBranchLengthMMul_ex = 0.03736;
float fBranchLengthPPyg_ex = 0.01183;
float fBranchLengthPPyg_in = 0.00638;
float fBranchLengthGGor_ex = 0.00592;
float fBranchLengthGGor_in = 0.00132;
float fBranchLengthHSap_ex = 0.00463;
float fBranchLengthHSap_in = 0.00330;
float fBranchLengthPTro_ex = 0.00130;
float fBranchLengthPPan_ex = 0.00132;
int n;
for( n = 0; n <= 2; ++n ) {
fMMul_ex_[n] = fGetKimuraProb( n, fBranchLengthMMul_ex );
fPPyg_ex_[n] = fGetKimuraProb( n, fBranchLengthPPyg_ex );
fPPyg_in_[n] = fGetKimuraProb( n, fBranchLengthPPyg_in );
fGGor_ex_[n] = fGetKimuraProb( n, fBranchLengthGGor_ex );
fGGor_in_[n] = fGetKimuraProb( n, fBranchLengthGGor_in );
fHSap_ex_[n] = fGetKimuraProb( n, fBranchLengthHSap_ex );
fHSap_in_[n] = fGetKimuraProb( n, fBranchLengthHSap_in );
fPTro_ex_[n] = fGetKimuraProb( n, fBranchLengthPTro_ex );
fPPan_ex_[n] = fGetKimuraProb( n, fBranchLengthPPan_ex );
}
for( n = 0; n <= 255; ++n ) {
n012_[n] = -666;
}
n012_[ 'n' ] = 0;
n012_[ 'i' ] = 1;
n012_[ 'v' ] = 2;
}
float kimuraBranchProb :: fGetKimuraProb(
const int n,
// 0 == no mutation
// 1 == transition
// 2 == transversion
const float fBranchLength ) {
if ( n == 0 ) {
return( 1.0 - fGetKimuraTransitionProb( fBranchLength ) -
fGetKimuraTransversionProb( fBranchLength ) );
}
else if ( n == 1 ) {
return( fGetKimuraTransitionProb( fBranchLength ) );
}
else if ( n == 2 ) {
return( fGetKimuraTransversionProb( fBranchLength ) );
}
else {
assert( false ); // program bug
}
}
float kimuraBranchProb :: fGetKimuraTransitionProb(
const float fBranchLength ) {
float fR = .5 * fTransitionTransversionRatio;
float fProb = 0.25 - 0.5*exp( - ((2.0*fR+1.0)/(fR+1.0))*fBranchLength) +
0.25*exp( -2.0*fBranchLength/(fR+1.0));
return( fProb );
}
float kimuraBranchProb :: fGetKimuraTransversionProb(
const float fBranchLength ) {
float fR = .5 * fTransitionTransversionRatio;
float fProb = 0.5 - 0.5 * exp( -2.0*fBranchLength/(fR+1.0) );
return( fProb );
}
float kimuraBranchProb :: fGetBranchProbability(
const int nIndexOfInitialBase,
const int nIndexOfNewBase,
const char cInitialBase,
const char cNewBase ) {
// is this a transition, a transversion, or no mutation?
// no mutation = 0
// transition = 1
// transversion = 2
char cTransitionOrTransversions =
cTransitionsTransversions[cInitialBase][cNewBase ];
int nTransitionOrTransversion =
n012_[ cTransitionOrTransversions ];
// which branch are we talking about?
RWCString soBranchName;
float fProb;
switch ( nIndexOfInitialBase ) {
case 0:
// branches are MMul_ex, PPyg_ex, and PPyg_in
switch (nIndexOfNewBase ) {
case 9:
soBranchName = "MMul_ex";
fProb = fMMul_ex_[ nTransitionOrTransversion ];
break;
case 8:
soBranchName = "PPyg_ex";
fProb = fPPyg_ex_[ nTransitionOrTransversion ];
break;
case 1:
soBranchName = "PPyg_in";
fProb = fPPyg_in_[ nTransitionOrTransversion ];
break;
default:
RWCString soError = "program error in fGetBranchProbability " +
RWCString( (long) nIndexOfInitialBase ) + " " +
RWCString( (long) nIndexOfNewBase );
THROW_ERROR( soError );
}
break;
case 1:
switch( nIndexOfNewBase ) {
case 7:
soBranchName = "GGor_ex";
fProb = fGGor_ex_[ nTransitionOrTransversion ];
break;
case 2:
soBranchName = "GGor_in";
fProb = fGGor_in_[ nTransitionOrTransversion ];
break;
default:
RWCString soError = "program error in fGetBranchProbability " +
RWCString( (long) nIndexOfInitialBase ) + " " +
RWCString( (long) nIndexOfNewBase );
THROW_ERROR( soError );
}
break;
case 2:
switch( nIndexOfNewBase ) {
case 6:
soBranchName = "HSap_ex";
fProb = fHSap_ex_[ nTransitionOrTransversion ];
break;
case 3:
soBranchName = "HSap_in";
fProb = fHSap_in_[ nTransitionOrTransversion ];
break;
default:
RWCString soError = "program error in fGetBranchProbability " +
RWCString( (long) nIndexOfInitialBase ) + " " +
RWCString( (long) nIndexOfNewBase );
THROW_ERROR( soError );
}
break;
case 3:
switch( nIndexOfNewBase ) {
case 5:
soBranchName = "PTro_ex";
fProb = fPTro_ex_[ nTransitionOrTransversion ];
break;
case 4:
soBranchName = "PPan_ex";
fProb = fPPan_ex_[ nTransitionOrTransversion ];
break;
default:
RWCString soError = "program error in fGetBranchProbability " +
RWCString( (long) nIndexOfInitialBase ) + " " +
RWCString( (long) nIndexOfNewBase );
THROW_ERROR( soError );
}
break;
}
return( fProb );
}