Rosalind – Transitions and Transversions

Mutações são classificadas de várias formas, se ela alteram alguma proteína, fazem alguma coisa mesmo, ou são silenciosas.

Mas um tipo especial de classificação é se elas são transições ou transversões

transitions-transversions

Transições são a troca de uma purina por outro purina ou uma pirimidina por outra pirimidina, ou seja, trocas dentro dos próprios grupos de bases nitrogenadas.

Transversões são trocas entre os grupos, purinas por pirimidinas e vice e versa.

Primeiro, podemos ver pela ilustração que temos mais possibilidades de transversões do que transições, se fosse tudo ao acaso, mas temos ainda implicações químicas também para elas ocorrem.

Agora se tiver duas sequências do mesmo tamanho, podemos contar cada uma desses coisas (transições e transversões) acontecendo, e fazer uma razão, que é dividir o número de um pelo número do outro.

Apesar de algo que parece bobo, isso pode ser uma medida para ver se estamos encontrando mutações corretamente num genoma sequenciado, ja que um valor muito diferente de algo esperado para a espécie, pode ser sinal de problemas. Podemos ver no manual do samtools que…

“A good measure of the callset quality is often ts/tv, the ratio of transitions and transversions.”

Mas chega de enrolação, dado que temos duas sequências.

>Rosalind_0209 GCAACGCACAACGAAAACCCTTAGGGACTGGATTATTTCGTGATCGTTGTAGTTATTGGA AGTACGGGCATCAACCCAGTT >Rosalind_2200 TTATCTGACAAAGAAAGCCGTCAACGGCTGGATAATTTCGCGATCGTGCTGGTTACTGGC GGTACGAGTGTTCCTTTGGGT

Temos que contar quando ocorre transições ou transversões, como elas tem o mesmo tamanho esta fácil, e não temos que se preocupar com inserções e deleções.

Então a gente tem que colocar essa sequência uma do lado da outra e olhar caso a caso

GCAACGCACAACGAAAACCCTTAGGGACTGGATTATTTCGTGATCGTTGTAGTTATTGGAAGTACGGGCATCAACCCAGTT | TTATCTGACAAAGAAAGCCGTCAACGGCTGGATAATTTCGCGATCGTGCTGGTTACTGGCGGTACGAGTGTTCCTTTGGGT

Esse primeiro caso, um G virou um T, Bem sabemos que:
Purinas = Adenina e Guanina
Pirimidinas = Citosina, Timina e Uracilo

Então um G virar um T trocamos uma purina (G) por uma pirimidina(T), veja que não interessa quem está em cima ou quem está em baixo, temos uma transversão.

GCAACGCACAACGAAAACCCTTAGGGACTGGATTATTTCGTGATCGTTGTAGTTATTGGAAGTACGGGCATCAACCCAGTT || TTATCTGACAAAGAAAGCCGTCAACGGCTGGATAATTTCGCGATCGTGCTGGTTACTGGCGGTACGAGTGTTCCTTTGGGT

Segunda base, trocamos C por T, outra transversão.

GCAACGCACAACGAAAACCCTTAGGGACTGGATTATTTCGTGATCGTTGTAGTTATTGGAAGTACGGGCATCAACCCAGTT ||| TTATCTGACAAAGAAAGCCGTCAACGGCTGGATAATTTCGCGATCGTGCTGGTTACTGGCGGTACGAGTGTTCCTTTGGGT

Na terceira nada acontece, não contamos nada

GCAACGCACAACGAAAACCCTTAGGGACTGGATTATTTCGTGATCGTTGTAGTTATTGGAAGTACGGGCATCAACCCAGTT ||||||||||||||||| TTATCTGACAAAGAAAGCCGTCAACGGCTGGATAATTTCGCGATCGTGCTGGTTACTGGCGGTACGAGTGTTCCTTTGGGT

Na base 17, temos uma transição, troca de purinas, A e G.

Então já deu para ver que temos uma atividade simples, é só percorrer os dois vetores, e quando as bases são diferentes, a gente conta mais um para transição ou transversão, dependendo do caso.

Isso pode ser feito com uma função simples.

1
2
3
4
5
6
7
def razao_tstv(s1, s2):
    transicoes = set([('A', 'G'), ('G', 'A'), ('C', 'T'), ('T', 'C')])
    razao = {True: 0.0, False: 0.0}
    for p in zip(s1, s2):
        if p[0] != p[1]:
            razao[p in transicoes] += 1
    return razao[True] / razao[False]

