Your name:

Your email address:

MrBayes by example: Identification of sites under positive selection in a protein

Professor Walter M. Fitch and assistant research biologist Robin M. Bush of UCI's Department of Ecology and Evolutionary Biology, working with researchers at the Centers for Disease Control and Prevention, studied the evolution of a prevalent form of the influenza A virus during an 11-year period from 1986 to 1997. They discovered that viruses having mutations in certain parts of an important viral surface protein were more likely than other strains to spawn future influenza lineages. Human susceptibility to infection depends on immunity gained during past bouts of influenza; thus, new viral mutations are required for new epidemics to occur. Knowing which currently circulating mutant strains are more likely to have successful offspring potentially may help in vaccine strain selection. The researchers' findings appear in the Dec. 3 issue of Science magazine.

Fitch and his fellow researchers followed the evolutionary pattern of the influenza virus, one that involves a never-ending battle between the virus and its host. The human body fights the invading virus by making antibodies against it. The antibodies recognize the shape of proteins on the viral surface. Previous infections only prepare the body to fight viruses with recognizable shapes. Thus, only those viruses that have undergone mutations that change their shape can cause disease. Over time, new strains of the virus continually emerge, spread and produce offspring lineages that undergo further mutations. This process is called antigenic drift. "The cycle goes on and on-new antibodies, new mutants," Fitch said.

The research into the virus' genetic data focused on the evolution of the hemagglutinin gene-the gene that codes for the major influenza surface protein. Fitch and fellow researchers constructed "family trees" for viral strains from 11 consecutive flu seasons. Each branch on the tree represents a new mutant strain of the virus. They found that the viral strains undergoing the greatest number of amino acid changes in specified positions of the hemagglutinin gene were most closely related to future influenza lineages in nine of the 11 flu seasons tested.

By studying the family trees of various flu strains, Fitch said, researchers can attempt to predict the evolution of an influenza virus and thus potentially aid in the development of more effective influenza vaccines.

The research team is currently expanding its work to include all three groups of circulating influenza viruses, hoping that contrasting their evolutionary strategies may lend more insight into the evolution of influenza.

Along with Fitch and Bush, Catherine A. Bender, Kanta Subbarao and Nancy J. Cox of the Centers for Disease Control and Prevention participated in the study.

A talk by Walter Fitch (slides and sound)
is here

The goal of this exercise is to detect sites in hemagglutinin that are under positive selection.

Since the analysis takes a very long time to run (several days), here
are the saved results of the MrBayes run: Fitch_HA.nex.p.txt,
Fitch_HA.nex.t.txt .

The original data file is flu_data.paup
. The dataset is obtained from an
article by Yang et al, 2000 . The File used for MrBayes is here

The MrBayes block used to obtain results above is:

set autoclose=yes;

lset nst=2 rates=gamma nucmodel=codon omegavar=Ny98;

mcmcp samplefreq=500 printfreq=500;

mcmc ngen=500000;

sump burnin=50;

sumt burnin=50; end;

Selecting a nucmodel=codon with Omegavar=Ny98 specifies a model in which for every codon the ratio of the rate of non-synonymous to synonymous substitutions is considered. This ratio is called OMEGA. The Ny98 model considers three different omegas, one equal to 1 (no selection, this site is neutral); the second with omega < 1, these sites are under purifying selection; and the third with Omega >1, i.e. these sites are under positive or diversifying selection. (The problem of this model is that only three distinct omegas are estimated, and for each site the probability to fall into one of these three classes. If the omega>1 is estimated to be very large, because one site has a large omega, the other sites might not have a high probability to have the same omega, even though they might also be under positive selection. This leads to the site with largest omega to be identified with confidence, the others have more moderate probabilities to be under positive selection).

Note : Version 2.0 of Mr Bayes has a model that estimates omega for each site individually, the new version only allows the Ny98 model as described above..

- First, you need to detect how many generations to burn in (meaning
the number of samples you will have to discard). Open the file Fitch_HA.nex.p.txt
with Excel and plot # of generations versus -LnL values. Determine
after how many generations the graph becomes "stationary" (hint:
change the Y-axis bounds to "zoom in", e.g., -3300 min to -3200 max).
The burnin value is that number of generations divided by 50 (since
only every 50th generation was sampled; i.e. the burnin value roughly
is equal to the number of rows - not quite because there is a header).
To more accurately determine the burnin, you need to rescale the
Y-axis (click at the Y-axis -- if you aim accurately, you'll get a box
that allows re-scaling).

The result (scatterplot of LogL versus generation) might look like this:

(You also could load the file into tracer, but I am not sure how well this can work to find the site under positive selection.)

