Rosalind – Introduction to Set Operations

Como faz tempo que não coloco nada aqui no blog, tem um problema que resolvi a algum tempo do Rosalind, bem simples que nem coloquei aqui.

O problema é bem simples, basta realizar alguma operações de conjunto, apenas para ver como funciona.

Então a entrada é assim:

10
{1, 2, 3, 4, 5}
{2, 8, 5, 10}

O primeiro valor da primeira linha, 10 quer dizer que o espaço inteiro dos elementos é de 1 até 10. As duas linhas depois são dois conjuntos, com os elementos entre as chaves e separados por virgulas.

A partir dessa entrada a gente tem que produzir 6 conjuntos.

A \cup B , A \cap B , A - B , B - A , A^c e B^c

União e interseção e diferenças a gente usa somente os conjuntos enquanto para o complemento, a gente precisa saber todos os elementos do espaço dos conjuntos.

Podemos definir uma função para fazer as operações a partir da entrada

1
2
3
4
5
6
7
8
9
def conjuntos(n, a, b):
    return [
        a | b,
        a & b,
        a - b,
        b - a,
        set(range(1, n + 1)) - a,
        set(range(1, n + 1)) - b
    ]

E dai é so fazer o parse da entrada, que é bem simples, são sempre 3 linhas, so pegar o primeiro elemento e depois os elementos dentro das chaves, e temos que arrumar a saída também, para sair na formação que a pergunta pede. Eu usei sys.stdout.write pois achei mais fácil para formatar com ele do que com a função print do python. Essa função é o motivo do import sys no começo do arquivo.

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 e até mais.

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
import sys
 
def conjuntos(n, a, b):
    return [
        a | b,
        a & b,
        a - b,
        b - a,
        set(range(1, n + 1)) - a,
        set(range(1, n + 1)) - b
    ]
 
 
'''
n=10
a=set([1, 2, 3, 4, 5])
b=set([2, 8, 5, 10])
 
saida=conjuntos(n,a,b)
'''
dados=open("/home/augusto/Downloads/rosalind_seto.txt",'r').read()
 
n=int(dados[0:dados.find('\n')])
 
chave1=dados.find('{')
chave2=dados.find('{',chave1+1)
 
a=dados[chave1+1:dados.find('}')]
a=a.split(', ')
a=map(int,a)
a=set(a)
 
b=dados[chave2+1:dados.find('}',chave2)]
b=b.split(', ')
b=map(int,b)
b=set(b)
 
saida=conjuntos(n,a,b)
 
for conjunto in saida:
    conjunto=list(conjunto)
    tamanho=len(conjunto)
    print "{",
    for i in range(0,tamanho):
        if i == tamanho-1:
            sys.stdout.write("%d" % conjunto[i])
        else:
            sys.stdout.write("%d, " % conjunto[i])
    print"}"

Programação multi-thread em python

Esse talvez possa parecer um post meio maluco, para o que comumente discutimos aqui nesse blog, mas é muito interessante e curioso, para saber o que está acontecendo sobre o capo do carro, ou do computador.

Uma das coisas mais diferentes que já vi, foi quando fiz a disciplina de sistemas operacionais, isso porque tudo que eu achava que sabia sobre computador, vi que estava errado. Nem tudo é tão simples como um script continuo que fazemos aqui no R, e uma prova disso é a programação Multi-thread.

thread

Vamos pensar na figura acima, tomando como exemplo os MCMC que vemos comumente aqui, veja alguma exemplos de implementação aqui de um gibbs samples e aqui do mais simples metropolis hastings.

Nesses dois casos fizemos algo como em single-threaded, isso porque fizemos uma cadeia e fomos realizando vários cálculos e salvando resultados. Acontece que, todas as analises de dados que fazemos, não usamos apenas uma cadeia, usamos várias cadeias, no último post mesmo, usando o openbugs, usamos três cadeias, que podem ser processadas separadamente, podem ser multithreaded.
Os valores das amostras, parâmetros, modelos, tudo é igual, então não compensa gravar isso varias vezes na memoria, é melhor compartilhar, mas para uma cadeia de markov, a gente sempre precisa do valor anterior para continuar ela, então as cadeias não precisam dividir os cálculos que estão fazendo entre si, apenas os dados e parâmetros e modelos, agora olha no modelinho ali em cima, não é exatamente isso que diz respeito a programação multi-thread?
Veja que com multi-thread, cada thread tem sua pilha, seus registradores, seus dados particulares, que é sua cadeia de markov que está computando, mas também divide dados, que são as amostras, modelo, parâmetros entre todas as threads.