Veja que a gente define somente as transições, que são menos definições, faz um dicionário com verdadeiro e falso como chave apontando para as respectivas contagens, e faz um loop com a função zip.
Zip é meio estranho, mas ele vai fazer um iterador para tudo que estivermos colocando nele, assim para cada iteração, p vais ser uma lista que vai ter a primeira letra em s1 e s2, como uma lista, na segunda iteração a segunda letra, e assim por diante. Mas o zip vai apenas apontar, e não criar um novo objeto, então se alterássemos algo, se eu não estiver errado, esse algo seria alterado em s1 e s2, mas como eles não são mais usados nesse escopo, não interessa muito.
Dai é fácil, cada iteração que fazemos que começamos la em cima, olhamos se tem diferença entre as bases, se tiver olhamos se é uma transição, senão com certeza é transversão, e depois e só dividir as duas contagens, e como já começamos os números como floats, a divisão já vem como divisão entre floats, não precisamos fazer nenhum tipo de casting.

O resto é fácil, é so abrir o fasta, como sempre fazemos, dar um parse nele e mandar, que é fácil, so tem duas sequências, e mandar para a função.

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
#!/usr/bin/env python
'''
 
'''
 
def razao_tstv(s1, s2):
    transitions = set([('A', 'G'), ('G', 'A'), ('C', 'T'), ('T', 'C')])
    ratio = {True: 0.0, False: 0.0}
    for p in zip(s1, s2):
        if p[0] != p[1]:
            ratio[p in transitions] += 1
    return ratio[True] / ratio[False]
 
arquivo = open("/home/augusto/Downloads/rosalind_tran.txt")
 
sequencia=[]
i=-1
 
for linha in arquivo:
	if linha.find('>') != -1:		
		i+=1
		sequencia.append('')
	else:		
		sequencia[i]=sequencia[i]+linha.rstrip()
 
print razao_tstv(sequencia[0], sequencia[1])

Rosalind – RNA Splicing

A gente já passou DNA para RNA e para aminoácidos, alias, uma fita de dna é copiada em uma fita de rna durante a transcrição, o que eu acho meio confuso, porque um lado do dna que realmente codifica uma proteína, e o outro lado fica la aberto? Bem as coisas são mais complicadas que isso, e por isso sobram mais dúvidas que certezas.

No núcleo, uma enzima chamada RNA polymerase (RNAP) inicia a transcrição quebrando as ligações que ligam as bases de DNA, então criam uma molécula precursora do mRNA,chamada de pre-mRNA, transcrevendo uma das fitas do dna como base, trocando as timinas por uracilas.

Ou seja, a gente pega uma das fitas e vai complementando ela trocando T por U para virar o RNA e ja sabemos que proteína teremos? As coisas não são tão simples.Como RNA é construído baseado na complementariedade, a segunda fita de DNA, chamada fita de codificação (“coding strand”), é quem identifica a nova fita de RNA, exceto pela troca de timina por uracila. Depois que o RNAP criou vários nucleotídeos de RNA, a fita complementar junta de volta a fita original, fechando aquele zíper que é o dna, abrindo para fazer pre-mRNA rapidinho e fechando de novo.

Agora, ainda antes do virar proteína, o pre-RNA é cortado em segmentos menos chamados de introns e exons e para virar proteína, os introns são jogados fora, e os exons grudados juntos sequencialmente para produzir a fita final de mRNA. Esse processo e corta e cola é chamado de splicing, e é facilidade por uma coleção de RNA e proteínas chamadas spliceosome. Como o spliceosome é feito de RNA e proteína é algo do tipo, quem veio primeiro, o ovo ou a galinha e ainda não é bem claro.

Olhando para a fita de dna somente, as partes que geram os exons são a parte codificante do dna. E nesse problema a gente tem que simplesmente codificar uma proteína, levando em conta esse mecanismo.

A nossa entrada vai ser o seguinte:

>Rosalind_10 ATGGTCTACATAGCTGACAAACAGCACGTAGCAATCGGTCGAATCTCGAGAGGCATATGGTCACATGATCGGTCGAGCGTGTTTCAAAGTTTGCGCCTAG >Rosalind_12 ATCGGTCGAA >Rosalind_15 ATCGGTCGAGCGTGT

A primeira sequencia é a dna, e depois teremos uma série de introns, dai temos que pegar esses introns, cortar eles fora, como no spliceing, e ver que proteína será codificada.

Então pegamos o dna

