Random Web Surfer, como o google classifica suas páginas

Eu fiz algumas aulas no coursera, e em computação, as que mais me chamaram atenção foram as do Robert Sedgewick, muito famoso acho que pela tese dele sobre a análise de algorítimos do Quicksort, que ele mostra que é em média mais rápido que muitos outros, e acho que o mais comum em uso hoje. Bem também depois da aula eu comecei a olhar os livros deles, que são bem legais, entre eles o Introduction to Programming in Java, que é legal principalmente pelos exercícios que ele propõem.

Um que eu achei muito legal foi o Random Web Surfer, que eu vou colocar aqui, adaptando o código para python.

Nos primórdios da internet, a gente tinha os listão de página, para conseguir acessar, veja o post do pi-hole para mais alguma informação, mas, então, como a gente tinha uma lista, chegou a hora de ordenar ela, ordem alfabética não era muito útil, talvez classificar por tipo de páginas, mas ainda sim fica a pergunta, qual seria a melhor, qual consultar primeiro? Os primeiros web search engines usaram o critérios de quanto mais links apontavam para uma página, melhor ela era, porque todo mundo apontava para ela, logo a gente sempre ia parar nele, logo ela era muito boa. Porém todo mundo começou a burlar isso, se alguém queria ficar pra cima nesse critério, era só pagar para outras páginas colocarem um link para sua página, e rapidamente você estava no topo, ou ainda você poderia criar mil páginas mequetrefes, e colocar o link que quiser nelas, e assim controlar quem está no topo da lista de uma busca da internet, e logo esse artificio começou dificultar a vida dos buscadores de página. Mas o google usou uma forma diferente para classificar páginas.

Bem como funciona a internet, temos uma conjunto de páginas, que tem coneções umas para as outras.

O primeiro passo é como representar essa estrutura, que para nossa sorte pode ser representado facilmente usando um grafo, onde as páginas são vértices e os link entre elas são arestas, um grafo direcional.

Então a entrada pode ser um arquivo onde na primeira linha a gente tem o número de páginas que vamos trabalhar, e depois todos as arestas existentes, qual página tem um link que leva para qual página.

5 0 1 1 2 1 2 1 3 1 3 1 4 2 3 3 0 4 0 4 2

Bem então temos uma representação da internet, o que fazer com ela? De forma a não cair na armadilha anterior, que é contar o quantos links chegam a você?
Uma forma é simular como uma pessoa navega, uma pessoa entra numa página, ai pode clicar num link para outra página, ou digitar o endereço de uma nova página. Mas esses eventos tendem a acontecer com uma distribuição de pareto, a maior parte do tempo, você segue links para as próximas páginas que vai visitar, mas de vez em quanto você simplismente entra numa página por vontade própria, isso numa proporção mais ou menos de 90% do tempo seguindo links e 10% do tempo entrando em novas páginas. Seguindo essa lógica, podemos fazer uma matriz de transição entre páginas, para realizar uma simulação de como o navegador se comporta, e contabilizar quais páginas ele visita mais frequentemente.

Diferente do que uso aqui, numa programação mais em um script continuo realizando as coisas passo a passo, vamos tentar usar um pouco de orientação a objetos, que nada mais é um paradigma de programação para deixar as coisas mais organizadas.

Bem vamos começar criando uma classe em python, e vamos ter um script que vai usar essa classe.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
 
#!/usr/bin/env python
 
import sys
import random
 
class Randomwebsurfer(object):
    #Aqui vai métodos(funções da minha classe)
 
#####################################
## Random Surfer
#####################################
 
if __name__ == '__main__':
    #Aqui eu faço o script, como normalmente fazemos.

Alguns detalhes é que essa linha aqui:

1
#!/usr/bin/env python

Avisa o sistema operacional, quem é o interpretador desse programa, ou seja, se a gente já avisou quem é o interpretador, ao inves de digitar python programa.py, a gente so precisa digitar ./programa.py no terminal.

