Propriedades em rede de interações, como descrever as comunidades de plantas e seus polinizadores?

Eu ja falei uma vez sobre como fazer gráficos de rede aqui, mas podemos explorar também várias propriedades de redes interessantes para a biologia.

Primeiro vamos simular aqui quatro redes para ajudar a pensar, como eu não sei muito bem como começar, vamos usar o modelo de Erdos para a simulação da rede de interações. Alias Paul Erdos era um cara muito interessante, diria até fascinante quanto ao seus estilo de vida, vale a pena uma olhadinha na história dele.

Mas por enquanto vamos olhar as redes que criamos onde bolinhas verdades são espécies de plantas e quadradinhos azuis são espécies de polinizadores.

01

Bem, eu fiz 4 redes, todas tem 15 polinizadores (linhas) e 15 espécies de plantas (colunas). As plantas dependem dos polinizadores para reprodução, e os polinizadores também dependem das flores para viver (de forma geral, existem excessões). Isso nos leva a conclusão óbvia que todo mundo que a gente ta vendo ai tem pelo menos 1 interação, cada polinizador visita pelo menos um espécie de flor.

Mas além disso podemos pensar em alguma possibilidades pensando nos polinizadores. Uma possibilidade é o polinizador ser generalista, visitar várias flores. Isso é bom já que ele deve encontrar flores mais facilmente, porém ele tem que estar preparado para ser um “pau para toda obra” já que cada flor tem uma forma (morfologia), as vezes disponibiliza seus recursos de uma forma diferente umas das outras, o que pode tornar esse polinizador meia sola em muitas espécies de plantas.

Outra estratégia que pode ser interessante se pensarmos bem é ser especializado. Ja que melhor que pegar só um pouquinho de comida de porta em porta que se bate, é ganhar um montão de comida em uma porta. Mas diferente do generalista, o especialista só consegue o que quer de uma única espécie de planta.

Em termo de energia ou tanto de comida que você consegue coletar, as duas estratégias podem ser equivalente, já que passar bastante tempo voando para pegar bastante pólen em uma única flor ou voar menos e ir pegando um pouquinho de pólen em cada flor pode dar na mesma no final das contas, mas existe outra coisa a se levar em conta, um polinizador generalista pode chegar em novos ambientes facilmente, já que ele consegue comida em varias espécies, enquanto o especialista só vai chegar no ambiente depois que sua “fonte de comida” chegou.

Até agora só consideramos comida mas podem ter flores outras utilidades para os polinizadores e existem outros tipos de interações. Existem ainda os interesses das plantas em todo esse jogo relações.

Mas talvez antes de ficar pensando em possibilidades e mais possibilidades, casos e mais casos, vamos olhar como é essa rede de interações como um todo.

Vamos olhar primeiro o exemplo 1 da perspectiva dos polinizadores.
Nele a gente vê bastante espécies de polinizadores especialista, que só visitam uma espécie de planta, como 2 casos de exclusividade polinizador e planta no meio do gráfico, 3 casos de grupinho envolvendo 2 polinizadores e 2 plantas e um grupão maior no canto inferior direito mais ou menos.

Mas a gente pode contar quantas plantas cada polinizador visita e fazer o seguinte gráfico;

02

Então vamos la, no eixo x temos o número de interações e no eixo y temos a contagem de polinizadores. Ou seja exitem 15 polinizadores, e como todo mundo interage com pelo menos uma planta (o obvio, senão nem na matriz de interação eles estariam), nos temos que com 1 planta os 15 polinizadores interagem, desses 15 polinizadores, apenas 7 interagem com 2 espécies de plantas e um apenas interagiu com 3 plantas.

Certo mas esse gráfico vai ser dependente do número de polinizadores no sistema, mas podemos resolver esse problema transformando tudo em porcentagem, ou seja dividir pelo número de polinizadores total. Então ficamos com…

03

Exatamente o mesmo gráfico, mas notem que agora o eixo X vai de 0 a 1, ou seja é independente do número de polinizadores, isso facilita se quisermos comparar características que vemos nesse gráfico com outras comunidades, independente do número de espécies de polinizadores. Mas ao invéz de barrinhas de contagem, podemos simplesmente usar bolinhas e ficamos com isso.

04

Agora sim, bolinha é bom para ajustar modelos, mas na verdade ainda estamos olhando a mesma coisa do primeiro gráfico de barras, apenas alguns ajustes aqui e ali.

Nesse caso um dos modelos que ajuda a descrever essa comunidade segundo o Jordano (2003) é o exponencial, é um modelo com dois parâmetros que ele chama de b e gamma. Gamma é a letra grega  gamma

E o modelo é assim:

 p(k) = b cdot exp(-gamma cdot k)

Mas notem que no artigo o Jordano (2003) não cita o b, que não fala muita coisa, mas ele não usa o símbolo de igual, ele usa o símbolo de aproximadamente, acho que para não citar o b, não entendo muito o porque disso, mas vamos em frente, é vivendo e aprendendo.

Bem podemos ajustar esse modelo não linear no R usando a função nls(). E ajustando temos…

05

Agora que porra é essa? Bem esse modelo (a linha vermelha) ta nos dizendo como muda a distribuição acumulada de interações, ou seja, numa comunidade onde todo mundo é especialista como essa do exemplo 1, as barrinhas dos gráficos la em cima descem muito rápido, formam uma escada muito íngreme, de degraus altos, enquanto tivessem mais generalistas, a escadinha seria bem menos íngreme e subir essa escadinha seria fácil. E é isso que o gamma do modelo ta falando.
Podemos fazer o mesmo procedimento para os 4 exemplos de comunidade e ver o gamma neles.

coef(m1.exe1) gamma b 0.9372338 2.5877177 coef(m1.exe2) gamma b 0.5020859 1.6452600 coef(m1.exe3) gamma b 0.2602011 1.3744055 coef(m1.exe4) gamma b 0.1417438 1.2410299

E sem olhar para o b por enquanto, mas somente o gamma que nos interessa agora, sabemos que a comunidade 1 é a que tem mais especialistas e menos generalistas, enquanto a comunidade do exemplo 4 tinha muitos mais generalistas e poucos especialistas, olhem o emaranhado que são as interações no exemplo 4.

Isso se reflete nos valores de gamma, olhem para o exemplo 4, ele é bem próximo de zero, bem baixinho 0.14, mas o exemplo 1, cheio de especialistas e quase nenhum generalista forte temos um valor de gamma de 0.93, quase 1.

Ou seja valor alto de gamma é uma escada íngreme naquele gráfico de barras, sinal de especialistas na comunidade, muitos especialistas.
Valor baixo de gamma é sinal de mais generalistas.
Agora temos um pouco de intuição sobre o que o modelo esta ajustando não?

