Poisson Point Process

A gente ja falou aqui sobre a distribuição dos pontos nos espaço, mas existe mais coisa a se pensar.

Eu comecei a ler o novo livro do Marc Kery, Applied hierarchical modeling in Ecology, eu ja gostava muito das coisas que ele escrevia desde os dois primeiros livros dele, Introduction to WinBUGS for Ecologists e Bayesian population analysis using WinBUGS, bem se o livro do Ben Bolker Ecological Models and Data in R abria os olhos para um monte de coisa, o Marc Kery detalhava bem modelos bayesianos desde de um começo simples, e me parece que com esse novo livro, vamos um pouquinho mais para a frente, ao mesmo estilo.

Bem uma coisa que comum no livro dele, que fala principalmente de modelos populacionais relacionados a abundância e ocupação, é o processo de pontos de poisson, ou Poisson Point Process. A gente ja falou da distribuição de poisson, que aliás é comum a muitas disciplinas, não so a ecologia. Mas vamos pensar um pouco mais sobre ela.

O Poisson point process (me perdoem usar o termo em inglês, mas as vezes é mais fácil ja ver o termo em inglês, tanto porque é mais fácil para procurar no google, achar na literatura, como muitas vezes eu traduzo as coisas erradas com meu inglês nórdico e depois vira uma confusão) tem a propriedade que cada ponto é estocasticamente independente de todos os outros pontos no processo. Apesar disso, ele é amplamente usado para modelar fenomemos representados por pontos, como fazemos direto ao usar o glm com distribuição de poisson, mas isso implica que ele nao descreve adequadamente fenomenos que tem interação suficientemente forte entre os pontos, o que é bem comum em biologia, e que explica porque a gente ve coisas como modelos quasi-poisson no glm e porque nos preocupamos com overdispersion.

O processo tem seu nome do matemático Siméon Denis Poisson e vem do fato que se uma coleção de pontos aleatórios no espaço formam um processo de poisson, então o número de pontos em uma região finita é uma variável aleatória com uma distribuição de poisson. O processo depende de um único parâmetro, que, dependendo do contexto, pode ser constante ou uma função localmente integrável (na maioria dos nossos casos a+bx ou uma função com várias médias, como numa anova). Esse parâmetro representa o valor esperado da densidade de pontos em um processo de Poisson, também conhecido como intensidade e normalmente representado pela letra grega lambda (\lambda).

A função massa de probabilidade é dada por:

p(n)=\frac{\lambda^n e^{-\lambda}}{n!}

Bem, a integral de p(n) da 1, como toda função massa de probabilidade (isso parece obvio hoje, mas eu demorei so me liguei disso depois de uma disciplina no coursera).

Para um \lambda grande, teremos o equivalente a uma distribuição normal.

Lembrando, a distribuição normal é:

p(x)=\frac{1}{\sqrt{2\sigma^2\pi}}e^{\frac{(x-\mu)^2}{2\sigma^2}}

Tudo isso basta dar ?rnorm e ?rpois, que você vai ver que ambas as distribuições são implementadas assim no R. Mas agora olha so que mágico.

Se você olhar uma parada chamada aproximação de stirling, que é uma formula fechada que mais ou menos da o resultado de um número fatorial para valores grandes, nos temos que

x! \approx \sqrt{2\pi x}e^{-x}x^x \ quando \ x \rightarrow \infty

Então, poisson é uma distribuição discreta e a normal continua, mas seja x=n=\lambda(1+\delta) onde \lambda \geq 10 (essa relação so vale para números grandes) e \delta \leq 1 (poisson é discreto, normal é contínuo). Dai substituindo temos…

p(x)=\frac{\lambda^{\lambda(1+\delta)} e^{-\lambda}}{\sqrt{2\pi \lambda(1+\delta)}e^{-\lambda(1+\delta)}(\lambda(1+\delta))^{\lambda(1+\delta)}}

Primeiro trazemos a potência do e do divisor para cima

p(x)=\frac{ e^{\lambda \delta}}{\sqrt{2\pi \lambda(1+\delta)}(\lambda(1+\delta))^{\lambda(1+\delta) - \lambda(1+\delta)}}

Depois descemos a potência do lambda la de cima

p(x)=\frac{ e^{\lambda \delta}}{\sqrt{2\pi \lambda(1+\delta)}(\lambda(1+\delta))}

simplificamos

p(x)=\frac{ e^{\lambda \delta}}{\sqrt{2\pi \lambda}(\lambda(1+\delta)) (1+\delta)^{1/2}}

ai tiramos o delta de dentro da raiz