ATGGTCTACATAGCTGACAAACAGCACGTAGCAATCGGTCGAATCTCGAGAGGCATATGGTCACATGATCGGTCGAGCGTGTTTCAAAGTTTGCG

Achamos o primeiro intron

ATCGGTCGAA

Acho ele no dna

ATGGTCTACATAGCTGACAAACAGCACGTAGCAATCGGTCGAATCTCGAGAGGCATATGGTCACATGATCGGTCGAGCGTGTTTCAAAGTTTGCG

E corto ele fora

ATGGTCTACATAGCTGACAAACAGCACGTAGCATCTCGAGAGGCATATGGTCACATGATCGGTCGAGCGTGTTTCAAAGTTTGCG

Próximo intron

>Rosalind_15 ATCGGTCGAGCGTGT

achamos ele

ATGGTCTACATAGCTGACAAACAGCACGTAGCATCTCGAGAGGCATATGGTCACATGATCGGTCGAGCGTGTTTCAAAGTTTGCG

Cortamos fora

ATGGTCTACATAGCTGACAAACAGCACGTAGCATCTCGAGAGGCATATGGTCACATGTTCAAAGTTTGCG

E fazemos isso até tirar todos os introns, depois a gente traduz esse dna para RNA e depois para a sequencia de aminoácidos que vai fazer uma proteína, mas como a gente olha muito mais para o dna, ja facilitaram a vida fazendo a tabela que faz direto o DNA->RNA->Aminoacido.

O código fica bem simples, vou comentar as linhas aqui, mas no final deixo o código completo

#!/usr/bin/env python

Essa linhas é apenas para falar para o bash quem é o interpretador, para poder usar o script direto, como um programa, isso é legal que funciona para o R também.

import re

Usamos a biblioteca de expressões regulares.

RNA_CODON_TABLE = {
'ATA':'I', 'ATC':'I', 'ATT':'I', 'ATG':'M',
'ACA':'T', 'ACC':'T', 'ACG':'T', 'ACT':'T',
'AAC':'N', 'AAT':'N', 'AAA':'K', 'AAG':'K',
'AGC':'S', 'AGT':'S', 'AGA':'R', 'AGG':'R',
'CTA':'L', 'CTC':'L', 'CTG':'L', 'CTT':'L',
'CCA':'P', 'CCC':'P', 'CCG':'P', 'CCT':'P',
'CAC':'H', 'CAT':'H', 'CAA':'Q', 'CAG':'Q',
'CGA':'R', 'CGC':'R', 'CGG':'R', 'CGT':'R',
'GTA':'V', 'GTC':'V', 'GTG':'V', 'GTT':'V',
'GCA':'A', 'GCC':'A', 'GCG':'A', 'GCT':'A',
'GAC':'D', 'GAT':'D', 'GAA':'E', 'GAG':'E',
'GGA':'G', 'GGC':'G', 'GGG':'G', 'GGT':'G',
'TCA':'S', 'TCC':'S', 'TCG':'S', 'TCT':'S',
'TTC':'F', 'TTT':'F', 'TTA':'L', 'TTG':'L',
'TAC':'Y', 'TAT':'Y', 'TAA':'_', 'TAG':'_',
'TGC':'C', 'TGT':'C', 'TGA':'_', 'TGG':'W',
}

Ai temos nosso dicionario de dna para aminoácido

def protein_string(mrna):
    result = ''
 
    for i in range(0, len(mrna), 3):
        symbol = RNA_CODON_TABLE[mrna[i:i+3]]
        if symbol == '_':
            break
        result += symbol
 
    return result

Criamos uma função que recebe um dna e o codifica para a proteína, que basicamente vai ler um codon, o string que é o dna de 3 em 3 letras, por isso o range é 3, vai pegar no symbol quais 3 letras são essas e usar como a chave do dicionario que retorna o aminoácido e vai concatenando no resultado, lembrando que o “_” está como o stop codon, então quanto ele aparece, nos quebramos o laço, ai o retorno é a string que representa a sequencia de aminoácidos.

arquivo = open("/home/augusto/Downloads/rosalind_splc.txt")

A partir dai é so abrir o arquivo, que a gente faz download do rosalind

sequencia=[]
i=-1

Criar uma lista para receber o arquivo fasta, esse -1 é apenar para sempre começar do começo no find.

for linha in arquivo:
	if linha.find('>') != -1:		
		i+=1
		sequencia.append('')
	else:		
		sequencia[i]=sequencia[i]+linha.rstrip()

