A simple phylogeny workflow in R

I often have to do some routine phylogenetic analyses, where one has to find a 'good' model of substitution using some sort of model testing framework, then reconstruct a tree using maximum likelihood. This can all be done in R with the ape and phangorn libraries. Here, I'm going to reconstruct a phylogeny from 71 hepatitis C sequences of the core/E1 region from infected individuals in Egypt, studied by Ray et al. 2000. You can download the sequences here.

First of all, we need to load the libraries.

library(ape)
library(phangorn)

Next, we set the filename of a multiple sequence alignment in FASTA format.

myalignment.filename <- "ray2000.fas"
myalignment <- read.dna(myalignment.filename, format = "fasta", as.matrix = TRUE)

Use of the as.matrix=TRUE option means that an error will be thrown if sequences are of unequal length (which means they may not be aligned), and also allows an image of the alignment to be displayed.

image(myalignment, show.labels = FALSE)

In order to obtain a good-fitting model, I need a rough tree. I use neighbour joining (NJ), assuming a TN93 distance model (just because I like it). NJ takes a distance matrix as input, so I generate this matrix first.

myalignment.dist.tn93 <- dist.dna(myalignment, model = "TN93", as.matrix = TRUE)

It's always sensible to plot out the resulting distances.

hist(myalignment.dist.tn93[lower.tri(myalignment.dist.tn93)], xlab = "Distance", 
    ylab = "Frequency", main = "", col = "grey")

Neighbour joining is but a single command away.

myalignment.tn93.nj <- nj(myalignment.dist.tn93)

Here is a plot of the resulting tree.

plot(myalignment.tn93.nj, "unrooted", cex = 0.5, show.tip.label = FALSE)

In order to fit models using maximum likelihood in phangorn, I have to change the format of the data from ape's DNAbin format, to the phyDat structure used in phangorn. I can then run a selection of (nucleotide) substitution models.

myalignment.phydat <- as.phyDat(myalignment)  # Convert to a format modeltest understand
myalignment.modeltest <- modelTest(myalignment.phydat, tree = myalignment.tn93.nj, 
    model = c("JC", "F81", "K80", "HKY", "SYM", "GTR"), G = TRUE, I = TRUE, 
    k = 4, control = pml.control(epsilon = 1e-08, maxit = 3, trace = 1), multicore = FALSE)

## Warning: negative edges length changed to 0!

myalignment.modeltest

##      Model  df logLik   AIC   BIC
## 1       JC 139  -9094 18465 19024
## 2     JC+I 140  -8618 17516 18079
## 3     JC+G 140  -8368 17017 17579
## 4   JC+G+I 141  -8352 16985 17552
## 5      F81 142  -9006 18296 18866
## 6    F81+I 143  -8525 17337 17911
## 7    F81+G 143  -8252 16790 17364
## 8  F81+G+I 144  -8239 16766 17344
## 9      K80 140  -8483 17246 17809
## 10   K80+I 141  -7980 16243 16810
## 11   K80+G 141  -7686 15653 16220
## 12 K80+G+I 142  -7668 15619 16190
## 13     HKY 143  -8441 17167 17742
## 14   HKY+I 144  -7916 16120 16699
## 15   HKY+G 144  -7636 15561 16139
## 16 HKY+G+I 145  -7617 15524 16106
## 17     SYM 144  -8445 17179 17758
## 18   SYM+I 145  -7943 16177 16759
## 19   SYM+G 145  -7650 15589 16172
## 20 SYM+G+I 146  -7632 15555 16142
## 21     GTR 147  -8436 17166 17757
## 22   GTR+I 148  -7903 16103 16698
## 23   GTR+G 148  -7622 15541 16136
## 24 GTR+G+I 149  -7602 15502 16100

It can be difficult to spot the 'best' model, so I write a function that pulls out the row of the dataframe with the lowest AIC or BIC.

findBestModel <- function(mt, ic = "AIC") {
    if (ic == "AIC") {
        minic <- min(mt$AIC)
        bestmodel <- mt[mt$AIC == minic, ]
    } else if (ic == "BIC") {
        minic <- min(mt$BIC)
        bestmodel <- mt[mt$BIC == minic, ]
    }
    bestmodel
}
myalignment.bestmodel <- findBestModel(myalignment.modeltest)
myalignment.bestmodel

##      Model  df logLik   AIC   BIC
## 24 GTR+G+I 149  -7602 15502 16100

As a first step in getting a maximum likelihood tree, I generate a phylogenetic maximum likelihood fit using pml. Note that this doesn't optimise the parameters; to do that, I have to run optim.pml. In this example, I do this in several stages, to show how different aspects of the tree - the substitution model, the branch lengths, and the topology - can be optimised separately from one another.

myalignment.mltree <- pml(myalignment.tn93.nj, myalignment.phydat, model = as.character(myalignment.bestmodel$Model), 
    k = 4)

## Warning: negative edges length changed to 0!

# optimise model parameters without fitting edges
myalignment.mltree <- optim.pml(myalignment.mltree, optNni = FALSE, optBf = TRUE, 
    optQ = TRUE, optInv = TRUE, optGamma = TRUE, optEdge = FALSE, optRate = TRUE, 
    control = pml.control(epsilon = 1e-08, maxit = 10, trace = 0))
# Now fix substitution model and fit tree
myalignment.mltree <- optim.pml(myalignment.mltree, optNni = TRUE, optBf = FALSE, 
    optQ = FALSE, optInv = FALSE, optGamma = FALSE, optEdge = TRUE, control = pml.control(epsilon = 1e-08, 
        maxit = 10, trace = 0))
