Calculating dNdS for NGS datasets

• Post by: Richard Orton
• November 16, 2014
• 22 Comments

vNvS

Our upcoming tool vNvS calculates the dN/dS ratio at each site, codon and also for the sample as a whole, here is an explanation of the theory behind it. vNvS is currently in development – for more information email Richard.Orton@glasgow.ac.uk

dN/dS

dN/dS is the ratio of the number of nonsynonymous substitutions per non-synonymous site (p_N) to the number of synonymous substitutions per synonymous site (p_S), which can be used as an indicator of selective pressure acting on a protein coding gene. A ratio greater than one implies positive or darwinian selection, less than one implies purifying selection, and a ratio of one indicates neutral selection.

A synonymous substitution is a nucleotide mutation that does not alter the amino acid (AA) of a protein; also called a silent mutation. However, synonymous substitutions can still potentially affect transcription, splicing, mRNA transport or translation, any of which could alter the phenotype and render the synonymous mutation non-silent.

A nonsynonymous substitution is a nucleotide mutation that alters the AA of a protein; also called a missense mutation. However, nonsynonymous substitutions are called nonsense mutations if the codon changes into a (premature) stop codon (i.e. TAA, TAG, TGA).

We will now consider two example AA codons and evaluate the types of mutations for each: Tryptophan which is mutationally fragile and Arginine which is mutationally robust. A codon consists of three nucleotides (the so called triplet code), and as there are four possible nucleotides (A, C, G, T), there are a total of 9 possible mutations to consider for each codon.

[table id=1 /]

Any mutation at any of the three codon positions for Tryptophan results in a nonsynonymous mutation. This emphasises why the expected number of nonsynonymous mutations  – or the number of nonsynonymous sites – should be taken into account in dN/dS analyses, as observed counts alone would lead to biased results. We will come to the N-Sites column later.

[table id=2 /]

The Arginine codon CGG has 4 possible synonymous mutations and 5 nonsynonymous mutations. It can withstand any mutation at the 3^rd codon position and still code for Arginine and also withstand one substitution at the 1^st codon position. However, any mutation at the 2^nd codon position is nonsynonymous as it results in an AA substitution.

This is true for all AA codons – any mutation at the 2^nd codon position will result in an AA substitution, whilst mutations at the 3^rd codon position either cause a silent mutation or a nonsynonymous mutation that typically maintains the hydrophilicity/hydrophobicity of the AA. This is one reason why the 3^rd codon position is most liable to mutation.

Expected Sites for a Reference Sequence

Central to the calculation of the dN/dS ratio is calculating the expected number of nonsynonymous sites for each codon – this is why the dN/dS ratio is described in part as “nonsynonymous substitutions per non-synonymous site”.

For any particular codon, we denote f_i to be the fraction of  substitutions that are nonsynonymous at the i^th codon position (i = 1, 2, 3). The number of nonsynonymous (n) and synonymous (s) sites for a codon is then calculated with the formulas:

Using our previous examples, and the values already populated in the N-Sites column of the tables:

• TGG – Tryptophan – Trp – W
• n = (1/3 + 1/3 + 1/3) + (1/3 + 1/3 + 1/3) + (1/3 + 1/3 + 1/3) = 3
• s = 3 – n = 3 – 3 = 0
• The number of nonsynonymous sites for Tryptophan TGG is 3, whilst the number of synonymous sites is 0.

• CGG – Arginine – Arg – R
• n = (0 + 1/3 + 1/3) + (1/3 + 1/3 + 1/3) + (0 + 0 + 0) = 1.66666
• s = 3 – n = 3 – 1.66666 = 1.33333
• The number of nonsynonymous sites for Arginine CGG is 1.66666, whilst the number of synonymous sites is 1.33333.

For a protein encoding DNA sequence of length r codons, the total number of nonsynonymous (N) and synonymous (S) sites for the whole sequence can therefore be calculated by summing up the individual n and s values from all codons using the formulas:

As an example, consider the hypothetical reference sequence ATG – AAA – CCC – GGG – TTT – TAA. Given this sequence, N and S are calculated as:

[table id=3 /]

Observed Substitutions for a Sample Sequence

Now consider that we have a new sequence, our sample sequence, and we want to compare it to our original reference sequence, so we must count the number of observed nonsynonymous (n_d) and synonymous (s_d) mutations that our sample sequence has at each codon compared to the reference.

