Models of Lineage-Specific Rate Variation

December 27, 2013. Unfortunately, because of the conjunction of Revels and winter break, I will not be meeting with Dr. Miller for almost a month! However, most of our work now involves in reading papers online and using computer programs, so the long break time won't be too big of a problem. For a little recap, our goal is to find the divergence time of mycobacterium strains identified in our samples. We are going to construct a tree and calculate the divergence time by using compute program BEAUTi (Bayesian Evolutionary Analysis Utility) and BEAST (Bayesian Evolutionary Analysis Sampling Trees) (very creative names). BEAUti is a program with a graphical user interface for creating an input file for BEAST, which muct be written in XML language. The part that really takes time is deciding which models to use in estimating the divergence time of bacteria. Therefore, what we are doing now is reading journals and previous studies on similar experiments to help us decide our models.

The first type of models introduced here is lineage-specific rate variation. A tutorial mentioned five common models:

• global molecular clock
• local molecular clock
• compound poisson process
• autocorrelated rate
• uncorrelated rate

Sebastián Duchêne's  power point provides simplified explanation of those molecular clock, in which I found very helpful.

Different color represents different rate of substitution. Relaxed clocks include auto- and uncorrelated rates. Credit to S.  Duchene (page 8).

After reading several papers and websites, I decided to try global and local clocks because bacteria evolve quickly in a short period of time. Therefore, the lineage-specific rates variation of different strains are more likely to be same or similar to one and another.

While I am still absorbing more info, I wish everyone happy holidays! I can't believe another year has passed and that I will be graduating in less than 6 months! Nonetheless, I hope to have significant contribution to the bird project before I graduate:)

Happy Holiday!

Today (November 22, 2013), we had a short talk on a phone with the zoo people whom we will be working with from San Diego. They are a group of very vibrant and intelligent people even without seeing their faces. Although most of the time I was listening to the conversation as there was a lot to digest, I was nonetheless very excited for this project! After several follow-ups emails, my goal was narrowed down to a single, clear task: What is the divergence time for these avium micrbacterium strains?

Dr. Miller and I looking at data on our computer. Most of our work can be achieved by computer programs, so no fancy gears or species to show B-)

We won't be meeting next Friday because it's Thanksgiving break (Black Friday shopping!!)! Happy Holiday everyone!:)

Bayesian Inference vs. Maximum Likelihood

Last Friday (November 8, 2013), my mentor and I discussed a bit about the Bayesian Inference of Trees and the program BEAST.

Like Maximum Likelihood (ML), Bayesian Inference is a character-based tree method, and they both generate several trees and use some criterion to decide which tree is the best. However, BI differs from ML in that BI "seeks the tree that is most likely given the data and the chosen substitution model" whereas ML "seeks the tree that makes the data the most likely" (Hall, 140). Sounds the same, I know. However, after a short period of intense research and seeing a lot of alien equations, I interpreted the major distinction like this (please correct me if I am wrong, clarification needed!!):

• ML - make  different trees and calculate the likelihood of each tree
• BI -  finds the tree that constitute the X (observation) with given knowledge (i.e. the descendants or a substitution model).

Mathematically, according to Bayes' Theorem:

In Bayesian inference, the event B is fixed (the "given knowledge") in the discussion, and we wish to consider the impact of its having been observed on our belief in various possible events A (how the tree split). In such a situation the denominator of the last expression, the probability of the given evidence B, is fixed; what we want to vary is A. Thus, the posterior probabilities are proportional to the numerator: 

• P(B|A) is the prior probability, also the likelihood. It can be interpreted as P(data | hypothesis). Prior probability indicates our state of knowledge about the truth of a hypothesis before we have observed the data. 
• P(A|B) is the posterior probability. We can rewrite as (hypothesis | data). It shows how well our models agree with the observed data --> Bayesian

I drew a picture and hopefully it helps with visualization. Suppose we know there are species A, B, C, E, and F. For ML, we are finding in what way are the species located on the tree that P(A)xP(B|A) is the largest. For BI, we know from data that A, D, E, B are in the order they are. Given this knowledge, we are finding where C and F locate so that P(A|B) is the greatest.