Mas o Jordano citou três modelos no artigo dele, ajustar o modelo desse jeito vai falhar em um caso que imagino. No caso de todo mundo ser generalista. Isso eu imagino pois a gente pode tentar olhar o que acontece com as curvas de acordo com a mudança de gamma.
Como olhando a formula não me traz muitos insights pois sou péssimo de matemática, eu resolvo ela para tudo que é valor de gamma e fico olhando o que acontece.

06

E olha ai, vemos como quanto maior o valor de gamma, mais íngreme é a subida, e biologicamente isso é mais gente com 1 interação apenas, ou seja, especialistas dominam na comunidade. Mas o gamma vai ser no máximo 1, quando ele for 1 o modelo vai ser uma reta inclinada, e no caso de somente haver generalistas, precisaríamos de uma reta não inclinada, ou seja esse modelo ia falhar. Mas o Jordano foi legal deixou outras opções. Mas investigaremos elas depois 🙂

Tudo muito legal mas foram um milhão de comandos, eu não sei fazer tudo isso denovo. Tire a preocupação do seu coração que o senhor Carsten foi muito gentil em resumir todo esse processo numa função do pacote bipartite. La é só digitar degreedistr(sua comunidade) que toda essa magia se realiza para os polinizadores e para as plantas de uma vez so, apesar que não abordamos plantas aqui mas a idéia continua a mesma.

O Resultado vai sair assim e com um gráfico que podemos comparar com o que fizemos até aqui, mas lembrando que aqui estamos olhando para o exemplo 4, a função do bipartite não vai te deixar ajustar o modelo se vc tiver somente 3 barrinhas como nos temos no exemplo 1, o que é fato ja que são poucos graus de liberdade, mas que dificultaria nossa didática aqui.

07

degreedistr(exe4, plot.it=T) $`lower trophic level dd fits` Estimate Std. Error Pr(>|t|) R2 AIC exponential 0.1417439 0.01486731 1.210787e-05 0.9717753 -18.632825 power law 0.4722450 0.09654721 1.206915e-03 0.8745991 -4.366322 truncated power law -0.3072618 0.12245915 4.045036e-02 0.9856348 -23.467338 $`higher trophic level dd fits` Estimate Std. Error Pr(>|t|) R2 AIC exponential 0.1377984 0.01249272 1.573018e-06 0.9746622 -23.366621 power law 0.4814095 0.08664691 3.536921e-04 0.8791545 -7.092832 truncated power law -0.2837947 0.11314399 3.646958e-02 0.9863076 -28.219390

E lembrando que polinizadores que eram nossas linhas são o “lower trophic level” nesse caso, e olhem que o modelo exponential tem um estimate de 0.14…

coef(m1.exe4) gamma b 0.1417438 1.2410299

E nosso estimate é o mesmo valor, 0.14, ohhhhh, será que porque eu copiei o procedimento dele? 🙂

Bem é isso, fica no script todo o procedimento que fiz para gerar as figuras e a estimativa, depois eu preparo mais posts com os outros modelos, principalmente que essa propriedade tem outro grande interesse em biologia, se um gamma maior significa mais especialistas, um gamma maior também traz algo a respeito da susceptibilidade da comunidade a extinção em cascata.

#Semente para os resultados serem reproduziveis
set.seed(123)
###################################################
#Função para gerar rede de interação baseado no modelo de Erdős–Rényi
#Creditos a Edwin Grappin, vi o post no R-bloggers
#http://probaperception.blogspot.com.br/2012/12/have-you-tried-to-understand-your.html
###################################################
generateER = function(n = 100, p = 0.5){
map = diag(rep(1, n))
link = rbinom(n*(n-1)/2, 1,p)
t = 1
for(j in 2:n){
for(i in 1:(j-1)){
map[i,j] = link[t]
t = t + 1
}
}
dimnames(map)<-list(paste("Pol",1:n),paste("Sp",1:n))
return(map)
}
###########################
#Minha função para transformar matriz em lista para o igraph
##############################
mat2lis<-function(matriz) {
lista<-dimnames(matriz)
dados<-data.frame(Pol=rep(lista[[1]],length(lista[[1]])),
SP=rep(lista[[2]],each=length(lista[[2]])),
int=c(matriz))
return(dados[dados$int==1,1:2])
}
############################
library(igraph)
 
#Gerando exemplos
exe1<-generateER(n=15,p=0.1)
exe2<-generateER(n=15,p=0.2)
exe3<-generateER(n=15,p=0.4)
exe4<-generateER(n=15,p=0.8)
 
#Transformando na forma de dados que o igraph usa
exe1.l<-mat2lis(exe1)
exe1.g <-graph.data.frame(exe1.l, directed=F)
 
exe2.l<-mat2lis(exe2)
exe2.g <-graph.data.frame(exe2.l, directed=F)
 
exe3.l<-mat2lis(exe3)
exe3.g <-graph.data.frame(exe2.l, directed=F)
 
exe4.l<-mat2lis(exe4)
exe4.g <-graph.data.frame(exe4.l, directed=F)
 
#Plots das redes
par(mfrow=c(2,2))
 
plot(exe1.g,layout=layout.fruchterman.reingold,
vertex.color=rep(c("blue","green"),each=15),
vertex.label=c(paste("P",1:15),paste("F",1:15)),
vertex.shape=rep(c("square","circle"),each=15),
edge.color="black",edge.width=3,vertex.label.cex=0.7,
vertex.label.color="black"
)
title("Exemplo 1")
 
plot(exe2.g,layout=layout.fruchterman.reingold,
vertex.color=rep(c("blue","green"),each=15),
vertex.label=c(paste("P",1:15),paste("F",1:15)),
vertex.shape=rep(c("square","circle"),each=15),
edge.color="black",edge.width=3,vertex.label.cex=0.7,
vertex.label.color="black"
)
title("Exemplo 2")
 
plot(exe3.g,layout=layout.fruchterman.reingold,
vertex.color=rep(c("blue","green"),each=15),
vertex.label=c(paste("P",1:15),paste("F",1:15)),
vertex.shape=rep(c("square","circle"),each=15),
edge.color="black",edge.width=3,vertex.label.cex=0.7,
vertex.label.color="black"
)
title("Exemplo 3")
 
plot(exe4.g,layout=layout.fruchterman.reingold,
vertex.color=rep(c("blue","green"),each=15),
vertex.label=c(paste("P",1:15),paste("F",1:15)),
vertex.shape=rep(c("square","circle"),each=15),
edge.color="black",edge.width=3,vertex.label.cex=0.5,
vertex.label.color="black"
)
title("Exemplo 4")
 
#Ajustando o modelo exponencial para polinizadores
#Jordano, P., Bascompte, J. and Olesen, J. M. (2003)
#Invariant properties in coevolutionary networks of
#plant-animal interactions. Ecology Letters 6, 69–81
 
