Reconstrução Ancestral Contínua no R

O script abaixo faz uma reconstrução ancestral de uma variável contínua (numérica) sobre uma filogenia, usando a função Ace do Pacote Ape do R.
O tutorial usa os seguintes arquivos:

• rubiaceaephymlgtr_rps16.tree - arvore filogenética em formato newick [usar read.nexus() ao invés de read.tree() se tua árvore estiver em formato nexus].
• rubiaceaetraitsinventados.csv - tabela com duas variáveis, uma categórica e outra quantitativa usadas nos exemplos. Essas variáveis são inventadas e não representam nada.
• Este script do R - gera a Figura Abaixo.

#chama o pacote filogenético
library(ape)
#remove qualquer objeto
rm(list=ls())
#lista arquivos na pasta de trabalho
dir()
#le a arvore que aparece nesta lista
arv = read.tree("rubiaceaePhyMLgtr_rps16.tree")
#le uma matriz de atributos (o arquivo esta separado por \t tabulacao, a primeira linha tem nome das colunas e as colunas estao entre aspas==quote)
traits = read.table("rubiaceaeTraitsInventados.csv",sep='\t',as.is=T,header=T,quote="\"")
#adiciona os tips como nome de linhas aos traits
rownames(traits) = traits[,"TipLabel"]
#faz a reconstrução ancestral da variável quantitativa
acequant = ace(traits$Quantitativo,arv, type = "continuous", method = "REML",CI = TRUE,model = "BM")
#obtive um erro:
#Error in ace(traits$Quantitativo, arv, type = "continuous", method = "REML",  : "phy" is not rooted AND fully dichotomous.
#isso quer dizer que há politomias na filogenia e a reconstrução não é possível
#no caso do exemplo, basta remover 1 taxon do grupo externo
novaarv = drop.tip(arv,tip="Fockea_apensis_KP245405.1")
#refaz a reconstrução ancestral da variável categórica
rn = novaarv$tip.label
variavel = traits[rn,"Quantitativo"]
acequant = ace(variavel,novaarv, type = "continuous", method = "REML",CI = TRUE,model = "BM")
#Há muitos parâmetros e você deve entendê-los. Veja também as funcões anc.ML() e anc.Bayes() do pacote Phytools que são outras alternativas 
#veja o resultado da reconstrucao
acequant
#pegue a reconstrucao para cada nó da filogenia, na mesma ordem que aparecem no objeto novaarv
reconstrucao = acequant$ace
#isso é um vetor onde cada valor é a reconstrucao  da variavel quantitativa para cada ancestral (nó)
reconstrucao
 
#vamos gerar um pdf com essa reconstrução
d1 = Sys.Date()
arq = paste(d1,'reconstrucaoContinua_rubiaceae_rps16_GTR.pdf',sep='')
#abre o pdf
pdf(file=arq,height=10,width=8)
#muda a margem
par(mar=c(5,2,2,1))
#plota a filogenia
plot.phylo(novaarv,type='phylogram',use.edge.length=T,show.tip.label=T,show.node.label=F,cex=0.8,label.offset =0.01)
#note o argumento label.offset, distanciando um pouco a legenda dos tips, para poder plotar um símbolo
#plota o estado do caracter para os terminais
#define cores para a variavel numérica
#vamos usar um gradiente de cores
#juntas os valores atuais e ancestrais para isso
osval = sort(unique(c(reconstrucao,variavel)))
#corestodas = rainbow(n=length(osval), end=0.8)
corestodas = colorRampPalette(c("green","yellow","red"))(length(osval))
names(corestodas) = osval
coresdostips = corestodas[as.character(variavel)]
#no exemplo poderia ser apenas coresdostipos=variavel
#plota um circulo (21) na posicao dos terminais
tiplabels(text=round(variavel,2),frame='circle',bg=coresdostips)
#plota os valores nos nós como pizzas
#define cores para os estados
cores = corestodas[as.character(reconstrucao)] #no exemplo é o próprio nome dos estados de variacao
#no 
nodelabels(text=round(reconstrucao,2),frame='circle',cex=0.5,bg=cores)
#adiciona uma legenda
ss = seq(1,length(corestodas),by=6)
ss = unique(c(ss,length(corestodas)))
legcor = corestodas[ss]
legtxt = round(as.numeric(names(corestodas)[ss]),2)
legend('topright',legend=legtxt,pch=21,pt.bg=legcor,pt.cex=2,inset=0.05,title ='Variável quantitativa',box.lwd = 0.5)
#adiciona uma legenda para a figura
mtext(text="Árvore de Máxima Verossimilhança (PhyML) do marcador \nde cloroplasto rps16 intron1 para algumas Rubiaceae, com a reconstrução ancestral de um caractere quantitativo hipotético",side=1,line=3,cex=0.8)
#fecha o pdf
dev.off()