Getting started

R

We assume you have some basic understanding of programming and R. You will not have to write a lot of code yourself, but you should be able to read and run the scripts that we provide you with. If you have never seen any R code before, it might be useful for you to work through a tutorial, such as the first three pages of this tutorial by Bart de Boer. Of course, you can always ask us for help if you are stuck on a particular piece of code!

RStudio

We strongly recommend you use a particular environment for running your R code, called RStudio. If you haven’t installed RStudio yet, you should do so by finding the appropriate download of RStudio Desktop (Open Source License) for your operating system on rstudio.com and following the instructions until you can open the program.

Console & Scripts

When you open RStudio, you will find your screen divided into a few different parts. The most important ones that allow you to interact with R are on the left. The bottom left panel, called the Console, is where you can enter commands directly, and where the results of your command will be shown immediately after pressing Enter. When R is ready to accept commands, the Console shows a > prompt. In the Console, you can retrieve commands that you entered previously by using the up arrow key. However, when you close RStudio, all the commands you entered in the Console will be forgotten. When you are using a sequence of multiple commands and/or would like to run your code more often than once, it is better to save them in a Script (a plain text file that contains your code). The top left panel in your RStudio window is the Script editor. Here you can enter commands that will not disappear after you run them (you can run a selected line of code from the Script editor by pressing Ctrl-Enter). The tab at the top of the panel will show the name of your script. R script filenames should end in .R. For the commands given in the labs of this course, we will always tell you if you should enter them in the Console or save them to your Script.

Working directory

To run a full script in R, you type source("scriptname") in the Console. To be able to run a script, it should be available to R. In RStudio, the easiest way to do this is to put them both in the same folder and set this folder as the working directory. In the labs for this course, you should always set your working directory to the folder with materials we provide for that lab (for this lab, that folder is lab-evolution-reconstruction). Look at this website to find out how to set the working directory for your version of RStudio.

Simulated Evolution

The goals of this section are:

• to better understand the concepts of genotype, genotype space, fitness, fitness landscape, selection, mutation, selection-mutation balance, frequency dependent selection;

• to see how these concepts can be formalised in a computer program;

• to appreciate both the power and the limits of natural selection.

You need the following R files in the folder: and auxiliary_functions.R and simulated_evolution.R.

Creating a population

In the first part of this computer lab, we will use R to simulate the evolution of (DNA) strings under a particular fitness function. We will represent DNA strings (the genotype) as a sequence of the letters ‘A’, ‘G’, ‘C’ and ‘U’. In R, you can generate a random sequence of 10 of these letters using the following command (try it out in the Console):

sample(c('A','G','C','U'), size=10, replace=TRUE)

To store the output in a variable (for instance x), you type:

x <- sample(c('A','G','C','U'), size=10, replace=TRUE)

You can then view the contents of a particular variable by simply typing its name (here, x) and pressing Enter. Let’s experiment a bit with this.

Now, let’s create a population of DNA strings. To do this, we will make 100 genotype strings. We will store our population in a matrix (the population matrix), where each member of our population is represented as a row of the matrix.

population <- matrix(rep(0, 100), 100, 50))

population_size <- 100
for (i in 1:population_size) {
    population[i,] <- sample(c('A','G','C','U'),
      size=50, replace=TRUE)
}

Defining a fitness function

Now we need to define a fitness function that computes the fitness of the individual members of our population. Imagine, for instance, that the string ‘CAC’ codes for some very useful aminoacid. The more CAC’s in the genome, the higher the expected number of offspring. In our simulation of evolution, let’s define the fitness as the number of times the substring ‘CAC’ appears in the genotype string (without overlap, so the string ‘CACAC’ contains one copy of ‘CAC’).

To keep track of the fitness of all members in our population (which are represented as rows in the population matrix), we create a vector where each element of the vector represents the fitness value for one member of the population.

fitness <- rep(0, population_size)

# loop over population size
for (i in 1:population_size) {      
  # generate string representation
  member <- paste(population[i,], collapse='') 
  
  # compute fitness member   
  fitness_member <- str_count(member, "CAC")   
  
  # store in fitness vector
  fitness[i] <- fitness_member
}

Implementing selection

Now we will generate the next generation. We assume that each member of the next generation inherits the genome of one of the members of the previous generation. The probability of inheriting each genome is proportional to the genome’s fitness: a child is most likely to inherit the genome of the fittest member of the previous population. This simulates selection.

Evolution with mutation

In the previous simulation, we looked at selection without mutation. Let’s now look at the case where every child’s nucleotide has a probability $\mu$ to change into a random other nucleotide.

Simulated evolution with ancestry

In the previous section we implemented computer simulations of evolution in R. In this section, we extend this simulation, and think about the patterns of genetic variation that the evolutionary process leaves in a population. We look at methods for phylogenetic tree reconstruction. These methods use the genetic variation in the most recent (current) generation to reconstruct the evolutionary history of a population or a set of species. We use a simple clustering algorithm⁴ for phylogenetic tree reconstruction.

The goal of this section is to learn about the possibilities and difficulties that are involved in phylogenetic tree reconstruction.

You need these two R files in the lab-evolution-reconstruction folder: auxiliary_functions.R and simulated_evolution_ancestry.R. To plot the phylogenetic trees with R, you will also need to have the package ggtree. Install it by typing the following into the Console:

source("https://bioconductor.org/biocLite.R")
biocLite("ggtree")

Tracking ancestry