veja que vou iterando linha por linha do fasta, até achar as linhas que tem o nome da sequencia, e vou juntando os dados abaixo dessa linha, começando com um nada, para criar o sequencia, e depois concatenando linhas posteriores.

principal=sequencia[0]
sequencia.pop(0)

Depois eu separo a sequencia principal, isso não é necessário na verdade, so para facilitar, bastaria usar os índices de maneira correta, mas de toda forma…

for seq in sequencia:
	principal=re.sub(seq, '', principal)

Agora eu vou apagando intron por intron usando o sub das expressões regulares

print protein_string(principal)

E imprimir o aminoacido do que sobrar.

Bem é isso ai, acho que esse problema é mais para ilustrar o processo, que é legal e ver que não basta olhar pro dna diretamente para ver as coisas. 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:

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
#!/usr/bin/env python
'''
 
'''
 
import re
 
RNA_CODON_TABLE = {
'ATA':'I', 'ATC':'I', 'ATT':'I', 'ATG':'M',
'ACA':'T', 'ACC':'T', 'ACG':'T', 'ACT':'T',
'AAC':'N', 'AAT':'N', 'AAA':'K', 'AAG':'K',
'AGC':'S', 'AGT':'S', 'AGA':'R', 'AGG':'R',
'CTA':'L', 'CTC':'L', 'CTG':'L', 'CTT':'L',
'CCA':'P', 'CCC':'P', 'CCG':'P', 'CCT':'P',
'CAC':'H', 'CAT':'H', 'CAA':'Q', 'CAG':'Q',
'CGA':'R', 'CGC':'R', 'CGG':'R', 'CGT':'R',
'GTA':'V', 'GTC':'V', 'GTG':'V', 'GTT':'V',
'GCA':'A', 'GCC':'A', 'GCG':'A', 'GCT':'A',
'GAC':'D', 'GAT':'D', 'GAA':'E', 'GAG':'E',
'GGA':'G', 'GGC':'G', 'GGG':'G', 'GGT':'G',
'TCA':'S', 'TCC':'S', 'TCG':'S', 'TCT':'S',
'TTC':'F', 'TTT':'F', 'TTA':'L', 'TTG':'L',
'TAC':'Y', 'TAT':'Y', 'TAA':'_', 'TAG':'_',
'TGC':'C', 'TGT':'C', 'TGA':'_', 'TGG':'W',
}
 
 
def protein_string(mrna):
    result = ''
 
    for i in range(0, len(mrna), 3):
        symbol = RNA_CODON_TABLE[mrna[i:i+3]]
        if symbol == '_':
            break
        result += symbol
 
    return result
 
 
 
arquivo = open("/home/augusto/Downloads/rosalind_splc.txt")
 
sequencia=[]
i=-1
 
for linha in arquivo:
	if linha.find('>') != -1:		
		i+=1
		sequencia.append('')
	else:		
		sequencia[i]=sequencia[i]+linha.rstrip()
 
 
principal=sequencia[0]
sequencia.pop(0)
 
 
for seq in sequencia:
	principal=re.sub(seq, '', principal)
 
print protein_string(principal)

Rosalind – Introduction to Random Strings

Nesse problema, a gente tem que simular strings baseado na quantidade de GC dela. Uma das relevâncias do conteúdo de GC, se eu não estiver falando muita bobeira, é que G e C se ligam fazendo 3 pontes de hidrogênio, o que confere estabilidade a Molécula de DNA, mas tem muitos motivos diversos para medir a quantidade de GC.

Ok, mas sem devagar do problema, nossa temos uma sequência de DNA e um vetor de valores entre 0 e 1, que são proporções do conteúdo de GC de uma molécula de DNA. E temos que calcular qual a probabilidade de, ao acaso, a sequência acima ter sido gerada a partir desse conteúdo de GC.
Isso é muito parecido com sortear caras e coroas a partir de um likelihood (aqui, um post de 2012, esse é velho hehe).

Nossa entrada é:

ACGATACAA 0.129 0.287 0.423 0.476 0.641 0.742 0.783

Como são quatro bases, temos 4 possibilidades para cada posição, então temos 4\cdot4\cdot4\cdot4\cdot4\cdot4\cdot4\cdot4\cdot4=4^9=262144 possibilidades de moléculas de DNA. Agora suponto que a gente atribua uma probabilidade p de 0.129 para sair C ou G e 1-p, 0.871 para A,T, qual a chance de sair uma molécula com a mesma quantidade de GC que estamos vendo.