p(x)=\frac{ e^{\lambda \delta} \lambda(1+\delta)^{-1/2}}{\sqrt{2\pi \lambda}(1+\delta)}

subimos ele

p(x)=\frac{ e^{\lambda \delta} \lambda^{-1/2}}{\sqrt{2\pi \lambda}}

cortamos o um mais delta

p(x)=\frac{ \lambda e^{-(\lambda \delta)/2}}{\sqrt{2\pi \lambda}}

e arrumamos as potências

E agora tem uma parte que eu não entendo perfeitamente, mas lembre-se que a média e a variância da distribuição de poisson são iguais, então \lambda=\mu=\sigma^2.
Assim, a parte debaixo já está pronta, e a parte de cima da equação, já é praticamente a pdf da distribuição normal, só não sei exatamente o que fazer agora, mas como vemos a relação dela com a distribuição normal nesse diagrama do blog do John Cook, é algo por ae, se alguém quiser comentar o que fazer daqui para frente, eu ficaria feliz 🙂

distribution_chart

Tudo bem que eu não sou muito bom em matemática, mas se a gente fizer um gráfico das pdfs, assim:

1
2
3
plot(1:30,dpois(1:30,10),type="l",lty=3,lwd=3,col="red",frame=F,xlab="x",ylab="Probabilidade")
points(1:30,dnorm(1:30,10,sqrt(10)),type="l",lty=3,lwd=3,col="blue")
legend("topright",legend=c("Poisson","Normal"),lty=3,lwd=3,col=c("red","blue"),bty="n")

A gente vai ver o seguinte:

01

Bem é isso ai, sem script hoje, mas olhe o repositório recologia de qualquer forma, e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail, que para esse post em especial, eu tenho certeza que devo ter falado algo errado!!! Então correções são bem vindas.

Referência:

John D. Cook – Diagram of distribution relationships
Marc Kéry & Andy Royle 2015 Applied hierarchical modeling in Ecology Elsevier 183pp

Rosalind – Suboptimal Local Alignment

Diferentes regiões do genoma evoluíram de maneira diferentes. Regiões codificantes, onde qualquer mudança pode ser o fim da vida, a maioria da variação são mutações simples, point mutations. Agora em regiões não codificantes (introns ou espaços entre genes), nos vemos uma situação diferente, intervalos inteiros podem ser facilmente duplicados, inseridos, deletados ou invertidos. Alias, a hora que todo mundo achava que esse era o DNA lixo, veio a utilidade dele para os testes de paternidade, que é feito contando repetições no “junk dna”.

Durante os anos de 1940 e 1950, Barbara McClintock uma botanica descobriu o efeito da transposição genética no milho e demostrou que muitos dos rearranjos do genoma são causados por regiões chamadas de transposons, ou intervalos pequenos de dna que podem mudar sua localização dentro do genoma. Descoberta que rendeu um nobel em Medicina, alias a única mulher a receber um nobel sozinha. Mas a gente so precisa começar a se preocupar com isso, quando o genoma é muitooo grande, o bom o velho 16s da mitocôndria não sofre esse mal, eu acho.

Mas então vamos ver o problema aqui. Nos vamos receber duas sequências:

>Rosalind_12 GACTCCTTTGTTTGCCTTAAATAGATACATATTTACTCTTGACTCTTTTGTTGGCCTTAAATAGATACATATTTGTGCGACTCCACGAGTGATTCGTA >Rosalind_37 ATGGACTCCTTTGTTTGCCTTAAATAGATACATATTCAACAAGTGTGCACTTAGCCTTGCCGACTCCTTTGTTTGCCTTAAATAGATACATATTTG

Que vamos chamar de s e t, e s e t tem uma substring r não exata de 32 a 40 pares. Devemos então contar qual a maior substring r tem o maior número de ocorrências.

Vamos olhar um exemplo então se eu pegar uma substring de t, a segunda sequencia, e comparar com s

Pegamos a segunda sequencia:

>Rosalind_37 ATGGACTCCTTTGTTTGCCTTAAATAGATACATATTCAACAAGTGTGCACTTAGCCTTGCCGACTCCTTTGTTTGCCTTAAATAGATACATATTTG

Ai a gente encaixa esse pedaço, na sequência 1

>Rosalind_12 GACTCCTTTGTTTGCCTTAAATAGATACATATTTACTCTTGACTCTTTTGTTGGCCTTAAATAGATACATATTTGTGCGACTCCACGAGTGATTCGTA x|||||x||||||x|||||||||||||||||| GGACTCCTTTGTTTGCCTTAAATAGATACATA

