################ LIBRARIES:
library (ape)
library (phytools)
## Loading required package: maps
library(vegan)
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.2-1
############### INPUT PARAMETERS:
N.trees = 1000 # No. trees to process
N.burnin = 3000 # No. of trees to discard (burn-in set in Bayesian analyses)
#
############### PACo FUNCTION:
PACo <- function (H.dist, P.dist, HP.bin)
{
HP.bin <- which(HP.bin > 0, arr.in=TRUE)
H.PCo <- pcoa(H.dist, correction="cailliez")$vectors #Performs PCo of Host distances
P.PCo <- pcoa(P.dist, correction="cailliez")$vectors #Performs PCo of Parasite distances
H.PCo <- H.PCo[HP.bin[,1],] #adjust Host PCo vectors
P.PCo <- P.PCo[HP.bin[,2],] ##adjust Parasite PCo vectors
list (H.PCo = H.PCo, P.PCo = P.PCo)
}
####### Read trees and compress for memory efficiency::
treeH <- .compressTipLabel(read.nexus(file= "C:/Users/balbuena/Dropbox/MCRALink/1K_30A/30A_10/30A_10a_Bayesian/infile.nex.run1.t")) # change input files as needed
treeP <- .compressTipLabel(read.nexus(file= "C:/Users/balbuena/Dropbox/MCRALink/1K_30A/30A_10/30A_10a_Bayesian/infile.nex.run2.t")) # id.
HP <- as.matrix(read.table("C:/Users/balbuena/Dropbox/HP.txt", header=TRUE))
# consensus trees:
treeCH <- read.nexus(file= "C:/Users/balbuena/Dropbox/MCRALink/1K_30A/30A_10/30A_10a_Bayesian/infile.nex.con.tre") # change input files as needed
treeCP <- read.nexus(file= "C:/Users/balbuena/Dropbox/MCRALink/1K_30A/30A_10/30A_10a_Bayesian/infile.nex.con.tre") # id.
#
# apply PACo to consensus trees:
host.D <- cophenetic (treeCH)
para.D <- cophenetic (treeCP)
PACo.fit <- PACo(host.D, para.D, HP)
HP.proc <- procrustes(PACo.fit$H.PCo, PACo.fit$P.PCo) #Procrustes Ordination
#
# Drop burn-in set:
if (length(treeH) >= N.burnin+2 & length(treeP) >= N.burnin+2) {
treeH <- treeH[(N.burnin+1): length(treeH)]
treeP <- treeP[(N.burnin+1): length(treeP)]
} else stop("Insufficient number of trees to perform the analysis. Check your input files")
# Select a random saple of trees of size N.trees
if (length(treeH) >= N.trees & length(treeP) >= N.trees) {
treeH <- sample(treeH, N.trees)
treeP <- sample(treeP, N.trees)
} else stop("Insufficient number of trees to perform the analysis. Check your input files")
#
# D.wrapper reads 1 nexus tree, executes PACo and returns squared residuals
D.wrapper <- function(n) {
DH.add <- cophenetic(treeH[[n]]) #
DP.add <- cophenetic(treeP[[n]]) #
PACo.add <- PACo(DH.add, DP.add, HP)
Proc.add <- procrustes(PACo.add$H.PCo, PACo.add$P.PCo) #Procrustes Ordination
add.res <- residuals(Proc.add)^2
}
# Matrix of squared residuals:
SQres <- suppressWarnings(t(sapply(1:N.trees, D.wrapper)))
colnames (SQres) <- paste(rownames(HP.proc$X),rownames(HP.proc$Yrot), sep="-") #colnames identify the H-P links
medianSQ <- apply(SQres, 2, median) # median values
meanSQ <- apply(SQres, 2, mean) # mean values
upperCI <- apply(SQres, 2, quantile, probs=0.975) # upper 95% confidence interval
lowerCI <- apply(SQres, 2, quantile, probs=0.025) # lower CI
# Plot
pat.bar <- barplot(medianSQ, names.arg = " ", space = 0.25, col="white", xlab= "Host-parasite links",
ylab= "Squared residuals", ylim=c(0, max(upperCI)), cex.lab=1.2)
text(pat.bar, par("usr")[3] - 0.001, srt = 330, adj = 0, labels = colnames(SQres), xpd = TRUE, font = 1, cex=0.6)
arrows(pat.bar, lowerCI, pat.bar, upperCI, length= 0.05, angle=90, code=3)
abline(a=HP.proc$ss/sum(HP), b=0, lty=2, xpd=FALSE) #draws a line across m**2/number of H-P links
# Alt Plot
plot(medianSQ, ylim=c(min(lowerCI),max(upperCI)), xaxt="n", xlab= "Host-parasite links",
ylab= "Squared residuals" )
axis(1, 1:sum(HP), labels= FALSE,mgp=c(1,0.5,0))
text(1:sum(HP), par("usr")[3] - 0.07, labels = colnames(SQres), srt = 90, adj = c(1,0.5),
cex=0.6, xpd = TRUE)
arrows(1:sum(HP), upperCI, 1:sum(HP),lowerCI, code=3, length=0.02,angle=90,col='blue')
abline(a=HP.proc$ss/sum(HP), b=0, lty=2, xpd=FALSE, col ="red")
#
#optional
points(1:sum(HP), meanSQ, pch=20)
points(1:sum(HP), residuals(HP.proc)^2, pch=4)
