[Mattice-commits] r37 - pkg/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Nov 19 16:43:08 CET 2008
Author: andrew_hipp
Date: 2008-11-19 16:43:08 +0100 (Wed, 19 Nov 2008)
New Revision: 37
Modified:
pkg/R/regimes.R
pkg/R/treeTraversal.R
Log:
multitrees
Modified: pkg/R/regimes.R
===================================================================
--- pkg/R/regimes.R 2008-11-19 15:20:52 UTC (rev 36)
+++ pkg/R/regimes.R 2008-11-19 15:43:08 UTC (rev 37)
@@ -1,3 +1,14 @@
+# ---------------------------------------------------
+# FUNCTIONS FOR PAINTING REGIMES ON AN S4 OUCH TREE #
+# ---------------------------------------------------
+
+# Modified from functions used in Hipp 2007 Evolution paper
+# Initially written for ouch v 1.2-4
+# updated to ouch >= 2.4-2 Nov 2008
+
+# FINISH REGIME MAKER
+
+
paintBranches <-
# Paints branches with regimes changing at nodes specified
# arguments
@@ -117,8 +128,9 @@
regimeMaker <- function(ouchTrees, regMatrix, nodeMembers) {
## supplants the old 'allPossibleRegimes'
+## takes a list of ouchtree objects, a regimeMatrix ouput, and a list of nodeMembers (the taxa definining each node of interest)
## Value:
-## regList = a list of regimes for each tree (i.e., a list of lists)
+## regList = a list of nodes defining the change points for each tree (i.e., a list of lists)
## nodeMatrix = a matrix of trees (rows) by nodes (columns) indicating whether the node is present in each tree
# set up variables
@@ -126,15 +138,15 @@
numNodes <- length(nodeMembers)
if(numNodes != dim(regMatrix)[1] stop('Number of nodes (columns) in regMatrix must equal number of items in nodeMembers list')
nodeMatrix <- matrix(NA, nrow = numTrees, ncol = numNodes, dimnames = list(seq(numTrees), dimnames(regMatrix)[2]))
+ changeNodes <- list(numTrees)
regList <- list(numTrees)
- # fill nodeMatrix
+ # fill outData
for(i in seq(numNodes)) nodeMatrix[, i] <- unlist(lapply(ouchTrees, isMonophyletic, taxa = nodeMembers[[i]]))
-
- # fill regList
- for(i in seq(numTrees)) {
- regList <- XXX
- STOPPED HERE
+ for(i in seq(numTrees)) changeNodes[[i]] <- unlist(lapply(nodeMembers[as.logical(nodeMatrix[i, ], mrcaOUCH, tree = ouchTrees[[i]])]))
+ FILL UP regList .. should be a list (trees) of lists (regimes)
+ outdata <- list(regList = regList, nodeMatrix = nodeMatrix)
+ return(outdata)
}
@@ -178,10 +190,12 @@
# Generates the list of painted branches representing all possible selective regimes for OU analyses, taking as argument
# species vectors that describe the clades at the bases of which regimes are specified to change.
# Arguments:
-# "node" "ancestor" "times" "species" = the standard tree specification vectors of the OUCH-style tree
+# "tree" = the standard tree specification vectors of the OUCH-style tree
# "cladeMembersList" = list of vectors containing names of the members of each clade (except for the root of the tree)
# Value: list of vectors that can each be plugged directly into OU analysis as the "regimes" argument
-function(tree, cladeMembersList, maxNodes = NULL) {
+# 19 nov 08: changing to accept a list of trees
+
+function(ouchTrees, cladeMembersList, maxNodes = NULL) {
## ------------------ begin ouchtree block -----------------
## check to see if tree inherits 'ouchtree'
if (!is(tree,'ouchtree'))
@@ -199,14 +213,14 @@
#changeNodesVector = vector("character", length(changeNodesList))
#for (i in 1:length(changeNodesList)) # Changing cladeMemberList to a 1-d vector
# {changeNodesVector[i] = changeNodesList[[i]]}
- apr = allPossibleRegimes(changeNodesVector, maxNodes)
- allRegimes <- apr$regimeList
- regimeMatrix <- apr$regimeMatrix
+ regMatrix <- CALL REG MATRIX
+ apr = regimeMaker(xxx) ## HOLD IT! NOW REGIME MAKER WORKS ON ALL TREES AT ONCE... RETHINK THIS
+ allRegimes <- regimesList
regimePaintings = vector("list", length(allRegimes))
for (i in 1:length(allRegimes)) {
allRegimes[[i]] <- c("1", allRegimes[[i]])
regimePaintings[[i]] <- as.factor(paintBranches(tree, allRegimes[[i]], as.character(allRegimes[[i]])))
names(regimePaintings[[i]]) <- tree at nodes
message(paste('Created regime',i))}
- outdata <- list(regimeList = regimePaintings, regimeMatrix = regimeMatrix)
+ outdata <- list(regimeList = regimePaintings, regimeMatrix = regMatrix)
return(outdata) }
Modified: pkg/R/treeTraversal.R
===================================================================
--- pkg/R/treeTraversal.R 2008-11-19 15:20:52 UTC (rev 36)
+++ pkg/R/treeTraversal.R 2008-11-19 15:43:08 UTC (rev 37)
@@ -1,6 +1,6 @@
-# ---------------------------------------------------------------
-# FUNCTIONS FOR TRAVERSING AN S4 OUCH TREE AND PAINTING REGIMES #
-# ---------------------------------------------------------------
+# ------------------------------------------
+# FUNCTIONS FOR TRAVERSING AN S4 OUCH TREE #
+# ------------------------------------------
# Modified from functions used in Hipp 2007 Evolution paper
# Initially written for ouch v 1.2-4
@@ -13,67 +13,17 @@
# 4. allPossibleRegimes
# 5. regimeVectors
+# 19 nov 08: moved all regime-involved routines to regimes.R
+
# To do:
# 2. make allPossibleRegimes more efficient when maxNodes < length(nodes)
-paintBranches <-
-# Paints branches with regimes changing at nodes specified
-# arguments
-# "node" "ancestor" "times" = the standard tree specification vectors of the OUCH-style tree
-# "regimeShiftNodes" = a vector of nodes at which selective regimes shift: root must be included, but tips are meaningless in this context
-# "regimeTitles" = a vector of titles for the regimes that begin at the root and at the nodes indicated in "regimeShiftNodes",
-# in order of description in "regimeShiftNodes", except that the root is listed first in "regimeTitles"
-# but not at all in "regimeShiftNodes"... defaults to "node[x]regime
-# Value: a vector of regimes that can be plopped right into an OUCH-style tree data frame
-function(tree, regimeShiftNodes, regimeTitles) {
- ## ------------------ begin ouchtree block -----------------
- ## check to see if tree inherits 'ouchtree'
- if (!is(tree,'ouchtree'))
- stop(paste('This function has been rewritten to use the new S4 ', sQuote('ouchtree'), ' class.',
- '\nYou can generate a tree of this class by calling ', sQuote('ouchtree()'), '.', sep = ""))
- ## get the vectors we need:
- ancestor <- tree at ancestors # class = "character"
- node <- tree at nodes # class = "character"
- species <- tree at nodelabels # class = "character" -- note that nodelabels is more general than this indicates and the name should be changed throughout at some point
- times <- tree at times # class = "numeric"
- ## ------------------ end ouchtree block -------------------
-
- names(regimeTitles) = as.character(regimeShiftNodes)
- colorsVector = character(length(node))
- for (i in 1:length(ancestor)) {
- # First three lines fill up the vector for nodes that are hit in order
- if (is.na(ancestor[i])) {
- colorsVector[i] = regimeTitles["1"]
- next }
- if (as.character(ancestor[i]) %in% as.character(regimeShiftNodes)) {
- colorsVector[i] = regimeTitles[as.character(ancestor[i])]
- next }
- if (colorsVector[as.integer(ancestor[i])] != "") {
- colorsVector[i] = colorsVector[as.integer(ancestor[i])]
- next }
- # These lines fill up the vector for nodes run reached before their immediate ancestor
- nodeQ = integer(length(node))
- ii = i
- repeat {
- nodeQ = c(ii, nodeQ)
- ii = as.numeric(ancestor[ii])
- if (as.character(ancestor[ii]) %in% as.character(regimeShiftNodes)) {
- colorsVector[ii] = colorsVector[as.integer(ancestor[ii])]
- break}
- if (colorsVector[as.integer(ancestor[ii])] != "") {
- colorsVector[ii] = colorsVector[as.integer(ancestor[ii])]
- break} }
+isMonophyletic <- function(tree, taxa) {
+# returns T or F on whether a group of taxa is monophyletic in an ouch tree
+ identical(sort(taxa), sort(nodeDescendents(tree, mrcaOUCH(taxa, tree))))
+}
- for(j in nodeQ) {
- colorsVector[j] = colorsVector[as.integer(ancestor[j])] }
-
- } # closes for(i in 1:length(ancestor)) loop
-
- # a little hack to fix a problem I don't understand... with the undesired side effect that it colors the stem of some subtrees rather than the crown as originally written
- for(i in 1:length(colorsVector)) if(colorsVector[i] == "") colorsVector[i] <- as.character(i)
- return(colorsVector) }
-
nodeDescendents <- function(tree, startNode) {
## Recursive function to find all the descendents of a node on an 'ouchtree' object
startNode <- as.character(startNode) # just to be safe
@@ -85,11 +35,6 @@
return(nodeNames[!is.na(nodeNames)])
}
-isMonophyletic <- function(tree, taxa) {
-# returns T or F on whether a group of taxa is monophyletic in an ouch tree
- identical(sort(taxa), sort(nodeDescendents(tree, mrcaOUCH(taxa, tree))))
-}
-
mrcaOUCH <-
# Finds most recent common ancestor for a vector of tips by:
# 1. Creating a vector of ancestral nodes leading to each tip
@@ -132,6 +77,7 @@
# "tip" = the tip node to trace back
# Value: a vector of nodes leading from a tip to the root
# 10 nov 08: changed to just grab tree at lineages and make a vector that fits the old code
+# this should really be exluded at some point
function(tip, tree) {
## ------------------ begin ouchtree block -----------------
## check to see if tree inherits 'ouchtree'
@@ -154,147 +100,4 @@
# tip = ancestor[tip] }
nodesVector <- c(as.character(tree at lineages[[tip]][2:length(tree at lineages[[tip]])]), NA)
return(nodesVector)
- }
-
-allPossibleRegimes <-
-# Generates a list of vectors of all possible 2^n regimes for a given list of n ancestor nodes, assuming that each node
-# can either be a change node or not, and that all combinations are meaningful.
-# NOTE: This routine returns all possible permutations, ignoring the fact that the root must be included as a changeNode
-# for paintBranches to work properly. Exclude the root when calling this routine.
-# Arguments:
-# "changeNodes" = vector of all change nodes in all possible scenarios, or number of regimes assumed if nodeMatrix = T
-# "maxNodes" = single number indicating the maximum number of nodes at which a regime can change
-# "nodeMatrix" = indicates whether to return a binary table for interp or a list of changeNode vectors for analysis
-# Value:
-# regimeList = a list of changeNode vectors (assumes type = "character"), one for each possible scenario
-# regimeMatrix = : a binary table indicating whether a regime node is present or absent based on allPossibleRegimes output;
-# presumes nodes are labeled 1:n
-# 10 nov 08: this function now takes over regimeNodes
-function(changeNodes, maxNodes = NULL, nodeMatrix = F) {
- if(!identical(maxNodes, NULL) && maxNodes > length(changeNodes)) warning(paste(sQuote('maxNodes'), 'cannot be larger than the number of nodes; maxNodes ignored'))
- numberOfRegimes = ifelse(length(changeNodes) == 1, 2, 2^length(changeNodes))
- regime = vector("list", numberOfRegimes)
- for (i in 1:(numberOfRegimes - 1)) {
- remainder = i
- n = NULL
- for (j in as.integer(log2(i)):0) {
- if (2^j > remainder) n[j+1] = NA
- else {
- n[j+1] = changeNodes[j+1]
- remainder = remainder %% 2^j
- }
- }
- regime[[i]] = sort(n[!is.na(n)])
- }
- regime[[numberOfRegimes]] = rep("0", times = as.integer(log2(i)) + 1)
-
- ## make node matrix
- regimesNameMatrix = matrix(
- data = NA, nrow = numberOfRegimes, ncol = length(changeNodes), dimnames = list(
- as.character(seq(numberOfRegimes)), as.character(seq(length(changeNodes)))
- )
- )
- for (i in seq(numberOfRegimes)) {
- for (j in seq(length(changeNodes))) {
- if (is.na(match(changeNodes[j],regime[[i]]))) regimesNameMatrix[i,j] = 0 # changed this so that j indexes changeNodes
- else regimesNameMatrix[i,j] = 1
- }
- }
- outmatrix <- regimesNameMatrix
- if(!identical(maxNodes, NULL)) {
- outmatrix <- outmatrix[apply(outmatrix,1,sum) <= maxNodes, ]
- dimnames(outmatrix)[[1]] = as.character(seq(dim(outmatrix)[1]))
- }
-
- ## prune list
- outlist <- regime
- if(!identical(maxNodes, NULL)) {
- outlist <- outlist[sapply(outlist, length) <= maxNodes]
- outlist[[length(outlist) + 1]] <- rep("0", length(changeNodes))
- }
- # }
- outdata <- list(regimeList = outlist, regimeMatrix = outmatrix)
- return(outdata) }
-
-regimeMaker <- function(ouchTrees, regMatrix, nodeMembers) {
-## supplants the old 'allPossibleRegimes'
-## Value:
-## regList = a list of regimes for each tree (i.e., a list of lists)
-## nodeMatrix = a matrix of trees (rows) by nodes (columns) indicating whether the node is present in each tree
- nodeMatrix <- lapply(isMonophyletic
-
-}
-
-
-regimeMatrix <- function(n = NULL, nodeNames = NULL, regimeNames = NULL, maxNodes = NULL) {
- if(identical(n, NULL) && identical(nodeNames, NULL)) stop("You have to give regimeMatrix the number of nodes, a vector of node names, or both")
- if(identical(nodeNames, NULL)) nodeNames <- as.character(seq(n))
- else n <- length(nodeNames)
- numberOfRegimes <- ifelse(n == 1, 2, 2^n)
- outmatrix <- matrix(NA, nrow = numberOfRegimes, ncol = n, dimnames = list(regimeNames, nodeNames))
- for(i in 1:(numberOfRegimes - 1)) outmatrix[i, ] <- as.binary(i, digits = n)
- outmatrix[numberOfRegimes, ] <- as.binary(0, digits = n)
- if(!identical(maxNodes, NULL)) {
- outmatrix <- outmatrix[apply(outmatrix,1,sum) <= maxNodes, ]
- dimnames(outmatrix)[[1]] = as.character(seq(dim(outmatrix)[1]))
- }
- return(outmatrix)
-}
-
-as.binary <- function(n, base = 2, r = FALSE, digits = NULL)
-# Robin Hankin <initialDOTsurname at soc.soton.ac.uk (edit in obvious way; spam precaution)>
-# submitted to R listserv Thu Apr 15 12:27:39 CEST 2004
-# AH added 'digits' to make it work with regimeMatrix
-# https://stat.ethz.ch/pipermail/r-help/2004-April/049419.html
-
-{
- out <- NULL
- while(n > 0) {
- if(r) {
- out <- c(out , n%%base)
- } else {
- out <- c(n%%base , out)
- }
- n <- n %/% base
- }
- if(!identical(digits, NULL) && !r) out <- c(rep(0, digits-length(out)), out)
- if(!identical(digits, NULL) && r) out <- c(out, rep(0, digits-length(out)))
- return(out)
-}
-
-regimeVectors <-
-# Generates the list of painted branches representing all possible selective regimes for OU analyses, taking as argument
-# species vectors that describe the clades at the bases of which regimes are specified to change.
-# Arguments:
-# "node" "ancestor" "times" "species" = the standard tree specification vectors of the OUCH-style tree
-# "cladeMembersList" = list of vectors containing names of the members of each clade (except for the root of the tree)
-# Value: list of vectors that can each be plugged directly into OU analysis as the "regimes" argument
-function(tree, cladeMembersList, maxNodes = NULL) {
- ## ------------------ begin ouchtree block -----------------
- ## check to see if tree inherits 'ouchtree'
- if (!is(tree,'ouchtree'))
- stop(paste('This function has been rewritten to use the new S4 ', sQuote('ouchtree'), ' class.',
- '\nYou can generate a tree of this class by calling ', sQuote('ouchtree()'), '.', sep = ""))
- ## get the vectors we need:
- ancestor <- tree at ancestors # class = "character"
- node <- tree at nodes # class = "character"
- species <- tree at nodelabels # class = "character" -- note that nodelabels is more general than this indicates and the name should be changed throughout at some point
- times <- tree at times # class = "numeric"
- ## ------------------ end ouchtree block -------------------
-
- changeNodesList <- lapply(cladeMembersList, mrcaOUCH, tree = tree) #Returns a list of length-1 character vectors, each containing a single changeNode -- the fact that this is a list causes problems in paintBranches if not changed to a 1-d vector
- changeNodesVector <- unlist(changeNodesList)
- #changeNodesVector = vector("character", length(changeNodesList))
- #for (i in 1:length(changeNodesList)) # Changing cladeMemberList to a 1-d vector
- # {changeNodesVector[i] = changeNodesList[[i]]}
- apr = allPossibleRegimes(changeNodesVector, maxNodes)
- allRegimes <- apr$regimeList
- regimeMatrix <- apr$regimeMatrix
- regimePaintings = vector("list", length(allRegimes))
- for (i in 1:length(allRegimes)) {
- allRegimes[[i]] <- c("1", allRegimes[[i]])
- regimePaintings[[i]] <- as.factor(paintBranches(tree, allRegimes[[i]], as.character(allRegimes[[i]])))
- names(regimePaintings[[i]]) <- tree at nodes
- message(paste('Created regime',i))}
- outdata <- list(regimeList = regimePaintings, regimeMatrix = regimeMatrix)
- return(outdata) }
+ }
\ No newline at end of file
More information about the Mattice-commits
mailing list