Rosalind – Genome Assembly as Shortest Superstring

Após mais algum tempo sem nada, vamos ver mais um problema do Rosalind sobre montar genomas. Dado que a tecnologia mais em conta agora é o sequenciamento High-throughput, a gente sempre tem um conjunto de reads para montar.

Reads são esses pedacinhos de dna, que vem da molécula original, eles tem partes repetidas que permitem a gente “montar” de volta a molécula original. Bem só lendo como funciona o processo para entender melhor, mas aqui a gente tem uma versão bem simples para começar a entender o que tem que ser feito.

Mas então, a gente tem um conjunto de vários pedacinhos de dna (reads) no formato fasta derivados de uma única molécula e temos que montar a menor sequência de dna possível (Shortest Superstring) que pode ter derivado todos esses pedacinhos.

Agora é precisa ver o seguinte, para esses dados, é garantido que existe uma única sequência original, não tem como montar de mais de um jeito, e existe uma sobreposição maior que 50% do reads, ou seja, a gente tem que encaixar no mínimo metade de um read no outro.

Vamos olhar os dados, que talvez facilite entender melhor o que pode ser feito. Esse é o conjunto de entrada.

>Rosalind_56 ATTAGACCTG >Rosalind_57 CCTGCCGGAA >Rosalind_58 AGACCTGCCG >Rosalind_59 GCCGGAATAC

E queremos fazer a menor sequência que tem todos esses reads como uma substring, então a gente pode começar a tentar encaixar, pegando os dois primeiros…

ATTAGACCTG xxx|x|xxx CCTGCCGGAA

Mas a gente pode tentar encaixar pelo outro lado também.

ATTAGACCTG x|xxxxxxx CCTGCCGGAA

A gente começa movendo só uma posição para tentar fazer o menor substring, ai como não deu, a gente diminuir em dois o tamanho do encaixe,

ATTAGACCTG xxxx||x| CCTGCCGGAA

Veja que o read 56 encaixa no 58 por 7 bases, e como os reads tem 10 bases, é mais da metade.

ATTAGACCTG ||||||| AGACCTGCCG

E encaixando tudo, a gente monta a sequência inicial

>Rosalind_56 ATTAGACCTG >Rosalind_58 |||AGACCTGCCG >Rosalind_57 ||||||CCTGCCGGAA >Rosalind_59 |||||||||GCCGGAATAC ATTAGACCTGCCGGAATAC

A gente pode fazer uma função recursiva no python.

def long(reads, genoma=''):
	if len(reads) == 0:
		return genoma
	if len(genoma)==0:
		genoma = reads.pop(0)
		return long(reads, genoma)  
	for i in range(len(reads)):
		a = reads[i]
		l = len(a)
		for p in range(l / 2):			
			q = l - p
			if genoma.startswith(a[p:]):
				reads.pop(i)
				return long(reads, a[:p] + genoma)
			if genoma.endswith(a[:q]):
				reads.pop(i)
				return long(reads, genoma + a[q:])

A função é assim, primeiro o nosso caso de parada, que é quando não temos mais reads para montar

	if len(reads) == 0:
		return genoma

Quando estamos começando, só tem um jeito de montar, logo a gente pode começar do primeiro read

	if len(genoma)==0:
		genoma = reads.pop(0)
		return long(reads, genoma)

Ai , para todos os reads

	for i in range(len(reads)):
		a = reads[i]
		l = len(a)

Para pelo menos metade do tamanho dos reads

		for p in range(l / 2):

A gente começa diminuindo uma base para encaixar, depois 2, e assim por diante

			q = l - p

Ai tentamos encaixar na frente

			if genoma.startswith(a[p:]):
				reads.pop(i)
				return long(reads, a[:p] + genoma)

E encaixar atrás

			if genoma.endswith(a[:q]):
				reads.pop(i)
				return long(reads, genoma + a[q:])

Ai veja que só existe uma maneira de montar, logo a gente vai entrar nessa recursão somente uma vez, porque uma vez que o read seja unido, a gente tira ele da lista de reads e já junta ao genoma, e outra coisa, é que so tem um read para cada ponta.

Bem é isso ai, o script vai estar la no repositório recologia, e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
#!/usr/bin/env python
'''
 
'''
 
def long(reads, genoma=''):
	if len(reads) == 0:
		return genoma
	if len(genoma)==0:
		genoma = reads.pop(0)
		return long(reads, genoma) 
	for i in range(len(reads)):
		a = reads[i]
		l = len(a)
		for p in range(l / 2):			
			q = l - p
 
			if genoma.startswith(a[p:]):
				reads.pop(i)
				return long(reads, a[:p] + genoma)
 
			if genoma.endswith(a[:q]):
				reads.pop(i)
				return long(reads, genoma + a[q:])
