Determinando a distribuição espacial a partir da Distribuição do vizinho mais próximo

Quando a gente estava medindo as distâncias ao vizinho mais próximo aqui, se a gente olhar a distribuição desses, a gente vai ver algo interessante aparecer.

Primeiro, vamos transformar em funções aquelas operações.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
##Distância entre pontos
distancia <- function(x1,y1,x2,y2){
    sqrt((x2 - x1)^2 + (y2 - y1)^2)
}
 
#Calcula distancias do vizinho mais próximo
calcula_r<- function(x,y){
    r <- vector()
    nn <- vector()
    d <- vector()
    for (i in 1:100) {
        d <- distancia(x[i],y[i],x,y)
        r[i] <- min(d[-i])
    }
    return(r)
}

O que nos interessa, é calcular a média das distâncias ao vizinho mais próximo.

1
2
3
4
5
6
7
n<-1000
dist_vmp<-vector()
 
for(i in 1:n){
    r<-calcula_r(runif(100),runif(100))
    dist_vmp[i]<-mean(r)
}

Fazemos isso muitas vezes, uma mil, e podemos dar uma olhada em como fica a distribuição dessas mil médias. Veja que estamos gerando uma distribuição de pontos ao acaso no espaço, e para cada um desses conjuntos, a gente calcula a média dos vizinhos mais próximos.

01

Veja que é assim que fica a média de vizinhos mais próximos, quando os pontos estão distribuídos ao acaso no espaço. Veja que se a gente pegar outro exemplo, e vamos ilustrar os pontos aqui.

02

Em relação a nossa distribuição, veja onde fica o média de vizinhos mais próximos.

03

Bem, até aqui, nada de novo, a gente sempre está gerando os pontos ao acaso, então temos essa distribuição média. Agora vamos considerar um segundo caso.

04

Agora nos temos os pontos agregados em alguns pontos. Agregações na natureza ocorrem por muitos motivos, pode ser mais seguro ficar perto de outros pontinhos, ali pode ser o local onde os melhores recursos estão estocados, enfim várias coisas ponde acontecer. Mas isso tem algum impacto na média do vizinhos mais próximo? Bem, todo mundo vai tender a estar mais pertinho, e assim teremos uma média mais baixa. Para esse exemplo, podemos calcular a média e comparar com a distribuição dos mil casos aleatórios.

05

E o que acontece é exatamente isso, o valor de média fica bem menor, bem antes da distribuição, já que todo mundo mais pertinho que o acaso. Porque algum motivo atrai todo os pontinho para o mesmo lugar.

Agora existe outra possibilidade, que é algo como, melhor sozinho do que acompanhado, ou seja, pense num caso de competição extrema, então a gente tenta ficar o mais longe possível de outros competidores.

07

Que gera um padrão mais ou menos como o acima. Agora o que vai acontecer com a média de distância do vizinho mais próximo? Agora ela fica um valor alto, porque sempre temos valores altos.

08

O que estamos vendo aqui? Baseado na média da distancia do vizinho mais próximo, podemo ter uma noção da distribuição espacial dos pontos, sabendo se ela é aleatória ou não, e não sendo aleatória, por esse valor, temos ainda uma noção se tudo está agregado ou espaçado. Sendo que existe um gradiente entre agregado e espaço, e o aleatório está no meio desse gradiente.

09

Teoricamente, a gente nem precisa calcular a distribuição aleatória dos pontos, porque temos uma expectativa para esse valor, que é

E(r)=\frac{\sqrt{p}}{2}

Onde r é a média das distancias do vizinho mais próximo, E é de esperança, a média para distribuições e p é a densidade de pontos, quantos pontos temos por unidade de espaço.

Se a gente olhar a média da nossa simulação

> mean(dist_vmp) [1] 0.05220953

E sabendo que temos 100 pontos, para uma área de 10000 cm2, área aqui é uma unidade meio vaga ja que estamos numa simulação.

> E_r<-sqrt(100/10000)/2 > E_r [1] 0.05