This is at least what I've got so far. Nonetheless, BI can be burdensome if no obvious prior distribution for a parameter exists. The researcher to ensure that the prior selected is not inadvertently influencing the posterior distribution of parameters of interest. I clarified my thoughts every time I wrote it out. Before we get into the bird data, we try to ask ourselves some fundamental questions about the phylogenetic trees. Next Friday, I will be joining a talk with the data holders to discuss our participation in the Avian mycobacteria project. 

More Phylo Talks!

Last Friday (November 1st, 2013), My mentor and I discussed more about the bootstrapping methods in terms of how it work. He gave me a paper, in which I will be reading this week. We also talked about the derivation of the Jukes Cantor model, which is the simplest evolutionary model used to predict the rate at which nucleotide substitutions occur during evolution. It has two main assumptions:

• Equal frequencies of the four bases
• The probability of changing from one state to a different state is always equal, i.e. A-->G is as likely to happen as G-->C

Q Matrix of Jukes Cantor Model

f(t) represents the mutation rate, which is equal for all if a nucleotide is substituted with a different nucleotide. The probability of changing from A to A is 1- f(t) since the sum of the probabilities should add to 1.

However, there are many modifications for this model later published, including Felenstein 81 (F81), Kimura 2-parameter (K2P), HKY85, TN93, GTR, etc.

While we are still waiting for the bird data for my project, Dr, Miller sent me a practice data, and he would like me to explore his website on my own first (http://www.phylo.org/). So I looked at the demo, uploaded the data, and let the data run by accepting the default setting. I found the website pretty user-friendly, yet I don't really know what those result means. Hopefully I will learn more about what I can get from these results this Friday.

Neighbor Joining Trees

October 25, 2013. For the past two weeks, I had been reading about the introduction of several major methods for estimating phylogenetic trees. There are primarily two approaches to tree estimation:

• Algorithmic - uses an algorithm to estimate a tree from the data. Advantages of using this method include fast speed and yielding only a single tree. Algorithmic method includes Neighbor Joining, UPGMA
• Tree-searching - estimates many trees and then uses some criterion to decide which is the best tree of all.

Another way to categorize those methods are distance vs. character-based methods.

• distance - Neighbor Joining, UPGMA
• character-based - Parsimony, Maximum Likelihood, Bayesian Interference

Our discussion today focused on Neighbor Joining Method (NJ). Nonetheless, it is very important to acknowledge that the 'right' tree does not exist since we cannot know what exactly happened in the past. All of these methods only allow us to deduce the the order in which existing taxa (sequences) diverged from a hypothetical common ancestor, and to calculate the amount of changes along the branches between the diverging events.

NJ is one of the most popular distance algorithmic method. It produces a single, strictly bifurcating tree (meaning that each internal node has exactly two branched descending from it). I downloaded another file for practice, in which I am going to show you here.

First I opened up the file LargeData.meg from MEGA. The window shows DNA sequence alignment. The program only show a base when it is different from that of the first sequence. Otherwise it'll just show "."

Then we calculate the average Jukes-Cantor (JC) distance for our data. The data are not suitable for NJ if JC >1.0 In this case, I got 0.811 for this data.

Next, I constructed the tree by choosing Construct/Test Neighbor Joining Tree from Phylogeny menu. A window of Option Summary will pop up. While it is good to stay with default setting (Maximum Composite Likelihood), my mentor and I discussed several other models for the Model/Method section (for simplicity, I won't get into detail here). 

• Jukes-Cantor model
• Kimura 2-Parameter model (K2P)
• Tamura-Nei model

In order to test the reliability of our tree, we also entered "bootstrap method" for Phylogeny Test with a replication of 1000 times. Bootstrapping is, in essence, a method to test how likely a certain "split" is going to occur by using computational simulation.

As we hit compute...

The number assigned to each internal node is the bootstrap value (bootstrap percentage). The high number shows that the split is not random. We can resize the tree a little bit to get a better view by clicking Display Only Topology.

Each internal node only splits into two branches (strictly bifurcating)

The bootstrap value less than 50% means we really have no idea what the branching order is. Scientists usually collapse those branches into polytomies - nodes from which more than two branches descend.

Majority rule tree

Finally, we can do various things to change the appearance of the tree to make it best represent our data...

A circular tree is an unrooted tree

This chapter was quite hard for me because it involved in a lot of equations and concept. Nonetheless, I am happy that I constructed a descent-looking phylogenetic tree in the end! Next week I will discuss with my mentor more about the bootstrapping method and hopefully it will make more sense to me!:)

More MEGA 5 with Protein Sequences

Last Friday (October 4th, 2013) I learned more about using MEGA 5 and Blast for constructing phylogenetic trees. This time, however, we used protein blast instead of nucleotide blast because searching protein sequences can detect much more distantly related homologs than searching nucleotide sequences. If two sequences are "homologous," we assume that they descended from a common ancestor (need to be distinguished from "similarity"). Recall that an amino acid is coded by a codon, and that the same amino acid can be coded by several codons, a silent substitution will not cause the amino acid sequence to change. For DNA sequences, there are only 4 possible states (A,T,C,G) of each characters, so when the sequences diverge greatly that there are only 25% identical, the program would classified them as not closely related even though they may have very similar amino acid sequences if you translate them. Thus, the solution id to use protein sequence as query. Proteins have 20 possible states so the lower limit of detectable homology drops to about 5% (instead of 25%). 

For this exercise we used EbgC protein sequence as a query and use blastp.

Those protein sequences are identified as closely related to our query sequences. If you click on the first entry...

The first protein sequences actually include 1087 DNA sequences. Those DNA sequences all code for the same protein sequences even though their DNA sequences vary (silent mutation). Therefore, all those 1087 sequences are very closely related. If we originally did a blastn instead of blastp, the result would only show 100 of these 1087 sequences, and we wouldn't even know other distantly related sequences. The second entry (evolved beta-galactosidase subunit beta [Escherichia coli]) is an example of distantly related sequences (even though they are still pretty close).

Later I chose several protein sequences and translated them back to DNA sequences to align them on MEGA and established a tree. The whole process took so long! Nonetheless, the main idea here is that protein sequences give us a bigger picture regarding the relatedness between species. We start from DNA sequence --> protein sequence --> protein structure --> functions. The bigger we look at , the more distantly related species we are incorporating.

I will be learning some computer code for the command line for the next few weeks. However, I will not meet with Dr. Miller for the next two weeks because of the schedule conflicts. In the meantime, I will continue reading the book and get myself more familiar with blast and MEGA:)

