Distograma filogenético

Walkthrough do código — Aula 3: Estimativa de estado ancestral e sinal filogenético

Walkthrough do código usado em sala. Os comandos não são executados pelo renderizador — copie e cole no seu R local. Você pode baixar o script .R original e os datasets em Aula 3.

# Generates either a distogram or an autocorrelogram based on Hardy, O. J., and S. Pavoine. 2012. Evolution 66:2614-2621 and Diniz-Filho, J. A. F. 2001. Evolution 55:1104-1109
# code is modified from Hardy and Pavoine 2012
# includes a comparison to Brownian motion model (simulations scaled by the original data) and to a case of no phylogenetic structure (permutations of the original data)
# Input is data matrix (for multivariate data) or vector (for univariate data) and tree as a phylo object
# dcfun specifies whether the Euclidean distance ("ED") is calculated or Moran's I ("IM")
# nsims specifies the number of simulations for the BM simulations and permutations
# default (plot=TRUE) generates the plot itself (either distogram or autocorrelogram), otherwise you can set plot=FALSE and plot it yourself based on the output
# output is a matrix in which every row is a class of patristic distance over which the distance/similarity is averaged. columns are the mean patristic distance for that class ("pdcm"), the observed value of distance/similarity ("obs"), the mean ("Pm"), lower ("Pl"), and upper ("Pu") bounds of the 95% CI of permutations, and the mean ("BMm"), lower ("BMl"), and upper ("BMu") bounds of the 95% CI of BM simulations

require(abind); require(geiger)

distcor <- function(X, tree, dcfun=c("ED","IM"), nsims = 100, plot=TRUE) {
    X <- as.matrix(X)
    pd <- cophenetic.phylo(tree)
    pd <- pd/max(pd)
    class <-c(0, (1:10)/10)
    if (dcfun=="ED") class <- class^2
    pdcm <- c()
    for(j in 2:length(class)) {pdcm <- c(pdcm, mean(pd[pd>class[j-1] & pd<=class[j]]))}
    XP <- array(apply(array(X, dim=c(nrow(X),ncol(X),nsims)), 3, function(X){X[sample(1:nrow(X),nrow(X)),]}), dim=c(nrow(X),ncol(X),nsims)) #permutations of species to simulate no phylogenetic structure
    XS <- sim.char(tree, var(X), nsim=nsims) # simulations of BM
    XX <- abind(XP, XS, X, along=3)
    if (dcfun=="ED") dfun <- function(X){as.matrix(dist(X))} else if (dcfun=="IM") dfun <- function(X){X%*%t(X)/sum(diag(var(X))) + 1/(nrow(X)-1)}
    cM0 <- c()
    for (h in 1:dim(XX)[3]) {
        X <- as.matrix(XX[,,h])
        D <- dfun(X)
        cm <- c()
        for(j in 2:length(class)) {cm <- c(cm, mean(D[pd>class[j-1] & pd<=class[j]]))}
        cM0 <- cbind(cM0, cm)           
        }
    na <- which(is.na(pdcm))
    cMP <- cM0[-na,1:nsims]
    cMS <- cM0[-na,(nsims+1):(nsims*2)]
    cM <- cbind("pdcm"=pdcm[-na], "obs"=cM0[-na,ncol(cM0)], "Pm"=rowMeans(cMP), "Pl"=apply(cMP, 1, quantile, 0.025), "Pu"=apply(cMP, 1, quantile, 0.975), "BMm"=rowMeans(cMS), "BMl"=apply(cMS, 1, quantile, 0.025), "BMu"=apply(cMS, 1, quantile, 0.975))
    if (plot==TRUE) {
        ltys <- c(NA,1,3,3,3,1,1,1)
        cols <- c(NA,"Black","black","grey","grey","black","grey","grey")
        if(dcfun=="ED") {
            cM[,"pdcm"] <- sqrt(cM[,"pdcm"])
            legpos <- "bottomright"
            yttl <- "Euclidean distance"
            xttl <- "Square root of patristic distance"
            } 
        else if(dcfun=="IM") {
            legpos <- "topright"
            yttl <- "Moran's I (trait similarity)"
            xttl <- "patristic distance"
            }
        plot(cM[,"pdcm"], cM[,"obs"], ylim=c(min(0,min(cM[,-1])),max(cM[,-1])), pch=19, xlab=xttl, ylab=yttl)
        for (j in 2:ncol(cM)) {lines(cM[,"pdcm"], cM[,j], lty=ltys[j], col=cols[j])}
        legend(legpos, legend=c("Observed","Permutated","Brownian"), cex=0.8, lty=c(1,3,1), pch=c(19,NA,NA))
        }
    cM
    }