Programação multi-thread em python

Esse talvez possa parecer um post meio maluco, para o que comumente discutimos aqui nesse blog, mas é muito interessante e curioso, para saber o que está acontecendo sobre o capo do carro, ou do computador.

Uma das coisas mais diferentes que já vi, foi quando fiz a disciplina de sistemas operacionais, isso porque tudo que eu achava que sabia sobre computador, vi que estava errado. Nem tudo é tão simples como um script continuo que fazemos aqui no R, e uma prova disso é a programação Multi-thread.

thread

Vamos pensar na figura acima, tomando como exemplo os MCMC que vemos comumente aqui, veja alguma exemplos de implementação aqui de um gibbs samples e aqui do mais simples metropolis hastings.

Nesses dois casos fizemos algo como em single-threaded, isso porque fizemos uma cadeia e fomos realizando vários cálculos e salvando resultados. Acontece que, todas as analises de dados que fazemos, não usamos apenas uma cadeia, usamos várias cadeias, no último post mesmo, usando o openbugs, usamos três cadeias, que podem ser processadas separadamente, podem ser multithreaded.
Os valores das amostras, parâmetros, modelos, tudo é igual, então não compensa gravar isso varias vezes na memoria, é melhor compartilhar, mas para uma cadeia de markov, a gente sempre precisa do valor anterior para continuar ela, então as cadeias não precisam dividir os cálculos que estão fazendo entre si, apenas os dados e parâmetros e modelos, agora olha no modelinho ali em cima, não é exatamente isso que diz respeito a programação multi-thread?
Veja que com multi-thread, cada thread tem sua pilha, seus registradores, seus dados particulares, que é sua cadeia de markov que está computando, mas também divide dados, que são as amostras, modelo, parâmetros entre todas as threads.

E qual a relevância disso? Ai que entre aqueles processadores de muitos núcleos, sem isso, sem paralelizar, dividir a computação, a gente teria que fazer uma cadeia no processador, terminar, fazer outra, terminar, pense que leva um minuto por cadeia, com três então levamos três minutos para fazer tudo, mas se fizermos algo multi-thread, da para mandar cada cadeia para cada núcleo do processador (mais ou menos isso, simplificando, sei que não é assim), como cada um leva 1 minuto e da para fazer em paralelo, ao invés de esperar 3 minutos pelo resultado, esperamos um.

image0021225707175745

Podemos fazer um exemplo em python que é bem simples, usando a biblioteca threading.

Para esse exemplo simples, vamos criar uma classe bem bobinha

1
2
3
4
5
6
7
8
9
10
class minhaThread (threading.Thread):
    def __init__(self, threadID, nome, contador):
        threading.Thread.__init__(self)
        self.threadID = threadID
        self.nome = nome
        self.contador = contador
    def run(self):
        print "Iniciando thread %s com %d processos" % (self.name,self.contador)
        processo(self.nome, self.contador)
        print "Finalizando " + self.nome

A gente ja falou de orientação a objetos em R aqui e aqui, uma das coisas que a orientação a objetos permite é a herança, que é herdar tudo que a classe pai tem, aqui estamos herdando a classe Thread definida na biblioteca threading, isso é o que esta escrito aqui:

1
class minhaThread (threading.Thread):

Por isso o parenteses, lembro que ja usei classe em algum problema do rosalind aqui, mas não lembro em qual post agora, mas beleza, o resto é bem besta a classe vai ter 3 atributos, um ID, um nome e um contador, que seria algo como os parâmetros para fazer um mcmc, vamos abstrair de uma atividade útil agora, vamos apenas pensar em termos de computação.

1
2
3
4
5
    def __init__(self, threadID, nome, contador):
        threading.Thread.__init__(self)
        self.threadID = threadID
        self.nome = nome
        self.contador = contador

Então nos construtor da thread definimos esses atributos e ela vai ter apenas um método, uma função chamada run.

1
2
3
4
    def run(self):
        print "Iniciando thread %s com %d processos" % (self.name,self.contador)
        processo(self.nome, self.contador)
        print "Finalizando " + self.nome

Tem que chamar run, por causa que é assim que funciona o pacote, mas o run vai so imprimir que estamos começando, processar algo, pensando que processar pode ser qualquer coisa e depois avisar que terminamos.

Aqui o processo definimos assim

1
2
3
4
def processo(nome, contador):
    while contador:
        print "Thread %s fazendo o processo %d" % (nome, contador)
        contador -= 1

Ele é uma função externa, não é da classe, mas ele apenas imprimi uma mensagem, falando que esta processando, e diminui o contador em 1, até zerar, e terminar.

Depois disso está tudo pronto, so precisamos criar nossas threads

1
2
3
# Criando as threads
thread1 = minhaThread(1, "Alice", 8)
thread2 = minhaThread(2, "Bob", 8)

Criamos elas, simples assim, aqui no caso são duas, chamadas de alice e bob, hehe, alice e bob são exemplos mais comuns em física, criptografia e teoria dos jogos, ao invés de falar A e B, os caras falam Alice e Bob, para ficar mais fácil de entender.

Ok, mas apos criado, nada foi processado ainda, dai então começamos elas.

1
2
3
# Comecando novas Threads
thread1.start()
thread2.start()

O resto do programa é meramente para fazer a thread mãe esperar alice e bob para terminar, não convêm explicar agora, mas essa parte so faz isso, esperar alice e bob terminarem

1
2
3
4
5
6
threads = []
threads.append(thread1)
threads.append(thread2)
 
for t in threads:
    t.join()

Agora é legal quando a gente roda esse programa.

As vezes ele roda uma thread seguida da outra

Screenshot from 2015-07-23 16:21:42

E as vezes tudo vira uma loucura.

