Desenhando árvores com mr bayes

Mr Bayes é um software para inferência bayesiana e escolha de modelos para uma grande variedade de modelos filogenéticos e evolutivos. Ele usa Markov Chain Monte Carlo (MCMC) para estimar a distribuição posterior dos parâmetros do modelo.
Basicamente ele tem um interface de linha de comando, que é similar a sintaxe usada por outro programa muito comum, o Paup, mas o Paup é um software proprietário e o Mr é um Software Livre.

Certo, bem ele pode ser conseguido aqui, e depois de baixar, no linux é só rodar um ./configure, sendo que temos que ter a biblioteca beagle instalada, senão temos que avisar o ./configure para ele não pedir o beagle, feito isso é só compilar com o make. Pro windows ele vem como um instalador, o que torna as coisas mais fácil, mas uma vez com ele instalado, podemos começar a dar uma olhada na sintaxe básica, ja que para como compilar, tem um capítulo inteiro no manual do Mr Bayes so de como compilar do fonte.

Bem, como entrada, usamos um alinhamento no formato nexus, permite guardar mais informações além das bases em si.

A cara dele é assim:

Captura de tela de 2016-05-31 15-38-30

Beleza, vamos usar o arquivo primatemitDNA.nexus, eu vou deixar ele disponível no meu repositório do github, mas basicamente é um alinhamento de 5 sequências de dna mitocondrial de primatas. É o arquivo do exercício do curso de Evolução Molecular que fiz no coursera.

Para abrir o arquivo no Mr Bayes, a gente usa o comando execute.

MrBayes > execute primatemitDNA.nexus Executing file "primatemitDNA.nexus" UNIX line termination Longest line length = 911 Parsing file Expecting NEXUS formatted file Reading data block Allocated taxon set Allocated matrix Defining new matrix with 5 taxa and 896 characters Data is Dna Gaps coded as - Taxon 1 -> Human Taxon 2 -> Chimpanzee Taxon 3 -> Gorilla Taxon 4 -> Orangutan Taxon 5 -> Gibbon Successfully read matrix Setting default partition (does not divide up characters) Setting model defaults Seed (for generating default start values) = 1464723572 Setting output file names to "primatemitDNA.nexus.run." Exiting data block Reached end of file MrBayes >

Veja que a gente abre as sequencias, de 911 bases, da para ver o nome das sequencias e algumas outras informações. A gente pode conferir ainda se os dados estão ok usando o comando ShowMatrix.

MrBayes > Showmatrix 1 Human AAGCTTCACCGGCGCAGTCATTCTCATAATCGCCCACGGACTTACATCCT Chimpanzee AAGCTTCACCGGCGCAATTATCCTCATAATCGCCCACGGACTTACATCCT Gorilla AAGCTTCACCGGCGCAGTTGTTCTTATAATTGCCCACGGACTTACATCAT Orangutan AAGCTTCACCGGCGCAACCACCCTCATGATTGCCCATGGACTCACATCCT Gibbon AAGCTTTACAGGTGCAACCGTCCTCATAATCGCCCACGGACTAACCTCTT . . . 851 Human ACTGACACTGAGCCACAACCCAAACAACCCAGCTCTCCCTAAGCTT Chimpanzee ACTGGCACTGAGCAACAACCCAAACAACCCAGCTCTCCCTAAGCTT Gorilla GCTGACACTGAGCAACAACCCAAACAATTCAACTCTCCCTAAGCTT Orangutan ACTGATGCTGAACAACCACCCAGACACTACAACTCTCACTAAGCTT Gibbon ACTGACACTGAACTGCAACCCAAACGCTAGAACTCTCCCTAAGCTT

Eu suprimi parte da saída, mas assim a gente olha os dados que foram carregados.
Agora vamos definir um “outgroup” para a árvore, para termos uma árvore com raiz

MrBayes > outgroup Gibbon Setting outgroup to taxon "Gibbon" e vamos definor o modelo que vamos usar MrBayes > lset nst=2 rates=gamma Setting Nst to 2 Setting Rates to Gamma Successfully set likelihood model parameters

Aqui estamos especificando que modelo evolutivo vamos usar. Nst é o número de parâmetros do modelo evolutivo, dizer que queremos 2 parâmetros aqui, dois tipos de substituições, usando o modelo de kimura, eu acho pelo menos, que distinguimos entre transições e transversões. Além disso, rates=gamma significa que nos queremos usar a distribuição gamma para representar diferentes taxas em diferentes partes da sequência.

Podemos mudar outras coisas, mas por enquanto vamos tentar fazer o básico, que é desenhar a árvore, podemos ver como é o modelo que está definido até agora com Showmodel, mas antes de iniciar a análise essa informação é mostrada, então não precisamos nos preocupar.

Para iniciar o MCMC, usamos o comando mcmc, com os parâmetros

ngen=100000 : tamanho da cadeia
samplefreq=100 : de quanto em quanto vamos salvar uma árvore e todos os parâmetros (amostras do distribuição posterior)

nchains=3 : número de cadeias paralelas, mas isso significa que uma é “fria”, onde as amostras são feitas e duas são “quentes” (“heated”) que se movem no espaço do parâmetros mais rapidamente para char picos na distribuição de probabilidade.

diagnfreq=5000 tem haver com ver se o MrBayes está rodando corretamente. Resumidamente, MrBayes vai começar duas sessões independentes, começando de diferentes árvores aleatórios. No começo, as duas sessões vão ter amostras bem diferentes, de árvores diferentes, mas quando elas convergirem (a distribuição posterior começar a ficar igual), as duas árvores devem ser bem parecidas, toda diagonfreq gerações, o MrBayes mede o quão similar as amostras das árvores estão e como uma regra geral (rule of thumb), nos vamos esperar esse valor ficar menor que 0.01.

Mas vamos la, vamos rodar o modelo.

