[[cursos:bot99:aulas:praticas:pratica06catego]]

Reconstrução Ancestral Discreta no R 01

O script abaixo faz uma reconstrução ancestral de uma variável discreta (categórica) sobre uma filogenia, usando o método de Máxima Verossimilhança implementado na função Ace do Pacote Ape do R.
O tutorial tem usa os seguintes arquivos:

#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 categórica
acecateg = ace(traits$Categorico,arv, type = "discrete", method = "ML",CI = TRUE,model = "ER", ip = 0.1)
#obtive um erro:
#Error in ace(traits$Categorico, arv, type = "discrete", method = "ML",  : "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,"Categorico"]
acecateg = ace(variavel,novaarv, type = "discrete", method = "ML",CI = TRUE,model = "ER", ip = 0.1)
#veja o resultado da reconstrucao
acecateg
#pegue a reconstrucao para cada nó da filogenia, na mesma ordem que aparecem no objeto novaarv
reconstrucao = acecateg$lik.anc
#isso é uma tabela onde cada linha é um nó da filogenia e cada coluna é a probabilidade do estado do caracter para o nó
head(reconstrucao,3)
#portanto, a soma dos valores das linhas deve ser igual a 1 para todos os nós:
rowSums(reconstrucao)
#vamos gerar um pdf com essa reconstrução
d1 = Sys.Date()
arq = paste(d1,'reconstrucaoDiscretaML_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 os tips (desnecessario no exemplo, porque a variavel tem valores de cores que o R entende)
ascores = c("blue","green","red","yellow") #tantas cores quantos estados de caracter
coresdostips = ascores[as.numeric(as.factor(variavel))]
#no exemplo poderia ser apenas coresdostipos=variavel
#plota um circulo (21) na posicao dos terminais
tiplabels(pch=21,bg=coresdostips,cex=1.5)
 
#plota os valores nos nós como pizzas
#define cores para os estados
cores = ascores[as.numeric(as.factor(colnames(reconstrucao)))] #no exemplo é o próprio nome dos estados de variacao
#no 
nodelabels(text=NULL,frame='none',cex=0.8,pie=reconstrucao,piecol=cores)
#adiciona uma legenda
texto = colnames(reconstrucao)
legend('topright',legend=texto,pch=21,pt.bg=cores,pt.cex=2,inset=0.05,title ='Variável categórica',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 categórico hipotético",side=1,line=3,cex=0.8)
#fecha o pdf
dev.off()

seguinte figura

  • cursos/bot99/aulas/praticas/pratica06catego.txt
  • Última modificação: 25/38/2015 23:38
  • por labotam_admin