phyclust (version 0.1-9)

ms: Generating Samples under a Wright-Fisher Neutral Model of Genetic Variation

Description

This function modifies the original standalone code of ms developed by Hudson (2002) for generating samples/coalescent trees under a Wright-Fisher neutral model.

Usage

ms(nsam = NULL, nreps = 1, opts = NULL)

Arguments

nsam

number of samples/coalescent trees, usually greater than 2.

nreps

number of replications.

opts

options as the standalone version

Value

• This function returns a vector, and each element stores one line of STDOUT of ms separated by newline, and store as a class ms. The details of output format can found on the website http://home.uchicago.edu/~rhudson1/source.html and its manual.

Details

This function directly reuses the C code of ms by arguments as input from the STDIN. The options opts is followed from the original ms except nsam and nreps.

For examples, options commonly used in phyclust are:

• "-T": generate trees in a neutral model.
• "-G": generate trees with a population growth rate, e.g. "-G 0.5".

These will return trees in a NEWICK format which can be read by the read.tree of ape and passed to seqgen to generate sequences.

References

Phylogenetic Clustering Website: http://thirteen-01.stat.iastate.edu/snoweye/phyclust/

Hudson, R.R. (2002) Generating Samples under a Wright-Fisher Neutral Model of Genetic Variation, Bioinformatics, 18, 337-338. http://home.uchicago.edu/~rhudson1/source.html

See Also

print.ms, read.tree, bind.tree, seqgen.

Examples

ms()

# an ancestral tree
set.seed(1234)
(ret.ms <- ms(nsam = 3, opts = "-T -G 0.1"))
(tree.anc <- read.tree(text = ret.ms[3]))
tree.anc$tip.label <- paste("a", 1:K, sep = "")

# adjacent descendant trees to the ancestral tree
K <- 3
N <- 12
N.k <- c(3, 4, 5)
ms.dec <- NULL         # a list to store trees of ms
tree.dec <- NULL       # a list to store the trees in phylo class
tree.joint <- tree.anc
for(k in 1:K){
  ms.dec[[k]] <- ms(N.k[k], opts = "-T -G 1.0")
  tree.dec[[k]] <- read.tree(text = ms.dec[[k]][3])
  tree.dec[[k]]$tip.label <- paste("d", k, ".", 1:N.k[k], sep = "")
  tree.joint <- bind.tree(tree.joint, tree.dec[[k]],
                          where = which(tree.joint$tip.label ==
                                        paste("a", k, sep = "")))
}
str(tree.joint)

# plot trees
par(mfrow = c(2, 3))
plot(tree.anc, main = paste("anc (", K, ")", sep = ""))
axis(1)
for(k in 1:K){
  plot(tree.dec[[k]], main = paste("dec", k, " (", N.k[k], ")", sep = ""))
  axis(1)
}
plot(tree.joint, main = paste("joint (", N, ")", sep = ""))
axis(1)