Primeiro vamos colocar esses dados no python

1
2
3
dna='ACGATACAA'
prob_gc='0.129 0.287 0.423 0.476 0.641 0.742 0.783'       
prob_gc = map(float, prob_gc.split())

A função map é para transformar string da leitura em um tipo que nos interessa, nesse caso float. Dai contamos a quantidade de GC

1
2
3
4
5
6
7
conteudo_gc = 0
conteudo_at = 0
for codon in dna:
	if codon in ['C', 'G']:
		conteudo_gc += 1
        else:
                conteudo_at += 1

Agora, pensando em uma quantidade apenas, para cada base, ela pode ser GC ou não, então é a chance da base ser GC ou não, coincidir com o que esperamos e a próxima também, e a próxima também. E em probabilidade, um ‘e’ significa que multiplicamos probabilidades.

1
2
3
4
5
6
7
gc=0.129
total=1.0
for base in dna:
        if base in ["G","C"]:
                total*=  gc*0.5
        else:
                total*=  (1-gc)*0.5

Mas veja que como falamos de uma base no conjunto ATCG, estamos falando da metade da probabilidade, dae multiplicamos por meio, porque a outra metade não ia resolver, certo? Mas então somamos tudo isso e é só tirar o log10

1
print log10(total)

Mas ai entramos em um problema de computação, que é o underflow, veja que a gente ta multiplicando e multiplicando um número, e ele ta ficando pequeno e pequeno até não ser possível mais representá-lo, e é ai que o log começa a valer a pena.

Veja que logs tem a propriedade de

log(x\cdot y) = log(x)+log(y)

E somando, a gente não vai encolhendo os números, então resolve o problema de precisão. Facilita bem pelo menos. Então podemos modificar o código acima para

1
2
3
4
5
6
7
8
total=0.0
for base in dna:
        if base in ["G","C"]:
                total+=  log10(gc*0.5)
        else:
                total+=  log10((1-gc)*0.5)
 
print total

Agora é só fazer isso, conteúdo de GC por conteúdo de GC

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
from math import log10
 
with open('/home/augusto/Downloads/rosalind_prob.txt') as input_data:
    dna, prob_gc = input_data.readlines()
 
prob_gc = map(float, prob_gc.split())
 
#Conta o número de GC e AT
conteudo_gc = 0
conteudo_at = 0
for codon in dna:
	if codon in ['C', 'G']:
		conteudo_gc += 1
        else:
                conteudo_at += 1
 
prob_lista  = []
for prob in prob_gc:        
	prob_final = conteudo_gc*log10(0.5*prob) + (len(dna)-conteudo_gc) * log10(0.5*(1-prob))
	prob_lista.append(str(prob_final))
 
print ' '.join(prob_lista)

Outra coisa, é que como estamos falando de quantidade de GC, não interessa a ordem das bases, e por último, temos a intuição, que o maior likelihhod vai ser no conteúdo de GC mais próximo do valor encontrado na sequência mesmo. E o join no string final é só para não imprimir a resposta desorganizada, como quando a gente imprime uma lista no python.

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:
Rosalind

Rosalind – Partial Permutations

Um exercício do Rosalind, para matar a saudade de fazer probleminhas.

Chama-se partial permutations.

Quando eu li esse nome, eu fui no google e fiz uma busca por partial permutations. E tem uma página do Wikipédia com esse título aqui. Veja que essa não é exatamente a solução. O que o problema realmente quer é o que está no fim dessa página, que é fazer os Restricted partial permutations.

O que a gente quer, por exemplo, é de um conjunto de 1 a 8, combinando de 3 em 3, quantas possibilidades temos, sem repetir os números.

Então alguns exemplos de conjuntos seriam:

(1,2,3) (1,2,4) (1,2,5) ... (7,8,9)

Então P(n,k) é quantos casos desses podemos escrever, onde n é as coisas que usamos para formar os conjuntos e k é quanto “espaço” temos. Isso serve para pensar por exemplo, quantos possibilidades temos de organizar os genes no genoma.

Bem, isso tem uma formula fechada, so consultar aqui

E assim perdeu a graça, em python, a função fatorial tem na biblioteca math, então:

1
2
3
4
import math
n=21
k=7
print (math.factorial(n)/math.factorial(n-k))%1000000

E temos