#Rede do Exemplo 1
#Número de interações
k<-sum(exe1)
#Número de interações por espécie de polinizador
pol.n<-rowSums(exe1)
#Número de polinizadores por quantidade de interação
P.pol<-sapply(sort(unique(pol.n)), function(x) sum(pol.n>=x))
#Juntando os dados
dados.pol1 <- cbind.data.frame(k = sort(unique(pol.n)), P = P.pol/max(P.pol))
#Ajustando modelo exponencial aos dados
m1.exe1<-nls(P~b*exp(-gamma*k),start=list(gamma = 0, b = 1),data=dados.pol1)
 
#Grafico 2
barplot(P.pol,names.arg=c(1,2,3),
xlab="Número de espécies de flores que visita (K)",
ylab="Contagem de polinizadores")
 
#Grafico 3
barplot(dados.pol1[,2],names.arg=c(1,2,3),
xlab="Número de espécies de flores que visita (K)",
ylab="Contagem de polinizadores")
 
#Grafico 4
plot(dados.pol1[,2]~dados.pol1[,1],pch=19,cex=2,ylab="Distribuição Acumulada",
xlab="K",frame.plot=F,ylim=c(0,1),xlim=c(1,5))
abline(h=1,lty=2)
 
#Grafico 5
plot(dados.pol1[,2]~dados.pol1[,1],pch=19,cex=2,ylab="Distribuição Acumulada",
xlab="K",frame.plot=F,ylim=c(0,1),xlim=c(1,5))
abline(h=1,lty=2)
curve(coef(m1.exe1)[2]*exp(-coef(m1.exe1)[1]*x),1,5,add=T,lwd=2,lty=2,
cex=1.5,col="red")
 
#Exemplo 2
k<-sum(exe2)
S<-sum(dim(exe2))
pol.n<-rowSums(exe2)
P.pol<-sapply(sort(unique(pol.n)), function(x) sum(pol.n>=x))
dados.pol2 <- cbind.data.frame(k = sort(unique(pol.n)), P = P.pol/max(P.pol))
m1.exe2<-nls(P~b*exp(-gamma*k),start=list(gamma = 0, b = 1),data=dados.pol2)
 
#Exemplo 3
k<-sum(exe3)
pol.n<-rowSums(exe3)
P.pol<-sapply(sort(unique(pol.n)), function(x) sum(pol.n>=x))
dados.pol3 <- cbind.data.frame(k = sort(unique(pol.n)), P = P.pol/max(P.pol))
m1.exe3<-nls(P~b*exp(-gamma*k),start=list(gamma = 0, b = 1),data=dados.pol3)
 
#Exemplo 4
k<-sum(exe4)
pol.n<-rowSums(exe4)
P.pol<-sapply(sort(unique(pol.n)), function(x) sum(pol.n>=x))
dados.pol4 <- cbind.data.frame(k = sort(unique(pol.n)), P = P.pol/max(P.pol))
m1.exe4<-nls(P~b*exp(-gamma*k),start=list(gamma = 0, b = 1),data=dados.pol4)
 
#Extraindo coeficientes
coef(m1.exe1)
coef(m1.exe2)
coef(m1.exe3)
coef(m1.exe4)
 
#Niveis de Gamma
graus<-seq(0.1,0.9,0.1)
 
#GRafico 6
jpeg("06.jpg")
plot(0,0,type="n",xlim=c(0,10),ylim=c(0,1),frame.plot=F,
xlab="Número de interações",
ylab="Distribuição Cumulativa do número de interações")
for(i in 1:length(graus)) {
curve(1*exp(-graus[i]*x),0,10,add=T,col=i,lty=2,lwd=2)
}
legend("topright",lty=2,lwd=2,col=c(1:length(graus)),legend=paste("Gamma é ",graus))
dev.off()
 
#Grafico e ajuste usando o bipartite para o exemplo 4
library(bipartite)
degreedistr(exe4, plot.it=T)

Modelo de crescimento populacional

A famosa equação de crescimento populacional, ou crescimento densidade dependente, ou crescimento logístico.

  {dN\over dt} = r\cdot N\cdot (1-{N\over K})

Ela foi descrita a primeira vez por Pierre François Verhulst, mas depois redescrita pelo Lotka, que o nome normalmente está associada a ela. Mas o interessante é que o Lotka foi quem começou a expandir ela, como no modelo de predador presa dos coelhos e dos linces.

Mas antes gente devia pensar, hoje, porque ainda vale a pena perder tempo olhando para essa equação? se a gente pode fazer modelos infinitamente mais complexos.

modelo

O problema que apesar da gente ser sim capaz de por exemplo usar o computador para produzir modelos melhores com muitas constantes e uma excelente capacidade de predição, um r² alto, eles serão difíceis de nos fornecer insights, de pensar o que esta acontecendo o que vai acontecer justamente porque são complexos. O modelo de crescimento populacional pode ser “ruim”, ter uma capacidade de previsão pior, mas rapidamente a gente consegue insights de como o mundo funciona olhando para ele.

Então vamos olhar melhor para essa equação diferencial, alias chama diferencial que ela fala a diferença de uma quantidade (N, tamanho populacional) dada uma diferença de tempo (t, tempo).
Além desses dois itens temos o r que é a taxa de crescimento.

Mas olhando ela a gente pode separar ela em duas partes.

Essa:

  r\cdot N

Que nos diz o que? A população cresce numa taxa R, e quanto maior o N, mais ela cresce, como o N cresce constantemente o crescimento é exponencial, considerando so esse pedacinho. Mas ele multiplica outra parte que é…

  (1-{N\over K})

Agora olha ai, é 1- N dividido por k, k é um numero fixo, a capacidade suporte e N é o tamanho da população. Aqui a gente pode pensar em termos de 3 momentos.
Primeiro quando a população N for bem pequenininha, vamos ter um número pequenininho divido por k, um numero grande, logo algo pequeno dividido por algo grande vai ser um númerozinho bemmmm pequeno, e 1 menos um numero pequeno, quase nada, da praticamente 1, logo o crescimento é rápido, praticamente exponencial.

Com o tempo o N vai aumentando, até o N ficar grande, quando o N ficar grande, vai ser um numero grande dividido por um numero grande, que vai dar 0, logo toda a equação vai ser multiplicada por 0 e o crescimento vai dar zero, assim a população se mantém estável no valor de k.

E por ultimo e muito interessante, é que se N for um numero maior que o K, o que vai acontecer?, a divisão vai ser um numero negativo maior que 1, logo 1 menos um numero negativo maior que 1 vai dar um numero negativo, e o crescimento vai multiplicar um número negativo, então vai ser negativo, e isso até ele voltar a população voltar a k.

Magnifico não?

podemos fazer uma função simples no R para calcular o crescimento logístico e ficar brincando com ela para ver se essas predições ocorrem. Alias isso é bem legal, ficar mudando os valores para ver se todas essas coisas funcionam assim mesmo.

#Função
dlogistic <- function(K= 100, rd = 1, N0 = 2, t = 15) {
N <- c(N0, numeric(t))
for (i in 1:t) N[i + 1] <- {
N[i] + rd * N[i] * (1 -N[i]/K)
}
return(N)
}
 