Bem toda classe tem seus atributos, ou seja a informação dela, nesse caso tem que ser o grafo que representa a relação entre as páginas, que podemos representar como uma matriz, que em python vai ser uma lista de listas, e para não ficar olhando toda hora a matriz para saber quantas páginas temos, vamos anotar o numero de páginas também, como python não precisa realmente declarar os atributos da classe, eu coloquei apenas para não me perder do que tem na classe.

1
2
3
4
5
6
7
8
9
10
11
12
13
class Randomwebsurfer(object):
    n=None
    matriz=None
 
Além disso, toda classe tem que ter um construtor, uma função que vai rodar toda vez que um objeto dessa classe for instanciado, que em classes de python a gente define o __init__
 
class Randomwebsurfer(object):
    n=None
    matriz=None
 
    def __init__(self,n):
        self.n=n
        self.matriz=self.inicia_matriz(n)

Agora toda vez que eu instancio um objeto dessa classe, por exemplo

1
objeto=Randomsurfer(5)

Eu tenho que rodar esse método (para classes a gente chama funções de métodos, mas é a mesma coisa).

Mas o que eu quero fazer, eu quero que alguém me diga quantas páginas existem, para fazer uma matriz quadrada de zeros, para anotar nela os links.

1
2
3
4
5
6
7
8
9
10
11
class Randomwebsurfer(object):
    n=None
    matriz=None
 
    def __init__(self,n):
        self.n=n
        self.matriz=self.inicia_matriz(n)
 
    def inicia_matriz(self,n):
        matriz = [[0 for i in xrange(n)] for j in xrange(n)]
        return matriz

Para tanto eu fiz esse método chamado inicia_matriz, que inicia a matriz de zeros, então dentro da classe eu chamo um método que so existe nela mesma, por isso o “self.”, assim o python sabe onde procurar, nesse caso nela mesma, na classe.

Para ver o que fizemos, podemos fazer um método para imprimir essa matriz, para ver como ela é.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
class Randomwebsurfer(object):
    n=None
    matriz=None
 
    def __init__(self,n):
        self.n=n
        self.matriz=self.inicia_matriz(n)
 
    def inicia_matriz(self,n):
        matriz = [[0 for i in xrange(n)] for j in xrange(n)]
        return matriz
 
    def imprime_matriz(self):
        for linha in self.matriz:
            print ''
            for item in linha:
                print item,
        print ''

Agora temos que ter um método para adiconar links também nessa matriz, porque ela não pode ser zerada. Então fazemos mais um método.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
class Randomwebsurfer(object):
    n=None
    matriz=None
 
    def __init__(self,n):
        self.n=n
        self.matriz=self.inicia_matriz(n)
 
    def inicia_matriz(self,n):
        matriz = [[0 for i in xrange(n)] for j in xrange(n)]
        return matriz
 
    def adiciona_link(self,linha,coluna):
        self.matriz[linha][coluna]+=1
 
    def imprime_matriz(self):
        for linha in self.matriz:
            print ''
            for item in linha:
                print item,
        print ''

Assim adicionamos o método adiciona_link, que recebe o link, da onde ele vem, que é a linha e para onde ele vai, que é uma coluna, ai adicionamos a esse item 1, que é um link, uma aresta, uma conexão entre páginas.

Agora nos vamos precisar calcular a matriz de transição, para conseguir realizar a simulação, a matriz de transição vai ter que incorporar aquelas porcentagens de 10% de ir parar por acaso em qualquer página e 90% de seguir um link na página.

As porcentagens ao acaso é fácil, basta dividir esse 10% pelo número de páginas, que fica uma pouquinho de porcentagem para cada um, se temos 5 páginas, 2% para cara uma. Já os 90% temos que dividir de forma ponderada, pelo número de vezes que outra página é linkada a ela. Para isso podemos calcular o grau de saída da página, dividimos 90% por esse número, que é dividir 90% por quantos links tem nessa página, e ai multiplicamos pela linha da matriz, o valor que esta na linha, se é zero, não há chance, se tem 2 links a chance é maior que uma que tem 1 link.