A gente vê que a substring encaixa la, mas com 3 erros, mas então essa substring teria uma contagem de 1. Alias, essa distancia se chama hamming distance.

No problema é sugerido que a gente use o programa Lalign, mas porque não tentar uma força bruta primeiro hehehe.

Primeiro, precisamos de uma função para calcular a hamming distance.

1
2
3
4
5
6
7
8
def hamming_dist(a,b,limite_erros=5):
    erros = 0
    for i in range(len(a)):
        if a[i] != b[i]:
            erros += 1
            if erros>limite_erros:
            	return erros
    return erros

Que basicamente é receber duas strings, e ir olhando quantas diferenças de caracteres tem neles. Aqui eu coloquei um limite de erros, apenas para diminuir um pouco o processamento, senão também demora demais, ter que sempre andar as strings inteiras. Dai vem a força bruta…

1
2
3
4
5
6
7
8
9
10
11
12
def subo(a, b):
	contagem_final=0
	for tamanho in range(32,41):
		for i in range(len(a)-tamanho):
			contagem=0
			for j in range(len(b)-tamanho):
				if hamming_dist(a[i:i+tamanho], b[j:j+tamanho])<=3:
					contagem+=1
					print a[i:i+tamanho]
					if contagem>contagem_final:
						contagem_final=contagem
	return contagem_final

Então a gente recebe duas sequências, ai vamos guardar o número de vezes que contaremos a substring, porque esse é o valor que precisamos, dai vamos fazer 3 laços, primeiro, para cada tamanho de substring valido, que é de 32 a quarenta, lembrando que o range é para string, então ele itera até o valor do segundo argumento -1, que funciona para strings, mas não no nosso caso, isso é importante, que se ficar 40, a gente não conta as substring de tamanho 40 ai responde errado. Bem dai a gente vai começar a olhar a substring a, zeramos a contagem para esse tamanho e vamos percorrer a substring b.

Então veja que se a gente tem

>A GACTCCTTTGTTTGCCTTAAATAGATACATATTTACTCTT >B ATGGACTCCTTTGTTTGCCTTAAATAGATACATATTCAAC

Tamanho vai começar valendo 32, o i 0 e o j 0, então a gente vai mandar para o função
do hamming distance

>A GACTCCTTTGTTTGCCTTAAATAGATACATATTTACTCTT >B ATGGACTCCTTTGTTTGCCTTAAATAGATACATATTCAAC

Terminando, a gente incrementa o j, então comparamos

>A GACTCCTTTGTTTGCCTTAAATAGATACATATTTACTCTT >B ATGGACTCCTTTGTTTGCCTTAAATAGATACATATTCAAC

Ai rodamos toda a string B, até acabar, dai andamos um base no A, até acabar e repetimos tudo para o substring de tamanho 33.

Cada match, nos aumentamos 1 na contagem, e toda horas que andamos o i, nos reiniciamos essa contagem, mas vamos sempre guardando o maior registro para retornar como contagem_final.

Isso é forma forma péssima de fazer, primeiro porque fazemos muitas comparações, segundo que eu acho que poderíamos usar a comparação da substring de 32 para a substring de 33, acho que mudando a ordem dos laços, e sendo um pouco mais inteligente, isso pouparia muito tempo, mas como diz o ditado…

itnotstupid

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:
Suboptimal Local Alignment

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
#!/usr/bin/env python
'''
 
'''
 
def hamming_dist(a,b,limite_erros=5):
    erros = 0
    for i in range(len(a)):
        if a[i] != b[i]:
            erros += 1
            if erros>limite_erros:
            	return erros
    return erros 
 
 
arquivo = open("/home/augusto/Downloads/rosalind_subo.txt")
 
sequencia=[]
i=-1
 
for linha in arquivo:
	if linha.find('>') != -1:		
		i+=1
		sequencia.append('')
	else:		
		sequencia[i]=sequencia[i]+linha.rstrip()
 
 
def subo(a, b):
	contagem_final=0
	for tamanho in range(32,41):
		for i in range(len(a)-tamanho):
			contagem=0
			for j in range(len(b)-tamanho):
				if hamming_dist(a[i:i+tamanho], b[j:j+tamanho])<=3:
					contagem+=1
					print a[i:i+tamanho]
					if contagem>contagem_final:
						contagem_final=contagem
	return contagem_final
print subo(sequencia[1], sequencia[0])
print subo(sequencia[0], sequencia[1])

Mapas aleatórios e o número de componentes

