libsequence  1.9.5
ThreeSubs.cc
1 /*
2 
3  Copyright (C) 2003-2009 Kevin Thornton, krthornt[]@[]uci.edu
4 
5  Remove the brackets to email me.
6 
7  This file is part of libsequence.
8 
9  libsequence is free software: you can redistribute it and/or modify
10  it under the terms of the GNU General Public License as published by
11  the Free Software Foundation, either version 3 of the License, or
12  (at your option) any later version.
13 
14  libsequence is distributed in the hope that it will be useful,
15  but WITHOUT ANY WARRANTY; without even the implied warranty of
16  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17  GNU General Public License for more details.
18 
19  You should have received a copy of the GNU General Public License
20  long with libsequence. If not, see <http://www.gnu.org/licenses/>.
21 
22 */
23 
24 #include <cassert>
25 #include <Sequence/Grantham.hpp>
29 #include <Sequence/SingleSub.hpp>
30 #include <Sequence/ThreeSubs.hpp>
31 
32 using std::string;
33 //handle cases where codons are fully substituted
34 namespace Sequence
35 {
36  struct ThreeSubs::ThreeSubsImpl
37  {
38  double p0, p2S, p2V, p4, q0, q2S, q2V, q4;
39  void Calculate (const RedundancyCom95 & sitesObj,
40  const Inter3_t & intermediates,
41  const std::string &codon1, const std::string &codon2,
42  double w_path1,double w_path2, double w_path3,
43  double w_path4,double w_path5,double w_path6);
44  ThreeSubsImpl():p0(0.), p2S(0.), p2V(0.), p4(0.), q0(0.), q2S(0.), q2V(0.), q4(0.)
45  {
46  }
47  };
48 
49  ThreeSubs::ThreeSubs() : impl(std::unique_ptr<ThreeSubsImpl>(new ThreeSubsImpl()))
50  {
51  }
52 
57  {
58  }
59 
60  void ThreeSubs::operator() (const RedundancyCom95 & sitesObj,
61  const std::string &codon1, const std::string &codon2,
62  const Sequence::WeightingScheme3 *weights3)
70  {
71  assert(codon1.length() == 3 && codon2.length() == 3);
72  auto intermediates =Intermediates3(codon1,codon2);
73  impl->p0 = impl->p2S = impl->p2V = impl->p4 = impl->q0 = impl->q2S = impl->q2V = impl->q4 = 0.0;
74 
75  auto weights = weights3->operator()(codon1,codon2,sitesObj.gencode());
76  impl->Calculate (sitesObj, intermediates, codon1, codon2, weights[0],
77  weights[1], weights[2], weights[3], weights[4], weights[5]);
78  }
79 
80  void
81  ThreeSubs::ThreeSubsImpl::Calculate (const RedundancyCom95 & sitesObj,
82  const Inter3_t & intermediates,
83  const std::string & codon1, const std::string & codon2,
84  double w_path1, double w_path2, double w_path3,
85  double w_path4, double w_path5, double w_path6)
89  {
90  //arrays to store transition / transversion values per branch
91  double p0_b[15], p2S_b[15], p2V_b[15], p4_b[15];
92  double q0_b[15], q2S_b[15], q2V_b[15], q4_b[15];
93 
94  //initialuze the SingleSub class for the first path
95  SingleSub Single;
96  //there are 15 branches to iterate over--you'll need to draw a picture...
97  //there is a picture in the documentation
98  for (int i = 0; i <= 14; ++i)
99  {
100  switch (i)
101  {
102  case 0:
103  //already initialized above
104  // *Single = SingleSub (sitesObj, codon1,intermediates[0]);
105  Single(sitesObj, codon1,intermediates[0]);
106  break;
107  case 1:
108  Single (sitesObj, intermediates[0],
109  intermediates[1]);
110  break;
111  case 2:
112  Single (sitesObj, intermediates[1],
113  codon2);
114  break;
115  case 3:
116  Single (sitesObj, intermediates[0],
117  intermediates[2]);
118  break;
119  case 4:
120  Single (sitesObj, intermediates[2],
121  codon2);
122  break;
123  case 5:
124  Single (sitesObj, codon1,
125  intermediates[3]);
126  break;
127  case 6:
128  Single (sitesObj, intermediates[3],
129  intermediates[4]);
130  break;
131  case 7:
132  Single (sitesObj, intermediates[4],
133  codon2);
134  break;
135  case 8:
136  Single (sitesObj, intermediates[3],
137  intermediates[5]);
138  break;
139  case 9:
140  Single (sitesObj, intermediates[5],
141  codon2);
142  break;
143  case 10:
144  Single (sitesObj, codon1,
145  intermediates[6]);
146  break;
147  case 11:
148  Single (sitesObj, intermediates[6],
149  intermediates[7]);
150  break;
151  case 12:
152  Single (sitesObj, intermediates[7],
153  codon2);
154  break;
155  case 13:
156  Single (sitesObj, intermediates[6],
157  intermediates[8]);
158  break;
159  case 14:
160  Single (sitesObj, intermediates[8],
161  codon2);
162  break;
163  }
164  p0_b[i] = Single.P0();
165  p2S_b[i] =Single.P2S();
166  p2V_b[i] =Single.P2V();
167  p4_b[i] = Single.P4();
168  q0_b[i] = Single.Q0();
169  q2S_b[i] =Single.Q2S();
170  q2V_b[i] =Single.Q2V();
171  q4_b[i] = Single.Q4();
172  }
173  //sum up changes along each branch, weighting by the
174  //weight factor for each path
175  p0 = (p0_b[0] + p0_b[1] + p0_b[2]) * w_path1
176  + (p0_b[0] + p0_b[3] + p0_b[4]) * w_path2
177  + (p0_b[5] + p0_b[6] + p0_b[7]) * w_path3
178  + (p0_b[5] + p0_b[8] + p0_b[9]) * w_path4
179  + (p0_b[10] + p0_b[11] + p0_b[12]) * w_path5
180  + (p0_b[10] + p0_b[13] + p0_b[14]) * w_path6;
181 
182  p2S = (p2S_b[0] + p2S_b[1] + p2S_b[2]) * w_path1
183  + (p2S_b[0] + p2S_b[3] + p2S_b[4]) * w_path2
184  + (p2S_b[5] + p2S_b[6] + p2S_b[7]) * w_path3
185  + (p2S_b[5] + p2S_b[8] + p2S_b[9]) * w_path4
186  + (p2S_b[10] + p2S_b[11] + p2S_b[12]) * w_path5
187  + (p2S_b[10] + p2S_b[13] + p2S_b[14]) * w_path6;
188 
189  p2V = (p2V_b[0] + p2V_b[1] + p2V_b[2]) * w_path1
190  + (p2V_b[0] + p2V_b[3] + p2V_b[4]) * w_path2
191  + (p2V_b[5] + p2V_b[6] + p2V_b[7]) * w_path3
192  + (p2V_b[5] + p2V_b[8] + p2V_b[9]) * w_path4
193  + (p2V_b[10] + p2V_b[11] + p2V_b[12]) * w_path5
194  + (p2V_b[10] + p2V_b[13] + p2V_b[14]) * w_path6;
195 
196  p4 = (p4_b[0] + p4_b[1] + p4_b[2]) * w_path1
197  + (p4_b[0] + p4_b[3] + p4_b[4]) * w_path2
198  + (p4_b[5] + p4_b[6] + p4_b[7]) * w_path3
199  + (p4_b[5] + p4_b[8] + p4_b[9]) * w_path4
200  + (p4_b[10] + p4_b[11] + p4_b[12]) * w_path5
201  + (p4_b[10] + p4_b[13] + p4_b[14]) * w_path6;
202 
203  q0 = (q0_b[0] + q0_b[1] + q0_b[2]) * w_path1
204  + (q0_b[0] + q0_b[3] + q0_b[4]) * w_path2
205  + (q0_b[5] + q0_b[6] + q0_b[7]) * w_path3
206  + (q0_b[5] + q0_b[8] + q0_b[9]) * w_path4
207  + (q0_b[10] + q0_b[11] + q0_b[12]) * w_path5
208  + (q0_b[10] + q0_b[13] + q0_b[14]) * w_path6;
209 
210  q2S = (q2S_b[0] + q2S_b[1] + q2S_b[2]) * w_path1
211  + (q2S_b[0] + q2S_b[3] + q2S_b[4]) * w_path2
212  + (q2S_b[5] + q2S_b[6] + q2S_b[7]) * w_path3
213  + (q2S_b[5] + q2S_b[8] + q2S_b[9]) * w_path4
214  + (q2S_b[10] +q2S_b[11] + q2S_b[12]) * w_path5
215  + (q2S_b[10] + q2S_b[13] + q2S_b[14]) * w_path6;
216 
217  q2V = (q2V_b[0] + q2V_b[1] + q2V_b[2]) * w_path1
218  + (q2V_b[0] + q2V_b[3] + q2V_b[4]) * w_path2
219  + (q2V_b[5] + q2V_b[6] + q2V_b[7]) * w_path3
220  + (q2V_b[5] + q2V_b[8] + q2V_b[9]) * w_path4
221  + (q2V_b[10] + q2V_b[11] + q2V_b[12]) * w_path5
222  + (q2V_b[10] + q2V_b[13] + q2V_b[14]) * w_path6;
223 
224  q4 = (q4_b[0] + q4_b[1] + q4_b[2]) * w_path1
225  + (q4_b[0] + q4_b[3] + q4_b[4]) * w_path2
226  + (q4_b[5] + q4_b[6] + q4_b[7]) * w_path3
227  + (q4_b[5] + q4_b[8] + q4_b[9]) * w_path4
228  + (q4_b[10] + q4_b[11] + q4_b[12]) * w_path5
229  + (q4_b[10] + q4_b[13] + q4_b[14]) * w_path6;
230  }
231 
232  double ThreeSubs::P0 (void) const
233  {
234  return impl->p0;
235  }
236 
237  double ThreeSubs::P2S (void) const
238  {
239  return impl->p2S;
240  }
241 
242  double ThreeSubs::P2V (void) const
243  {
244  return impl->p2V;
245  }
246 
247  double ThreeSubs::P4 (void) const
248  {
249  return impl->p4;
250  }
251 
252  double ThreeSubs::Q0 (void) const
253  {
254  return impl->q0;
255  }
256 
257  double ThreeSubs::Q2S (void) const
258  {
259  return impl->q2S;
260  }
261 
262  double ThreeSubs::Q2V (void) const
263  {
264  return impl->q2V;
265  }
266 
267  double ThreeSubs::Q4 (void) const
268  {
269  return impl->q4;
270  }
271 }
272 
double P2S(void) const
Definition: SingleSub.cc:78
double Q0(void) const
Definition: SingleSub.cc:104
double Q2S(void) const
Definition: SingleSub.cc:113
double Q2V(void) const
Definition: SingleSub.cc:122
double P2S(void) const
Definition: ThreeSubs.cc:237
double P4(void) const
Definition: ThreeSubs.cc:247
STL namespace.
The namespace in which this library resides.
Deal with codons differing at 1 position.
Definition: SingleSub.hpp:45
class Sequence::RedundancyCom95
used by Sequence::Comeron95, class Sequence::ThreeSubs calculates divergence between codons that diff...
Grantham&#39;s distances (Sequence::Grantham)
void operator()(const RedundancyCom95 &sitesObj, const std::string &codon1, const std::string &codon2, const Sequence::WeightingScheme3 *weights3)
Definition: ThreeSubs.cc:60
used by Sequence::Comeron95, class Sequence::SingleSub calculates divergence between codons that diff...
abstract interface to weighting schemes when codons differ at 3 positions
declarations of Sequence::Intermediates2 and Sequence::Intermediates3
double Q0(void) const
Definition: ThreeSubs.cc:252
double P0(void) const
Definition: SingleSub.cc:69
double Q2S(void) const
Definition: ThreeSubs.cc:257
double P4(void) const
Definition: SingleSub.cc:96
double Q4(void) const
Definition: SingleSub.cc:131
abstract interface to weighting schemes when codons differ at 2 positions
double P2V(void) const
Definition: SingleSub.cc:87
double Q2V(void) const
Definition: ThreeSubs.cc:262
Inter3_t Intermediates3(const std::string &codon1, const std::string &codon2)
Calculate the intermediate codons between a pair of codons diverged at 3 positions.
Calculate redundancy of a genetic code using Comeron&#39;s counting scheme.
double P2V(void) const
Definition: ThreeSubs.cc:242
double Q4(void) const
Definition: ThreeSubs.cc:267
double P0(void) const
Definition: ThreeSubs.cc:232