Python 2.7.10 (default, Oct 14 2015, 16:09:02) [GCC 5.2.1 20151010] on linux2 Type "help", "copyright", "credits" or "license" for more information. >>> >>> >>> >>> 51200

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:
Rosalind – Partial Permutations

Grafos – 101

Grafos, isso é algo que quando a gente choaqualha muito esse blog, começa a aparecer aqui e ali. Bem grafos são uma representação matemática que pode ser utilizado para resolver muitos problemas da no mundo real.

Em biologia a gente constantemente vê grafos em árvores filogenéticas, redes de interações, entre outras coisas, mas um monte de problemas de computação são modelados como grafos.

Mas o que é um grafo?

A definição que tem no livro (a maioria das figura desse post peguei desse livro, ta na referência) do Robert Sedgewick é

Definition. A graph is a set of vertices and a collection of edges that each connect a pair of vertices.

Simples assim, dois conjuntos, um de vértices e outro de arestas, vértices são as bolinhas e as linhas.

Grafos apareceram no mundo após um artigo publicado pelo Euler, o mesmo cara que ezinho do exponencial, ele estudava o problema das pontes de Königsberg, e para resolver esse problema, ele começou pensar no problema como arestas e vértices, e assim começou a teoria dos grafos.

Bem pra começar, vamos so tentar criar uma classe em python, para guardar grafos, para começar a tentar fazer alguma coisa, para entender melhor como os grafos funcionam e tentar usar eles melhores em biologia.

Existem varias formas de representar grafos, por exemplo a gente pode guardar as arestas e dizer quantos vértices temos.

O exemplo que tem no livro de algorítimos do Sedgewick é assim:

13 13 0 5 4 3 0 1 9 12 6 4 5 4 0 2 11 12 9 10 0 6 7 8 9 11 5 3

O primeiro número, é número de vértices, o segundo é o numero de arestas, nesse caso, depois disso estamos listando todas as arestas. Alguns detalhes são que temos 13 vértices, mas veja que eles começam no zero, assim o maior vértice é o 12. Também não estamos dando nomes pra eles ainda, mas isso é simples, futuramente é só usar um dicionario. E estamos tratando de um tipo especifico de grafo, os não não direcionais, ou seja, a gente não se preocupa com a direção da aresta.

Mas então. Veja que uma representação visual seria assim:

01

Veja que todas as arestas ali em cima estão nessa figura. Agora a primeira coisa que complica, é que não existe uma única representação gráfica possível, veja que a gente pode desenhar o mesmo grafo assim:

02

Veja que esse grafo é a mesma coisa, e o fato de poder desenhar ele de mais de uma forma, é a mesma coisa que confunde em árvores filogenéticas por exemplo, que podem ser “rotacionadas” quando aos seus cluster.

Certo, mas como podemos guardar isso no computador? Qual estrutura de dados pode ser usada. Uma forma é simplesmente guardar o conjunto de arestas como visto ali em cima, mas isso é ruim, porque para extrair informação do grafo depois, tudo fica muito difícil, pense, como contar quantos vértices temos, ou avaliar o grau de cada vértice, isso se torna difícil, ja que constantemente temos que olhar todos o vértices para responder essas perguntas.

Uma segunda forma é guardar os dados em matrizes de adjacência. Essa matriz tem o tamanho do número de vértices ao quadrado e e colocamos um 1 na célula ij quando existe uma aresta de i para j. Mas essa também uma uma solução ruim, porque o tamanho da matriz cresce ao quadrado do número de vértices, ou seja, a cada novo vértice, precisamos de muito mais espaço para guardar o grafo, e logo isso se torna proibitivo.

Uma forma mais amena, que deixa os dados de forma a serem acessado mais rapidamente e não ocupa tanto espaço é uma lista de listas, um vetor de ponteiros no java como no livro, mas aqui podemos usar a lista do python.

03

Basicamente a gente faz uma lista que tem o tamanho do número de vértices, e cada posição dessa lista guarda outra lista, que é aonde podemos chegar a a partir desse vértice. Assim, podemos começar nossa classe no python, já que já temos como guardar os dados.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
class Graph:
	#Contrutor para grafo vazio ou com um argumento, o numero de vertices
	def __init__(self,V=None):
 
		if V is not None:
			self.V=V
			self.E=0
			self.adj_Lista()
 
		else:
			self.V=int(raw_input())
			self.E=0
			self.adj_Lista()		
 
			e=int(raw_input())
 
			for i in range(e):
				entrada = raw_input()
				entrada = entrada.split(' ')
				v=int(entrada[0])
				w=int(entrada[1])
				self.addEdge(v,w)