Então no código fazemos assim, um método para calcular o grau:

1
2
3
4
5
6
    def calcula_grau(self):
        self.grau=[0]*self.n
 
        for i in range(self.n):
            for j in range(self.n):
                self.grau[i]+=self.matriz[i][j]

E um método para iniciar a matriz de transição, onde calculamos o grau, fazemos uma matriz de pulo ao acaso, que é dividir o 10% entre todas as páginas, uma de link que faz a conta que eu disse, divide o 90% pelo grau da página, todos os links que ela tem, e dai multiplica pela quantidade de link para cada local, feito essas duas, a gente simplesmente soma essa duas matrizes que temos a matriz de transição.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
    def inicia_matriz_transisao(self):
        self.calcula_grau()        
        pulo=[[0 for i in xrange(n)] for j in xrange(n)]
        link=[[0 for i in xrange(n)] for j in xrange(n)]
        self.transisao=[[0 for i in xrange(n)] for j in xrange(n)]
 
        for i in range(self.n):
            for j in range(self.n):
                pulo[i][j]=0.1/float(self.n)
 
        for i in range(self.n):
            for j in range(self.n):
                link[i][j]=self.matriz[i][j] * (0.9/float(self.grau[i]))
 
        for i in range(self.n):
            for j in range(self.n):
                self.transisao[i][j]=pulo[i][j]+link[i][j]

Agora veja que ela é a essência da simulação, então ela tem que ser um atributo também, além é claro, do resultado da simulação, então já vamos deixar um frequencia, para guardar o resultado da simulação.

1
2
3
4
5
6
class Randomwebsurfer(object):
    n=None
    matriz=None
    grau=None
    transisao=None
    frequencia=None

Agora vamos finalmente fazer a simulação, que vai ser assim:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
 def simulacao(self,passos):
        self.inicia_matriz_transisao()
 
        if self.frequencia==None:
            self.frequencia=[0]*self.n
 
        pagina=0
        i=0
        while i < passos:
            sorteio=random.random()
            soma=0.0
 
            for j in xrange(self.n):
                soma+=self.transisao[pagina][j]
                if sorteio < soma:
                    pagina=j
                    break
 
            self.frequencia[pagina]+=1
            i+=1

Primeiro, a gente não vai navegar para sempre, então número de passos é o número de páginas que vamos simular, precisamos desse valor para realizar a simulação.
Nesse momento o modelo da internet, as páginas já tem que estar inserido, acho que poderíamos colocar um erro aqui, caso esse não seja o caso, mas depois eu faço as exceções.

A gente usa nosso proprio método de iniciar a matriz de transisão, depois, se não temos um local para guardar resultado, nos o iniciamos, eu fiz dessa forma testando se o resultado ja existia porque imaginei que poderiamos continuar a simulação depois, mas ok.

E temos que ter um ponto de partida, que nesse caso é a página zero, que sempre será a inicial.

Certo então começamos a simulação no while, nele a gente sorteia um numero de zero a um, depois ve em que página caiu. Que é essa parte em espeficio

1
2
3
4
5
            for j in xrange(self.n):
                soma+=self.transisao[pagina][j]
                if sorteio < soma:
                    pagina=j
                    break

suponha que temos três páginas com as porcentagens [20% 30% 50%], eu sorteio um número ao acaso entre zero e um, se caiu 0.97, qual página ele é? A porcentagem pode ser dividida assim:

20% 30% 50% [----------0.2------------0.5------------------------------------1] Página 1 Página 2 Página 3

Então temos intervalos entre zero e um, cada um de tamanho correspondente a porcentagem dada para cada página, e dependendo onde o ponto cai, é qual página a pessoa foi parar.

Por isso o for vai somando as porcentagens e comparando com o ponto, até encontrar a página.

