procura de sequencias DNA

[Guilherme Loss]

Olá,

estou com um pouco de dificuldade em fazer um pequeno contador (acumulador)...

tenho 2 arquivos:

ambos no formato FASTA

ou seja

>seq1

AAAAAAAAAAAA

>seq2

TTTTTTTTTTTTT

gostaria de saber como contar o numero de ocorrências (match perfeito) das sequencias de um arquivo no segundo...

o segundo arquivo geralmente é bastante grande (na casa dos Mb)...

tenho uma implementação desse contador utilizando dicionários...

mas é muuuuuuuuuuito lenta, e tem o problema inerente de dicionários não fornecerem a mesma ordem do input..

o que eu gostaria na verdade seria algo tipo um 'grep -c' que funcionasse em loop para cada sequencia...

andei olhando e tem o módulo subprocess, o qual acredito que teria como chamar o grep...

Mas não sei ao certo se é possivel utilizar isso....

Alguém conhece alguma outra forma de fazer esse contador de maneira mais eficiente??

vlw

abraço

--
Guilherme Loss de Morais M.Sc.
Ph.D. student in Cellular and Molecular Biology
CBiot - UFRGS - Brazil
http://www.ufrgs.br/rnai/LGPP.htm

[Laerte M. Rodrigues]

ola guilherme,

tarbalho com analise de sequencias, e seu problema eh simples, mas o volume de dados é grande, para isto, lhe aconselho utilizar threads ou a biblioteca mpi4py para implementação,

mas vc pode fazer alinhamento de genes, existem dois algoritmos, com alinhamento local ou globa, no caso, o custo de espaço computacional seria 

O(len(a)*len(b)) [tamanho da matriz de a * a matriz de b]

--
------------------------------------
Grupo Python-Brasil
http://www.python.org.br/wiki/AntesDePerguntar
 
<*> Para visitar o site do grupo na web, acesse:
http://groups.google.com/group/python-brasil
 
<*> Para sair deste grupo, envie um e-mail para:
python-brasi...@googlegroups.com

--
Grato,

Laerte Mateus Rodrigues
Mestrando em informática (PUC Minas)

http://laertemateus.com.br

[Christian S. Perone]

Se você não estiver fazendo isto por aprendizado, deveria estar usando frameworks já bem estabelecidos (aka. biopython), não tem porque reinventar a roda sem motivos.

[1] http://biopython.org/DIST/docs/tutorial/Tutorial.html

2012/5/23 Guilherme Loss <guilher...@gmail.com>

--
------------------------------------
Grupo Python-Brasil
http://www.python.org.br/wiki/AntesDePerguntar
 
<*> Para visitar o site do grupo na web, acesse:
http://groups.google.com/group/python-brasil
 
<*> Para sair deste grupo, envie um e-mail para:
python-brasi...@googlegroups.com

--
"Forgive, O Lord, my little jokes on Thee, and I'll forgive Thy great big joke on me."
http://pyevolve.sourceforge.net/wordpress/

[Guilherme Loss]

Oi Laerte, 

é, esse problema é simples, mas escala (O) bem mal mesmo.. ehehe

quanto a utilizar algoritmos de busca local e global.. ja tentei utilizar implementação de ambos como o fasta e o BLAST...

mas o propósito destes é outro ... no caso para esses eles permitem gaps, e mismatchs entre as sequencias 'iscas' e o banco de dados...

e como resultado ele não me dá o número de ocorrencias, e permite variaçãoes entre a isca e o banco de dados

quanto a usar mpi, nunca tentei.. vou dar uma estudada...

vlw abraço

 

2012/5/23 Laerte M. Rodrigues <laerte...@gmail.com>

[Laerte M. Rodrigues]

tenta entao distancia de levenshtein (acho q eh isso msm) onde é encontrado o maior subsequencia comum entre duas strings, é a base dos algoritmos de alinhamento, mas ele nao possui penalidade, apenas verifica o mair matching possivel

[Guilherme Loss]

não é para apreendizado...

utilizo biopython já a algum tempo...

vlw, já conheço o tutorial...

mas mesmo utiliando biopython... caio no mesmo problema de escala...

 vlw

2012/5/23 Christian S. Perone <christia...@gmail.com>

[Guilherme Loss]

ok

me recordo desse algorítimo.. da distancia de levenstein, vou dar uma olhada...

vlw

[Eder Souza]

Infelizmente Leveintein nao vai ajudar a ganhar desempenho ele apenas vai lhe dar as sequencias mais parecidas ou exatas, neste caso você vai precisar das sequencias exatas portanto descarte este algoritmo ele possui um intensivo processamento, o seu caso esta mais para brute-force, existe uma saída simples para isso !

1º Eu acho que os seus arquivos presas devem ser grandes e os iscas nem tanto certo ?