#Gerando a população com os argumentos padrão
Nts <- dlogistic()
 
#Grafico
t <- 15
k <- 100
plot(0:t, Nts,xlab="Tempo",ylab="Tamanho da População")
abline(h = k, lty = 3)
text(5,k-5,"Capacidade Suporte")

E temos nosso gráfico

01

Outra coisa ainda é pensar como diabos as pessoas pensam nessas formulas, o professor Gotteli demonstra como da para chegar na equação de crescimento logístico das taxas de natalidade e mortalidade.

Começamos que o crescimento é a natalidade menos a mortalidade

 {dN\over dt} = (b\cdot N-d\cdot N)*N

A taxa de natalidade de um momento vai ser o que nasceu em um ano menos o que nasceu no outro, lembrando que a natalidade não é constante, começa a falta alimento, espaço ela diminui, ou seja como esta mudando a natalidade, a mortalidade vai ser o mesmo esquema

 bN = b_{0}-b_{1}\cdot N

Substituimos isso para a natalidade e a mortalidade

 {dN\over dt} = ((b_{0}-b_{1}\cdot N) - (d_{0}-d_{1}\cdot N))\cdot N

Com algum algebrismo mexemos agora as coisa para ca, para la. Se você como eu é pessimo de algebra, pode estudar junto comigo aqui 🙂

 {dN\over dt} = (b_{0}-d_{0})\cdot N\cdot (1-{(b_{1}-d_{1})\over (b_{0}-d_{0})}\cdot N)

E no final temos a equação de crescimento logístico, igualzinho la em cima, so trocar esses termos por r e k.

  {dN\over dt} = r\cdot N\cdot (1-{N\over K})

Bem, vamos tentar falar um pouco mais disso, principalmente porque dessa equação que a gente chega em metapopulações, que é um conjunto de populações, ou se preferir um conjunto dessas equações heheh.

Referencias:

Soetaert, K, Herman, P M. J. – 2008 – A Practical Guide to Ecological Modelling Springer

Stevens, M. H. – 2011 – A Primer of Ecology with R Srpinger

Rosalind problema 10 – Finding a Shared Motif LCS

E la vamos nos em mais um problema de informática. Dado um n de sequências, achar a maior sub-sequência comum a todas as sequências.

Primeiro eu fui dar uma olhada no google e a primeira parada foi essa pagina da wikipedia, parece que esse problema não é especifico a bioinformática. Ele ensina um esquema muito louco usando matriz. Mas após pensar muito eu imaginei três caminhos a seguir. O esquema com matriz do wikipedia, o esquema com árvore que também é comentado no wikipedia e a boa e velha força bruta.

Mas vamos la.

Então o exemplo do problema fornece 3 sequências:

GATTACA TAGACCA ATACA

Aqui qual a maior sub-sequência comum a todo mundo? Nesse caso a maior so tem 2 pares de base. Mas existe mais de uma.

GAT--TA--CA TA--GACCA A--TA--CA

Todo mundo tem TA, mas outras sequências são comuns a todo mundo também como:

GATTA--CA TAGAC--CA CA ATA--CA GATT--AC--A TAG--AC--CA AC AT--AC--A

Certo, esse problema é fácil de entender, mas solucionar é meio complicado para mim pelo menos, o primeiro esquema com matriz é assim.

Vamos supor outras string so para ver uma sequência comum maior.

"ACACACACACGTAC" "GTGTACACACACAAA"

Então primeiro comparando so essas 2 string a gente faz o seguinte, usa uma como nomes das linhas e outra como nome das colunas de uma matriz.

Dai vamos comparar a primera letra da sequência 2 com a sequência 1, e vemos que o G da sequência 2 ocorre na 11 base da sequência 1

A C A C A C A C A C G T A C G 0 0 0 0 0 0 0 0 0 0 1 0 0 0

Dai vamos para a segunda linha.

A C A C A C A C A C G T A C G 0 0 0 0 0 0 0 0 0 0 1 0 0 0 T 0 0 0 0 0 0 0 0 0 0 0 2 0 0

Agora um T da sequência 2, vemos que na 12 base ele ocorreu, mas opa na 11 ja tinha um 1, so olhar a diagonal ali, então eu preencho com um 2 ja que é a segunda ocorrência em linha. Ai preenchendo toda a matriz eu tenho:

A C A C A C A C A C G T A C G 0 0 0 0 0 0 0 0 0 0 1 0 0 0 T 0 0 0 0 0 0 0 0 0 0 0 2 0 0 G 0 0 0 0 0 0 0 0 0 0 1 0 0 0 T 0 0 0 0 0 0 0 0 0 0 0 2 0 0 A 1 0 1 0 1 0 1 0 1 0 0 0 3 0 C 0 2 0 2 0 2 0 2 0 2 0 0 0 4 A 1 0 3 0 3 0 3 0 3 0 0 0 1 0 C 0 2 0 4 0 4 0 4 0 4 0 0 0 2 A 1 0 3 0 5 0 5 0 5 0 0 0 1 0 C 0 2 0 4 0 6 0 6 0 6 0 0 0 2 A 1 0 3 0 5 0 7 0 7 0 0 0 1 0 C 0 2 0 4 0 6 0 8 0 8 0 0 0 2 A 1 0 3 0 5 0 7 0 9 0 0 0 1 0 A 1 0 1 0 1 0 1 0 1 0 0 0 1 0 A 1 0 1 0 1 0 1 0 1 0 0 0 1 0

Dai quanto eu tirar o maior numero dessa matriz, ele será a maior sequência comum, e a sequência será a diagonal, ou o nome das colunas da diagonal que ele ocorre. Mas como podem ocorrer várias sub-sequências, eu pego e retiro todas as subsequências por ordem em uma lista.

lista.cs [[1]] [1] "AC" "AC" "AC" "AC" "AC" "AC" "AC" "AC" "GT" "GT" "AC" "AC" "AC" [[2]] [1] "ACA" "ACA" "ACA" "ACA" "ACA" "ACA" "ACA" "GTA" [[3]] [1] "ACAC" "ACAC" "ACAC" "ACAC" "ACAC" "ACAC" "GTAC" [[4]] [1] "ACACA" "ACACA" "ACACA" "ACACA" "ACACA" [[5]] [1] "ACACAC" "ACACAC" "ACACAC" "ACACAC" [[6]] [1] "ACACACA" "ACACACA" "ACACACA" [[7]] [1] "ACACACAC" "ACACACAC" [[8]] [1] "ACACACACA"

E pimba, eu tenho de baixo para cima as maiores sequências comuns de acordo com o número de bases. Minha ideia original era dai usar essa lista, e procurar ela na outras sequências, sistematicamente, da maior para a menor. Mas como eu usei muito loop de for no R, sem vetorizar as operações, meu script ficou muitooooo lento. Para desenhar essa matriz com muitas bases demorava muito. O código funcionava perfeito, mas mais de uma 500 bases para cima ficava ruim. Eu pensei em algumas coisas para melhorar mais ai mudei de idéia.

