Quanto č utile/interessante questa discussione:
Autore |
Discussione |
|
HnACP
Nuovo Arrivato
Prov.: MI
Cittā: Milano
78 Messaggi |
Inserito il - 13 maggio 2009 : 09:46:07
|
ciao a tutti! ho scritto un codice per trovare siti di restrizione in una sequenza fasta, ma nn riesco a dare l'output, nel senso che quando lo faccio girare nn mi da niente. Copio il codice qua: qualcuno mi sa aiutare? Se funziona magari puo' essere di aiuto a qualcuno!
from Bio import SeqIO
from string import *
handle = open("unr.fasta")
for sequnr in SeqIO.parse(handle, "fasta") :
def digest(sequnr, enz) :
EcoRI = 'GAATTC'
BamHI = 'GGATCC'
HindIII = 'AAGCTT'
enz = [ EcoRI, BamHI, HindIII ]
res = []
site = sequnr.find(enz)
while site != -1 :
res.append(site)
site = sequnr.find(enz, site + 1)
return res
handle.close()
|
|
|
dallolio_gm
Moderatore
Prov.: Bo!
Cittā: Barcelona/Bologna
2445 Messaggi |
Inserito il - 13 maggio 2009 : 10:04:51
|
ciao, forse l'indentazione é sbagliata, o si é persa qualche riga nel post, perķ non capisco bene se stai definendo la funzione digest per ogni esecuzione del ciclo for. Se il codice é quello, non ti funziona perché non stai chiamando la funzione digest(), la stai solo definendo.
Innanzitutto, una buona regola quando si scrivono programmi é quella di pensare subito a un test case per verificare che i risultati siano corretti. Per esempio, dovresti scrivere un paio di sequenze casuali e verificare che il tuo script funzioni correttamente.
Secondo, su biopython esiste un modulo chiamato Restriction, per lavorare su enzimi di restrizione. Perķ il modulo sembra essere piuttosto vecchio e almeno a me, non funziona. Ti sconsiglio di utilizzarlo, e sarebbe il caso di avvertire gli sviluppatori (posso farlo io se vuoi).
Terzo, per lavorare su motivi e matrici pssm, ti consiglio questo modulo: - http://fraenkel.mit.edu/TAMO/TAMO_tutorial.html E' ottimo e se dai una occhiata al tutorial, puoi vedere che esiste un metodo .scan() che permette di cercare le occorrenze del motivo su una sequenza. In ogni caso, ti sconsiglio di utilizzare il semplice .find, piuttosto dai una occhiata alle espressioni regolari (modulo re).
ok, adesso aspetto una tua risposta :) |
Il mio blog di bioinformatics (inglese): BioinfoBlog Sono un po' lento a rispondere, posso tardare anche qualche giorno... ma abbiate fede! :-) |
|
|
dallolio_gm
Moderatore
Prov.: Bo!
Cittā: Barcelona/Bologna
2445 Messaggi |
Inserito il - 13 maggio 2009 : 10:24:13
|
from Bio import SeqIO
from TAMO import MotifTools
# from string import * # non necessario
handle = open("unr.fasta", 'r')
EcoRI = MotifTools.Motif_from_text('GAATTC')
BamHI = MotifTools.Motif_from_text('GGATCC')
HindIII = MotifTools.Motif_from_text('AAGCTT')
enzimes = [ EcoRI, BamHI, HindIII ]
for sequnr in SeqIO.parse(handle, "fasta"):
for enz in enzimes:
digest(sequnr, enz)
def digest(sequnr, enz) :
"""
trova tutte le occorrenze di un motivo in una sequenza. (corretto?)
esempio:
>>> from StringIO import StringIO
>>> from TAMO import MotifTools
>>> from Bio import SeqIO
>>> fastafile = StringIO('>seq1\nGGATcCCGGATGCATGCTA')
>>> fastafile.seek(0)
>>> enz = MotifTools.Motif_from_text('GGATCC')
>>> print digest(SeqIO.parse(fastafile, 'fasta'), enz)
(['GGATcC'], [(0, 5)], [11.685195480635185])
"""
occurrencies = enz.scan(sequnr.seq.tostring())
return occurencies
Probabilmente puoi spostare alcune chiamate all'interno della funzione, come la creazione del motivo e il parsing del file. Scrivi prima il doctest e poi ti sará facile decidere come impostarlo.
ultima cosa: non usare il metodo find di biopython, non é case insensitive e non distingue le N come basi speciali. |
Il mio blog di bioinformatics (inglese): BioinfoBlog Sono un po' lento a rispondere, posso tardare anche qualche giorno... ma abbiate fede! :-) |
|
|
HnACP
Nuovo Arrivato
Prov.: MI
Cittā: Milano
78 Messaggi |
Inserito il - 13 maggio 2009 : 23:12:26
|
ciao! quel codice, un po arrangiato da me, l'avevo trovato sul tutorial di biopython dell'istituto pasteaur, solo che č cosi stringato che nn si capisce niente...
in pratica volevo, da una sequenza in formato fasta, vedere se dei dati siti di restrizione erano presenti in quella sequenza e quanti erano. Con un altro script sono riuscito solo che mi indicava solo il primo sito e poi si fermava! o cmq se nn sono siti di restrizione magari cercare una data sequenza di un promotore, ma penso che la procedura sia la stessa.
quel modulo TAMO sembra molto interesseante, solo che su linux nn riesco ad installarlo...
sul tutorial usa quasi sempre find... che poi ha funzionato solo una volta con me..
alla fine io mi basavo solo sul tutorial ufficiale di biopython e quello francese.. ma vedo che hanno cose datate... potrebbero aggiornarsi un po, no?
visto che nn riesco ad installare il modulo TAMO nn c'č il modo di fare lo script in puro python? |
|
|
dallolio_gm
Moderatore
Prov.: Bo!
Cittā: Barcelona/Bologna
2445 Messaggi |
Inserito il - 14 maggio 2009 : 00:36:02
|
Che problemi hai ad installare TAMO? Basta scompattare il file e dare 'sudo python setup.py install'. E' richiesto il modulo numpy... lo puoi installare tramite easy_install, o come pacchetto nella tua distro.
Se non vuoi utilizzare moduli esterni (male... mai riscrivere quello che giā esiste), usa le espressioni regolari. Qualcosa del genere:
>>> import re
>>> seq1 = 'ACGATGGCTACGTACGTACGTACGA'
>>> seq2 = 'AAAAAAAACCAAAAAAAAAAAAAAA'
>>> re.findall('AC', seq1)
['AC']
>>> for match in re.finditer('AA', seq2):
print match.start(), match.end(), match.group()
# nota come anche con le regex il comportamento é incorretto, perché non si tiene conto che una N rappresenta qualsiasi nucleotide:
>>> seq3 = 'ANAA'
>>> re.findall('AA', seq3)
|
Il mio blog di bioinformatics (inglese): BioinfoBlog Sono un po' lento a rispondere, posso tardare anche qualche giorno... ma abbiate fede! :-) |
|
|
HnACP
Nuovo Arrivato
Prov.: MI
Cittā: Milano
78 Messaggi |
Inserito il - 14 maggio 2009 : 23:05:51
|
ciao! si ho provato ad installarlo ma ad un certo punto mi da errore:
creating build/temp.linux-i686-2.6/TAMO/MD/MDsupport_source gcc -pthread -fno-strict-aliasing -DNDEBUG -g -fwrapv -O2 -Wall -Wstrict-prototypes -fPIC -I/usr/include/python2.6 -c TAMO/MD/MDsupport_source/MDsupport.cxx -o build/temp.linux-i686-2.6/TAMO/MD/MDsupport_source/MDsupport.o gcc: error trying to exec 'cc1plus': execvp: No such file or directory error: command 'gcc' failed with exit status 1
|
|
|
|
Discussione |
|
|
|
Quanto č utile/interessante questa discussione:
MolecularLab.it |
© 2003-18 MolecularLab.it |
|
|
|