if __name__ == "__main__":
 
	arquivo = open("/home/augusto/Downloads/rosalind_long.txt")
	sequencia=[]
	i=-1
	for linha in arquivo:
		if linha.find('>') != -1:		
			i+=1
			sequencia.append('')
		else:		
			sequencia[i]=sequencia[i]+linha.rstrip()
	print long(sequencia)

Poisson Point Process

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

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

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

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

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

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

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

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

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

Lembrando, a distribuição normal é:

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

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

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

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

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

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

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

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

Depois descemos a potência do lambda la de cima

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

simplificamos

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

ai tiramos o delta de dentro da raiz

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

subimos ele

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

cortamos o um mais delta

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

e arrumamos as potências

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

distribution_chart

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

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

A gente vai ver o seguinte:

01

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

Referência:

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

Rosalind – Suboptimal Local Alignment

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

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

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

>Rosalind_12 GACTCCTTTGTTTGCCTTAAATAGATACATATTTACTCTTGACTCTTTTGTTGGCCTTAAATAGATACATATTTGTGCGACTCCACGAGTGATTCGTA >Rosalind_37 ATGGACTCCTTTGTTTGCCTTAAATAGATACATATTCAACAAGTGTGCACTTAGCCTTGCCGACTCCTTTGTTTGCCTTAAATAGATACATATTTG

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

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

Pegamos a segunda sequencia:

>Rosalind_37 ATGGACTCCTTTGTTTGCCTTAAATAGATACATATTCAACAAGTGTGCACTTAGCCTTGCCGACTCCTTTGTTTGCCTTAAATAGATACATATTTG

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

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

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

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

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

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

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

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

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

Então veja que se a gente tem

>A GACTCCTTTGTTTGCCTTAAATAGATACATATTTACTCTT >B ATGGACTCCTTTGTTTGCCTTAAATAGATACATATTCAAC

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

>A GACTCCTTTGTTTGCCTTAAATAGATACATATTTACTCTT >B ATGGACTCCTTTGTTTGCCTTAAATAGATACATATTCAAC

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

>A GACTCCTTTGTTTGCCTTAAATAGATACATATTTACTCTT >B ATGGACTCCTTTGTTTGCCTTAAATAGATACATATTCAAC

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

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

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

itnotstupid

Bem é isso ai, o script vai estar la no repositório recologia, e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail.

Referência:
Suboptimal Local Alignment

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
#!/usr/bin/env python
'''
 
'''
 
def hamming_dist(a,b,limite_erros=5):
    erros = 0
    for i in range(len(a)):
        if a[i] != b[i]:
            erros += 1
            if erros>limite_erros:
            	return erros
    return erros 
 
 
arquivo = open("/home/augusto/Downloads/rosalind_subo.txt")
 
sequencia=[]
i=-1
 
for linha in arquivo:
	if linha.find('>') != -1:		
		i+=1
		sequencia.append('')
	else:		
		sequencia[i]=sequencia[i]+linha.rstrip()
 
 
def subo(a, b):
	contagem_final=0
	for tamanho in range(32,41):
		for i in range(len(a)-tamanho):
			contagem=0
			for j in range(len(b)-tamanho):
				if hamming_dist(a[i:i+tamanho], b[j:j+tamanho])<=3:
					contagem+=1
					print a[i:i+tamanho]
					if contagem>contagem_final:
						contagem_final=contagem
	return contagem_final
print subo(sequencia[1], sequencia[0])
print subo(sequencia[0], sequencia[1])

Mapas aleatórios e o número de componentes

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

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

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

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

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

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

01

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

02

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

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

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

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

Entao temos um pixel e 8 pixels em volta dele

1 2 3 4 x 5 6 7 8

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

03

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

(linha,coluna)

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

Usando a função do igraph, teremos:

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

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

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

04

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

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

Bem é isso ai, apesar dos mapas serem gerados aleatoriamente, ainda sim temos algumas certezas, o que permite alguma previsibilidade, mesmo o script vai estar la no repositório recologia, e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail.

Referência:
Robert Sedgewick Kevin Wayne 2011 Algorithms 4th ed Addison-Wesley Professional 992 pp

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
#####################################
## Conectividade
#####################################
library(igraph)
 
##Função para gerar mapas quadrados aleatorios, n é o tamanho e p é a probabilidade de threshold
gera_mapa<-function(n,p) {
    matriz<-matrix(0,n,n)
    for(i in 1:nrow(matriz)){
        for(j in 1:ncol(matriz)){
            if(runif(1)<=p){
                matriz[i,j]<-1
            }
        }
    }
    return(matriz)
}
 