MrBayes > mcmc ngen=100000 samplefreq=100 nchains=3 diagnfreq=5000 Setting number of generations to 100000 Setting sample frequency to 100 Setting number of chains to 3 Setting diagnosing frequency to 5000 Running Markov chain MCMC stamp = 4955053902 Seed = 557265027 Swapseed = 1464723572 Model settings: Data not partitioned -- Datatype = DNA Nucmodel = 4by4 Nst = 2 Transition and transversion rates, expressed as proportions of the rate sum, have a Beta(1.00,1.00) prior Covarion = No # States = 4 State frequencies have a Dirichlet prior (1.00,1.00,1.00,1.00) Rates = Gamma The distribution is approximated using 4 categories. Likelihood summarized over all rate categories in each generation. Shape parameter is exponentially distributed with parameter (1.00). Active parameters: Parameters --------------------- Tratio 1 Statefreq 2 Shape 3 Ratemultiplier 4 Topology 5 Brlens 6 --------------------- 1 -- Parameter = Tratio Type = Transition and transversion rates Prior = Beta(1.00,1.00) 2 -- Parameter = Pi Type = Stationary state frequencies Prior = Dirichlet 3 -- Parameter = Alpha Type = Shape of scaled gamma distribution of site rates Prior = Exponential(1.00) 4 -- Parameter = Ratemultiplier Type = Partition-specific rate multiplier Prior = Fixed(1.0) 5 -- Parameter = Tau Type = Topology Prior = All topologies equally probable a priori Subparam. = V 6 -- Parameter = V Type = Branch lengths Prior = Unconstrained:GammaDir(1.0,0.1000,1.0,1.0) The MCMC sampler will use the following moves: With prob. Chain will use move 1.89 % Dirichlet(Tratio) 0.94 % Dirichlet(Pi) 0.94 % Slider(Pi) 1.89 % Multiplier(Alpha) 9.43 % ExtSPR(Tau,V) 9.43 % ExtTBR(Tau,V) 9.43 % NNI(Tau,V) 9.43 % ParsSPR(Tau,V) 37.74 % Multiplier(V) 13.21 % Nodeslider(V) 5.66 % TLMultiplier(V) Division 1 has 87 unique site patterns Initializing conditional likelihoods Using standard SSE likelihood calculator for division 1 (single-precision) Initial log likelihoods and log prior probs for run 1: Chain 1 -- -3213.775230 -- 14.143053 Chain 2 -- -3181.879527 -- 14.143053 Chain 3 -- -3086.072235 -- 14.143053 Initial log likelihoods and log prior probs for run 2: Chain 1 -- -3177.836224 -- 14.143053 Chain 2 -- -3222.902892 -- 14.143053 Chain 3 -- -3230.891466 -- 14.143053 There are results from a previous run saved using the same filename(s). Do you want to overwrite these results? (yes/no): yes Overwriting file "primatemitDNA.nexus.mcmc" Overwriting file "primatemitDNA.nexus.run1.p" Overwriting file "primatemitDNA.nexus.run1.t" Overwriting file "primatemitDNA.nexus.run2.p" Overwriting file "primatemitDNA.nexus.run2.t" Using a relative burnin of 25.0 % for diagnostics Chain results (100000 generations requested): 0 -- [-3213.775] (-3181.880) (-3086.072) * [-3177.836] (-3222.903) (-3230.891) 1000 -- (-2648.611) (-2710.304) [-2622.754] * (-2681.408) (-2664.530) [-2620.152] -- 0:00:00 2000 -- [-2617.298] (-2620.320) (-2615.762) * [-2617.340] (-2627.230) (-2617.470) -- 0:00:00 3000 -- (-2618.457) (-2620.569) [-2620.056] * (-2618.272) (-2618.647) [-2616.704] -- 0:00:00 4000 -- [-2623.450] (-2617.093) (-2620.776) * (-2624.657) (-2617.662) [-2618.030] -- 0:00:24 5000 -- (-2615.517) [-2616.710] (-2626.709) * [-2617.718] (-2615.441) (-2618.128) -- 0:00:19 . . . Average standard deviation of split frequencies: 0.005447 96000 -- [-2616.134] (-2616.792) (-2615.722) * (-2621.091) (-2624.519) [-2616.766] -- 0:00:00 97000 -- (-2617.755) [-2614.612] (-2617.002) * (-2621.728) [-2614.617] (-2616.776) -- 0:00:00 98000 -- [-2613.363] (-2624.533) (-2620.102) * (-2616.026) (-2619.449) [-2615.646] -- 0:00:00 99000 -- (-2617.593) [-2614.617] (-2617.432) * (-2624.591) [-2617.346] (-2620.014) -- 0:00:00 100000 -- [-2620.524] (-2620.783) (-2616.498) * (-2618.557) [-2616.113] (-2623.220) -- 0:00:00 Average standard deviation of split frequencies: 0.006120 Continue with analysis? (yes/no): no Analysis completed in 13 seconds Analysis used 13.26 seconds of CPU time Likelihood of best state for "cold" chain of run 1 was -2611.84 Likelihood of best state for "cold" chain of run 2 was -2611.87 Acceptance rates for the moves in the "cold" chain of run 1: With prob. (last 100) chain accepted proposals by move 22.7 % ( 18 %) Dirichlet(Tratio) 22.4 % ( 34 %) Dirichlet(Pi) 25.9 % ( 27 %) Slider(Pi) 35.1 % ( 34 %) Multiplier(Alpha) 1.2 % ( 0 %) ExtSPR(Tau,V) 1.8 % ( 3 %) ExtTBR(Tau,V) 1.8 % ( 3 %) NNI(Tau,V) 1.6 % ( 3 %) ParsSPR(Tau,V) 25.2 % ( 20 %) Multiplier(V) 26.4 % ( 22 %) Nodeslider(V) 25.7 % ( 26 %) TLMultiplier(V) Acceptance rates for the moves in the "cold" chain of run 2: With prob. (last 100) chain accepted proposals by move 22.7 % ( 23 %) Dirichlet(Tratio) 20.9 % ( 21 %) Dirichlet(Pi) 27.1 % ( 26 %) Slider(Pi) 34.2 % ( 28 %) Multiplier(Alpha) 1.4 % ( 0 %) ExtSPR(Tau,V) 1.3 % ( 4 %) ExtTBR(Tau,V) 1.7 % ( 0 %) NNI(Tau,V) 1.1 % ( 1 %) ParsSPR(Tau,V) 25.7 % ( 17 %) Multiplier(V) 25.8 % ( 25 %) Nodeslider(V) 25.1 % ( 28 %) TLMultiplier(V) Chain swap information for run 1: 1 2 3 ----------------------- 1 | 0.86 0.73 2 | 33001 0.87 3 | 33411 33588 Chain swap information for run 2: 1 2 3 ----------------------- 1 | 0.86 0.73 2 | 33298 0.87 3 | 33395 33307 Upper diagonal: Proportion of successful state exchanges between chains Lower diagonal: Number of attempted state exchanges between chains Chain information: ID -- Heat ----------- 1 -- 1.00 (cold chain) 2 -- 0.91 3 -- 0.83 Heat = 1 / (1 + T * (ID - 1)) (where T = 0.10 is the temperature and ID is the chain number) MrBayes >

Primeiro temos um report de todo o modelo que estamos usando, inicializamos as cadeias, e particularmente aqui, como é a segunda vez que estava rodando, ja tem um resultado pronto, mas eu peço para rodar denovo, e temos os passos do mcmc, e depois do número de passos definidos, ele pergunta se queremos continuar, como estamos com uma pequena diferença entre as duas sessões falamos que nao precisa, e pronto, ja temos a cadeia do mcmc. Arquivos foram gerados na pasta que estamos trabalhando, mas ainda não temos que sumarizar os resultados. E podemos fazer isso com o comando sumt

