Competição entre duas espécies.

Agora vamos ver como expandir o modelo de crescimento logístico para considerar mais de uma espécie.

Eu demorei um pouco para entender que tudo na verdade é apenas uma pecinha a mais no modelo de crescimento populacional.

Bem no início a gente começou assim:

 N_{t+1}= N_t+N_t\cdot(r)

Que quer dizer, a cada unidade de tempo, conforme o tempo passa, essa população dessa espécie cresce a uma taxa r, onde r é um parâmetro fixo no nosso modelinho de crescimento populacional, ou seja ele é uma característica da população, dessa população que estamos olhando. O  N_{t+1} é o mesmo que o tamanho N da população no tempo + 1 enquanto o  N_{t} é o tamanho N da população agora, neste momento.
Se a gente lembrar, isso vai dar um gráfico assim:

01

Um crescimento exponencial, mas ainda temos mais, se você como eu acha aquela formula complicada, podemos ainda pensar num sisteminha esquemático, aonde o que temos aqui é um parâmetro r (taxa de crescimento) da população da espécie que estamos trabalhando, vamos chamar de espécie 1, então temos um sisteminha assim:

01

A, mas no que me ajuda um gráfico que é uma bolinha? Calma, que conforme progredimos acaba ficando mais fácil pensar com bolinhas e risquinhos que com formulas.

Então continuamos nossa escalada a um modelo que reflete melhor o que acontece na natureza, e pensamos, poxa uma espécie não simplesmente cresce assim indiscriminadamente, conforme a população aumenta, essa população interage com ela mesmo, existe uma competição intra-específica, quanto mais gente, começa a falta espaço, comida e muitas coisas que podemos listar aqui dependendo da espécie, mas uma coisa é certa, as populações interagem com elas mesmas, e assim chegamos ao modelo de crescimento logístico.

Como era mesmo o modelo de crescimento logístico?

 N_{t+1}= N_t+N_t\cdot(r)\cdot (1-{N_t\cdot \alpha})

Basicamente a mesma coisa do modelo anterior, a gente só ta multiplicando um termo a mais, que é o que? A interação da espécie com ela mesma.

Lembrando que alpha é igual a 1/k.

Aquilo la é a mesma coisa disso aqui:

 N_{t+1}= N_t+N_t\cdot(r)\cdot (1-{N_t\over K})

E o gráfico disso é o seguinte:

01

Mas a gente pode olhar esse mesmo sisteminha como um grafo, como a figurinha de bolinha que fizemos ali em cima, e o que adicionamos aqui é uma interação da espécie com ela mesma. E isso da:

02

Temos agora a população da espécie 1 e uma interação dela, nela mesma, um alpha da população 1 na população 1.

Muito bem, até aqui ja exploramos um bocado por aqui. Mas o mundo não é feito de populações isoladas, o mundo é um monte de população de várias espécies emboladas numa pedra grandona (conhecida como planeta terra, pelo menos por nos, alienígenas devem ter outro nome hehe).
Então vamos colocar outra espécie na jogada.

Primeira coisa, vamos pensar quantas interações teremos?

Bem que a espécie 1 interage com ela mesma e a espécie 2 interage com ela mesma,isso é certo, a interação intra-específica, se não consideramos isso estamos diminuindo a escala de complexidade que estamos tentado subir aqui.

Mas o que mais? Bem a espécie 1 tem que interagir com a espécie 2, senão vai continuar a mesma coisa do sistema anterior, mas se a espécie 1 interage com a espécie 2, então a espécie 2 tem que interagir com a espécie 1, isso é a interação inter-específica, uma espécie interagindo com a outra.

E onde isso nos leva? A isso:

03

Olha que bonito, e mais legal, se você pensar como nas interações, quanto mais espécies colocamos, sempre teremos uma matriz quadrada. Nesse caso não temos quatro alphas?
Podemos colocar eles organizadinhos dentro de uma matriz assim:

 \left(\begin{array}{cc}  \alpha_{11} & \alpha_{12} \\  \alpha_{21} & \alpha_{22} \end{array}\right)

Viu, a diagonal da matriz é exatamente a interação da espécie com ela mesmo, ai um lado da matriz é a interação que uma espécie causa na outra, e o outro lado é a interação que ela sofre.

E como isso fica na formulinha?

 \begin{cases}  N_{1,t+1}= N_{1,t}+ N_{1,t}\cdot(r_{1})\cdot (1-{N_{1,t}\cdot \alpha_{11}} - {N_{2,t}\cdot \alpha_{21}}) \\  N_{2,t+1}= N_{2,t}+N_{2,t}\cdot(r_{2})\cdot (1-{N_{2,t}\cdot \alpha_{22}} - {N_{2,t}\cdot \alpha_{12}})  \end{cases}

Ai você pode pensar, nossa que difícil, mas não, a gente so colocou mais uma coisinha, que é essa parte aqui:

 - {N_{2,t}\cdot \alpha_{21}}

E é so pensar, poxa, se a interação da espécie com ela mesma a gente descontava aqui nesse parenteses, qualquer outra interação, a gente desconta o efeito dela nesse mesmo parenteses, e isso vale para todas as espécies que estamos levando em conta, aqui são duas, mas poderiam ser mais se você quiser, mas depois vemos como lidar com mais espécies.

Agora é facil implementar isso no R, é so a gente escrever essa formula ai no R assim:

dlvcomp2 <- function(N, alpha, rd=c(1,1)) {
  N1.t1 <- N[1] + rd[1] * N[1] * (1 - alpha[1,1]*N[1] - alpha[1,2]*N[2])
  N2.t1 <- N[2] + rd[2] * N[2] * (1 - alpha[2,1]*N[1] - alpha[2,2]*N[2])
  c(N1.t1, N2.t1)
}

Notem que alpha a gente implementa como uma matriz, exatamente a matriz que fizemos ali em cima.

So que ai vamos dar os valores que quisermos.

Por exemplo, para esses valores:

[,1] [,2] [1,] 0.010 0.005 [2,] 0.008 0.010

A população iniciando com tamanho 10 para cada espécie e os r iguais, de 1 para cada espécie, veremos o seguinte:

04

Muito legal né, agora as espécies estão interagindo, mas podemos ainda implementar esse mesmo sistema de equações de forma continua.

