Once the splits package was installed, you are ready to run the gmyc to delimit species.
First, load the package, and load your tree,
>tr <- read.tree(“./youtree.tre”)
The read.tree function reads a tree file written in Newick format. If your tree is in Nexus format, use read.nexus instead of read.tree.
Your tree must be ultrametric and fully bifurcating. A tree with polytomy can be easily converted to a dichotomous tree by the multi2di function. Making ultrametric trees can be done by any dating programs like BEAST or r8s. The method using the relaxed clock is better. However, the GMYC is not hyper-sensitive to the dating methods. Even dating with the molecular clock often returns the result similar to the relaxed clock (for example, see Monaghan et al. 2009).
Calling gmyc runs the gmyc analysis. The default setting is the single threshold method.
>result <- gmyc(tr)
You will see outputs like below after a few seconds (or minutes/hours, depending on tree size).
node T loglik
2 -0.954603 863.3993
3 -0.9010095 864.1611
148 -0.00011231 861.0263
149 -4.374e-05 860.83
Mon Apr 22 12:45:51 2013
The analysis is done. Now, you can check the number of delimited species and the significance of the GMYC model by the summary function.
Result of GMYC species delimitation
likelihood of null model: 858.9326
maximum likelihood of GMYC model: 872.3755
likelihood ratio: 26.88581
result of LR test: 6.220952e-06***
number of ML clusters: 33
confidence interval: 31-34
number of ML entities: 46
confidence interval: 42-56
threshold time: -0.05174417
The “number of ML clusters” is the number of delimited “clusters”, groups including more than 2 samples. The “number of ML entities” is the number of all delimited species including singletons (species represented by a single sample). The GMYC is usually robust to a moderate number of singletons (10-20% of singletons or more). Rare species will be reasonably delimited.
You can plot a lineage through time plot, a likelihood surface and a tree with delimited species by calling plot.
The red vertical line indicates the threshold time between inter – intraspecific branching. The red branches on the gene tree are delimited species.
If you want detailed associations of samples to species, use spec.list.
1 1 spec1.5
2 1 spec1.4
3 1 spec1.3
4 1 spec1.2
This function returns a dataframe showing the sample-species associations. You may save this result of delimitation by calling the write.table function and use it for other analyses. (spec.list is available in the splits >= 1.0-15.)