E qual a relevância disso? Ai que entre aqueles processadores de muitos núcleos, sem isso, sem paralelizar, dividir a computação, a gente teria que fazer uma cadeia no processador, terminar, fazer outra, terminar, pense que leva um minuto por cadeia, com três então levamos três minutos para fazer tudo, mas se fizermos algo multi-thread, da para mandar cada cadeia para cada núcleo do processador (mais ou menos isso, simplificando, sei que não é assim), como cada um leva 1 minuto e da para fazer em paralelo, ao invés de esperar 3 minutos pelo resultado, esperamos um.

image0021225707175745

Podemos fazer um exemplo em python que é bem simples, usando a biblioteca threading.

Para esse exemplo simples, vamos criar uma classe bem bobinha

1
2
3
4
5
6
7
8
9
10
class minhaThread (threading.Thread):
    def __init__(self, threadID, nome, contador):
        threading.Thread.__init__(self)
        self.threadID = threadID
        self.nome = nome
        self.contador = contador
    def run(self):
        print "Iniciando thread %s com %d processos" % (self.name,self.contador)
        processo(self.nome, self.contador)
        print "Finalizando " + self.nome

A gente ja falou de orientação a objetos em R aqui e aqui, uma das coisas que a orientação a objetos permite é a herança, que é herdar tudo que a classe pai tem, aqui estamos herdando a classe Thread definida na biblioteca threading, isso é o que esta escrito aqui:

1
class minhaThread (threading.Thread):

Por isso o parenteses, lembro que ja usei classe em algum problema do rosalind aqui, mas não lembro em qual post agora, mas beleza, o resto é bem besta a classe vai ter 3 atributos, um ID, um nome e um contador, que seria algo como os parâmetros para fazer um mcmc, vamos abstrair de uma atividade útil agora, vamos apenas pensar em termos de computação.

1
2
3
4
5
    def __init__(self, threadID, nome, contador):
        threading.Thread.__init__(self)
        self.threadID = threadID
        self.nome = nome
        self.contador = contador

Então nos construtor da thread definimos esses atributos e ela vai ter apenas um método, uma função chamada run.

1
2
3
4
    def run(self):
        print "Iniciando thread %s com %d processos" % (self.name,self.contador)
        processo(self.nome, self.contador)
        print "Finalizando " + self.nome

Tem que chamar run, por causa que é assim que funciona o pacote, mas o run vai so imprimir que estamos começando, processar algo, pensando que processar pode ser qualquer coisa e depois avisar que terminamos.

Aqui o processo definimos assim

1
2
3
4
def processo(nome, contador):
    while contador:
        print "Thread %s fazendo o processo %d" % (nome, contador)
        contador -= 1

Ele é uma função externa, não é da classe, mas ele apenas imprimi uma mensagem, falando que esta processando, e diminui o contador em 1, até zerar, e terminar.

Depois disso está tudo pronto, so precisamos criar nossas threads

1
2
3
# Criando as threads
thread1 = minhaThread(1, "Alice", 8)
thread2 = minhaThread(2, "Bob", 8)

Criamos elas, simples assim, aqui no caso são duas, chamadas de alice e bob, hehe, alice e bob são exemplos mais comuns em física, criptografia e teoria dos jogos, ao invés de falar A e B, os caras falam Alice e Bob, para ficar mais fácil de entender.

Ok, mas apos criado, nada foi processado ainda, dai então começamos elas.

1
2
3
# Comecando novas Threads
thread1.start()
thread2.start()

O resto do programa é meramente para fazer a thread mãe esperar alice e bob para terminar, não convêm explicar agora, mas essa parte so faz isso, esperar alice e bob terminarem

1
2
3
4
5
6
threads = []
threads.append(thread1)
threads.append(thread2)
 
for t in threads:
    t.join()

Agora é legal quando a gente roda esse programa.