Uma forma de entender o ambiente melhor, é tentar ver se apesar de construirmos mapas aleatórios, ainda sim conseguimos enxergar algum padrão nessa aleatoriedade, como o número e tamanho de manchas formados

Um jeito de produzir um mapa aleatório, é simplesmente percorrer uma matriz, e sortear preencher cada célula com 0 ou 1 com uma certa probabilidade.

Podemos fazer uma função simples de ler da seguinte forma:

1
2
3
4
5
6
7
8
9
10
11
gera_mapa<-function(n,p) {
    matriz<-matrix(0,n,n)
    for(i in 1:nrow(matriz)){
        for(j in 1:ncol(matriz)){
            if(runif(1)<=p){
                matriz[i,j]<-1
            }
        }
    }
    return(matriz)
}

Então, pensando em construir mapas quadrados (menos número de linhas e colunas), nos geramos uma matriz de zeros de um tamanho n, então para cada coluna de cada linha (i,j) nos sorteamos um número entre 0 e 1, com distribuição uniforme, e se for menor que o nosso threshold, nos trocamos o zero por 1, senão passamos pro próximo.

Veja que podemos fazer uma imagem com essa matriz, onde quadradinhos pretos são 1 e brancos são 0.

01

Nesse caso, usamos um threshold de 0.5, mas se pegarmos um valor menor, tipo 0.2, temos menos quadradinhos pretos.

02

Agora, agora veja que isso é uma imagem dessa matriz:

> mapa [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [1,] 0 0 0 0 0 0 0 0 0 1 [2,] 0 0 1 0 0 1 0 0 0 0 [3,] 0 0 0 0 0 0 1 1 1 0 [4,] 0 0 0 0 0 0 0 0 0 0 [5,] 0 0 1 1 0 0 0 1 0 0 [6,] 0 0 0 0 1 0 0 0 0 0 [7,] 0 0 0 0 0 0 0 0 0 0 [8,] 0 1 0 0 0 0 0 0 0 0 [9,] 0 0 0 0 0 0 0 0 0 0 [10,] 0 0 0 0 0 0 0 0 1 0

Veja que da para ver a semelhança dos pixels pretos com os 1, em relação a localização. Até ai tudo bem. Agora o que seria uma mancha num mapa? Uma mancha seria a união de vários 1 adjacentes.

Mais ou menos assim, pegando um pixel qualquer, se houver algum pixel nas 8 posições em volta desse pixel escolhido, esses 2 pixels são u conjunto, são uma mancha.

Entao temos um pixel e 8 pixels em volta dele

1 2 3 4 x 5 6 7 8

Se em alguma posição de 1 a 8 tem outro 1, unimos os dois, como? Ue, só pensar em um pixel como um vértice e a união então é uma aresta. E assim fazemos nosso grafo.

03

Veja que nessa figura, eu dei nome para os vértices da seguinte forma:

(linha,coluna)

Ou seja, o número de vértices é igual ao número de 1 da matriz, agora, tendo o grafo é facil, basta a gente ver os componentes do grafo. Que no pacote igraph ja tem uma função implementada para fazer isso, mas poderíamos usar o algorítimo union-find como descrito no livro do Robert Sedgewick

Usando a função do igraph, teremos:

> components(grafo) $membership (8,2) (2,3) (5,3) (5,4) (6,5) (2,6) (3,7) (3,8) (5,8) (3,9) (10,9) (1,10) 1 2 3 3 3 4 4 4 5 4 6 7 $csize [1] 1 1 3 4 1 1 1 $no [1] 7

E veja que aqui conseguimos ver tanto o número de componentes formados, como o tamanho de cada componente em vértices, que é a quantidade de pixels do elemento, ou seja, o tamanho da mancha.

Certo agora estamos olhando isso para um mapinha que geramos ao acaso, e se olharmos para muitos em diferentes probabilidades de geração. Bem podemos fazer isso, com 10 repetições para cada probabilidade de geração e em média veremos isso:

04

Primeiro, podemos ver que conforme aumentamos o threshold de geração do mapa, temos menos componente, a não ser quando ele é muito baixo. Mas é so pensar, se for 0, não teremos nenhum componente, e se for 1 teremos exatamente 1 componente, cobrindo todo o mapa. Agora o meio termo é onde as coisas acontecem, rapidamente a gente começa a ter muitos componentes, até um pico la pelo 0.2 e 0.3, variável dependendo de quando fazemos a simulação, mas sempre para por aqui. Depois disso, a quantidade de componentes so diminui.

Agora, avaliando em relação ao mesmo valor de threshold, não é difícil imaginar que, quanto maior o threshold, mais arestas aparecem, mais os pixels acabam unidos, então temos manchas, em média, cada vez maiores. Mas isso demora um tempo para acontecer, tanto em probabilidades de 0.1 como 0.2, temos manchas pequenas, isso porque, a relação não é exatamente linear, so a partir de um momento, que começamos a ter um aumento de probabilidade resultando num aumento do tamanho das manchas.

Bem é isso ai, apesar dos mapas serem gerados aleatoriamente, ainda sim temos algumas certezas, o que permite alguma previsibilidade, mesmo 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:
Robert Sedgewick Kevin Wayne 2011 Algorithms 4th ed Addison-Wesley Professional 992 pp

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
#####################################
## Conectividade
#####################################
library(igraph)
 
##Função para gerar mapas quadrados aleatorios, n é o tamanho e p é a probabilidade de threshold
gera_mapa<-function(n,p) {
    matriz<-matrix(0,n,n)
    for(i in 1:nrow(matriz)){
        for(j in 1:ncol(matriz)){
            if(runif(1)<=p){
                matriz[i,j]<-1
            }
        }
    }
    return(matriz)
}
 
##Função para transformar um mapa em um grafo, com pixels com vértices e pontos adjacentes ligados por arestas
matriz_para_grafo<-function(matriz){
    aresta<-1
    grafo<-data.frame(origem=vector(),destino=vector())
    for(i in 1:nrow(matriz)){
        for(j in 1:ncol(matriz)){        
            if(matriz[i,j]!=0){            
                for(m in -1:1){
                    for(n in -1:1){
                        if(i+m>0 & i+m<=nrow(matriz) &
                           j+n>0 & j+n<=ncol(matriz)){
                            if(matriz[i+m,j+n]!=0){
                                grafo[aresta,]<-c(paste("(",i,",",j,")",sep=""),paste("(",i+m,",",j+n,")",sep=""))
                                aresta<-aresta+1
                            }
                        }
                    }
                }
            }
        }
    }
    indice<-!grafo[,1]==grafo[,2]
    indices<-which(mapa!=0,arr.ind = T)    
    grafo<-simplify(graph.data.frame(grafo[indice,],directed=F,vertices=paste("(",indices[,"row"],",",indices[,"col"],")",sep="")))  
    return(grafo)
}
 
##Função para fazer uma imagem da matriz como vemos ela na tela
plota_mapa<-function(mapa){
    image(t(mapa)[,ncol(mapa):1],col=0:1)
    }
 
 
#####################################
## Mapas aleatorios
#####################################
set.seed(123)
 
##Gera um mapa com threshold de 0.5
mapa<-gera_mapa(10,0.5)
mapa
#jpeg("01.jpg")
plota_mapa(mapa)
#dev.off()
 
##Threshold de 0.2
mapa<-gera_mapa(10,0.2)
#jpeg("02.jpg")
plota_mapa(mapa)
#dev.off()
 
##Transformando o mapa em grafo
grafo<-matriz_para_grafo(mapa)
 
#jpeg("03.jpg")
plot(grafo,vertex.size=20)
#dev.off()
 
##Usando a função do igraph
components(grafo)
 
##Realizando uma simulção com 10 repetições para probabilidades de 0.1 a 0.9
repeticoes<-10
probs<-seq(0.1,0.9,0.05)
tamanho<-15
 
##Matrizes para receber os resultados da simulação
n_compomentes<-matrix(0,ncol=repeticoes,nrow=length(probs),dimnames = list(paste("p=",probs,sep=""),paste("Amostra",1:10)))
maior_componente<-matrix(0,ncol=repeticoes,nrow=length(probs),dimnames = list(paste("p=",probs,sep=""),paste("Amostra",1:10)))
 
##Realizando a simulação
for(i in 1:length(probs)){
    for(j in 1:repeticoes){
        mapa<-gera_mapa(tamanho,probs[i])
        grafo<-matriz_para_grafo(mapa)
        n_compomentes[i,j]<-components(grafo)$no
        maior_componente[i,j]<-max(components(grafo)$csize)
        print(paste("Concluido i=",i," j=",j,sep=""))
    }
}
 
##Avaliando o resultado
#jpeg("04.jpg")
par(mfrow=c(2,1))
plot(probs,rowMeans(n_compomentes),type="b",xlim=c(0,1),ylab="Número médio de componentes",xlab="Threshold geração mapa",frame=F)
plot(probs,rowMeans(maior_componente),type="b",xlim=c(0,1),ylab="Tamanho do maior componente",xlab="Threshold geração mapa",frame=F)
#dev.off()

Rosalind – Data Formats

Como ninguém é uma ilha e existem mais dados do que podemos guardar em planilhas, o que alias é algo que eu acho que é bem ruim na biologia, sempre guardarmos as coisas em tabelas, e planilhas, o que faz com que, quando a base de dados começa a crescer, começamos a sofrer com organização e inconsistências, mesmo com poucas pessoas manipulando os dados, mas essa discussão fica para depois.

Esse problema é sobre como os dados são guardados, e como acessar um banco de dados.
O Genbank é um banco de dados de nucleotídeos (DNA) e possui uma interface para acesso remota para a maioria das linguagens de computador comuns, entre elas R e python.
O python tem um projeto chamado biopython, que traz um pacote com muita coisa para bioinformática, entre elas o pacote Entrez.

Aqui a gente so tem que aprender a olhar, e basicamente ler a documentação.

A gente recebe uma lista de ids

FJ817486 JX069768 JX469983

E tem que imprimir o registro em formato fasta do que tiver a menor sequência nessa lista.

Veja que se a gente digitar FJ817486 na entrade de busca do genbank, vamos cair aqui.

Dai é so contar o tamanho da sequência la em baixo, ou mesmo olhar ali na descrição que tem a quantidade de pares de base sequenciados, e ver quem ta lista tem a menor e imprimir no formato fasta assim:

>gi|408690371|gb|JX469983.1| Zea mays subsp. mays clone UT3343 G2-like transcription factor mRNA, partial cds ATGATGTATCATGCGAAGAATTTTTCTGTGCCCTTTGCTCCGCAGAGGGCACAGGATAATGAGCATGCAA GTAATATTGGAGGTATTGGTGGACCCAACATAAGCAACCCTGCTAATCCTGTAGGAAGTGGGAAACAACG GCTACGGTGGACATCGGATCTTCATAATCGCTTTGTGGATGCCATCGCCCAGCTTGGTGGACCAGACAGA GCTACACCTAAAGGGGTTCTCACTGTGATGGGTGTACCAGGGATCACAATTTATCATGTGAAGAGCCATC TGCAGAAGTATCGCCTTGCAAAGTATATACCCGACTCTCCTGCTGAAGGTTCCAAGGACGAAAAGAAAGA TTCGAGTGATTCCCTCTCGAACACGGATTCGGCACCAGGATTGCAAATCAATGAGGCACTAAAGATGCAA ATGGAGGTTCAGAAGCGACTACATGAGCAACTCGAGGTTCAAAGACAACTGCAACTAAGAATTGAAGCAC AAGGAAGATACTTGCAGATGATCATTGAGGAGCAACAAAAGCTTGGTGGATCAATTAAGGCTTCTGAGGA TCAGAAGCTTTCTGATTCACCTCCAAGCTTAGATGACTACCCAGAGAGCATGCAACCTTCTCCCAAGAAA CCAAGGATAGACGCATTATCACCAGATTCAGAGCGCGATACAACACAACCTGAATTCGAATCCCATTTGA TCGGTCCGTGGGATCACGGCATTGCATTCCCAGTGGAGGAGTTCAAAGCAGGCCCTGCTATGAGCAAGTC A

Claro que faremos isso usando um script do python, que vai ser assim:

from Bio import Entrez
from Bio import SeqIO

Primeiro importamos os pacotes que precisamos, o Entrez é para acessar o banco de dados e o SeqIO é para fazer o parse na saída do query que o Entrez faz.

arquivo = open('/home/augusto/Downloads/rosalind_frmt.txt')
 
ids=[]
for linha in arquivo:
	ids.extend(linha.rstrip().split())

Abrimos o arquivo que tem as chaves de acesso, para buscar no genbank, e colocamos ela dentro de uma lista, veja que o extend porque estamos adicionando na lista ids, se usássemos o append, teríamos listas dentro de listas.

Entrez.email = 'ribas.aca@gmail.com'
handle = Entrez.efetch(db='nucleotide', id=ids, rettype='fasta')
records = list(SeqIO.parse(handle, 'fasta'))

Então com al ista de ids, pedimos pelo Entrez, para buscar no banco de nucleotídeos, a lista de ids, e retornar tudo no formato fasta, depois damos um parse para deixar num formato de lista, melhor para trabalhar no python.

menor=10**6
indice=-1
for i in range(len(records)):
	if len(records[i].seq) < menor:
		menor=len(records[i].seq)
		indice=i

Agora achamos o que tem o menor tamanho, veja que records é uma lista, que tem um atributo .seq para cada instancia, que é uma string, e com len vemos o tamanho dela, assumimos um número muito grande para o menor, cuidado que potencia no python é com ** e nao ^ como no R, eu errei aqui :), e depois é so achar o menor, podíamos ter considerado o primeiro da lista como o menor também, mas tanto faz.