lvcomp2 <- function(t, n, parms) {
  with(as.list(parms), {
    dn1dt <- r1*n[1]*(1-a11*n[1] - a12*n[2])
    dn2dt <- r2*n[2]*(1-a22*n[2] - a21*n[1])
    list(c(dn1dt, dn2dt))
  } )
}

Agora aqui podemos usar o pacote deSolve para resolver essa equação, vamos olhar um exemplo com outros valores de parâmetros.

05

Bom, parece que funciona bem. Agora é so sair resolvendo essas equações para tudo que é combinação de parâmetros possível e ver o que acontece, que acabamos por pegar alguns insights interessantes.

Mas para começar a gente pode ir jogando valores para tentar repetir algumas situações reais, como as que tem nos exemplos do livro de ecologia do Begon, Townsend e Harper, esse exemplo aqui:
begonecology

Pense, como deveriam ser os alphas e os rs para isso ai acontecer? Conseguindo fazer isso é um sinal que estamos entendo bem o que ta acontecendo com essas populações.

Bem é isso, apos alguns dias sem posts novos. Mas agora podemos começar a explorar alguns conceitos de competição entre espécies também.

Bem o script esta la no repositório recologia do github e aqui em baixo. Bom fim de semana a todos.

#########################################################################
#Competição Lotka Volterra
#########################################################################
library(igraph)
library(deSolve)
 
#Figura 1
plot(graph.formula("Sp1" ),vertex.size=50,vertex.color=NA,vertex.label.color="black",vertex.label.cex=1.5,
     vertex.label=expression(r[1]))
 
sistema<-matrix(c("Sp1", "Sp1"),nc=2,byrow=TRUE)
sistema
rede.sistema<-graph.data.frame(sistema, directed=TRUE, vertices=NULL)
#Figura 2
plot(rede.sistema,edge.curved=seq(0.3, 0.3, length = ecount(rede.sistema)),
     edge.color="blue",edge.width=1,vertex.size=50,vertex.color=NA,vertex.label.color="black",vertex.label.cex=1.5,
     edge.label=c(expression(alpha[11])),edge.label.cex=2,edge.label.color="black",vertex.label=expression(r[1]))
 
 
sistema<-matrix(c("Sp1", "Sp1","Sp1", "Sp2","Sp2", "Sp2","Sp2", "Sp1"),nc=2,byrow=TRUE)
sistema
rede.sistema<-graph.data.frame(sistema, directed=T, vertices=NULL)
#Figura 3
plot(rede.sistema,layout=layout.spring,edge.curved=seq(0.3, 0.3, length = ecount(rede.sistema)),
     edge.color="blue",edge.width=1,vertex.size=50,vertex.color=NA,vertex.label.color="black",vertex.label.cex=1.5,
     edge.label=c(expression(alpha[11]),expression(alpha[12]),expression(alpha[22]),expression(alpha[21])),
     edge.label.cex=2,edge.label.color="black",vertex.label=c(expression(r[1]),expression(r[2])))
 
#Competição discreta:
dlvcomp2 <- function(N, alpha, rd=c(1,1)) {
  N1.t1 <- N[1] + rd[1] * N[1] * (1 - alpha[1,1]*N[1] - alpha[1,2]*N[2])
  N2.t1 <- N[2] + rd[2] * N[2] * (1 - alpha[2,1]*N[1] - alpha[2,2]*N[2])
  c(N1.t1, N2.t1)
}
 
#Parametros
alphs<-matrix(c(0.01,0.005,0.008,0.01),ncol=2,byrow=TRUE)
t<-20
 
N <- matrix(NA, nrow=t+1, ncol=2)
N[1,] <- c(10,10)
 
for(i in 1:t) {
  N[i+1,]<-dlvcomp2(N[i,],alphs)
}
 
 
#Gráfico
#Figura 4
matplot(0:t, N, type='l', col=1, ylim=c(0,110),frame.plot=F,xlab="Tempo",
        ylab="Tamanho da População",main="Crescimento logístico discreto para duas espécies")
abline(h=1/alphs[1,1], lty=3)
text(0, 1/alphs[1,1], "K", adj=c(0,0))
legend("right", c(expression("Espécie 1 "*(alpha[21]==0.008)),
                  expression("Espécie 2 "*(alpha[12]==0.005))),
       lty=1:2,bty='n')
dev.off()
 
#Crescimento continuo
lvcomp2 <- function(t, n, parms) {
  with(as.list(parms), {
    dn1dt <- r1*n[1]*(1-a11*n[1] - a12*n[2])
    dn2dt <- r2*n[2]*(1-a22*n[2] - a21*n[1])
    list(c(dn1dt, dn2dt))
  } )
}
 
 
#resolvendo continuamente com desolve
 
parms <- c(r1=1,r2=0.1,a11=0.2,a21=0.1,a22=0.02,a12=0.01);
initialN<-c(2,1)
out<-ode(y=initialN, times=1:100, func=lvcomp2, parms=parms)
 
#Figura 5
matplot(out[,1],out[,-1],type='l',xlab="Tempo",ylab="Tamanho populacional",
        frame.plot=F,main="Crescimento logístico continuo para duas espécies")
legend("left",c(expression("Espécie 1 "*(alpha[21]==0.1)),
                expression("Espécie 2 "*(alpha[12]==0.01))),
       lty=1:2,bty='n',col=c("black","red"))

Update:

Notem que isso…

 N_{t+1}= N_t+N_t\cdot(r)

é o mesmo que…

 N_{t+1}= N_t \cdot ( 1 + r)

É so por o N_t em evidência. Caso alguém esteja se perguntando porque essa formula pode parecer diferente da de alguns posts anteriores.

Teste t Binomial – GLM com prevalências.

A gente ja usou um bocado estatística bayesiana usando Openbugs. Os modelos feitos no Openbugs são facilmente replicados no Winbugs, Jags e com um pouquinho de dedicação no Stan. Mas a maioria dos posts anteriores era apenas usando a distribuição normal para descrever os dados.

Mas as vezes não é possível descrever o mundo assim. Um exemplo é quando a gente trabalha com algo como a prevalência de um parasita.
Vamos pegar um exemplo aqui, passarinhos sendo parasitados por Haemoproteus columbae.

Haemoproteus é um genero de protozoa que parasita passarinhos, repteis e anfíbios. O seu nome é derivado do grego, “Haima” = Sangue e
“Proteus” = Deus dos mares que tem o poder de assumir varias formas diferentes. Dai da para sacar que que ele vive no sangue e tem mais de uma forma.