As vezes ele roda uma thread seguida da outra

Screenshot from 2015-07-23 16:21:42

E as vezes tudo vira uma loucura.

Screenshot from 2015-07-23 16:23:00

Isso porque uma vez que você inicio o processamento, o sistema operacional que vai escalonar as threads, dai a gente não manda mais, só o sistema operacional, então se o processador ta liberado, ele manda elas para la para fazer conta, mas as vezes ta ocupado, então a thread tem que esperar sua vez, agora uma thread pode ir para um núcleo e ir até o final la, enquanto a outra entrou num outro e teve que esperar alguém terminar de processar, como no segundo exemplo, Alice começou, o Bob esperou um pouco, Alice ja tinha processado 2 vezes, ai o Bob começou a trabalhar, mas terminou depois.

Um exemplo bem besta, mas não é preciso saber fazer muitas contas complexas para imaginar que isso é melhor que primeiro alice fazer todo o trabalho dela e depois bob fazer todo o trabalho dela. Pena que isso é extremamente complexo, ainda mais quando a gente tem que compartilhar algum dado e mexer nele, porque ai a vez de quem mexeu pode ser importante, mas isso é outra historia.

Claro que programação não é a meta de muitas pessoas que visitam aqui, mas ter uma ideia que isso existe, pode ajudar a entender comentários de manuais de programa entre outras coisas.

Por exemplo, um programa comum em biologia é o mrBayes, basta olhar o manual que você vai notar alguma notas de instalação, veja essa parte aqui:

instructions for compiling and running the parallel version of the program

Ou seja, se você instala as coisas so dando dois cliques, não le o manual, pode estar perdendo essa possibilidade de usar multi-thread e algo que pode rodar em uma hora vai levar um dia. Ja que ainda no exemplo do mrBayes, da para usar até o processador da placa de video para ajudar a fazer cálculos e agilizar o processamento. E pode abrir o olho, que tudo que você ler MCMC é um forte candidato a usar multi-thread.

O R em si não é multi-thread, mas existem vários pacotes que tentam incorporar isso, como visto no taskview de high performance. Além de que o R possui interfaces que ligam eles a muitos outros programas, como é o caso do Openbugs mesmos, que os calculos são feitos neles, mas depois devolvemos os dados no R para fazer figuras e outras analises, mas temos muitos outros exemplo.

Apesar que na maioria das vezes isso não é muito relevante, para fazer regressões e muitos modelos, tudo é tão rápido que a gente nem sente diferença, mas algumas coisas mais avançadas (talvez “avançando” seja um termo ruim aqui, coisas que precisam fazer mais contas, comum em evolução), isso começa a ser importante.

Bem é isso ai, a primeira vez que ouvi sobre esse negocio de multi-thread, a primeira coisa que veio na minha cara é MCMC, acho que ter ouvido falar de MCMC antes de fazer Sistemas Operacionais foi um diferencial para tornar tudo mais interessante e belo, mas tive um excelente professor em sistemas operacionais também para ajudar. O script vai estar la no repositório recologia, tem mais alguma exemplos da minha aula de sistemas operacionais no github em outro repositório aqui mas em linguagem C, e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail e até mais.

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
import threading
 
class minhaThread (threading.Thread):
    def __init__(self, threadID, nome, contador):
        threading.Thread.__init__(self)
        self.threadID = threadID
        self.nome = nome
        self.contador = contador
    def run(self):
        print "Iniciando thread %s com %d processos" % (self.name,self.contador)
        processo(self.nome, self.contador)
        print "Finalizando " + self.nome
 
def processo(nome, contador):
    while contador:
        print "Thread %s fazendo o processo %d" % (nome, contador)
        contador -= 1
 
# Criando as threads
thread1 = minhaThread(1, "Alice", 8)
thread2 = minhaThread(2, "Bob", 8)
 
# Comecando novas Threads
thread1.start()
thread2.start()
 
threads = []
threads.append(thread1)
threads.append(thread2)
 
for t in threads:
    t.join()
 
print "Saindo da main"

Problema – Rosalind – Completing a Tree – TREE

Ja faz um tempo que não olhava mais o Rosalind, aqui vai mais um problema.

Nesse problema a gente tem que completar um grafo do tipo arvore.

