Simulating sequences and phylogenetic trees

Using the seqgen() function available in the phyclust package [1], data were simulated by mutating DNA sequences along phylogenetic trees. The topology of the trees were generated at two stages via the ms program [1, 2]. At the first stage we randomly generated a root sequence and a star-like ancestor tree having equal branch lengths and tips equal to the number of desired clusters. The tip label ag denotes the ancestor node, or recent common ancestor, of the gth cluster. The second stage generates a descendant tree for each cluster which is rooted at ag. The number of tips in the gth descendant tree represents the number of members in cluster g. The complete tree is obtained by attaching the G descendant trees to the tips of the ancestor tree. For instance, Figure S1 shows the complete tree (right panel) created by binding six randomly generating descendent trees (example in middle panel) to the corresponding tips on the ancestor tree (left panel). Standard phylogenetic trees typically involve bifurcating branches, i.e., exactly two diverging sequences per node. It may be the case, however, that the rapid diversification of HIV is better modeled using multifurcating trees [3]. For this reason, our simulations were created using star-like ancestor trees which take on a “pitchfork” structure (cf. Figure S1). In our context, this topology indicates that patients have simultaneously (or very rapidly) diverged from a common root [4]. The height of the ancestor tree, denoted by rA, determines the diversity between the ancestor node sequences, while the height of the descendant trees, rD, determines the diversity between sequences belonging to the same group. The proportion of the ancestor to descendent height is determined by the growth rate ratio rAD = rA/(rA + rD). Herein, complete trees are scaled to have total height one; consequently, rAD = rA. Based on the topology of the tree, we generated DNA sequences with seqgen(). Sequences were mutated according to a GTR model which assumes that sites evolve according to random variable R ∼ Γ(α, β) (parameterized such that E[R] = α/β). To avoid too many parameters, it is commonplace to assume a Gamma distribution parameterized by a single shape parameter α, where α = β. In addition, a proportion of sites (p) are assumed to be invariable. Together, these specifications are referred to as the GTR + I + Γ model. In accordance with [5] we set these parameters to values typical for HIV-1 Subtype B on the pol gene; namely, we set α = 0.7589, the proportion of invariant sites to p = 0.4817 and specify the rate matrix to: