Rosalind – Genome Assembly as Shortest Superstring

Após mais algum tempo sem nada, vamos ver mais um problema do Rosalind sobre montar genomas. Dado que a tecnologia mais em conta agora é o sequenciamento High-throughput, a gente sempre tem um conjunto de reads para montar.

Reads são esses pedacinhos de dna, que vem da molécula original, eles tem partes repetidas que permitem a gente “montar” de volta a molécula original. Bem só lendo como funciona o processo para entender melhor, mas aqui a gente tem uma versão bem simples para começar a entender o que tem que ser feito.

Mas então, a gente tem um conjunto de vários pedacinhos de dna (reads) no formato fasta derivados de uma única molécula e temos que montar a menor sequência de dna possível (Shortest Superstring) que pode ter derivado todos esses pedacinhos.

Agora é precisa ver o seguinte, para esses dados, é garantido que existe uma única sequência original, não tem como montar de mais de um jeito, e existe uma sobreposição maior que 50% do reads, ou seja, a gente tem que encaixar no mínimo metade de um read no outro.

Vamos olhar os dados, que talvez facilite entender melhor o que pode ser feito. Esse é o conjunto de entrada.

>Rosalind_56 ATTAGACCTG >Rosalind_57 CCTGCCGGAA >Rosalind_58 AGACCTGCCG >Rosalind_59 GCCGGAATAC

E queremos fazer a menor sequência que tem todos esses reads como uma substring, então a gente pode começar a tentar encaixar, pegando os dois primeiros…

ATTAGACCTG xxx|x|xxx CCTGCCGGAA

Mas a gente pode tentar encaixar pelo outro lado também.

ATTAGACCTG x|xxxxxxx CCTGCCGGAA

A gente começa movendo só uma posição para tentar fazer o menor substring, ai como não deu, a gente diminuir em dois o tamanho do encaixe,

ATTAGACCTG xxxx||x| CCTGCCGGAA

Veja que o read 56 encaixa no 58 por 7 bases, e como os reads tem 10 bases, é mais da metade.

ATTAGACCTG ||||||| AGACCTGCCG

E encaixando tudo, a gente monta a sequência inicial

>Rosalind_56 ATTAGACCTG >Rosalind_58 |||AGACCTGCCG >Rosalind_57 ||||||CCTGCCGGAA >Rosalind_59 |||||||||GCCGGAATAC ATTAGACCTGCCGGAATAC

A gente pode fazer uma função recursiva no python.

def long(reads, genoma=''):
	if len(reads) == 0:
		return genoma
	if len(genoma)==0:
		genoma = reads.pop(0)
		return long(reads, genoma)  
	for i in range(len(reads)):
		a = reads[i]
		l = len(a)
		for p in range(l / 2):			
			q = l - p
			if genoma.startswith(a[p:]):
				reads.pop(i)
				return long(reads, a[:p] + genoma)
			if genoma.endswith(a[:q]):
				reads.pop(i)
				return long(reads, genoma + a[q:])

A função é assim, primeiro o nosso caso de parada, que é quando não temos mais reads para montar

	if len(reads) == 0:
		return genoma

Quando estamos começando, só tem um jeito de montar, logo a gente pode começar do primeiro read

	if len(genoma)==0:
		genoma = reads.pop(0)
		return long(reads, genoma)

Ai , para todos os reads

	for i in range(len(reads)):
		a = reads[i]
		l = len(a)

Para pelo menos metade do tamanho dos reads

		for p in range(l / 2):

A gente começa diminuindo uma base para encaixar, depois 2, e assim por diante

			q = l - p

Ai tentamos encaixar na frente

			if genoma.startswith(a[p:]):
				reads.pop(i)
				return long(reads, a[:p] + genoma)

E encaixar atrás

			if genoma.endswith(a[:q]):
				reads.pop(i)
				return long(reads, genoma + a[q:])

Ai veja que só existe uma maneira de montar, logo a gente vai entrar nessa recursão somente uma vez, porque uma vez que o read seja unido, a gente tira ele da lista de reads e já junta ao genoma, e outra coisa, é que so tem um read para cada ponta.

Bem é isso ai, 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.

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
#!/usr/bin/env python
'''
 
'''
 
def long(reads, genoma=''):
	if len(reads) == 0:
		return genoma
	if len(genoma)==0:
		genoma = reads.pop(0)
		return long(reads, genoma) 
	for i in range(len(reads)):
		a = reads[i]
		l = len(a)
		for p in range(l / 2):			
			q = l - p
 
			if genoma.startswith(a[p:]):
				reads.pop(i)
				return long(reads, a[:p] + genoma)
 
			if genoma.endswith(a[:q]):
				reads.pop(i)
				return long(reads, genoma + a[q:])