MrBayes > sumt contype=halfcompat conformat=simple relburnin=yes burninfrac=0.25 showtreeprobs=yes Setting sumt contype to Halfcompat Setting sumt conformat to Simple Using relative burnin (a fraction of samples discarded). Setting burnin fraction to 0.25 Setting showtreeprobs to yes Summarizing trees in files "primatemitDNA.nexus.run1.t" and "primatemitDNA.nexus.run2.t" Using relative burnin ('relburnin=yes'), discarding the first 25 % of sampled trees Writing statistics to files primatemitDNA.nexus. Examining first file ... Found one tree block in file "primatemitDNA.nexus.run1.t" with 1001 trees in last block Expecting the same number of trees in the last tree block of all files Tree reading status: 0 10 20 30 40 50 60 70 80 90 100 v-------v-------v-------v-------v-------v-------v-------v-------v-------v-------v ********************************************************************************* Read a total of 2002 trees in 2 files (sampling 1502 of them) (Each file contained 1001 trees of which 751 were sampled) Overwriting file "primatemitDNA.nexus.parts" Overwriting file "primatemitDNA.nexus.tstat" Overwriting file "primatemitDNA.nexus.vstat" Overwriting file "primatemitDNA.nexus.con.tre" Overwriting file "primatemitDNA.nexus.trprobs" General explanation: In an unrooted tree, a taxon bipartition (split) is specified by removing a branch, thereby dividing the species into those to the left and those to the right of the branch. Here, taxa to one side of the removed branch are denoted '.' and those to the other side are denoted '*'. Specifically, the '.' symbol is used for the taxa on the same side as the outgroup. In a rooted or clock tree, the tree is rooted using the model and not by reference to an outgroup. Each bipartition therefore corresponds to a clade, that is, a group that includes all the descendants of a particular branch in the tree. Taxa that are included in each clade are denoted using '*', and taxa that are not included are denoted using the '.' symbol. The output first includes a key to all the bipartitions with frequency larger or equual to (Minpartfreq) in at least one run. Minpartfreq is a parameter to sumt command and currently it is set to 0.10. This is followed by a table with statistics for the informative bipartitions (those including at least two taxa), sorted from highest to lowest probability. For each bipartition, the table gives the number of times the partition or split was observed in all runs (#obs) and the posterior probability of the bipartition (Probab.), which is the same as the split frequency. If several runs are summarized, this is followed by the minimum split frequency (Min(s)), the maximum frequency (Max(s)), and the standard deviation of frequencies (Stddev(s)) across runs. The latter value should approach 0 for all bipartitions as MCMC runs converge. This is followed by a table summarizing branch lengths, node heights (if a clock model was used) and relaxed clock parameters (if a relaxed clock model was used). The mean, variance, and 95 % credible interval are given for each of these parameters. If several runs are summarized, the potential scale reduction factor (PSRF) is also given; it should approach 1 as runs converge. Node heights will take calibration points into account, if such points were used in the analysis. Note that Stddev may be unreliable if the partition is not present in all runs (the last column indicates the number of runs that sampled the partition if more than one run is summarized). The PSRF is not calculated at all if the partition is not present in all runs.The PSRF is also sensitive to small sample sizes and it should only be considered a rough guide to convergence since some of the assumptions allowing one to interpret it as a true potential scale reduction factor are violated in MrBayes. List of taxa in bipartitions: 1 -- Human 2 -- Chimpanzee 3 -- Gorilla 4 -- Orangutan 5 -- Gibbon Key to taxon bipartitions (saved to file "primatemitDNA.nexus.parts"): ID -- Partition ----------- 1 -- *.... 2 -- .*... 3 -- ..*.. 4 -- ...*. 5 -- ****. 6 -- ***.. 7 -- **... ----------- Summary statistics for informative taxon bipartitions (saved to file "primatemitDNA.nexus.tstat"): ID #obs Probab. Sd(s)+ Min(s) Max(s) Nruns ---------------------------------------------------------------- 6 1501 0.999334 0.000942 0.998668 1.000000 2 7 1464 0.974700 0.011299 0.966711 0.982690 2 ---------------------------------------------------------------- + Convergence diagnostic (standard deviation of split frequencies) should approach 0.0 as runs converge. Summary statistics for branch and node parameters (saved to file "primatemitDNA.nexus.vstat"): 95% HPD Interval -------------------- Parameter Mean Variance Lower Upper Median PSRF+ Nruns ------------------------------------------------------------------------------------- length[1] 0.059671 0.000347 0.026276 0.095918 0.058513 1.001 2 length[2] 0.080990 0.000482 0.041008 0.126153 0.078997 1.000 2 length[3] 0.075218 0.000727 0.027528 0.127487 0.071678 1.000 2 length[4] 0.344477 0.015647 0.154275 0.596846 0.322574 1.002 2 length[5] 0.524957 0.031039 0.246341 0.882645 0.490579 1.001 2 length[6] 0.175448 0.005446 0.054393 0.326580 0.163428 0.999 2 length[7] 0.043979 0.000544 0.005098 0.088742 0.040356 1.001 2 ------------------------------------------------------------------------------------- + Convergence diagnostic (PSRF = Potential Scale Reduction Factor; Gelman and Rubin, 1992) should approach 1.0 as runs converge. NA is reported when deviation of parameter values within all runs is 0 or when a parameter value (a branch length, for instance) is not sampled in all runs. Summary statistics for partitions with frequency >= 0.10 in at least one run: Average standard deviation of split frequencies = 0.006120 Maximum standard deviation of split frequencies = 0.011299 Average PSRF for parameter values (excluding NA and >10.0) = 1.001 Maximum PSRF for parameter values = 1.002 Clade credibility values: /---------------------------------------------------------------- Gibbon (5) | |---------------------------------------------------------------- Orangutan (4) | + /--------------------- Human (1) | /----------97---------+ | | \--------------------- Chimpanzee (2) \---------100--------+ \------------------------------------------- Gorilla (3) Phylogram (based on average branch lengths): /-------------------------------------------------------------------- Gibbon (5) | |--------------------------------------------- Orangutan (4) | + /-------- Human (1) | /----+ | | \----------- Chimpanzee (2) \----------------------+ \---------- Gorilla (3) |------------| 0.100 expected changes per site Calculating tree probabilities... Tree 1 (p = 0.974, P = 0.974): /---------------------------------------------------------- | | /------------------ | /-------------------+ +-------------------+ \------------------ | | | \-------------------------------------- | \---------------------------------------------------------- Tree 2 (p = 0.015, P = 0.989): /------------------ /-------------------+ /-------------------+ \------------------ | | | \-------------------------------------- | +---------------------------------------------------------- | \---------------------------------------------------------- Tree 3 (p = 0.010, P = 0.999): /---------------------------- /-----------------------------+ | \---------------------------- | | /---------------------------- +-----------------------------+ | \---------------------------- | \---------------------------------------------------------- Tree 4 (p = 0.001, P = 1.000): /---------------------------------------------------------- | | /-------------------------------------- +-------------------+ | | /------------------ | \-------------------+ | \------------------ | \---------------------------------------------------------- Credible sets of trees (4 trees sampled): 99 % credible set contains 3 trees MrBayes >

Esse comando da um sumario das árvores estimadas. A opção contype=halfcompat faz com que a a arvore consensus seja calculada a partir das árvores que não foram descartadas no burnin. A arvore consensus é primeiro plotada na tela, abaixo está o cladograma, um phylograma consensual. As distâncias das arestas é a media sobre as árvores que tiveram essas arestas presentes. O cladograma também tem um valor de credibilidade.

O que nos interessa é a lista de árvores que está depois do phylogram. Essas são nomeadas como Tree 1, Tree 2, etc, e estão ordenadas apartir da sua probabilidade posterior, que é indicada pelo p minisculo depois de cada árvore. O P maisculo é a probabilidade acumulada, diferente de desenhar árvores com ML, aqui temos uma lista de árvores e suas probabilidades associadas.

Bem é isso ai, arquivo do alinhamento vai estar la no repositório recologia na pasta dados, e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail, e hoje não tem script, são só meia duzia de comandos.

Referência:

Anders Gorm Pedersen – Computational Molecular Evolution

Least Squares para árvores filogenéticas.

Uma coisa que sempre achei complicado de entender, era como árvores filogenéticas eram feitas com verosimilhança.
Pensando em um teste t ou uma regressão linear, sabemos que temos um conjunto de amostras, e ajustamos um modelo a eles, ajustando os parâmetros do modelo de forma a minimizar o erro em todas as amostras.

Talvez a parte complexa é imaginar o que diabos é equivalente a amostra, para usar na construção de uma árvore. Então pensando no método dos mínimos quadrados, temos o seguinte.

Primeiro, imagine que temos quatro sequências qualquer, de quatro linhagens de alguma espécie. Mas vamos criá-las no formato fasta, que é um formato padrão para guardar sequências.

>A AAAAA >B CAAAA >C CGGGG >D GGGGG

Então temos as sequências A,B,C e D, e aqui, mesmo sem análise, da para chutar fácil quem é mais parecido com quem.

Mas tudo bem, o modelo para gerar uma árvore vai ser assim.
Veja que podemos calcular a uma distância entre sequências:

> ex.dna <- read.dna("exemplo.fasta", format = "fasta") > dist.dna(ex.dna,model="raw") A B C B 0.2 C 1.0 0.8 D 1.0 1.0 0.2

O pacote ape tem a função read.dna, que nos permite ler o DNA de várias formas, então como temos um arquivos fasta, lemos o dna no formato fasta, e calculamos a distância entre as sequências, usando um modelo “raw”, que e só o número de bases diferentes dividido pelo tamanho da sequência, esse nem é um modelo usado mesmo na biologia molecular, mas depois vemos isso.
A para B é 0.2 diferente, uma base em cinco, A e C é 1 diferente, porque todas as bases são diferentes, e assim vai.

Agora pense em uma topologia qualquer, uma organização numa árvore, com tamanhos nas arestas, com uma quantidade para a diferença. Vamos criar uma.

arvore1 <- "((A:1,B:1):1,(C:1,D:1):1);" arvore1.tre <- read.tree(text=arvore1) ###Como é árvore plot(arvore1.tre)

Veja que esse formato é relativamente simples, basicamente, uma virgula divide dois indivíduos, até um ancestral.

A:1,B:1

O dois pontos, é o tamanho desse cara até o ancestral comum, ai como temos que fazer mais uniões, colocamos ele em parenteses

(A:1,B:1)

Agora esse é como se fosse um indivíduo, é um ancestral de a com b, mas é um indivíduos, veja que unimos o C e o D

(C:1,D:1)

e depois unimos os ancestrais desses dois caras

((A:1,B:1):1,(C:1,D:1):1)

O ponto e vírgula é para indicar que a arvore terminou.
E essa árvore é a seguinte:

01

Beleza, então temos uma árvore e uma topologia, agora é o seguinte.

Essa topologia é um grafo, onde as folhas são as sequências e a arestas internas são ancestrais, só para esclarecer, na teoria dos grafos, uma árvore é um grafo não direcionado e em que quaisquer dois vértices são conectados por exatamente um caminho. Em outras palavras, qualquer gráfico acíclico conectado é uma árvore. Não existem ciclos.

Captura de tela de 2016-03-20 00-28-56

Agora essas sequências tem uma distância entre elas. Veja que cada aresta tem um valor, e a distância que aparece na matriz de distâncias da sequência tem que ser a soma dessas distâncias. Assim…

Captura de tela de 2016-03-20 00-29-15

Nesse exemplo, a distância filogenética entre a sequência 1 e 2 é o percurso de a->b->c, por isso a distância é a+b+c. Mas a não ser que essa seja uma filogenia perfeita, essa conta nunca fecha. Sempre falta ou sobra um pouquinho.

Captura de tela de 2016-03-20 00-29-39

E que pouquinho que sobra é exatamente o que permite calcular o mínimo quadrado, ou seja, quais valores a gente atribui para a,b,c,d,e que fazem com que tenhamos o mínimo possível de sobra? E veja que essa conta vai ser feita para todos as distâncias, todos os pares de distâncias, logo é minimizar o valor para todos os casos, como se fossem “amostras”.

E assim o mínimo quadrado é aplicado a filogenia.

No pacote phangorn, a gente tem esse método implementado já, e por enquanto vamos usar a função de la, mas pegando aquelas nossas sequências, podemos usar a função pml, que calcula o “Likelihood of a tree”. Então com a árvore e nossas sequências, assumindo um modelo de evolução calculamos um likelihhod, logo, assim como na parcimônia podemos decidir qual a melhor árvore, mas temos um valor de verosimilhança, que diz o quão melhor ela é.

Então pensando numa segunda árvore.

02

Veja que essa árvore é muito ruim, ela junta as sequências A e D, que são completamente diferentes, e usando a função pml vemos o seguinte:

03

Agora se a gente quer fazer uma topologia é so calcular todas as topologia possíveis, com todas as combinações de tamanho de arestas possíveis, e encontramos a melhor. Veja que o valor da melhor topologia é maior (menos negativo). Lembrando que estamos trabalhando com log aqui, mas podemos voltar para os número de probabilidade com exponencial

> format(exp(ajuste1$logLik),scientific = F)
[1] “0.000000000003288564”
> format(exp(ajuste2$logLik),scientific = F)
[1] “0.000000000000593519”

E é isso.

Bem é isso ai, tentando apenas fazer um link entre verossimilhança e árvores, o script vai estar la no repositório recologia, e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail.

Referência:
O curso Computational Molecular Evolution

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
rm(list=ls())
##install.packages("phangorn")
library(phangorn)
 
 
###Gerando um arquivo com sequencias
cat(">A",
    "AAAAA",
    ">B",
    "CAAAA",
    ">C",
    "CGGGG",
    ">D",
    "GGGGG",
file = "exemplo.fasta", sep = "\n")
 
###Lendo sequencias e transformando na classe phyDat
ex.dna <- read.dna("exemplo.fasta", format = "fasta")
dist.dna(ex.dna,model="raw")
 
ex.phyDat<-as.phyDat(ex.dna)
 
 
#####################################
## Topologia 1
#####################################
###Criando uma árvore
arvore1 <- "((A:1,B:1):1,(C:1,D:1):1);"
arvore1.tre <- read.tree(text=arvore1)
 
###Como é árvore
plot(arvore1.tre)
###Calculando o ajuste
ajuste1<-pml(arvore1.tre, data=ex.phyDat,model="JC")
ajuste1
 
#####################################
## Topologia 2
#####################################
###Agora uma topologia ruim
arvore2 <- "((C:1,B:1):1,(A:1,D:1):1);"
arvore2.tre <- read.tree(text=arvore2)
 
###Como é árvore
plot(arvore2.tre)
ajuste2<-pml(arvore2.tre, data=ex.phyDat,model="JC")
 
ex.dna <- read.dna("exemplo.fasta", format = "fasta",as.character=T,as.matrix=F)
 
arvore1.tre$tip.label<-sapply(ex.dna,paste,collapse = "")[match(arvore1.tre$tip.label,names(ex.dna))]
arvore2.tre$tip.label<-sapply(ex.dna,paste,collapse = "")[match(arvore2.tre$tip.label,names(ex.dna))]
 
par(mfrow=c(1,2),mar=c(0,0,2,0))
plot(arvore1.tre,main=paste("Verossimilhança:",round(ajuste1$logLik,2)),type="unrooted")
plot(arvore2.tre,main=paste("Verossimilhança:",round(ajuste2$logLik,2)),type="unrooted")
 
format(exp(ajuste2$logLik),scientific = F)

Métodos comparativos filogenéticos – Observando a correlação

Bem, como sou brasileiro, eu não desisto nunca. Vamos continuar olhando esse esquema de levar em conta a filogenia na análise de dados.

Vamos fazer uma árvore pequena, de 6 espécies, para entender melhor o que está acontecendo.

01

Agora lembrando quando fizemos uma simulação do movimento browniano, em uma dimensão, pensando na variação adquirida por uma característica, vamos fazer isso denovo para essa filogenia.

02

Tudo bem, so que esse é o resultado de uma simulação, usando a função fastBM do pacote phytools, vamos fazer 500 vezes isso, que vamos ver algo bem legal acontecer.

A função fastBM faz isso automático pra gente, basta não salvarmos as o resultado para as espécies internas (espécies ancestrais)…

1
simulacao <- fastBM(arvore,sig2=0.01,internal = TRUE)

Que é a opção padrão da função na verdade (não salvar, para salvar temos que colocar o TRUE), estamos declarando que queríamos essa informação ao falar internal=TRUE, a outra coisa é que agora temos que falar que queremos um n maior, o padrão é 1, então usamos a função assim:

1
saida <- fastBM(arvore, nsim = 500)

Veja que temos uma matriz de 6 linhas e 500 colunas

> str(saida) num [1:6, 1:500] -0.69 -1.438 -0.852 -1.243 -0.294 ... - attr(*, "dimnames")=List of 2 ..$ : chr [1:6] "Sp 3" "Sp 4" "Sp 5" "Sp 6" ... ..$ : NULL

Onde cada coluna é uma simulação, e a linha mostra o resultado final para cada espécie.

> saida[,1:5] [,1] [,2] [,3] [,4] [,5] Sp 3 -0.6898868 -1.3435907 -0.6903719 -0.69457835 0.2401247 Sp 4 -1.4382537 0.1451417 -0.1127161 -1.14320129 -0.2096551 Sp 5 -0.8523258 0.1719962 -0.4519925 -0.60225251 -0.5355870 Sp 6 -1.2430154 -0.6474081 -0.9506413 0.08653215 0.2347673 Sp 2 -0.2935698 -0.2171037 -0.9515205 0.19553651 0.2630972 Sp 1 -3.1330723 0.3918438 1.5757986 0.61583625 -0.6656815

Agora vamos colocar esse resultado numa figura usando a função pairs do R.

03

Veja que para usar o modelo browniano aqui, precisamos de uma árvore com os tamanhos das arestas, o tamanho dos ramos. Nessa figura ai em cima, temos na diagonal o histograma do resultado da 500 simulações de cada espécie, que parecem apresentar uma distribuição normal. Mas quando comparamos os resultado de simulação por simulação entre as espécies, por exemplo da espécie 5 com a 6, veja que temos uma correlação muito alta. Isso porque se você olhar na árvore la em cima, é exatamente as espécies mais próximas, elas “viram” 2 espécies a partir da espécie ancestral delas muito recentemente, e como vimos no modelo browniano, que depende do tempo para acumular variabilidade, não deu tempo quase nenhum nessa árvore de acumular variabilidade.

Agora olha que legal a espécie 1, ela é a espécie mais externa da árvore, essa amiguinha sim, teve todo o tempo evolutivo destinado a essa árvore para diferenciar de todo mundo, logo a correlação dela é muito baixa com todo mundo.

Mas estamos falando de correlação aqui. Vamos dar uma olhada o que estamos calculando com o vcv, que é a função do ape que faz matriz de covariância ou correlação a partir de uma árvore.

Se a gente usar ela…

> vcv(arvore,model = "Brownian", corr = T) Sp 3 Sp 4 Sp 5 Sp 6 Sp 2 Sp 1 Sp 3 1.0000000 0.4586521 0.4586521 0.4586521 0.4474836 0 Sp 4 0.4586521 1.0000000 0.7478113 0.7478113 0.4474836 0 Sp 5 0.4586521 0.7478113 1.0000000 0.8497196 0.4474836 0 Sp 6 0.4586521 0.7478113 0.8497196 1.0000000 0.4474836 0 Sp 2 0.4474836 0.4474836 0.4474836 0.4474836 1.0000000 0 Sp 1 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 1

Veja a correlação entre a espécie 5 com a 6, é a mesma que está na figura la em cima, é a correlação, claro que esses caras não fazem simulações, tem um pulo do gato nas contas que eu não consegui entender ainda, mas estamos caminhando. Mas esse é o caso da correlação, a opção default usa a variância-covariância.

> vcv(arvore,model = "Brownian", corr = F) Sp 3 Sp 4 Sp 5 Sp 6 Sp 2 Sp 1 Sp 3 0.9424450 0.4322544 0.4322544 0.4322544 0.4217286 0.000000 Sp 4 0.4322544 0.9424450 0.7047710 0.7047710 0.4217286 0.000000 Sp 5 0.4322544 0.7047710 0.9424450 0.8008140 0.4217286 0.000000 Sp 6 0.4322544 0.7047710 0.8008140 0.9424450 0.4217286 0.000000 Sp 2 0.4217286 0.4217286 0.4217286 0.4217286 0.9424450 0.000000 Sp 1 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.942445

Que de modo geral vai representar a mesma coisa, mas veja que agora a diagonal não é um. O jeito pra tentar entender como essa conta está sendo feita, é olhar o código da função. (Veja esse post aqui, se boiou no código a seguir)

> vcv function (phy, ...) UseMethod("vcv") 'environment: namespace:ape'

Ra, mas o vcv é uma função genérica, ela tem mais de um método.

> methods(vcv) [1] vcv.corPhyl vcv.phylo

Aqui os métodos usados são do pacote ape, e estamos usando esse vcv.phylo, que tem um código até razoavelmente compreensivo.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
> vcv.phylo
function (phy, model = "Brownian", corr = FALSE, ...) 
{
    if (is.null(phy$edge.length)) 
        stop("the tree has no branch lengths")
    pp <- prop.part(phy)
    phy <- reorder(phy, "postorder")
    n <- length(phy$tip.label)
    e1 <- phy$edge[, 1]
    e2 <- phy$edge[, 2]
    EL <- phy$edge.length
    xx <- numeric(n + phy$Nnode)
    vcv <- matrix(0, n, n)
    for (i in length(e1):1) {
        var.cur.node <- xx[e1[i]]
        xx[e2[i]] <- var.cur.node + EL[i]
        j <- i - 1L
        while (e1[j] == e1[i] && j > 0) {
            left <- if (e2[j] > n) 
                pp[[e2[j] - n]]
            else e2[j]
            right <- if (e2[i] > n) 
                pp[[e2[i] - n]]
            else e2[i]
            vcv[left, right] <- vcv[right, left] <- var.cur.node
            j <- j - 1L
        }
    }
    diag.elts <- 1 + 0:(n - 1) * (n + 1)
    vcv[diag.elts] <- xx[1:n]
    if (corr) {
        Is <- sqrt(1/vcv[diag.elts])
        vcv <- Is * vcv * rep(Is, each = n)
        vcv[diag.elts] <- 1
    }
    dimnames(vcv)[1:2] <- list(phy$tip.label)
    vcv
}
'environment: namespace:ape'

Veja que ele vai pega aonde a estão as partições aqui:

1
2
    e1 <- phy$edge[, 1]
    e2 <- phy$edge[, 2]

E como são as partições aqui:

1
    pp <- prop.part(phy)

Esse prop.part é outra função do pacote ape que vai dizer onde estão as partições nas árvore, aonde, no sentido de quem está com quem, so usar ela que da para pegar o feeling da coisa.

> prop.part(arvore) 1: Sp 3 2: Sp 4 3: Sp 5 4: Sp 6 5: Sp 2 6: Sp 1 ==> 1 time(s):[1] 1 2 3 4 5 6 ==> 1 time(s):[1] 1 2 3 4 5 ==> 1 time(s):[1] 1 2 3 4 ==> 1 time(s):[1] 2 3 4 ==> 1 time(s):[1] 3 4

Veja que primeiro temos uma lista que relação de um número com o nome da espécie, depois temos 1 partição com 3 4, que é a espécie 5 com a 6, as espécies mais próxima na árvore, ai uma vez temos a espécies 5, 6 e 4, e bem, é a árvore ai, acho que deu para entender o que essa função está calculando. Então essa informação mais os tamanhos das arestas é o que é determinado para calcular a variância-covariância ou correlação. Que como tem uma relação linear, não precisa de simulação na verdade.

Pegando a explicação do artigo do artigo citado na função do vcv, esse aqui, pense que estamos pegando o extremo da árvore onde não estão as folhas, e somando a variação de todos os ancestrais até ele. No nosso caso, entre as espécies 5 e 6, existe apenas 1 ancestral, então somamos pouquinha variação, logo a correlação é alta.
entre a espécie 5 e 1, temos que percorrer 4 ancestrais até chegar nela, então a variação é alta, e a correlação muito baixa.

Pra ver essa conta ao invés de so falar, vamos pegar as distâncias:

> arvore$edge.length [1] 0.42172863 0.01052579 0.51019057 0.27251660 0.23767397 0.09604295 [7] 0.14163102 0.14163102 0.52071636 0.94244499

E vamos colocar essa informação na forma de grafo, usando o igraph, pra ficar mais fácil ver, os vértices internos vamos chamar de 1 a 5, e as folhas, que são as espécies, damos os nomes delas mesmo, visualize agora a árvore la de cima como um grafo que tem a distancia mostrada nas arestas (não consegui manipular essas distancias nas arestas em si)

04

O vértice 1 é o ancestral lá de traz, pra espécie 1 chegar nele, ela percorre o 0.94, para a espécie 2 chegar nele, ela percorre 0.52 para chegar no 2 e depois 0.42 pra chegar no 1, que da 0.94, veja que essa é a diagonal da matriz la do vcv, como ja falamos.
A variância-covariância entre a espécie 5 e 6 por exemplo, é o caminho que o ancestral delas percorre até o vértice 1, que vai ser, 0.1 do 5 (ancestral comum entre a espécie 5 e 6) pro 4, 0.27 do 4 pro 3, 0.01 do 3 pro 2 e 0.42 do 2 pro 1, que da 0.8, o valor que ta na matriz[4,3].

Acho que com isso ja da para começar a entender bem como está sendo estabelecida essa correlação entre as espécies. No artigo do Garland e Ives, tem o uso dele na regressão, ou no método de mínimos quadrados (GLS), que vai ficar para o próximo post.

E é isso ai, agora estamos um pouco mais próximos de considerar a evolução em análises de dados, o script vai estar la no repositório recologia, e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo, ou mande um e-mail e até mais.

Referência:
Garland, T. Jr. and Ives, A. R. (2000) Using the past to predict the present: confidence intervals for regression equations in phylogenetic comparative methods. American Naturalist, 155, 346–364.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
set.seed(123)
library(phytools)
# Simulando uma arvore de 5 espécies
arvore <- pbtree(n = 6)
arvore$tip.label <- paste("Sp",sub("t","",arvore$tip.label))
jpeg("01.jpg")
plot(arvore)
dev.off()
 
## 1 simulação de movimento browniano
simulacao <- fastBM(arvore,sig2=0.01,internal = TRUE)
jpeg("02.jpg")
phenogram(arvore,simulacao, spread.labels = TRUE, spread.cost = c(1, 0))
dev.off()
 
# 500 simulações de movimento browniano
saida <- fastBM(arvore, nsim = 500)
 
str(saida)
saida[,1:5]
 
#funções auxiliares para o pairs
#para plotar um histograma na diagonal do gráfico
panel.hist <- function(x, ...) {
    usr <- par("usr"); on.exit(par(usr))
    par(usr = c(usr[1:2], 0, 1.5) )
    h <- hist(x, plot = FALSE)
    breaks <- h$breaks; nB <- length(breaks)
    y <- h$counts; y <- y/max(y)
    rect(breaks[-nB], 0, breaks[-1], y, col = "cyan", ...)
}
# Para colocar as correlações, com tamanho do texto proporcional ao valor.
panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...) {
    usr <- par("usr"); on.exit(par(usr))
    par(usr = c(0, 1, 0, 1))
    r <- abs(cor(x, y))
    txt <- format(c(r, 0.123456789), digits = digits)[1]
    txt <- paste0(prefix, txt)
    if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
    text(0.5, 0.5, txt, cex = cex.cor * r + 0.5)
}
 