Fui para a força bruta. O que fiz foi o seguinte, imaginemos uma sequência. A primeira do exemplo la em cima, se a gente aninhar 2 loops de for, da para retirar todas as sub-sequências possíveis.

> for(i in 1:nchar(sequencia)){ + for(j in 1:nchar(data[1])) { + print(substr(data[1],i,i+j)) + } + } [1] "GA" [1] "GAT" [1] "GATT" [1] "GATTA" [1] "GATTAC" [1] "GATTACA" [1] "GATTACA" … [1] "CA"

A partir daqui é só usar uma função ajudante para testar uma a uma dessas sequências se são sub-string das outras, no R com expressão regular, usando o grepl para um resultado de TRUE or FALSE, a partir dai é só fazer a sequência, a gente começa com um sub-string com 0 carácter, e sempre que algo for sub-string de todo mundo, e tiver mais caracteres que o sub-string atual, a gente sobrescreve ele, senão não, assim no final a gente vai ter um dos maiores sub-strings possíveis que todo mundo compartilha. Isso ficou um pouco mais rápido devido a minha incompetência para calcular a matriz ali em cima, mas no final fica a mesma coisa.

Então eu usei essa função no R, a is.substr é auxiliar, ela so testa se é sub-string de todo mundo.

#Função Auxiliar
is_substr<-function(find,data) {
if(length(data)<1 & nchar(find)<1) { T } else {
resp<-logical()
for(i in 1:length(data)) {
if(grepl(find,data[i])) {resp<-TRUE} else {resp<-FALSE;break}
}
resp
}
}
 
#Exemplo de como vai estar os dados
data<-c("GATTACA","TAGACCA","ATACA")
 
#Função para achar o menor sub-string comum
long_substr<-function(data) {
sub.m<-""
if(length(data) > 1 & nchar(data[1]) > 0) {
for(i in 1:nchar(data[1])){
for(j in 1:nchar(data[1])) {
if(is_substr(substr(data[1],i,i+j),data) &
nchar(sub.m)<nchar(substr(data[1],i,i+j))) {
sub.m<-substr(data[1],i,i+j)
}
}
}
}
sub.m
}

E ainda sim ficou meio lento, meio eu to sendo gentil, muito lento na verdade. Mas foi preguiça também, pois é possível comparar tudo de forma mais eficiente se compararmos do maior substring para o menor, e parar no primeiro que encontrar dai, mas ja estava triste com o primeiro código e parei.

E em python o código ficou idêntico ao do R. Que na verdade eu tinha visto esse código aqui em python, mas ele é bem simples, tanto que ficou a mesma coisa em R.

def is_substr(find, data):
    if len(data) < 1 and len(find) < 1:
        return False
    for i in range(len(data)):
        if find not in data[i]:
            return False
    return True
 
def long_substr(data):
    substr = ''
    if len(data) < 1 and len(data[0]) < 0:
        for i in range(len(data[0])):
            for j in range(len(data[0])-i+1):
                if j < len(substr) and is_substr(data[0][i:i+j], data):
                    substr = data[0][i:i+j]
    return substr
 
f = open('/home/augusto/Downloads/rosalind_lcs.txt', 'r')
lista = f.read().splitlines()
f.close()
 
print long_substr(lista)

So olhando que agora eu aprendi a abrir arquivos em python, pelo menos para tirar só as sequências. Eu abro com open, dai leio com read e separo em uma lista com splitlines e depois fecho o arquivo, que ele para de ocupar memoria. Mas a solução em python é bem mais rapida, ja que nele não precisa vectoriza as coisas.

E segue o código usado, esse problema foi complicado, principalmente por conta de fazer um método eficiente. No fórum de soluções as pessoas colocam o tempo que seus scripts levam, mas pelo menos deu para encontrar uma solução. Otimização de código ainda não esta em pauta por enquanto, segue todo o código que usei caso alguém queria dar uma olhada.

########
#Função para checar se algo é comum a todos os strings de um vetor
is_substr<-function(find,data) {
if(length(data)<1 & nchar(find)<1) { T } else {
resp<-logical()
for(i in 1:length(data)) {
if(grepl(find,data[i])) {resp<-TRUE} else {resp<-FALSE;break}
}
resp
}
}
 
#testando
is_substr("A",c("GATTACA","TAGACCA","ATACA"))
is_substr("G",c("GATTACA","TAGACCA","ATACA"))
 
#Agora fornecendo varios strings em um vetor
data<-c("GATTACA","TAGACCA","ATACA")
 
#função para encontrar o maior sub-string dos strings de um vetor
long_substr<-function(data) {
sub.m<-""
if(length(data) > 1 & nchar(data[1]) > 0) {
for(i in 1:nchar(data[1])){
for(j in 1:nchar(data[1])) {
if(is_substr(substr(data[1],i,i+j),data) &
nchar(sub.m)<nchar(substr(data[1],i,i+j))) {
sub.m<-substr(data[1],i,i+j)
}
}
}
}
sub.m
}
 
#Testando, mas fica bem lento conforme se aumenta
#o número e tamanho dos strings
long_substr(data)
 
#Abrindo os dados do Rosalind e transformando num Vetor
strings<-readLines("/home/augusto/Downloads/rosalind_lcs.txt")
strings<-as.vector(strings)
 
#Usando a função
long_substr(strings)
 
#Exemplo com matriz do wikipedia
string1<-"ACACACACACGTAC"
string2<-"GTGTACACACACAAA"
string3<-"ATACA"
 
#picando os 2 primeiros strings
str1<-strsplit(string1,"")[[1]]
str2<-strsplit(string2,"")[[1]]
 
#olhando o tamanho de cada um
m<-length(str1)
n<-length(str2)
 
#fazendo uma matriz em que os strings são os nomes da linhas e colunas
matriz<-matrix(0,ncol=m,nrow=n,dimnames=list(str2,str1))
matriz
 
#Preenchendo com 1 onde ocorre a mesma letra
for(i in 1:m) {
for(j in 1:n) {
if(str1[i]==str2[j]) {matriz[j,i]<-1}
}
}
 
matriz
 
#Colocando a sequencia, aqui precisa de melhoria para
#ficar mais rapido
menor<-min(c(m,n))
 
for(k in 1:(menor-1)) {
for(i in 1:(m-1)){
for(j in 1:(n-1)) {
if(matriz[j,i]==k & matriz[j+1,i+1]>=1) {matriz[j+1,i+1]<-k+1}
}
}
}
 
matriz
 
#Extraindo uma lista com os sub-strings
lcs.n<-max(matriz)
 
lista.cs<-vector("list",lcs.n-1)
lista.cs
 
for(j in (lcs.n-1):1) {
posi<-which(matriz==j+1,arr.ind=T)
for (i in 1:nrow(posi)) {
lista.cs[[j]][i]<-paste(str1[(posi[i,2]-(j+1)+1):posi[i,2]],collapse="")
}
}
 