Screenshot from 2015-07-23 16:23:00

Isso porque uma vez que você inicio o processamento, o sistema operacional que vai escalonar as threads, dai a gente não manda mais, só o sistema operacional, então se o processador ta liberado, ele manda elas para la para fazer conta, mas as vezes ta ocupado, então a thread tem que esperar sua vez, agora uma thread pode ir para um núcleo e ir até o final la, enquanto a outra entrou num outro e teve que esperar alguém terminar de processar, como no segundo exemplo, Alice começou, o Bob esperou um pouco, Alice ja tinha processado 2 vezes, ai o Bob começou a trabalhar, mas terminou depois.

Um exemplo bem besta, mas não é preciso saber fazer muitas contas complexas para imaginar que isso é melhor que primeiro alice fazer todo o trabalho dela e depois bob fazer todo o trabalho dela. Pena que isso é extremamente complexo, ainda mais quando a gente tem que compartilhar algum dado e mexer nele, porque ai a vez de quem mexeu pode ser importante, mas isso é outra historia.

Claro que programação não é a meta de muitas pessoas que visitam aqui, mas ter uma ideia que isso existe, pode ajudar a entender comentários de manuais de programa entre outras coisas.

Por exemplo, um programa comum em biologia é o mrBayes, basta olhar o manual que você vai notar alguma notas de instalação, veja essa parte aqui:

instructions for compiling and running the parallel version of the program

Ou seja, se você instala as coisas so dando dois cliques, não le o manual, pode estar perdendo essa possibilidade de usar multi-thread e algo que pode rodar em uma hora vai levar um dia. Ja que ainda no exemplo do mrBayes, da para usar até o processador da placa de video para ajudar a fazer cálculos e agilizar o processamento. E pode abrir o olho, que tudo que você ler MCMC é um forte candidato a usar multi-thread.

O R em si não é multi-thread, mas existem vários pacotes que tentam incorporar isso, como visto no taskview de high performance. Além de que o R possui interfaces que ligam eles a muitos outros programas, como é o caso do Openbugs mesmos, que os calculos são feitos neles, mas depois devolvemos os dados no R para fazer figuras e outras analises, mas temos muitos outros exemplo.

Apesar que na maioria das vezes isso não é muito relevante, para fazer regressões e muitos modelos, tudo é tão rápido que a gente nem sente diferença, mas algumas coisas mais avançadas (talvez “avançando” seja um termo ruim aqui, coisas que precisam fazer mais contas, comum em evolução), isso começa a ser importante.

Bem é isso ai, a primeira vez que ouvi sobre esse negocio de multi-thread, a primeira coisa que veio na minha cara é MCMC, acho que ter ouvido falar de MCMC antes de fazer Sistemas Operacionais foi um diferencial para tornar tudo mais interessante e belo, mas tive um excelente professor em sistemas operacionais também para ajudar. O script vai estar la no repositório recologia, tem mais alguma exemplos da minha aula de sistemas operacionais no github em outro repositório aqui mas em linguagem C, e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail e até mais.

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
import threading
 
class minhaThread (threading.Thread):
    def __init__(self, threadID, nome, contador):
        threading.Thread.__init__(self)
        self.threadID = threadID
        self.nome = nome
        self.contador = contador
    def run(self):
        print "Iniciando thread %s com %d processos" % (self.name,self.contador)
        processo(self.nome, self.contador)
        print "Finalizando " + self.nome
 
def processo(nome, contador):
    while contador:
        print "Thread %s fazendo o processo %d" % (nome, contador)
        contador -= 1
 
# Criando as threads
thread1 = minhaThread(1, "Alice", 8)
thread2 = minhaThread(2, "Bob", 8)
 
# Comecando novas Threads
thread1.start()
thread2.start()
 
threads = []
threads.append(thread1)
threads.append(thread2)
 
for t in threads:
    t.join()
 
print "Saindo da main"

Data Augmentation para dados de Captura e Recaptura com efeito do tempo

Ano passado, a gente olhou aqui a ideia de data augumentation. Que é muito interessante. Ai eu queria fazer posts sobre todos esses modelos, que são muito legais, mas fiquei enrolando, mas agora sai mais um e sabe-se la quando o resto.

Mas vamos la, esse modelo assume que a probabilidade de detecção p varia por ocasião de coleta, se estivermos coletando passarinhos, talvez a chuva faça com que eles não saiam para voar e cair na rede, ou para pequenos mamíferos, as condições do clima estejam ruins, o dia para ficar na cama abaixo das cobertas, ou ainda diferentes armadilhas ou mecanismos de detecção foram usados.

Então vamos gerar dados, pensando em passarinhos caindo numa rede, que colocamos no campo, e conseguimos ir pegando os passarinhos e anilhando eles, conforme caem na rede, então identificamos eles perfeitamente.

Então imagine:

1
2
3
4
n = 100
probabilidade_detecao = 0.25
tempo = 5
variacao_Tempo = runif(tempo,-2,2)

Cem passarinhos (n), coletados em 5 campanha diferentes (tempo), mas além desses parâmetros, sabemos que… , bem nos não sabemos, mas a mãe natureza estabelece uma probabilidade desse passarinho cair na rede (probabilidade_detecao), mas essa probabilidade não é constante ao longo das campanhas, as vezes ela é mais um pouquinho, quando o clima ta bom, sem vento, e as vezes menos um pouquinho (esses valores estão no vetor variacao_Tempo).

Ok, agora vamos fazer a matriz de dados que a mãe natureza vê, ela vê todos dos 100 passarinhos que la residem, e ela observa as 5 campanhas também, ou seja uma matriz de n pelo tempo.

1
yreal<- matrix(NA, nrow=n,ncol=tempo)

