# Difference between revisions of "Phylogenetics: HyPhy Lab"

Paul Lewis (Talk | contribs) (→Define the HKY85 rate matrix) |
Paul Lewis (Talk | contribs) (→Fixing up edge lengths) |
||

Line 234: | Line 234: | ||

=== Fixing up edge lengths === | === Fixing up edge lengths === | ||

− | The problem is that applying a new model to the 4 edge lengths caused the edge length information to be erased! We need to now set the <tt> | + | The problem is that applying a new model to the 4 edge lengths caused the edge length information to be erased! We need to now set the <tt>betat</tt> parameter for each of those edges to be compatible with the edge lengths that were originally there. This process is a bit more tedious than I would like, so you'll have to bear with me through this next section. The <tt>AUTOMATICALLY_CONVERT_BRANCH_LENGTHS = 1</tt> that we placed at the top of the file saves us from having to go through this procedure for all the <tt>hky1</tt> edges, but now that we've erased the models from 4 edges, we need to do a bit of work to get the edge lengths back. |

− | Recall that, under the JC69 model, the edge length equals <tt>v = 3*beta*t</tt>. The <tt>betat</tt> parameter that appears in our <tt>HKY85RateMatrix</tt> is the <tt>beta*t</tt> part, so to set the <tt>betat</tt> parameter we need to divide the desired edge length by 3: <tt>betat = v/3</tt>. The scaling factor 3 becomes more complex for the HKY85 model, but the principle is the same. Our first goal is thus to compute the scaling factor we need to convert edge lengths to <tt>betat</tt> values. | + | Recall that, under the JC69 model, the edge length equals <tt>v = 3*beta*t</tt>. The <tt>betat</tt> parameter that appears in our <tt>HKY85RateMatrix</tt> is the <tt>beta*t</tt> part, so to set the <tt>betat</tt> parameter we need to divide the desired edge length by 3: that is, <tt>betat = v/3</tt>. The scaling factor 3 becomes more complex for the HKY85 model, but the principle is the same. Our first goal is thus to compute the scaling factor we need to convert edge lengths to <tt>betat</tt> values. It is not necessary to create a function for this, but I will do it that way so in order to illlustrate how functions are defined in HBL: |

function computeScalingFactor(rateMatrix, baseFreqs) { | function computeScalingFactor(rateMatrix, baseFreqs) { | ||

Line 250: | Line 250: | ||

} | } | ||

scalingFactor = computeScalingFactor(HKY85RateMatrix, highATfreqs); | scalingFactor = computeScalingFactor(HKY85RateMatrix, highATfreqs); | ||

+ | |||

+ | The last line calls the function, passing the ''arguments'' <tt>HKY85RateMatrix</tt> and <tt>highATfreqs</tt> into the function. The function's ''parameters'' (<tt>rateMatrix</tt> and <tt>baseFreqs</tt>) are arbitrary names that are used within the function body. Thus, within the function, <tt>rateMatrix</tt> is used for the actual rate matrix <tt>HKY85RateMatrix</tt> that was passed into the function. (Defining a function would be more useful if we had to compute scaling factors for many different models.) The function body goes through two nested loops that visit every off-diagonal element of the rate matrix, multiply it by the base frequency of that row and the base frequency of that column and add this product of 3 terms to the growing sum <tt>sf</tt>. The function returns the value <tt>sf</tt>, which is used to set the value of the variable <tt>scalingFactor</tt>. Review your lecture notes if you don't remember why the scaling factor is computed in this way. | ||

+ | |||

+ | Now we simply need to set the <tt>betat</tt> values for the four edges of interest: | ||

+ | constrainedTree.microbats.beta := 0.097851/scalingFactor; | ||

+ | constrainedTree.Tonatia_bidens.beta := 0.008252/scalingFactor; | ||

+ | constrainedTree.Tonatia_silvicola.beta := 0.0001/scalingFactor; | ||

+ | constrainedTree.Pteropus.beta := 0.104663/scalingFactor; | ||

+ | |||

+ | Place the loop that shows edge lengths after these 4 lines in order to check and make sure the edge lengths have been set correctly. | ||

[[Category: Phylogenetics]] | [[Category: Phylogenetics]] |

## Revision as of 19:23, 24 February 2020

EEB 349: Phylogenetics | |

## Contents

## Goal

The goal of this lab exercise is to show you how to use the HyPhy program for data exploration and hypothesis testing within a maximum likelihood framework. Although much can be done with PAUP* and IQ-TREE, HyPhy lets you to do some interesting and useful things that these programs cannot, such as allowing the model of evolution to change across a tree.

## Login to Xanadu

Login to Xanadu and request a machine as usual:

srun --pty -p mcbstudent --qos=mcbstudent bash

## Loading modules needed

Load the module needed for this exercise:

module load hyphy/2.3.11 module load paup/4.0a-166

## The goal of this lab

In this lab, we will recreate the very interesting parametric bootstrapping analysis performed in this paper:

RA Vandenbussche, RJ Baker, JP Huelsenbeck, and DM Hillis. 1998. Base compositional bias and phylogenetic analyses: a test of the "Flying DNA" hypothesis. Molecular Phylogenetics and Evolution 13(3): 408-416. [DOI:10.1006/mpev.1998.0531 https://doi.org/10.1006/mpev.1998.0531]

In short, this paper demonstrated that the "Flying DNA" hypothesis proposed earlier by Pettigrew (1994. Flying DNA. Curr. Biol. 4: 277–280) was not viable. The Flying DNA hypothesis proposed that microbats and megabats are actually unrelated, but appear as a monophyletic group in phylogenetic trees due to the fact that both have high AT bias in the genes used to reconstruct phylogeny. The idea is that this strong nucleotide composition bias makes convergence much more probable, as there are effectively only two states (A and T) rather than four (A, C, G, T), and parsimony is mistaking such convergence as historical relatedness.

The Vandenbussche et al. paper simulated data under the null hypothesis (micro- and mega-bats are unrelated) but added various amounts of AT bias when simulating the bat lineages. If Pettigrew was correct, trees reconstructed from such data should show bats monophyletic, even though they were not together in the true tree used for the simulation.

This type of simulation required ad hoc software in 1998 because most software that can carry out simulations on phylogenetic trees assumes that the model is the same across the tree. Fortunately, these days we have HyPhy, which offers a way to simulate (and analyze) under pretty much any model you can imagine.

Vandenbussche et al. used the K80 (kappa=4) model across most of the tree, but the lineage leading to Pteropus (the lone megabat in the analysis) and the lineages within the microbat clade (Tonatia bidens, Tonatia silvicola, and their stem lineage) used HKY85 with kappa=4 but nucleotide frequencies that are AT-biased (e.g. piA=0.4, piC=0.1, piG=0.1, piT=0.4). The question is: how much AT-bias does one need to put into the simulation in order to see the convergence that Pettigrew claimed was happening. Our goal will be to recreate the parsimony part of figure 4 from the Vandenbussche paper. We will use HyPhy to simulate the data, and PAUP* to do the parsimony analyses.

## Creating the HyPhy batch file

HyPhy has its own scripting language known as the HyPhy Batch Language (HBL). HyPhy can be run from the command line to carry out phylogenetic analyses that are scripted in HBL, much like running Python to interpret a Python script. Start by creating a new file named `bats.bf` in a directory called `bats` (you can of course use any name you like, but these are the names I will be using). Also create a directory `simdata` inside your `bats` directory.

Download the data we will use for this analysis using curl:

cd ~/bats curl -LO http://hydrodictyon.eeb.uconn.edu/people/plewis/courses/phylogenetics/labs/irbp.nex

This is what your directory structure should look like at this point:

bats/ simdata/ bats.bf irbp.nex

Open `bats.bf` in your favorite text editor (e.g. nano) and enter the following text:

AUTOMATICALLY_CONVERT_BRANCH_LENGTHS = 1; ACCEPT_ROOTED_TREES = TRUE; ASSUME_REVERSIBLE_MODELS = -1; /***************************************************************************************** | Read in the data and store the result in the variable nucleotideSequences. */ DataSet nucleotideSequences = ReadDataFile("irbp.nex"); /***************************************************************************************** | Filter the data, specifying which sites are to be used. The first 1 means treat each | site separately (3 would cause the data to be interpreted as codons). The quoted | string (last argument) specifies which sites (where first site = 0) to use (we are | excluding sites with a lot of missing data and or alignment issues). This leaves us with | 978 sites rather than the 935 used by Vandenbussche, but it is impossible to determine | exactly which sites were excluded in the original study. */ DataSetFilter filteredData = CreateFilter(nucleotideSequences,1,"106-183,190-1089"); /***************************************************************************************** | Store empirical nucleotide frequencies in the variable observedFreqs. The 1,1,1 means | unit=1, atom=1, position-specific=1. These settings create one global set of | frequencies (setting, for example, unit=3, atom=3, would tally 64 codon frequencies, | which is not what we need because we will not be using a codon model). */ HarvestFrequencies(observedFreqs, filteredData, 1, 1, 1); fprintf(stdout, observedFreqs);

### Run bats.bf

Run your `bats.bf` file as follows:

hyphy bats.bf

You should see the empirical nucleotide frequencies displayed as follows:

{ {0.1936054742803209} {0.3104845052697813} {0.3106418121755545} {0.1852682082743432} }

I have provided copious comments in the batch file to explain what each command is doing. Comments in HyPhy batch files follow the C programming convention: either surround the comment with slash-asterisk delimiters (`/* comment */`) or begin the line with two slashes (`// comment`). The three lines at the beginning will be explained when each becomes relevant.

The `fprintf` command also mimics the C programming language. It allows you to print objects to `stdout` (standard output). This is a useful tool for performing sanity checks, such as checking to ensure that the frequencies were indeed harvested and stored in the variable `observedFreqs`.

### Define the HKY85 rate matrix

Add the following to the bottom of your `bats.bf` file:

/***************************************************************************************** | Define the KHY substitution matrix. '*' is used for the diagonal elements that can be |.computed automatically by HyPhy. The transition-transversion rate ratio (kappa) is | declared to be global, meaning it is shared by all edges. */ global kappa = 4.224618; HKY85RateMatrix = {{ * , betat , betat*kappa , betat } { betat , * , betat , betat*kappa } { betat*kappa , betat , * , betat } { betat , betat*kappa , betat , * }}; fprintf(stdout, HKY85RateMatrix);

Run `bats.bf` in HyPhy. Did the `HKY85RateMatrix` variable have the value you expected? Why or why not? (Hint: the variable `betat` was initialized to zero.) Set `betat = 1;` before the `fprintf` statement and run the batch file again to see the effect. Does the output make sense now?

### Combine frequencies with rate matrix to create an HKY85 model

/***************************************************************************************** | Define the HKY85 model, by combining the substitution matrix with the vector of | empirical frequencies. */ Model hky1 = (HKY85RateMatrix, observedFreqs);

Now we have a model variable (`hky1`) that can be applied to each edge of a tree.

### Create a tree representing the null hypothesis

The next step is to create the model tree that we will use for the simulations. The tree topology is from Figure 2a in the Vandenbussche et al. paper. I have estimated the edge lengths and transition/transversion rate ratio (kappa) using PAUP*.

/***************************************************************************************** | Define the tree variable, using the tree description read from the data file. | By default, the last defined model (hky1) is assigned to all edges of the tree. */ Tree constrainedTree = "((((((Homo:0.077544,(Tarsius:0.084863,Galago:0.075292):0.009462):0.026367,((Cynocephalus:0.067955,Tupaia:0.093035):0.016468,(Oryctolagus:0.093866,Mus:0.143079):0.013506):0.017052,Pteropus:0.102675):0.008768,Bos:0.099273):0.007976,(Tonatia_bidens:0.008137,Tonatia_silvicola:0)microbats:0.096022):0.013987,Felis:0.044428):0.043248,Didelphis:0.247617)";;

### Creating the second model

Hopefully everything we've done so far in HyPhy makes sense. Now comes the tricky part! We need to create a second HKY85 model that has elevated A and T frequencies, and that model needs to be applied to only certain edges in the tree. Right now, the `hky1` model has been applied to every node in the tree. You can verify this using the following HBL code:

GetInformation(modelMap, constrainedTree); fprintf(stdout, modelMap);

Add these two lines to the bottom of your `bats.bf` file and run it. You should see this near the bottom of your output:

{ "Bos":"hky1", "Cynocephalus":"hky1", "Didelphis":"hky1", "Felis":"hky1", "Galago":"hky1", "Homo":"hky1", "Mus":"hky1", "Node1":"hky1", "Node10":"hky1", "Node11":"hky1", "Node14":"hky1", "Node19":"hky1", "Node2":"hky1", "Node3":"hky1", "Node4":"hky1", "Node5":"hky1", "Node7":"hky1", "Oryctolagus":"hky1", "Pteropus":"hky1", "Tarsius":"hky1", "Tonatia_bidens":"hky1", "Tonatia_silvicola":"hky1", "Tupaia":"hky1" }

First, create a second model named `hky2` and apply it to four specific edges in the tree:

/***************************************************************************************** | Define a second AT-rich HKY85 model named hky2 and apply it to selected edges. */ highATfreqs = {{.45}{.05}{.05}{.45}}; Model hky2 = (HKY85RateMatrix, highATfreqs); SetParameter(constrainedTree.microbats, MODEL, hky2); SetParameter(constrainedTree.Tonatia_bidens, MODEL, hky2); SetParameter(constrainedTree.Tonatia_silvicola, MODEL, hky2); SetParameter(constrainedTree.Pteropus, MODEL, hky2);

If you move your `GetInformation` call and its accompanying `fprintf` after these 6 lines, then you should see this in the output:

{ "Bos":"hky1", "Cynocephalus":"hky1", "Didelphis":"hky1", "Felis":"hky1", "Galago":"hky1", "Homo":"hky1", "Mus":"hky1", "Node1":"hky1", "Node10":"hky1", "Node11":"hky1", "Node14":"hky1", "Node2":"hky1", "Node3":"hky1", "Node4":"hky1", "Node5":"hky1", "Node7":"hky1", "Oryctolagus":"hky1", "Pteropus":"hky2", "Tarsius":"hky1", "Tonatia_bidens":"hky2", "Tonatia_silvicola":"hky2", "Tupaia":"hky1", "microbats":"hky2" }

### Don't get too comfortable just yet

It is great to sit back and admire your work thus far, but there is a small problem looming just below the surface. Let's print out all the edge lengths in the tree:

fprintf(stdout, "\n\nCurrent edge lengths:\n"); edgeNames = BranchName(constrainedTree,-1); edgeLengths = BranchLength(constrainedTree,-1); for (k = 0; k < Columns(edgeNames) - 1; k = k + 1) { fprintf(stdout, Format(edgeLengths[k],10,5), " ", edgeNames[k], "\n"); }

The first line simply skips a couple of lines (`\n\n`) and prints a header announcing that current edge lengths will follow.

The second and third lines ask HyPhy for the names of all branches and the lengths of all branches. The `-1` in these functions is somewhat obscure, but means "give me all of them". (If you used `5` rather than `-1`, it would give you the name of the edge having index 5.)

The last 3 lines are a loop in which the variable `k` ranges from 0 to the number of edges minus 1. For each value of `k`, the `fprintf` statement prints out the length of edge `k` (the `Format` command causes it to use exactly 10 spaces and 5 decimal places) followed by a couple of spaces and then the name of edge `k`. The newline character (`"\n"`) at the end of the `fprintf` statement causes a carriage return so that the edge lengths and names do not all end up on the same line of output.

After running `bats.bf`, what possibly important detail do you notice about the lengths of the edges to which we attached the `hky2` model?

### Fixing up edge lengths

The problem is that applying a new model to the 4 edge lengths caused the edge length information to be erased! We need to now set the `betat` parameter for each of those edges to be compatible with the edge lengths that were originally there. This process is a bit more tedious than I would like, so you'll have to bear with me through this next section. The `AUTOMATICALLY_CONVERT_BRANCH_LENGTHS = 1` that we placed at the top of the file saves us from having to go through this procedure for all the `hky1` edges, but now that we've erased the models from 4 edges, we need to do a bit of work to get the edge lengths back.

Recall that, under the JC69 model, the edge length equals `v = 3*beta*t`. The `betat` parameter that appears in our `HKY85RateMatrix` is the `beta*t` part, so to set the `betat` parameter we need to divide the desired edge length by 3: that is, `betat = v/3`. The scaling factor 3 becomes more complex for the HKY85 model, but the principle is the same. Our first goal is thus to compute the scaling factor we need to convert edge lengths to `betat` values. It is not necessary to create a function for this, but I will do it that way so in order to illlustrate how functions are defined in HBL:

function computeScalingFactor(rateMatrix, baseFreqs) { sf = 0; for (n1 = 0; n1 < 4; n1 = n1+1) { for (n2 = 0; n2 < 4; n2 = n2+1) { if (n2!=n1) { sf = sf + baseFreqs[n1]*baseFreqs[n2]*rateMatrix[n1][n2]; } } } return sf; } scalingFactor = computeScalingFactor(HKY85RateMatrix, highATfreqs);

The last line calls the function, passing the *arguments* `HKY85RateMatrix` and `highATfreqs` into the function. The function's *parameters* (`rateMatrix` and `baseFreqs`) are arbitrary names that are used within the function body. Thus, within the function, `rateMatrix` is used for the actual rate matrix `HKY85RateMatrix` that was passed into the function. (Defining a function would be more useful if we had to compute scaling factors for many different models.) The function body goes through two nested loops that visit every off-diagonal element of the rate matrix, multiply it by the base frequency of that row and the base frequency of that column and add this product of 3 terms to the growing sum `sf`. The function returns the value `sf`, which is used to set the value of the variable `scalingFactor`. Review your lecture notes if you don't remember why the scaling factor is computed in this way.

Now we simply need to set the `betat` values for the four edges of interest:

constrainedTree.microbats.beta := 0.097851/scalingFactor; constrainedTree.Tonatia_bidens.beta := 0.008252/scalingFactor; constrainedTree.Tonatia_silvicola.beta := 0.0001/scalingFactor; constrainedTree.Pteropus.beta := 0.104663/scalingFactor;

Place the loop that shows edge lengths after these 4 lines in order to check and make sure the edge lengths have been set correctly.