#Conferindo com o resto dos strings, não cheguei a investir muito tempo
#nessa parte
for (i in 1:length(unlist(lista.cs))) {
if(grepl(rev(unlist(lista.cs))[i],string3)) {
cat(rev(unlist(lista.cs))[i],"presente","n")
break
} else {cat(rev(unlist(lista.cs))[i],"ausente","n")}
}

Estatística F e a idéia da anova.

Algo que eu sempre me fez sofrer, é o fato de as vezes a gente “aprender” coisas, mas por não entender perfeitamente, ficar sem nenhuma intuição de o que está acontecendo.

Isso é o caso da estatística F, tão comum.

Aquelas tabelas de anova para mim eram um amontoado de números que não faziam o menor sentido por muito tempo. Eu só fazia a análise e olhava o valor p e depois sofria um pouco para escrever o resultado, pois nunca sabia o que deveria escrever no papel.

Mas vamos pensar melhor sobre o como funciona a estatística F do tio Fisher.

Todo modelo probabilístico é composto normalmente de uma parte determinística e uma parte probabilística, parte probabilística é o que a gente usa as distribuições estatísticas para representar, que no caso da anova é a distribuição normal. A parte determinística é uma função, que no caso da anova é as médias das coisas.

Mas vamos pensar aqui num exemplo, vamos pensar que medimos o tamanho de 10 aranhas em três locais (A,B e C), totalizando 30 aranhas.
Ai obtivemos isso:

01

Então qual o modelo mais simples possível? É termos a mesma média para os 3 locais. Para facilitar, vamos fazer um gráfico com as 30 aranhas, representando elas com a sua respectiva letrinha da onde veio. E no eixo X, vamos colocar na ordem, da primeira a trigésima aranha e vemos isso:

02

Agora isso seria o caso mais simples possível, outra possibilidade seria que cada local tenha sua própria média, que refazendo a figura acima fica assim.

03

Agora como falamos, média é 1 número, é a parte determinística, o que sobre é o resíduo. O resíduo é o que a parte determinística num deu conta de prever. No caso do modelo simples podemos representar isso como a distância da medida até a média, colocando um pauzinho no gráfico assim.

04

Agora olha so, o que fazemos com a média geral, a gente pode repetir com as médias para cada local, vendo o que restou.

05

A partir daqui já deu para ter uma idéia do que vamos comparar, vamos comparar o tamanho desses pauzinhos. Então vamos colocar os pauzinhos do modelo mais simples (só uma media, cor preta) com o do mais complexo (3 média, coloridos) lado a lado e dar uma olhada.

06

E pimba, os pauzinhos do modelo mais simples são quase sempre muito maiores que os do modelo mais complexo. E quanto menor os pauzinhos (resíduo) melhor é o modelo de certa forma, ja que precisou de menos “acaso” para explicar onde ta a medida. Ou seja aqui a gente já fica inclinado a acreditar que o modelo com 3 médias (pauzinhos coloridos, uma cor para cada local) é bem melhor.

Mas ao invés de ficar olhando assim, a gente poderia fazer o seguinte, vamos somar todos os tamanhos de pauzinhos. Para medir um pauzinho desse é só pegar o valor da medida e diminuir da média, mas como alguns vão dar valores positivos e outros negativos, a gente eleva ao quadrado.

Então somando todos os pauzinhos pretos da:

> sse<-sum((medidas-mean(medidas))^2)
> sse
[1] 112.5723

Um pauzinho ao quadrado é um quadradinho, mas o que interessa é que quanto maior os pauzinhos, maior vai ser esse número, quanto menor, menor o número.
Agora fazemos a mesma coisa com os pauzinhos coloridos:

> ssy<-sum((medidas-rep(medias[,2],each=10))^2)
> ssy
[1] 19.52451

E os pauzinho coloridos não eram muito menor? Logo o número da soma dos quadrados de todos eles da um número pequenininho uai.
Agora a estatística F é so diminuir esses 2 números.

> ssm<-sse-ssy
> ssm
[1] 93.0478

E temos o que, a diferença do modelo mais simples para o mais complexo. Dai temos que fazer mais duas coisas. Como qualquer coisa que você some mais e mais números so vai crescer mais e mais, a gente pega e divide pelo número de coisas que a gente somou. Isso são os graus de liberdade. E para ficar uma coisa onde o mais simples possível é 0 a gente ainda divide pelo soma dos pauzinhos do modelo mais complexo.

Agora ficou confuso, mas nem tanto. vejamos:

Primeiro quantas coisas somamos? Aquela diferença é a diferença entre o modelo simples com uma média e o modelo complexo que tem 3 médias diferentes, ai o grau de liberdade sempre é n-1, que não vamos discutir agora o porque desse -1. Mas então o grau de liberdade da diferença a gente define como 2, quanto mais médias, maior o grau de liberdade do numerador da divisão, a parte de cima.

O denominador, são 30 aranhas, a gente tira -1 do n de cada local e ficamos com 27. Sei que isso parece meio complicado, mas não é aqui a intuição que queria chegar ainda.

Então entramos com esses graus de liberdade:

> gld<-2
> gln<-27

Dividimos os valores das somas dos pauzinho pelos graus de liberdade

> mssm<-ssm/gld
> mssy<-ssy/gln

Ai é dividir isso que temos o valor F

> calculado<-mssm/mssy
> calculado
[1] 64.33685

Ohhhh, e para conferir, vamos fazer uma anova no R e ver se ele fala a mesma coisa ou nos perdemos nos meios dos cálculos.

Analysis of Variance Table Response: medidas Df Sum Sq Mean Sq F value Pr(>F) local 2 93.048 46.524 64.337 5.352e-11 *** Residuals 27 19.525 0.723 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Exatamente assim que se calcula o valor F, e olhem como foi significativo.
Mas porque foi significativo? Vamos pensar agora, estamos dividindo a diferença entre o modelo mais complexo pelo modelo mais simples pelo valor da soma dos pauzinhos do modelo mais complexo.

Primeira coisa a pensar é o seguinte, se o modelo mais complexo tiver pauzinhos iguais ao modelo mais simples, digo de mesmo tamanho, a diferença entre eles vai ser zero. Quanto mais diferente, maior vai ser o número da diferença, isso vai fazer o valor F crescer, ja que o numerador da divisão ta crescendo. So que tudo isso é ponderado pelo tamanho dos pauzinhos do modelo mais complexo, a soma deles.

E o que é a distribuição F, é vc tirar a atribuição de local das aranhas, embaralhar elas e ver se essa diferença fica igual. Então a gente embaralha as medidas e faz todas essa conta denovo e vê o que acontece.

07

