Levando em conta a correlação filogenética em uma análise de dados.

A cada ano mais e mais filogenias são publicadas, esses dados estão cada vez mais fáceis de serem obtidos, mas ainda talvez não tão utilizados como deveriam.
No R existe o pacote ape, que traz várias ferramentas para trabalhar com filogenia.
Mas uma em particular eu achei bem legal. E vamos tentar usar ela em um exemplo.

Bem pensar num exemplo, vamos comparar o tamanho de 10 espécies, 5 de uma ilha, 5 de outra ilha e ver se existe diferença entre o tamanho delas.
Então fomos em uma ilha, medimos o tamanho médio de 5 passarinhos e depois fomos a outra ilha e medimos o tamanho de 5 outras espécies de passarinho.

Com os dados em mão vamos tentar demonstrar que o tamanho dessas espécies são diferentes entre as ilhas.
Primeiro fazemos um gráfico dos dados:

01

Hummm, tudo parece bem diferente, ai fazemos o famoso teste T entre dois grupos e lindo, as ilhas são diferentes, é estrelinha de significativo para todo lado.

> summary(lm(tam~pop)) Call: lm(formula = tam ~ pop) Residuals: Min 1Q Median 3Q Max -0.82754 -0.18734 0.07906 0.29716 0.55995 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 2.3190 0.1994 11.631 2.72e-06 *** popB 1.5201 0.2820 5.391 0.000653 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.4458 on 8 degrees of freedom Multiple R-squared: 0.7842, Adjusted R-squared: 0.7572 F-statistic: 29.06 on 1 and 8 DF, p-value: 0.0006529

Mas nos dias de hoje facilmente temos acesso a filogenia, talvez não tão facilmente assim, mas com um pouco de perseverança sim, então vamos olhar a filogenia dessas espécies em relação ao seus tamanhos em dois exemplos de filogenias diferentes.

02

So para entender melhor essa figura, eu representei a filogenia e o tamanho das bolinhas é relativo ao tamanho da espécie. Eu também coloquei duas cores, uma cor para cada ilha.
A questão é que aqui a gente vê a diferença entre as espécies(tamanho das bolinhas) e vê que nos vários ramos existem tanto bolinhas grandes como pequenas, assim como nos vários ramos existem tanto bolinhas verdes como vermelha. Ou seja muitas oportunidades devem ter existido de colonizar essas ilhas, e apesar de aberta a todas as espécies, alguma coisa deve ter feito somente espécies grandes se sobressair em uma ilha, enquanto na outra ilha somente espécies pequenas. Ja vimos uma diferença entre os tamanhos, mas com a filogenia conseguimos pensar um pouco além, que todo mundo pode ter chegado a ilha, mas a ilha talvez tenha selecionado os tamanhos. Podemos então colocar essa filogenia na analise, e controlar o quanto acreditamos na diferenças entre espécies vem por causa do “ambiente” e o retirando o parentesco (efeito filogenético).
Então levando a filogenia em conta temos o seguinte resultado:

> compar.gee(tam~pop, phy =phylo1) Beginning Cgee S-function, @(#) geeformula.q 4.13 98/01/27 running glm to get initial regression estimate (Intercept) popB 2.319003 1.520127 Call: compar.gee(formula = tam ~ pop, phy = phylo1) Number of observations: 10 Model: Link: identity Variance to Mean Relation: gaussian QIC: 10.59292 Summary of Residuals: Min 1Q Median 3Q Max -0.7536925 -0.1102418 0.1550745 0.3753477 0.6337977 Coefficients: Estimate S.E. t Pr(T > |t|) (Intercept) 2.245158 0.2963257 7.576656 0.01027012 popB 1.515783 0.3384849 4.478140 0.03350552 Estimated Scale Parameter: 0.2059916 "Phylogenetic" df (dfP): 4.373501

Bem, chegamos a mesma conclusão (p=0.034), existe uma diferença entre o tamanho, não esperaríamos nada diferente disso. Mas agora sabemos que deve ter algo nas ilhas levando a essa diferença.
Agora vamos pensar em outro exemplo:

03

Bem agora comparando esse caso com o primeiro (salvo que minha árvore ta meia mal feita), a questão é que todas as espécies pequenas estão em um ramo, enquanto todas as espécies grandes estão em outro ramo, e essa ainda por cima é exatamente a divisão das ilhas. Bem não é nada sobrenatural, podemos pensar que 2 ancestrais colonizaram aquelas ilhas, um grande foi para uma ilha, um pequeno para outra. O problema é que se simplesmente fazemos um teste T e falamos que existe uma distribuição normal para descrever uma ilha, isso não vai dar certo, porque deve haver uma dependência dessa ancestralidade. Não deve ter como mudar muito em poucas gerações, e agora podemos pensar que essa diferença não deve ter haver com algo na ilha, mas a ancestralidade. Então existe uma forte correlção entre o tamanho das espécies, e deveríamos levar isso em conta no teste T. Podemos fazer a analise seguindo a ideia do Paradis e temos:

> compar.gee(tam~pop, phy =phylo2) Beginning Cgee S-function, @(#) geeformula.q 4.13 98/01/27 running glm to get initial regression estimate (Intercept) popB 2.319003 1.520127 Call: compar.gee(formula = tam ~ pop, phy = phylo2) Number of observations: 10 Model: Link: identity Variance to Mean Relation: gaussian QIC: 14.66343 Summary of Residuals: Min 1Q Median 3Q Max -0.87844244 -0.19787730 0.05506749 0.29373377 0.50904782 Coefficients: Estimate S.E. t Pr(T > |t|) (Intercept) 2.369908 0.321294 7.376135 0.008810888 popB 1.466297 0.509500 2.877914 0.077381531 Estimated Scale Parameter: 0.2003873 "Phylogenetic" df (dfP): 4.545455

Opa, agora não acreditamos mais que existe uma diferença significativa (p=0.077) entre as ilha, oras é só olhar a filogenia que nossa crença sobre o que esta acontecendo, quais as forças que impulsionam as espécies a serem desse jeito muda totalmente no exemplo 2.
A ancestralidade é muito forte e ela tem que ser levada em conta nos nossos intervalos de confiança. Esse esquema que o Paradis fez foi levar em conta a filogenia, como uma matriz de correlação entre as espécies junto com gee. Um tipo de correção de auto-correlação nos graus de liberdade. Salvo a amostra pequena (n=10 espécies) e a árvore que fiz que ficou meio feia no exemplo 2, a do exemplo 1 foi feita ao acaso, logo não mudou a conclusão, da para entender a idéia.

Essa abordagem pode ajudar muito na interpretação dos dados e enriquecer muito a discussão. E os ingredientes necessários são os dados e uma filogenia, nada mais. Como os dados a gente vai ter, é so procurar ou fazer uma árvore filogenética. Talvez seja difícil fazer uma árvore, mas existem muitas por ai que podem ser usadas. Existem muitas abordagens para incluir filogenia nas análises, mas essa é de fácil acesso (so ter o pacote ape e os ingredientes certo) e pode ser aplicada em qualquer caso que se usaria um glm. As contas podem ser difíceis de imaginar, mas a interpretação pode ser mais fácil que quando não levado em conta a filogenia como deve ter sido nesse exemplo.

Para entender melhor isso tudo da para ler o artigo do Paradis que explica essa analise, além de ele ser uma pessoa bastante acessível e legal na lista de filogenia do R.

Paradis, E. and Claude J. (2002) Analysis of comparative data using generalized estimating equations. Journal of theoretical Biology, 218, 175–185.

########################
#Gerando dados
library(ape)
set.seed(21)
pop<-factor(rep(c("A","B"),each=5))
tam<-c(rnorm(5,2,0.4),rnorm(5,4,0.4))
ind<-1:10
phylo1<-rtree(10,tip.label=ind)
cat("(((1:1,2:1):1,(3:1,4:1):1,5:2):1,((6:1,7:1):1,(8:1,9:1):1,10:1):4);",
file = "ex.tre", sep = "n")
phylo2 <- read.tree("ex.tre")
 
#Grafico e Teste T
boxplot(tam~pop,frame.plot=F,xlab="Locais",ylab="Tamanho")
 
summary(lm(tam~pop))
 
#Grafico e Analise com Exemplo 1
plot.phylo(phylo1,direction="downwards",use.edge.length=F,show.tip.label=F)
points(phylo1$tip.label,rep(0,10),cex=tam,pch=19,col=as.numeric(pop)+1)
title("Exemplo 1")
 
compar.gee(tam~pop, phy =phylo1)
 
#Grafico e analise com Exemplo 2
plot.phylo(phylo2,direction="downwards",use.edge.length=F,show.tip.label=F)
points(phylo2$tip.label,rep(0,10),cex=tam,pch=19,col=as.numeric(pop)+1)
title("Exemplo 2")
 
compar.gee(tam~pop, phy =phylo2)
##################################

Olhando o MCMC mais de perto…

Algumas vezes pode ser frustrante tentar usar estatística bayesiana.
Mas vamos continuar usando o Metropolis-Hasting para um modelo de regressão linear. E para tanto simulamos mais alguns dados

01

Primeiro vamos ver o seguinte, eu alterei o prior para sempre usar uma distribuição normal agora

prior <- function(param){
alfa = param[1]
beta = param[2]
sd = param[3]
alfaprior = dnorm(alfa, sd=2, log = T)
betaprior = dnorm(beta, sd=1, log = T)
sdprior = dnorm(sd, sd=5, log = T)
return(alfaprior+betaprior+sdprior)
}

E alterei um pouquinho a função que propõe um novo passo.

proposalfunction <- function(param){
return(rnorm(3,mean = param, sd= c(1,0.1,1)))
}

Aqui reside a seguinte decisão, se você quer passos grandes ou passos pequenos. Como assim? A cada novo passo, uma condição é testada, se vamos ficar no mesmo lugar ou vamos para esse novo lugar, se esse novo lugar tem um likelehood ruim, ficamos onde estamos, lembre-se que tem um if que faz essa restrição, sempre é sorteado uma chance de aceitar ou não o novo passo, mas quanto melhor posicionado estamos perto daquilo que deve ser a verdade, menor a chance de aceitar um passo para muito longo.
O ponto também são que usamos alfa e beta, os parâmetros do modelo, o anterior da cadeia, então ali nos desvios, quando estamos falando sd=número, é o quão longe podemos sair da onde estamos.

metropolis-hastings

Mas pagaremos um preço se dermos passos curtos, que é de ter que caminhar mais, meio óbvio não.

Mas vamos olhar num plano, como mudam os valores de intercepto e inclinação ao longo da cadeia. Então o modelo sempre tem um ponto que ele esta , e vamos unindo os pontos com linhas, então a linha vermelha é a união de como esta indo a cadeia. e vamos apagar os pontos mais antigos para facilitar a visualização, clareando eles com o tempo. Como é uma simulação, vamos deixar um pontinho azul que representa a “verdade”, que infelizmente nunca saberemos numa análise com dados reais, mas como aqui o intuito é entender melhor como a análise ta se comportando, temos esse privilegio de saber a verdade.

E vemos o seguinte

rw

Se alguém lembrou da figurinha do random walk não é por acaso, só que aqui a probabilidade da uma direção, e o random walk é em torno dos parâmetros reais.

Olhe como começamos de um ponto bem longe e a cadeia caminha para perto do ponto azul, que são os parâmetros reais usados na simulação, e la começa a ficar rodando. Usamos uma cadeia de 50 mil passos aqui, e fui colocando exponencialmente os passos no gráfico, cada vez mais e mais passos de uma vez.
Outra coisa é que não importa da onde começamos, a cadeia sempre vai convergir para a “verdade”. Por isso vai ser comum ao invés de usar um ponto inicial fico, começar de um ponto aleatório, já que se o ponto inicial influenciar no resultado final, algo não deve estar certo.

Mas vamos olhar como esta a distribuição desses 50 mil pontos ai nesse plano dos dois parâmetros (lembrando que a cadeia tem 3 partes, existe o desvio que não esta ai representado).

03

Fazendo curvas de nível para representar as densidades, igual curvas de altitude em mapinhas, vemos que o pico esta bem próximo da verdade.

Outra coisa que a gente pode ver é como parece que um parâmetro compensa o outro.
Veja a cadeias em relação aos parâmetros reais.

02

Veja como os valores de inclinação ficam um pouco abaixo da linha vermelha (o parâmetro real) e o intercepto faz o contrario, fica acima, mas quando um muda outro muda junto, vemos isso na figura anterior também, já que a figura parece uma bolinha achatada, mas que tem uma inclinação.

Por isso é importante olhar como ficam as cadeias no final de uma análise e não somente computar alguma estatística do tipo aceitação (acceptance), que apesar de ficar alta (alto ou baixo sempre são conceitos abstratos sem algo para comparar), quase 55%, ainda podemos ver esse pequeno desviozinho.

Realmente todo esse processo é bem difícil de fazer, e ficar mudando os algoritmos para qualquer alteraçãozinha, mas agora nos sabemos (apesar que bugs usa um esquema ligeiramente diferente, o gibbs sampler) o que o BUGS vai fazer com os dados, ele pula toda essa parte sofrida de fazer o algoritmos, e só vamos falar o likelehood e o prior e dar como queremos as cadeias (aqui só estávamos fazendo uma, mas o ideal é fazer mais de uma) e ele faz essa parte chata, e nos preocupamos em interpretar os resultados.

Bem por último, para aqueles que tem dificuldade de entender as curvas de nível, podemos montar aquele gráfico em perspectiva, onde a montanha é o empilhamento dos pontos e quanto mais alto a montanha, maior era a densidade de pontos ali.
Gráficos assim são ruins de visualizar, principalmente representar numeração em 3d é, mas fica tão legal ele girando que não quis jogar ele fora.

per

Agora podemos voltar a usar o Openbugs ou Winbugs com um pouco mais de conforto quanto ao que diabos esta acontecendo. Ou mesmo o stan 🙂

Aqui também estão o blog original onde vi a ideia de como fazer no R o Metropolis-hasting, o Theoretical Ecology, e se alguém quer ver um gibbs sampler, olhem o Darren Wilkinson’s research blog. É isso 🙂

###################################
# Metropolis-Hastings MCMC Denovo #
###################################
#dados de exemplo
set.seed(5)
 
alfa.verdadeiro <- 4
beta.verdadeiro <- 8
desv.verdadeiro <- 6
tamanho.amostra <- 30
 
# Criando valores do preditor
x <-runif(30,0,10)
# Criando valores dependentes a partir de "beta*x + alfa + N(0,sd)"
y <-alfa.verdadeiro + beta.verdadeiro*x + rnorm(
n=tamanho.amostra,mean=0,sd=desv.verdadeiro)
 
#Grafico
plot(x,y, main="Dados Simulados",xlab="Preditor",ylab="Resposta",
frame.plot=F,pch=19)
abline(lm(y~x))
 
#Likelihood
likelihood <- function(param){
alfa = param[1]
beta = param[2]
sd = param[3]
pred = alfa + beta*x
singlelikelihoods = dnorm(y, mean = pred, sd = sd, log = T)
sumll = sum(singlelikelihoods)
return(sumll)
}
 
#Prior distribution (Agora com distribuição normal para todos os priors)
prior <- function(param){
alfa = param[1]
beta = param[2]
sd = param[3]
alfaprior = dnorm(alfa, sd=2, log = T)
betaprior = dnorm(beta, sd=1, log = T)
sdprior = dnorm(sd, sd=5, log = T)
return(alfaprior+betaprior+sdprior)
}
 
posterior <- function(param){
return (likelihood(param) + prior(param))
}
 
######## Metropolis algorithm ################
#Note que aqui voce vai configurar o tamanho dos passos.
proposalfunction <- function(param){
return(rnorm(3,mean = param, sd= c(1,0.1,1)))
}
 
run_metropolis_MCMC <- function(startvalue, iterations){
chain = array(dim = c(iterations+1,3))
chain[1,] = startvalue
for (i in 1:iterations){
proposal = proposalfunction(chain[i,])
probab = exp(posterior(proposal) – posterior(chain[i,]))
if (runif(1) < probab){
chain[i+1,] = proposal
}else{
chain[i+1,] = chain[i,]
}
}
return(chain)
}
 
#vamos começar de algum lugar aleatorio agora
#mas cuidado, um start value ruim, pode não ser bom
startvalue <- runif(3,0,10)
chain <- run_metropolis_MCMC(startvalue, 50000)
burnIn <- 25000
 
#Agora a cadeia é maior, mas temos um acceptance menor
#Mas lembre-se que o acceptance vai depender do prior
#e do proposal
acceptance <- 1-mean(duplicated(chain[-(1:burnIn),]))
acceptance
 
#olhando estabilidade das cadeias
par(mfrow = c(3,1))
plot(chain[-(1:burnIn),1], type = "l",
xlab="" , main = "Cadeia dos valores de intercepto", )
abline(h=alfa.verdadeiro,lwd=2,col="red")
plot(chain[-(1:burnIn),2], type = "l",
xlab="" , main = "Cadeia dos valores da inclinação", )
abline(h=beta.verdadeiro,lwd=2,col="red")
plot(chain[-(1:burnIn),3], type = "l",
xlab="" , main = "Cadeia dos valores do desvio", )
abline(h=desv.verdadeiro,lwd=1,col="red")
 
#Grafico animado do MCMC para os parametros
#intercepto e inclinação
passos<-c(round(seq(1,223,length=200)^2),50000)
 
for(i in 1:200) {
jpeg(sprintf("RW%03d.jpg",i), width = 300, height = 300)
plot(0,0,type="n",xlab="Intercepto",ylab="Inclinação",
xlim=c(-2,12),ylim=c(0,10),main=paste("Passo",passos[i+1]))
cores<-rev(heat.colors(i, alpha =0.5))
for(j in 1:i) {
lines(chain[passos[j]:passos[j+1],1:2],type="l",col=cores[j])
}
points(alfa.verdadeiro,beta.verdadeiro,pch=19,col="blue")
dev.off()
}
 
system("convert RW*.jpg -delay 5 rw.gif")
system("rm RW*.jpg")
 
#preparando a densidade dos pontos
library(MASS)
z<-kde2d(chain[,1],chain[,2])
 
#contorno da densidade de pontos (passos)
plot(chain[,1],chain[,2],col=heat.colors(nrow(chain),alpha=0.4))
contour(z,lwd=3,add=T)
points(alfa.verdadeiro,beta.verdadeiro,pch=19,col="blue")
 
#perspectiva animada da densidade de pontos
theta<-seq(0,360,2)
for(i in 1:length(theta)) {
jpeg(sprintf("per%03d.jpg",i), width = 350, height = 350)
persp(z,phi = 15,theta=theta[i],xlab="Intercepto",
ylab="Inclinação",zlab="Densidade")
dev.off()
}
system("convert per*.jpg -delay 30 per.gif")
system("rm per*.jpg")
 
#lembre-se que trabalhar em outro diretorio apenas para os graficos
#pode evitar transtornos.

Rosalind Problema 6 – Finding a Motif in DNA (SUBS)

Aqui também é fácil, depois de suar a camisa para ler arquivos fasta, esses últimos dois problemas são simples. Agora a atividade é achar um motif em uma sequência. Bem motif basicamente é um pedacinho de DNA, e a gente quer encontrar onde ele esta em uma sequência, ocorra quantas vezes for. Isso deve útil para entender o que são os BLAST e o Primer da vida, mas vamos la
Então os ingredientes são uma sequência e um motif, e temos que responder onde esse motif começa na sequência.

Então vemos o seguinte:

ACGTACGTACGTACGT GTA

A primeira linha é a sequência e a segunda o motif. Bom o jeito mais simples no R seria usando Expressão Regular, só procurar esse padrão na sequência.

> gregexpr(motif,sequencia) [[1]] [1] 3 7 11 attr(,"match.length") [1] 3 3 3 attr(,"useBytes") [1] TRUE

E ali esta a resposta, no primeiro item da lista, o motif ocorre 3 vezes, e começa nas bases 3, 7 e 11.
Mas se quisermos continuar quebrando as bases com strsplit e fingir que expressões regulares não existem e reinventar a roda, poderíamos fazer o seguinte.

Picotar a sequência e o motif, dai o numero de comparações seria o tamanho da sequência – o tamanho do motif, já que temos que comparar ele inteiro, a última comparação são as ultimas 3 bases aqui.
Feito isso é só comparar o motif com a sequência, separamos as 3 primeiras bases da sequência, de 1:3 ACG e comparamos com o motif, depois denovo de 2:4 CGT e comparamos com o motif, e quanto as 3 bases forem iguais, registramos a comparação, então é um loop com um if dentro registrando quanto encontra o motif, 3 bases (o tamanho do motif) iguais.

E Tudo fica assim:

[sourcecode language=”r”]

sequencia<-"ACGTACGTACGTACGT"
motif<-"GTA"

seq.l<-strsplit(sequencia,"")[[1]]
mot.l<-strsplit(motif,"")[[1]]
comps<-length(seq.l)-length(mot.l)

local<-vector()

for(i in 1:comps){
if(sum(seq.l[i:(i+(length(mot.l)-1))]==mot.l)==length(mot.l)) {
local<-c(local,i)
}
}
[/sourcecode]

E temos no final:

> local [1] 3 7 11

Bem vamos tentar com python agora. Como em python é mais difícil reinventar a roda, vamos usar um pacote que vem com o python, o re que trabalha com expressões regulares. Daria para usar o sequencia.find(mofit), mas ele so da a primeira posição do motif, então vamos com o match e finditer. É so ir digitando esses comandos no google que a gente vai achando o que cada um faz, alias, essa solução veio da procura no google que caiu aqui.
Feito isso é só lembrar que a lista começa do zero, então é preciso somar 1 nas resposta, que fiz da pior forma possível, eu fiz uma função só para somar um no resultado e reutilizar a função map.
Espero um dia aprender direito python, mas por enquanto novamente a solução funcionou então vamos levando 🙂

[sourcecode language=”python”]
import re

sequencia = "ACGTACGTACGTACGT"
motif = "GTA"

resultado = re.search(motif,sequencia)

inicio = [match.start() for match in re.finditer(re.escape(motif), sequencia)]
inicio = map(int, inicio)

def soma1(x):
y=x+1
return y

inicio = map(soma1, inicio)
print inicio
[/sourcecode]

E a resposta é:

>>> [3, 7, 11]

Essencialmente, todos os modelos estão errados, mas alguns são úteis


 

“Essencialmente, todos os modelos estão errados, mas alguns são úteis”, ta ai uma frase que é repetida a exaustão, do original:

“essentially, all models are wrong, but some are useful”

É legal que o cara que falou isso caiu de para-quedas nessa área antes de fazer suas importantes contribuições.

O senhor George E. P. Box fazia faculdade de química, parou quando foi chamado para servir na segunda guerra mundial. La, como conhecia algo de química foi mandado pro laboratório onde testavam o efeito de armas químicas. Após alguns testes ele chegou no superior dele e disse que precisavam de um estatístico para entender aqueles resultados, ai o superior perguntou o que ele sabia de estatística, ele respondeu que tinha feito uma disciplina no primeiro ano da faculdade e o superior dele respondeu. “Está claro que você entende melhor de estatística do que eu, você analisa esses dados”, talvez não exatamente com essas palavras, mas ai ele começou a ler alguns livros, estudar, e hoje provavelmente o nome dele aparece em muitos livros de estatística.

Ele ainda foi aluno do Ronald Aylmer Fisher, e não satisfeito se casou com a filha dele.

Random Walk ou a caminhada aleatória

rw

Random walk, caminhada ao acaso ou talvez caminhada aleatória é um termo recorrente em biologia, mas que eu achei bem difícil entender.

Ele aparece e reaparece dentro da biologia sempre não importa para onde vamos.

A primeira vez que ouvi menção sobre ele foi ligado a teoria unificada da biodiversidade (UNTB) que diz que dentro dos níveis tróficos, as espécies são equivalentes, e que mudanças na abundância, composição e riqueza de espécies variam ao acaso, como um random walk. Se a gente vê uma espécie ali é porque ela teve sorte e não porque é melhor adaptada. Isso pode parecer bater de frente com a evolução, mas talvez tenha uma ponta de verdade já que como a evolução esta sempre em ação, realmente quem a gente esta vendo agora é porque sobreviveu a seleção natural, ou seja pelo menos adaptado ao momento de agora a espécie esta até certo ponto, então assumir alguma equivalência nessa adaptação entre espécies me parece razoável. Mas caso a idéia seja olhar melhor para a UNTB, o R tem o pacote UNTB para melhor entender como ela funciona e suas previsões.

Outro lugar onde a palavra Random Walk aparece é na genética, na teoria neutra, que alias diferente da UNTB, não tem muita discussão sobre sua ocorrência. Então temos alguns pedaços do DNA, sequências que não influenciam na nossa sobrevivência, logo não estão sobre influencia da seleção natural, logo elas ficam mudando e acumulando mudança, caminhando aleatoriamente sobre o plano das possibilidades. Proposta pelo senhor Kimura, tem desdobramentos muito importantes, como a teoria da coalescência e os relógios moleculares.

Mas então, Random Walks são somente uma coleção de passos ao acaso, sejam esses passos as mudanças das bases em uma sequência de DNA, mudanças da abundância das espécies em uma floresta ou pecinhas andando em um joguinho de ludo.

Podemos simular isso fácil no R simplesmente sorteando valores para coordenadas. Por exemplo, pra desenhar um ponto precisamos de cordenadas no eixo x e y, para o eixo X, a cada geração ou passo a gente sorteia entres valores 1,0 ou -1, a gente pode ir para frente, para traz ou ficar parado, fazemos o mesmo para o eixo Y e temos uma caminhada em duas dimensões. Mas ao invés de eixos e esses valores poderíamos colocar pares de bases, abundância das espécies ou ainda a diferenciação entre sequências e ver o que acontece.

Vamos simular aqui 6 populações, 6 caminhadas diferentes e ver o que acontece, vamos começar todas do centro do gráfico e para todas e andar mil passos.

01

Após caminhados 1000 passos, cada caminhada acabou em um lugar. Eu coloquei cada população com uma cor, é legal como a vermelha e a preta acabaram perto, mas cada uma das outras acabaram em um extremo do gráfico.

Eu vi a ideia desse script pela primeira vez no livro do Michael Crawley e sempre achei legal poder pensar como com o passar do tempo é assim que podem estar acontecendo as mudanças no nosso DNA e nas espécies das florestas que vemos aqui e ali.

Nesse script também tem o exemplo do inicio do post, o gráfico animado e como ele foi feito.
Eu não sei exatamente o ponto que gostaria de chegar nesse post, mas é legal ver a animação do random walk.

#Random Walk
 
#criando uma lista para depositar o resultado
random.walks<-list()
 
#simulação, o j são quantas populações eu quero, aqui são 6
#dai eu ponho no item j da lista uma matrix para receber o resultado
#ai no segundo loop interno, eu simulo o random walk, aqui para mil passos
#lembrando que tem 1001 linhas na matrix, pois para todo mundo a linha 1 é 50,50
#o meio do grafico
for (j in 1:6) {
random.walks[[j]]<-matrix(NA,ncol=2,nrow=1001,dimnames=list(1:1001,c("x","y")))
random.walks[[j]][1,]<-c(50,50)
for (i in 1:1000) {
xi<-sample(c(1,0,-1),1)
yi<-sample(c(1,0,-1),1)
random.walks[[j]][i+1,"x"]<-random.walks[[j]][i,"x"]+xi
random.walks[[j]][i+1,"y"]<-random.walks[[j]][i,"y"]+yi
}
}
 
#depois faço um grafico basico
plot(0:100,0:100,type="n",xlab="",ylab="")
 
#coloco as 6 linhas
for(i in 1:6){
lines(random.walks[[i]],col=i)
}
#uma bolinha no ultimo ponto para identificarmos onde acabou a caminhada
#apos os mil passos
for(i in 1:6){
points(random.walks[[i]][1001,"x"],random.walks[[i]][1001,"y"],col=i,pch=19,cex=1.5)
}
 
#Fazendo 200 figura jpg, organizadas por numero
for(i in 1:200){
jpeg(sprintf("RW%03d.jpg",i), width = 300, height = 300)
plot(0:100,0:100,type="n",xlab="",ylab="",main=paste("Passo",i*5))
lines(random.walks[[1]][1:(i*5),])
dev.off()
}
 
#depois é so converter para um gif animado usando o comando
#convert do imagemagik aqui no linux
system("convert RW*.jpg -delay 5 rw.gif")
#e apagar as 200 imagens
system("rm RW*.jpg")
#por segurança, é bom fazer isso em um diretorio separado