Feito isso a gente faz um update de qual página estamos e somamos 1 a frequência com que a página é visitada, ou seja, é essa simulação se comporta como uma cadeia de markov, é um random walk sem memoria das páginas de internet.

Agora que temos nossa classe arrumadinha, podemos utilizá-la da seguinte forma

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
if __name__ == '__main__':
 
    print "Lendo entrada..."
    arquivo = open(sys.argv[1], 'r')
    n=int(arquivo.readline())
    print "Numero de paginas: " +  str(n)
 
    minharede=Randomwebsurfer(n)
 
    for linha in arquivo:
        vetor=linha.strip('\n')
        vetor=vetor.split(' ')        
        while '' in vetor: vetor.remove('')
 
        for i in range(0,len(vetor),2):
            minharede.adiciona_link(int(vetor[i]),int(vetor[i+1]))
 
 
    print "Matriz de entrada:",
    minharede.imprime_matriz()
 
    print ''
    print 'Rodando simulacao...'
    minharede.simulacao(100)
 
    print ''
    print 'Resultado:'
    minharede.imprime_resultado_simulacao()

A gente faz uma método main, que é o que o python realmente vai rodar, é o nosso script. Bem a entrada vem de um arquivo, naquel formato que falamos, usandos o sys.argv para ler o nome do arquivo que a pessoa fornece ao utilizar o programa.

1
arquivo = open(sys.argv[1], 'r')

Na primeira linha temos o número de páginas

1
2
n=int(arquivo.readline())
print "Numero de paginas: " +  str(n)

A gente le e imprime ele, para ver se ta fazendo correto, lembre-se que lemos tudo como string, então temos que converter para números.
Iniciamos um objeto da nossa classe recem criada

1
minharede=Randomwebsurfer(n)

E vamos depois lendo os pares de inteiro que são os links das páginas e adicionando ao modelo com o método adiciona_link

1
2
3
4
5
6
7
for linha in arquivo:
        vetor=linha.strip('\n')
        vetor=vetor.split(' ')        
        while '' in vetor: vetor.remove('')
 
        for i in range(0,len(vetor),2):
            minharede.adiciona_link(int(vetor[i]),int(vetor[i+1]))

Com isso pronto, a gente pode imprimir a matriz que representa o modelo, pra ver se está correta

1
2
print "Matriz de entrada:",
minharede.imprime_matriz()

Realizamos a simulação, com 100 passos para experimentar

1
2
print 'Rodando simulacao...'
minharede.simulacao(100)

E imprimimos o resultado para olhar

1
2
print 'Resultado:'
minharede.imprime_resultado_simulacao()

Então utilizando nosso programinha, vemos o seguinte:

augusto@augusto-xps:~/git/recologia$ ./randomwebsurfer.py tiny.txt Lendo entrada... Numero de paginas: 5 Matriz de entrada: 0 1 0 0 0 0 0 2 2 1 0 0 0 1 0 1 0 0 0 0 1 0 1 0 0 Rodando simulacao... Resultado: 25 30 14 24 7

E pronto, temos uma pequena implementação do Random Surfer.

Acho que depois eu faço um post discutindo esse resultado, porque esse post ta ficando muito grande já, mas é legal programar um pouquinho, eu ando muito preguiçoso e a gente vai esquecendo como fazer as coisas se não treina, mas depois eu continuo esse post. Bem é isso ai, o script vai estar la no repositório recologia, vejam o exercicio orignal aqui , provavelmente com uma explicação melhor que a minha, mas eu queria apensar programar alguma coisa aqui para passar algum tempo, e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail.

Referência:

Robert Sedgewick e Kevin Wayne 2007 An Introduction to Programming in Java: An Interdisciplinary Approach – Pearson Addison Wesley 736pp

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
142
143
144
145
146
#!/usr/bin/env python
 
import sys
import random
 
class Randomwebsurfer(object):
    n=None
    matriz=None
    grau=None
    transisao=None
    frequencia=None
 
    def __init__(self,n):
        self.n=n
        self.matriz=self.inicia_matriz(n)
 
 
    def simulacao(self,passos):
        '''
        Esse método realiza a simulacao de navegação por um usuário.
 
        '''
        self.inicia_matriz_transisao()
 
        if self.frequencia==None:
            self.frequencia=[0]*self.n
 
        pagina=0
        i=0
        while i < passos:
            sorteio=random.random()
            soma=0.0
 
            for j in xrange(self.n):
                soma+=self.transisao[pagina][j]
                if sorteio < soma:
                    pagina=j
                    break
 
            self.frequencia[pagina]+=1
            i+=1
 
 
    def inicia_matriz(self,n):
        '''
        Esse método inicia a matriz que representa a rede, mas sem nenhuma conecao entre paginas,
        recebe um argumento que é o numero de paginas do sistema
 
        '''
        matriz = [[0 for i in xrange(n)] for j in xrange(n)]
        return matriz
 
    def adiciona_link(self,linha,coluna):
        '''
        Esse método adiciona uma conecao entre duas paginas, sendo o primeiro argumento a origem e o segundo o destino
 
        '''
        self.matriz[linha][coluna]+=1
 
    def calcula_grau(self):
        '''
        Esse método calcula o grau de cada pagina.
 
        '''
        self.grau=[0]*self.n
 
        for i in range(self.n):
            for j in range(self.n):
                self.grau[i]+=self.matriz[i][j]
 
    def inicia_matriz_transisao(self):
        '''
        Esse método calcula a matriz de transicao, a matriz e links precisam estar todos adicionados
 
        '''
        self.calcula_grau()        
        pulo=[[0 for i in xrange(n)] for j in xrange(n)]
        link=[[0 for i in xrange(n)] for j in xrange(n)]
        self.transisao=[[0 for i in xrange(n)] for j in xrange(n)]
 
        for i in range(self.n):
            for j in range(self.n):
                pulo[i][j]=0.1/float(self.n)
 
        for i in range(self.n):
            for j in range(self.n):
                link[i][j]=self.matriz[i][j] * (0.9/float(self.grau[i]))
 
        for i in range(self.n):
            for j in range(self.n):
                self.transisao[i][j]=pulo[i][j]+link[i][j]       
 
    def imprime_matriz(self):
        '''
        Esse método imprime a matriz de dados que presenta o sistema
 
        '''
        for linha in self.matriz:
            print ''
            for item in linha:
                print item,
                print ''
 
    def imprime_resultado_simulacao(self):
        '''
        Esse método imprime o resultado da simulação
 
        '''
        for i in self.frequencia:
            print i, 
 
 
 
 
#####################################
## Random Surfer
#####################################
 
if __name__ == '__main__':
 
    print "Lendo entrada..."
    arquivo = open(sys.argv[1], 'r')
    n=int(arquivo.readline())
    print "Numero de paginas: " +  str(n)
 
    minharede=Randomwebsurfer(n)
 
    for linha in arquivo:
        vetor=linha.strip('\n')
        vetor=vetor.split(' ')        
        while '' in vetor: vetor.remove('')
 
        for i in range(0,len(vetor),2):
            minharede.adiciona_link(int(vetor[i]),int(vetor[i+1]))
 
 
    print "Matriz de entrada:",
    minharede.imprime_matriz()
 
    print ''
    print 'Rodando simulacao...'
    minharede.simulacao(100)
 
    print ''
    print 'Resultado:'
    minharede.imprime_resultado_simulacao()

Rosalind – Enumerating k-mers Lexicographically

Esse problema é muito simples, dada um alfabeto, uma coleção com um certo número de símbolos e um tamanho, temos que fazer todas as combinações possíveis organizados lexicograficamente, na “ordem alfabética”.

Primeiro, a quantidade de possibilidades será o número de símbolos elevado pelo tamanho da palavra que vamos formar.