Bem, a primeira pergunta é o que diabos é um grafo, um grafo é um conjunto de pontos (vértices ou node em inglês) e suas ligações (as arestas ou edges em inglês).

Existe uma disciplina chamada teoria dos grafos que estuda propriedades e características de grafos. Dentro dessa teoria eles são divididos em vários tipos, que apresentam características em comum, um desses tipos é o tipo arvore ou Tree citado no nome do problema.

O grafo de árvore é um grafo conexo (existe caminho entre quaisquer dois de seus vértices) e acíclico (não possui ciclos). Basicamente esses caras são as arvores filogeneticas.
Ele é conexo porque aceitamos, na hora de reconstruir uma arvore qualquer, cada par de especies tem pelo menos um ancestral comum, e as especies, que aqui são as folhas, não se ligam umas nas outras, folhas no nosso caso são as espécies. Para ter uma ideia de uma olhada aqui.

Mas após ver tudo isso, é so olhar as propriedades dos grafos tipo arvore para matar esse problema, todo grafo do tipo arvore terá n-1 arestas, sendo n o número de folhas ou espécies.

Agora o problema fornece as seguintes informações

10 1 2 2 8 4 10 5 9 6 10 7 9

10 é o número vértices, e essas duplas de números são as arestas.

Usando o R e o pacote igraph, a gente pode tentar construir esse grafo de curiosidade da seguinte forma

library(igraph)
plot(graph.formula( 1--2, 2--8, 4--10, 5--9,6--10,7--9,3 ))

Eu coloquei o 3 ali sozinho no final, porque ele não aparece na lista de arestas, então ele não apareceria.

01

Mas é isso ai. Um monte de partezinha perdidas. Para ele ser uma árvore, ele precisa n-1 arestas, olha ai temos 7 linhas, então ta faltando 3 linhas, então se adicionarmos essas linhas respeitando as outras propriedades do grafo para ser uma uma árvore, temos 7 arestas já, 10-1 da 9, temos que ter nove arestas, então ta faltando 3, adicionamos elas e ficamos com o seguinte.

plot(graph.formula( 1--2, 2--8, 4--10, 5--9,6--10,7--9,2--3,10--3,9--3 ))

02

Agora a visualização pode parecer estranha, isso porque para arvores filogenéticas, estamos mais acostumados com algo mais ou menos assim:

arvore-eps

Uma arvorezinha com raiz. você pode arrumar a disposição dos vértices usando o comando tkplot assim:

tkplot(graph.formula( 1--2, 2--8, 4--10, 5--9,6--10,7--9,2--3,10--3,9--3 ))

Então a conta a fazer é bem simples se temos que ter n-1 arestas sempre, então contamos quantas arestas ja temos e diminuímos o número de vértices do número de arestas ja existentes -1. A não ser que as arestas estejam zuadas, teremos uma árvore.

Então em python fazemos a função:

def minarv(n, edges):
    return n - len(edges) - 1

E para os dados do exemplo é só usar a função:

n=10
edges=[(1, 2), (2, 8), (4, 10), (5, 9), (6, 10), (7, 9)]
 
print minarv(n,edges)
>>>3

A diferença para o teste é que temos muito mais vértices, e muito mais arestas, então não da para escrever a mão, logo temos que ler o arquivo txt de dados, basicamente a primeira linha é o número de vértices o todas as outras linhas são arestas.

file = open( "/home/augusto/Downloads/rosalind_tree.txt" )
n = file.readline()
edges = []
 
for linha in file:
    edges.append(linha.split())
 
file.close()
 
print int(n)
print len(edges)
 
print minarv(int(n),edges)

Aqui é só tomar cuidado que quando a gente le o arquivo em python, tudo é lido como string, então para o número de vértices temos que usar int para transformar ele em um número inteiro, para poder fazer continhas.
Como as arestas apenas temos que contar por enquanto, então não precisamos se preocupar por enquanto.
E voala, concluímos mais um problema do Rosalind.

Agora entrando nesses problemas de árvores filogenéticas deve ficar bem legal, essa é uma área legal.

Boa semana 🙂

Problema – Rosalind – Open Reading Frames – ORF

DNA_ORF

