Ka and Ks, also known as dn and ds, computation
Ks and Ka are, respectively, the number of substitutions per synonymous site and per non-synonymous site between two protein-coding genes. They are also denoted as ds and dn in the literature. The ratio of nonsynonymous (Ka) to synonymous (Ks) nucleotide substitution rates is an indicator of selective pressures on genes. A ratio significantly greater than 1 indicates positive selective pressure. A ratio around 1 indicates either neutral evolution at the protein level or an averaging of sites under positive and negative selective pressures. A ratio less than 1 indicates pressures to conserve protein sequence (i.e. purifying selection). This function estimates the Ka and Ks values for a set of aligned sequences using the method published by Li (1993) and gives the associated variance matrix.
kaks(x, verbose = FALSE, debug = FALSE, forceUpperCase = TRUE, rmgap = TRUE)
x |
An object of class |
verbose |
If TRUE add to the results the value of L0, L2, L4 (respectively the frequency of non-synonymous sites, of 2-fold synonymous sites, of 4-fold synonymous sites), A0, A2, A4 (respectively the number of transitional changes at non-synonymous, 2-fold, and 4-fold synonymous sites ) and B0, B2, B4 (respectively the number of transversional changes at non-synonymous, 2-fold, and 4-fold synonymous sites). |
debug |
If TRUE turns debug mode on. |
forceUpperCase |
If TRUE, the default value, all character in sequences are forced to the upper case if at least one 'a', 'c', 'g', or 't' is found in the sequences. Turning it to FALSE if the sequences are already in upper case will save time. |
rmgap |
If TRUE all positions with at least one gap are removed. If FALSE only positions with nothing else than gaps are removed. |
ks |
matrix of Ks values |
ka |
matrix of Ka values |
vks |
variance matrix of Ks |
vka |
variance matrix of Ka |
Computing Ka and Ks makes sense for coding sequences that have been aligned at the amino-acid level before retro-translating the alignement at the nucleic acid level to ensure that sequences are compared on a codon-by-codon basis. Function reverse.align may help for this.
As from seqinR 2.0-3, when there is at least one non ACGT base in a codon, this codon is considered as a gap-codon (---). This makes the computation more robust with respect to alignments with out-of-frame gaps, see example section.
Gap-codons (---) are not used for computations.
When the alignment does not contain enough information (i.e. close to saturation), the Ka and Ks values are forced to 10 (more exactly to 9.999999).
Negative values indicate that Ka and Ks can not be computed.
According to Li (1993) and Pamilo and Bianchi (1993), the rate of synonymous substitutions Ks is computed as: Ks = (L2.A2 + L4.A4) / (L2 + L4) + B4
and the rate of non-synonymous substitutions Ka is computed as: Ka = A0 + (L0.B0 + L2.B2) / (L0 + L2)
D. Charif, J.R. Lobry
Li, W.-H., Wu, C.-I., Luo, C.-C. (1985) A new method for estimating synonymous and nonsynonymous rates of nucleotide substitution considering the relative likelihood of nucleotide and codon changes. Mol. Biol. Evol, 2:150-174
Li, W.-H. (1993) Unbiased estimation of the rates of synonymous and nonsynonymous substitution. J. Mol. Evol., 36:96-99.
Pamilo, P., Bianchi, N.O. (1993) Evolution of the Zfx and Zfy genes: Rates and interdependence between genes. Mol. Biol. Evol, 10:271-281
Hurst, L.D. (2002) The Ka/Ks ratio: diagnosing the form of sequence evolution.
Trends Genet., 18:486-486.
The C programm implementing this method was provided by Manolo Gouy. More info is
needed here to trace back the original C source so as to credit correct source.
The original FORTRAN-77 code by Chung-I Wu modified by Ken Wolfe was available
here http://wolfe.gen.tcd.ie/lab/pub/li93/ but this is no more true as 2017-07-01.
For a more recent discussion about the estimation of Ka and Ks see:
Tzeng, Y.H., Pan, R., Li, W.-H. (2004) Comparison of three methods for estimating
rates of synonymous and nonsynonymous nucleotide substitutions.
Mol. Biol. Evol, 21:2290-2298.
The method implemented here is noted LWL85 in the above paper.
The cite this package in a publication, as any R package, try something as citation("seqinr")
at your R prompt.
read.alignment to import alignments from files, reverse.align to align CDS at the aa level,
kaksTorture for test on one-codon CDS.
#
# Simple Toy example:
#
s <- read.alignment(file = system.file("sequences/test.phylip", package = "seqinr"),
format = "phylip")
kaks(s)
#
# Check numeric results on an simple test example:
#
data(AnoukResult)
Anouk <- read.alignment(file = system.file("sequences/Anouk.fasta", package = "seqinr"),
format = "fasta")
if( ! all.equal(kaks(Anouk), AnoukResult) ) {
warning("Poor numeric results with respect to AnoukResult standard")
} else {
print("Results are consistent with AnoukResult standard")
}
#
# As from seqinR 2.0-3 the following alignment with out-of-frame gaps
# should return a zero Ka value.
#
# >Reference
# ATGTGGTCGAGATATCGAAAGCTAGGGATATCGATTATATATAGCAAGATCGATAGAGGA
# TCGATGATCGATCGGGATCGACAGCTG
# >With out-of-frame gaps
# AT-TGGTCCAGGTATCGTAAGCTAGGGATATCGATTATATATAGCAAGATCGATAGGGGA
# TCGATGATCGATCGGGA--GACAGCTG
#
# This test example provided by Darren Obbard is now used as a routine check:
#
Darren <- read.alignment(file = system.file("sequences/DarrenObbard.fasta", package = "seqinr"),
format = "fasta")
stopifnot( all.equal(kaks(Darren)$ka[1], 0) )
#
# As from seqinR 3.4-0, non-finite values should never be returned for
# Ka and Ks even for small sequences. The following test checks that this
# is true for an alignement of the 64 codons, so that we compute Ka and
# Ks for all possible pairs of codons.
#
wrd <- as.alignment(nb = 64, nam = words(), seq = words())
res <- kaks(wrd)
if(any(!is.finite(res$ka))) stop("Non finite value returned for Ka")
if(any(!is.finite(res$ks))) stop("Non finite value returned for Ks")Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.