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()
Figura da Reconstrução