Mas chegamos a um valor bem próximo de média. E como temos uma expectativa, isso facilita a tentar interpretar esse valor como uma razão, se dividirmos o valor que encontramos por esse valor esperado, temos que:

Nos casos de distribuição aleatória, ele estará próximo de 1

> r_ale/E_r [1] 1.111665

Nos casos de agregamento, ele estára abaixo de 1

> r_agre/E_r [1] 0.4468558

E no caso de distribuição uniforme, ele será maior que 1

> r_uni/E_r [1] 1.810889

Mas apesar de tudo essa técnica leva a um outro problema, que é o seguinte. Para conseguirmos ter uma ideia do que é esperado para o local, precisamos saber a densidade de pontos do local, mas isso a gente só sabe com certeza, se tivermos certeza que estamos vendo todos os pontos do local, e como já vimos varias vezes, nossas contagens são normalmente tudo menos precisa. Logo isso é um problema para ter um valor esperado para r, o que dificulta bastante o uso desse método. Mas na teoria ele funciona bem.

Bem é isso ai, brincando mais um pouquinho com as distancias do vizinho mais próximo, 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:
Michael J Crawley 2007 – R Book 2nd. Editora Wiley 1072pp .

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
##Distância entre pontos
distancia <- function(x1,y1,x2,y2){
    sqrt((x2 - x1)^2 + (y2 - y1)^2)
}
 
#Calcula distancias do vizinho mais próximo
calcula_r<- function(x,y){
    r <- vector()
    nn <- vector()
    d <- vector()
    for (i in 1:100) {
        d <- distancia(x[i],y[i],x,y)
        r[i] <- min(d[-i])
    }
    return(r)
}
 
##Simular média do vizinho mais próximo
n<-1000
dist_vmp<-vector()
 
for(i in 1:n){
    r<-calcula_r(runif(100),runif(100))
    dist_vmp[i]<-mean(r)
}
 
##Resultado da simulação
hist(dist_vmp,main="Distribuição de médias do vizinho mais próximo")
 
##Avaliando alguns casos:
##Caso aléatorio
x<-runif(100)
y<-runif(100)
 
plot(x,y)
 
r_ale<-mean(calcula_r(x,y))
hist(dist_vmp)
abline(v=r_ale,lwd=2,col="blue",lty=2)
legend("topright",lwd=2,col="blue",lty=2,legend="Aleatório",bty="n")
 
##Caso Agregado
x<-c(rnorm(50,0.25,0.05),rnorm(50,0.6,0.1))
y<-c(rnorm(50,0.25,0.1),rnorm(50,0.7,0.05))
 
plot(x,y)
 
r_agre<-mean(calcula_r(x,y))
hist(dist_vmp,xlim=c(min(dist_vmp,r_agre),max(dist_vmp)))
abline(v=r_agre,lwd=2,col="red",lty=2)
legend("top",lwd=2,col="red",lty=2,legend="Agregado",bty="n")
 
##Caso Uniforme
x<-rep(seq(0.01,0.99,length=10),10)+runif(100,-0.002,0.05)
y<-rep(seq(0.01,0.99,length=10),each=10)+runif(100,-0.002,0.05)
 
plot(x,y)
 
r_uni<-mean(calcula_r(x,y))
hist(dist_vmp,xlim=c(min(dist_vmp),max(dist_vmp,r_uni)))
abline(v=r_uni,lwd=2,col="green",lty=2)
legend("top",lwd=2,col="green",lty=2,legend="Uniforme",bty="n")
 
hist(dist_vmp,xlim=c(min(dist_vmp,r_agre),max(dist_vmp,r_uni)))
abline(v=r_ale,lwd=2,col="blue",lty=2)
abline(v=r_agre,lwd=2,col="red",lty=2)
abline(v=r_uni,lwd=2,col="green",lty=2)
legend("topright",lwd=2,col=c("blue","red","green"),lty=2,legend=c("Aleatório","Agregado","Uniforme"),bty="n")
 
##Observado para uma distribuição ao acaso
mean(dist_vmp)
 
##Valor Predito
E_r<-sqrt(100/10000)/2
E_r
 
#Razão observado por esperado
r_ale/E_r
r_agre/E_r
r_uni/E_r