Para simular os dados fazemos o seguinte

1
2
3
4
5
6
7
#Num loop para todas as campanha
for (j in 1:tempo){
    #Qual a probabilidade de detecção nessa campanha j?
    p <- plogis(log(probabilidade_detecao / (1-probabilidade_detecao)) + variacao_Tempo[j])
    #Ai simulamos todos os 100 individuos de uma vez
    yreal[,j] <- rbinom(n = n, size = 1, prob = p)
}

E pronto, a simulação é isso, gerar números a partir de uma distribuição binomial, na verdade.
Só que quem a gente nunca viu em nenhuma campanha, quem nunca foi anilhado não existe pra gente, só foi anilhado quem foi pego pelo menos uma vez

1
2
detectado_pelo_menos_uma_vez <- apply(yreal, 1, max)
C <- sum(detectado_pelo_menos_uma_vez)

E assim esses são os cara que observamos

1
2
yobservado <- matrix(NA, nrow=n,ncol=tempo)
yobservado <- yreal[detectado_pelo_menos_uma_vez == 1,]

Agora basta a gente estabelecer nosso modelo para a linguagem bugs

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
 model {
        # Priors
        omega ~ dunif(0, 1)
        for (i in 1:tempo){
            p[i] ~ dunif(0, 1)
        }
        # Likelihood
        for (i in 1:amostras){
            z[i] ~ dbern(omega)
            for (j in 1:tempo){
                yaug[i,j] ~ dbern(p.eff[i,j])
                p.eff[i,j] <- z[i] * p[j]
                } #j
            } #i
        # Quantidade derivadas
        N <- sum(z[])
        p.medio <- mean(p[])
    }

E a principio o likelihood pode parecer idêntico ao do post anterior, mas veja que a diferença está em usar um vetor para as probabilidade de detecção da campanha

1
p.eff[i,j] <- z[i] * p[j]

Veja que no prior, inicializamos um vetor agora, e não um único valor, para contabilizar essa variabilidade entre campanha

1
2
3
for (i in 1:tempo){
     p[i] ~ dunif(0, 1)
}

Ai agora é só estimar os parâmetros, temos que colocar os dados no formado do openbugs que estamos usando aqui, lembrando das outras opções aqui, estabelecer valores iniciais para as cadeia do MCMC, seu parâmetros como tamanho, valores iniciais a serem queimados, número de cadeia e rodar.

Aqui o resultado fica o seguinte:

> print(out_tempo, dig = 3) Inference for Bugs model at "modelo_tempo.txt", Current: 3 chains, each with 2500 iterations (first 500 discarded), n.thin = 2 Cumulative: n.sims = 6000 iterations saved mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff N 99.604 4.646 92.000 96.000 99.000 102.000 110.000 1.001 6000 p[1] 0.542 0.055 0.436 0.505 0.542 0.579 0.649 1.001 6000 p[2] 0.108 0.031 0.055 0.086 0.105 0.127 0.179 1.001 6000 p[3] 0.138 0.034 0.078 0.113 0.135 0.160 0.211 1.001 6000 p[4] 0.672 0.056 0.560 0.635 0.673 0.711 0.776 1.001 6000 p[5] 0.099 0.030 0.049 0.078 0.096 0.117 0.164 1.001 6000 p.medio 0.312 0.022 0.268 0.297 0.312 0.327 0.354 1.001 6000 omega 0.418 0.037 0.349 0.392 0.417 0.443 0.495 1.002 2700 deviance 470.030 20.741 434.300 455.500 468.300 483.100 514.302 1.001 6000 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 = var(deviance)/2) pD = 215.1 and DIC = 685.2 DIC is an estimate of expected predictive error (lower deviance is better).

Bem existem mais ou menos entre 92 a 110 indivíduos por ali, o intervalo de confiança, o que cobre a verdade (já que sabemos os valores usados pela mãe natureza foi 100), e temos um valor de detecção por campanha e duas quantidades derivadas, o p médio e o omega, que falamos depois.

Podemos tentar olhar esse resultado numa figura das nossas estimativas.

01

Beleza, mas será que valeu a pena o esforço, o que aconteceria se a gente ajusta-se o modelo considerando um único valor de detecção. Bem, podemos ajustar e ver o que acontece, por curiosidade.

E a saída é o seguinte

> print(out, dig = 3) Inference for Bugs model at "modelo.txt", Current: 3 chains, each with 2500 iterations (first 500 discarded), n.thin = 2 Cumulative: n.sims = 6000 iterations saved mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff N 112.268 7.386 100.000 107.000 112.000 117.000 129.000 1.002 1700 p 0.274 0.026 0.224 0.256 0.274 0.292 0.326 1.002 1700 omega 0.471 0.044 0.390 0.440 0.469 0.499 0.561 1.002 1400 deviance 657.828 23.208 616.100 641.300 656.900 672.425 707.800 1.002 1700 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 = var(deviance)/2) pD = 269.1 and DIC = 926.9 DIC is an estimate of expected predictive error (lower deviance is better).

Aqui o intervalo de confiança ainda pegou o valor real, que era de 100, o tamanho populacional teve uma estimativa de 100 a 129 indivíduos, o que está ok, com um p de 274, a probabilidade média de captura. Agora o que isso quer dizer?
Vamos olhar os valores da mãe natureza.

Veja que os valores de detecção foram

> plogis(log(probabilidade_detecao/(1-probabilidade_detecao))+variacao_Tempo) [1] 0.50132242 0.09139408 0.12948952 0.69882099 0.09683113

Da média original, na primeira campanha foi bem mais, na segunda e terceira bem menos, na quarta a mais, e na quinta a menos. se fizermos a média disso da