Esse problema do Rosalind será o com maior número de funções reutilizado desde o primeiro.
Mas vamos la, primeiro temos que entender do que diabos esse problema se trata, eu tive bastante problemas nessa fase, antes de começar a fazer qualquer coisa.

A gente ja viu aqui que podemos, a partir de uma sequência de RNA, prever qual proteína será produzida quando temos em mãos uma tabela que diz qual aminoácido um determinado códon produz. Inclusive vimos algumas propriedade que essa previsão proporciona aqui.

Certo, até aqui tudo bem, mas o RNA vem do DNA, logo a gente pode transcrever DNA em RNA como visto a muito tempo, aqui. E assim usar esse RNA e ver qual proteína ele codifica.

Mas para simplificar as coisas, temos tabelas, como para o RNA, dizendo o que um códon, três letrinhas de DNA geram de aminoácido, sabendo que antes o DNA é transcrito em RNA e depois vira uma proteína, ou algo da sequência de aminoácidos, posso estar sendo muito simplista ou errado nas nomenclaturas aqui.

Essas tabelas tem no wikipedia até, a do DNA para aminoácido aqui e do RNA para aminoácido aqui.

Agora vamos entrar no problema em si.

Acontece que nem todo DNA é transcrito para RNA, o tal do DNA lixo que as vezes a gente ouve falar. Então de uma sequência de DNA, em algum lugar pode começar o pedaço que vai ser transcrito, vai virar RNA, mas a gente pode não saber exatamente onde.

Como assim?

Vamos olhar a sequência de exemplo:

>Rosalind_99 AGCCATGTAGCTAACTCAGGTTACATGGGGATGACCCCGCGACTTGGATTAGAGTCTCTTTTGGAATAAGCCTGAATGATCCGAGTAGCATCTCAG

Certo, a gente tem essa sequência ai, que vai virar RNA, e depois o RNA em códons, que são sequências de letras vão virar aminoácidos que juntos formam proteínas.

Então se a gente separar os tripletos, os grupinhos de três letras aqui, a gente pode prever o aminoácido que vai ser formado.

Mas podemos começar aqui:

AGC - CAT - GTA - GCT - AAC - ...

Virão, começando do início da sequência acima.
Mas derrepente o RNA começa a ser transcrito da segunda base, ou seja assim:

(A) GCC - ATG - TAG - CTA - ACT - CAG - GTT - ...

Notem que os códon, os tripletos começam na segunda base, a primeira, sobrou, ficou ali como “lixo”, ela não influência em nada.
Aqui ficamos com duas possibilidades ja.

E existe uma terceira ainda…

(AG) CCA - TGT - AGC - TAA - CTC - AGG - TTA - ...

Aqui eu deixei em parenteses as duas primeiras bases, AG, elas não entraram no códon. São lixinho e agora formamos três possibilidades.
Será que daria para ter uma quarta?

(AGC) CAT - GTA - GCT - AAC - TCA - ...

Bem se fossemos pensar nessa quarta, agora o primeiro códon é igual ao segundo da nossa primeira opção, alias a sequência toda, so pulou o primeiro códon, por isso da quase na mesma.

Então para cada sequência de DNA, temos 3 possibilidades de Reading Frames, que é o nome chique para esse esquema de dependendo de qual ponto você começa, temos outra sequência de códons, ja que a quarta é igual a primeira pulando o primeiro códon, fica mais fácil pensar nos três primeiros casos so, ja que o quarto caso é o primeiro com menos códon, mas não paramos por aqui ainda.

Certo, só que o DNA não simplesmente faz uma proteína gigante, ele faz muitos tipos de proteínas e outras coisas. Acontece que tem um códon que inicia a proteína. No caso dos Eucariotos é comumente a Metionina (M), que é codificado por “ATG” no DNA e “AUG” no RNA.
Então naqueles Reading Frames ali em cima, das possibilidades, quando aparece um “ATG” a sequência de aminoácido começa a ser feita, e so para na hora que chegar no num códon de parada, so olhar na tabela quais são as paradas.

Mas para complicar mais ainda, além dessas possibilidades, a gente tem duas fitas, o complemento do DNA. Como a gente viu aqui, podemos ver o complemento dessa sequência e ai tem mais três Reading Frames para fazer proteína, é um bocado de possibilidades né.

