c_cpp/etc/KaKs_Calculator/src/NG86.cpp

   1 /*********************************************************
   2 * Copyright (C) 2005, BGI of Chinese Academy of Sciences
   3 * All rights reserved.
   4 *
   5 * Filename: NG86.cpp
   6 * Abstract: Definition of NG86 class.
   7 *
   8 * Version: 1.0
   9 * Author: Zhang Zhang (zhanghzhang@genomics.org.cn)
  10 * Date: Feb.21, 2005
  11
  12   Reference: Nei M, Gojobori T  (1986)  Simple methods for
  13   estimating the numbers of synonymous and nonsynonymous
  14   nucleotide substitutions. Mol Biol Evol 3:418-426.
  15 **********************************************************/
  16
  17 #include "NG86.h"
  18
  19 NG86::NG86() {
  20         name = "NG";
  21 }
  22
  23
  24 //Count synonymous(S) and nonsynonymous(N) sites
  25 void NG86::getCondonSite(string codon) {
  26
  27         int i,j,stop;
  28         string temp="";
  29         double syn = 0.0;
  30
  31         if (getAminoAcid(codon)=='!')
  32                 return;
  33
  34         /* Synonymous sites only occur at first and third position in a codon */
  35         for(i=0, stop=0; i<3; i+=2) {
  36                 for(j=0; j<4; j++) {
  37                         temp = codon;
  38                         if (j!=convertChar(temp[i])) {
  39                                 temp[i] = convertInt(j);
  40                                 if (getAminoAcid(temp)=='!') {
  41                                         stop++;
  42                                 }
  43                                 else {
  44                                         if (getAminoAcid(temp)==getAminoAcid(codon))
  45                                                 syn++;
  46                                 }
  47                         }
  48                 }
  49         }
  50         S += (syn/3.0);
  51         N += (3-stop/3.0-syn/3.0);
  52 }
  53
  54 //Count synonymous(Sd) and nonsynonymous(Nd) differences
  55 void NG86::getCondonDifference(string codon1, string codon2) {
  56
  57         int i,j,k,diff[CODONLENGTH];
  58         int num = 0;
  59         int stop = 0;
  60         int path = 1;
  61         double sd_temp = 0.0;
  62         double nd_temp = 0.0;
  63         string temp1, temp2;
  64
  65
  66         if (getAminoAcid(codon1)=='!' || getAminoAcid(codon2)=='!')
  67                 return;
  68
  69         for (i=0; i<CODONLENGTH; i++) {
  70                 diff[i] = -1;
  71                 if (codon1[i]!=codon2[i])
  72                         diff[num++] = i;
  73         }
  74
  75         //two codons are same
  76         if (num==0) return;
  77
  78         snp += num;
  79
  80         //Pathway of evolution from the differences
  81         for (i=1; i<=num; i++)
  82                 path *= i;
  83
  84         //Only one difference between two codons
  85         if (num==1) {
  86                 if (getAminoAcid(codon1)==getAminoAcid(codon2))
  87                         sd_temp++;
  88                 else
  89                         nd_temp++;
  90         }
  91
  92         //Two differences between two codons
  93         if (num==2) {
  94                 for(i=0;i<num;i++)
  95                         for(j=0;j<num;j++)
  96                                 if(i!=j) {
  97                                         temp1 = codon1;
  98                                         temp1[diff[i]] = codon2[diff[i]];
  99                                         if (getAminoAcid(temp1)!='!') {
 100
 101                                                 //codon1<->temp1
 102                                                 if (getAminoAcid(temp1)==getAminoAcid(codon1)) sd_temp++;
 103                                                 else nd_temp++;
 104
 105                                                 //temp1<->codon2
 106                                                 if (getAminoAcid(temp1)==getAminoAcid(codon2)) sd_temp++;
 107                                                 else nd_temp++;
 108                                         }
 109                                         else
 110                                                 stop++;
 111                                 }
 112         }
 113         //Three differences between two codons
 114         if(num==3) {
 115                 for (i=0;i<3;i++) {
 116                 for (j=0;j<3;j++) {
 117                         for (k=0;k<3;k++) {
 118                                         if ((i!=j) && (i!=k) && (j!=k)) {
 119                                                 temp1 = codon1;
 120                                                 temp1[diff[i]] = codon2[diff[i]];
 121                                                 temp2 = temp1;
 122                                                 temp2[diff[j]] = codon2[diff[j]];
 123                                                 if (getAminoAcid(temp1)!='!' && getAminoAcid(temp2)!='!') {
 124                                                         //codon1<->temp1
 125                                                         if (getAminoAcid(temp1)==getAminoAcid(codon1)) sd_temp++;
 126                                                         else nd_temp++;
 127
 128                                                         //temp1<->temp2
 129                                                         if (getAminoAcid(temp2)==getAminoAcid(temp1)) sd_temp++;
 130                                                         else nd_temp++;
 131
 132                                                         //temp2<->codon2
 133                                                         if (getAminoAcid(codon2)==getAminoAcid(temp2)) sd_temp++;
 134                                                         else nd_temp++;
 135
 136                                                 }
 137                                                 else
 138                                                         stop++;
 139                                         }
 140                                 }
 141                         }
 142                 }
 143         }
 144         if (path==stop) {
 145                 //All pathways are through stop codons
 146                 if (num==2) {
 147                         Sd+=0.5;        Nd+=1.5;
 148                 }
 149                 else {
 150                         Sd+=1.0;        Nd+=2.0;
 151                 }
 152         }
 153         else {
 154                 Sd += (double)(sd_temp/(path-stop));
 155                 Nd += (double)(nd_temp/(path-stop));
 156         }
 157
 158 }
 159
 160 void NG86::PreProcess(string seq1, string seq2) {
 161
 162         long i;
 163
 164         //Count sites and differences
 165         for(i=0; i<seq1.length(); i=i+3) {
 166                 getCondonSite(seq1.substr(i,3));
 167                 getCondonSite(seq2.substr(i,3));
 168                 getCondonDifference(seq1.substr(i,3), seq2.substr(i,3));
 169         }
 170
 171         S/=2.0;
 172         N/=2.0;
 173
 174         //Scale the sum of S+N to the length of sequence.
 175         double y=seq1.length()/(S+N);
 176         S*=y;
 177         N*=y;
 178 }
 179
 180 /* Jukes & Cantor's one-parameter formula for correction */
 181 double NG86::kaks_formula(double p) {
 182         //Equation (3) in the reference of NG86
 183         double d = 1-(4*p)/3;
 184         if (d<0.0) {
 185                 d = NA;
 186         }
 187         else {
 188                 d = log(d);
 189                 if (d>0.0)
 190                         d = NA;
 191                 else
 192                         d = (-3.0)*d/4.0;
 193         }
 194
 195         return d;
 196 }
 197
 198 string NG86::Run(string seq1, string seq2) {
 199
 200         PreProcess(seq1, seq2);
 201
 202         Ks = kaks_formula(Sd/S);
 203         Ka = kaks_formula(Nd/N);
 204
 205         t = (S*Ks+N*Ka)/(S+N);
 206
 207         return  parseOutput();
 208 }
 209
 210
 211
 212 /***********************************************************
 213   NONE: an in-house algorithm in BGI for testing Ka and Ks.
 214   NONE is NG86 without correction for multiple substitutions.
 215 ************************************************************/
 216 NONE::NONE() {
 217         name = "NONE";
 218 }
 219
 220 string NONE::Run(string seq1, string seq2) {
 221
 222         PreProcess(seq1, seq2);
 223
 224         Ks = Sd/S;
 225         Ka = Nd/N;
 226
 227         t = (S*Ks+N*Ka)/(S+N);
 228
 229         return  parseOutput();
 230 }