Pick one of three gene sets to work on, MHC-class1, CYP2, or
SOD1.
MHC-class1 - one of the two proteins that binds and presents
antigenic peptides for the immune system (this is the same gene as
in the Swanson paper, but not the same species).
CYP2 - a family of P450 enzymes critical for degredation of
xenobiotics (including drugs). Thomas, 2007.
SOD1 - superoxide dismutase enzyme, important for oxidant
protection.
Save the following three files (or the equivalent files for CYP2
or SOD1):
MHC-class1_codons.paln - a phylip format codon alignment of MHC Class I genes from 19 mammals, named by species
MHC-class1_tree.newick - a Newick format maximum-likelihood tree for the same genes
MHC.ctl - the file that
codeml uses to set parameters and I/O files for a run
[CYP2_codons.paln, CYP2_tree.newick, CYP2.ctl,
SOD1_codons.paln, SOD1_tree.newick, SOD.ctl]
Each of the sets is coding sequences from various mammals. They
have been modified a bit to make viewing and interpretation
simpler (alignment degapped and names simplified).
XXX below depends on which data set you chose (MHC-class1 or CYP2 or SOD1).
1) Open the XXX.ctl file and look at the parameter settings,
which are commented and most of which you should now understand.
(the file name doesn't matter, I use the ".ctl" extension to
suggest "control" file)
2) Put those three files in the same directory as the codeml program.
3) Open a command prompt, navigate to your directory and type
"codeml XXX.ctl"
4) Start the runs as soon as possible - they take a LONG
time to finish. While they run I will make some comments and
explanations.
When the run has finished, open the file "XXX_codeml_out.txt" and
look through it.
Find the relevant log likelihoods (model 7 and model 8). They are
kind of buried in a lot of other stuff, but you can find them if
you look hard - the label is lnL and the likelihood is the large
negative number following the colon. Also look at the omega
classes, following the line "dN/dS (w) for site classes (K=11)" -
remember that the discretized distribution has equal class
frequencies (p) but is allowed to vary in omega value (w).
Open the file "rst" (it is a text file), scroll down to the part
that starts with "Bayes Empirical Bayes (BEB) probabilities for 11
classes" and interpret what you are seeing. Note - you have no
control over this file name and if you make another codeml run in
the same directory it will be overwritten, along with the other
mysterious files that appeared during your run.
Write up an interpretation of 1) the statistical test comparing
models 7 and 8 using chi-square, two degrees of freedom (the two
degrees of freedom are the 11th omega class and its estimated
omega value). 2) what the omega classes tell you about the
fraction of codons with varying degrees of negative selection. 3)
what the BEB part of the rst file is telling you.
If you get ambitious (totally optional), map the positively
selected sites onto the crystal structure using Cn3D (go to the
NCBI structure database to find structures and Cn3D,
https://www.ncbi.nlm.nih.gov/Structure/). Be warned that you will
have to figure out how the amino acids in your data set correspond
to the protein used for the crystal structure (get the crystal
structure amino acid sequence, translate the codons from one of
the sequences from your data set and align them).