(* ::Package:: *) (* This function prints the version number for the installed verison of this package. *) PhylogeneticsVersion[]:=Print["PollyPhylogenetics 2.0\n(c) P. David Polly, 23 July 2012\n"]; PhylogeneticsVersion[] (* This functions reads one or more Newick format trees from a text file. Usage: ReadNewick[FileName] where FilePath is the name (possibly including path) to the tree file. Created by P. David Polly, 21 March 2012 *) ReadNewick[TreeFile_]:=Module[{Tree,NumTrees,NumNodes,NumBranchLengths,NumTips}, Tree=Import[TreeFile,"Text"]; Tree=StringReplace[Tree,{RegularExpression["[^\\x21-\\x7E]"]->""}]; If[StringMatchQ[StringTake[Tree,-1],";"],Tree=Tree,Tree=Tree<>";"]; NumTrees=Length[StringCases[Tree,";"]]; NumNodes=Length[StringCases[Tree,"("]]/NumTrees; NumBranchLengths=Length[StringCases[Tree,":"]]/NumTrees; NumTips=(Length[StringCases[Tree,","]]/NumTrees)+1; If[NumNodes>NumBranchLengths,Return["Your tree does not have lengths for all branches"]]; If[NumTrees>1,(Tree=#<>";"&/@StringSplit[Tree,";"])]; Return[Tree]; ]; (* This function transforms a Newick format tree into a table where each row represents a branch in the tree, the first column giving the name of the branch end point, the second column giving the name of the branch ancestor, the third column giving the length of the branch, and the fourth column indicating whether the branch ends in a tip (1) or a node (0). The table is returned with a header line that must be discarded for functions that use the table as input. Usage: TreeToTable[tree] where tree is a Newick format tree in a single string. The tree must have branch lengths and must end in a semicolon. Example Tree: "((A:1,B:1):1,C:1);" Created by P. David Polly, 21 March 2012 *) TreeToTable[Tree_]:=Module[{LastNode,SplitTree,Pos,TreeTable,Tip,Node,ParseTree}, (* This subfunction processes a tip, reading the name and length and inserting it into TreeTable *) Tip[AncNum_]:=Block[{MyName, MyNum, MyAnc, MyLength}, MyAnc="Node "<>ToString[AncNum]; MyName=""; MyLength=""; While[StringMatchQ[SplitTree[[Pos]], RegularExpression["[^:]"]], MyName=MyName<>SplitTree[[Pos]]; Pos=Pos+1; ]; Pos=Pos+1; While[StringMatchQ[SplitTree[[Pos]], RegularExpression["[^,);]"]], MyLength=MyLength<>SplitTree[[Pos]]; Pos=Pos+1; ]; MyLength=ToExpression[MyLength]; TreeTable=Append[TreeTable,{MyName,MyAnc,MyLength, 1}]; Pos=Pos-1; Return[]; ]; (* This subfunction processes a node, calling new node functions or tips as they are encountered within the node *) Node[AncNum_]:=Block[{MyName, MyAnc,MyLength,MyNum}, MyLength=""; MyNum=LastNode; LastNode=LastNode+1; MyName="Node "<>ToString[MyNum]; MyAnc="Node "<>ToString[AncNum]; While[StringMatchQ[SplitTree[[Pos]],RegularExpression["[^);]+"]], Pos=Pos+1; If[StringMatchQ[SplitTree[[Pos]],"("],Node[MyNum], If[StringMatchQ[SplitTree[[Pos]], RegularExpression["[A-Za-z0-9]"]],Tip[MyNum], If[StringMatchQ[SplitTree[[Pos]],RegularExpression["[);]"]],Break[]]]]; ]; Pos=Pos+1; If[StringMatchQ[SplitTree[[Pos]],":"],Pos=Pos+1]; While[StringMatchQ[SplitTree[[Pos]], RegularExpression["[^,);]+"]], MyLength=MyLength<>SplitTree[[Pos]]; Pos=Pos+1; ]; MyLength=ToExpression[MyLength]; TreeTable=Append[TreeTable,{MyName,MyAnc,MyLength,0}]; Return[]; ]; (* main code *) TreeTable={{"Descendant","Ancestor","Branch Length", "Tip?"}}; Pos=1; LastNode=0; SplitTree=StringSplit[Tree,""]; While[True, If[StringMatchQ[SplitTree[[Pos]],"("],Node[LastNode]]; If[StringMatchQ[SplitTree[[Pos]],";"],Break[]]; ]; Return[TreeTable[[1;;-2]]]; ]; (* This function calculates the three matrices used by Martins and Hansen (1997) for phylogenetic regressions. The matrices are: (1) Var[Y], which is the matrix of shared ancestry of tree tips; (2) Var[AY], which is the matrix of shared ancestry of tips; and (3) Var[A], which is the matrix of shared ancestry of tree nodes; and and nodes. This function takes a tree table as its argument, in the form produced by the function TreeToTable[] (but without the row of column labels!). The three matrices are returned with column and row labels. Usage: {varY,varAY,varA}=PhylogeneticMatrices[tree]; where varY, varAY, and varA are the three output matrices described above and tree is a string containing a tree with branch lengths in Newick format or a tree table of the sort produced by TreeToTable[]. Created by P. David Polly, updated 20 May 2012 *) PhylogeneticMatrices[tree_]:=Module[{FindLineage,TipEntries,NodeEntries,TipLineages,NodeLineages,varY,varA,varAY,i,j,myTable}, If[StringQ[tree], myTable=TreeToTable[tree][[2;;]], If[tree[[1]]=={"Descendant","Ancestor","Branch Length","Tip?"}, myTable=tree[[2;;]], If[Length[tree[[1]]]==4,myTable=tree,Print["Tree variable must be a Newick tree string or a tree table."];];]; ]; FindLineage[Point_]:=Module[{Lineage,NewBranch}, Lineage={Point}; While[!StringMatchQ[Lineage[[-1,2]],"Node 0"], NewBranch=Flatten[myTable[[Flatten[Position[myTable[[1;;,1]],Lineage[[-1,2]]]]]]]; Lineage=Append[Lineage,NewBranch]; ]; Return[Lineage]; ]; TipEntries=myTable[[Flatten[Position[StringCases[myTable[[1;;,1]],"Node "~~__],{}]]]]; NodeEntries=myTable[[Flatten[Position[myTable[[1;;,1]],#]&/@Flatten[StringCases[myTable[[1;;,1]],"Node "~~__]]]]]; TipLineages=FindLineage[#]&/@TipEntries; NodeLineages=FindLineage[#]&/@NodeEntries; varY=Table[Table[Plus@@(Intersection[TipLineages[[i]],TipLineages[[j]]][[1;;,3]]),{i,Length[TipLineages]}],{j,Length[TipLineages]}]; varA=Table[Table[Plus@@(Intersection[NodeLineages[[i]],NodeLineages[[j]]][[1;;,3]]),{i,Length[NodeLineages]}],{j,Length[NodeLineages]}]; varA=Transpose[Prepend[Transpose[Prepend[varA,Table[1,{Length[varA[[1]]]}]]],Table[1,{Length[varA[[1]]]+1}]]]; varAY=Table[Table[Plus@@(Intersection[TipLineages[[i]],NodeLineages[[j]]][[1;;,3]]),{i,Length[TipLineages]}],{j,Length[NodeLineages]}]; varAY=Prepend[varAY,Table[0,{Length[varAY[[1]]]}]]; Return[{Prepend[Transpose[Prepend[varY,TipEntries[[1;;,1]]]],Prepend[TipEntries[[1;;,1]],"Var Y"]],Transpose[Prepend[Transpose[Prepend[varAY,TipEntries[[1;;,1]]]],Prepend[Prepend[NodeEntries[[1;;,1]],"Node 0"],"Var AY"]]],Prepend[Transpose[Prepend[varA,Prepend[NodeEntries[[1;;,1]],"Node 0"]]],Prepend[Prepend[NodeEntries[[1;;,1]],"Node 0"],"Var A"]]}]; ]; (* This function reconstructs ancestral quantitative values for traits on a phylogenetic tree following Martins & Hansen, 1997, Am. Nat., 149:646-667. Estimate of the standard error at the base of the tree is from Rohlf, 2001, Evolution, 55: 2143-2160. Usage: {rate, nodes, se} = NodeReconstruct[tree, trait] where tree is a string containing a tree with branch lengths in Newick format or a tree table of the sort produced by TreeToTable[]. The function returns the estimated per-step rate of the trait on the tree, the ancestral values of Y at each node in the tree assuming a Brownian motion, and the uncertainty in the ancestral values due to the Brownian motion evolutionary process expressed as 1 standard deviation of the variance due to the evolutionary process (the "standard error" of some authors). Created by P. David Polly, updated 20 May 2012 *) ReconstructNodes[tree_,trait_]:=Module[{J,GrandMean,Yprime,Nodes,Rate,x,contrasts, i,j,SDNode,X,SyGLS,varY,varAY,varA,TipLabels,NodeLabels,Y,VarY,VarAY,VarA,myTable,RateSquared,VarYNew}, If[StringQ[tree],myTable=TreeToTable[tree][[2;;]], If[tree[[1]]=={"Descendant","Ancestor","Branch Length","Tip?"},myTable=tree[[2;;]], If[Length[tree[[1]]]==4,myTable=tree,Print["Tree variable must be a Newick tree string or a tree table."];]; ]; ]; {VarY,VarAY,VarA}=PhylogeneticMatrices[myTable]; varY=VarY[[2;;,2;;]];varAY=VarAY[[2;;,2;;]]; varA=VarA[[2;;,2;;]]; TipLabels=VarY[[2;;,1]];NodeLabels=VarA[[2;;,1]]; If[Length[trait[[1]]]!=2,Return["Trait must be a two-column matrix with tip labels in the first column and trait values in the second column"]]; If[Sort[TipLabels]!=Sort[trait[[1;;,1]]],Return["Traits labels do not match tip labels from Var Y."]]; Y=trait[[1;;,2]][[Flatten[Position[trait[[1;;,1]],#]&/@TipLabels]]]; J=Table[1,{Length[varY]}]; GrandMean=(J.Inverse[varY].Y)/(J.Inverse[varY].J); Yprime=#-GrandMean&/@Y; Nodes=#+GrandMean&/@(varAY.Inverse[varY].Yprime); contrasts=IndependentContrasts[tree,trait]; RateSquared=Mean[contrasts^2]; Rate=Sqrt[RateSquared]; SDNode={}; Do[ VarYNew=PhylogeneticMatrices[RerootTree[myTable,NodeLabels[[i]]]][[1,2;;,2;;]]; SyGLS=Sqrt[RateSquared/(Plus@@Plus@@Inverse[VarYNew])]//N; SDNode=Append[SDNode,{NodeLabels[[i]],SyGLS}] ,{i,Length[NodeLabels]}]; Return[{{"Rate",Rate},Transpose[Prepend[{Nodes},#&/@NodeLabels]],SDNode}]; ]; (* This function draws a Newick format tree and labels the nodes and tips, making use of helpful suggestions given by Felsenstein (2004). Usage: DrawNewickTree[tree] where tree tree is a string containing a tree with branch lengths in Newick format or a tree table of the sort produced by TreeToTable[], The function requires that taxon names be unique. It handles polytomies, but fails when the tree format is not specified correctly or when branch lengths are not given. Created by P. David Polly, 5 April 2012 *) DrawNewickTree[tree_]:=Module[{TreeTable,myTable,ExtNodeNames,ImmediateDescendants,FindLineage,TipEntries,NodeEntries,NodeNames,TipNames,TipY,TipLineages,NodeLineages,Descendants,NodeY,TreeX,TreeY,TreeXY,x}, FindLineage[Point_]:=Module[{Lineage,NewBranch}, Lineage={Point}; While[!StringMatchQ[Lineage[[-1,2]],"Node 0"], NewBranch=Flatten[TreeTable[[Flatten[Position[TreeTable[[1;;,1]],Lineage[[-1,2]]]]]]]; Lineage=Append[Lineage,NewBranch]; ]; Return[Lineage]; ]; If[StringQ[tree], myTable=TreeToTable[tree][[2;;]], If[tree[[1]]=={"Descendant","Ancestor","Branch Length","Tip?"}, myTable=tree[[2;;]], If[Length[tree[[1]]]==4,myTable=tree,Print["Tree variable must be a Newick tree string or a tree table."];];]; ]; TreeTable=myTable; Clear[myTable]; TipEntries=TreeTable[[Flatten[Position[StringCases[TreeTable[[1;;,1]],"Node "~~__],{}]]]]; NodeEntries=TreeTable[[Flatten[Position[TreeTable[[1;;,1]],#]&/@Flatten[StringCases[TreeTable[[1;;,1]],"Node "~~__]]]]]; NodeNames=NodeEntries[[1;;,1]]; TipNames=TipEntries[[1;;,1]]; TipY=Transpose[{TipNames,Table[x,{x,Length[TipNames],1,-1}]}]; TipLineages=FindLineage[#]&/@TipEntries; NodeLineages=FindLineage[#]&/@NodeEntries; Descendants=Table[Cases[If[MemberQ[#,{__,NodeNames[[x]],__}],#[[1,1]]]&/@TipLineages,Except[Null]],{x,Length[NodeNames]}]; NodeY=Transpose[{NodeNames,Table[Mean[TipY[[Flatten[Position[TipY[[1;;,1]],#]&/@Descendants[[x]]],2]]]//N,{x,Length[Descendants]}]}]; TreeX=Union[Sort[Partition[Flatten[Table[Transpose[{TipLineages[[x,1;;,1]],Reverse[Drop[FoldList[Plus,0,Reverse[TipLineages[[x,1;;,3]]]],1]]}],{x,Length[TipLineages]}]],2]]]; TreeY=Sort[Flatten[{TipY,NodeY},1]]; TreeXY=Transpose[{TreeX[[1;;,1]],TreeX[[1;;,2]],TreeY[[1;;,2]]}]; TreeXY=Sort[Append[TreeXY,{"Node 0",0,Mean[TreeY[[1;;,2]]]}]]; ExtNodeNames=Reverse[Sort[Append[NodeNames,"Node 0"]]]; ImmediateDescendants=TreeTable[[Flatten[Position[TreeTable[[1;;,2]],#]],1]]&/@ExtNodeNames; Do[Do[ TreeXY[[Flatten[Position[TreeXY[[1;;,1]],ExtNodeNames[[x]]]],3]]=Mean[TreeXY[[Flatten[Position[TreeXY[[1;;,1]],#]&/@ImmediateDescendants[[x]]],3]]]//N, {x,Length[ExtNodeNames]}],{2}]; Return[ Graphics[ { (* Horizontal lines *) {Thick,Line[{Flatten[TreeXY[[Flatten[Position[TreeXY[[1;;,1]],#[[1]]]],{2,3}]]],Flatten[{Flatten[TreeXY[[Flatten[Position[TreeXY[[1;;,1]],#[[2]]]],{2}]]],Flatten[TreeXY[[Flatten[Position[TreeXY[[1;;,1]],#[[1]]]],{3}]]]}]}]&/@TreeTable}, (* Vertical lines *) {Thick,Line[{Flatten[{Flatten[TreeXY[[Flatten[Position[TreeXY[[1;;,1]],#[[2]]]],{2}]]],Flatten[TreeXY[[Flatten[Position[TreeXY[[1;;,1]],#[[1]]]],{3}]]]}],Flatten[TreeXY[[Flatten[Position[TreeXY[[1;;,1]],#[[2]]]],{2,3}]]]}]&/@TreeTable}, (* tip labels *) Table[{FontFamily->"Arial",FontSize->12,Text[TipNames[[x]],Flatten[TreeXY[[Flatten[Position[TreeXY[[1;;,1]],TipNames[[x]]]],{2,3}]]],{-1.3,0}]},{x,Length[TipNames]}], (* node labels *) Table[{FontFamily->"Arial",FontSize->10,Text[NodeNames[[x]],Flatten[TreeXY[[Flatten[Position[TreeXY[[1;;,1]],NodeNames[[x]]]],{2,3}]]],{-1.3,0}]},{x,Length[NodeNames]}], (* Base node and label *) {PointSize[0.02], Point[Flatten[TreeXY[[Flatten[Position[TreeXY[[1;;,1]],"Node 0"]],{2,3}]]]]}, Text["Node 0",Flatten[TreeXY[[Flatten[Position[TreeXY[[1;;,1]],"Node 0"]],{2,3}]]],{-1.3,0}] } ,Axes->{True,False},AxesOrigin->{0,0*Length[TipNames]},BaseStyle->{FontFamily->"Arial",FontSize->10},AspectRatio->1/GoldenRatio] ]; ]; (* This function takes a list of string labels and converts them to PHYLIP format by truncating labels greater than 10 characters or padding ones less than 10 characters. Usage: phyliplabels = MakePHYLIPLabels[labels] where labels is a list of strings. Created by P. David Polly, 5 April 2012 *) MakePHYLIPLabels[labels_]:=Module[{phyliplabels,x}, phyliplabels=labels; Return[Table[If[StringLength[phyliplabels[[x]]]<= 10, phyliplabels[[1]]=phyliplabels[[x]]<>StringJoin[Table[" ",{10-StringLength[phyliplabels[[x]]]}]],StringDrop[phyliplabels[[x]], -1*(StringLength[phyliplabels[[x]]]-10)]],{x,Length[phyliplabels]}]]; ] (* This function performs a brownian motion random walk for n generations with a step rate of i. Usage: walk = RandomWalk[100, 1]; Created by P. David Polly, 13 April 2012 *) RandomWalk[n_,i_]:=NestList[#+Random[NormalDistribution[0,i]]&,0,n]; (* This function returns standardized Independent Contrasts (Felsenstein, 1985) for a continuous trait. Usage: contrasts = IndependentContrasts[tree, trait]; where tree is a phylogentic tree with branch lengths in Newick format and trait is a a two-column matrix with taxon labels in the first column (which must match the taxon names in tree) and trait values in the second column. Created by P. David Polly, Updated 19 May 2012 *) IndependentContrasts[tree_,trait_]:=Module[{myTable, contrasts,myData0,x,stdcontrasts,FindDescendants,ProcessTip,TraverseNode}, FindDescendants[me1_]:=Module[{}, Return[Select[myTable,#[[2]]==me1&]] ]; ProcessTip[me2_]:=Module[{myLength,myID,myValue}, myLength=Select[myTable,#[[1]]==me2&][[1,3]]; myValue=Select[trait,#[[1]]==me2&][[1,2]]; Return[{myValue,myLength}]; ]; TraverseNode[me3_]:=Module[{myDescendants,myLength,myValue,myData}, myLength=Select[myTable,#[[1]]==me3&][[1,3]]; myDescendants=FindDescendants[me3]; myData=If[#[[4]]==0,TraverseNode[#[[1]]],ProcessTip[#[[1]]]]&/@myDescendants; Do[ contrasts=Append[contrasts,{myData[[x,1]]-myData[[x+1,1]],myData[[x,2]]+myData[[x+1,2]]}]; ,{x,Length[myData]-1}]; myValue=(Plus@@((1/#[[2]])*#[[1]]&/@myData))/(Plus@@(1/#[[2]]&/@myData)); myLength=myLength+(Fold[Times, 1, myData[[1;;,2]]]/(Plus@@myData[[1;;,2]])); (* Return[{myValue,myLength}]; *) Return[{myValue,myLength}]; ]; If[StringQ[tree], myTable=TreeToTable[tree][[2;;]], If[tree[[1]]=={"Descendant","Ancestor","Branch Length","Tip?"}, myTable=tree[[2;;]], If[Length[tree[[1]]]==4,myTable=tree,Return["Tree variable must be a Newick tree string or a tree table."];];]; ]; If[Length[Dimensions[trait]]!= 2, Return["Trait variable must have two columns"]]; If[Sort[trait[[1;;,1]]]!= Sort[Union[Select[myTable,#[[4]]==1&][[1;;,1]]]], Return["Trait labels must equal tree names"]]; contrasts={}; myData0=If[#[[4]]==0,TraverseNode[#[[1]]],ProcessTip[#[[1]]]]&/@FindDescendants["Node 0"]; Do[ contrasts=Append[contrasts,{myData0[[x,1]]-myData0[[x+1,1]],myData0[[x,2]]+myData0[[x+1,2]]}//N]; ,{x,Length[myData0]-1}]; stdcontrasts=#[[1]]/Sqrt[#[[2]]]&/@contrasts; Return[stdcontrasts]; ]; (* This function converts a tree in tabular form to a Newick format string. . Usage: tree = TableToTree[treetable]; where treetable is a table in the same format as produced by the TreeToTable[] function. Created by P. David Polly, 3 May 2012 *) TableToTree[myTable_]:=Module[{contrasts,myData0,x,stdcontrasts,newick,myValue0,myLength0,myEntry0,FindDescendants,ProcessTip,TraverseNode}, FindDescendants[me1_]:=Module[{}, Return[Select[myTable,#[[2]]==me1&]] ]; ProcessTip[me2_]:=Module[{myLength,myEntry}, myLength=Select[myTable,#[[1]]==me2&][[1,3]]; myEntry=me2<>":"<>ToString[myLength]; Return[myEntry]; ]; TraverseNode[me3_]:=Module[{myDescendants,myLength,myData,myEntry}, myLength=Select[myTable,#[[1]]==me3&][[1,3]]; myDescendants=FindDescendants[me3]; myData=If[#[[4]]==0,TraverseNode[#[[1]]],ProcessTip[#[[1]]]]&/@myDescendants; myEntry="("<>StringJoin[Riffle[myData,","]]<>"):"<>ToString[myLength]; Return[myEntry]; ]; myData0=If[#[[4]]==0,TraverseNode[#[[1]]],ProcessTip[#[[1]]]]&/@FindDescendants["Node 0"]; myEntry0="("<>StringJoin[Riffle[myData0,","]]<>");"; Return[myEntry0]; ]; (* This function reroots a Newick tree at the specified node. Usage: newtree = RerootTree[tree, node]; where tree is a string containing a tree with branch lengths in Newick format or a tree table of the sort produced by TreeToTable[], and node is the name of the node where the new tree will be rooted. The name of the node can be found by plotting the tree with DrawNewickTree[tree]. Created by P. David Polly, 10 May 2012 *) RerootTree[tree_,root_]:=Module[{myData0,x,myLength0,myEntry0,TablePosition,myLine,myTable,FindLineage, FindDescendants,ProcessTip,TraverseNode}, FindLineage[Point_]:=Module[{Lineage,NewBranch}, Lineage={Point}; While[!StringMatchQ[Lineage[[-1,2]],"Node 0"], TablePosition=Append[TablePosition,Flatten[Position[myTable[[1;;,1]],Lineage[[-1,2]]]]]; NewBranch=Flatten[myTable[[TablePosition[[-1]]]]]; Lineage=Append[Lineage,NewBranch]; ]; Return[{Lineage,Flatten[TablePosition]}]; ]; FindDescendants[me1_]:=Module[{}, Return[Select[myTable,#[[2]]==me1&]] ]; ProcessTip[me2_]:=Module[{myLength,myEntry}, myLength=Select[myTable,#[[1]]==me2&][[1,3]]; myEntry=me2<>":"<>ToString[myLength]; Return[myEntry]; ]; TraverseNode[me3_]:=Module[{myDescendants,myLength,myData,myEntry}, myLength=Select[myTable,#[[1]]==me3&][[1,3]]; myDescendants=FindDescendants[me3]; myData=If[#[[4]]==0,TraverseNode[#[[1]]],ProcessTip[#[[1]]]]&/@myDescendants; myEntry="("<>StringJoin[Riffle[myData,","]]<>"):"<>ToString[myLength]; Return[myEntry]; ]; If[StringQ[tree], myTable=TreeToTable[tree][[2;;]], If[tree[[1]]=={"Descendant","Ancestor","Branch Length","Tip?"}, myTable=tree[[2;;]], If[Length[tree[[1]]]==4,myTable=tree,Print["Tree variable must be a Newick tree string or a tree table."];];]; ]; TablePosition={}; myLine=FindLineage[myTable[[ Flatten[Position[myTable[[1;;,2]],root]][[1]] ]] ]; Do[ myTable[[TablePosition[[x]],{1,2}]]=myTable[[TablePosition[[x]],{2,1}]], {x,Length[TablePosition]}]; myData0=If[#[[4]]==0,TraverseNode[#[[1]]],ProcessTip[#[[1]]]]&/@FindDescendants[root]; myEntry0="("<>StringJoin[Riffle[myData0,","]]<>");"; Return[myEntry0]; ]; (* This function simulates the evolution of a continuous trait on a tree using a Brownian motion model of evolution. By default the trait starts with a value of 0.0 at the root of the tree (this can be changed by specifying a different value as the third input parameter of the function). Evolution is simulated along each branch by Brownian motion using the rate specified by the second input parameter. The function returns a list of node and tip labels followed by the simulated trait value. Usage: SimulateContinuousTraitOnTree[tree, rate (,startvalue, includenodes)] where tree is a string containing a tree with branch lengths in Newick format or a tree table of the sort produced by TreeToTable[], rate is the amount of trait change per unit of branch length (as specified in the tree), startvalue is the optional starting value of the trait at the root of the tree, and includenodes is the optional parameter that specifies whether simulated node values should be returned (1) or not returned (0, default). Created by P. David Polly, 20 April 2012 *) SimulateContinuousTraitOnTree[tree_,Rate_,StartValue_:0,IncludeNodes_:0]:=Module[{myTable,x,myData0,TraitValues,FindDescendants,ProcessTip,TraverseNode}, FindDescendants[me_]:=Module[{}, Return[Select[myTable,#[[2]]==me&]] ]; ProcessTip[me1_,myAnc1_,myLength1_,myStart1_]:=Module[{myValue}, myValue=Random[NormalDistribution[myStart1,Sqrt[myLength1]*Rate]]; TraitValues=Append[TraitValues,{me1,myValue}]; Return[]; ]; TraverseNode[me2_,myAnc2_,myLength2_,myStart2_]:=Module[{myDescendants,myValue}, myValue=Random[NormalDistribution[myStart2,Sqrt[myLength2]*Rate]]; TraitValues=Append[TraitValues,{me2,myValue}]; myDescendants=FindDescendants[me2]; If[#[[4]]==0,TraverseNode[#[[1]],#[[2]],#[[3]],myValue],ProcessTip[#[[1]],#[[2]],#[[3]],myValue]]&/@myDescendants; Return[]; ]; If[StringQ[tree], myTable=TreeToTable[tree][[2;;]], If[tree[[1]]=={"Descendant","Ancestor","Branch Length","Tip?"}, myTable=tree[[2;;]], If[Length[tree[[1]]]==4,myTable=tree,Print["Tree variable must be a Newick tree string or a tree table."];];]; ]; TraitValues={{"Node 0",StartValue}}; If[#[[4]]==0,TraverseNode[#[[1]],#[[2]],#[[3]],StartValue],ProcessTip[#[[1]],#[[2]],#[[3]],StartValue]]&/@FindDescendants["Node 0"]; If[IncludeNodes==0,TraitValues=TraitValues[[Flatten[Position[TraitValues[[1;;,1]],#]&/@Complement[TraitValues[[1;;,1]],Union[Flatten[StringCases[TraitValues[[1;;,1]],"Node "~~DigitCharacter..]]]]]]]]; Return[TraitValues]; ]; (* This function simulates the evolution of correlated continuous traits on a tree using a Brownian motion model of evolution. By default the traits each start with a value of 0.0 at the root of the tree (this can be changed by specifying a vector of values as the third input parameter of the function). Evolution is simulated along each branch by Brownian motion using the rate specified by the second input parameter with covariance between the traits as specified in the matrix. The number of traits is determined from the structure of the covariance matrix. The rates of evolution are derived from the trace of the covariance matrix, with the rate at each step of the simulation being equivalent to a random step whose variance is equal to the variance of that trait in the covariance matrix. The function returns a list of node and tip labels followed by the simulated trait values. Usage: SimulateCorrelatedTraitsOnTree[tree, covariancematrix (,startvector, includenodes)] where tree is a string containing a tree with branch lengths in Newick format or a tree table of the sort produced by TreeToTable[], covariancematrix is the covariance matrix of the traits. per unit of branch length (as specified in the tree), startvector is the optional starting values of the traits at the root of the tree, and includenodes is the optional parameter that specifies whether simulated node values should be returned (1) or not returned (0, default). Created by P. David Polly, 23 July 2012 *) SimulateCorrelatedTraitsOnTree[tree_,P_,StartVector_:0,IncludeNodes_:0]:=Module[{myTable,x,myData0,TraitValues,FindDescendants,ProcessTip,TraverseNode,Chol}, FindDescendants[me_]:=Module[{}, Return[Select[myTable,#[[2]]==me&]] ]; ProcessTip[me1_,myAnc1_,myLength1_,myStart1_]:=Module[{myValue}, myValue=Table[Random[NormalDistribution[0,1]],{Length[Chol]}].(Sqrt[myLength1]*Chol)+myStart1; TraitValues=Append[TraitValues,{me1,myValue}]; Return[]; ]; TraverseNode[me2_,myAnc2_,myLength2_,myStart2_]:=Module[{myDescendants,myValue}, myValue=Table[Random[NormalDistribution[0,1]],{Length[Chol]}].(Sqrt[myLength2]*Chol)+myStart2; TraitValues=Append[TraitValues,{me2,myValue}]; myDescendants=FindDescendants[me2]; If[#[[4]]==0,TraverseNode[#[[1]],#[[2]],#[[3]],myValue],ProcessTip[#[[1]],#[[2]],#[[3]],myValue]]&/@myDescendants; Return[]; ]; If[StringQ[tree], myTable=TreeToTable[tree][[2;;]], If[tree[[1]]=={"Descendant","Ancestor","Branch Length","Tip?"}, myTable=tree[[2;;]], If[Length[tree[[1]]]==4,myTable=tree,Print["Tree variable must be a Newick tree string or a tree table."];];]; ]; Chol=CholeskyDecomposition[P]; TraitValues={{"Node 0",If[StartVector==0,Table[0,{Length[Chol]}],StartVector]}}; If[#[[4]]==0,TraverseNode[#[[1]],#[[2]],#[[3]],StartVector],ProcessTip[#[[1]],#[[2]],#[[3]],StartVector]]&/@FindDescendants["Node 0"]; If[IncludeNodes==0,TraitValues=TraitValues[[Flatten[Position[TraitValues[[1;;,1]],#]&/@Complement[TraitValues[[1;;,1]],Union[Flatten[StringCases[TraitValues[[1;;,1]],"Node "~~DigitCharacter..]]]]]]]]; Return[TraitValues]; ];