O exemplo de dados é:

A C G T 2

Então o alfabeto é o ‘A’, ‘C’, ‘G’ e ‘T’ e o tamanho é 2, ou seja temos 2 posições

_ _

E cada posição tem quatro possibilidades.

A primeira solução que eu pensaria é fazer um loop, usando um for para cada posição, algo como:

1
2
3
4
5
alfabeto=['a','c','g','t']
 
for i in alfabeto:
    for j in alfabeto:
        print i+j

Mas ai temos o problema de como criaríamos for de acordo com o tamanho de entrada, já que ele é variável. Mas podemos criar uma solução recursiva para esse problema.

1
2
3
4
5
6
7
8
def fn_re(n):
 
    if n>1:        
        fn_re(n-1)
 
    print n
 
fn_re(len(alfabeto))
>>>4 3 2 1

Ou seja, enquanto n for maior que o tamanho for maior que zero, a gente faz algo, até a condição de parada, que é n chegar a zero, não ter o que fazer, ai a gente só retorna o processamento. Então o que a gente precisa é mandar o alfabeto, e mandar o tamanho da palavra que vamos formar, ai para cada letra do alfabeto, a gente chama de novo a função.

1
2
3
4
5
6
7
def combina(alfabeto,n,kmer,lista):
    if n==0:
        lista.append(kmer)
    else:
        for letra in alfabeto:
            combina(alfabeto,n-1,kmer+letra,lista)
    return lista

Um jeito legal de entender melhor o que está acontecendo, é imprimir o kmer e a lista, no inicio da função, assim:

1
2
3
4
5
6
7
8
9
def combina(alfabeto,n,kmer,lista):
    print kmer
    print lista
    if n==0:
        lista.append(kmer)
    else:
        for letra in alfabeto:
            combina(alfabeto,n-1,kmer+letra,lista)
    return lista

Ai quando a gente usa ela

1
2
3
alfabeto=['a','c','g','t']
n=2
resposta= combina(alfabeto,n,'',[])

A gente vê o seguinte:

[] a [] aa [] ac ['aa'] ag ['aa', 'ac'] at ['aa', 'ac', 'ag'] c ['aa', 'ac', 'ag', 'at'] ca ['aa', 'ac', 'ag', 'at'] cc ['aa', 'ac', 'ag', 'at', 'ca'] . . .

A primeira chamada, o kmer é ” e a lista ta vazia.
O n é diferente de zero, então a gente continua e então chamamos a função para todo o alfabeto, no primeira chamada recursiva, vamos diminuir um no n, que vai ficar 1 e veja que vamos concatenar no for a uma letra ao kmer, a primeira letra do nosso alfabeto é o a, então na segunda vez que imprimimos na tela, vemos o a, porque fizemos kmer+letra na chamada e temos o n sendo 1, então chamamos recursivo de novo, só que agora, na recursão o kmer é ‘a’, e no for vamos fazer kmer+letra de novo, e chamar recursivo, mas agora temos ‘aa’, agora a gente vai cair dentro do if, e ao invés de fazer a chamada recursiva, a gente so coloca o ‘aa’ na lista, e assim seguimos.

Bem agora é só ler o arquivo, pegar os parâmetros e pimba, 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
def combina(alfabeto,n,kmer,lista):
    print kmer
    print lista
    if n==0:
        lista.append(kmer)
    else:
        for letra in alfabeto:
            combina(alfabeto,n-1,kmer+letra,lista)
    return lista
 
"""
file = open( "/mnt/B21AA1BD1AA17ECB/downloads_chromium_ubuntu/rosalind_lexf.txt" )
alfabeto = file.readline().split()
n = int(file.readline())
file.close()
"""
 
###Exemplo
alfabeto=['a','c','g','t']
n=2
 
 
print "Número de possibilidades: "+str(len(alfabeto)**n)
 
resposta= combina(alfabeto,n,'',[])
 
for i in resposta:
    print i

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)

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])

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