Agora vamos ver quais são essas possibilidades, como eu disse, a gente vai reciclar um monte de código.

Primeiro, como já falamos ali em cima, como vamos ter que olhar a sequência e seu complemento, temos que definir uma função para ver o reverso do complemente, e isso ja esta pronto aqui:

def revc(sequencia):
    intab = "ACTG"
    outtab = "TGAC"
    trantab = maketrans(intab, outtab)
    return(sequencia.translate(trantab)[::-1])

Lembrando que para usar a função maketrans, tem que abrir o library string (“from string import maketrans”)

A tabela códon -> aminoácido o exercício da, e tem no wikipedia, então sem grandes problemas.

Antes de começar, temos ainda outra atividade que vai sempre ser refeita, que é traduzir os códons, transformar a sequência de códon em aminoácidos, de três em três. Basicamente essa tarefa aqui, so mudar a tabela da de RNA para DNA.

def traduzir_codon(codon):
    amino = ""
    if len(codon) == 3 and DNA_CODON_TABLE.has_key(codon):
        amino = DNA_CODON_TABLE[codon]
    return amino

Agora, aqui tem um problema que eu pedalei muito, ao fazer uma função para passar o DNA para proteína, a gente tem que ter o cuidado de so fazer isso quando tem um códon de stop. Pelo que vi de duvidas no question, foi muita gente além de mim que pedalou aqui.

Como assim?
Bem, pode acontecer de quando a sequência estiver em algum ponto, como por exemplo perto do fim, a gente tem um codon “ATG” para começar, beleza, mas vamos adicionando aminoácidos depois dele, e a sequência acaba e não vimos um códon de stop, ou seja, talvez esteja mais para a frente na sequência e o primer cortou no meio essa sequência, não pegamos um pedaço completinho, sei la, mas se não tem stop, só começa, eu retornei um string vaziu, ta la o ultimo return “” da função, quando chega em um stop, a proteína é retornada. Assim depois eu so vou contabilizar como possibilidade o que não for “”, vazio.

def dna2prot(seqdna,inicio):
    result = ""
    for i in range(inicio, len(seqdna), 3):
        aminoacido = traduzir_codon(seqdna[i:i+3])
        if aminoacido == 'Stop':
            return result
        else:
            result += aminoacido
    return ""

Agora, com tudo definido, é só usar esse monte de função para chegar as possibilidades:

def possibilidades(sequencia):
    possibilidades = []
    inicio = sequencia.find("ATG",0)
    while inicio != -1:
        resul = dna2prot(sequencia,inicio)
        inicio = sequencia.find("ATG",inicio+1)
        if resul != "":
            possibilidades.append(resul)
    return possibilidades

O que estramos fazendo aqui? Começamos uma lista vazia, achamos o primeiro “ATG”, dai a partir dele vemos que proteina é codificada com a função dna2prot, terminando isso, achamos o próximo “ATG” do string, note que usamos a função find aqui, e depois que achamos o primeiro caso, começamos a procurar do inicio+1, para achar o próximo inicio, o próximo “ATG”, so que como eu to lendo letra por letra, eu estou contemplando os três reading frames.
Mas voltando, eu criei uma lista vazia, achei o primeiro “ATG”, comecei um loop, testo se ainda tem “ATG”, se não tiver o find vai dar um resultado de -1, se tiver, ele vai dizer onde começa. Como tem, o loop começa, traduzo a primeira proteína, acho o próximo “ATG” e olho a proteína que saiu, se ela não for vazia, não for “”, ela ponho ela na lista, e começo o loop denovo, e fico nessa até acabarem os “ATG”.

Agora eu tenho todas as possibilidades para essa sequência, mas, eu ainda preciso das possibilidades do complemento, mais como a gente ja definiu uma função para fazer o complemento da sequência, é baba, é so usar ela e depois ver as possibilidades no complemento.

Para não repetir possibilidades ao somar as possibilidades da sequencia e do complemento dela, usamos a função set do python para termos somente termos únicos.

Como fazia tempo que eu não me preocupava em abrir arquivos txt no python, eu sofri um pouco aqui. Mas eu fiz o seguinte:

arquivo = open('/home/augusto/Downloads/rosalind_orf.txt')
 
dados = arquivo.readline()
dataset=""
 
for line in arquivo:
        dataset+=line.strip()

Primeiro com open a gente abre o arquivo, readline() le uma linha do arquivo, como a primeira não me interessa eu coloco ela em algum lugar, as linhas são lidas de forma sequencial com readline(), então tudo que sobrou é a sequência nas linhas quebradas, a primeira linha era a identificação da sequencia que não me interessava, dai então para as linhas no arquivo, eu coloco elas no meu dataset, usando o .strip depois para tirar espaços, que não é o caso, mas para tirar o \n que quebram linhas, que se isso entrar nas funções la em cima vira uma merda.

Dai o dataset fica um string gigante de toda a sequência.

Bem uma observação legal, é que esse esquema de usar find, para achar o próximo caso, num loop igualzinho aquele la, eu vi assistindo as aulinhas de introdução a computação do udacity. La era ensinado como fazer um web-crawler. Um programinha que lia uma página da internet, o html, dai ia olhando o texto até achar um link, guardava o link e ia procurar o próximo, e assim até acabar o página salva. Ora bolas, não é exatamente o que estamos fazendo aqui, so que ao invés de links, paginas de internet são proteínas.
Eu acho legal como o mesmo problema pode de certa forma se repetir em áreas que pra mim parecem tão distintas.

Outra coisa legal, desse exercício é pensar como uma mutação em alguns lugares vai fazer um fuzue danado, pode mudar muito as possibilidades de proteínas que obtemos, enquanto em outros lugares não vai influencias em nada, vai ser completamente neutra, um random walk.

Bem mas é isso, o script vai estar la no repositório do Rosalind além de aqui em baixo.
Foi um problema bem legal, apesar de ter apanhando muito para resolver, principalmente no final onde o problema estava em somente ler os dados de um txt, e não em resolver o problema, é muito frustante conseguir resolver certinho o exemplo e depois sempre falhar com os dados do teste, mas nada que um pouco de perseverança não resolva.

A saída do script abaixo é:

>>> MLLGSFRLIPKETLIQVAGSSPCNLS M MGMTPRLGLESLLE MTPRLGLESLLE

E sem mais enrolar, el script:

from string import maketrans
 
DNA_CODON_TABLE = {
    'TTT': 'F',     'CTT': 'L',     'ATT': 'I',     'GTT': 'V',
    'TTC': 'F',     'CTC': 'L',     'ATC': 'I',     'GTC': 'V',
    'TTA': 'L',     'CTA': 'L',     'ATA': 'I',     'GTA': 'V',
    'TTG': 'L',     'CTG': 'L',     'ATG': 'M',     'GTG': 'V',
    'TCT': 'S',     'CCT': 'P',     'ACT': 'T',     'GCT': 'A',
    'TCC': 'S',     'CCC': 'P',     'ACC': 'T',     'GCC': 'A',
    'TCA': 'S',     'CCA': 'P',     'ACA': 'T',     'GCA': 'A',
    'TCG': 'S',     'CCG': 'P',     'ACG': 'T',     'GCG': 'A',
    'TAT': 'Y',     'CAT': 'H',     'AAT': 'N',     'GAT': 'D',
    'TAC': 'Y',     'CAC': 'H',     'AAC': 'N',     'GAC': 'D',
    'TAA': 'Stop',  'CAA': 'Q',     'AAA': 'K',     'GAA': 'E',
    'TAG': 'Stop',  'CAG': 'Q',     'AAG': 'K',     'GAG': 'E',
    'TGT': 'C',     'CGT': 'R',     'AGT': 'S',     'GGT': 'G',
    'TGC': 'C',     'CGC': 'R',     'AGC': 'S',     'GGC': 'G',
    'TGA': 'Stop',  'CGA': 'R',     'AGA': 'R',     'GGA': 'G',
    'TGG': 'W',     'CGG': 'R',     'AGG': 'R',     'GGG': 'G'
}
 
def revc(sequencia):
    intab = "ACTG"
    outtab = "TGAC"
    trantab = maketrans(intab, outtab)
    return(sequencia.translate(trantab)[::-1])
 