if __name__ == "__main__":
 
	arquivo = open("/home/augusto/Downloads/rosalind_long.txt")
	sequencia=[]
	i=-1
	for linha in arquivo:
		if linha.find('>') != -1:		
			i+=1
			sequencia.append('')
		else:		
			sequencia[i]=sequencia[i]+linha.rstrip()
	print long(sequencia)

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)

Teoria da Coalescência

Teoria da coalescência, essa faz um tempo que estou tentando entender. Tai algo que você le, parece que da pra entender o que as pessoas falam, mas aquelas fórmulas, é complicado de fazer sentido delas, ainda mais quando vem umas umas figuras assim:

fig

Essa figura é do livro Bayesian Evolutionary Analysis with Beast, mas vamos la, se a gente olhar o o artigo no wikipedia, ou em qualquer livro, sempre a gente vai cair na seguinte equação:

Pr\{t\}=(\frac{1}{N}) \cdot  (1-\frac{1}{N})^{t-1}

Essa função da a chance para coalescer no tempo t, de dois indivíduos dentro de uma população de tamanho N.

No wikipedia está:

Pr\{t\}=(\frac{1}{2N}) \cdot  (1-\frac{1}{2N})^{t-1}

Mas vamos ignorar esse 2 por enquanto, pensando em bactérias, onde não precisamos de pai e mãe, todo mundo é haploide, então a fórmula que vamos usar vai ser essa:

Pr\{t\}=(\frac{1}{N}) \cdot  (1-\frac{1}{N})^{t-1}

Primeiro, vamos passar isso para o R, criando uma função.

1
2
3
coalescent<- function(n,t) {
    return((1/n) * (1-(1/4))^(t-1))
}

Agora o que diabos é calculado ae? Bem é calculado a chance de coalescer, que nada mais é, que a chance de ter mães em comum. Pense em dois indivíduos, qual a chance desses coisa caras terem a chance de ter a mesma mãe, se na geração anterior só tinha um indivíduo?
Podemos usar nossa função, onde primeiro falamos qual o tamanho da população e depois a quantas gerações a atrás queremos saber a probabilidade.

> coalescent(1,1) [1] 1

A resposta é 1, porque como o tamanho da população é 1, só tem como esse único individuo ser a mãe desses dois indivíduos, oras bolas.

Agora vamos pensar numa população de 2 indivíduos, veja que o N ser dois significa que a geração das mães tinha 2 indivíduos, aqui estamos sempre pensando para dois filhos.

> coalescent(2,1) [1] 0.5

Isso da meio, isso porque ou você é filho um dos indivíduos ou do outro? Esse N dividindo ali faz parecer algo desse tipo, mas não é bem assim. Vamos pensar aqui, existe duas possíveis mães, bactérias originarias. Vamos chamá-las de “a” e “b”.

> mae<-letters[1:2] > mae [1] "a" "b"

Agora veja que se temos duas mães, duas bactérias mãe, e temos dois filhos, podemos explorar as possibilidades, que são 4

> possibilidades<-expand.grid(individuo1=mae,individuo2=mae) > possibilidades individuo1 individuo2 1 a a 2 b a 3 a b 4 b b

Agora coalescer, é ter a mesma origem, a mesma bactéria mãe,

> possibilidades$individuo1==possibilidades$individuo2 [1] TRUE FALSE FALSE TRUE

E tem dois casos em que isso acontece, quando “a” é a mãe dos dois filhos ou “b” é a mãe dos dois filhos. Só que isso vem de 4 possíveis combinações como vimos ali em cima, ou podemos representar num grafo:

01

Logo, existe 0.5, 50% de chance que é esses dois casos onde a mãe é igual para todos os filhos dividido pelo total de casos, que são 4, 2/4 da 0.5

> sum(possibilidades$individuo1==possibilidades$individuo2)/nrow(possibilidades) [1] 0.5

Então veja que essa conta é a mesma da formula la, lembrando que no tempo 1, teremos aquele elemento que é elevado pelo tempo-1 como zero, e todo número elevado a zero da um, ai multiplicando o 1/n da 1/n

É claro que conforme a gente vai aumentando a população

> mae<-letters[1:3] > mae [1] "a" "b" "c"

as possibilidades aumentam

> possibilidades<-expand.grid(individuo1=mae,individuo2=mae) > possibilidades individuo1 individuo2 1 a a 2 b a 3 c a 4 a b 5 b b 6 c b 7 a c 8 b c 9 c c

mas a continha continua a mesma

> possibilidades$individuo1==possibilidades$individuo2 [1] TRUE FALSE FALSE FALSE TRUE FALSE FALSE FALSE TRUE > sum(possibilidades$individuo1==possibilidades$individuo2)/nrow(possibilidades) [1] 0.3333333

assim como nossa formula funciona bem.

> coalescent(3,1) [1] 0.3333333

E vamos aumentar mais um a população