jpeg("03.jpg")
pairs(t(saida),diag.panel = panel.hist,lower.panel = panel.smooth,upper.panel = panel.cor)
dev.off()
 
#vcv para correlação
vcv(arvore,model = "Brownian", corr = T)
 
#vcv para variancia-covariancia
vcv(arvore,model = "Brownian", corr = F)
 
#olhando o código do método vcv
vcv
methods(vcv)
 
vcv.phylo
 
#usando o prop.part
prop.part(arvore)
 
#Grafo
library(igraph)
dados.g <-graph.data.frame(arvore$edge, directed=F)
jpeg("04.jpg")
plot.igraph(dados.g,edge.label=round(arvore$edge.length,2),vertex.label=c(1:5,arvore$tip.label))
dev.off()

Crescimento Populacional e a Seleção Natural

A teoria da evolução é talvez a teoria mais importante para a biologia. Ela unificou muitos ramos científicos atribuídos a biologia que antes dela eram praticamente áreas distintas, como fisiologia, taxonomia, comportamento, etc.

Temos a frase celebre do nosso colega Theodosius Dobzhansky.

Nada em biologia faz sentido se não for à luz da evolução

Mas será que a gente realmente entende como a evolução funciona?

Quero dizer, provavelmente, todo mundo tem uma noção básica. O que a gente aprende na escola é tranquilo, mas o problema é que muitas vezes nos lembramos do exemplo básico e fazemos generalizações muito amplas, ou pior ainda, generalizações curtas e não conseguimos ver o processo como ele realmente ocorre.

