library(ape)
library(phytools)
library(geiger)Roteiro prático — A árvore e os dados
Aula 1 · Métodos Filogenéticos Comparativos
0. Pacotes e o básico do R
Os objetos que mais vamos usar: vetores com nomes (um valor por espécie), tabelas (data.frame) e matrizes.
x <- c(sp1 = 2.51, sp2 = 2.80, sp3 = 2.23)
x["sp2"]sp2
2.8 names(x)[1] "sp1" "sp2" "sp3"d <- data.frame(sp = names(x), A = x)
d[d$A > 2.3, ] sp A
sp1 sp1 2.51
sp2 sp2 2.801. O objeto phylo
Uma árvore no R é uma lista da classe phylo. Vamos ler a árvore-guia do curso a partir do texto Newick.
guia <- read.tree(text = "(((sp1:0.18,sp2:0.18):1.06,(sp3:0.02,sp4:0.02):1.22):0.20,(sp5:1.37,(sp6:0.23,sp7:0.23):1.14):0.07);")
guia
Phylogenetic tree with 7 tips and 6 internal nodes.
Tip labels:
sp1, sp2, sp3, sp4, sp5, sp6, ...
Rooted; includes branch length(s).str(guia)List of 4
$ edge : int [1:12, 1:2] 8 9 10 10 9 11 11 8 12 12 ...
$ edge.length: num [1:12] 0.2 1.06 0.18 0.18 1.22 0.02 0.02 0.07 1.37 1.14 ...
$ Nnode : int 6
$ tip.label : chr [1:7] "sp1" "sp2" "sp3" "sp4" ...
- attr(*, "class")= chr "phylo"
- attr(*, "order")= chr "cladewise"As peças: edge (cada linha é um ramo, do nó pai ao filho), tip.label (as pontas são os números 1 a n), Nnode e edge.length, na mesma ordem de edge. Os nós internos são numerados de n + 1 em diante, e a raiz é n + 1.
cbind(guia$edge, comprimento = guia$edge.length) comprimento
[1,] 8 9 0.20
[2,] 9 10 1.06
[3,] 10 1 0.18
[4,] 10 2 0.18
[5,] 9 11 1.22
[6,] 11 3 0.02
[7,] 11 4 0.02
[8,] 8 12 0.07
[9,] 12 5 1.37
[10,] 12 13 1.14
[11,] 13 6 0.23
[12,] 13 7 0.23plot(guia, label.offset = .03)
nodelabels(frame = "circle", bg = "white", cex = .8)
edgelabels(round(guia$edge.length, 2), frame = "none", adj = c(.5, -.4), cex = .75)
2. Escrevendo uma árvore: Newick passo a passo
Os mesmos cinco passos dos slides. Rode um por vez e compare o texto com o desenho.
passos <- c("(A,B);", "(A:2.75,B:3.76);", "((A:2.75,B:3.76):2.0,C:2.6);",
"((A:2.75,B:3.76)N1:2.0,C:2.6);", "((A:2.75,B:3.76)N1:2.0,C:2.6)N2;")
t5 <- read.tree(text = passos[5])
t5$node.label[1] "N2" "N1"plot(t5, show.node.label = TRUE, cex = 1.3, x.lim = c(0, 7))
E o caminho inverso, da árvore para o texto:
write.tree(t5)[1] "((A:2.75,B:3.76)N1:2,C:2.6)N2;"3. Lendo arquivos: Nexus e Newick
barb <- read.nexus("BarbetTree.nex") # Nexus (salvo pelo Mesquite)
cent <- read.tree("Centrarchidae.tre") # Newick
barb
Phylogenetic tree with 42 tips and 41 internal nodes.
Tip labels:
M_asiatica_asiatica, M_lineata_hodgsoni, M_australis_australis, M_haemacephala_intermedia, M_eximia, M_viridis, ...
Rooted; includes branch length(s).cent
Phylogenetic tree with 28 tips and 27 internal nodes.
Tip labels:
Acantharchus_pomotis, Lepomis_gibbosus, Lepomis_microlophus, Lepomis_punctatus, Lepomis_miniatus, Lepomis_auritus, ...
Rooted; includes branch length(s).Uma árvore sem comprimentos de ramo (só a topologia) e a mesma árvore datada, dos anuros de Gomes:
gomes <- read.tree("Gomes_phylo.txt")
gomes_dat <- read.tree("Gomes_tree_datada.txt")
is.null(gomes$edge.length) # TRUE: só a topologia[1] TRUEis.ultrametric(gomes_dat)[1] FALSEO is.ultrametric() usa uma tolerância numérica: arredondamentos no arquivo podem fazer uma árvore datada “falhar”. Compare a maior e a menor distância raiz–ponta:
alturas <- node.depth.edgelength(gomes_dat)[1:Ntip(gomes_dat)]
range(alturas)[1] 154.1 154.1diff(range(alturas)) # diferença de arredondamento, não de biologia[1] 5e-06is.ultrametric(gomes_dat, tol = 1e-4)[1] TRUEplot(gomes_dat, cex = .7); axisPhylo()
4. Enraizamento, politomias e comprimentos
is.rooted(barb); is.binary(barb); is.ultrametric(barb)[1] TRUE[1] TRUE[1] TRUEDesenraizar e enraizar de novo muda a leitura da árvore (slides, bloco 1):
u <- unroot(read.tree(text = "((A:1,B:1):1,(C:1,D:1):1);"))
r1 <- root(u, "A", resolve.root = TRUE)
par(mfrow = c(1, 3))
plot(u, type = "unrooted", main = "não enraizada", cex = 1.6, lab4ut = "axial")
plot(r1, main = "raiz no ramo de A", cex = 1.6)
plot(read.tree(text = "((A,B),(C,D));"), main = "raiz no ramo central", cex = 1.6)
Uma politomia e como resolvê-la (aleatoriamente, com ramos de comprimento zero) para os métodos que exigem árvore binária:
p <- read.tree(text = "(A,(B,C,D));")
is.binary(p)[1] FALSEset.seed(1)
is.binary(multi2di(p))[1] TRUEQuando só há topologia, compute.brlen() atribui comprimentos (por exemplo, pelo método de Grafen). Os comprimentos são inventados: servem para desenhar ou explorar, não como tempo.
g1 <- compute.brlen(gomes, method = "Grafen")
is.ultrametric(g1)[1] TRUEplot(g1, cex = .7); axisPhylo()
5. A árvore e a tabela de atributos
Os atributos ficam numa planilha; para os métodos comparativos, quase sempre precisamos de um vetor com os nomes das espécies.
obj <- read.csv("Centrarchidae.csv", row.names = 1)
head(obj) feeding.mode gape.width buccal.length
Acantharchus_pomotis pisc 0.114 -0.009
Lepomis_gibbosus non -0.133 -0.009
Lepomis_microlophus non -0.151 0.012
Lepomis_punctatus non -0.103 -0.019
Lepomis_miniatus non -0.134 0.001
Lepomis_auritus non -0.222 -0.039gape <- setNames(obj$gape.width, rownames(obj))
head(gape)Acantharchus_pomotis Lepomis_gibbosus Lepomis_microlophus
0.114 -0.133 -0.151
Lepomis_punctatus Lepomis_miniatus Lepomis_auritus
-0.103 -0.134 -0.222 Agora o passo que mais dá problema: os nomes da árvore e da tabela batem?
dados_barb <- read.csv("Barbetdata.csv", row.names = 1)
nc <- name.check(barb, dados_barb)
nc$tree_not_data
[1] "M_asiatica_chersonesus" "M_australis_australis"
[3] "M_haemacephala_celestinoi" "M_haemacephala_delica"
[5] "M_haemacephala_intermedia" "M_haemacephala_rosea"
[7] "M_lineata_lineata" "M_oorti_faber"
[9] "M_oorti_sini"
$data_not_tree
character(0)Nove subespécies da árvore não têm dados. Duas saídas: podar a árvore (drop.tip()) ou deixar o treedata() casar as duas coisas de uma vez.
casado <- treedata(barb, dados_barb, warnings = FALSE)
Ntip(casado$phy); nrow(casado$data)[1] 33[1] 33all(casado$phy$tip.label == rownames(casado$data))[1] FALSEAtenção à última linha: por padrão, o treedata() não põe as linhas da tabela na ordem das pontas. Funções que usam vetores com nomes (como phylosig() e fitContinuous()) casam pelos nomes e não se importam; as que usam a ordem das linhas (como gls() sem form = ~sp) dão resultados errados em silêncio. Peça a ordenação, ou ordene você mesmo:
casado <- treedata(barb, dados_barb, sort = TRUE, warnings = FALSE)
all(casado$phy$tip.label == rownames(casado$data))[1] TRUE# ou: dados_barb[casado$phy$tip.label, ]O treedata() devolve uma matriz de texto quando a tabela tem colunas de tipos diferentes; converta os atributos numéricos antes de usar.
asa <- setNames(as.numeric(casado$data[, "wing"]), rownames(casado$data))
plotTree.barplot(casado$phy, asa, args.plotTree = list(fsize = .5),
args.barplot = list(xlab = "comprimento da asa (log)"))
6. Por que a árvore importa: dois atributos independentes
A simulação dos slides (bloco 0), passo a passo. Simulamos uma árvore e dois atributos que evoluem por movimento browniano sem nenhuma relação entre si.
set.seed(21)
arv <- rcoal(n = 100)
x <- fastBM(arv); y <- fastBM(arv)
ajuste_ols <- lm(y ~ x)
anova(ajuste_ols)Analysis of Variance Table
Response: y
Df Sum Sq Mean Sq F value Pr(>F)
x 1 3.275 3.2745 8.7871 0.003809 **
Residuals 98 36.520 0.3726
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1phylomorphospace(arv, cbind(x, y), label = "off", node.size = c(.4, .6))
abline(ajuste_ols, lwd = 2, lty = "dashed", col = "firebrick")
O filomorfoespaço mostra o motivo: espécies próximas estão próximas no gráfico porque herdaram valores parecidos. Com os contrastes independentes (Aula 2), a “relação” desaparece:
ix <- pic(x, arv); iy <- pic(y, arv)
anova(lm(iy ~ ix - 1))Analysis of Variance Table
Response: iy
Df Sum Sq Mean Sq F value Pr(>F)
ix 1 0.182 0.18238 0.1988 0.6567
Residuals 98 89.917 0.91752 Rode de novo com outras sementes (set.seed(2), set.seed(3)…): quantas vezes o OLS dá p < 0,05? E os contrastes?
7. Suporte dos nós: um bootstrap de verdade
A árvore com bootstrap dos slides (bloco 1), com as 15 sequências de camundongos-do-mato do ape.
data(woodmouse)
f <- function(s) root(nj(dist.dna(s, model = "K80")), 1, resolve.root = TRUE)
nj_wm <- f(woodmouse)
set.seed(103)
bs <- boot.phylo(nj_wm, woodmouse, f, B = 100, quiet = TRUE, rooted = TRUE)
bs [1] 100 31 31 55 48 44 33 12 43 70 84 83 46 9plot(nj_wm, cex = .8)
nodelabels(bs, frame = "none", adj = c(1.2, -.4), cex = .75)
8. A forma da árvore
O índice de Colless soma, em cada nó, a diferença entre o número de pontas dos dois lados. Escrever a função é um bom exercício de navegação no objeto phylo:
colless <- function(tr) {
n <- Ntip(tr)
dif <- sapply((n + 1):(n + tr$Nnode), \(no) {
filhos <- tr$edge[tr$edge[, 1] == no, 2]
pontas <- sapply(filhos, \(f) if (f <= n) 1 else sum(getDescendants(tr, f) <= n))
abs(diff(pontas))
})
sum(dif)
}
c(balanceada = colless(stree(8, "balanced")), pente = colless(stree(8, "left")))balanceada pente
0 21 O γ de Pybus compara a posição dos nós no tempo com a esperada sob um processo de nascimento puro: γ < 0, nós concentrados perto da raiz.
gammaStat(cent)[1] -0.28378649. Para fora do R: obter árvores
Estes pacotes não estão instalados nos computadores do curso (precisam de internet); o código fica de referência, sem rodar aqui.
# Open Tree of Life: árvore sintética para uma lista de nomes
library(rotl)
nomes <- tnrs_match_names(c("Rhinella icterica", "Physalaemus cuvieri", "Scinax hayii"))
arv_otl <- tol_induced_subtree(ott_ids = ott_id(nomes))
# phyndr: trocar espécies da árvore sem dados por equivalentes com dados
# (ver o roteiro de 2020 e Pennell et al. 2016)Resumo
- Uma árvore é um objeto
phylo:edge,tip.label,Nnode,edge.length. - Newick e Nexus leem com
read.tree()eread.nexus(); confira sempre se a árvore é enraizada, binária e ultramétrica. - Antes de qualquer análise:
name.check()etreedata(), e confira se a ordem das linhas é a das pontas (sort = TRUE). - Dados de espécies aparentadas não são independentes: o OLS pode achar relação onde não há.
O roteiro da edição de 2020 traz mais exemplos de manipulação de dados e o phyndr.