So far, we have simulated the evolution of strings of symbols, and looked at the effect of different fitness functions. Now we will repeat this simulation, but during the evolutionary process we will keep track of ancestry, so that we can reconstruct family trees of different individuals. We will start with a very simple simulation. The script simulated_evolution_ancestry.R runs the same simulation we saw in the previous part of this lab. This time, however, the script also generates a matrix called parent_matrix that specifies the parent of each member in each generation (where the parent is the individual of the previous generation whose genetic material was inherited). At the end of the simulation, a plot illustrating the development of both the average population fitness and the diversity of the population is generated.

source("simulated_evolution_ancestry.R")

population_size <- 100
simulation_length <- 1000

Plotting the family tree

Printing the parent matrix for such large simulations is not very helpful, because the network is too dense to properly visualise (you may try if you want). Rather than looking the parent matrix, we will use the parent matrix to reconstruct a family tree for only the last generation (that is, we only look at the members of previous generations whose offspring appears in the last generation.) We will then generate a visual representation of the tree.

tree <- reconstruct_tree(parent_matrix)
plot_tree(tree)

Phylogenetic reconstruction of biological data

In evolutionary research, the elaborate ancestry information represented by the parent matrix is usually not available. To reconstruct family trees we have to resort to different methods. Information about when species branched off (‘speciated’) can be deduced from genetic variation in the current population. For instance, horses are genetically more similar to donkeys than to, say, frogs. So, the last common ancestor of horses and frogs most likely lived much further in the past than the last common ancestor of horses and donkeys. In other words, the branch that would eventually evolve into frogs split off from the branch that would eventually evolve into horses earlier than the branch that would eventually evolve into donkeys. This type of analysis is called phylogenetic reconstruction. It’s based on genetic similarity between members of the current generation, which we measure using a distance measure. There are R packages that can automatically perform this reconstruction. Lets start with installing these packages:

install.packages("ape")
install.packages("phangorn")

library(ape)
library(phangorn)

The phangorn package comes with a dataset that contains real genetic data (i.e., RNA samples) from many different species. You can load this dataset by typing the following (add it to your Script):

data(Laurasiatherian)

To show a summary of the data you can type str(Laurasiatherian) in the Console. The data originates from the (now closed) Allan Wilson Centre in New Zealand. To find out more about this data, have a look at their website.

We will try to reconstruct a phylogenetic tree for these species. That is, we will try to reconstruct when different species branched off from each other, based only on genetic information of the current population (the last generation). The first step is to measure ‘genetic distance’ between the genetic samples for each species. For simplicity, we assume that (1) all species ultimately originate from a single common ancestor (an uncontroversial assumption in evolutionary biology), and (2) that species have diverged genetically by picking up mutations at a roughly constant rate (a more problematic assumption). (Try and convince yourself that the phylogenetic tree reconstruction method we described requires the second assumption — and that this assumption is problematic when considering evolution in the real world.)

The distance between strings of DNA or RNA is typically measured by counting the number of mutations required to change one into the other. Because of the second assumption, the genetic distance between two species is proportional to the time that has passed since their last common ancestor.

mysubset <- subset(Laurasiatherian, subset=c(19,20,28,29,30)

str(mysubset) 

distance_matrix <- dist.ml(mysubset)
print(distance_matrix)

Phylogenetic trees with hierarchical clustering

UPGMA

We can use a simple method called ‘hierarchical clustering’ to build the phylogenetic trees such as in the previous sections automatically. Hierarchical clustering can be done with the UPGMA algorithm. (If you need a refresher on how exactly it worked, look up the UPGMA chapter from the readings again!) In short, UPGMA follows these steps:

• Treat each datapoint (for example, a RNA sample of a species) as a separate “cluster” containing just one datapoint;
• Compute the distances between all clusters (using some distance measure; for example, genetic distance);
• Merge the two clusters that are nearest to each other into a new cluster;
• Repeat steps 2 and 3 until only all datapoints are in cluster.

Constructing a tree

To construct a phylogenetic tree, we can think of each merging of clusters as the joining of two branches. In the simplest version of this algorithm, we define ‘distance’ between a cluster A and a cluster B as the average distance between any datapoint in A and any datapoint in B — this is what UPGMA does. A slightly more complicated method, Ward’s clustering, uses the square root of the average of the squared point-to-point distances.

The phangorn package we installed earlier provides pre-defined functions implementing different hierarchical clustering methods.

tree <- upgma(distance_matrix, method='average')
plot(tree)

We will now investigate what happens if we perform phylogenetic analysis on the population resulting from our own simulated evolution. Remember that, since this is a simulation over which we have full control, we can reconstructed the actual phylogenetic tree using the information that we stored in the parent matrix.

source("simulated_evolution_ancestry.R")

distance_matrix <- compute_distance_matrix(population)

gold_standard_tree <- reconstruct_tree(parent_matrix)
plot_tree(gold_standard_tree)

Phylogenetic reconstruction of languages

In the tutorial you have learned about features of language and music used in comparative research, and about variation across languages and musics. By studying some specific examples, you have seen that both languages and musical traditions are transmitted culturally, and are subject to a process of cultural evolution. In this final section of the computer lab, you will use the methods that we saw for phylogenetic tree reconstruction in the previous sections to reconstruct the cultural evolutionary history of languages.

You will need the file language_data.Rdata from the lab-evolution-reconstruction folder, and you will also need to have the packages ape and phangorn installed.

load('language_data.Rdata')

library(ape)
library(phangorn)

mysubset <- subset(mydata, c(40:42,58,60))

distance_matrix <- dist.hamming(mysubset)

tree <- upgma(distance_matrix, method='ward.D')

plot(tree, use.edge.length=FALSE, cex=2)