Na verdade a teoria da evolução é muito mais do que vamos ver aqui, mas vamos olhar como funciona o negócio de variação no número de descendentes produzidos. Quem não lembra do exemplo da mariposa e da revolução industrial. A Biston betularia. Essa ai é aquela mariposa que existia de dois tipo, ou dois fenótipos, como preferir.

Tinhamos a versão, ou fenótipo mais claro.
claro

E a versão, ou fenótipo, mais escuro.
escuro

Ai antes da poluição, a mais clarinha era comum, e a escura era muito rara, com o acumulo de fuligem nos troncos das arvores, a clarinha começou a não ser muito boa para se camuflar, a escura, da cor da fuligem, começou a se dar bem e o cenário se inverteu, com a clarinhas passando a ser a mais rara e a escura sendo a mais comum.

Daqui pra frente vamos chamar de fenótipos sempre. Se a gente fosse la no campo e conta-se todas as mariposas existentes, poderíamos dizer que existe uma quantidade N de mariposas, essa quantidade sendo independente do fenótipo, que é a cor para as mariposas.

E elas vão vivendo, se reproduzindo a uma certa taxa, e bem isso nos faz lembrar da equação de crescimento populacional.

N_t = N_0 \cdot r^t

Lembre-se que estamos falando do crescimento discreto ai em cima, não confunda com a formula N_t = N_0 \cdot e^{r \cdot t} , a gente já discutiu sobre crescimento discreto versus crescimento contínuo aqui, de uma olhada para relembrar e ver a diferença.