O que acontece é que quase sempre o valor sai baixinho, olhe nos observamos um valor F de 64, e dificilmente ao acaso veríamos um valor F assim. A linha vermelha é a distribuição F, que é uma função que a gente pode plotar no gráfico. Note que para conseguir a distribuição F eu tive que ficar aleatorizando os valores, que com computadores é fácil, mas a algumas décadas atrás essa atividade deveria se,r se não impossível, pelo menos árdua. A grande sacada das distribuições é que você não precisa ficar sorteando, assumindo algumas coisas como independência das amostras e etc, a distribuição ja mostra como vai ser esse histrograma, antes de a gente fazer ele, tudo isso com calculo. Eu sempre achei que caras como o Fisher tinham uns dados e moedas na gaveta ou uma rodinha de bingo no armário, mas as distribuições nulas e valores p envolvem calculo e não jogar moedinhas ou tirar bolinhas de bingo e anotar.

Antes vamos olhar agora um exemplo de tamanhos de pauzinho para dados com diferenças significativas de médias e sem diferenças.

08

Olha la, quando não há diferença os pauzinhos são meio que do mesmo tamanho, não exatamente, mas bem parecido, enquanto as diferenças de tamanho quando existe uma diferença de médias é grande.

Ai é calcular essa diferenças, fazendo a razão entre a diferenças entre grupos dividido pelo caso mais complexo, ai comparar com a distribuição F para ver o quão intenso é isso.

O que a gente tem que falar depois, nesse caso é isso (f2,27=64.337,p<0.001). Da onde saiu o 2? É o grau de liberdade dos parâmetros, o tanto de locais que coletamos. Da onde saiu o 27? É o grau de liberdade do tanto de aranhas que coletamos, o nosso número de amostras. E o valor é a razão F, como vimos acima, que vai de 0 até o infinito, e quanto maior ela é, maior é a diferença entre o modelo mais simples e mais complexo. E o valor p é a comparação com as aleatorizações.E porque eu escrevi assim? Porque espaço em texto custa dinheiro, muito mais quando livros e revistas eram sempre impressos, minha conclusão fica a mesma eu escrevendo p=0.00000000005352 ou que esse valor ae é menor que o 0.001, p<0.001. A diferença entre escrever dos dois jeito é a economia de espaço, ao invés de quase ocupar uma linha eu ocupo só o espacinho de uma palavra.

Essa é a idéia da coisa. So que não precisamos trabalhar somente com médias, já que embaixo das repressões lineares também tem um valor F. Mas essa fica para a próxima.

#Entendendo estatistica F
#Criando dados de exemplo
set.seed(15)
local<-factor(rep(letters[1:3],each=10))
medidas<-rnorm(30,rep(c(2,4,6),each=10))
 
#Figura 01
boxplot(medidas~local,frame.plot=F,ylim=c(0,10),ylab="Medidas",xlab="Local")
 
#Figura 02
plot(medidas,pch=as.character(local),col=as.numeric(local)+1,frame.plot=F)
abline(h=mean(medidas),lwd=2)
text(25,mean(medidas)+0.2,"Média Geral")
 
#Figura 03
plot(medidas,pch=as.character(local),col=as.numeric(local)+1,frame.plot=F)
#Calculando as médias por local
medias<-aggregate(medidas,list(local),mean)
 
abline(h=medias[1,2],lwd=2,col=2)
text(25,medias[1,2]+0.2,"Média a",col=2)
abline(h=medias[2,2],lwd=2,col=3)
text(25,medias[2,2]+0.2,"Média b",col=3)
abline(h=medias[3,2],lwd=2,col=4)
text(15,medias[3,2]+0.2,"Média c",col=4)
 
#Figura 4
plot(medidas,pch=as.character(local),col=as.numeric(local)+1,frame.plot=F,
type="n")
abline(h=mean(medidas),lwd=2)
text(25,mean(medidas)-0.2,"Média Geral")
segments(1:30, mean(medidas),1:30,medidas,lwd=3)
points(medidas,pch=as.character(local),col=as.numeric(local)+1,cex=1.2)
 
#Figura 5
plot(medidas,pch=as.character(local),col=as.numeric(local)+1,frame.plot=F,
type="n")
abline(h=medias[1,2],lwd=2,col=2)
abline(h=medias[2,2],lwd=2,col=3)
abline(h=medias[3,2],lwd=2,col=4)
segments(1:30, rep(medias[,2],each=10),1:30,medidas,lwd=3,
col=as.numeric(local)+1)
text(25,medias[1,2]+0.2,"Média a",col=2)
text(25,medias[2,2]+0.2,"Média b",col=3)
text(15,medias[3,2]+0.2,"Média c",col=4)
points(medidas,pch=as.character(local),col=as.numeric(local)+1,cex=1.2)
 
#Figura 6
plot(medidas,pch=as.character(local),col=as.numeric(local)+1,frame.plot=F,
type="n")
segments(1:30, rep(medias[,2],each=10),1:30,medidas,lwd=3,
col=as.numeric(local)+1)
segments(1.5:30.5, mean(medidas),1.5:30.5,medidas,lwd=3)
 
#Tabela de Anova
sse<-sum((medidas-mean(medidas))^2)
sse
ssy<-sum((medidas-rep(medias[,2],each=10))^2)
ssy
ssm<-sse-ssy
ssm
gld<-2
gln<-27
 
mssm<-ssm/gld
mssy<-ssy/gln
 
calculado<-mssm/mssy
calculado
 
anova(lm(medidas~local))
#Criar um vetor para guardar o resultado
saida<-vector()
 
#Fazer aleatorizações
for(i in 1:5000) {
#Aqui eu resorteio sem reposição os valores de medidas
sorteio<-sample(medidas)
#E daqui em diante eu so repito tudo que fizemos acima
medias<-aggregate(sorteio,list(local),mean)
sse<-sum((sorteio-mean(sorteio))^2)
ssy<-sum((sorteio-rep(medias[,2],each=10))^2)
ssm<-sse-ssy
gld<-2
gln<-27
mssm<-ssm/gld
mssy<-ssy/27
#E vou guardando os valores F no vetor de saida
saida[i]<-mssm/mssy
}
 
#Figura 7
hist(saida,probability =T,main="Distribuição F do exemplo")
curve(df(x,2,27),0,20,add=T,col="red",lwd=2)
 
medidas2<-rnorm(30,4)
medias2<-aggregate(medidas2,list(local),mean)
 
#Figura 8
par(mfrow=c(2,1))
plot(medidas,pch=as.character(local),col=as.numeric(local)+1,frame.plot=F,
type="n")
medias<-aggregate(medidas,list(local),mean)
segments(1:30,medidas ,1:30,rep(medias[,2],each=10),lwd=3,
col=as.numeric(local)+1)
segments(1.5:30.5, mean(medidas),1.5:30.5,medidas,lwd=3)
title("Efeito Significativo")
 
plot(medidas2,pch=as.character(local),col=as.numeric(local)+1,frame.plot=F,
type="n")
segments(1:30, rep(medias2[,2],each=10),1:30,medidas2,lwd=3,
col=as.numeric(local)+1)
segments(1.5:30.5, mean(medidas2),1.5:30.5,medidas2,lwd=3)
title("Efeito não significativo")