# Now fine tune
myalignment.mltree <- optim.pml(myalignment.mltree, optNni = TRUE, optBf = TRUE, 
    optQ = TRUE, optInv = TRUE, optGamma = TRUE, optEdge = TRUE, optRate = FALSE, 
    control = pml.control(epsilon = 1e-08, maxit = 10, trace = 0))

Finally, I display the tree and save it as a Newick-formatted string to a file.

myalignment.mltree

## 
##  loglikelihood: -7529 
## 
## unconstrained loglikelihood: -1991 
## Proportion of invariant sites: 0.2699 
## Discrete gamma model
## Number of rate categories: 4 
## Shape parameter: 0.9088 
## 
## Rate matrix:
##       a       c      g      t
## a 0.000  1.3620 8.0664  1.791
## c 1.362  0.0000 0.7302 11.374
## g 8.066  0.7302 0.0000  1.000
## t 1.791 11.3741 1.0000  0.000
## 
## Base frequencies:  
## 0.1749 0.3425 0.2533 0.2293

myalignment.mltree$tree

## 
## Phylogenetic tree with 71 tips and 69 internal nodes.
## 
## Tip labels:
##  AF271819, AF271820, AF271821, AF271823, AF271824, AF271822, ...
## 
## Unrooted; includes branch lengths.

plot(myalignment.mltree$tree, "unrooted", show.tip.label = FALSE)

write.tree(myalignment.mltree$tree, "mytree.tre")

An R markdown file with the workflow can be found here.

Comparing the performance of R and Rcpp for a simple epidemiological model

While R is great for statistics, like many dynamic languages, it's terribly slow in loops. This can make dynamical models, like those commonly used in epidemiology, rather slow. However, thanks to the wonderful Rcpp package, such simulations can be sped up considerably. Let's take the 'standard' susceptible-infected-recovered (SIR) model. A deterministic version of this model is as follows.

\[\begin{align} \frac{dS(t)}{dt} & = -\beta S(t) I(t) \cr \frac{dI(t)}{dt} & = \beta S(t) I(t)- \gamma I(t) \cr \frac{dR(t)}{dt} & = \gamma I(t) \end{align}\]
Let's consider a stochastic version of the SIR model.

\[ \begin{align} {\rm Transition} & \quad {\rm Rate} \cr S \rightarrow S-1,\; I \rightarrow I+1 & \quad \beta S(t) I(t) \cr I \rightarrow I-1,\; R \rightarrow R+1 & \quad \gamma I(t) \end{align}\]
We can simulate this model using the Gillespie method, also known as the stochastic simulation algorithm, or SSA. In plain R, the stochastic model is as follows.

sir <- function(beta, gamma, N, S0, I0, R0, tf) {
    time <- 0
    S <- S0
    I <- I0
    R <- R0
    ta <- numeric(0)
    Sa <- numeric(0)
    Ia <- numeric(0)
    Ra <- numeric(0)
    while (time < tf) {
        ta <- c(ta, time)
        Sa <- c(Sa, S)
        Ia <- c(Ia, I)
        Ra <- c(Ra, R)
        pf1 <- beta * S * I
        pf2 <- gamma * I
        pf <- pf1 + pf2
        dt <- rexp(1, rate = pf)
        time <- time + dt
        if (time > tf) {
            break
        }
        ru <- runif(1)
        if (ru < (pf1/pf)) {
            S <- S - 1
            I <- I + 1
        } else {
            I <- I - 1
            R <- R + 1
        }
        if (I == 0) {
            break
        }
    }
    results <- data.frame(time = ta, S = Sa, I = Ia, R = Ra)
    # print(head(results))
    return(results)
}

Here is the version written using Rcpp, inlining the C++ code.

library(Rcpp)
cppFunction('
  List sirc(double beta, double gamma, double N, double S0, double I0, double R0, double tf){
    double t = 0;
    double S = S0;
    double I = I0;
    double R = R0;
    std::vector<double> ta;
    std::vector<double> Sa;
    std::vector<double> Ia;
    std::vector<double> Ra;
    do{
      ta.push_back(t);
      Sa.push_back(S);
      Ia.push_back(I);
      Ra.push_back(R);
      double pf1 = beta*S*I;
      double pf2 = gamma*I;
      double pf = pf1+pf2;
      double dt = rexp(1,pf)[0];
      t += dt;
      double r = runif(1)[0];
      if(r<pf1/pf){
        S--;
        I++;
      }else{
        I--;
        R++;
      }
      if(I==0){break;}
    } while (t<=tf && (I>0));
  return List::create(_["time"] = ta, _["S"] = Sa, _["I"] = Ia, _["R"]=Ra);
  }'
)

Now we run both the vanilla R and Rcpp versions of the same model 100 times, with parameter values \( \beta \) = 0.1/10000 and \( \gamma \) = 0.05, initial conditions \( S(0) \) =9999, \( I(0) \)=1, \( R(0) \)=0, and a simulation time of 1000.

N <- 100  # Run 100 times
set.seed(4)
sirsim <- sum(replicate(N, system.time(x <- sir(0.1/10000, 0.05, 10000, 9999, 
    1, 0, 1000))["elapsed"]), trim = 0.05)
set.seed(4)
sircsim <- sum(replicate(N, system.time(x <- sirc(0.1/10000, 0.05, 10000, 9999, 
    1, 0, 1000))["elapsed"]), trim = 0.05)

As you can see, the Rcpp version is considerably faster.

paste("Plain R:", sirsim)

## [1] "Plain R: 114.178"

paste("Rcpp version:", sircsim)

## [1] "Rcpp version: 0.175000000000028"

paste("R took", sirsim/sircsim, "times longer")

## [1] "R took 652.445714285608 times longer"