[table id=4 /]

In the above, there are three codons that have mutations:

Codon 3: CCC (pro) -> CGC (arg): This has a single nucleotide mutation (distance = 1) resulting in a nonsynonymous AA substitution, and therefore n_d = 1 and s_d = 0 for this codon.

Codon 4: GGG (gly) -> GGC (gly): This has a single nucleotide mutation (distance = 1) resulting in a synonymous AA substitution, and therefore n_d = 0 and s_d = 1 for this codon.

Codon 5: TTT (phe) -> TAC (tyr): This has two nucleotide mutations (distance = 2) and is a nonsynonymous AA substitution. However, due to the double mutation it is a more complex situation, as one must consider the two mutation pathways that could have led to this state (you can not assume that as the AA is nonsynonymous then both mutations are nonsynonymous):

1. TTT (phe) –> TAT (tyr) -> TAC (tyr): 1 nonsynonymous and 1 synonymous mutation
2. TTT (phe) –> TTC (phe) -> TAC (tyr): 1 synonymous and 1 nonsynonymous mutation

• As we consider the two above pathways to occur with equal probability, for this codon n_d = 1 and s_d = 1 (sums to two as we have two mutations in the codon).

Hypothetically, let us expand the example above to a three mutation (distance = 3) situation from TTT (phe) to GAC (asp). A three mutation situation results in 6 possible mutation pathways:

1. TTT (phe) -> TTC (phe) -> TAC (tyr) -> GAC (asp): 2n and 1s
2. TTT (phe) -> TTC (phe) -> GTC (val) -> GAC (asp): 2n and 1s
3. TTT (phe) -> TAT (tyr) -> TAC (val) -> GAC (asp): 3n
4. TTT (phe) -> TAT (tyr) -> GAT (val) -> GAC (asp): 3n
5. TTT (phe) -> GTT (val) -> GTC (val) -> GAC (asp): 2n and 1s
6. TTT (phe) -> GTT (val) -> GAT (val) -> GAC (asp): 2n and 1s

• As we consider the six above pathways to occur with equal probability, for this codon n_d = 2.333 and s_d = 0.666 (sums to three as we have three mutations in the codon).

So, for a protein encoding DNA sequence of length r codons, the total number of observed nonsynonymous (N_d) and synonymous (S_d) mutations between two sequences can therefore be calculated by summing up the individual n_d and s_d values from all codons using the formulas:

dN/dS ratio

One then calculates the proportion of nonsynonymous (p_n) and synonymous (p_s) differences with the following equations:

Then to estimate the number of nonsynonymous substitutions (d_N) and synonymous substitutions (d_S) per site, and the dN/dS ratio itself, we use the formulas:

So using our reference and sample sequence examples above we now have the following values:

• N = 14.666
• S = 3.333
• N_d = 2
• S_d = 2
• p_N = 0.1364
• p_S = 0.6001
• d_N = 0.1505
• d_S = 1.2074
• dN/dS = 0.1247

NGS Datasets

The expansion of all of the above to NGS datasets is relatively straightforward. The calculation of the number of nonsynonymous and synonymous sites in the reference sequences proceeds the same as above.

One must then calculate the number of observed nonsynonymous and synonymous mutations in the reads when compared to the reference (this will need information on where the open reading frame [ORF] starts and stops). There are two options when calculating the observed numbers:

1. Consider all observed mutations in the reads covering any part of the codon.
2. Consider only those mutations where the read fully covers the codon the mutation occurs in, i.e. ignore partially covered codons at the read ends, or partially covered codons due to indels.

To adapt to NGS datasets, one must consider read coverage. The approach taken by Morelli et al (2013) adjusts the formula for p_N (and likewise p_S) as follows to take into account the read coverage (c) at each codon:

 

Essentially, for each read (c) that covers a particular codon, the observed number of nonsynonymous mutations in the read compared to the reference codon is calculated, and divided by the expected number. The values for all reads at the codon are then summed and averaged. Then the value for all codons is summed to give a single value for the whole ORF. This value of p_N and p_S can they be plugged in to the same dN and dS formulas to calculate the dN/dS ratio.

References

Simple methods for estimating the numbers of synonymous and nonsynonymous nucleotide substitutions

Evolution of foot-and-mouth disease virus intra-sample sequence diversity during serial transmission in bovine hosts