> mean(plogis(log(probabilidade_detecao/(1-probabilidade_detecao))+variacao_Tempo)) [1] 0.3035716

O valor que aquele p realmente queria chegar era 0.304, mas ele ficou em 0.274, isso quer dizer, que na segunda, terceira e quinta campanha, a gente deveria ter pensado em achar menos passarinhos, mas o modelo fez simulação com valores altos, sempre achando mais, por isso a gente super-estima muito, veja que comparando as duas distribuições, nesse exemplo o modelo sem levar o tempo foi la para frente.

figura2

Veja que isso é a quantidade derivada la do modelo com o tempo, que ficou próximo desse valor real.

p.medio 0.312 0.022 0.268 0.297 0.312 0.327 0.354 1.001 6000

Um pouquinho maior, mas próximo.

Veja que isso vai depender claro do número de campanhas, pensando que a não ser que haja uma tendência nesses valores de detecção, por exemplo, a não ser que eles sempre melhorem, ou piorem, se eles simplesmente forem um processo estocástico com média zero, com muitas campanhas, os modelo sem o tempo conseguiria um valor bom de estimativa populacional, mas basta pensar um pouco para saber que isso não vai acontecer, porque?

Porque campanhas de coletas são caras, dificilmente o dinheiro vai dar para 30 ou 90 campanhas, e com poucas campanhas, a chance de você considerar um valor fixo de detecção e dar errado é maior, chance no sentido de que o modelo sempre vai ser pior, então talvez começar com um modelo assim seja mais seguro, e é fácil você simplificar se for possível.

Veja o deviance do modelo sem o tempo é muito pior:

deviance 657.828 23.208 616.100 641.300 656.900 672.425 707.800 1.002 1700

Compare como o deviance do modelo com o tempo é muito melhor

deviance 470.030 20.741 434.300 455.500 468.300 483.100 514.302 1.001 6000

Veja o quanto, quando comparamos as estimativas, o desvio padrão das estimativas para o modelo com tempo é menor.

figura3

A importante contribuição da incerteza sobre o tamanho da população vem da incerteza sobre a probabilidade p. Só que a gente paga em parâmetros a mais que precisamos para estimar p. Mas ignorar esses parâmetros levam a erros, que podem ser muito ruins. Talvez se a variação fosse pequena, esse modelo nem vale-se a pena, mas na dúvida, compensa começar do modelo mais complexo, e ir simplificando, vendo se é possível simplificar, afinal, se todas as probabilidades de detecção fossem iguais, daria para ver nos parâmetros estimados e não precisaríamos de um vetor de probabilidade de detecção, mas pero que si, pero que no, compensa tentar todas as possibilidades.

Nesse modelo, nos assumimos que os efeitos do tempo são fixos, mas poderíamos ainda ajustar um modelo estocástico também, pensando que ao invés de valores fixos de p, eles viriam de uma distribuição como a normal ou uniforme, que é o caso dos dados que geramos aqui.

Bem é isso ai, o script vai estar la no repositório recologia, falta uns dois modelos a mais sobre esse tema no livro muito legal do Kery, mas um dia chegamos la, isso tem no capítulo 6. Se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail e até mais.

Referência:
Kéry e Shaub 2011 – Bayesian population analysis using WinBUGS. A hierarchical perspective – Academic Press 554pp

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
#
library(R2OpenBUGS)
 
##################################
## Gerando dados
##################################
set.seed(51)
n = 100
probabilidade_detecao = 0.25
tempo = 5
variacao_Tempo = runif(tempo,-2,2)
 
yreal<- matrix(NA, nrow=n,ncol=tempo)
 
for (j in 1:tempo){
    p <- plogis(log(probabilidade_detecao / (1-probabilidade_detecao)) + variacao_Tempo[j])
    yreal[,j] <- rbinom(n = n, size = 1, prob = p)
}
 
detectado_pelo_menos_uma_vez <- apply(yreal, 1, max)
C <- sum(detectado_pelo_menos_uma_vez)
 
yobservado <- matrix(NA, nrow=n,ncol=tempo)
yobservado <- yreal[detectado_pelo_menos_uma_vez == 1,]
 
##################################
## Modelo com tempo
##################################
 
# Data Augumentation
nz <- 150
yaug <- rbind(yobservado, array(0, dim = c(nz, tempo)))
 