vol102(3)_foto_capa2

Então temos passarinhos, e eles podem ou não estar contaminado. Mas como descobrimos?

Bem existem dois jeitos, ou a gente faz uma lâmina do sangue, prepara ela e olha no microscópio e fica procurando o Haemoproteus ou a gente prepara o sangue, faz uma pcr colocando primers que replicam o DNA do haemoproteus mas não do passarinho hospedeiro (devem vir varias células, logo vários DNAS diferentes num pouquinho de sangue) e se houver bastante DNA amplificado no final, bingo, o parasita tem que estar ali.

Mas no final das contas o que vamos falar é o seguinte, vamos ter sangue de vários hospedeiros, e nos dois casos vamos apenas falar se o indivíduo, o passarinho esta ou não parasitado pelo Haemoproteus.

Individuo rep.dados.Sp..each...15. Parasita 1 1 Columbina talpacoti Sim 2 2 Columbina talpacoti Sim 3 3 Columbina talpacoti Sim 4 4 Columbina talpacoti Sim 5 5 Columbina talpacoti Não . . . 25 10 Scardafella squammata Sim 26 11 Scardafella squammata Não 27 12 Scardafella squammata Não 28 13 Scardafella squammata Sim 29 14 Scardafella squammata Não 30 15 Scardafella squammata Não . . .

Acho que da para entender como estão os dados, temos uma identidade para cada indivíduo, cada amostra, ou seja um passarinho é uma amostra. Cada amostra pertence a uma espécie, e sabemos se ele esta ou não parasitada, temos algo como Sim ou Não, parasitado ou não parasitado, como quiser chamar.

Bem como são sim ou não, normalmente as pessoas resumem essa informação em contagens de quantos passarinhos estão parasitados e quantos não estão parasitados.

> tabela Parasita Sp Sim Não Columbina talpacoti 32 30 Scardafella squammata 10 47

Aqui mantemos a informação de quantos passarinhos foram amostrados (soma de todos as células da matriz), quantos de cada espécie, quantos parasitados ou não por espécie, notem que com essa matriz podemos reconstruir a planilha la em cima, isso se não nos importamos com a identidade de cada passarinho para outras informações, mas temos aqui a mesma informação que la em cima.

Virou uma tabelinha assim, a primeira coisa que vem a cabeça é teste de chi-quadrado para ver se temos um desvio nas proporções. Como assim? Uai se o parasita ocorre numa frequência constante, independente da espécie, então esperaríamos uma proporção igual entre as linhas e colunas, agora se algo acontece, uma espécie é mais parasitada ou, sei la, a prevalência é baixinha nas espécies, ou a gente pode observar se o desvio das contagens é devido ao acaso ou há razões para desconfiarmos que algo a mais acontece.

Então realizamos o teste de chi-quadrado:

> teste.chi <- chisq.test(tabela) > teste.chi Pearson's Chi-squared test with Yates' continuity correction data: tabela X-squared = 13.6387, df = 1, p-value = 0.0002216 > teste.chi$observed # Contagens observadas Parasita Sp Sim Não Columbina talpacoti 32 30 Scardafella squammata 10 47 > teste.chi$expected # Contagens esperada sobre a hipotese nula Parasita Sp Sim Não Columbina talpacoti 21.88235 40.11765 Scardafella squammata 20.11765 36.88235

E temos o clássico p muito menor que 0.05, e mais um monte de valor. Bem, como não fornecemos as probabilidades para o esperado, ele calcula automaticamente, mas podemos ver o que esperávamos no resultado olhando o $expected do modelo ajustado. Bem eu ja fiquei perdido nesses números. Quando a gente olha por exemplo o resultado de cruzamentos de Drosophila melanogaster em relação a algumas mutações, é fácil estabelecer mentalmente o que é o esperado, 1:1 em caso do sexo, 3:1 no caso de uma mutação dominante e assim vai, mas essa pergunta pode ser mais complicada para alguns casos. Então a pergunta, “o que você espera de resultado” é sempre uma pergunta interessante. Mas de qualquer forma aqui a gente vê que algo acontece, agora vamos olhar mais a frente.

Outra forma de olhar os dados é calcular a prevalência, que é, para uma espécie nesse exemplo, o quanto de casos de parasitados tivemos do total de dados.

> dados Sp N Prevalence simul 1 Columbina talpacoti 62 0.516 32 2 Scardafella squammata 57 0.193 10

E temos a prevalência, que é mais ou menos a chance de um indivíduo estar parasitado. E Para nossa surpresa, pelo menos pra minha, dividir o número de indivíduos parasitados pelo total de indivíduos amostrados da a máxima verosimilhança (Maximum likelihood) da distribuição Binomial, mas no próximo post vamos olhar isso.
Mas antes de continuar, vamos ver o que estamos querendo falar ao dizer que 0.516, ou 52% é a prevalência do Haemoproteus na espécie Columbina talpacoti

01

Bem, são dados de contagem aqui, mas temos que ter em mente que essa contagem tem um limite, que no primeiro caso é 62 e no segundo 57. Deixamos riscado no gráfico esse limite.
Agora para o gráfico 1, onde esta a bolinha mais alta? Oras, no 32, o maximum likelihood, mas existe um erro inerente a qualquer amostragem, mas para essa mesma prevalência, ainda seria razoável pensar nela como o valor real de prevalência, a real chance de um passarinho ser contaminado mesmo com 31 passarinhos doentes. Agora pense no valor de 30 ou 29, ja estamos ficando com metade da altura que tem o maximum likelihood, aqui talvez não seja muito legal pensar em 51.6% como prevalência, e pensar que a prevalência é um pouco mais baixa. O mesmo vale para outra especie.

Mas então temos crenças diferentes para diferentes valores de prevalência para cada uma das especies. Mas além de estabelecer quanto é a prevalência de cada espécie, podemos tentar estimar a diferença entre essas prevalência. Que seria o que? Algo como o quanto deve ser a diferença entre elas.

So que aqui tem um esquema, meramente falando em diferença, tudo fica meio confuso, a gente não esta meramente diminuindo 0.516 de 0.193. A gente quer saber a distribuição, assim como essas que estimamos para cada especie, a distribuição das diferenças, porque não vamos simplismente falar um valor, ele pode estar errado, a gente tem que ver todas as possibilidades de diferenças.

