---
title: "Workshop 1: Analysis of low throughput sequence data and phylogeny"
author: ""
output: 
  html_document:
    depth: 1
    theme: cerulean
    toc: yes
    toc_float:
      collapsed: no
      smooth_scroll: yes
  word_document: default
  pdf_document:
    toc: yes
---

[Link back to main course page](http://www-users.york.ac.uk/~dj757/popgenomics/popgenomics_index.html)

This workshop was originally designed by Francois Balloux (UCL). Daniel Jeffares and Ville Friman (York) edited and updated it.

##Introduction

**Phylogenetics** is a relatively complex field, with a variety of tools adapted to very specific questions and analyses. It is also a field which has not clear-cut boundaries with related disciplines such as population genetics, statistics and bioinformatics. 

As such, it is not reasonable to expect an in depth treatment of any aspect of phylogenetics within an afternoon.  This practical has been devised to give a glimpse of the vast publicly available sequence resources and illustrate the kind of research questions that can be addressed. While the practical is (hopefully) easy to follow as a cookbook recipe, the program should be sufficiently light to leave some time to explore the databases and the **MEGA software**. 


#Part 1: Molecular Change in the Influenza Virus.

##Getting data

[Genbank](https://www.ncbi.nlm.nih.gov/genbank/) is one of three main databases that store DNA sequence data. GenBank is part of the International Nucleotide Sequence Database Collaboration, which comprises the DNA DataBank of Japan (DDBJ), the European Nucleotide Archive (ENA), and GenBank at NCBI. These three organizations exchange data on a daily basis.

There is a wealth of freely available sequence data available on Genbank. The downside is that most data is often poorly annotated and misses the critical associated information (e.g. date and place of collection). Here we will take advantage of the new [influenza virus resource database](https://www.ncbi.nlm.nih.gov/genomes/FLU/Database/nph-select.cgi?go=database), which is relatively well curated.

***


##Generating a dataset
If you go to the influenza virus resource database by clicking on the link [here](https://www.ncbi.nlm.nih.gov/genomes/FLU/Database/nph-select.cgi?go=database).

You should see a page like this:

![](pics/Picture1.png)

Figure 1: Influenza virus resource website frontpage

***

As an exercise, we can download all the human full-length H3N2 hemagglutinin (HA) sequences from Wellington (New Zealand). To do so you have to fill the query site as in figure 2. 

![](pics/Picture2.png)
Figure 2. First, click Nucleotide and then select following Type: A, Host: Human, Country/Region: any, Segment: HA, Subtype: H3N2. Finally, tick full length only.

***

If you now click on “Show results”, you will open a new page with a list of sequences (figure 3). You should get 72 strains from 1985 to 2005. First, sort these by year.

![](pics/Picture3.png)
Figure 3. List of queried sequences

***


**You can now align your sequences.** This is usually done with dedicated programs but the online tool on the influenza virus resource site is remarkably accurate and fast. Thus, click the “Do multiple alignment” tab. About a minute later or so, the alignment should be finished. Download it and give it a meaningful name (e.g. H3N2_Wellington_seq.fa). This is a Fasta file (.fa extension), which is one of the most standard format for sequence data. We can load this aligned sequence file into MEGA.

***

##Visualising and analysing data

***

###MEGA software

MEGA is a tool for phylogenetic and other evoltionary analysis of DNA. It also produces multiple sequence alignments. It is reasonably neat package, even if somewhat counterintuitive at first. It uses various methods to make phylogenetic trees, analyse and display the trees. It offers three distance-based tree reconstruction methods (Neighbour Joining, Minimum Evolution and UPGMA) and also does Maximum Parsimony.

***

###Download and install MEGA
We will use MEGA X. First, download MEGA X (Windows 32 Bit version) from here: http://www.megasoftware.net/.

***

###First steps with MEGA

If you start the Mega program, you will get the following window. 

![](pics/Picture4.png)
Figure 4. Mega starting window


**To explore its features**, we recommend you first go to the Tutorial (under Help Docs, bottom left of screen). Skip the “Sequence alignment section”, but try to go briefly through the “Building Trees from Distance Data” and the “Computing Statistical Quantities for Nucleotide Sequences” tutorials by using the inbuilt examples.

***

###Running MEGA: loading data

Once you got somewhat accustomed to the idiosyncrasy of the multiple windows, load the H3N2_Wellington_seq.fa file through the **Data > Open a File/Session** Menu (Ctrl+O). 


Then click **Analyse** (not align), click OK to **nucleotide sequences** to **Protein coding** and to **Standard Code**. 

You will get the window displayed in figure 6. This datafile comprises 72 complete HA sequences from influenza H3N2 sampled between 1985 to 2005. Sequences have also been labelled according to their year of collection.


![](pics/Picture5.png)


Figure 5. Wellington set of 72 strains from 1985 to 2005.


***

###Running MEGA: grouping strains (samples)


**Now click on the “TA” window at top left** to see the Sequence Data Explorer. You can scroll around and resize this window.

Now let’s **group the strains by year**. This will be fairly easy because we sorted the data by year before we downloaded it. Samples have the year embedded in the same. 

*For example, this one is from 1985*

Influenza A virus (A/Wellington/4/**1985**(H3N2)). 

*And this is from 2004* 

Influenza A virus (A/Wellington/34/**2004**(H3N2))

***

**To group strains**, select one or more strains from a year by clicking on the name of the first strain in the grey part of the window, and shift+click to select the last one, as below.

Then choose the Groups menu, then Add/Edit Group name. In this case we will group strains by year.

![](pics/Picture6.png)

Figure 6. Grouping strains in MEGA.

***

##Building a tree in MEGA

There are a large number of possible options for analysis in MEGA. Feel free to explore the various tools and methods. 

We will build one tree along a single, robust methodology. 

* Go to **Phylogeny** menu on the main menu. 

* Choose **Construct/Test Neighbor-Joining Tree**’, say **Yes** to use currently active data.

* Retain **Bootstrap method** with 1000 replications (Bootstrap represents statistical support for individual clades, values in excess of 70-80% are considered as indicative of well supported).

* Choose a **Kimura 2-parameter model** of evolution in **Model/Method** - otherwise leave all other default options. 

* After less than a minute, you will get a phylogenetic tree. *Have a careful look at it. Discuss it with your group and a tutor*.

* The different menus allow you to change the presentation. *Try some of these out*.

*** 

###Drawing a clearer tree

* It is possible to define the root of the tree. The clade with the sequences from 1985 at the bottom makes a biologically reasonable root. To define a root, simply right click with your mouse on the chosen node and select **Place root**. 

* One obvious problem with the graphical representation is that the number of taxa is very large. Nicer trees can be generated switching off the tames of the strains, and just showing the group names. 

* You can do this by clicking on **Toggle the display of taxa names**, the lowest button the left hand button bar. You should see a tree like Figure 7.


![](pics/Picture7.png)

Figure 7. Years shown for Wellington H3N3 HA sequences.

***

<span style="color:purple">**DISCUSS WITH YOUR GROUP**</span>

<span style="color:purple">What pattern do you notice about this tree?</span>

***

###Measuring molecular change over time

An important question in phylogenetics is whether the accumulation of mutation (substitution rate) is constant over time. One way to test for this would be to count the number of mutations from the root to each tip. Estimating the most likely root is not completely straightforward and counting the number of mutations from root to tip would require some scripting well beyond the scope of this practical. However, we can test whether genetic distances increase linearly with time.

***

<span style="color:purple">**DISCUSS WITH YOUR GROUP**</span>

<span style="color:purple">Why might genetic distances increase linearly with time?</span>

<span style="color:purple">In what situations might substitution rates *not* be constant over time?</span>


***

**To address the issue of linearity between substitution rates and time**, we can first compute genetic distances between the sequences grouped by years. To do this, click **Compute Between Group Means**, in the Distance menu (light blue master window). Keep Kimura 2 parameter model and click Compute. You will get the matrix shown in figure 8.


![](pics/Picture8.png)
Figure 8. Genetic distance matrix between different years 


Click on the **File** menu, then **Export/Print Distances**, choosing the **unformatted text** and **column** (not matrix) options as Figure 9, below.


![](pics/Picture9.png)
Figure 9. Output genetic distances from MEGA.

***

###Genetic distance and the molecular clock

In this part, we will test if influenza genomes change gradually with time in a [Molecular clock](https://en.wikipedia.org/wiki/Molecular_clock). If you are nor familiar with the molecular clock read about it [here](https://en.wikipedia.org/wiki/Molecular_clock).

The remainder of this part of the workshop is we’ll carry out in **Microsoft Excel.** So start Excel now.

No load the the genetic distance values, using the **Data**, **Get External Data**, **Import Text File...** menus, as below.

![](pics/Picture10.png)
Figure 10. Importing data in to Excel.

You'll need to remove some header lines, and lines at the end of the file, until you have only the column headers, and the data.


Then highlight the data, then choose **Data**, **Text from columns**, choose **Delimited**, and delimit by **space**. And thern calculate the difference between years, as below. Then copy this calculation cell to all the columns below.

![](pics/Picture11.png)

Finally, make a scatter plot, using the *Genetic Distance* and *Time difference (years)* columns.

***

<span style="color:purple">**DISCUSS WITH YOUR GROUP**</span>

<span style="color:purple">Does the scatter plot confirm to expectations about a molecular clock?</span>

<span style="color:purple">What does this mean about the *speed of molecular change* in this virus?</span>


***


#Part 2: Does the influenza virus spread across the Tasman Sea?

##Introduction
Another use of virus genetic data is to determine how fast pathogens spread. Here, we’ll show an example where we examine whether people in New Zealand share the same strains as their Australian neighbours across the Tasman Sea.

***

![](pics/Picture12.png)
Figure 12. There are frequent flights between Syndey and Auckland. Do people carry the influenza virus with them?

**If you like, you can choose any two countries, and test these.**

***
##Getting data, from both sides of the Tasman Sea.

To do this, go back to the influenza database [here](https://www.ncbi.nlm.nih.gov/genomes/FLU/Database/nph-select.cgi?go=database).

**Then:**

**1.** Search for strains from Southern Temperate region, HA, H3N2 from year 2000 to year 2000. This will find strains from Australia and New Zealand.

**2.**	So that you don’t give MEGA too much data, choose about 100 nucleotide sequences, from this subset in some way (eg: choose all from females). 

**3.**	Sort the samples by virus name, which happens to sort by country and region. We’ll use this to make groups in MEGA. 

**4.**	Align, as before

**5.**	Output fasta file, as before. And save the file with a sensible name (eg: **tasman.fa**).

***

Then open MEGA, and:

**1.**	Add groups to sets of sequences from the same regions within Australia and NZ. Be sure to name groups to that they include country prefixes, like NZ- and AU-. We’ll use these later.

**2.** Compute inter-group distances, as before.

**3.** Output the distance file, as before.


***

#Part 3: Using R Studio

We will now start to use the statistical package **R** to analye the tree data. Although  **R** can be challenging to learn at first, it is very powerful. R can conduct any sophicated analysis and create beautiful plots that extend far beyond the power of Excel. We will use a popular R graphic user interface (GUI) called **R Studio**. Both R and R Studio are free.

First, we will use R Studio to re-analyse the correlation of genetic distance and time distance that we conducted with Excel.

***

##First steps with R Studio

###Make a folder
The words ‘folder’ and ‘directory’ mean the same thing (for most purposes). Folder is mainly a Windows term; directory is more classical and widespread. It is best to work in one folder named for the particular task. I recommend creating a folder called Population Genomics (or similar) containing folders called ‘data’ and ‘scripts’.

Then start RStudio from the Start menu

In RStudio, set your working directory to your scripts folder, as below.

![](pics/Setwd.png)


Then make a new script file to carry out the rest of the work. Give your script a sensible name. It is traditional to give R scripts the suffix **.R**, like **genomics2019.R**, for example.

![](pics/Newscript.png)

##Make a plot and test a correlation

Now we will load some data into R, and make a plot. Copy or type the commands below into your script then select the **run** command at the top of the script window. You may need to change the *wellington-distances.txt2* file name in your script.

Note: Any lines in your script that start with a **hash** symbol (#) are comments, that are not run as code.


```{r,eval=F}
#read your distances file into a data frame called 'flu'
#PS: your distance file will probably have a different name
#So edit the command to have that name
flu <- read.table("wellington-distances.txt2",h=T)
```


Now we will calculate the difference between the two dates. This is the same as we did with excel. But easier, once you know how. Your columns will need to be named species1 and species2 (like the excel file) for this to work. 

Of course in our case the genetic distances are between strains names by dates and they are not species. 

```{r,eval=F}
#calculate how many years between two groups
flu$time.passed = flu$species1 - flu$species1
```

<br/> 
<br/> 
We can now make a scatter plot, comparing genetic distance and time distance.

```{r,eval=F}
#see if distance increases with time, as mutations accumulate
plot(flu$time.passed,flu$distance,
     xlab ="time (years)",
     ylab ="genetic distance")

#add a line of best fit:
abline(lm(flu$distance ~ flu$time.passed),col=2,lty=2)
```

<br/> 
<br/> 
Finally, we can use a test to determine if the correlation is stattically significant. We will us the default test in R, Pearson's product moment correlation.

```{r,eval=F}
#test the correlation
cor.test(flu$time.passed,flu$distance)
```


<br/> 
<br/> 

**Congratulations, you have taken your first steps in R!**

***




##Examine the between-country differences

Now we will use R Studio to examine whether New Zealand (NZ) and Australia share the same population of influenza viruses.

We can test this by examining whether within-country genetic distances are any different from between-country genetic distances.

<br/> 
<span style="color:purple">**DISCUSS WITH YOUR GROUP**</span>

<span style="color:purple">What would we expect if NZ and Australia share influenza viruses?</span>

<span style="color:purple">What would we expect if they did **not** share influenza viruses?</span>
<br/> 

***

First, load the data into R.
```{r,eval=F}
#open the tasman data
tas <- read.table("tasman",h=T)
```
<br/> 

Then add country labels for each comparison.

```{r,eval=F}
#label regions, using grep
#set up AUS as a default for both regions
tas$country1 = "AU"
tas$country2 = "AU"

#now label the NZ rows, using grep
tas[grep("NZ", tas$region1),]$country1 = "NZ"
tas[grep("NZ", tas$region2),]$country2 = "NZ"
```
<br/> 

Now create new data frames called **within** and **between** from subsets of the initial **tas** data frame.
These contain all the within-country genetic distances, and all between-country genetic distances.

```{r,eval=F}
#get the within country data
within <- subset(tas, country1 == country2)
between <- subset(tas, country1 != country2)
```
</br>


Now we will amke a plot to show the within-country genetic distances, and all between-country genetic distances.
```{r,eval=F}
#make a plot to see if within-country groups are more closely relared than between-country groups
stripchart(
  list(within$distance,between$distance),
  vert=T,
  method="jitter",
  pch=1,
  group.names=c("within","between"),
  ylab="distance"
)
```
</br>


Does it look to you like within-country distances, are very different from between-country distances? We can test your intuition using a statistical test. Usually we could use Student's t-test to determine if the means of two sets of data are significantly different from each other. 

But in thia case, we will be careful and use  Wilcoxon rank sum test, which does not assume that the data are normally distrubuted. (This test is sometimes called a Mann-Whitney U test).


Here is how we run a Wilcoxon rank sum test:


```{r,eval=F}
wilcox.test(within$distance,between$distance)
```
</br>

***

<span style="color:purple">**DISCUSS WITH YOUR GROUP**</span>

<span style="color:purple">Is there a significant difference?</span>

<span style="color:purple">What does this mean?</span>


***

##The End
That is it for today. We hope you found the practical interesting and useful.

**REMEMBER: Store your script files and plots. You will need to present your results on Friday with your group.**

***
[Link back to main course page](http://www-users.york.ac.uk/~dj757/popgenomics/popgenomics_index.html)





