c_cpp/etc/KaKs_Calculator/src/LPB93.cpp

   1 /*********************************************************
   2 * Copyright (C) 2005, BGI of Chinese Academy of Sciences
   3 * All rights reserved.
   4 *
   5 * Filename: LPB93.cpp
   6 * Abstract: Definition of LPB93 and MLPB93 class.
   7 *
   8 * Version: 1.0
   9 * Author: Zhang Zhang (zhanghzhang@genomics.org.cn)
  10 * Date: Jan.21, 2005
  11
  12   References:
  13   Li WH  (1993)  Unbiased estimation of the Rates of synonymous
  14   and nonsynonymous substitution. J. Mol. Evol. 36:96-99.
  15
  16   Pamilo P, Bianchi NO  (1993)  Evolution of the Zfx and Zfy
  17   genes: rates and interdependence between the genes. Mol. Biol.
  18   Evol. 10:271-281.
  19
  20   Tzeng Y-H, Pan R, Li W-H  (2004)  Comparison of Three Methods
  21   for Estimating Rates of Synonymous and Nonsynonymous Nucleotide
  22   Substitutions. Mol. Biol. Evol. 21:2290-2298.
  23 **********************************************************/
  24
  25
  26 #include "LPB93.h"
  27
  28 LPB93::LPB93() {
  29         name = "LPB";
  30 }
  31
  32 /* Similar to LWL85 except the formulas for calculating ka and ks*/
  33 string LPB93::Run(string seq1, string seq2) {
  34
  35         preProcess(seq1, seq2);
  36
  37         Ks = B[4] + (L[2]*A[2] + L[4]*A[4])/(L[2] + L[4]);
  38         Ka = A[0] + (L[0]*B[0] + L[2]*B[2])/(L[0] + L[2]);
  39
  40         Sd = L[2]*A[2] + L[4]*K[4];
  41         Nd = L[2]*B[2] + L[0]*K[0];
  42
  43         S = Sd/Ks;
  44         N = Nd/Ka;
  45
  46         t = (L[0]*K[0]+L[2]*K[2]+L[4]*K[4])/(L[0]+L[2]+L[4]);
  47
  48         return parseOutput();
  49 }
  50
  51
  52 /*The difference between LPB93 and MLPB93 focuses on the definition of transition & transversion*/
  53 MLPB93::MLPB93() {
  54         name = "MLPB";
  55 }
  56
  57 /*For more detail see reference: Tzeng Y-H, Pan R, Li W-H  (2004)  Mol. Biol. Evol.*/
  58 int MLPB93::TransitionTransversion(string codon1, string codon2, int pos) {
  59
  60         //1:synonymous, 0:nonsynonymous, -1:uncalculate
  61         int isSyn=-1;
  62
  63         //Ile: ATT, ATC, ATA; Met: ATA
  64         if ((codon1=="ATA" && codon2=="ATG" && pos==2)||(codon1=="ATG" && codon2=="ATA" && pos==2)) {
  65                 isSyn = 0;
  66         }
  67         if ((codon1=="ATA" && (codon2=="ATC"||codon2=="ATT") && pos==2) || ( (codon1=="ATC"||codon1=="ATT") && codon2=="ATA" && pos==2)) {
  68                 isSyn = 1;
  69         }
  70
  71         //Arg: CGA, CGG, AGA, AGG
  72         if ((codon1=="CGA" && codon2=="AGA" && pos==0)||(codon1=="AGA" && codon2=="CGA" && pos==0)) {
  73                 isSyn = 1;
  74         }
  75         if ((codon1=="CGG" && codon2=="AGG" && pos==0)||(codon1=="AGG" && codon2=="CGG" && pos==0)) {
  76                 isSyn = 1;
  77         }
  78
  79         //Synonymous: A<->G, C<->T
  80         //Normal situation
  81         if (isSyn==-1) {
  82                 int sum = convertChar(codon1[pos]) + convertChar(codon2[pos]);
  83                 if (sum==5 || sum==1)
  84                         isSyn = 1;
  85                 else
  86                         isSyn = 0;
  87         }
  88
  89         int class1=getCodonClass(codon1,pos);
  90         int class2=getCodonClass(codon2,pos);
  91         if (isSyn==1) {
  92                 Si_temp[class1] += 0.5;
  93                 Si_temp[class2] += 0.5;
  94         }
  95         if (isSyn==0) {
  96                 Vi_temp[class1] += 0.5;
  97                 Vi_temp[class2] += 0.5;
  98         }
  99
 100         return 0;
 101 }
 102
 103
 104
 105