E a gente pode fazer isso como similar a forma como fizemos o teste t aqui, so que para usar a distribuição binomial, precisamos de uma função de ligação, que nesse caso é a função logistica.

Vamos ver primeiro como é a função na linguagem Bugs:

sink("Binomial.t.test.txt")
cat("
model {
 
# Priors
 alpha ~ dnorm(0,0.01)
 beta ~ dnorm(0,0.01)
 
# Likelihood
 for (i in 1:n) {
    C[i] ~ dbin(p[i], N[i])
    logit(p[i]) <- alpha + beta *species[i]
 }
# Derived quantities
 prev.C_talpacoti <- exp(alpha) / (1 + exp(alpha))
 prev.S_squammata <- exp(alpha + beta) / (1 + exp(alpha + beta))
 prev.Diff <- prev.C_talpacoti - prev.S_squammata	# Teste
}
",fill=TRUE)
sink()

A gente ja fez quase isso aqui, sem usar MCMC. Mas agora a diferenças que a gente tem, ao invés de apenas 1 parâmetro, um valor de probabilidade (que aqui a gente le prevalência), a gente tem dois, uma para uma espécie e outro para outra espécie.

Então fazemos aqui temos que fazer isso para as duas espécies.

Então o likelihood tem um loop que vai de 1:2, aqui o n é dois quando a gente preenche os dados.

# Likelihood
 for (i in 1:n) {
    "Meu modelo aqui"
 }

Dentro desse loop a gente vai por nosso modelo.
E nosso modelo é como?
Tem uma parte estocástica e uma parte determinística.

Como é a parte estocástica?

C[i] ~ dbin(p[i], N[i])

Oras, a parte estocástica é uma probabilidade, que depende do N, o número de amostras, é uma distribuição binomial. Mas esse p é determinado pela nossa função determinística, e como é ela?

logit(p[i]) <- alpha + beta *species[i]

Igual foi em todo modelo linear, a+bx, intercepto + inclinação * espécie. Vezes a espécie porque eu vou dar um valor para cada espécie. Aqui que ta o macete, eu vou atribuir a probabilidade de uma espécie como alpha, e o da outra como alpha main beta. so que beta é multiplicado por um valor, que é a hora que eu mando junto dos dados o 0 ou 1, ou seja, vamos supor um exemplo simples, eu tenho espécies A e B, A tem prevalência de 20 e B tem prevalência de 30
Nesse caso, para esse modelo, o intercepto é 20, a inclinação é 10 e o beta da especie A é 0 e o beta da especie B é 1.
Assim ficamos com o seguinte:

Para a espécie A A = interceopto + inclinação * especie A = 20 + 10 * 0 A = 20

Como qualquer número multiplicado por zero é zero, esse 10 some para a especie 1, e ela é meramente o valor do intercepto.
Ja pra a espécie B, a espécie 1, temos o seguinte:

Para a espécie B B = interceopto + inclinação * especie B = 20 + 10 * 1 B = 30

Oras, 10 vezes 1 é 10 mesmo, que somado com o 20 da 30, a prevalência da segunda espécie.
Isso é só pra facilitar a visualização do que ta dentro do likelihood ali, mas observe que o valor de p(de prevalência) esta dentro da função logit, a função logística, usado para generalizar para linearizar a relação e assim usar a distribuição binomial. Por isso o a+bx da equação vai sair com um número esquisito de olhar, que não é porcentagem, não é um número entre 0 e 1, entre 0% e 100%. Para olhar ele depois temos que destransformar, e fazemos isso usando as quantidades derivadas no Openbugs.

# Derived quantities
 prev.C_talpacoti <- exp(alpha) / (1 + exp(alpha))
 prev.S_squammata <- exp(alpha + beta) / (1 + exp(alpha + beta))
 prev.Diff <- prev.C_talpacoti - prev.S_squammata	# Teste

Agora olha essa equação ai e lembra da formula da equação logística, poxa aqui a gente so ta destransformando a continha que falamos acima, veja a formula no wikipedia, é exatamente isso a formula da função logística, so lembrar que o inverso de log é exponencial. So que além de acessar os valores de prevalência das duas espécies, olhamos também a diferença entre eles, e ja teremos uma saída direto.

Agora vamos rodar o modelo e ver o que acontece.

> out <- bugs(data=bugs.dados, inits=inits, parameters.to.save=params,model.file="Binomial.t.test.txt", + n.thin=nt, n.chains=nc, n.burnin=nb,n.iter=ni) > print(out, dig = 3) Inference for Bugs model at "Binomial.t.test.txt", Current: 3 chains, each with 1200 iterations (first 200 discarded), n.thin = 2 Cumulative: n.sims = 3000 iterations saved mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff alpha 0.065 0.253 -0.432 -0.113 0.068 0.243 0.539 1.001 3000 beta -1.636 0.424 -2.452 -1.929 -1.637 -1.342 -0.807 1.001 2400 prev.C_talpacoti 0.516 0.062 0.394 0.472 0.517 0.561 0.632 1.001 3000 prev.S_squammata 0.177 0.049 0.095 0.142 0.174 0.208 0.285 1.001 3000 prev.Diff 0.339 0.079 0.180 0.284 0.342 0.395 0.481 1.003 2800 deviance 10.488 1.855 8.609 9.127 9.914 11.250 15.331 1.001 2100 For each parameter, n.eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor (at convergence, Rhat=1). DIC info (using the rule, pD = Dbar-Dhat) pD = 1.9 and DIC = 12.4 DIC is an estimate of expected predictive error (lower deviance is better).

Primeiro olhamos o Rhat para ver a convergência dos valores, senão nem adianta olhar o resultado, como todos estão ok a gente procede para olhar os resultados, se você não lembra o que é o Rhat, pode dar uma conferida aqui.

Mas então o alpha e beta, que é o a+bx da equação da reta ta ae, mas são números que estão transformados pela função logística. Vamos olhar direto a prevalência das espécies, que ja destransformamos os dados diretamente no bugs, e la temos, dado que simulamos dados com essas mesmas características, recuperamos muito bem os valores de prevalência, como recuperamos bem, vamos olhar a diferença, e a diferença de prevalências tem média de 0.339 e desvio de 0.079. Ta mas mais legal que isso vamos olhar qual o intervalo de confiança de 95% dele, é algo que vai de 0.180 a 0.481, ou seja dentro desse intervalo não tem nenhum valor de 0, ou negativo, isso quer dizer que no minimo existe uma diferença de 0.18 e no máximo 0.481 para as prevalência, então se pegarmos por exemplo a média de C. talpacoti, que é 0.516 de prevalência, nesse caso o S. squammata tem no minimo 0.516-0.180 e no maximo 0.516-0.481.
Essa conta é esse teste aqui:

prev.Diff <- prev.C_talpacoti - prev.S_squammata	# Teste

Mas fazemos ele para cada valor do MCMC, por isso temos a cadeiazinha dele e assim uma distribuição.
Vamos olhar esses três valores das quantidades derivadas num gráfico.

02

O primeiro e o segundo, são as distribuições de qual a prevalência, e nossa confiança no valor. O risco vermelho é a media da distribuições, bem próximo das distribuições que tinhamos originalmente certo?

So que aqui também calculamos a distribuição da diferença, e vemos que ela é quase que certamente maior que zero, a gente não vê a distribuição nem encostando no zero, e a diferença ser zero é o mesmo que não haver diferença, se você é zero cm mais alto que eu, eu posso afirmar que temos a mesma altura, mas se você é no minimo três centímetros mais alto que eu, eu não posso afirmar que temos a mesma altura, simples assim.

Então nesse caso podemos afirmar que “Existe uma diferença na prevalência de Haemoproteus nessas duas espécies”, a mais ainda, podemos afirmar que “Existe uma diferença, e ela deve ser algo entre 18% e 34% menos prevalente em S. squammata do que em C. talpacoti“, agora porque isso ocorre, precisamos ler mais do que eu li para escrever esse post.

Bem, para constar, podemos fazer o mesmo teste usando estatística frequentista assim.

> modelo.glm<-glm(cbind(dados$simul, dados$N - dados$simul) ~ dados$Sp, family = binomial) > summary(modelo.glm) Call: glm(formula = cbind(dados$simul, dados$N - dados$simul) ~ dados$Sp, family = binomial) Deviance Residuals: [1] 0 0 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.06454 0.25413 0.254 0.799529 dados$SpScardafella squammata -1.61210 0.43111 -3.739 0.000184 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 1.5693e+01 on 1 degrees of freedom Residual deviance: -3.5527e-15 on 0 degrees of freedom AIC: 12.551 Number of Fisher Scoring iterations: 3

Notem que temos os mesmos valores para alpha e beta praticamente, não esperamos nada diferente, já que o prior que usamos foi pouco informativo, mas como fizemos aqui a analise de três formas diferentes, chi-quadrado, glm num framework bayesiano e glm num framework frequentistca, e tudo convergiu para a mesma resposta, existe uma diferença entre a espécies quanto a sua prevalência de Haemoproteus, devemos estar no caminho certo.

Bem acho que acabei falando um monte de coisa aqui, mas acho que agora com o tempo da para ir destrinchando o glm melhor nos próximos posts. Mas glms são muito legais, ja que da para usar, nesse mesmo esquema, muitos tipos de distribuições, cobrindo uma gama muito maior de tipos de dados que com os modelo lineares mas usando todas as estrategias que ja tinhamos para cobrir os modelos lineares, aqueles em que os dados tem uma distribuição normal.

Referencias
Adriano E.A. & Cordeiro N.S. 2001 – Prevalence and Intensity of Haemoproteus columbae in Three Species of Wild Doves from Brazil – Mem. Inst. Oswaldo Cruz 96(2) : 175-178

Marc Kéry 2010 – Introduction to WinBUGS for Ecologists.

set.seed(123)
### Gerando dados
dados<-data.frame(Sp=c("Columbina talpacoti","Scardafella squammata"),N=c(62,57),
                  Prevalence=c(0.516,0.193))
 
dados$simul<-0
 
for(i in 1:2) {
    dados[i,"simul"]<-rbinom(1, dados[i,"N"], prob = dados[i,"Prevalence"])
    }
 
layout(matrix(1:2,ncol=1,nrow=2))
for(i in 1:2) {
    plot(dbinom(x=0:dados$N[i],size=dados$N[i],prob=dados$Prevalence[i])~c(0:dados$N[i]),frame=F,pch=19,
         main=paste(dados$Sp[i]," N=",dados$N[i],sep=""),ylab="Frequência",xlim=c(0,70),xlab="",ylim=c(0,0.15))
    legend("topright",legend=paste("Prevalência =",dados$Prevalence[i]*100,"%"),bty="n")
    lines(c(dados$N[i],dados$N[i]),c(0.05,0),col="blue",lty=2,lwd=3)
    text(dados$N[i]+5,0.05,paste("Limite N =",dados$N[i]),cex=0.7)
    }
 
tabela<- as.table(cbind(dados$simul, dados$N - dados$simul))
dimnames(tabela) <- list(Sp = c("Columbina talpacoti","Scardafella squammata"),Parasita = c("Sim","Não"))
tabela
 
teste.chi <- chisq.test(tabela)
 
str(teste.chi)
teste.chi$observed   # Contagens observadas
teste.chi$expected   # Contagens esperada sobre a hipotese nula
teste.chi$residuals  # Residuo de Pearson
teste.chi$stdres     # Residos "standardized"
 
### Análise usando Bugs
# Definindo Modelo
 
sink("Binomial.t.test.txt")
cat("
model {
 
# Priors
 alpha ~ dnorm(0,0.01)
 beta ~ dnorm(0,0.01)
 
# Likelihood
 for (i in 1:n) {
    C[i] ~ dbin(p[i], N[i]) 
    logit(p[i]) <- alpha + beta *species[i]
 }
# Derived quantities
 prev.C_talpacoti <- exp(alpha) / (1 + exp(alpha))
 prev.S_squammata <- exp(alpha + beta) / (1 + exp(alpha + beta))
 prev.Diff <- prev.C_talpacoti - prev.S_squammata	# Teste
}
",fill=TRUE)
sink()
 
# Junte os dados numa lista
bugs.dados <- list(C = dados$simul, N = dados$N, species = c(0,1), n = length(dados$Sp))
 
# Função para iniciar os dados
inits <- function(){ list(alpha=rlnorm(1), beta=rlnorm(1))}
 
# Parametros para acompanhar
params <- c("alpha", "beta", "prev.C_talpacoti", "prev.S_squammata", "prev.Diff")
 
# Configurações do MCMC
nc <- 3
ni <- 1200
nb <- 200
nt <- 2
 
# Iniciando o Gibbs sampling
library(R2OpenBUGS)
out <- bugs(data=bugs.dados, inits=inits, parameters.to.save=params,model.file="Binomial.t.test.txt",
            n.thin=nt, n.chains=nc, n.burnin=nb,n.iter=ni)
 
print(out, dig = 3)
 
out$summary
 
#Figura 2
str(out)
par(mfrow = c(3,1))
hist(out$sims.list$prev.C_talpacoti, col = "grey", xlim = c(0,1), xlab="",
	main = "Prevalência de C. talpacoti", breaks = 30,ylim=c(0,300))
abline(v = out$mean$prev.C_talpacoti, lwd = 3, col = "red")
hist(out$sims.list$prev.S_squammata, col = "grey", xlim = c(0,1), xlab= "",
	main = "Prevalência de S. squammata", breaks = 30,ylim=c(0,300))
abline(v = out$mean$prev.S_squammata, lwd = 3, col = "red")
hist(out$sims.list$prev.Diff, col = "grey", xlim = c(0,1), xlab = "",
	main = "Diferença nas Prevalências", breaks = 30,ylim=c(0,300))
abline(v = 0, lwd = 3, col = "red")
 
### Analise usando GLM
modelo.glm<-glm(cbind(dados$simul, dados$N - dados$simul) ~ dados$Sp, family = binomial)
summary(modelo.glm)

Arduino – Pisca luzes binomial

DSC04419

E mais um pequeno projetinho bem bobo com Arduino.

Domingo, sem nada para fazer, eu tentei ver como faz para ler dados do Arduino no computador.
Na linguagem dele existe um comando Serial e o modelo Uno tem uma porta serial para transmitir ou receber dados do computador, que da para fazer via o cabo usb.

Bem eu fiz um pequeno programinha, ainda de piscar luzes, mas agora eu fiz elas piscarem de forma aleatória. Usando o gerador do números aleatórios que vem no c, o random(), eu gerava um número, ai testava se ele era par ou impar, se for par eu acendo uma luz azul, senão eu acendo a luz amarela no caso de impar. Ele vai ficar no fim do post.

DSC04423

Até ai não tem muita coisa, mas o legal é que agora eu criei duas variáveis, para manter uma contagem de quantas vezes cada luz ficava acesa, depois eu imprimia no monitor do ide do Arduino.

Captura de tela de 2013-06-09 13:53:29

Assim eu conseguia receber um dado direto do Arduino.

Falta aprender melhor como fazer um esquema para o computador ir lendo os dados e salvando num arquivo, isso seria bem útil, e parece que não será uma tarefa muito árdua ja que muita gente quer fazer isso.

Em python, eu vi aqui que existe uma library chamada pyserial exatamente para fazer isso, so que não especifica para Arduino. Na verdade essa é uma pergunta bem recorrente como vemos nos posts aqui e aqui por exemplo, para iniciar.

Bom é isso, bom domingo.

//Apenas para usar um define para aprender, agora toda pausa vale 700
//So mudar o valor aqui para pausas diferentes
#define pausa 700L
 
//Iniciando as variáveis que serão usadas
int carapin = 13 , coroapin = 12 ;
int cara = 0 , coroa = 0 , n = 100 , numero , i;
 
void setup() {
//Aqui a gente inicia tantos os pinos como a porta serial
//para comunicar com o computador
  Serial.begin(9600);
  pinMode(carapin, OUTPUT);
  pinMode(coroapin, OUTPUT);
  randomSeed(analogRead(0));
}
 
// Aqui é o loop que vai ficar em modo infinito.
void loop() {
 
  //Vamos fazer contagens de 100 vezes apenas, para ver o final
  for(i=1;i<=n;i++) {
    //Aqui geramos um número pseudo-aleatório
    numero = random(0,1000);
 
// Se ele for par, a gente acende a luz por 700 milisegundos de pausa
// Mas veja que antes eu incremento o contador de cara
// Imprimo o número de caras e pulo uma linha.
// Aqui eu apanhai um pouco porque achei que tudo era igual o printf do C
    if((numero % 2) == 0) {
      cara++;
      Serial.print("Cara = ");
      Serial.print(cara);
      Serial.print("\n");
      digitalWrite(carapin, HIGH);
      delay(pausa);
      digitalWrite(carapin, LOW);
      delay(250);
    } else {
      coroa++;
      Serial.print("Coroa = ");
      Serial.print(coroa);
      Serial.print("\n");
      digitalWrite(coroapin, HIGH);
      delay(pausa);
      digitalWrite(coroapin, LOW);
      delay(250);
    }
  }
 
//Depois de 100 vezes no loop acendendo e apagando luzes e contando
// Eu imprimo a contagem final de cara e coroa e reinicio as variáveis
// Antes de começar a contar denovo. 
  Serial.print("Final");
  Serial.print("\n");
  Serial.print("Cara = ");
  Serial.print(cara);
  Serial.print(" e Coroa =");
  Serial.print(coroa);
  Serial.print("\n");
  Serial.print("--------------------------");
  Serial.print("\n");
 
  cara = 0;
  coroa= 0;
  delay(1000);
 
}

Problema – Rosalind – Open Reading Frames – ORF

DNA_ORF

Esse problema do Rosalind será o com maior número de funções reutilizado desde o primeiro.
Mas vamos la, primeiro temos que entender do que diabos esse problema se trata, eu tive bastante problemas nessa fase, antes de começar a fazer qualquer coisa.

A gente ja viu aqui que podemos, a partir de uma sequência de RNA, prever qual proteína será produzida quando temos em mãos uma tabela que diz qual aminoácido um determinado códon produz. Inclusive vimos algumas propriedade que essa previsão proporciona aqui.

Certo, até aqui tudo bem, mas o RNA vem do DNA, logo a gente pode transcrever DNA em RNA como visto a muito tempo, aqui. E assim usar esse RNA e ver qual proteína ele codifica.

Mas para simplificar as coisas, temos tabelas, como para o RNA, dizendo o que um códon, três letrinhas de DNA geram de aminoácido, sabendo que antes o DNA é transcrito em RNA e depois vira uma proteína, ou algo da sequência de aminoácidos, posso estar sendo muito simplista ou errado nas nomenclaturas aqui.

Essas tabelas tem no wikipedia até, a do DNA para aminoácido aqui e do RNA para aminoácido aqui.

Agora vamos entrar no problema em si.

Acontece que nem todo DNA é transcrito para RNA, o tal do DNA lixo que as vezes a gente ouve falar. Então de uma sequência de DNA, em algum lugar pode começar o pedaço que vai ser transcrito, vai virar RNA, mas a gente pode não saber exatamente onde.

Como assim?

Vamos olhar a sequência de exemplo:

>Rosalind_99 AGCCATGTAGCTAACTCAGGTTACATGGGGATGACCCCGCGACTTGGATTAGAGTCTCTTTTGGAATAAGCCTGAATGATCCGAGTAGCATCTCAG

Certo, a gente tem essa sequência ai, que vai virar RNA, e depois o RNA em códons, que são sequências de letras vão virar aminoácidos que juntos formam proteínas.

Então se a gente separar os tripletos, os grupinhos de três letras aqui, a gente pode prever o aminoácido que vai ser formado.

Mas podemos começar aqui:

AGC - CAT - GTA - GCT - AAC - ...

Virão, começando do início da sequência acima.
Mas derrepente o RNA começa a ser transcrito da segunda base, ou seja assim:

(A) GCC - ATG - TAG - CTA - ACT - CAG - GTT - ...

Notem que os códon, os tripletos começam na segunda base, a primeira, sobrou, ficou ali como “lixo”, ela não influência em nada.
Aqui ficamos com duas possibilidades ja.

E existe uma terceira ainda…

(AG) CCA - TGT - AGC - TAA - CTC - AGG - TTA - ...

Aqui eu deixei em parenteses as duas primeiras bases, AG, elas não entraram no códon. São lixinho e agora formamos três possibilidades.
Será que daria para ter uma quarta?

(AGC) CAT - GTA - GCT - AAC - TCA - ...

Bem se fossemos pensar nessa quarta, agora o primeiro códon é igual ao segundo da nossa primeira opção, alias a sequência toda, so pulou o primeiro códon, por isso da quase na mesma.

Então para cada sequência de DNA, temos 3 possibilidades de Reading Frames, que é o nome chique para esse esquema de dependendo de qual ponto você começa, temos outra sequência de códons, ja que a quarta é igual a primeira pulando o primeiro códon, fica mais fácil pensar nos três primeiros casos so, ja que o quarto caso é o primeiro com menos códon, mas não paramos por aqui ainda.

Certo, só que o DNA não simplesmente faz uma proteína gigante, ele faz muitos tipos de proteínas e outras coisas. Acontece que tem um códon que inicia a proteína. No caso dos Eucariotos é comumente a Metionina (M), que é codificado por “ATG” no DNA e “AUG” no RNA.
Então naqueles Reading Frames ali em cima, das possibilidades, quando aparece um “ATG” a sequência de aminoácido começa a ser feita, e so para na hora que chegar no num códon de parada, so olhar na tabela quais são as paradas.

Mas para complicar mais ainda, além dessas possibilidades, a gente tem duas fitas, o complemento do DNA. Como a gente viu aqui, podemos ver o complemento dessa sequência e ai tem mais três Reading Frames para fazer proteína, é um bocado de possibilidades né.

Agora vamos ver quais são essas possibilidades, como eu disse, a gente vai reciclar um monte de código.

Primeiro, como já falamos ali em cima, como vamos ter que olhar a sequência e seu complemento, temos que definir uma função para ver o reverso do complemente, e isso ja esta pronto aqui:

def revc(sequencia):
    intab = "ACTG"
    outtab = "TGAC"
    trantab = maketrans(intab, outtab)
    return(sequencia.translate(trantab)[::-1])

Lembrando que para usar a função maketrans, tem que abrir o library string (“from string import maketrans”)

A tabela códon -> aminoácido o exercício da, e tem no wikipedia, então sem grandes problemas.

Antes de começar, temos ainda outra atividade que vai sempre ser refeita, que é traduzir os códons, transformar a sequência de códon em aminoácidos, de três em três. Basicamente essa tarefa aqui, so mudar a tabela da de RNA para DNA.

def traduzir_codon(codon):
    amino = ""
    if len(codon) == 3 and DNA_CODON_TABLE.has_key(codon):
        amino = DNA_CODON_TABLE[codon]
    return amino

Agora, aqui tem um problema que eu pedalei muito, ao fazer uma função para passar o DNA para proteína, a gente tem que ter o cuidado de so fazer isso quando tem um códon de stop. Pelo que vi de duvidas no question, foi muita gente além de mim que pedalou aqui.

Como assim?
Bem, pode acontecer de quando a sequência estiver em algum ponto, como por exemplo perto do fim, a gente tem um codon “ATG” para começar, beleza, mas vamos adicionando aminoácidos depois dele, e a sequência acaba e não vimos um códon de stop, ou seja, talvez esteja mais para a frente na sequência e o primer cortou no meio essa sequência, não pegamos um pedaço completinho, sei la, mas se não tem stop, só começa, eu retornei um string vaziu, ta la o ultimo return “” da função, quando chega em um stop, a proteína é retornada. Assim depois eu so vou contabilizar como possibilidade o que não for “”, vazio.

def dna2prot(seqdna,inicio):
    result = ""
    for i in range(inicio, len(seqdna), 3):
        aminoacido = traduzir_codon(seqdna[i:i+3])
        if aminoacido == 'Stop':
            return result
        else:
            result += aminoacido
    return ""

Agora, com tudo definido, é só usar esse monte de função para chegar as possibilidades:

def possibilidades(sequencia):
    possibilidades = []
    inicio = sequencia.find("ATG",0)
    while inicio != -1:
        resul = dna2prot(sequencia,inicio)
        inicio = sequencia.find("ATG",inicio+1)
        if resul != "":
            possibilidades.append(resul)
    return possibilidades

O que estramos fazendo aqui? Começamos uma lista vazia, achamos o primeiro “ATG”, dai a partir dele vemos que proteina é codificada com a função dna2prot, terminando isso, achamos o próximo “ATG” do string, note que usamos a função find aqui, e depois que achamos o primeiro caso, começamos a procurar do inicio+1, para achar o próximo inicio, o próximo “ATG”, so que como eu to lendo letra por letra, eu estou contemplando os três reading frames.
Mas voltando, eu criei uma lista vazia, achei o primeiro “ATG”, comecei um loop, testo se ainda tem “ATG”, se não tiver o find vai dar um resultado de -1, se tiver, ele vai dizer onde começa. Como tem, o loop começa, traduzo a primeira proteína, acho o próximo “ATG” e olho a proteína que saiu, se ela não for vazia, não for “”, ela ponho ela na lista, e começo o loop denovo, e fico nessa até acabarem os “ATG”.

Agora eu tenho todas as possibilidades para essa sequência, mas, eu ainda preciso das possibilidades do complemento, mais como a gente ja definiu uma função para fazer o complemento da sequência, é baba, é so usar ela e depois ver as possibilidades no complemento.

Para não repetir possibilidades ao somar as possibilidades da sequencia e do complemento dela, usamos a função set do python para termos somente termos únicos.

Como fazia tempo que eu não me preocupava em abrir arquivos txt no python, eu sofri um pouco aqui. Mas eu fiz o seguinte:

arquivo = open('/home/augusto/Downloads/rosalind_orf.txt')
 
dados = arquivo.readline()
dataset=""
 
for line in arquivo:
        dataset+=line.strip()

Primeiro com open a gente abre o arquivo, readline() le uma linha do arquivo, como a primeira não me interessa eu coloco ela em algum lugar, as linhas são lidas de forma sequencial com readline(), então tudo que sobrou é a sequência nas linhas quebradas, a primeira linha era a identificação da sequencia que não me interessava, dai então para as linhas no arquivo, eu coloco elas no meu dataset, usando o .strip depois para tirar espaços, que não é o caso, mas para tirar o \n que quebram linhas, que se isso entrar nas funções la em cima vira uma merda.

Dai o dataset fica um string gigante de toda a sequência.

Bem uma observação legal, é que esse esquema de usar find, para achar o próximo caso, num loop igualzinho aquele la, eu vi assistindo as aulinhas de introdução a computação do udacity. La era ensinado como fazer um web-crawler. Um programinha que lia uma página da internet, o html, dai ia olhando o texto até achar um link, guardava o link e ia procurar o próximo, e assim até acabar o página salva. Ora bolas, não é exatamente o que estamos fazendo aqui, so que ao invés de links, paginas de internet são proteínas.
Eu acho legal como o mesmo problema pode de certa forma se repetir em áreas que pra mim parecem tão distintas.

Outra coisa legal, desse exercício é pensar como uma mutação em alguns lugares vai fazer um fuzue danado, pode mudar muito as possibilidades de proteínas que obtemos, enquanto em outros lugares não vai influencias em nada, vai ser completamente neutra, um random walk.

Bem mas é isso, o script vai estar la no repositório do Rosalind além de aqui em baixo.
Foi um problema bem legal, apesar de ter apanhando muito para resolver, principalmente no final onde o problema estava em somente ler os dados de um txt, e não em resolver o problema, é muito frustante conseguir resolver certinho o exemplo e depois sempre falhar com os dados do teste, mas nada que um pouco de perseverança não resolva.

A saída do script abaixo é:

>>> MLLGSFRLIPKETLIQVAGSSPCNLS M MGMTPRLGLESLLE MTPRLGLESLLE

E sem mais enrolar, el script:

from string import maketrans
 
DNA_CODON_TABLE = {
    'TTT': 'F',     'CTT': 'L',     'ATT': 'I',     'GTT': 'V',
    'TTC': 'F',     'CTC': 'L',     'ATC': 'I',     'GTC': 'V',
    'TTA': 'L',     'CTA': 'L',     'ATA': 'I',     'GTA': 'V',
    'TTG': 'L',     'CTG': 'L',     'ATG': 'M',     'GTG': 'V',
    'TCT': 'S',     'CCT': 'P',     'ACT': 'T',     'GCT': 'A',
    'TCC': 'S',     'CCC': 'P',     'ACC': 'T',     'GCC': 'A',
    'TCA': 'S',     'CCA': 'P',     'ACA': 'T',     'GCA': 'A',
    'TCG': 'S',     'CCG': 'P',     'ACG': 'T',     'GCG': 'A',
    'TAT': 'Y',     'CAT': 'H',     'AAT': 'N',     'GAT': 'D',
    'TAC': 'Y',     'CAC': 'H',     'AAC': 'N',     'GAC': 'D',
    'TAA': 'Stop',  'CAA': 'Q',     'AAA': 'K',     'GAA': 'E',
    'TAG': 'Stop',  'CAG': 'Q',     'AAG': 'K',     'GAG': 'E',
    'TGT': 'C',     'CGT': 'R',     'AGT': 'S',     'GGT': 'G',
    'TGC': 'C',     'CGC': 'R',     'AGC': 'S',     'GGC': 'G',
    'TGA': 'Stop',  'CGA': 'R',     'AGA': 'R',     'GGA': 'G',
    'TGG': 'W',     'CGG': 'R',     'AGG': 'R',     'GGG': 'G'
}
 
def revc(sequencia):
    intab = "ACTG"
    outtab = "TGAC"
    trantab = maketrans(intab, outtab)
    return(sequencia.translate(trantab)[::-1])
 
def traduzir_codon(codon):
    amino = ""
    if len(codon) == 3 and DNA_CODON_TABLE.has_key(codon):
        amino = DNA_CODON_TABLE[codon]
    return amino
 
def dna2prot(seqdna,inicio):
    result = ""
    for i in range(inicio, len(seqdna), 3):
        aminoacido = traduzir_codon(seqdna[i:i+3])
        if aminoacido == 'Stop':
            return result
        else:
            result += aminoacido
    return ""
 
def possibilidades(sequencia):
    possibilidades = []
    inicio = sequencia.find("ATG",0)
    while inicio != -1:
        resul = dna2prot(sequencia,inicio)
        inicio = sequencia.find("ATG",inicio+1)
        if resul != "":
            possibilidades.append(resul)
    return possibilidades
 
exemplo = "AGCCATGTAGCTAACTCAGGTTACATGGGGATGACCCCGCGACTTGGATTAGAGTCTCTTTTGGAATAAGCCTGAATGATCCGAGTAGCATCTCAG"
#arquivo = open('/home/augusto/Downloads/rosalind_orf.txt')
 
#dados = arquivo.readline()
#dataset=""
 
#for line in arquivo:
#        dataset+=line.strip()
 
#print str(dataset)
 
 
A = possibilidades(exemplo)
B = possibilidades(revc(exemplo))
 
print "\n".join(set(A + B))