Mas continuando, então temos nosso crescimento populacional, vamos ficar com ele por enquanto para facilitar o raciocinio, se começar a meter aquele número e de Euler vai ficar tudo mais difícil.

N_t = N_0 \cdot r^t

So que essa equação é valida para toda a população dessa determinada espécie, todos os indivíduos de uma região que compõem uma população de uma certa espécie, nossa mariposinha por exemplo.

Vamos pensar nos fenótipos claro e escuro como fenótipo A e fenótipo a, cada um deles então vai ter uma frequência na população, que chamaremos de f_A e f_a e cada um deles vai ter sua taxa de crescimento R, r_A e r_a, que também é conhecido como o fitness absoluto (Absolute fitness), que a princípio vamos considerar como igual a taxa de crescimento da população como um todo, escrevendo isso de forma mais clara, vamos assumir que r_A = r_a = R.

Primeiro, vamos pensar que o número inicial de indivíduos A seria igual a sua frequência no total inicial da população.

N_{A,0} = N_0 \cdot f_A

Certo, se a população inicial é N_0, sabemos que uma fração dela é de mariposas claras, a fração é f_A, por isso se multiplicar o total pela fração N_0 \cdot f_A da N_{A,0}

O mesmo vai ser valido para a população de a,

N_{a,0} = N_0 \cdot f_a

Lembrando que f_A + f_a = 1, tem que ser o total da população.

Vamos em frente, A e a são a população total dessa espécie.

N_0 = N_0 \cdot f_A + N_0 \cdot f_a

Sem muita mágica até aqui.
Certo, agora como o crescimento vai funcionar aqui?
Muito simples, para a próxima geração (vamos pensar apenas na diferença de uma geração para a próxima) temos a quantidade inicial vezes a taxa de crescimento, so que por cada um dos fenótipos.

Para A temos

N_{A,1} = N_0 \cdot f_A \cdot r_A

E para a temos

N_{a,1} = N_0 \cdot f_a \cdot r_a

Mas como r_A = r_a = r, o fitness dos dois são iguais, então poderíamos simplesmente dizer.

N_{A,1} = N_0 \cdot f_A \cdot r

e

N_{a,1} = N_0 \cdot f_a \cdot r

Então podemos ver que a população total no tempo 1 N_1 vai ser

N_1 = N_0 \cdot f_A \cdot r + N_0 \cdot f_a \cdot r

Veja que o termo N_0, que é a população total inicial e o termo rocorrem nas duas multiplicações, então podemos simplificar isso.

N_1 = N_0 \cdot r \cdot (f_A + f_a)

Lembre-se que f_A + f_a = 1 já que são frequência, e esses dois A e a tem que ser o todo.

Então apos uma geração a frequência de A será a quantidade de indivíduos A pelo total de indivíduos da população, a frequência f_A será o seguinte.

{f_{A,1}} = {{N_0 \cdot f_{A,0} \cdot r}\over{N_0 \cdot r}}

Como é tudo uma multiplicação podemos cortar os N_0 e r e assim ficamos somente com.

f_{A,1}= f_{A,0}

A frequência na próxima geração não mudou nada.

Mas é isso que esperamos não? Porque nesse exemplo r_A = r_a = r, ou seja os dois fenótipos eram igualmente bons, mas e se r_A \ne r_a, um deles pode ser melhor que o outro, sei la r_A > r_a ou r_A < r_a, mas ambos os casos r_A \ne r_a, como fica esse mesmo esquema?

A principio, a próxima geração vai ter a mesma formula do crescimento.

N_{A,1} = N_0 \cdot f_A \cdot r_A

E para a temos

N_{a,1} = N_0 \cdot f_a \cdot r_a,

Ou ainda, generalizando para o tempo t, temos

N_{A,t} = N_0 \cdot f_A \cdot {r_A}^t

E para a temos

N_{a,t} = N_0 \cdot f_a \cdot {r_a}^t

Ou seja a população N_t será

N_{t} = N_0 \cdot f_A \cdot {r_A}^t + N_0 \cdot f_a \cdot {r_a}^t

Agora vamos fazer aquela conta denovo, após t gerações, a frequência f_A vai ser a quantidade de individuos A sobre o total N

{f_{A,t}} = {N_0 \cdot f_A \cdot {r_A}^t \over{N_0 \cdot f_A \cdot {r_A}^t + N_0 \cdot f_a \cdot {r_a}^t}}

Agora nos podemos avaliar a frequência de a depois de t gerações assim como fizemos antes, mas podemos simplificar um pouco essa equação. Podemos começar tirando esse N_0 em cima e em baixo na divisão.

A frequência f_A vai ser o total dela sobre o tamanho total da população.
{f_{A,t}} = {f_A \cdot {r_A}^t \over{f_A \cdot {r_A}^t + f_a \cdot {r_a}^t}}

Agora podemos dividir em cima e em baixo por r_A e ficamos com o seguinte

{f_{A,t}} = {{f_A}\over{{f_A} + {f_a \cdot {{r_a}^t\over{{r_A}^t}}}}}

E finalmente temos a nossa equação para acompanhar como a frequência gênica muda ao longo das gerações.
Se a gente pensar um pouco, quando a razão {r_a}^t\over{{r_A}^t} é 1, ambos os fenótipos são igualmente bons, essa razão da um, e teremos  {f_A}\over{f_A + f_a \cdot 1}, como  f_A + f_a = 1, teremos  f_A dividido por 1 e nada nunca muda. Agora qualquer desvio na razão {r_a}^t\over{{r_A}^t} para mais ou para menos que 1, teremos mudanças.

Mas de qualquer forma, vamos jogar alguns números e fazer uma simulação no R.

Precisamos estabelecer algumas condições iniciais, e pensar por quantas gerações queremos simular.

# ************************************************** # Condições iniciais # ************************************************** N0 <- 1000 # População inicial rate_A <- 1.2 # Taxa de crescimento do alelo "A" rate_a <- 1.2 # Taxa de crescimento do alelo "a" fA <- 0.3 # Frequência do alelo "A" max_gen <- 20 # Número de gerações a simular

Veja que baseado naquelas quantidades, o resto dos ingredientes são derivados

# ************************************************** # Calculando variáveis derivadas # ************************************************** fa <- 1.0 - fA # Frequência do alelo "a" NA_0 <- N0 * fA # População inicial do alelo "A" Na_0 <- N0 * fa # População inicial do alelo "a"

Mas vamos la, criamos uma função para com aquelas condições iniciais e usando as equações que estabelecemos acima para ver o que acontece.

evosimu<-function(N0,rate_A,rate_a,fA,max_gen) {
    fa <-1.0-fA
    NA_0<- N0 * fA
    Na_0<- N0 * fa
 
    resultado<-matrix(NA,ncol=5,nrow=max_gen+1)
    colnames(resultado)<-c("Nt","NA_t","Na_t","fA_t","fa_t")
    resultado[1,]<-c(N0,NA_0,Na_0,fA,fa)
 
    for(t in 2:(max_gen+1)) {
        resultado[t,"NA_t"] = NA_0 * (rate_A ^ t)
        resultado[t,"Na_t"] = Na_0 * (rate_a ^ t)
        resultado[t,"Nt"] =  resultado[t,"NA_t"] + resultado[t,"Na_t"]
        resultado[t,"fA_t"] =  resultado[t,"NA_t"] / resultado[t,"Nt"]
        resultado[t,"fa_t"] = resultado[t,"Na_t"] / resultado[t,"Nt"]
    }
    return(resultado)
}

E testamos nossa funçãozinha.