Começamos nossa classe definindo um construtor pra ela. Algo que é um pouco diferente do java, c++ e outras linguagens, é que no python, a gente não precisa declarar quais os atributos da classe, a gente declara isso logo dentro do construtor, gerando o atributo.
Outra coisa é que no python a gente só faz 1 construtor, não podemos ter vários com assinaturas diferentes como no c++ ou java. Aqui eu quis fazer 2 construtores, um para iniciar um grafo vazio, que no caso deveríamos apenas dizer o número de vértices que ele tem, e outro em que leriamos o arquivo de exemplo do livro, para ler o arquivo de exemplo, não mandaríamos nenhum argumento, mas leriamos do terminal, ou de um arquivo mandado diretamente, ou seja se V, que é o número de vértices é fornecido, fazemos um grafo vazio com V vértices, senão a gente le o padrão de input do livro.
Não vou explicar tudo, senão ficarei até amanha falando de código aqui, mas veja que eu fiz uma função chamada adj_Lista(), ela é que cria a lista de listas vazias

1
2
3
4
	def adj_Lista(self):
		self.adj=[]
		for _ in range(self.V):
			self.adj.append([])

Aqui a gente so cria uma lista, e para o número de vértices do grafo, criamos listas vazias dentro dessa lista, essa é a estrutura para guarda o grafo em si.

Assim, adicionar vértices é fácil

1
2
3
	def addEdge(self,v,w):
		self.adj[v].append(w)
		self.E+=1

Basta colocar o vértice o valor w, pra onde a aresta vai, na posição da lista v que é onde o vértice ta saindo.

Agora para tentar fazer alguma coisa, podemos ver como é o algorítimo usado para testar se o grafo é conectado, o que é ser conectado para um grafo? Conectado é quando a partir de qualquer vértice do grafo, podemos chegar a qualquer outro vértice, ou seja, não tem vértices isolados.

Eu fiz mais algumas funções, que estão disponíveis no código, particularmente uma para gerar uma visualização do grafo, onde eu uso o igraph, que é o mesmo igraph do R, ele tem uma implementação em python também, e a função deep first search para testar se um grafo é conectado.

Primeiro vamos ver dois exemplos do uso aqui, para ver como ficam os grafos conectados e não conectados.

Primeiro o não conectado:

Captura de tela de 2015-11-04 20:04:54

E o conectado:

Captura de tela de 2015-11-04 20:05:24

Então o programa recebe os dados de um arquivo txt, monta o grafo e usa o dfs, deep search first para avaliar se o grafo é conectado ou não.

Pra facilitar entender como isso é feito, vamos olhar como é o código.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
	#o grafo e conectado
	def connected(self):
		self.visitado = [False]*self.V
		self.count=0
 
		self.dfs(0)
 
		if self.count == self.V:
			return "Conectado"
		else:
			return "Nao conectado"
 
	#deep first searsh para a funcao conectado
	def dfs(self, origem):
		self.visitado[origem] = True
		self.count+=1
 
		for i in self.adj[origem]:
			if self.visitado[i] == False:
				self.dfs(i)

Ele está dividido em duas funções:
A primeira:

1
2
3
4
5
6
7
8
9
10
11
	#o grafo e conectado
	def connected(self):
		self.visitado = [False]*self.V
		self.count=0
 
		self.dfs(0)
 
		if self.count == self.V:
			return "Conectado"
		else:
			return "Nao conectado"

Essa função faz o seguinte, veja que vc criou um objeto dessa classe, logo ele é um grafo, então ela é um método da classe que responde se o grafo é conectado ou não, depois podemos trocar isso para uma resposta do tipo verdadeiro ou falso.

Veja que ela vai criar uma lista do tamanho do número de vértices, com falso em todas as posições, e uma contagem, iniciando em zero. Feito isso a gente chama o método dfs.

1
2
3
4
5
6
7
8
	#deep first searsh para a funcao conectado
	def dfs(self, origem):
		self.visitado[origem] = True
		self.count+=1
 
		for i in self.adj[origem]:
			if self.visitado[i] == False:
				self.dfs(i)