> coalescent(4,1) [1] 0.25

tudo continua funcionando.

> mae<-letters[1:4] > possibilidades<-expand.grid(individuo1=mae,individuo2=mae) > possibilidades$individuo1==possibilidades$individuo2 [1] TRUE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE TRUE FALSE FALSE [14] FALSE FALSE TRUE > sum(possibilidades$individuo1==possibilidades$individuo2)/nrow(possibilidades) [1] 0.25

Agora vamos pensar onde entra o tempo. Se a gente olhar a formula no tempo 2, a gente vai ver o seguinte:

> coalescent(4,2) [1] 0.1875

Veja que o cenário continua o mesmo, logo nossa população é constante em 4 indivíduos.

> mae<-letters[1:4] > mae [1] "a" "b" "c" "d"

As possibilidades se mantém

> possibilidades<-expand.grid(individuo1=mae,individuo2=mae) > possibilidades individuo1 individuo2 1 a a 2 b a 3 c a 4 d a 5 a b 6 b b 7 c b 8 d b 9 a c 10 b c 11 c c 12 d c 13 a d 14 b d 15 c d 16 d d

Mas para você coalescer a 2 gerações atrás, precisamos de duas coisas, primeiro você não pode coalescer na geração anterior

> não_coalecer<-sum(!possibilidades$individuo1==possibilidades$individuo2)/nrow(possibilidades) > não_coalecer [1] 0.75

E depois você tem que coalescer na próxima geração, a 2 gerações atrás

> coalecer<-sum(possibilidades$individuo1==possibilidades$individuo2)/(nrow(possibilidades)) > coalecer [1] 0.25

So que as duas coisas tem que acontecer, logo é a chance de uma coisa acontecer e outra também, então multiplicamos as probabilidades

> não_coalecer*coalecer [1] 0.1875

e assim temos aquele número, veja que essa chance sempre vai diminuir por causa disso.

Então se a gente olhar isso ao longo do tempo, numa figura, agora que a conta ta mais clara, temos o seguinte.

02

E é isso ai, é mais provável que a gente coalesça logo, mas essa chance vai caindo. E mais uma coisa, se a gente somar essas probabilidades.

03

Que que ta acontecendo? Hora, uma hora a população tem que ter começado, então em algum momento a gente tem que coalescer, a gente tem que ter tido uma mãe em comum.

> sum(coalescent(4,1:10)) [1] 0.9436865 > sum(coalescent(4,1:20)) [1] 0.9968288 > sum(coalescent(4,1:40)) [1] 0.9999899

Veja que quanto mais para trás no tempo a gente anda, mais a chance da coalescência já ter ocorrido para quaisquer dois indivíduos. E bem, acho que é assim que a gente começa a estudar esse negocio de teoria da coalescência.

Claro que podemos expandir esse resultado para gerações não discretas como estamos olhando aqui, além de muitas outras coisas. Por exemplo, algo interessante é que vimos que depende do tempo e do tamanho da população, então podemos estimar ao longo do tempo, se uma população oscilou de tamanho, ao comparar quanto tempo dois indivíduos deveriam coalescer, com sua divergência dado um relógio molecular, se esses valores divergem, a gente pode pensar em um N variável ao longo do tempo. Mas isso vem bem depois.

Bem é isso ai, só um gostinho da teoria da coalescência, e essa formula que é jogada na nossa cara, e não sei se sou só eu que sou fraquinho para demorar a entender isso, mas 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:

Jotun Hein, Mikkel H. Schierup , Carsten Wiuf 2005 Gene Genealogies, Variation and Evolution: A Primer in Coalescent Theory 290pp

Alexei Drummond, Remco Bouckaert 2015 Bayesian Evolutionary Analysis with Beast 260pp

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
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
###
library(igraph)
 
###Definindo a função de probabilidade de coalescencia
coalescent<- function(n,t) {
    return((1/n) * (1-(1/4))^(t-1))
}
 
###Testando o modelo
coalescent(1,1)
 
coalescent(2,1)
 
###Avaliando o resultado
mae<-letters[1:2]
mae
possibilidades<-expand.grid(individuo1=mae,individuo2=mae)
 
possibilidades$individuo1==possibilidades$individuo2
 
sum(possibilidades$individuo1==possibilidades$individuo2)/nrow(possibilidades)
 
###Representação grafica das possibilidades
grafo1<-graph_from_literal(F1-a,F2-a,b)
grafo2<-graph_from_literal(F1-b,F2-a)
grafo3<-graph_from_literal(F1-a,F2-b)
grafo4<-graph_from_literal(F1-b,F2-b,a)
 
local<-matrix(c(0,0,1,0,0,1,1,1),ncol = 2,nrow=4,dimnames = list(c("F1","F2","a","b"),c("x","y")),byrow=T)
 