sequencia=records[indice].seq
lista=[]
while sequencia:
    lista.append(sequencia[:70])
    sequencia = sequencia[70:]

Uma vez que achamos, lembre-se que no records temos uma mega string, então temos que formatar de novo para fasta, então basta quebrar em linhas de 70 caracteres, aqui uma lista com esses pedaços

print '>'+records[indice].description
for linha in lista:
	print linha

E agora é so imprimir, detalhe que eu errei umas 3 vezes esse problema, porque esqueci de adicionar o “>” antes do nome da sequencia, os detalhes são sempre importantes, para fazer as coisas certas afinal.

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
'''
 
'''
 
from Bio import Entrez
from Bio import SeqIO
 
arquivo = open('/home/augusto/Downloads/rosalind_frmt.txt')
 
ids=[]
for linha in arquivo:
	ids.extend(linha.rstrip().split())
 
Entrez.email = 'ribas.aca@gmail.com'
handle = Entrez.efetch(db='nucleotide', id=ids, rettype='fasta')
records = list(SeqIO.parse(handle, 'fasta'))
 
menor=10**6
indice=-1
for i in range(len(records)):
	if len(records[i].seq) < menor:
		menor=len(records[i].seq)
		indice=i
 
print indice
 
sequencia=records[indice].seq
lista=[]
while sequencia:
    lista.append(sequencia[:70])
    sequencia = sequencia[70:]
 