> evosimu(N0=N0,rate_A=rate_A,rate_a=rate_a,fA=fA,max_gen=10) Nt NA_t Na_t fA_t fa_t [1,] 1000.000 300.0000 700.000 0.3 0.7 [2,] 1440.000 432.0000 1008.000 0.3 0.7 [3,] 1728.000 518.4000 1209.600 0.3 0.7 [4,] 2073.600 622.0800 1451.520 0.3 0.7 [5,] 2488.320 746.4960 1741.824 0.3 0.7 [6,] 2985.984 895.7952 2090.189 0.3 0.7 [7,] 3583.181 1074.9542 2508.227 0.3 0.7 [8,] 4299.817 1289.9451 3009.872 0.3 0.7 [9,] 5159.780 1547.9341 3611.846 0.3 0.7 [10,] 6191.736 1857.5209 4334.215 0.3 0.7 [11,] 7430.084 2229.0251 5201.059 0.3 0.7

Temos como entradas as quantidades iniciais, que repetimos na primeira linha, a geração zero, e temos na saída o tamanho da população total como primeira coluna, como ela vai mudando sendo cada geração uma linha, como cada fenótipo muda, bem como suas frequências ao longo das gerações, mas como olhar para uma matriz de números é meio chato, vamos representar esses dados em dois gráficos, um de como está a população e outro de como estão as frequências.

01

Certo, bem simples não, a população está crescendo exponencialmente, e as frequências genicas estão constantes, como previmos, porque as taxas de crescimento iniciais que colocamos para essa simulação são iguais.

rate_A <- 1.2 # Taxa de crescimento do alelo "A" rate_a <- 1.2 # Taxa de crescimento do alelo "a"

Agora lembre-se que nem tudo na vida são flores, se tivermos taxas de crescimento abaixo de 1, a população estará diminuindo.

Vamos pegar esse caso aqui.

#N0 = 1000 rate_A = 0.7 rate_a = 0.7 fA = 0.3

02

Mas ainda temos as frequências dos fenótipos iguais, aqui ta todo mundo se ferrando igual. Agora quando a gente diz que está havendo evolução é quando as frequência começam a mudar, quando aparece um novo mutante que derrepente se da melhor, se alguém tem a taxa de crescimento maior que sua contraparte, vemos o seguinte.

#N0 = 1000 rate_A = 2.0 rate_a = 1.5 fA = 0.02

03

Veja que aqui, nenhum dos fenótipos está necessáriamente se ferrando, só temos um fenótipo que é melhor que o outro, mas ninguém está diminuindo em quantidade ao longo do tempo, a população está crescendo como um todo, todos os fenótipos estão crescendo, apenas um mais rápido que o outro, mas isso faz com que a frequência dos fenótipos mude. Quando eu digo que as vezes a evolução pode não ser tão simples, porque como no exemplo das mariposas, a gente não pensa nesse caso assim, a gente sempre pensa que se alguém é pior, ele tem que estar se ferrando, como no próximo exemplo

#N0 = 1000 rate_A = 1.2 rate_a = 0.9 fA = 0.02

04

Essa é a situação que a gente tem na cabeça para as mariposas, pelo menos a primeira que vinha na minha, se a gente pensar na linha azul como as mariposas brancas e vermelho as mariposas escuras, a gente quer que alguém seja muito comum e alguém raro, ai esse cara mais comum, as mariposas brancas, começam a se dar mal por causa da revolução industrial, as mariposas escuras, que antes eram raras e ruim, começam a se dar bem, e temos uma inversão de quem é mais raro e quem é mais comum, mas vemos que a população decaiu um pouco e depois voltou a crescer. Mas como vimos, a evolução pode acontecer mesmo quando todo mundo está se dando bem, com a população crescendo, e além disso ainda, podemos pensar em outra situação.

#N0 = 10000 rate_A = 0.8 rate_a = 0.6 fA = 0.02

05

E é isso ai, ambos os fenótipos estão se ferrando, ambos estão diminuindo sua frequência. Cada vez que você olha a mata, você ve menos e menos daquela espécia que você observa, mas isso não quer dizer que a evolução não possa estar agindo ali. Veja que aqui não é quem se ferra e quem se da bem, aqui a situação é quem se ferra menos, quem pode perdurar aos tempos difíceis, talvez essa espécie seja extinta, mas talvez não, e quando os tempos mudarem, ela pode voltar a crescer, mas as frequências dos fenótipos estão alteradas já na próxima largada para o crescimento populacional.

Mas, não é a toa que nos humanos somos péssimos geradores de números aleatórios, pelo menos eu, sempre tendia a pensar em evolução, no caso de uma espécie que tem um fenótipo se tornando mais raro, porque a população dele está diminuindo, e outro se tornando mais comum porque a população dele está aumentando, mas esteja a população estacionada, aumentando ou diminuindo, as frequências de alelos, fenótipos, snp, como queria pensar, podem estar mudando, e todos esses cenários tem que estar na nossas cabeças.

So que a primeira coisa que vem na minha cabeça é um fenótipo se dando bem e outro se dando mal. Mas existem muitos mais possibilidades. Então é preciso se policiar para não tomar conclusões precipitadas, e conduzir uma análise de dados que de possibilidade de todos os casos serem avaliados corretamente.

Alias, para teste em sistemas mais complexos, sempre esses teste envolvem simulações com parâmetros aleatórios, porque muitas vezes o que podemos pensar não contempla todas as possibilidades, pensa em quantas árvores podemos desenhar com n espécies, eu nunca pensei que existissem tantas possibilidades assim, e é fácil se perder nas contas. Mas ao tentar algumas simulações, passar tudo para um modelo e ir tentando várias possibilidades de parâmetros, podemos talvez ver aquilo que nunca tinha passado pela nossa cabeça.

Bem é isso ai, as equações la do começo do post eu peguei de um curso do coursera, o Computational Molecular Evolution, ótimo curso alias, uma excelente seleção de tópicos para a ementa e conceitos bem complicados são explicados de maneira bem simples, pena que a maior parte pratica é usando o paup, mas agora é ir adaptando tudo pro R com as ferramentas que tivermos a mãos. Scripts no repositório recologia e aqui em baixo, e até o próximo post.

# **************************************************
# Condições iniciais
# **************************************************
 
N0      <- 1000 # População inicial
rate_A  <- 1.2  # Taxa de crescimento do alelo "A"
rate_a  <- 1.2  # Taxa de crescimento do alelo "a"
fA      <- 0.3  # Frequência do alelo "A"
max_gen <- 20   # Número de gerações a simular
 
# **************************************************
# Calculando variáveis derivadas
# **************************************************
 
fa   <- 1.0 - fA  # Frequência do alelo "a"
NA_0 <- N0 * fA   # População inicial do alelo "A"
Na_0 <- N0 * fa   # Polulação inicial do alelo "a"
 
# **************************************************
# Simulação
# **************************************************
 
evosimu<-function(N0,rate_A,rate_a,fA,max_gen) {
    fa <-1.0-fA
    NA_0<- N0 * fA
    Na_0<- N0 * fa
 
    resultado<-matrix(NA,ncol=5,nrow=max_gen+1)
    colnames(resultado)<-c("Nt","NA_t","Na_t","fA_t","fa_t")
    resultado[1,]<-c(N0,NA_0,Na_0,fA,fa)
 
    for(t in 2:(max_gen+1)) {
        resultado[t,"NA_t"] = NA_0 * (rate_A ^ t)
        resultado[t,"Na_t"] = Na_0 * (rate_a ^ t)
        resultado[t,"Nt"] =  resultado[t,"NA_t"] + resultado[t,"Na_t"]
        resultado[t,"fA_t"] =  resultado[t,"NA_t"] / resultado[t,"Nt"]
        resultado[t,"fa_t"] = resultado[t,"Na_t"] / resultado[t,"Nt"]
    }
    return(resultado)
}
 
evosimu(N0=N0,rate_A=rate_A,rate_a=rate_a,fA=fA,max_gen=10)
 
saida<-evosimu(N0=N0,rate_A=rate_A,rate_a=rate_a,fA=fA,max_gen=max_gen)
 
layout(matrix(c(1,2),nrow=2,ncol=1))
plot(0:max_gen,saida[,"Nt"],frame=F,xlab="Geração",ylab="População",type="l",ylim=c(0,max(saida[,"Nt"])))
points(0:max_gen,saida[,"NA_t"],type="l",col="red")
points(0:max_gen,saida[,"Na_t"],type="l",col="blue")
legend("topleft",col=c("black","red","blue"),lty=1,legend=c("Total","A","a"),bty="n")
plot(0:max_gen,saida[,"fA_t"],frame=F,xlab="Geração",ylab="Frequência Gênica",type="l",ylim=c(0,1),col="red")
points(0:max_gen,saida[,"fa_t"],type="l",col="blue")
 