2º pq nao recriar todas as suas presas para o formato .Mat com o (scipy.io.savemat) você teria este trabalho inicial, mas ganharia mais tarde em desempenho, assim abriria todo o arquivo já como dicionario diretamente no python, cada sequencia seria sua chave contendo os valores de cada sequencia; ao comparar com o arquivo isca chame diretamente a chave necessária por exemplo: Presa[seq3] isso vai cai diretamente no conteúdo da seq3 sem passar por todo o resto, pegue os valores contidos lá e faça o match, vai ganhar mt em desempenho !

Documenação:

http://docs.scipy.org/doc/scipy/reference/generated/scipy.io.savemat.html

Eng Eder de Souza

2012/5/23 Guilherme Loss <guilher...@gmail.com>

[Guilherme Loss]

Vlw Eder,

é bem brute-force mesmo..

sim, as iscas são pequenas e o banco de dados extenso...

vou estudar sobre esse modulo...

outra alternativa que pensei.. foi pré-processar o meu banco de dados, utilizando collections 

onde após eu ficaria sequencia '\t' frequencia....ai fica tranquilo pegar com dicts...

ainda é brutal-force, mas melhora o desempenho (menos tempo pra procura)

abraço

2012/5/23 Eder Souza <ederw...@gmail.com>

[Pedro Werneck]

Que tal usar SQLite?


2012/5/23 Guilherme Loss <guilher...@gmail.com>:

---
Pedro Werneck

[Frederico Arnoldi]

Oi Guilherme

> gostaria de saber como contar o numero de ocorrências (match perfeito) das
> sequencias de um arquivo no segundo...

Se você precisa de um "match perfeito", você pode carregar as
sequencias como strings e usar o "find" que só acha "perfect matches"
não sobrepostos. Você só precisaria fazer um alinhamento e usar
distâncias caso fosse necessário buscar sequencias semelhantes.
Iguais, o "find" faz o serviço e bem rápido.

Se entendi seu problema minha sugestão tabajara é:

1 - Abra o seu arquivo 1 e carregue suas sequencias individualmente.
2 - Abra seu arquivo 2 com um open().read() (se esse for na ordem de
Mb você não vai ter problema de memória... se for algumas centenas de
Mb talvez seja interessante carrega-lo em blocos... mas cuidado para
abrir em blocos que não quebre uma sequencia no meio)
3 - Para cada sequencia do seu arquivo 1, use o find até o final do
arquivo 2 ou dos blocos. Isto é enquanto o valor retornado por find
for menor que a string do arquivo 2 você busca os matches.
4 - Salve a localização das sequencias apenas (o valor retornado pelo
find, e esse mesmo valor mais o comprimento do match e um id da seq do
arquivo 1) e não a sequencia propriamente dita. Use SQLite, como já
sugerido, ou em um arquivo, ou simplesmente deixe na memoria...
afinal, vc estara salvando, por exemplo, uma lista de tuplas com três
inteiros (o inicio e o fim do match no arquivo 2 e um id da sequência
do arquivo 1) na memoria.
5- Mapeie seu arquivo 2, de forma q vc saiba onde está o começo e o
fim de cada sequencia dele.
6- Tendo o mapa do arquivo 2 e a localização do matches.... você tem o
que quer

A não ser que seu arquivo seja realmente grande, acho que fara cada
seq em questão de segundos com o find.

Abs,
Fred

[Guilherme Loss]

obrigado pela ajuda Fred

abraço

2012/5/23 Frederico Arnoldi <fre...@hotmail.com>

--
------------------------------------
Grupo Python-Brasil
http://www.python.org.br/wiki/AntesDePerguntar

<*> Para visitar o site do grupo na web, acesse:
   http://groups.google.com/group/python-brasil

<*> Para sair deste grupo, envie um e-mail para:
   python-brasi...@googlegroups.com

[Eder Souza]

Bom fiz os testes agora escrevi o código em 10 minutos, Guilherme eu fiz os testes em escala pequena não sei como isso se comportaria ao usar os seus arquivos !

Aqui utilizando o método convencional abrindo os dois arquivos e fazendo o match!

python DNA.py

HSW7  encontrada 3 vezes

HSW5  encontrada 2 vezes

HSW9  encontrada 1 vezes

Tempo Gasto: 0.045 s

Aqui abrindo o arquivo gravado em .mat e fazendo o match por ele!

python.exe DNAMatchDict.py

HSW7  encontrada 3 vezes

HSW5  encontrada 2 vezes

HSW9  encontrada 1 vezes

Tempo Gasto: 0.008 s

Gustavo teste, aqui o source que converte seu arquivo presa e o coloca diretamente em um arquivo estruturado como dicionário: http://pastebin.com/H6rQXNYK

E Aqui como fazer o match abrindo o arquivo .mat: http://pastebin.com/KEjbrARC

Importante: use o scipy-0.9.0 ou superior me lembro que com versões menores tive problemas com String !

Eng Eder de Souza

2012/5/23 Frederico Arnoldi <fre...@hotmail.com>