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

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