# Modelo na linguagem bugs
sink("modelo_tempo.txt")
cat("
    model {
        # Priors
        omega ~ dunif(0, 1)
        for (i in 1:tempo){
            p[i] ~ dunif(0, 1)
        }
        # Likelihood
        for (i in 1:amostras){
            z[i] ~ dbern(omega)
            for (j in 1:tempo){
                yaug[i,j] ~ dbern(p.eff[i,j])
                p.eff[i,j] <- z[i] * p[j]
                } #j
            } #i
        # Quantidade derivadas
        N <- sum(z[])
        p.medio <- mean(p[])
    }
    ",fill = TRUE)
sink()
 
# Juntando os dados
bugs.data <- list(yaug = yaug, amostras = nrow(yaug), tempo = ncol(yaug))
 
# Valores iniciais
inits <- function() list(z = rep(1, nrow(yaug)), p = runif(tempo, 0, 1))
 
# Parametros para saida
params <- c("N", "p","p.medio", "omega")
 
# Configurações do MCMC
ni <- 2500
nt <- 2
nb <- 500
nc <- 3
 
#Rodando o modelo
out_tempo <- bugs(bugs.data, inits, params, "modelo_tempo.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb)
 
# Resultados
print(out_tempo, dig = 3)
 
#Figura1
hist(out_tempo$sims.list$N,nclass=30, col = "gray", main = "", xlab = "Tamanho Populacional", las = 1, xlim = c(80, 120))
abline(v = C, col = "black", lwd = 3)
abline(v = n, col = "red", lwd = 3)
legend("topright",c("Detectado","Real"),lty=1,lwd=3,col=c("black","red"),bty="n")
 
##################################
## Modelo sem variação de tempo
##################################
 
# Modelo na linguagem bugs
sink("modelo.txt")
cat("
    model {
        # Priors
        omega ~ dunif(0, 1)
        p ~ dunif(0, 1)
        # Likelihood
        for (i in 1:amostras){
            z[i] ~ dbern(omega)
            for (j in 1:tempo){
                yaug[i,j] ~ dbern(p.eff[i,j])
                p.eff[i,j] <- z[i] * p
                } #j
            } #i
        # Quantidade derivadas
        N <- sum(z[])
    }
    ",fill = TRUE)
sink()
 
# Valores iniciais
inits <- function() list(z = rep(1, nrow(yaug)), p = runif(1, 0, 1))
 
#
params <- c("N", "p", "omega")
 
#Rodando o modelo
out <- bugs(bugs.data, inits, params, "modelo.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb)
 
# Resultados
print(out, dig = 3)
print(out_tempo, dig = 3)
 
#Figura 2
hist(out_tempo$sims.list$N, nclass = 40, col = "darkgray", main = "", xlab = "Tamanho Populacional", las = 1, xlim = c(70, 150))
hist(out$sims.list$N, nclass = 40, col = "lightgray",add=T)
abline(v = mean(out_tempo$sims.list$N), col = "blue", lwd = 3)
abline(v = mean(out$sims.list$N), col = "red", lwd = 3)
legend("topright",c("Modelo com tempo","Modelo sem tempo"),lty=1,lwd=3,col=c("blue","red"),bty="n")
 
#
plogis(log(probabilidade_detecao/(1-probabilidade_detecao))+variacao_Tempo)
mean(plogis(log(probabilidade_detecao/(1-probabilidade_detecao))+variacao_Tempo))
 
#Figura 3
par(mfrow=c(2,1))
hist(out_tempo$sims.list$N, nclass = 40, col = "darkgray", main = "Modelo com variação no tempo",
     xlab = "Tamanho Populacional", las = 1, xlim = c(70, 150),ylim=c(0,600))
text(130,400,paste("Desvio padrão =",round(sd(out_tempo$sims.list$N),3)))
hist(out$sims.list$N, nclass = 40, col = "darkgray", main = "Modelo sem variação no tempo",
     xlab = "Tamanho Populacional", las = 1, xlim = c(70, 150),ylim=c(0,600))
text(130,400,paste("Desvio padrão =",round(sd(out$sims.list$N),3)))

Rosalind – Locating Restriction Sites

Na guerra invisível entre bactérias e seus agressores, os vírus bacteriófagos, ou “phages”, uma das armas das das bactérias são as enzimas de restrição.

O que o vírus faz? Ele introduz seu material genético na bactéria, insere seu dna no dna da bactéria e a partir dai ela copia esse dna, ou seja o vírus usa todo o maquinário da bactéria para funcionar.

E como a bactéria pode combater isso? Atacar o dna do vírus não é simples. Pensa que qualquer coisa que ataque DNA, ataca ela mesmo. Esse é o triunfo de vírus e doenças como câncer, falta de especificidade para conseguir atacar o problema.

Mas anos de evolução produziram enzimas de restrição eficientes e especificas o suficientes para ferrar com o vírus e não atrapalhar a bactéria. Ou seja, você deixa flutuando um dna um monte de enzima que corta dna, mas corta num local suficientemente especifico para ferrar somente o vírus, dai se phage adiciona algum dna la dentro da bactéria, ele pode ser incapacitado antes de fazer algum mal.

Ok, enzimas de restrição podem ser classificados em vários tipo, mas um dos tipos, o tipo 2 são os homodimers, que além de muitas características, são palíndromos no complemento, de 4 a 12 bases de tamanho.

Então nesse exercício temos que encontrar esses locais possíveis, em uma dada sequencia, que podemos assumir que será de um vírus certo? Mas apesar da especificidade do local do “corte”, podemos pensar que 4 pares de base é 0.25^4=0.00390625, pensando numa sequencia de Mycobacteirum bovis por exemplo que tem mais ou menos 4.300.000 de pares de base, partindo isso em 4, não contando que uma base que não foi do primeiro conjunto de 4 pode contribuir para o próximo, esperaríamos pelo menos (4300000/4)*(0.25^4)=4199 locais ao acaso para essa enzima de restrição funcionar, mas esse número é maior se fizermos as contas direto, mas a bactéria ainda tem o dna methylation para se salvar, que é uma das formas de controlar expressão genica, mas não vem ao caso.

Agora sim, ao problema. Temos a seguinte sequência.

>Rosalind_24 TCAATGCATGCGGGTCTATATGCAT

Veja que primeiro temos que saber fazer o complemento dela, que é

>Rosalind_24 TCAATGCATGCGGGTCTATATGCAT AGTTACGTACGCCCAGATATACGTA Complemento

E invertermos este complemento.

>Rosalind_24 TCAATGCATGCGGGTCTATATGCAT ATGCATATAGACCCGCATGCATTGA Complemento Reverso

Certo, ja tinha visto isso num problema anterior aqui, que diga-se de passagem, meu código era bem mais confuso hehe.

Agora, temos que encontrar matchs de 4 a 12 pares de bases nessa sequência
O match é onde elas são iguais, o primeiro acontece no quarto par de base

Original TCA ATGCAT GCGGGTCTATATGCAT Complemento TACGTA Reverso ATGCAT

Logo

Original TCA ATGCAT GCGGGTCTATATGCAT Reverso ATGCAT

Ohhh, que coincidência. Mas veja, não é para fazer o complemento de toda a sequência, como la em cima, mas o complemento de cada pedacinho.
Veja que o próximo caso, está “dentro” desse

Original TCAA TGCA TGCGGGTCTATATGCAT Complemento ACGT Reverso TGCA

Acho que como o dna é uma fita dupla, uma enzima de restrição consegue pegar logo as duas fitas de uma vez.

Agora como fazer isso por força bruta?

Pegue a sequência

TCAATGCATGCGGGTCTATATGCAT

Separe o primeiro candidato

TCAA TGCATGCGGGTCTATATGCAT

Faça o complemento reverso

OR:TCAA CR:TTGA

Pergunta se eles são iguais, se sim salve a posição e o tamanho, senão passa para o próximo caso, que vai ser

TCAAT GCATGCGGGTCTATATGCAT

Aumentar apenas uma base, fazemos isso até chegar ao tamanho máximo de 12

TCAATGCATGCG GGTCTATATGCAT

Dai pulamos para começar comparações a partir do segundo par de base

T CAAT GCATGCGGGTCTATATGCAT

E vamos aumentando de um em um par de base. Bem isso na linguagem python fica assim.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
def complemento_reverso(s):
    complementos = {'A':'T', 'T':'A', 'G':'C', 'C':'G'}
    return ''.join([complementos[c] for c in reversed(s)])
 
def palindromo_reverso(s):
    saida = []
 
    l = len(s)
 
    for i in range(0,l-4):
        for j in range(4, 13):
            if i + j < l:
                s1 = s[i:i+j]
                s2 = complemento_reverso(s1)
 
                if s1 == s2:
                    saida.append([i + 1, j])
 
    return saida

Então criamos uma lista para salvar a posição e tamanho dos palíndromos que encontrarmos chamado saída.

Usaremos dois loops aninhados, um até o tamanho da sequência, e um do tamanho mínimo ao máximo dos palíndromos que procuramos

    for i in range(0,l-4):
        for j in range(4, 13):

Isso quer dizer um algorítimo que cresce ao quadrado com o tamanho da entrada, logo espere algo lento para uma sequência muito longa. Na verdade até que não, ja que vai depender do tamanho da sequência vezes o tamanho dos palíndromos máximos, mas vai ser lento.

dai é aquele esquema, pegamos a sequência e seu complemento

                s1 = s[i:i+j]
                s2 = complemento_reverso(s1)

E se elas forem iguais, adicionamos a resposta

                if s1 == s2:
                    saida.append([i + 1, j])

Veja que para não fazer comparações ruins, isso é, eu vou no primeiro if até o tamanho da sequência menos o menor palíndromo possível, mas dependendo do maior possível, eu posso tentar acessar posições erradas no string, por isso temos que evitar isso com esse loop

            if i + j < l:

E é isso, agora é só receber os dados, fazer as devidas arrumações

dados = open('/home/augusto/Downloads/rosalind_revp.txt').read().strip()
dados = dados[dados.find('\n')+1:].replace('\n','')

Isso é, abrir o arquivo, retirar os cabeçalhos, retirar as quebras de linhas para ficar uma string grande.

E imprimir no formato correto a saída, que é uma lista de listas

 
saida = palindromo_reverso(sequencia)
print
for i in saida:
    print i[0],i[1]

E para os dados de exemplo temos

4 6 5 4 6 6 7 4 17 4 18 4 20 6 21 4

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 e até mais.

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
def complemento_reverso(s):
    complementos = {'A':'T', 'T':'A', 'G':'C', 'C':'G'}
    return ''.join([complementos[c] for c in reversed(s)])
 
def palindromo_reverso(s):
    saida = []
 
    l = len(s)
 
    for i in range(0,l-4):
        for j in range(4, 13):
            if i + j < l:
                s1 = s[i:i+j]
                s2 = complemento_reverso(s1)
 
                if s1 == s2:
                    #print "s1:",s1," s2:",s2
                    saida.append([i + 1, j])
 
    return saida
 
##################################
## Processamento
##################################
 
sequencia = "TCAATGCATGCGGGTCTATATGCAT"
#dados = open('/home/augusto/Downloads/rosalind_revp.txt').read().strip()
#dados = dados[dados.find('\n')+1:].replace('\n','')
 
saida = palindromo_reverso(sequencia)
 
print
for i in saida:
    print i[0],i[1]

Epidemiologia – O modelo SIR com dinâmica populacional

E continuando mais um pouquinho sobre o modelo SIR.

O modelo SIR assume um tamanho populacional constante, mas lembre-se que o grupo “resistente” podia ser considerado como o grupo daqueles que conseguem a última imunidade possível, la muerte, já que quem ta morto não fica mais doente.

De toda forma, nos podemos fazer um modelo mais complexo que inclua nascimentos e mortes independentes da doença que estamos tentando entender. Então nos temos que adicionar nascimentos como um parâmetro, o b aqui (de births). Nascimentos esses que podem vir de todos os tipos de indivíduos (S + I + R), mas nos assumimos que quem acabou de nascer é suscetível apenas. E nos adicionamos uma mortalidade m para cada grupo(mS,mI,mR).

\frac{dS}{dt} = b(S+I+R) -\beta S I - mS

\frac{dI}{dt} = \beta S I - \gamma I - mI

\frac{dR}{dt} = \gamma I - mR

E tentemos ver novamente na forma de uma figura feita com o igraph, que eu achei mais facil de entender que esse monte de equação.

01

Veja que nascimentos adicionam apenas ao grupo dos suscetíveis, mas ela adiciona indivíduos baseado no tamanho populacional atual (S+I+R, não apenas o S), enquanto a mortalidade é independente da densidade subtrai indivíduos de todos os grupos.

Nos podemos assumir (apenas para começar a entender o modelo) que se indivíduos infectados e resistentes podem contribuir para prole, para os nascimentos, então essa doença deve ser relativamente benigna. Assim podemos assumir também que a mortalidade é a mesma para todos os grupos (m_i=m). Por último, vamos assumir um tamanho populacional constante. Isso significa que a taxa de natalidade é igual a de mortalidade ou b=m.

Agora imagine uma grande cidade, com digamos, um milhão de pessoas. Vamos também assumir que nos começamos com a população virtualmente toda suscetível, mas introduzimos um único individuo infectado.

1
2
3
4
N <- 10^6
R <- 0
I <- 1
S <- N - I - R

Vamos também simular a ocorrência dessa doença ao longo de uns 14 dias, Lembre-se que \gamma é o inverso da duração de uma doença.

1
g<-1/(13/365)

Dada uma tamanho populacional constante e um crescimento exponencial, então tempo de vida médio é o inverso da taxa de natalidade. Vamos considerar que o tempo médio de vida é de 50 anos.

1
b <- 1/50

Para esse modelo, a força da infeção passa a ser R_0 = 1+\frac{1}{b+\alpha}, onde \alpha é a idade média do contagio da doença. Nos podemos então estimar \beta a partir de todos os outros parâmetros, incluindo o tamanho populacional, tempo médio de vida, idade média do início da doença e o tempo médio da duração da doença. Para o exemplo, imagine que temos uma doença de crianças, onde a doença se inicia comumente la pelos 5 anos, assim temos:

1
2
idade <- 5
R0 <- 1+1/(b*idade)

então \beta tem que ser

1
B <- R0 * (g+b)/N

Finalmente, nos podemos integrar a população e seus estados.

Acontece que devido a relativa dinâmica extrema, nos queremos que a função ode tome passos pequenos, para capturar propriamente a dinâmica, nos usamos o argumento hmax para ter certeza que o maior passo é relativamente pequeno.

1
2
3
4
5
6
7
8
9
parms <-c(B=B, #Taxa de transmissão
          g=g, #Taxa de criação de resistencia
          b=b, #Natalidade
          m=b) #Mortalidade
 
anos <- seq(0,30,by=0.1)
 
library(deSolve)
SIRbd.out <- data.frame(ode(c(S=S,I=I,R=R),anos,SIRbd,parms,hmax=0.01))

E vemos o seguinte:

02

Note que a população rapidamente se torna resistente. Note também que nos temos oscilações, no entanto elas são amortecidas. Um tratamento analítico do modelo, incluindo análise de valores eigen da matriz jacobiana pode nos mostrar precisamente a periodicidade. Ela depende da idade média de início da doença e da sua duração.

Bem é isso ai, daqui a pouco começaremos a ter um pouco de posts de estatística espacial, que ja passou da hora de começar a explorar esse tópico direito e aprender essas coisas que são muito legais e úteis. 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 e até mais.

Referência:
Stevens, M. H. – 2011 – A Primer of Ecology with R. Springer

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
##################################
## Grafo
##################################
library(igraph)
 
SIR_grafo<-graph.formula(S-+I,I-+R,Entrada-+S,Saida+-S,Saida+-I,Saida+-R)
E(SIR_grafo)$label<-c(expression(paste(beta,"IS")),
                      "m",
                      expression(paste(gamma,"I")),
                      "m",
                      "m",
                      "b")
loc<-matrix(c(0,1,1,1,2,1,0,1.5,1,0.5),ncol=2,byrow=T)
 
#jpeg("01.jpg")
plot(SIR_grafo,layout=loc,rescale=F,xlim=c(0,2.2),ylim=c(0,2),vertex.size=23)
#dev.off()
 
##################################
## Modelo
##################################
SIRbd <- function(t,y,p){
    S <- y[1]
    I <- y[2]
    R <- y[3]
    with(as.list(p), {
                    dS.dt <- b * (S+I+R) - B * I * S - m * S
                    dI.dt <- B * I * S - g * I - m * I
                    dR.dt <- g * I -m * R
                    return(list(c(dS.dt,dI.dt,dR.dt)))
                })
}
 
##################################
## Parametros
##################################
N <- 10^6
R <- 0
I <- 1
S <- N - I - R
g<-1/(13/365)
b <- 1/50
idade <- 5
R0 <- 1+1/(b*idade)
B <- R0 * (g+b)/N
 
##################################
## Figura
##################################
parms <-c(B=B, #Taxa de transmissão
          g=g, #Taxa de criação de resistencia
          b=b, #Natalidade
          m=b) #Mortalidade
 
anos <- seq(0,30,by=0.1)
 
library(deSolve)
SIRbd.out <- data.frame(ode(c(S=S,I=I,R=R),anos,SIRbd,parms,hmax=0.01))
 
#jpeg("02.jpg")
matplot(SIRbd.out[,1],sqrt(SIRbd.out[,-1]),type="l",col=1,lty=1:3,
        ylab="sqrt(No. of indivíduos)",xlab="Anos",frame=F)
legend("right",c("S","I","R"),lty=1:3,bty="n")
#dev.off()

Rosalind – Finding a Protein Motif

Assim como a gente vê varias partes do DNA que tem funções especiais, por exemplo com a parte que sofre metilação, a ilhas de CG, nas proteína a gente vê esse tipo de coisa também. Então a atividade de localizar essas estruturas de interesse, esses motifs, é sempre importante.

Proteínas são registradas de forma bem parecida com DNA, mas ao invés de 4 letras, 4 bases que formam o DNA, proteínas são formadas por 20 aminoácidos, que são codificados por aquelas trincas de DNA, os codons, na verdade o Dogma Central da Biologia Molecular é DNA -> RNA -> Proteínas, não perfeitamente isso, tem mais setinhas ai no meio, mas essa é a ideia.

Mas beleza, a ideia aqui então é que temos uma estrutura, um motif de proteína, que pode estar representado de várias formas, mas tem uma estrutura geral, a presentação é assim, [XY] significa “X ou Y” e {X} significa “qualquer aminoácido exceto X”. Por exemplo, o motif N-glycosylation está escrito como N{P}[ST]{P}, então é um N seguido de qualquer coisa menos P, depois um S ou T e depois qualquer coisa menos P denovo.

A tarefa é receber uma lista de nomes de proteínas, o id que estão depositados num banco de dados na verdade, e encontrar onde começa o N-glycosylation nelas, caso ocorra, em qualquer uma de suas formas.

O exemplo que temos é

A2Z669 B5ZC00 P07204_TRBM_HUMAN P20840_SAG1_YEAST

Então o primeiro passo é conseguir baixar os dados, aqui precisaremos da internet.

Mas basicamente é so pegar o link do banco de dados

http://www.uniprot.org/uniprot

e adicionar o id do deposito

http://www.uniprot.org/uniprot/A2Z669

Cole isso no navegador, que você vera a página do A2Z669. Então aqui acessar o banco de dados é simplesmente isso, não envolve query nem nada complicado. Então da pra usar o pacote requests do python para fazer o pedido http e receber o resultado.

Então é so montar o request e ele retorna um objeto da classe dele (), mais ou menos assim:

entradas =  arquivo.read().strip().split('\n')
print entradas
for entrada in entradas:
    r = requests.get('http://www.uniprot.org/uniprot/%s.fasta' % entrada)

Veja as entradas, e veja que depois, para cada uma estamos fazendo um request, e unindo a entrada no arquivo.
Outra coisa, não estamos pedindo a página completa, como la em cima, estamos pedindo a pelo nome do id.fasta, entre o seguinte link no navegador e veja o que é retornado.

http://www.uniprot.org/uniprot/A2Z669.fasta

Apenas texto, no formato fasta, que é mil vezes mais fácil de processar que a página html. A gente já mexeu nisso a muito tempo atrás aqui, usando o R.

Ok, feito isso, um dos atributos da classe request é o text, que é o texto que veio no pedido http.

No primeiro id, o texto é assim:

>sp|A2Z669|CSPLT_ORYSI CASP-like protein 5A2 OS=Oryza sativa subsp. indica GN=OsI_33147 PE=3 SV=1 MRASRPVVHPVEAPPPAALAVAAAAVAVEAGVGAGGGAAAHGGENAQPRGVRMKDPPGAP GTPGGLGLRLVQAFFAAAALAVMASTDDFPSVSAFCYLVAAAILQCLWSLSLAVVDIYAL LVKRSLRNPQAVCIFTIGDGITGTLTLGAACASAGITVLIGNDLNICANNHCASFETATA MAFISWFALAPSCVLNFWSMASR

Então a gente acha a primeira quebra de linha, e só usa o que está depois dela, que são é a proteína propriamente dita

enter = proteina.find('\n')
    proteina = proteina[(enter+1):]

Tira todos as quebras de linha, os ‘\n’

proteina = proteina.replace('\n','')

Depois usamos expressão regular para achar o motif, que a representação é basicamente uma expressão regular pronta, é até fácil converter

    busca = re.finditer('(?=(N[^P][ST][^P]))',proteina)
    saida=[]
    for i in busca:
        saida.append(i.start()+1)

Aqui só tem um pulo do gato que eu apanhei um monte para dar certo, a gente tem que achar todas as ocorrências, mas se você usar (N[^P][ST][^P]) como expressão regular apenas, ela vai achar uma ocorrência e depois vai começar a procurar depois do fim dessa ocorrência, assim você não vai achar todas as ocorrências, porque por exemplo, uma sequência XXXXXNNSTNXXXXX, tem duas ocorrências, NNST e NSTN, so que com aquela expressão regular a gente so pega uma, por isso da interrogação ali no começo, a interrogação vai fazer a gente voltar a procurar imediatamente após o início do primeiro motif encontrado. É só entrar no fórum de duvidas, que você vai ver um monte de gente, assim como eu diga-se de passagem, apanhando exatamente na mesma coisa.

Ok, então tudo certo aqui, achamos os motifs nos fastas de cada id

A2Z669 B5ZC00 NFSD NSSN NWTE NLSK NISA P07204_TRBM_HUMAN NASQ NNTS NTSY NQTS NQTA P20840_SAG1_YEAST NSSQ NTTF NFSD NDTN NTTY NTSA NRTT NITV NITN NATR NFTS

Agora é só formatar a resposta de forma correta, esse probleminha foi bem legal. 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 e até mais.

import re
import requests
 
#arquivo = open('/home/augusto/Downloads/rosalind_mprt.txt','r')
arquivo = open('mprt.in','r')
entradas =  arquivo.read().strip().split('\n')
 
#print entradas
 
for entrada in entradas:
    r = requests.get('http://www.uniprot.org/uniprot/%s.fasta' % entrada)
    proteina = r.text
    enter = proteina.find('\n')
    proteina = proteina[(enter+1):]
    proteina = proteina.replace('\n','')
    busca = re.finditer('(?=(N[^P][ST][^P]))',proteina)
    saida=[]
    for i in busca:
        saida.append(i.start()+1)
        #print i.start()+1,
        #print proteina[i.start():i.start()+4]
    if len(saida)>0:
        print entrada
        for i in saida:
            print i,
        print