- This file contains information for posterior probabilities for each
codon (columns) at each sampled generation (rows). Scroll to the right
to see these columns, starting with pr+(1,2,3), pr+(4,5,6), etc.
Calculate average posterior probability for each site of being under
positive selection (Do not forget to exclude
first N rows as a burnin; you should have detected value of N in the
first question of this exercise - to be clear on where the burnin
ends, you might want to highlight the
rows representing the burnin and select a different font color.
(Use the AVERAGE() function of Excel, enter the formula in a cell
below the values for the individual generations -- starting in column
pr+(1,2,3) -- copy the formula to all columns) (see slides from class
19)

- Plot average posterior probability
*vs.*site #. (select the row in which you calculated the averages, then click Graph, and select a bar graph). Write down the codon positions for a few sites with the highest posterior probability of being positively selected (the columns name pr+(1,2,3), pr+(4,5,6)....and so on. pr+(1,2,3) mean probability of codon #1 (nucleotide #1, #2 and #3) to be under positive selection))

- Determine the 95% credibility interval for the omega<1 value. To
do this you sort posterior probability column in ascending order
(Select data you want to analyze, copy them into a new sheet, then go
to Data->Sort... ). Again, do not forget to ignore the burnin ; the
easiest might be to actually delete it... After sorting, exclude 5% of
the data on the top and on the bottom. The range of the remaining data
gives you the 90% confidence interval. (Enter answer in box below!)

- The structure of hemagglutinin has been crystallized and is publicly
available through PDB. Examin the 2VIU.pdb file (here)
in chimera. Chain A of the PDB file corresponds to the
sequences we did our analysis with (color the molecule according to
chain). Below is a comparison of one of the sequences we used for
analyses with the sequence for which the structure was determined:

Using this alignment as a guide, map the site(s) which have the highest probability to belong to the class with omega>1. Where are these sites located in the protein? (Reminders: The position number in the nexus file corresponds to nucleotide sequence, the structure is based on the amino acid sequence - take the third codon position and divide by 3 to find the amino acid. You only want to be concerned with Chain A!)

What is the 95% crediibility interval for the omega smaller than 1? From your experience does this value indicate strong purifying selection? Which codon(s) showed signs of positive selection? Which position and which amino acid does this correspond to in the above alignment? Where is this aa located in the structure?

dN/dS ratios along a sequence. I ran the dataset from last weeks exercise (both the intein and the extein sequences) in MrBayes using the NY98 model (no partitions). The file is here. The model has a parameter for transition/transversion ratio and for the dN/dS ratio (called omega). The model uses and explores only two of these - one for sites under purifying selection, one for sites under diversifying selection. In addition for each codon the probability that the codon is in the omega+ group, and the omega value for each codon is estimated. This is the MrBayes block in the file:

begin mrbayes;

lset nst=2 rates=gamma nucmodel=codon omegavar=Ny98;

report possel = yes siteomega = yes;

mcmcp filename=analysisS;

mcmcp samplefreq=50 printfreq=50 diagnfreq=500;

mcmcp savebrlens=yes;

end;

This results in a lot of parameters, which
makes for a slow moving robot. The parameter files resulting form a 24h
run is here (already imported
into excel - a file with a sheet calculating the the expected omega
value along the sequence is here)

NOTE: the spreadsheet contains several sheets, these contain so many
rows (one for each generation) and columns (one for each parameter). The
scrol bar covers the sheet tabs. To switch from sheet to sheet press the
<ctrl ><Page Up> keys.|

Two sheets give the parameters for each run, the third contains the
combined values after the burnin, and the third one a bar graph for the
omega values along the sequence.

On the sheet giving the combined samples
after the burnin, scroll to the right to see the columns, starting with
pr+(1,2,3), pr+(4,5,6), and etc, and omega(1,2,3), omega(4,5,6) ... .
Calculate averages for each of these columns using the =AVERAGE()
function of Excel, enter the formula in a cell below the values for the
individual generations -- starting in column pr+(1,2,3) -- copy the
formula to all columns). You can do the same with the estimated omega
values (to the right of the pr(xyz columns). Plot the average values for
Omega for each column as a bar graph.

(Note: in my version of excel, this only worked when I specified the
cells from which to the average was to be calculated, i.e.,
=AVERAGE(BX2:BX1353) - when I highlighted the column, the resulting
formula would not copy correctly to the neighboring cells.)

**Note:** In some versions of
Excel the program refuses to plot more than 255values as a bargraph. To
solve this, copy the values you want to plot as bargraph, open a new
sheet, and do a special past of the **values only**.
To analyze the data for the extein and intein seperately: the
N-extein goes from nucleotide position 1 to 855 and the C-extein from
position 2686 to 3705. For codons or aa or columns in the bar graph, you
have to divide by three N-extein goes from position 1 to 285 and the
C-extein from position 896 to 1235.

Do you see different estimated P+ and Omega values for the extein and the Intein? If yes, which seems to be under stronger purifying selection? What are the average of the expected Omega values for sites in the N-extein What are the average of the expected Omega values for sites in the intein What are the average of the expected Omega values for sites in the C-Extein