print '>'+records[indice].description
for linha in lista:
	print linha

Duas espécies em metapopulações de Levins

Ja vimos metapopulações para uma espécie aqui, aqui e aqui, mas será que poderíamos adicionar mais uma espécie ao modelo de Levins?

Primeiro, vamos lembrar como é a metapopulação de levins.

\frac{dO}{dt}=cO(1-O)-eO

onde:

  • O – Ocupação
  • c – taxa de colonização
  • e – taxa de extinção

Então basicamente temos dois termos,

cN(1-N)

que faz a população crescer a uma taxa c, mas quanto maior a ocupação, mais propagulos são emitidos, mas veja que quanto maior a ocupação, menor o crescimento, que é a multiplicação por (1-N), conforme o N cresce, o valor tende a zero, porque cada vez teremos 1 menos alguma coisa, que da menos que um, estaremos multiplicando por uma fração menor.

eN

Cada mancha ocupada tem uma chance de extinção e, logo quanto mais manchas ocupadas, mais manchas veremos sendo extintas.

Agora vamos colocar outra espécie, ela vai ocupar as manchas vagas pela primeira espécie, ou seja, vagou, ela tem chance de entrar ocupar, mas a se a espécie um ocupar novamente, ela perde a competição.

Como são duas espécies, precisamos agora acompanhar a mudança na ocupação de cada uma, então precisamos de duas equações

E espécie 1 continua como a equação anterior
\frac{dO_1}{dt}=c_1O_1(1-O_1)-eO_1
mas a segunda fica
\frac{dO_2}{dt}=c_2O_2(1-O_1-O_2)-c_1*O_1*O_2-eO_2

Ou seja, a metapopulação 2 expande sua ocupação, diminuindo a expansão de acordo com o o quanto ja está ocupado pela ocupação da espécie 1 e dela mesmo, e tiramos as manchas que ela perde na competição, além do que ela perde com a extinção.

Agora, quem lembra do pacote desolve, vamos usar ele para ver como isso fica ao longo do tempo, então temos que transformar numa função para usar com a função ode.

1
2
3
4
5
6
7
8
9
10
11
metpop2sp <- function(t, y, p) {
    {
        o1 <-y[1] 
        o2 <-y[2]
    }
    with(as.list(p), {
        do1.dt <- c1*o1*(1-o1)-e*o1
        do2.dt <- c2*o2*(1-o1-o2)-c1*o1*o2-e*o2
        return(list(c(do1.dt, do2.dt)))
    })
}

Agora temos que determinar por quanto tempo queremos observar a ocupação e os paramtros definir os parametros de ocupação e extinção

1
2
3
4
5
6
tempo <- seq(0, 200, by =5)
c1 = 0.4
c2 = 0.5 
o1 = 0.1
o2 = 0.4 
e = 0.25

E vemos o seguinte

01

A espécie um rapidamente domina o ambiente, retirando a espécie 2 da jogada.

