Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Working similarity matrix code sensu Leinster & Cobbold 2012
require(reshape)
require(picante)
require(gtools)
#Function that lists all terminal tip labels descended from ancestral node
internal2tips.self = function (phy, int.node){
#require(picante); require(ape)
Ntaxa = length(phy$tip.label)
Nnode = phy$Nnode
if ((Ntaxa + Nnode - 1) != nrow(phy$edge)) {
print("tree structure error")
break
}
if (mode(int.node) == "character"){
nodes = which(phy$node.label == int.node) + Ntaxa
}else nodes = int.node
tips = c()
repeat {
nodes = phy$edge[which(phy$edge[, 1] %in% nodes), 2]
if (length(nodes) == 0)
break
tips = c(tips, nodes)
}
tips = tips[tips <= Ntaxa]
if( int.node <= Ntaxa & length(tips) == 0 ){
tips = int.node
}
tips = phy$tip.label[tips]
return(tips)
}
#Simulated ultrametric phylogeny with 3 species
phy <- read.tree(text="(3:3,(1:2,2:2):1);")
phy$tip.label <- as.numeric(phy$tip.label)
plot(phy)
nodelabels()
edge <- phy$edge[,2]
taxa <- length(phy$tip.label)
Bnode = taxa + 1
Lb = t(t(dist.nodes(phy)[1:taxa, Bnode]))
That <- sum(Lb * p)
p = rep((1/taxa), taxa)
S <- rep(1:length(phy$tip.label))
sp <- matrix(NA, nrow = length(S) + 1, ncol = length(edge) + 1)
sp[1,-1] <- sort(edge)
sp[-1,1] <- S
edges <- sp[1,]
for (i in (2:ncol(sp))) {
branch <- internal2tips.self(phy,edges[i])
for (j in (2:nrow(sp))){
sp[j,i] <- ifelse(sp[j,1] %in% branch, 1, 0)
}
}
#Number of historic species
Nhist <- sum(sp[-1,-1])
histsp <- sp[-1,-1]
colnames(histsp) <-sort(edge)
rownames(histsp) <- S
hs <- melt(histsp, varnames = c("species", "branch"))
hs$abund <- rep(p, length(edge))
hs$Lb <- rep(Lb, length(edge))
hs <- (subset(hs, hs$value != 0))
hs <- hs[order(hs$species, -hs$branch), ]
#Initialize Z similarity matrix
Z <- matrix(NA, Nhist, Nhist)
for(i in 1:Nhist){
branch <- internal2tips.self(phy,hs$branch[i])
for(j in 1:Nhist){
Z[i,j] <- ifelse(hs$species[j] %in% branch, That/Lb[hs$species[j]], 0)
}
}
#Print Similarity matrix
Z
pi = matrix((hs$Lb / That)*hs$abund)
Z %*% pi
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.