Tutorial: A Phylogenetic Tree with MEGA 5

Last Friday (September 27, 2013) instead of having an internship, I read the tutorial book my mentor gave me and made my first phylogenetic tree!

First of all, I had to download a program called MEGA 5, which allows us to align and compare related sequences with our sequence of interest. For this tutorial, we use a sequence of an alpha-glucuronidase gene from the bacterium Thermotoga petrophilia. We can obtained the sequence by clicking Do BLAST Search from the Align menu, and it will lead us to the following window:

MEGA 5 is linked to BLAST.

I entered the "accession numbers" and the "query subrange", which are given by the book, chose "Neucleotide collection" for the Database, and blastn for Program Selection.

The results showed up few seconds after I hit search. It gave all the sequences that produce similar alignment to our query sequence. The first entry was eventually my query sequence (because that's what I searched for!) Each entry consists of 8 elements: accession number, description, Max score, Total score, query coverage, E-value, and Max identity. For the sake of time, I will just briefly explained Max score, Total score, query coverage, and E-value. Based on the definitions given by the book,

• Max score - the score for the highest scoring segments of the subject sequence
• Total score - the sum of the scores for all the segments that aligned (including noncontiguous alignments)
• Query coverage (%) - shows how much of the query sequence aligns with the subject sequence.
• E-value - describes the number of hits with scores this high that one would "expect" to see by chance when searching a database of a particular size. The closer the E-value to 0 the better.

 The book asked us to include those subject sequences that have E-value <10^-3 and query coverage >60%. For that, I chose 9 sequences, and I clicked "Add to Alignment" for each (some sequenced require to be reversed before added to alignment).

Sequences of the nine samples

Then I clicked "Align DNA" and get:

This is the end portion of the sequences after being aligned

Finally, I construct a phylogenetic tree from my DNA alignment by choosing "Construct / Test Neighbor Joining Tree" and tada~

My first phylogenetic tree!

This Friday, I will discuss with Dr. Miller more about the assumptions of making phylogenetic tree and do more practices with MEGA 5. I can't wait!:)