Mesmo quando mudamos os parametros um pouco, deixando ambas com a mesma taxa de colonização, maior a taxa de extinção, lembrando que para uma espécie, quando a extinção ficava maior que a colonização, ela tendia a extinção, sendo a diferença o quão rápido isso ia acontecer.

1
2
3
4
5
c1 = 0.1 
c2 = 0.1
o1 = 0.05 
o2 = 0.05 
e = 0.05

02

Agora se mudarmos o jogo para a espécie 1 extinguir, a espécie 2 consegue dominar, salvo que ela também não tenha a colonização menor que a extinção

1
2
3
4
5
c1 = 0.05
c2 = 0.07 
o1 = 0.05 
o2 = 0.05
e = 0.06

03

E existem também os casos onde podemos ter uma coexistencia

1
2
3
4
5
c1 = 0.1
c2 = 0.5
o1 = 0.05 
o2 = 0.05
e = 0.05

04

Agora, as possibilidades de combinações de parâmetros são grandes, tornando mais complexo definir quando uma ou outra espécie ganha, ou quando temos a cooexistencia, mas podemos ter um feeling da coisa já. além disso, nossas equações garantem a vitoria da espécie 1, mas isso não precisa ser assim.

Bem é isso ai, lembrando um pouco de equações diferencias, que não falamos a um bom tempo aqui e vendo a parte essa expansão para a metapopulação, 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:
A. A. de Oliveira e P. I. Prado – Coexistência em Metapopulações – Roteiro no EcoVirtual
Robert May & Angela McLean 2007 – Theoretical Ecology: Principles and Applications. Oxford University 272 pp

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
##
##install.packages("deSolve")
library(deSolve)
rm(list=ls())
 
##Função
metpop2sp <- function(t, y, p) {
    {
        o1 <-y[1] 
        o2 <-y[2]
    }
    with(as.list(p), {
        do1.dt <- c1*o1*(1-o1)-e*o1
        do2.dt <- c2*o2*(1-o1-o2)-c1*o1*o2-e*o2
        return(list(c(do1.dt, do2.dt)))
    })
}
 
##Tempo
tempo <- seq(0, 200, by =5)
 
##Caso 1
c1 = 0.4
c2 = 0.5 
o1 = 0.1
o2 = 0.4 
e = 0.25
 
out.ode<-ode(c(o1,o2), tempo, metpop2sp, list(c1=c1,c2=c2,e=e))
 
##jpeg("01.jpg")
matplot(out.ode[,1], out.ode[, -1], type = "l", lty = 1:2,frame=F,xlab="Tempo",ylab="N",lwd=2,col=c("red","blue"))
legend("right", c("Sp1", "Sp2"), lty = 1:2, col = c("red","blue"), bty = "n",lwd=2)
##dev.off()
 
##Caso 2
c1 = 0.1 
c2 = 0.1
o1 = 0.05 
o2 = 0.05 
e = 0.05
 
out.ode<-ode(c(o1,o2), tempo, metpop2sp, list(c1=c1,c2=c2,e=e))
##jpeg("02.jpg")
matplot(out.ode[,1], out.ode[, -1], type = "l", lty = 1:2,frame=F,xlab="Tempo",ylab="N",lwd=2,col=c("red","blue"))
legend("right", c("Sp1", "Sp2"), lty = 1:2, col = c("red","blue"), bty = "n",lwd=2)
##dev.off()
 
##Caso 3
c1 = 0.05
c2 = 0.07 
o1 = 0.05 
o2 = 0.05
e = 0.06
 
out.ode<-ode(c(o1,o2), tempo, metpop2sp, list(c1=c1,c2=c2,e=e))
##jpeg("03.jpg")
matplot(out.ode[,1], out.ode[, -1], type = "l", lty = 1:2,frame=F,xlab="Tempo",ylab="N",lwd=2,col=c("red","blue"))
legend("right", c("Sp1", "Sp2"), lty = 1:2, col = c("red","blue"), bty = "n",lwd=2)
##dev.off()
 
##Caso 4
c1 = 0.1
c2 = 0.5
o1 = 0.05 
o2 = 0.05
e = 0.05
 
 
out.ode<-ode(c(o1,o2), tempo, metpop2sp, list(c1=c1,c2=c2,e=e))
##jpeg("04.jpg")
matplot(out.ode[,1], out.ode[, -1], type = "l", lty = 1:2,frame=F,xlab="Tempo",ylab="N",lwd=2,col=c("red","blue"))
legend("right", c("Sp1", "Sp2"), lty = 1:2, col = c("red","blue"), bty = "n",lwd=2)
##dev.off()