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)

Leave a Reply

Your email address will not be published. Required fields are marked *