jpeg("01.jpg")
par(mfrow=c(2,2))
local<-local[match(rownames(local),V(grafo1)$name),]
plot(grafo1,layout=local,vertex.size=50,vertex.color="lightblue",
     vertex.label.color="black",vertex.label.cex=1.8,edge.width=4,edge.color="black")
 
local<-local[match(rownames(local),V(grafo2)$name),]
plot(grafo2,layout=local,vertex.size=50,vertex.color="lightblue",
     vertex.label.color="black",vertex.label.cex=1.8,edge.width=4,edge.color="black")
 
local<-local[match(rownames(local),V(grafo3)$name),]
plot(grafo3,layout=local,vertex.size=50,vertex.color="lightblue",
     vertex.label.color="black",vertex.label.cex=1.8,edge.width=4,edge.color="black")
 
local<-local[match(rownames(local),V(grafo4)$name),]
plot(grafo4,layout=local,vertex.size=50,vertex.color="lightblue",
     vertex.label.color="black",vertex.label.cex=1.8,edge.width=4,edge.color="black")
dev.off()
 
###Mais testes
coalescent(3,1)
 
mae<-letters[1:3]
mae
possibilidades<-expand.grid(individuo1=mae,individuo2=mae)
possibilidades
possibilidades$individuo1==possibilidades$individuo2
sum(possibilidades$individuo1==possibilidades$individuo2)/nrow(possibilidades)
 
 
coalescent(4,1)
mae<-letters[1:4]
possibilidades<-expand.grid(individuo1=mae,individuo2=mae)
possibilidades$individuo1==possibilidades$individuo2
sum(possibilidades$individuo1==possibilidades$individuo2)/nrow(possibilidades)
 
###Alterando o tempo
 
coalescent(4,2)
 
mae<-letters[1:4]
mae
 
possibilidades<-expand.grid(individuo1=mae,individuo2=mae)
possibilidades
 
não_coalecer<-sum(!possibilidades$individuo1==possibilidades$individuo2)/nrow(possibilidades)
coalecer<-sum(possibilidades$individuo1==possibilidades$individuo2)/(nrow(possibilidades))
 
não_coalecer
coalecer
 
não_coalecer*coalecer
 
###Olhando ao longo do tempo
jpeg("02.jpg")
plot(coalescent(4,1:20),xlab="Geração para trás",ylab="Probabilidade",bty="n")
dev.off()
 
jpeg("03.jpg")
plot(cumsum(coalescent(4,1:20)),xlab="Geração para trás",ylab="Somatório das Probabilidades",bty="n")
dev.off()
 
###Somas ao longo do tempo
sum(coalescent(4,1:10))
sum(coalescent(4,1:20))
sum(coalescent(4,1:40))

Evolução – Fisherian Optimality Models – Parte 5

Descendo a quantidade de um post por mês quase, mas ainda tenho esperança de reverter esse quadro.

Agora, voltando a evolução básica de Fisher, lembrando que é básica para ele, pro Fisher, não que eu ache isso básico, de forma nenhuma, gostaria, mas não dou conta de achar isso básico ainda. Mas ok, até agora nos assumimos que a medida apropriada para o fitness é o sucesso reprodutivo ao longo da vida, R0. Apesar que essa medida pode ser apropriada para uma população estável, uma medida de fitness mais geral é o parâmetro maltusiano, r, que é igual a taxa de crescimento da população na distribuição de idades estável.

 \int_0^\infty e^{-rt} l(t)m(t)dt =1

 \sum\limits_{t=1}^\infty  e^{-rt}l_t m_t = 1

onde l(t), l_t são as probabilidades de sobrevivência na idade t e m(t), m_t são os nascimentos específicos da idade de fêmeas (= fecundidade/2, assumindo uma razão sexual de 1:1). A diferença da notação usada nas duas equações são geralmente de pouco ou nenhuma consequência ( a equação da diferença pode ser igualmente escrita como a equação da integral) e usada simplesmente para ilustrar como a diferença na notação não deve implicar em diferença na interpretação. Note que a equação da diferença (do somatório) começa no fim do primeiro período de tempo, desde que nos assumimos que a fecundidade é zero no nascimento (claro que poderíamos iniciar do zero se nos simplesmente colocarmos m_0=0). Apesar que essas equações não englobam diretamente a contribuição dos machos, nos poderíamos escrever de forma similar equações para os machos relacionando seu sucesso em encontrar fêmeas a taxa de crescimento. Esse pressuposto referente ao r é que qualquer mutação que aumenta o r devera aumentar em frequência na população. A algum tempo atras, escrevemos sobre isso aqui, esse post foi muito legal, para ver como o r é determinante ao longo do tempo, para ver como um fenótipo ou característica muda de proporção nas populações, mesmo quando por exemplo ta todo mundo decaindo. Outra coisa para ter claro na cabeça olhando esse modelo, que sempre estamos levando em conta o tamanho, é lembrar da seleção sexual. Da para dar uma olhada no artigo do Ronald Aylmer Fisher aqui para ter uma ideia de como é um cara que escreveu seu nome na historia rediz um artigo, alias esses artigos meio que escritos a maquina de datilografar, sem figuras, porque acho que o povo desenhava a mão, mas uma leitura uma vez na vida de uns paper desse não mata ninguém.

