[R-sig-phylo] correlation of binary and multistate
Hi all, Does anybody know of a package that does character correlations between binary and multistate characters (and even better with the possibility of missing data)? Thanks! ~John ___ R-sig-phylo mailing list - R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
[R-sig-phylo] correlation of characters from stochastic maps
Hi all, I’m curious to know whether any packages implement character correlations and tests for significance of character correlation using the outputs of stochastic character mapping, such as make.simmap(), allowing for for characters with different numbers of states. I know the package correlate does pairwise correlations from stochastic maps, but as far as I can tell does not do significance testing. Thanks! ~John ___ R-sig-phylo mailing list - R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
[R-sig-phylo] maximum agreement subtrees?
Hi all, Is there a function in an existing package for calculating maximum agreement subtrees (MASTs)? Thanks! ~John John S. S. Denton, Ph.D. Department of Vertebrate Paleontology American Museum of Natural History www.johnssdenton.com [[alternative HTML version deleted]] ___ R-sig-phylo mailing list - R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
[R-sig-phylo] evolve discrete characters on non-ultrametric tree?
Hi all, I'm curious to know whether there is a function for evolving discrete character data onto a non-ultrametric tree. I know sim.char() can technically produce output given a non-ultrametric tree, and the documentation does not specify that it is a requirement, but all the example trees I have seen applied with this function are ultrametric. What are the consequences of using sim.char() or similar functions on a tree generated by something like rtree(), rather than rcoal() or sim.bdtree()? Are there ways to make equivalencies, as through some kind of regression or transformation? Thanks! ~John John S. S. Denton, Ph.D. Department of Vertebrate Paleontology American Museum of Natural History www.johnssdenton.com [[alternative HTML version deleted]] ___ R-sig-phylo mailing list - R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
[R-sig-phylo] extracting from mrbayes
Hi folks, I need to do some simulation using the Mkv model, and am wondering if there is a way to extract Mk Q-matrix values from mrbayes. A look at report lists revmat as an output, but I do not see any parameters in the .p files when I finish the run. Thanks, ~John John S. S. Denton, Ph.D. Department of Vertebrate Paleontology American Museum of Natural History www.johnssdenton.com ___ R-sig-phylo mailing list - R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
[R-sig-phylo] mcmcglmm with different classes of predictors
Hi all,
I'm trying to analyze a dataset that has left-censored, right-censored, and
categorical predictor variables and a univariate response in mcmcglmm. However,
I am uncertain of how to specify the appropriate model command and priors.
The response (rates) is a net diversification rate derived from BAMM. The
right-censored data is body size (max.min, max=Inf), and size at sexual
maturity is left censored (f.min= -Inf, f.max). There are four binary variables
indicating presence/absence: M.sup, M.inf, F.sup, F.inf.
My code is as follows:
library(MCMCglmm)
dat <- read.table("appended3.txt", row.names=1, header=TRUE) ##Read data table
dat$f.min[which(is.na(dat$f.min))]<- -Inf ##Assign -Inf for left-censored
information
t <- read.tree("MCC_pruned.tre") ##Read phylogenetic tree
tnames <- t$tip.label[order(t$tip.label)] ##Get tip names
full <- as.data.frame(cbind(tnames, dat)) ##Bind tip names to existing data
frame
invT <- inverseA(t)$Ainv##Calculate inverse, assuming nodes="ALL", the
default
prior <- list(R=list(V=1,nu=0.002),
G=list(G1=list(V=1,
nu=0.002))) ##General
"get it running" prior I have seen in tutorials
model_test <- MCMCglmm(rates ~ cbind(max.min, max) + cbind(f.min, f.max)+ M.sup
+ M.inf + F.sup + F.inf,
random = ~tnames, prior=prior,
family=c("gaussian","cenexponential","cenexponential", "threshold","threshold",
"threshold","threshold"), ginverse=list(tnames=invT),
data=full, nitt=500,burnin=1000,thin=500)
##The above specification is attempting to assign gaussian for rates,
cenexponential for the left-and right-censored data, and threshold for the four
binary characters.
I get the error
Error in matrix(unlist(value, recursive = FALSE, use.names = FALSE), nrow = nr,
:
'data' must be of a vector type, was 'NULL'
traceback() gives the following:
3: matrix(unlist(value, recursive = FALSE, use.names = FALSE), nrow = nr,
dimnames = list(rn, cn))
2: Ops.data.frame(data[, which(names(data) == response.names[nt +
1])], data[, which(names(data) == response.names[nt])])
1: MCMCglmm(rates ~ cbind(max.min, max) + cbind(f.min, f.max) +
M.sup + M.inf + F.sup + F.inf, random = ~tnames, prior = prior,
family = c("cenexponential", "cenexponential", "threshold",
"threshold", "threshold", "threshold"), ginverse = list(tnames =
invT),
data = full, nitt = 5e+06, burnin = 1000, thin = 500)
and debug() stops after
debug: if (any(data[, which(names(data) == response.names[nt + 1])] <
data[, which(names(data) == response.names[nt])], na.rm = T)) {
stop("for censored traits left censoring point must be less than right
censoring point")
}
However, the lower bound for max is equal to or larger than the upper bound for
f.min.
How can I get this running?
Many thanks,
~John
John S. S. Denton, Ph.D.
Department of Vertebrate Paleontology
American Museum of Natural History
www.johnssdenton.com
___
R-sig-phylo mailing list - R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
[R-sig-phylo] add geoscale to histogram or density plot
Hi all, I'd like to add a geoscale to a histogram/density plot, in a manner similar to axisGeo() or add.geoscale() in phyloch. Is there a function that can do this? Thanks! ~John John S. S. Denton, Ph.D. Department of Vertebrate Paleontology American Museum of Natural History www.johnssdenton.com ___ R-sig-phylo mailing list - R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
[R-sig-phylo] generate liability/threshold plot from ancThresh
Hi folks, I'm trying to generate the thresholds and liabilities plot, including the posterior support values, from Revell (2014) Evolution 68(3), Figure 6B. Does anybody have a script for this? I see the liabilities output from ancThresh(), but no thresholds. Thanks! ~John John S. S. Denton, Ph.D. Department of Vertebrate Paleontology American Museum of Natural History www.johnssdenton.com ___ R-sig-phylo mailing list - R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
[R-sig-phylo] drop trees from multiphylo list
Hi folks,
I have a set of trees with tip states, generated by simulation, in a multiphylo
object. I'd like to drop trees from the list based on a criterion (some trees
do not exhibit all three tip states).
I've tried the following (which is admittedly clumsy; I am not very familiar
with multiPhylo lists):
trees.pool - apply(pars, 1, function(pars) trees(pars, type=musse, n=1,
max.taxa=250, max.t=40, include.extinct=FALSE, x0=1))
class(trees.pool) - multiPhylo
# Screen these trees for ntax size, all states present:
for (i in 1:length(trees.pool)) {
if (length(unique(trees.pool[[i]]$tip.state)) 3) {
trees.pool[[i]] - NULL
} else {
i = i+1
}
}
but I get the error
Error in `[[-.multiPhylo`(`*tmp*`, i, value = NULL) :
trying to assign an object not of class phylo into an object of class
multiPhylo.
Any ideas?
Thanks!
~John
John S. S. Denton, Ph.D.
Department of Vertebrate Paleontology
American Museum of Natural History
www.johnssdenton.com
___
R-sig-phylo mailing list - R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
[R-sig-phylo] FW: tps format problem reading in R
Hi folks, I'm trying to read a tps file in the geomorph v1.1-1 R package using readland.tps(file), but every time I try the above, I get the error Error in dimnames(coords)[[3]] - imageID : 'dimnames' must be a list In addition: Warning message: In readland.tps(Lter_mean.TPS) : NAs introduced by coercion I do not get the error when I read in some of my other tps files (all of which were produced using append files in tpsUtil). The file I'm trying to read in is the side-averaged Procrustes coordinates, generated in MorphoJ. I've checked the hidden formatting, the line breaks, the number of decimal places, and the related file image extension, but I can't seem to figure it out. The file is attached. ~John John S. S. Denton Ph.D. Candidate Department of Ichthyology and Richard Gilder Graduate School American Museum of Natural History www.johnssdenton.com Lter_mean.TPS Description: Lter_mean.TPS ___ R-sig-phylo mailing list - R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
Re: [R-sig-phylo] collapse descendants of a node to a polytomy
Hi Liam and David,
Thanks for the code! It will be very useful for the trees I'm working with. Two
clades are large, and I need to collapse them for diversitree.
Best,
John
John S. S. Denton
Ph.D. Candidate
Department of Ichthyology and Richard Gilder Graduate School
American Museum of Natural History
www.johnssdenton.com
From: Liam J. Revell [liam.rev...@umb.edu]
Sent: Saturday, June 22, 2013 8:09 AM
To: John Denton
Cc: David Bapst; r-sig-phylo@r-project.org
Subject: Re: [R-sig-phylo] collapse descendants of a node to a polytomy
Hi John.
Here is an alternative to David's second solution that also retains the
branch lengths everywhere else in the tree the height above the root
for the tips in the star subtree. It uses phytools (and its dependencies).
collapse.to.star-function(tree,node){
tt-splitTree(tree,split=list(node=node,bp=tree$edge.length[which(tree$edge[,2]==node)]))
ss-starTree(species=tt[[2]]$tip.label,branch.lengths=diag(vcv(tt[[2]])))
ss$root.edge-0
tree-paste.tree(tt[[1]],ss)
return(tree)
}
# e.g.,
library(phytools)
set.seed(1)
tree-pbtree(n=30)
plotTree(tree,node.numbers=T)
node-fastMRCA(tree,t15,t5)
tree2-collapse.to.star(tree,node)
plotTree(tree2)
Liam J. Revell, Assistant Professor of Biology
University of Massachusetts Boston
web: http://faculty.umb.edu/liam.revell/
email: liam.rev...@umb.edu
blog: http://blog.phytools.org
On 6/22/2013 5:21 AM, David Bapst wrote:
library(ape)
set.seed(444)
tree-rtree(10)
plot(tree)
node-mrca(tree)[t4,t2]
#collapse immediate branches
tree1-tree
tree1$edge.length-rep(1,
Nedge(tree1))#replace all edge lengths with 1
tree1$edge.length[tree1$edge[,1]==node]-0#replace descendent edge
lengths with 0
tree1-di2multi(tree1)
tree1$edge.length-NULL#get rid of branch lengths
plot(tree1)
#collapse all descendant branches
tree2-tree
library(phangorn)
descNodes-Descendants(tree2,node,all)
tree2$edge.length-rep(1,Nedge(tree2))#replace all edge lengths with 1
descEdges-sapply(tree2$edge[,2],function(x) any(x==descNodes))
tree2$edge.length[descEdges]-0
tree2-di2multi(tree2)
tree2$edge.length-NULL#get rid of branch lengths
plot(tree2)
___
R-sig-phylo mailing list - R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
[R-sig-phylo] collapse descendants of a node to a polytomy
Hi folks, I'd like to collapse the descendants of a node, identified using something like node - mrca(tree)[A, B]. I did not see a function in ape, geiger, phyloch, or picante to do something like collapse.descendants(node). Is there a package with a function like this? Thanks! ~John John S. S. Denton Ph.D. Candidate Department of Ichthyology and Richard Gilder Graduate School American Museum of Natural History www.johnssdenton.com ___ R-sig-phylo mailing list - R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
[R-sig-phylo] transform nexus SETS block into treefinder partition
Hi all, PartitionFinder selected a complex partitioning scheme that I would like to put into a TreeFinder analysis. As far as I know, TreeFinder ignores NEXUS SETS blocks, and so I would have to input the partitioning by hand. I do this normally for repetitive partitions, e.g. codon positions, but the partitioning is more tedious in this case: charset p1 = 1-373\3, 2-374\3, 3-375\3, 376-958\3, 377-959\3, 378-960\3, 961-1648\3, 962-1649\3, 1652-2495\3, 2497-3241\3, 2498-3242\3, 2499-3243\3, 3244-4057\3, 4061-4979\3 charset p2 = 963-1650\3, 1653-2496\3 charset p3 = 1651-2494\3, 3245-4058\3, 3246-4059\3, 4060-4978\3, 4062-4980\3 and TreeFinder files set partitions in a line at the top of the file, like partition 111233122111 taxon1 AACGATTTCTCT taxon2 AACGATCTTAT Has anybody used an R script for transforming NEXUS SETs into the treefinder partition? Thanks, ~John John S. S. Denton Ph.D. Candidate Department of Ichthyology and Richard Gilder Graduate School American Museum of Natural History www.johnssdenton.com ___ R-sig-phylo mailing list - R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
[R-sig-phylo] treefinder chronogram problems
Hi folks, I generated a quick chronogram in TreeFinder (attached) for an exploratory analysis in R. The tree appeared fine in the Treefinder plot window, but the branch length notation does not appear to be compatible with R plotting for a chronogram. The basal divergence time is 33.9 MYA (which appears in the Newick file), but FigTree is telling me the base is something around 22. Has anyone else encountered the problem with TreeFinder chronograms? Does anybody have a suggestion for how to rescale? I'd prefer to use this tree rather than simply make the original unrooted tree ultrametric. Thanks! ~John John S. S. Denton Ph.D. Candidate Department of Ichthyology and Richard Gilder Graduate School American Museum of Natural History www.johnssdenton.com chrono3.tre Description: chrono3.tre ___ R-sig-phylo mailing list R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
[R-sig-phylo] problems with assign(), paste(), and data.frame() for folders containing trees
Hi folks,
I am trying to recurse through several numbered subfolders in a directory. Each
folder has many trees that I want to display summary values for. I have been
expanding data frames using code with the structure name - rbind(name,
newvals) to produce a data frame with n rows equal to the number of files in
one of the folders, and n column equal to the number of values in the file.
I can loop over the values within a single subdirectory fine with, for example,
library(ape)
trees - list.files(pattern=*.tre)
iters=length(trees)
branchdata.5 - data.frame()
iterations - as.character(c(1:length(trees)))
for (i in 1:iters) {
tree - read.tree(trees[i])
iteration.edges.5 - as.numeric(tree$edge.length)
branchdata.5 - rbind(branchdata.5, iteration.edges.5)
}
The problem comes when I want to iterate through the numbered subdirectories
while also iterating through the files in a given directory. I want to
recursively assign these data frames as well, with something like
f - list.dirs(path = /.../.../etc, full.names = FALSE, recursive = FALSE)
for (j in 1:length(f)) {
setwd(paste(/.../.../.,j,sep=))
assign( paste(branchdata.5,j,sep=), data.frame() )
iterations - as.character(c(1:length(trees)))
for (i in 1:iters) {
tree - read.tree(trees[i])
assign(paste(iteration.edges.5,j,sep=), as.numeric(tree$edge.length) )
paste(branchdata.5,j,sep=) - rbind(paste(branchdata.5,j,sep=),
paste(iteration.edges.5,j,sep=))
}
names(iterations) - NULL
boxplot(t(paste(branchdata.5,j,sep=)) , horizontal=TRUE , names=iterations
, ylim=c(0,2), xlab=Branch Lengths , ylab=Iterations , main = )
}
The problem seems to be in the rbind() when using values with assign() and
paste(). I would love some help on this!
John S. S. Denton
Ph.D. Candidate
Department of Ichthyology and Richard Gilder Graduate School
American Museum of Natural History
www.johnssdenton.com
___
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
[R-sig-phylo] phylogenetic tree diameter
Hi all, I'm using the ape package, and was wondering if there was a method in the package or code for calculating the diameter (http://mathworld.wolfram.com/GraphDiameter.html) of phylogenetic tree objects. I had been looking in to exporting trees from ape into a graph theory package (igraph), but it looks like it would be more difficult to write a file conversion script than to try something in ape itself. Also, I have tried installing the phybase package from CRAN but keep getting errors. Does anyone have or know of a good script for calculating Robinson-Foulds differences for phy objects? Thanks! ~J -- ORGANIC LIFE beneath the shoreless waves Was born and nurs'd in ocean's pearly caves; First forms minute, unseen by spheric glass, Move on the mud, or pierce the watery mass; These, as successive generations bloom, New powers acquire and larger limbs assume; Whence countless groups of vegetation spring, And breathing realms of fin and feet and wing. ~Erasmus Darwin, The Temple of Nature Canto I.V ___ R-sig-phylo mailing list R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
[R-sig-phylo] outputting values during taxon pruning
Hi, all.
I'm working on a script that iteratively removes taxa from an initial
total tree and re-searches for optimal (pruned) topologies after
outputting a pruned alignment.
I'm trying to get a few things to happen:
First, the script iteratively checks to make sure that the taxon names
match in the alignments and trees, and terminates if they do not. I can
get it to print that the matches are okay as it checks each taxon:
for (i in totaltr$tip.label) {
if (all(totaltr$tip.label %in% rownames(alig)) all(rownames(alig)
%in%
totaltr$tip.label)) cat(OK: Taxon labels match in tree and alignment)
else stop(Taxon labels do not match. Script terminates.)
}
But I want each line of OK: Taxon labels match and Taxon labels do
not match to print the name of the taxon being checked, or the taxon for
which a match does not occur, something like
Taxon labels match in tree and alignment FOR TAXON i. and
Taxon labels do not match FOR TAXON i. Script terminates.
Suggestions?
Second, I want the script to store a list of the tree lengths it
calculates as it goes. I can calculate tree lengths via
sum(tree.name$edge.length), but I am having trouble, it seems,
initializing and storing the list. The code I have is
ts.values - list()
length(ts.values) - length(row.names(alig))
for (i in totaltr$tip.label) {
if (all(totaltr$tip.label %in% rownames(alig)) all(rownames(alig)
%in%
totaltr$tip.label)) cat(OK: Taxon labels match in tree and alignment)
else stop(Taxon labels do not match. Script terminates.)
}
Sys.sleep(3)
for (i in totaltr$tip.label) {
write.tree(drop.tip(totaltr, i), file = paste(minus_, i, .tre, sep =
) )
ts.values[[i]] - {sum(drop.tip(totaltr, i)$edge.length)}
}
So ts.values is the list I am trying to store the tree lengths in. After
the script finishes, though, I try to call ts.values via print(ts.values),
but it says the object does not exist.
I am tinkering with these two problems using a modified version of the
phymltest routine in the ape package, so the tree/alignment filenames
written are for phyml output.
Suggestions?
Thanks a lot.
~John
--
ORGANIC LIFE beneath the shoreless waves
Was born and nurs'd in ocean's pearly caves;
First forms minute, unseen by spheric glass,
Move on the mud, or pierce the watery mass;
These, as successive generations bloom,
New powers acquire and larger limbs assume;
Whence countless groups of vegetation spring,
And breathing realms of fin and feet and wing.
~Erasmus Darwin, The Temple of Nature Canto I.V
___
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
[R-sig-phylo] iterative pruning from alignment and tree
Hi folks, I'm starting to use APE, and am trying to write a script to iteratively prune taxa from an alignment / tree combination and print the reduced datasets to separate files. I was thinking of using some combination of drop.tip + DNAbin, but I am unsure of how to link removing the tip labels in the tree and the corresponding taxon labels in the alignment in such a way that the iteration completes without redundancy. Is there a function or loop you could recommend? ~John ___ R-sig-phylo mailing list R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