#N0 = 50     rate_A = 1.2  rate_a = 1.2  fA = 0.3
saida<-evosimu(N0=50,rate_A=1.2,rate_a=1.2,fA=0.3,max_gen=20)
 
layout(matrix(c(1,2),nrow=2,ncol=1))
plot(0:max_gen,saida[,"Nt"],frame=F,xlab="Geração",ylab="População",type="l",ylim=c(0,max(saida[,"Nt"])))
points(0:max_gen,saida[,"NA_t"],type="l",col="red")
points(0:max_gen,saida[,"Na_t"],type="l",col="blue")
legend("topleft",col=c("black","red","blue"),lty=1,legend=c("Total","A","a"),bty="n")
plot(0:max_gen,saida[,"fA_t"],frame=F,xlab="Geração",ylab="Frequência Gênica",type="l",ylim=c(0,1),col="red")
points(0:max_gen,saida[,"fa_t"],type="l",col="blue")
 
#N0 = 1000   rate_A = 0.7  rate_a = 0.7  fA = 0.3
saida<-evosimu(N0=100,rate_A=0.7,rate_a=0.7,fA=0.3,max_gen=max_gen)
 
layout(matrix(c(1,2),nrow=2,ncol=1))
plot(0:max_gen,saida[,"Nt"],frame=F,xlab="Geração",ylab="População",type="l",ylim=c(0,max(saida[,"Nt"])))
points(0:max_gen,saida[,"NA_t"],type="l",col="red")
points(0:max_gen,saida[,"Na_t"],type="l",col="blue")
legend("topright",col=c("black","red","blue"),lty=1,legend=c("Total","A","a"),bty="n")
plot(0:max_gen,saida[,"fA_t"],frame=F,xlab="Geração",ylab="Frequência Gênica",type="l",ylim=c(0,1),col="red")
points(0:max_gen,saida[,"fa_t"],type="l",col="blue")
 
#N0 = 1000   rate_A = 2.0  rate_a = 1.5  fA = 0.02
saida<-evosimu(N0=1000,rate_A=2.0,rate_a=1.5,fA=0.02,max_gen=max_gen)
 
layout(matrix(c(1,2),nrow=2,ncol=1))
plot(0:max_gen,saida[,"Nt"],frame=F,xlab="Geração",ylab="População",type="l",ylim=c(0,max(saida[,"Nt"])))
points(0:max_gen,saida[,"NA_t"],type="l",col="red")
points(0:max_gen,saida[,"Na_t"],type="l",col="blue")
legend("topleft",col=c("black","red","blue"),lty=1,legend=c("Total","A","a"),bty="n")
plot(0:max_gen,saida[,"fA_t"],frame=F,xlab="Geração",ylab="Frequência Gênica",type="l",ylim=c(0,1),col="red")
points(0:max_gen,saida[,"fa_t"],type="l",col="blue")
 
#N0 = 1000   rate_A = 1.2  rate_a = 0.9  fA = 0.02
saida<-evosimu(N0=1000,rate_A=1.2,rate_a=0.9,fA=0.02,max_gen=max_gen)
 
layout(matrix(c(1,2),nrow=2,ncol=1))
plot(0:max_gen,saida[,"Nt"],frame=F,xlab="Geração",ylab="População",type="l",ylim=c(0,max(saida[,"Nt"])))
points(0:max_gen,saida[,"NA_t"],type="l",col="red")
points(0:max_gen,saida[,"Na_t"],type="l",col="blue")
legend("topleft",col=c("black","red","blue"),lty=1,legend=c("Total","A","a"),bty="n")
plot(0:max_gen,saida[,"fA_t"],frame=F,xlab="Geração",ylab="Frequência Gênica",type="l",ylim=c(0,1),col="red")
points(0:max_gen,saida[,"fa_t"],type="l",col="blue")
 
#N0 = 10000  rate_A = 0.8  rate_a = 0.6  fA = 0.02
saida<-evosimu(N0=10000,rate_A=0.8,rate_a=0.6,fA=0.02,max_gen=max_gen)
 
layout(matrix(c(1,2),nrow=2,ncol=1))
plot(0:max_gen,saida[,"Nt"],frame=F,xlab="Geração",ylab="População",type="l",ylim=c(0,max(saida[,"Nt"])))
points(0:max_gen,saida[,"NA_t"],type="l",col="red")
points(0:max_gen,saida[,"Na_t"],type="l",col="blue")
legend("topleft",col=c("black","red","blue"),lty=1,legend=c("Total","A","a"),bty="n")
plot(0:max_gen,saida[,"fA_t"],frame=F,xlab="Geração",ylab="Frequência Gênica",type="l",ylim=c(0,1),col="red")
points(0:max_gen,saida[,"fa_t"],type="l",col="blue")

Problema – Rosalind – Counting Phylogenetic Ancestors – INOD

Esse problema é bem simples, a pergunta consiste em, dado um número n inteiro, que é o número de espécies que vamos avaliar, quantos ancestrais eles tem, quantos ancestrais precisamos para fazer uma árvore filogenética ligando eles.

A resposta aqui é bem simples, e nem precisa de um programa na verdade. A resposta é simplesmente n-2.

Mas para chegar até essa resposta, a gente tem que olhar algo legal, as árvores binárias sem raiz.

Essa é a estrutura de dados que é gerada quando a gente vai fazer a árvore genealógica de uma família por exemplo. Que é um tipo especial de árvore binária. Poxa, como teorias dos grafos é uma coisa legal, ja esbarramos neles aqui e aqui, mas essas eram estruturas e ideias bem diferentes da daqui, usando grafos bipartidos.

781px-Waldburg_Ahnentafel

Pela imagem , que é da família de Sigmund Christoph von Waldburg-Zeil-Trauchburg de 1776, da para ter a ideia, que mesmo a gente pensando isso hoje, essa estrutura de dados já esta por ai a muito tempo.

02

Em teoria dos grafos, esse grafo, que é uma estrutura de dados de árvore, é caracterizado pelo fato de cada vértice, ou node, ou bolinha se preferir, ter no máximo dois “filhos”, então ele vai ter um grau 3 no máximo, ou seja 3 arestas ou edge que saem da uma bolinha, no máximo, poderiam ser 2 como numa folha porém. Veja que cada bolinha tem no máximo 3 conexões. a gente fez um problema ja comentando sobre isso aqui.

Agora olhando o texto do wikipedia sobre árvores binárias sem raiz, a própria definição diz que para n folhas da árvore, que na biologia são as espécies, haverão n-2 vértices internos, ou ancestrais, e somente assim é possível construir essas árvores. Por isso o exemplo diz que para 4 espécies, precisamos de 2 ancestrais, e assim vai, n-2 vértices internos na árvore.

Agora outra coisa legal, que pelo menos eu sempre tive dificuldade de entender ou pensar, é como as pessoas veem coisas diferentes na mesma figura, comparando por exemplo computação e biologia.

Por exemplo, numa aula de grafos, uma árvore binaria sem raiz vai ser apresentada assim, talvez, num sei, eu imagino, apesar que nunca tive aulas sobres grafos:

01

Numa filogenia, a terminação ali é uma espécie, árvores binárias sem raiz não tem ciclos, ou seja, todas as pontas terminam numa extremidade solta, e além disse se você pensa na união dos vértices ali, todo mundo ta ligado a um ancestral.

O problema é que a primeira vista, isso não se parasse nada com o que vimos nesse post aqui não?

03

A gente quer ver uma filogenia mais desse jeito, uma árvore desse jeito.

Mas a partir daquela figura:

01

A gente tem que pensar como se pegássemos uma parte…

02

Ai a gente puxa essa parte…

03

Dai se apontarmos todas as espécies para o mesmo lado…

04

Fica com uma cara muito melhor de filogenia, ou árvore filogenética não? Da para dar mais uma melhoradinha ainda…

05

E olha ai, agora parece que da para entender, ver que é tudo igual, so depende da perspectiva. Depois de entender que o jeito que a gente visualiza as coisas é relativo, fica mais fácil entender da onde diabos vem o maximum likelihood para desenhar árvores, ja que a gente sempre pensa em mínimos quadrados da perspectiva de uma anova, como aqui onde somamos as distancias da média, mas agora que temos folhas, que são espécies ligadas por risquinhos, se a gente pensar que essa distância entre cada espécie, o caminho de uma folha a outra é o quão igual elas devem ser, da para começar a imaginar como se coloca o conceito de verossimilhança nesse problema. E além disso, se da para usar verissimilhança, da para ter um “pre-conceito”, uma noção de como a árvore deve ser e usar isso de prior para as distâncias das espécies, numa perspectiva bayesiana, mas isso fica para um próximo post.