Mas ok, vamos voltar ao que interessa, nos vamos considerar os mesmos pressupostos, com exceção que agora a medida de fitness é o r.

Pressupostos gerais.

1. O organismo é iteroparo.
2. Fecundidade, F, aumenta com o tamanho do corpo, x, o qual não muda depois da maturidade
3. Sobrevivência, S, diminui com o tamanho do corpo, x.
4. Fitness, W, é uma função da fecundidade e da sobrevivência.

Pressupostos matemáticos

1. Maturidade ocorre na idade 1 depois da qual, não há mais crescimento.
2. Fecundidade aumenta linearmente com o tamanho na maturidade, resultando na fecundidade sendo uma função uniforme da idade

 F_t = a_F + b_F x

3. A taxa instantânea de mortalidade aumenta linearmente com o tamanho do corpo obtido na idade 1 e é constante por unidade de tempo. Sobre esse pressuposto, a sobrevivência na idade t é dada por

 S_t = e^{-(a_S+b_S x) t}

Note que para fazer a sobrevivência declinar em função do tamanho do corpo, dado que temos uma função exponencial, nos trocamos o a_S-b_Sx por a_S+b_Sx, outra coisa que pode ser confuso aqui, é falar que é linear e a equação ser uma exponencial, veja que mesmo que tenhamos uma chance p de morrer no tempo 1, no tempo 2, vai ser o p da primeira e o p da segunda, no tempo 3, vai ser o p do tempo 1, o p do tempo 2 e o p do tempo 3, por isso aquele e elevado a equação ali.

4. O Fitness, W, é o parâmetro maltusiano r, pegando o r para ser a medida de fitness, a função de fitness é dada pela solução da equação característica:

 \int_1^\infty e^{-rt} (a_F+b_Fx)e^{-(a_S+b_Sx)t} dt =1

onde o valor inicial da integral é 1, ja que essa é a primeira idade reprodutiva.
Os dois expoentes podem ser absorvidos em um único termo, ja que os dois são “e” elevado a algo, dando

 \int_1^\infty (a_F+b_Fx)e^{-(a_S+b_Sx+r)t}dt=1

Agora, a equação acima tem a mesma forma geral da equação da da parte 4, e a integração dela vai ser mais ou menos a mesma coisa, a gente so esta adicionando a parte do crescimento populacional exponencial.

 1= (a_F+b_Fx) [\frac{-e^{-(a_S+b_Sx+r)t}}{a_S+b_Sx+r}]_1^\infty = 0 + \frac{a_F+b_Fx}{a_S+b_Sx+r} e^{-(a_S+b_Sx+r)}

Plotando a função do fitness

Para plotar r, nossa medida de fitness, como uma função de x, nos precisamos resolver a integral anterior ou a função acima para cada valor de x, ao gosto do cliente. Podemos primeiro examinar os métodos de estimar o r sem nos mesmos resolvermos a integral, e depois usando a integral resolvida.

Usando o R para fazer a integração.

Nos primeiro consideramos a integração numérica para estimar r. A estratégia é de iterar pelo espaço de x e para cada valor e calcular r.

1. integral produz a função a ser integrada, o  (a_F+b_Fx)e^{-(a_s+b_Sx+r)t}

1
2
3
4
5
6
7
8
9
10
11
 
# Função para retornar a função a ser integrada
integral <- function(idade,x,r) {
    #Parametros
    Af <- 0
    Bf <- 4*4 
    As <- 1 
    Bs <- 0.5
    #Função
    return ((Af+Bf*x)*exp(-(As+Bs*x+r)*idade))
}

2. cal_integral chama a função do R integrate para obter a integral, passando ela para o integral que acabamos de definir

1
2
3
4
cal_integral <- function(r,x) {
    #retornamos somente o valor
    return(1-integrate(integral,1,Inf,x,r)$value)
}

3. A função encontra_raiz usa a função do R uniroot para achar o valor de r que satisfaz a equação característica para um dado x.

1
2
3
encontra_raiz <- function(x) {
    uniroot(cal_integral, interval=c(1e-7,10),x)$root
}

4. Agora é só criar uma matriz de com uma única coluna valores de x e usa a função apply do R para usar encontra_raiz em cada valor de x (linha da matrix x) para pegar o valor requisito de r, guardando no vetor r.