def traduzir_codon(codon):
    amino = ""
    if len(codon) == 3 and DNA_CODON_TABLE.has_key(codon):
        amino = DNA_CODON_TABLE[codon]
    return amino
 
def dna2prot(seqdna,inicio):
    result = ""
    for i in range(inicio, len(seqdna), 3):
        aminoacido = traduzir_codon(seqdna[i:i+3])
        if aminoacido == 'Stop':
            return result
        else:
            result += aminoacido
    return ""
 
def possibilidades(sequencia):
    possibilidades = []
    inicio = sequencia.find("ATG",0)
    while inicio != -1:
        resul = dna2prot(sequencia,inicio)
        inicio = sequencia.find("ATG",inicio+1)
        if resul != "":
            possibilidades.append(resul)
    return possibilidades
 
exemplo = "AGCCATGTAGCTAACTCAGGTTACATGGGGATGACCCCGCGACTTGGATTAGAGTCTCTTTTGGAATAAGCCTGAATGATCCGAGTAGCATCTCAG"
#arquivo = open('/home/augusto/Downloads/rosalind_orf.txt')
 
#dados = arquivo.readline()
#dataset=""
 
#for line in arquivo:
#        dataset+=line.strip()
 
#print str(dataset)
 
 
A = possibilidades(exemplo)
B = possibilidades(revc(exemplo))
 
print "\n".join(set(A + B))

Problema – Rosalind – Calculating Protein Mass – PRTM

E aqui vamos nos com o problemas do projeto Rosalind. Bem, em um problema anterior, esse aqui mais especificamente, a gente viu como podemos a partir de uma sequencia, ver qual proteína será produzida. Agora uma coisa legal é que se a gente sabe o que vai ser produzido e algumas características dos tijolinhos que produzem, podemos prever varias características do que ta no final da linha de produção desse sistema todo ae.

O problema te fornece uma tabela com as massas para cada aminoácido, essa tabela tem la na wikipedia, assim como alguns pacotes que englobam funções voltadas para bioinformática tem essas e outras informações para previsões úteis.

Então o python para resolver o serviço, a idéia é bem simples, criamos um dicionario no python com o nome do aminoácido e sua respectiva massa. Agora eu entendi porque diabos as pessoas usam um dicionario e não uma lista para fazer isso. Se essa lista fosse muito longa, sempre para procurar um peso, na pior das hipóteses teríamos que rodar toda a lista procurando, em dados como um dicionario (usando {} chaves), o dicionário tem uma função ja feita nele de hash-table, isso basicamente é um índice remissivo, que torna a busca mais rápida, lembre-se que se queremos achar uma palavra, se você olhar o livro todo você ta ferrado, mas se você pega a primeira letra e vai procurar no índice remissivo no final do livro, serão menos páginas para olhar, logo a busca é mais rápida porque temos que olhar menos coisas, particularmente aqui isso não deve ter efeito nenhum, mas se a lista fosse grande teria 🙂

Certo, mas sabendo disso, nos optamos pelo dicionario ao invés de uma lista e dai é so somar letra a letra os pesos e temos as resposta.
Como so temos uma sequência, a leitura do arquivo é simples, eu so não sei se desse jeito o arquivo fica na memória, deve ficar la tomando espaço caso não usemos o close()

Mas então a função final fica assim:

mono_mass_table = {
    'A': 71.03711,
    'C': 103.00919,
    'D': 115.02694,
    'E': 129.04259,
    'F': 147.06841,
    'G': 57.02146,
    'H': 137.05891,
    'I': 113.08406,
    'K': 128.09496,
    'L': 113.08406,
    'M': 131.04049,
    'N': 114.04293,
    'P': 97.05276,
    'Q': 128.05858,
    'R': 156.10111,
    'S': 87.03203,
    'T': 101.04768,
    'V': 99.06841,
    'W': 186.07931,
    'Y': 163.06333,
}
 
def peso(seq_prot):
    n = 0
    for i in seq_prot:
        n += mono_mass_table[i]
    return n
 
dados = "SKADYEK"
#dados = open("/home/augusto/Downloads/rosalind_prtm.txt").read().strip()
 
print peso(dados)

Que vai nos dar o resultado de:

>>> 821.39192

E todos vivemos felizes para sempre.