EBImagem e a possibilidade de contar e medir células no R

A gente já viu como o R pode ajudar também a coletar dados em outro post. Algo que eu achei surpreendente foi usar o R para medir folhas como descrito nesse post ou ainda medir o crescimento de fungos nesse outro post. Essas são referencias muito melhor que a minha aqui, mas bem tenho que tentar também.

Fora a facilidade de bater uma foto e usar o computador para fazer medidas, uma característica bem legal de fazer as coisas assim é o fato da reproducibilidade dos seus resultados. Se guardamos as fotos que usamos mais o script, a qualquer momento o trabalho pode ser revisitado quanto as coletas de dados para esclarecer dúvidas, além de que quando alguém tem acesso a essas informações, repetir o trabalho se torna uma tarefa bem menos árdua.

Mas vamos la, o ponto aqui é o pacote EBImage, ele não esta no cran, ele foi feito pelo pessoal do bioconductor e ele foi feito na verdade para fazer diagnóstico de células, mas pode ser usado para muitas coisas.

Para instalar ele é muito importante ler o manual de instalação, la ele explica certinho como instalar no linux, mac e windows. O importante aqui é que ele tem dependencias fora do R, ou seja você tem que instalar coisas fora do R e essas coisas tem que estar disponível no computador para ele funcionar, mas é só seguir as instruções do manual que não tem erro.

Uma vez instalado temos que achar o que fazer com ele, fora os exemplos acima, uma amiga minha estava medindo o tamanho dos elementos de vasos dos cortes de plantinhas e eu pensei em tentar fazer isso com o R.

A imagem em questão é essa.

01

Nesse caso eu queria medir o tamanho (área branca, interna da célula) do elemento de vaso, que vemos alguns na imagem acima.

Após abrir a imagem no R, o que a gente faz é modificar a imagem e fazer filtros até separarmos as características que queremos.

Então primeiro eu usei a função thresh para deixar a imagem binária (preto e branco). Mas porque eu fiz isso? Pois é assim que esta no exemplo do pacote que eu estava tentando copiar. Devem ter outras abordagens, mas a imagem fica assim.

02

Depois existem muitas possibilidades de filtros, eu usei um chamado erode apenas, que ele faz o que? Erode o preto dentro do branco, so olhar o resultado, pontinhos brancas são todas cobertas por preto, e controlamos a intensidade disso com o makebrush, fazendo um kernel de o como essa transformação vai afetar a imagem.

03

Novamente no manual ele da varias idéias de o que as transformações fazem, nesse caso minha escolha é até meio ruim, que se comparamos com a primeira imagem, eu até diminui um pouco já a área do que quero medir. Mas aqui seria o caso de ficar testando várias possibilidades até achar uma que se adequa bem ao caso. Tem como evitar esse problema que o erode criou, mas como aqui é so um teste, não vou parar muito tempo nisso.

Então até agora deixamos a imagem pretro e branco e usamos filtros para separar bolinhas que queremos medir. Caso já estivéssemos satisfeito, agora é so separar as bolinhas (células) que queremos medir e medir.

Primeiro a gente usa a função bwlabel que vai separar e dar um “nome” para cada bolinha. E vemos:

> max(imagemlabel) [1] 112

112 bolinhas na imagem, o que é muito, estamos medindo mesmo outras células menores que não são de interesse. Com o comando computeFeatures.shape da para medir o tamanho dessas 112 bolinhas para obtermos:

> head(forma) s.area s.perimeter s.radius.mean s.radius.min s.radius.max 1 224 74 9.884910 3.1300635 18.495069 2 24 26 3.335892 0.6308161 6.511369 3 2129 178 26.563208 14.6193766 38.116709 4 1062 163 21.903856 3.7854004 38.100422 5 25 20 2.762887 0.8500000 5.071735 6 766 139 17.264785 0.9091544 28.619419

Da para ter uma ideia de como é o resultado. Daqui se fizermos um histograma de todos os tamanhos vemos que:

04

Aham, várias pequenas células estão sendo medidas, mas como a gente sabe que o nosso foco são as células grandes, é só a gente escolher um valor mínimo de tamanho (aqui o tamanho ta em pixels). Eu pego 5 mil e fazemos so um histograma com valores acima de 5 mil e temos

05

Para confirmar que não estamos medindo algo nada a ver, podemos marcar na imagem original o que medimos, só calcularmos os centroides das células que medimos e desenhar uma bolinha la. Aqui foi isso.

06

Daqui para frente ainda precisamos usar a escala na imagem para transformar essa medida de área em área em microns e todas as outras medidas que nos interessa como perímetro, raio , etc.

Por algum motivo algumas funções estão desabilitadas no pacote, mas ainda constam no manual, eu não sei porque isso, falta alguém perguntar para o autor, senão poderíamos escrever a medida do lado da célula. Numerar elas, etc.

Outra coisa legal é que para imagens muito grandes, a gente pode facilmente cortar pedaços da imagem e amostrar uma imagem, ou amostrar as medidas obtidas nessa imagem.

Bem aqui ainda faltou varias coisas e muita dedicação para usar eficientemente o pacote para essa atividade, mas a minha idéia era apenas tentar mostrar que é possível.

#Instalando e abrindo o pacote, lembre-se de ler o manual para instalar
#source("http://bioconductor.org/biocLite.R")
#biocLite("EBImage")
library(EBImage)
 
#Abrindo a imagem
imagem <- readImage("image.tif")
display(imagem,method ="raster")
 
#Tornando a imagem binaria (preto e branco)
imagem2 <- thresh(imagem[,,1], w=100, h=100, offset=0.25)
display(imagem2,method ="raster")
 
#Aplicando filtro (aqui precisa melhorar muito)
kern <- makeBrush(25, shape=’disc’)
imagem3 <- erode(imagem2, kern)
display(imagem3,method ="raster")
 
#Separando celulas
imagemlabel = bwlabel(imagem3)
max(imagemlabel)
 
#Calculando area e outras medidas
forma <- computeFeatures.shape(imagemlabel)
head(forma)
 
#histograma das areas
hist(forma[, "s.area"])
 
#contando celulas maior que 5000
sum(forma[, "s.area"]>5000)
 
#histograma das celulas maior que 5 mil
hist(forma[forma[, "s.area"]>5000, "s.area"])
 
#Calculando os centroides
centroides <- computeFeatures.moment(imagemlabel)
centroides[,1:2]
 
#Criando um vetor para separar so as celulas que queremos marcar
logico<-forma[, "s.area"]>5000
x<-centroides[logico,1]
y<-centroides[logico,2]
z<-rep(1,length(x))
 
#Desenhando um circulo nas celulas que medimos
circ<-drawCircle(imagem, x=x[1], y=y[1], radius=20, col="green", fill=T)
 
for(i in 1:length(x)) {
circ<-drawCircle(circ, x=x[i], y=y[i], radius=20, col="green", fill=T)
}
 
display(circ,method ="raster")