1
2
3
4
5
6
#Figura
#Valor de x a explorar
x <- matrix(seq(0.5,3, length=100)) 
length=100
# Calcular r dado x
r <- apply(x,1,encontra_raiz)

5.O vetor x e r são então usados para plotar a relação entre r(fitness) e x(tamanho do corpo). Rm sequencia, para um dado x, RCALC chama uniroot que chama INTEGRAL que então chama INTEGRAND

1
2
#Figura
plot(x,r,type="l",xlab="Tamanho do corpo, x",ylab="Fitness, r",las=1,lwd=4)

Usando nossa solução para a integral

Nos podemos usar a nossa solução de integral também, é mais simples no quesito de código, entretanto, como muitas funções não podem ser integradas analiticamente, ela é menos geral, além do que estamos assumindo que você consegue integrar bem as funções, eu sou péssimo com integrais, saindo daquilo que são os exemplos “triviais”.

#Fazendo a figura com nossa equação integrada

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
#Nossa função de integral
funcao <- function(r,x) {
    Af <- 0
    Bf <- 4*4 
    As <- 1 
    Bs <- 0.5
    S <- exp(-(r+As+Bs*x))*(Af+Bf*x)/(As+Bs*x+r)
    return(1-S)
}
 
# Função para achar r dado x
acha_raiz <- function(x){
    uniroot( FUNC, interval=c(1e-07,10),x)$root
}
 
#Explorando x
x <- matrix(seq(0.5,3, length=100))
r <- apply(x,1,acha_raiz)
#Figura
plot(x, r, type="l", xlab="Tamanho, x", ylab="Fitness, r",las=1,lwd=4)

E então temos nossa figura:

01

E temos ai, como o fitness varia em relação ao tamanho do corpo, podemos ver que temos aquele bom e velho ponto máximo, o melhor tamanho de corpo para se ter, agora é so iniciar a busca para saber qual valor é esse.

Encontrando o máximo usando cálculo.

Nossa medida de fitness não está mais em um lado da equação, e ela não pode ser simplificada para ser assim. Mas ainda podemos usar a diferenciação implícita ja que o r depende do x. Por “conveniência” (porque senão não anda), nos usamos logaritmos para converter a função em uma que é aditiva.

 0 = ln(a_F+b_Fx)-(a_S+b_Sx +r) - ln(a_S+b_Sx+r)

Veja que temos 3 termos separados por um sinal de menos

 0 = T1 - T2 - T3

Agora, trabalhando com cada termo separadamente temos:

T1 :  \frac{dr}{dx} = \frac{b_F}{a_F+b_Fx}
Esse é fácil, não tem r na parada

T2:  \frac{dr}{dx} = b_S + \frac{dr}{dx}
Aqui veja que temos o r e x

T3:  \frac{dr}{dx} = ( \frac{1}{a_S+b_Sx+r} )(\frac{dr}{dx}) + \frac{b_S}{a_S+b_S+r}

Juntando todos os termos, ficamos com

 \frac{dr}{dx} (1+ \frac{1}{a_S+B_Sx +r}) = \frac{b_F}{a_F+b_Fx} -b_S - \frac{b_S}{a_S+b_Sx+r}

Agora quando  \frac{dr}{dx}=0 , toda a parte da direita se iguala a zero, dado que o termo no parenteses do lado esquerdo não se iguala a zero também, o que geralmente não sera o caso. Assim podemos rearranjar a parte da direita para fazer o r uma função de x:

 r^*= \frac{b_S(a_F+b_Fx^*)}{b_F-b_Sa_F-b_Sb_Fx^*} -b_Sx^* - a_S

onde r^* é o máximo valor de r, obtido no x^*. Para achar x^* nos substituímos essa equação la na primeira com os logaritmos.

 0=ln(a_F+b_Fx^*)-(a_S+b_Sx^*+r^*)-ln(a_S+b_Sx^*+r^*)

 0=ln(a_F+b_Fx^*)-[a_S+b_Sx^*+\frac{b_S(a_F+b_Fx^*)}{b_F-b^Sa_F-b_Sb_Fx^*} -b_Sx^*-a_S]  -ln[a_S+b_Sx^*+\frac{b_S(a_F+b_Fx^*)}{b_F-b_Sa_F-b_Sb_Fx^*}-b_Sx^*-a_S]

que pode ser finalmente resolvido numericamente, mas deus me livre ter que fazer isso na ponta do lápis, acho que desde a parte 2, já não rola mais heheh. Mas o R ta ai pra isso.

Usando uniroot

A primeira abordagem é o uniroot. Note que os limites são setados ligeiramente perto dos valores requeridos. Se os limites estiverem muito longes um do outro () a função pode falhar, retornando uma mensagem de erro.

> uniroot(f=funcao,interval=c(1.2,3))$root Error in uniroot(f = funcao, interval = c(1.2, 3)) : f.upper = f(upper) é NA Além disso: Warning message: In log(As + Bs * x + r) : NaNs produzidos

