# phylo-ML.r # T. J. Finney, 2014-01-02. message("Maximum likelihood phylogenetic analysis...") # Set parameters. message("Input data matrix: ", input.file <- "../data/clean/Heb-UBS4.csv") message("Output directory: ", output.path <- "../phylo/ML/") # Set writing mode: [1] = write results; [2] = do not write results. message("Write mode: ", write <- c(TRUE, FALSE)[2]) # Read input. #input.frame <- read.csv(input.file, row.names=1, colClasses="character") # Do analysis. See vignette("Trees"), section 6. require(phangorn) require(parallel) data(Laurasiatherian) phy = Laurasiatherian dm = dist.ml(phy) # Alternative starting trees. treeUPGMA = upgma(dm) treeNJ = NJ(dm) message("Do ratchet...") treeRatchet <- pratchet(phy, trace=1) message("Do rearrangement...") treeRe = optim.parsimony(treeUPGMA, phy) #treeRe = optim.parsimony(treeNJ, phy) # Parsimony scores. pUPGMA = parsimony(treeUPGMA, phy) pNJ = parsimony(treeNJ, phy) pRatchet = parsimony(treeRatchet, phy) pRe = parsimony(treeRe, phy) message("Parsimony scores:") message("UPGMA: ", pUPGMA) message("NJ: ", pNJ) message("Ratchet: ", pRatchet) message("Rearranged : ", pRe) # Branch and bound algorithm (small datasets only). message("Do branch and bound...") (trees <- bab(subset(phy, 1:10))) # Maximum likelihood. message("Do maximum likelihood...") fit = pml(treeNJ, data=phy) fitJC = optim.pml(fit, TRUE) fitGTR = update(fit, k=4, inv=0.2) fitGTR = optim.pml(fitGTR, TRUE, TRUE, TRUE, TRUE, TRUE, control = pml.control(trace = 0)) # Try different models... #mt = modelTest(phy) #env <- attr(mt, "env") message("Do bootstrap...") #bs = bootstrap.pml(fitJC, bs=100, optNni=TRUE, control = pml.control(trace = 0)) # Do plots. par(bg="white") par(mar = c(.2,.1,.1,.1)) layout(matrix(1:6, 2, 3)) plot(treeUPGMA, main="UPGMA") plot(treeNJ, "unrooted", main="NJ") plot(treeRatchet, main="Ratchet") plot(treeRe, main="Rearranged") plot(fitJC, main="ML: JC") plot(fitGTR, main="ML: GTR") # Clean up. #rm(list=ls(all=TRUE))