Distância até o vizinho mais próximo.

A gente ja falou de point process aqui, agora vamos ver como é a distribuição do vizinho mais próximo, e o efeito de borda que podemos ter nele, isso é meio confuso, mas vamos la.

Suponha que nos temos o problema de desenhar linhas para unir os pares de vizinhos mais próximos em um dado conjunto de pontos (x,y) que estão mapeados em duas dimensões.

Então esse é o nosso conjunto:

01

Nos precisamos de três passos.

  • Calcular a distância para o vizinho mais próximo
  • Identificar o vizinho mais próximo de cada ponto
  • usar essas distâncias minimas para identificar todos os vizinhos mais pŕoximos

Para o primeiro passo, é só calcular a distância, usando o bom e velho triângulo retângulo.

triganguloretangulo

Então a gente implementa a fórmula d=\sqrt{(y_2-y_1)^2+(x_2-x_1)^2}

1
2
3
distancia <- function(x1,y1,x2,y2){
    sqrt((x2 - x1)^2 + (y2 - y1)^2)
}

Agora nos percorremos o pontos e calculamos um vetor de distâncias, que vamos chamar de d. Ai desse vetor a gente pega a menor distância, tirando o i porque não adianta a distância de um ponto para ele mesmo, e guardamos o índice também do ponto mais próximo.

1
2
3
4
5
6
7
for (i in 1:100) {
    for (j in 1:100){
        d[j] <- distancia(x[i],y[i],x[j],y[j])
    }
    r[i] <- min(d[-i])
    nn[i] <- which(d==r[i])
}

Agora podemos colocar no gráfico, uma linha unindo os pontos mais próximos.

02

Veja que alguns pontos tem a linha de mais próximo sobrepostos. Agora vamos identificar quem está mais próximo da borda. A borda inferior e da esquerda são zero, as bordas superiores e da direta são 1, a gente pode calcular a distância para todos os pontos, vendo o menor entre todas as bordas.

borda <- pmin(x,y,dist_borda_superior,dist_borda_direita) Veja que isso funciona, porque o valor mínimo é zero, ai temos a distancia para as bordas superiores e direta, e pegamos o valor mínimo desses conjuntos, como temos um quadrado, não tem porque pegar uma distância diagonal. Agora podemos ver quantos desses valores são menores que a distância para o vizinho mais próximo.