Esse erro mostra a importância do plot original, no “olhômetro”, da para ver quais seriam valores razoáveis e mais próximos do “ótimo”.

Então definimos a função que chegamos usando o calculo

1
2
3
4
5
6
7
8
funcao <- function(x) {
    Af <- 0;
    Bf<-4*4;
    As<-1;
    Bs<-0.5
    r <- Bs*(Af+Bf*x)/(Bf-Bs*Af-Bs*Bf*x)-Bs*x-As # r from eqn (2.40)
    return(log(Af+Bf*x)-(As+Bs*x+r)-log(As+Bs*x+r))
}

e usamos o uniroot

> uniroot(f=funcao,interval=c(1.2,1.8))$root [1] 1.389974

Usando o nlm

Uma forma alternativa é usar o nlm, usando o valor absoluto na função, que no caso o mínimo tem de ser zero.

1
2
3
4
5
6
7
8
funcao <- function(x) {
    Af<-0
    Bf<-4*4
    As<-1
    Bs<-0.5
    r <- Bs*(Af+Bf*x)/(Bf-Bs*Af-Bs*Bf*x)-Bs*x-As
    return(abs(log(Af+Bf*x)-(As+Bs*x+r)-log(As+Bs*x+r)))
}

Veja que o que mudou para a função anterior é somente o abs de valor absoluto ali na saída.
agora usamos a função nlm

> nlm(funcao, p=1.2)$estimate [1] 1.389943 Warning messages: 1: In log(As + Bs * x + r) : NaNs produzidos 2: In nlm(funcao, p = 1.2) : NA/Inf substituido pelo máximo valor positivo

Como x não é restrito a uma região, esse método não é satisfatório e apesar de termos a resposta correta, uma mensagem de aviso é gerada.

Usando optimize

Usar optimize é bem simples no entanto, basta usarmos a mesma função do nlm. Sendo que aqui definimos o mínimo e máximo, que é fácil estimar a partir da figura do fitness.

> optimize(f=funcao, interval = c(1.2,1.8),maximum= FALSE)$minimum [1] 1.389952

O intervalo de nlm de busca é muito grande, veja também que apesar de muito próximos, encontramos valores diferentes, diferentes a partir da quinta casa depois do zero, mas valores diferentes.

Usando métodos numéricos

Nos vamos ver duas abordagens aqui, a primeira, usando a equação integral resolvida e a segunda usando métodos númericos para a integração do modelo original. Esse é o caso mais geral usando métodos númericos mas também o que consome mais tempo. Mesmo que a integral possa ser resolvida, essa é uma boa prática para checar se a solução da integral está realmente correta.

Usando a função integrada.

Aqui nos usamos a função optimize, dando a ele a função que calcula a raiz (raiz_funcao) usando uniroot para achar a raiz da função especificada em “funcao”.

1
2
3
4
5
6
7
8
9
10
11
12
13
funcao <- function(r,x) {
    Af <- 0
    Bf <- 4*4
    As <- 1 
    Bs <- 0.5
    S <- exp(-(r+As+Bs*x))*(Af+Bf*x)/(As+Bs*x+r)
    return(1-S) #lembrando que aqui ja é a integral
}
 
# Achando r em função de x
raiz_funcao <- function(x){
    uniroot( FUNC, interval=c(1e-07,10),x)$root
}

E assim nos temos

> optimize(f = raiz_funcao, interval = c(.5,3),maximum = TRUE)$maximum [1] 1.389946

Usando integração númerica para a função.

Lembrando que a gente ja teve uma noção do que é integração númerica aqui e aqui, nos podemos pular a parte da integral feita a mão (caso lidemos com um modelo difícil) e integrar e depois achar o valor de r em função de x e depois maximizar o fitness. Como dito anteriormente, mesmo com uma função integrável, essa pode ser uma garantia que tudo está correto.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
funcao_integrar <- function(idade,x,r) {
    Af <- 0
    Bf <- 4*4
    As <- 1
    Bs <- 0.5
    return ((Af+Bf*x)*exp(-(As+Bs*x+r)*idade))
}
 
# Função que calcular a integral
integral <- function(r,x){
    1-integrate(funcao_integrar,1,Inf,x,r)$value
}
# Função para achar r dado x
raiz_integral <- function(x) {
    uniroot( integral, interval=c(1e-07,10),x)$root
}

E agora usamos o optimize para achar o máximo da função de fitness.

> optimize(f = raiz_integral, interval = c(1.2,1.8),maximum = TRUE)$maximum [1] 1.389934

Veja que aqui, estamos começando a esbarrar com problemas de precisão númerica, que não tem muito efeito por enquanto, mas pode vir a ser um problema.

