Probability profile

www.cardiff.ac.uk/biosi/research/biosoft/

Probability distribution Q

The 16S rRNA gene has an alternating pattern of conserved and hypervariable regions. Predicting where these regions lie, in any given sequence, is an important step towards distinguishing between natural intergene variation and a genuine sequence anomaly, when comparing two or more sequences. What is needed is a probability distribution for the 16S rRNA gene (with one probability allocated to each base position) which maps the position of the various variable regions.

Defining Q

Consider two, randomly selected, error-free sequences S_q and S_s and their resulting global alignment S_qs: let q_i be the probability that alignment position i contains a different base residue in sequence S_q to that in sequence S_s. Consider further Q = {q_i : q₁, q₂, ..., q_n}, a set of these probabilities, one for each position in the alignment. A knowledge of Q allows us to predict where in the alignment S_qs the two sequences should most vary from one another, and conversely where they should be most similar, assuming both sequences are error-free. Thus, if S_q and S_s are found to deviate from each other in a manner outside of that predicted by Q, then we have good evidence that one, or other (or both) of the sequences is in some way anomalous.

Calculating probability distribution Q

Q can be generated for any dataset of homologous gene sequences as follows:

Amass a collection of error-free, full-length gene sequences representing the taxon to be investigated.
Include a reference sequence to serve as a position guide*.
Align the sequences to create a multiple sequence alignment.
At each alignment position that coincides with a base in the reference sequence, determine the frequency of occurrence of the most common residue f, ignoring gap characters (Fig. 1), so generating the frequency distribution F = {f_i : f₁, f₂, ..., f_n}, where n is the total number of residues in the reference sequence.
Convert F to probability distribution P = {p_i: p₁, p₂, ..., p_n}, where p_i = (f_i - 0.25)/0.75^§ and is the probability that reference position i contains the same residue in each sequence of the alignment (that is, i is a conserved position).
Finally, convert P to Q = {q_i : q₁, q₂, ..., q_n}, where q_i = 1 - p_i and is the probability that reference position i does not contain the same residue in each sequence of the alignment (that is, i is a hypervariable position).

Figure 1. Calculating Q values. In this example, residue T occurs most often at reference sequence position i, so the relative frequency of this base (0.6) is used to calculate probability p_i and hence q_i (0.43). Similarly at position j, G is most common with a relative frequency of 0.89, ignoring gaps, hence q_j is 0.15.

Probability distribution Q as used by Pintail algorithm

The Pintail algorithm, as implemented in the Pintail and Mallard programs, uses the probability distribution Q generated from an aligned dataset of 4,383 full-length Bacteria type-strain sequences, as described in Ashelford et al. (2005).

A plot of Q from this data gives a visual representation of the conserved and hypervariable regions Q describes:

Figure 2. Plot of probability distribution Q generated from 4,383 full-length Bacteria type-strain sequences.

Smoothing the data enables a clear picture of the various regions. The peaks and troughs coincide closely with the known positions of the variable regions within the 16S rRNA gene.

Figure 3. Plot of probability distribution Q, after smoothing data by taking means with a sliding window of 50 bases, moving 1 base at a time. Peaks represent hypervariable regions whilst troughs represent conserved regions.

Note...

* In order to relate probability distribution Q to any other sequence alignment (e.g. during running of the Pintail algorithm) it is necessary to include a reference sequence to serve as a suitable position guide. Q probabilities are calculated with reference to this sequence, hence the size of Q corresponds to the length of the reference. The implementation of the Pintail algorithm used by the Pintail and Mallard programs uses a Q distribution generated with Escherichia coli K12 U00096 as reference.

^§ At any particular position i, the values for f_i will range from 0.25 to 1, inclusive; f_i can never be lower than 0.25 since the lowest possible frequency of the most common base occurs when all four bases (A, C, G, and T/U) occur in equal numbers (that is, are equally common). Consequently, to convert frequencies to probabilities (which should range from 0 to 1, inclusive), each frequency is converted with the equation p_i = (f_i - 0.25)/0.75.