No dfs, ela é chamado com um valor, o nome do vértice de origem. Primeiro registramos então que visitamos esse vértice e acrescentamos 1 a nossa contagem de vértices visitados, dai, para todos os vértices que podemos chegar a partir desse vértice e que ainda não foram visitados, chamamos recursivamente o método dfs.
E aqui está a mágica, dfs vai so vai ser chamada para quem não foi visitado ainda, ou seja partindo de um ponto, aqui sempre do vértice 0, a gente vai exatamente uma vez em cada vértice, guardando na memória se ja fomos ou não nos vértices, o mesmo algorítimo que teseu usou para matar o minotauro, so que ao invés de lista de verdadeiro ou falso, ele usava linha para marcar onde já passou.

Eu acho que deveria acabar por aqui esse post, que faz uma semana que está nos drafts e não é publicado nunca, porque não sei o que exatamente estou fazendo escrevendo ele. Na verdade, eu gostaria de aprender mais sobre grafos e nada melhor que ir testando os algorítimos. A classe em si, e como guardar o grafo é demasiadamente complicado já, e não quis ficar picando em mil posts sobre grafos, porque as coisas começam a ficar legal quando a gente começa a ver algorítimos que fazem alguma coisa, como testar se o grafo é conectado ou não.

No meu repositório recologia, vai ter pasta chamada python, e dentro dela tem outra pasta chamada grafos, o codigo fonte vai estar todo la, além da classe de grafos, muitas das outras funcionalidades são implementadas em outras classes no exemplos do Sedgewick, e espero passar tudo para python, um dia quem sabe.

Mas é isso, fica aqui esse post publicado, porque ele ja encheu saco na lista de draft, é legal para ver o código do python, é bem legal usar orientação a objetos, e bem diferente python de linguagens como c++ e java, logo so fazendo classes para aprender. Fico por aqui e espero deslanchar meu post sobre ggplot2 logo.

Referência:
Robert Sedgewick and Kevin Wayne 2011 Algorithms, 4th Edition

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
import igraph
 
class Graph:
 
	#Contrutor para grafo vazio ou com um argumento, o numero de vertices
	def __init__(self,V=None):
 
		if V is not None:
			self.V=V
			self.E=0
			self.adj_Lista()
 
		else:
			self.V=int(raw_input())
			self.E=0
			self.adj_Lista()		
 
			e=int(raw_input())
 
			for i in range(e):
				entrada = raw_input()
				entrada = entrada.split(' ')
				v=int(entrada[0])
				w=int(entrada[1])
				self.addEdge(v,w)
 
	#Funcao que cria a lista de adjacencia, usada no construtor
	def adj_Lista(self):
		self.adj=[]
		for _ in range(self.V):
			self.adj.append([])
 
	#Retorna o numero de vertices
	def vertices(self):
		return self.V
	#retorna o numero de arestas
	def arestas(self):
		return self.E
 
	#metodo para adicionar arestas
	def addEdge(self,v,w):
		self.adj[v].append(w)
		self.E+=1
 
	#representando como string
	def __str__(self):
		string= str(self.V) + " vertices e " + str(self.E) + " arestas"
		return string
 
	#imprimindo arestas
	def imprime_arestas(self):
		v=0
		n=0
		i=0
 
		while v < self.V:
			n = len(self.adj[v])
			if n>0:
				while i < n:
					print v,self.adj[v][i]
					i+=1
			v+=1
			i=0
 
	#grau de V
	def grau_de_v(self,v):
		return len(self.adj[v])
 
	#maior grau
	def maior_grau(self):
		lista = []
		for i in range(self.V):
			lista.append(self.grau_de_v(i))
		return max(lista)
 
	#grau medio
	def grau_medio(self):
		return 2*self.E/float(self.V)
 
	#conta auto-loops
	def numero_auto_loops(self):
		pass
 
	#o grafo e conectado
	def connected(self):
		self.visitado = [False]*self.V
		self.count=0
 
		self.dfs(0)
 
		if self.count == self.V:
			return "Conectado"
		else:
			return "Nao conectado"
 
	#deep first searsh para a funcao conectado
	def dfs(self, origem):
		self.visitado[origem] = True
		self.count+=1
 
		for i in self.adj[origem]:
			if self.visitado[i] == False:
				self.dfs(i)
 
	def plot(self):
		g = igraph.Graph()
		g.add_vertices(self.V)
 
		v=0
		n=0
		i=0
 
		while v < self.V:
			n = len(self.adj[v])
			if n>0:
				while i < n:
					g.add_edges((v,self.adj[v][i]))
					i+=1
			v+=1
			i=0
 
		layout = g.layout("kk")
		igraph.plot(g, layout = layout,vertex_label=range(self.V),vertex_size=30)