Bem é isso ai, 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. Esse livro do Derek Roff é realmente muito legal, estamos ainda no segundo capítulo para quem quiser acompanhar, e queria eu ter a oportunidade de desde a graduação, estar exposto a um conteúdo desses, mas nunca é tarde para estudar e aprender coisas novas. Nas próxima partes, os modelos vão começar a ficar realmente legais com a entrada de variação e mais variáveis a maximizar. Até a próxima espero eu que a próxima chegue mais rápido.

Referência:
Derek A. Roff 2010 – Modeling Evolution – an introduction to numerical methods. Oxford University Press 352pp.

R. A. Fisher 1915 The evolution of sexual preference Eugenics Reviews 7(3): 184–192.

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
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
#Plotando a função de fitness
 
 
# Função para retornar a função a ser integrada
integral <- function(idade,x,r) {
    #Parametros
    Af <- 0
    Bf <- 4*4 
    As <- 1 
    Bs <- 0.5
    #Função
    return ((Af+Bf*x)*exp(-(As+Bs*x+r)*idade))
}
 
# Função para integrar a equação e retornar os valores
cal_integral <- function(r,x) {
    #retornamos somente o valor
    return(1-integrate(integral,1,Inf,x,r)$value)
}
 
 
encontra_raiz <- function(x) {
    uniroot(cal_integral, interval=c(1e-7,10),x)$root
}
 
#Figura
#Valor de x a explorar
x <- matrix(seq(0.5,3, length=100)) 
length=100
# Calcular r dado x
r <- apply(x,1,encontra_raiz)
#Figura
jpeg("01.jpg")
plot(x,r,type="l",xlab="Tamanho do corpo, x",ylab="Fitness, r",las=1,lwd=4)
dev.off()
 
 
#Fazendo a figura com nossa equação integrada
 
funcao <- function(r,x) {
    Af <- 0
    Bf <- 4*4 
    As <- 1 
    Bs <- 0.5
    S <- exp(-(r+As+Bs*x))*(Af+Bf*x)/(As+Bs*x+r)
    return(1-S)
}
 
# Função para achar r dado x
acha_raiz <- function(x){
    uniroot( FUNC, interval=c(1e-07,10),x)$root
}
 
#Explorando x
x <- matrix(seq(0.5,3, length=100))
r <- apply(x,1,acha_raiz)
#Figura
plot(x, r, type="l", xlab="Tamanho, x", ylab="Fitness, r",las=1,lwd=4)
 
#####################################
## Achando o maximo com calculo
#####################################
 
##uniroot
funcao <- function(x) {
    Af <- 0;
    Bf<-4*4;
    As<-1;
    Bs<-0.5
    r <- Bs*(Af+Bf*x)/(Bf-Bs*Af-Bs*Bf*x)-Bs*x-As
    return(log(Af+Bf*x)-(As+Bs*x+r)-log(As+Bs*x+r))
}
#
uniroot(f=funcao,interval=c(1.2,1.8))$root
 
## Usando nlm
funcao <- function(x) {
    Af<-0
    Bf<-4*4
    As<-1
    Bs<-0.5
    r <- Bs*(Af+Bf*x)/(Bf-Bs*Af-Bs*Bf*x)-Bs*x-As
    return(abs(log(Af+Bf*x)-(As+Bs*x+r)-log(As+Bs*x+r)))
}
#
nlm(funcao, p=1.2)$estimate
 
## Optimize
optimize(f=funcao, interval = c(1.2,1.8),maximum= FALSE)$minimum
 
#####################################
## Usando métodos númericos
#####################################
 
funcao <- function(r,x) {
    Af <- 0
    Bf <- 4*4
    As <- 1 
    Bs <- 0.5
    S <- exp(-(r+As+Bs*x))*(Af+Bf*x)/(As+Bs*x+r)
    return(1-S) #lembrando que aqui ja é a integral
}
 
# Achando r em função de x
raiz_funcao <- function(x){
    uniroot( FUNC, interval=c(1e-07,10),x)$root
}
 
#Achando o maximo
optimize(f = raiz_funcao, interval = c(.5,3),maximum = TRUE)$maximum
 
#####################################
## Usando integração númerica
#####################################
 
funcao_integrar <- function(idade,x,r) {
    Af <- 0
    Bf <- 4*4
    As <- 1
    Bs <- 0.5
    return ((Af+Bf*x)*exp(-(As+Bs*x+r)*idade))
}
 
# Função que calcular a integral
integral <- function(r,x){
    1-integrate(funcao_integrar,1,Inf,x,r)$value
}
# Função para achar r dado x
raiz_integral <- function(x) {
    uniroot( integral, interval=c(1e-07,10),x)$root
}
 
#Achando o maximo
optimize(f = raiz_integral, interval = c(1.2,1.8),maximum = TRUE)$maximum