##Função para transformar um mapa em um grafo, com pixels com vértices e pontos adjacentes ligados por arestas
matriz_para_grafo<-function(matriz){
    aresta<-1
    grafo<-data.frame(origem=vector(),destino=vector())
    for(i in 1:nrow(matriz)){
        for(j in 1:ncol(matriz)){        
            if(matriz[i,j]!=0){            
                for(m in -1:1){
                    for(n in -1:1){
                        if(i+m>0 & i+m<=nrow(matriz) &
                           j+n>0 & j+n<=ncol(matriz)){
                            if(matriz[i+m,j+n]!=0){
                                grafo[aresta,]<-c(paste("(",i,",",j,")",sep=""),paste("(",i+m,",",j+n,")",sep=""))
                                aresta<-aresta+1
                            }
                        }
                    }
                }
            }
        }
    }
    indice<-!grafo[,1]==grafo[,2]
    indices<-which(mapa!=0,arr.ind = T)    
    grafo<-simplify(graph.data.frame(grafo[indice,],directed=F,vertices=paste("(",indices[,"row"],",",indices[,"col"],")",sep="")))  
    return(grafo)
}
 
##Função para fazer uma imagem da matriz como vemos ela na tela
plota_mapa<-function(mapa){
    image(t(mapa)[,ncol(mapa):1],col=0:1)
    }
 
 
#####################################
## Mapas aleatorios
#####################################
set.seed(123)
 
##Gera um mapa com threshold de 0.5
mapa<-gera_mapa(10,0.5)
mapa
#jpeg("01.jpg")
plota_mapa(mapa)
#dev.off()
 
##Threshold de 0.2
mapa<-gera_mapa(10,0.2)
#jpeg("02.jpg")
plota_mapa(mapa)
#dev.off()
 
##Transformando o mapa em grafo
grafo<-matriz_para_grafo(mapa)
 
#jpeg("03.jpg")
plot(grafo,vertex.size=20)
#dev.off()
 
##Usando a função do igraph
components(grafo)
 
##Realizando uma simulção com 10 repetições para probabilidades de 0.1 a 0.9
repeticoes<-10
probs<-seq(0.1,0.9,0.05)
tamanho<-15
 
##Matrizes para receber os resultados da simulação
n_compomentes<-matrix(0,ncol=repeticoes,nrow=length(probs),dimnames = list(paste("p=",probs,sep=""),paste("Amostra",1:10)))
maior_componente<-matrix(0,ncol=repeticoes,nrow=length(probs),dimnames = list(paste("p=",probs,sep=""),paste("Amostra",1:10)))
 
##Realizando a simulação
for(i in 1:length(probs)){
    for(j in 1:repeticoes){
        mapa<-gera_mapa(tamanho,probs[i])
        grafo<-matriz_para_grafo(mapa)
        n_compomentes[i,j]<-components(grafo)$no
        maior_componente[i,j]<-max(components(grafo)$csize)
        print(paste("Concluido i=",i," j=",j,sep=""))
    }
}
 
##Avaliando o resultado
#jpeg("04.jpg")
par(mfrow=c(2,1))
plot(probs,rowMeans(n_compomentes),type="b",xlim=c(0,1),ylab="Número médio de componentes",xlab="Threshold geração mapa",frame=F)
plot(probs,rowMeans(maior_componente),type="b",xlim=c(0,1),ylab="Tamanho do maior componente",xlab="Threshold geração mapa",frame=F)
#dev.off()

Rosalind – Data Formats

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

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

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

A gente recebe uma lista de ids

FJ817486 JX069768 JX469983

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

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

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

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

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

from Bio import Entrez
from Bio import SeqIO

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

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

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

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

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

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

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

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

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

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

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

Bem é isso ai, o script vai estar la no repositório recologia, e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
#!/usr/bin/env python
'''
 
'''
 
from Bio import Entrez
from Bio import SeqIO
 
arquivo = open('/home/augusto/Downloads/rosalind_frmt.txt')
 
ids=[]
for linha in arquivo:
	ids.extend(linha.rstrip().split())
 
Entrez.email = 'ribas.aca@gmail.com'
handle = Entrez.efetch(db='nucleotide', id=ids, rettype='fasta')
records = list(SeqIO.parse(handle, 'fasta'))
 
menor=10**6
indice=-1
for i in range(len(records)):
	if len(records[i].seq) < menor:
		menor=len(records[i].seq)
		indice=i
 
print indice
 
sequencia=records[indice].seq
lista=[]
while sequencia:
    lista.append(sequencia[:70])
    sequencia = sequencia[70:]
 
print '>'+records[indice].description
for linha in lista:
	print linha