> sum(borda

Podemos identificar esses caras no mapa

03

Agora olha o que acontece com esses valores, quando vemos a média deles com e sem esses indivíduos da borda.

04

Comumente, uma amostra é parte de um todo, então existe mais pontos, fora dessas bordas que não vemos, se pensarmos assim, veja que alguns pontos estão longe de outros pontos meramente porque não podemos ver o ponto mais próximo dele, que está fora do nosso campo de visão, fora da borda. Então eles terem recebido uma distância maior tem do vizinho mais próximo tem haver com isso na verdade. Veja que se identificamos esses caras perto demais da borda, podemos desconfiar deles, e se retiramos eles, a gente diminui a média de distâncias do vizinho mais próximo. Mas perdemos um pouco de pontos amostrados fazendo isso.

Agora pensando numa amostragem real, esse efeito pode ser muito maior se a gente usar um retângulos como área de amostra, o famoso transecto. Mas pode ser minimizado usando um circulo, que é mais difícil de implementar no campo. Então isso é algo a se levar em conta num delineamento, quando o espaço é importante, a distribuição espacial, a localização das amostras.

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.

Referência:
Michael J Crawley 2007 – R Book 2nd. Editora Wiley 1072pp .

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
set.seed(2)
 
##Distribuição de pontos aleatórios
x <- runif(100)
y <- runif(100)
 
##Num mapa
plot(x,y,pch=19,col="black",main="Mapa")
 
##Distância entre pontos
distancia <- function(x1,y1,x2,y2){
    sqrt((x2 - x1)^2 + (y2 - y1)^2)
}
 
##Menor Distância para cada ponto
 
##Vetor para guardar a menor distância
r <- vector()
##Vetor para guardar o indice do ponto mais próximo
nn <- vector()
#Vetor auxiliar para guardar todas as distâncias
d <- vector()
 
##Para cada ponto i
for (i in 1:100) {
    ##Calcular a distancia dele até o ponto j
    for (j in 1:100){
        d[j] <- distancia(x[i],y[i],x[j],y[j])
    }
    ##Qual a menor distância? Note que tiramos o i, porque a
    ##distância de um ponto pra ele mesmo não interessa, é 0
    ##Anotamos a distância
    r[i] <- min(d[-i])
    ##E o indice do ponto correspondente
    nn[i] <- which(d==r[i])
}
 
##Desenhando uma linha até o ponto mais próximo
plot(x,y,pch=19,col="black")
for (i in 1:100){
    lines(c(x[i],x[nn[i]]),c(y[i],y[nn[i]]),col="red",lty=3,lwd=2)
}
 
##Descobrindo quem está mais próximo da borda do que do vizinho mais próximo:
 
##bordas
dist_borda_superior <- 1-y
dist_borda_direita <- 1-x
 
##Distância para a borda
borda <- pmin(x,y,dist_borda_superior,dist_borda_direita)
 
##Quantos pontos estão mais próximos da borda
##Do que do vizinho mais próximo
sum(borda<r)
 
##
plot(x,y,pch=19)
indice <- which(borda<r)
points(x[indice],y[indice],col="red",cex=1.5,lwd=2)
 
##
layout(matrix(1:4,ncol=2,nrow=2))
hist(borda,xlim=c(0,0.5),ylim=c(0,20),breaks=seq(0,0.5,length.out = 20),main="Distância para borda")
hist(borda[-indice],xlim=c(0,0.5),ylim=c(0,20),breaks=seq(0,0.5,length.out = 20),main="Sem efeito de borda")
hist(r,xlim=c(0,0.15),ylim=c(0,20),breaks=seq(0,0.15,length.out = 20),main=paste("Vizinho Mais Próximo\n Média:",round(mean(r),4)))
hist(r[-indice],xlim=c(0,0.15),ylim=c(0,20),breaks=seq(0,0.15,length.out = 20),
     main=paste("Sem efeito de borda\n Média:",round(mean(r[-indice]),4)))

Existem infinitos números primos!

Essa é a pergunta que foi respondida pelo Euclides, cerca de 350 antes de Cristo, e hoje, até hoje a gente (pelo menos eu) apanha para entender.

Vamos olhar

Teorema: Existem infinitos números primos.
Suponha que exista uma quantidade finita de números primos. Seja a lista desses números p_1,p_2,p_3,\dots,p_n a lista de todos os números primos e um m um número tal que m = p_1 p_2 p_3 \dots p_n + 1. Note que m não é divisivel por p_1 ja que dividir m por p_1 vai dar um quociente de p_2 p_3 \dots p_n e um resto de 1. Similarmente, m não é divisivel por p_2,p_3,\dots,p_n.

Agora considerando o fato que todo número maior que 1 ou é um número primo, ou pode ser escrito como um produto de números primos. É claro que m é maior que 1, já que ele é alguma coisa mais 1, então m é primo ou deve ser um produto de primos.

Suponha que m seja primo, note que m é maior que todos os números na lista p_1,p_2,p_3,\dots,p_n, então nos encontramos um número que é primo que não está na lista, contradizendo nossa suposição original de que tinha a lista de todos os números primos.

Agora suponha que por outro lado m não é primo, então ele tem que ser um produto de primos. Suponha que q é um dos primos desse produto. Então m tem que ser divisível por q. Mas a gente já viu que m não é divisível por nenhum número na lista p_1 , p_2,\dots , p_n, então denovo nos temos uma contradição com a suposição de que temos uma lista de todos os números primos.
Já que a suposição de que existem uma quantidade finita de números primos leva a uma contradição, tanto considerando m como primo ou não primo, que são as duas únicas possibilidades, então deve existir uma quantidade infinita de números primos.

Se a nossa lista de primos for 2 e 3, então 2 vezes 3 da 6, mais 1 da 7 que também é primo. Se a lista fosse 2,3 e 5, multiplicando da 30, mais 1 da 31, que também e primo, ou seja, a gente pode ir gerando primos assim. A gente não gera o próximo primo da lista, mas gera um número primo, ou seja, a sempre tem mais um número primo.

Bem é isso ai, esse negócio de prova é muito divertido, 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:
Daniel J. Velleman 2006 – How to Prove It: A Structured Approach 2ed Cambridge University Press 384pp.

Biogeografia de ilhas

MacArthur e Wilson propuseram a teoria de biogeografia de ilhas onde o número de espécies de uma ilha oceânica seria uma função da taxa de imigração e extinção de espécies. O número de espécies em qualquer tempo é um equilíbrio dinâmico, resultando de ambos os lentos processos de extinção e continua chegada de novas espécies. Assim, é esperado que a composição de espécies pode mudar ao longo do tempo, gerando um “turnover” de espécies.

Vamos pensar da taxa de imigração, y, como uma função do número de espécies, x, que ja está em uma ilha. Essa relação tem uma inclinação negativa, porque conforme o número de espécies aumenta, a chance de um novo individuo chegar ser de uma nova espécie fica mais baixa, lembre-se de como é a distribuição das abundancias das espécies aqui. Assim, a imigração vai ter seu maior valor quando a ilha está vazia, quando o x é igual a zero, e vai cair para zero quando todas as espécies possíveis tiverem chegado la. Além disso, a inclinação deve desacelerar porque algumas espécies são muito mais prováveis que outras de imigrar. Isso significa que a taxa de imigração decai abruptamente assim que as espécies mais prováveis de imigrar chegam, e depois só aquelas que tem pouca chance de dispersar, ou baixas abundancias estão faltando chegar. Denovo pensando nas distribuições de abundancias de espécies, como Preston observou, algumas espécies são simplesmente mais abundantes que outras, produzindo mais indivíduos para dispersar, sendo “melhor dispersores”
Ou seja, esperamos algo desse jeito.

01

Agora vamos imaginar a taxa de extinção, y, como uma função do número de espécies, x. Essa relação vai ter uma relação positiva com o número de espécies, de tal forma que a probabilidade de extinção aumenta com o número de espécies, fácil pensar nisso, quando mais espécies, mais a gente vai adicionando a chance de extinção delas a taxa geral de extinção.

Essa relação também é predita de ter uma inclinação cada vez maior, mais ou menos pelos menos motivos do declínio cada vez maior da imigração. Algumas espécies são menos comuns que outras, e assim mais fácil de se tornarem extintas devido a estocasticidade demográfica e ambiental, tipo um drift genético, e segundo algumas espécies vão ter um menor fitness por muitas razões. Com o acumulo de espécies, maior sera o número de especies propicias a extinção (em relação a baixa abundancia e fitness) presentes, e que vão se extinguir. Lembrando que estamos falando de extinção localmente. Porque lembro que uma vez que apresentei um seminário sobre meta-populações,a palavra extinção soava agressiva para as pessoas, mas veja que extinção aqui não significa o fim da especie, só o fim local.
Então localmente, a taxa de extinção vai ser algo assim:

02

Assim, a taxa de mudança do número de especies em uma ilha vai ser \Delta R, que vai ser a diferença entre imigração I e extinção E, ou

\Delta R = I-E

Quando \Delta R=0, nos temos um equilíbrio. Se soubermos a forma quantitativa da imigração e extinção, nos podemos resolver para o equilíbrio. Esse equilíbrio vai ser um ponto no eixo x, onde as duas taxas se cruzam.

Na teoria do MacArthur e Wilson de biogeografia de ilhas, essas taxas podem ser determinadas pelo tamanho das ilhas, onde:

  • Ilhas maiores tem menor taxa de extinção por causa do maior tamanho médio das populações
  • Ilhas maiores tem maior taxa de colonização porque são alvos maiores para a dispersão

A distância entre as ilhas e as fontes de propágulos também influência essas taxas da seguinte forma:

  • Ilhas próximas do continente tem maiores taxas de colonização de novas espécies porque mais propágulos podem chegar, eles não morrem na jornada por exemplo
  • O efeito da área seria mais importante em ilhar longe do continente do que ilhas perto do continente, para reverter o continuo processo de extinção

Derivar uma curva de imigração ou extinção é um bocado complicado. Mas para ilustrar a teoria, a gente pode assumir que a taxa de imigração I, pode ser representada como uma função exponencial negativa e^{I_0-iR}, onde I_0 é a taxa de imigração em uma ilha vazia e i é o efeito negativo por espécies, que multiplica R que é quantas espécies tem la.

Então a gente pode criar um exemplo, com o número de espécies na ilha:

1
exemplo<-data.frame(especies=seq(0.1,50,0.2))

E ver a taxa de acordo com o número de espécies:

1
2
3
4
##Imigração
IO<-log(1)
b<-0.1
exemplo$imigracao<-exp(IO-b*exemplo$especies)

Agora para a extinção, E, está precisa ser zero se não houver espécies presentes. Imaginando que a extinção é uma função da densidade e que a densidade média declina conforme o número de espécies aumenta, então \bar{N}=\frac{1}{R}. então a gente vai ter algo mais ou menos assim

e^{dR}-1

Nos subtraimos esse 1 para ter certeiza que E=0 quando R=0.
O número de espécies, R, vai resultar de \Delta R = 0 = I-E, o ponto onde as linhas se cruzam.

I=e^{I_0-iR}
E=e^{dR}-1
\Delta R = 0 = I-E

03

O exato ponto onde essas linhas se cruzam podem ser encontrada de muitas formas, a gente tem visto isso em evolução aqui, mas vamos criar uma função no R para minimizar, nos vamos minimizar (I-E)^2, porque elevando ao quadrado a diferença nos temos uma quantidade com a conveniente propriedade que o mínimo ira ser aproximado tanto pelo lado positivo como o negativo.

1
2
3
deltaR<-function(R){
    (exp(IO-b*R) - (exp(d*R)-1))^2
}

Nos usamos essa função de um parâmetro no optimize, e especificamos que queremos saber o ótimo em um intervalo de 1 a 50.

> optimize(f=deltaR,interval = c(1,50)) $minimum [1] 16.91321 $objective [1] 2.813237e-13

A saída nos diz que o mínimo é algo em torno de 16.9, e podemos colocar isso na figura

04

Agora suponha uma ilha um pouco mais longe, do mesmo tamanho, mas mais longe do continente, veja que a taxa de imigração diminuirá, porque pensa, as espécies vão morrer antes de alcançar a ilha e isso vai resultar um outro ponto mínimo, outro ponto de equilíbrio.

05

A beleza dess teoria é que ela foca nossa atenção em processos a nível de paisagem, geralmente fora dos limites espaciais e temporais das nossas capacidades de amostragem. Ele demostra que qualquer fator que ajude a determinar taxas de imigração ou extinção, incluindo a area de uma ilha ou a sua proximidade a uma fonte de propágulos, vai alterar o equilíbrio do número de espécies. Mas temos que lembrar que a identidade das espécies pode mudar ao longo do tempo, isso é, um turnover de espécies, mas um turnover devagar, porque as espécies que tem uma grande chance de imigrar e menor chance de extinguir vão se manter no ambiente muito provavelmente.

Bem é isso ai, incrível como as teorias mais legais são extremamente simples, mas com um efeito profundo na forma como pensamos sobre as coisas, 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:
M.H.H. Stevens 2011, A Primer of Ecology with R

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
library(ggplot2)
 
exemplo<-data.frame(especies=seq(0.1,50,0.2))
##Imigração
IO<-log(1)
b<-0.1
exemplo$imigracao<-exp(IO-b*exemplo$especies)
 
ggplot(data=exemplo,aes(x=especies,y=imigracao))+geom_line()+
    xlab("Número de espécies") +
    ylab("Taxa de imigração")
ggsave("01.jpg")
##curve(exp(IO-b*x),0,50,xlab="n. de espécies",ylab="Taxa Imigração")
 
##Extinção
d<-0.01
exemplo$extincao<-exp(d*exemplo$especies)-1
ggplot(data=exemplo,aes(x=especies,y=extincao))+geom_line()+
    xlab("Número de espécies") +
    ylab("Taxa de extinção")
ggsave("02.jpg")
##curve(exp(d*x)-1,0,50,xlab="n. de espécies",ylab="Taxa Extinção")
 
 
##Os dois juntos
exemplo2<-stack(exemplo[,c("imigracao","extincao")])
colnames(exemplo2)<-c("valores","tipo")
exemplo2$especies<-rep(exemplo$especies,2)
 
ggplot(data=exemplo2,aes(x=especies,y=valores,group=tipo,colour=tipo))+geom_line()+
    scale_colour_discrete(name="Taxas",labels=c("Extinção","Imigração"))+
    xlab("Número de espécies") +
    ylab("Taxa")
ggsave("03.jpg")
##
deltaR<-function(R){
    (exp(IO-b*R) - (exp(d*R)-1))^2
}
estavel<-optimize(f=deltaR,interval = c(1,50))
segmento<-data.frame(x = c(estavel$minimum,estavel$minimum),y=c(0,max(apply(exemplo[,2:3],1,min))))
 
ggplot(data=exemplo2,aes(x=especies,y=valores,group=tipo,colour=tipo))+geom_line()+
    scale_colour_discrete(name="Taxas",labels=c("Extinção","Imigração"))+
    geom_line(aes(x = x, y = y,group="Estável", colour = "Estável"),data = segmento,colour="black")+
    geom_point(aes(x = x, y = y,group="Estável", colour = "Estável"), data = segmento[2,],colour="black")+
    xlab("Número de espécies") +
    ylab("Taxa")
ggsave("04.jpg")
 
exemplo<-data.frame(especies=seq(0.1,50,0.2))
exemplo$isp1<-exp(IO-0.1*exemplo$especies)
exemplo$isp2<-exp(IO-0.14*exemplo$especies)
exemplo$extincao<-exp(d*exemplo$especies)-1
exemplo2<-stack(exemplo[,c("isp1","isp2","extincao")])
colnames(exemplo2)<-c("valores","tipo")
exemplo2$especies<-rep(exemplo$especies,3)
 
##
deltaR<-function(R,b){
    (exp(IO-b*R) - (exp(d*R)-1))^2
}
b<-0.1
estavel<-optimize(f=deltaR,interval = c(1,50),b=b)
segmento1<-data.frame(x = c(estavel$minimum,estavel$minimum),y=c(0,max(apply(exemplo[,c(2,4)],1,min))))
b<-0.14
estavel<-optimize(f=deltaR,interval = c(1,50),b=b)
segmento2<-data.frame(x = c(estavel$minimum,estavel$minimum),y=c(0,max(apply(exemplo[,c(3,4)],1,min))))
 
ggplot(data=exemplo2,aes(x=especies,y=valores,group=tipo,colour=tipo))+geom_line()+
    scale_colour_discrete(name="Taxas",labels=c("Extinção","Imigração Sp1","Imigração Sp2"))+
    geom_line(aes(x = x, y = y,group="Estável", colour = "Estável"),data = segmento1,colour="black")+
    geom_point(aes(x = x, y = y,group="Estável", colour = "Estável"), data = segmento1[2,],colour="black")+
    geom_line(aes(x = x, y = y,group="Estável", colour = "Estável"),data = segmento2,colour="black")+
    geom_point(aes(x = x, y = y,group="Estável", colour = "Estável"), data = segmento2[2,],colour="black")+
    xlab("Número de espécies") +
    ylab("Taxa")
ggsave("05.jpg")

Rosalind – Introduction to Random Strings

Nesse problema, a gente tem que simular strings baseado na quantidade de GC dela. Uma das relevâncias do conteúdo de GC, se eu não estiver falando muita bobeira, é que G e C se ligam fazendo 3 pontes de hidrogênio, o que confere estabilidade a Molécula de DNA, mas tem muitos motivos diversos para medir a quantidade de GC.

Ok, mas sem devagar do problema, nossa temos uma sequência de DNA e um vetor de valores entre 0 e 1, que são proporções do conteúdo de GC de uma molécula de DNA. E temos que calcular qual a probabilidade de, ao acaso, a sequência acima ter sido gerada a partir desse conteúdo de GC.
Isso é muito parecido com sortear caras e coroas a partir de um likelihood (aqui, um post de 2012, esse é velho hehe).

Nossa entrada é:

ACGATACAA 0.129 0.287 0.423 0.476 0.641 0.742 0.783

Como são quatro bases, temos 4 possibilidades para cada posição, então temos 4\cdot4\cdot4\cdot4\cdot4\cdot4\cdot4\cdot4\cdot4=4^9=262144 possibilidades de moléculas de DNA. Agora suponto que a gente atribua uma probabilidade p de 0.129 para sair C ou G e 1-p, 0.871 para A,T, qual a chance de sair uma molécula com a mesma quantidade de GC que estamos vendo.

Primeiro vamos colocar esses dados no python

1
2
3
dna='ACGATACAA'
prob_gc='0.129 0.287 0.423 0.476 0.641 0.742 0.783'       
prob_gc = map(float, prob_gc.split())

A função map é para transformar string da leitura em um tipo que nos interessa, nesse caso float. Dai contamos a quantidade de GC

1
2
3
4
5
6
7
conteudo_gc = 0
conteudo_at = 0
for codon in dna:
	if codon in ['C', 'G']:
		conteudo_gc += 1
        else:
                conteudo_at += 1

Agora, pensando em uma quantidade apenas, para cada base, ela pode ser GC ou não, então é a chance da base ser GC ou não, coincidir com o que esperamos e a próxima também, e a próxima também. E em probabilidade, um ‘e’ significa que multiplicamos probabilidades.

1
2
3
4
5
6
7
gc=0.129
total=1.0
for base in dna:
        if base in ["G","C"]:
                total*=  gc*0.5
        else:
                total*=  (1-gc)*0.5

Mas veja que como falamos de uma base no conjunto ATCG, estamos falando da metade da probabilidade, dae multiplicamos por meio, porque a outra metade não ia resolver, certo? Mas então somamos tudo isso e é só tirar o log10

1
print log10(total)

Mas ai entramos em um problema de computação, que é o underflow, veja que a gente ta multiplicando e multiplicando um número, e ele ta ficando pequeno e pequeno até não ser possível mais representá-lo, e é ai que o log começa a valer a pena.

Veja que logs tem a propriedade de

log(x\cdot y) = log(x)+log(y)

E somando, a gente não vai encolhendo os números, então resolve o problema de precisão. Facilita bem pelo menos. Então podemos modificar o código acima para

1
2
3
4
5
6
7
8
total=0.0
for base in dna:
        if base in ["G","C"]:
                total+=  log10(gc*0.5)
        else:
                total+=  log10((1-gc)*0.5)
 
print total

Agora é só fazer isso, conteúdo de GC por conteúdo de GC

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
from math import log10
 
with open('/home/augusto/Downloads/rosalind_prob.txt') as input_data:
    dna, prob_gc = input_data.readlines()
 
prob_gc = map(float, prob_gc.split())
 
#Conta o número de GC e AT
conteudo_gc = 0
conteudo_at = 0
for codon in dna:
	if codon in ['C', 'G']:
		conteudo_gc += 1
        else:
                conteudo_at += 1
 
prob_lista  = []
for prob in prob_gc:        
	prob_final = conteudo_gc*log10(0.5*prob) + (len(dna)-conteudo_gc) * log10(0.5*(1-prob))
	prob_lista.append(str(prob_final))
 
print ' '.join(prob_lista)

Outra coisa, é que como estamos falando de quantidade de GC, não interessa a ordem das bases, e por último, temos a intuição, que o maior likelihhod vai ser no conteúdo de GC mais próximo do valor encontrado na sequência mesmo. E o join no string final é só para não imprimir a resposta desorganizada, como quando a gente imprime uma lista no python.

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.

Referência:
Rosalind