| Autore |
Discussione |
|
|
HnACP
Nuovo Arrivato
Prov.: MI
Città: Milano
78 Messaggi |
 Inserito il - 18 maggio 2009 : 15:52:12
|
ciao!
devo mappare delle seq su un genoma. Ho un file con tante sequenze in questo modo:
1||Count=1 TTCAATAGTTTAAGAAATTC R2 0 0 2
2||Count=1 CGAAAAAAAAAGTGTACGCG U2 0 0 1 3L 10367503 F
3||Count=1 CTGCACTTTCCTCTGGTTCG NM 0 0 0
4||Count=1 GTCTGAGTATTTGATGGTTT NM 0 0 0
5||Count=1 ATCTGATTTCAAACGAAGCG NM 0 0 0
6||Count=8 ATCTTAGT
Devo mappare solo le sequenze con la sigla NM, e ad ogni ciclo devo togliere una base di tali sequenze e rimappare sul genoma fino a fermarmi a 15 basi di lunghezza. Se count e' piu di 1 allora devo farlo per quel numero di volte.
Il genoma e' in formato .raw
qualcuno mi puo aiutare
grazie!
|
|
|
|
|
HnACP
Nuovo Arrivato
Prov.: MI
Città: Milano
78 Messaggi |
Inserito il - 19 maggio 2009 : 14:36:09
|
modifico un paio di cose:
-il file.txt ha le sequenze formattate in questo modo:
2||Count=1 CTATTCTTTTTTTTTAAATA U1 0 1 2R
3||Count=1 AACATTAATTGTTTTAATAA U1 0 1 3L
4||Count=1 CCCTCCGCCTTTTTTGCCAA NM 0 0
5||Count=1 CCCCCAGCGGCTGCACCAGA NM 0 0
6||Count=1 AAGGAAAAAAATTTTATTAT NM 0 0
7||Count=1 CCCGCCCTTCCCTTTATAAA NM 0 0
-il file del genoma e' .raw
-devo in pratica selezionare solo quelle sequenze che sono NM, ma nn riesco ad "estrarle dal file per poi usarle".
-una volta che le ho "selezionate" devo togliere una base alla fine e mapparle sul genoma che sta sul file genoma.raw.
-quelle sequenze che non mappano le butto via
-per quelle che rimangono tolgo un'altra base e rimappo sul genoma
-continuo con questa procedura fino a che le sequenze arrivano a 14 basi di lunghezza
-salvo i risultati (cioe' le sequenze che ho trovato che mappano sul genoma) in un file formato FASTA
so che devo fare dei cicli FOR pero' nn sono capace di combinare NR alla sequenza, togliere la base e infine aprire il file genoma.raw per mappare e continuare cosi...
spero qualcuno mi aiuti!
grazie mille!
ciao! |
 |
|
|
dallolio_gm
Moderatore
Prov.: Bo!
Città: Barcelona/Bologna
2445 Messaggi |
Inserito il - 20 maggio 2009 : 23:22:19
|
ciao,
l'ultima versione di biopython supporta sequenze in formato raw, ma non so a che punto sia lo sviluppo, ti conviene chiedere nella ml di biopython.
Stai cercando di assemblare i risultati di un sequenziamento? ci sono già delle pipeline per farlo, non ti so dire bene quali.
Citazione:
-devo in pratica selezionare solo quelle sequenze che sono NM, ma nn riesco ad "estrarle dal file per poi usarle".
in che senso? non puoi usare grep, o un equivalente in python?
Citazione:
-una volta che le ho "selezionate" devo togliere una base alla fine e mapparle sul genoma che sta sul file genoma.raw.
per lanciare allineamenti da linea di comando, io ho usato questo tool: www.ebi.ac.uk/~guy/exonerate
Citazione:
-quelle sequenze che non mappano le butto via
-per quelle che rimangono tolgo un'altra base e rimappo sul genoma
non vorrei ficcare il naso, ma non capisco bene il significato di questo procedimento... perché togliere una base alla fine e rimappare?
Citazione:
-salvo i risultati (cioe' le sequenze che ho trovato che mappano sul genoma) in un file formato FASTA
questo lo puoi fare con biopython, o anche manualmente visto che il formato fasta non é troppo complicato. |
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 - 21 maggio 2009 : 10:04:15
|
ciao grazie per la risposta!
le sequenze vengono da un sequenziatore, e devo verificare che quelle con NM, cioe' sequenze che nn risultano mappare sul genoma, togliendo man mano basi trovo delle corrispondenze. Lo so che nn ha molto senso ma e' solo una verifica di dati che si hanno gia'.
per leggere le righe forse ho trovato una soluzione.
ma io la mappatura la devo fare tramite codice che scrivo io, senza usare programmi esterni.
nn posso farlo manualmente visto che sono milioni di sequenze! |
|
|
|
dallolio_gm
Moderatore
Prov.: Bo!
Città: Barcelona/Bologna
2445 Messaggi |
Inserito il - 21 maggio 2009 : 10:28:42
|
Citazione:
le sequenze vengono da un sequenziatore, e devo verificare che quelle con NM, cioe' sequenze che nn risultano mappare sul genoma, togliendo man mano basi trovo delle corrispondenze. Lo so che nn ha molto senso ma e' solo una verifica di dati che si hanno gia'.
ok, ma non sarebbe piú comodo utilizzare un tool di allineamento, introducendo parametri via via piú permissivi?
Puo' anche darsi che una sequenza non si mappi correttamente perché contiene una base diversa (o in piu', in meno) rispetto al genoma con cui la stai confrontando, ricordati che ci possono essere errori o polimorfismi anche sullo stesso genoma standard.
|
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 - 21 maggio 2009 : 11:22:17
|
| si... ma lo scopo è che scrivi io lo script! |
|
|
|
dallolio_gm
Moderatore
Prov.: Bo!
Città: Barcelona/Bologna
2445 Messaggi |
Inserito il - 21 maggio 2009 : 13:49:03
|
per doverlo fare bene, dovresti riscrivere un algoritmo per l'allineamento. Sei sicuro che ne valga la pena?
Anche saper interfacciarsi a tool da riga di comando é una dote molto utile per un bioinformatico. Non dimenticare che ci sarà anche la fase di parsing dei risultati, di testing (che in genere porta via il doppio del tempo rispetto al resto), quella di interagire con diversi formati, etc..
Inoltre, ricordati che noi bioinformatici siamo pagati per ottenere risultati, non per scrivere programmi :) |
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 - 21 maggio 2009 : 16:27:50
|
| nn so che dire... me l'hanno detto loro di fare questo script... mi hanno detto che ci vogliono 3 min per scriverlo tramite un paio di cicli for.. |
|
|
|
Ra1D
Utente Junior
Prov.: Bergamo
Città: Treviglio
174 Messaggi |
Inserito il - 22 maggio 2009 : 01:31:37
|
Citazione:
Messaggio inserito da HnACP
nn so che dire... me l'hanno detto loro di fare questo script... mi hanno detto che ci vogliono 3 min per scriverlo tramite un paio di cicli for..
Magari non hanno ben chiaro il concetto di ottimizzazione degli script...magari ci vogliono veramente 3 minuti a scriverlo con un paio di cicli for, poi però ci mette 3 anni a scorrere su un milione di basi 
Cmq potresti usare le regex per controllare riga per riga la presenza di NM e dirgli di prendere come input la sequenza precendente (ora saprei farlo in perl, in python nn ricordo bene la sintassi)
Il secondo ciclo lo fai ripetere mappando sul genoma e togliendo una lettera finchè non arrivi a lunghezza 14....se usi un array puoi usare l'opzione [:-1] e fare il check che l'array sia più lungo di 14.
Non conosco però il formato .raw del genoma...e come output di mappatura cosa ti servirebbe? la posizione sul genoma in cui la tua sequenza combacia? |
...e vicinissimo al mio sapere si accampa la tenebra della mia ignoranza... |
|
|
|
HnACP
Nuovo Arrivato
Prov.: MI
Città: Milano
78 Messaggi |
Inserito il - 22 maggio 2009 : 10:06:00
|
ciao!
ho guardato le regex e sembra abbastanza comodo..
avevo calcolato di fare 2 array
come output ho bisogno di un file fasta cosi:
>1 sequenza |
|
|
|
dallolio_gm
Moderatore
Prov.: Bo!
Città: Barcelona/Bologna
2445 Messaggi |
Inserito il - 22 maggio 2009 : 10:35:50
|
il tuo problema é che ci potrebbero essere delle divergenze tra quello che é sul genoma e le sequenze che tu vuoi annotare.
Per esempio, immagina di avere un genoma così:
>genoma
CCCCCAAANAAACCCC
e di voler mappare una sequenza così:
>seq1
AAAAAAA
utilizzando un semplice confronto fra stringhe, non riuscirai mai a trovare il match giusto, perchè un N non viene riconosciuta come una A (e ci potrebbero essere numerose varianti).
Anche con le regex, non risolveresti il problema, a meno di non scrivere espressioni regolari abbastanza complesse (qualcosa del tipo: [AN]{7,7})
Un'altra soluzione potrebbe essere quella di usare il modulo TAMO che ti dicevo nell'altro topic, che nono soffre di queste limitazioni.
Però l'ideale sarebbe utilizzare un programma di allineamento, capace di tenere in conto sia gap che sostituzioni.
Se vuoi fare le cose per bene, scriviti dei test, ovvero dei casi in cui l'approccio più semplice fallirebbe (come quello che ti ho postato adesso), e questo ti aiuterà a scrivere un programma corretto. |
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 - 22 maggio 2009 : 10:49:46
|
ciao!
si quello l'ho fatto presente, e mi hanno detto di scartare le sequenze che hanno basi sconosciute N.
eh lo so che TAMO è utile, ma come dicevo nell'altro post mi viene un errore.. cmq adesso riproverò ad installarlo, perchè ieri cambiando layout della tastiera su ubuntu al successivo riavvio ha dato una serie di errori e nn si piu riavviato e quindi l'ho dovuto reinstallare. Quindi installo biopython, numpy e poi TAMO se nn sbaglio.
grazie mille dei consigli! |
|
|
|
HnACP
Nuovo Arrivato
Prov.: MI
Città: Milano
78 Messaggi |
Inserito il - 22 maggio 2009 : 15:37:42
|
| alla fine ce l'ho fatta, se volete copio il codice qui in modo da discuterne o in modo che sia di aiuto "didatticamente"... |
|
|
|
dallolio_gm
Moderatore
Prov.: Bo!
Città: Barcelona/Bologna
2445 Messaggi |
Inserito il - 22 maggio 2009 : 15:43:21
|
Grande!
Se vuoi posta pure il codice (puoi usare il tag [ code ] per impaginarlo meglio) |
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 - 22 maggio 2009 : 21:39:57
|
#-*- coding: utf-8 -*-
import string, re, sys
class Parser(object):
filetxtname = "female_unr_seq.txt"
filerawname = "genome.raw"
filefastaname = "find_seq.fasta"
filetxtptr = None
filerawstr = ""
arraystrmatched = []
def __init__(self):
name = raw_input("insert txt filename (if blank will be used 'female_unr_seq.txt'): ")
if name != "":
self.filetxtname = name
name = raw_input("insert raw filename (if blank will be used 'genome.raw'): ")
if name != "":
self.filerawname = name
try:
self.filetxtptr = open(self.filetxtname,'r')
print "opened "+self.filetxtname+"!"
except:
print "unable to open txt file: "+self.filetxtname
sys.exit(0)
try:
filerawptr = open(self.filerawname,'r')
print "opened "+self.filerawname+"!"
try:
self.filerawstr = filerawptr.read()
self.filerawstr = self.filerawstr.replace ( "\n", "" )
filerawptr.close()
except:
print "unable to read raw file: "+self.filerawname
filerawptr.close()
sys.exit(0)
except:
print "unable to open raw file: "+self.filerawname
sys.exit(0)
return
def parse(self):
try:
for line in self.filetxtptr:
linetoken = string.split(line)
if linetoken[2] == "NM":
print "compatible NM sequence: "+line[:-1]
#controlliamo se la sequenza contiene "N"
if string.find(linetoken[1], "N") == -1:
print "compatible sequence without N factor found"
sequence = linetoken[1][:-1]
while len(sequence) > 15:
matches = re.search(sequence,self.filerawstr,re.IGNORECASE)
if matches is None:
sequence = sequence[:-1]
print "scaling to "+sequence
else:
print "sequence found!"
self.arraystrmatched.append(sequence)
break;
if len(self.arraystrmatched) != 0:
print "\nbuilding .FASTA entity...",
try:
fastaptr = open(self.filefastaname,"w")
for index in range(len(self.arraystrmatched)):
fastaptr.write(">"+str(index+1)+" "+self.arraystrmatched[index]+"\n")
fastaptr.close()
print "done!"
except:
print "error!"
else:
print "no matches found!"
finally:
self.filetxtptr.close()
def main():
print "beginning sequential txt parser"
parser = Parser()
parser.parse()
return
#entry-point
main() |
|
|
|
Ra1D
Utente Junior
Prov.: Bergamo
Città: Treviglio
174 Messaggi |
Inserito il - 22 maggio 2009 : 21:44:46
|
Posta posta.....o se non puoi scrivere tutto il codice per qualche cavillo burocratico, mettici solo le parti fondamentali 
Cmq penso che la questione fosse di fargli scrivere uno script a livello didattico, per il pattern matching, senza affrontare il discorso di omologia o varianti biologiche, nè il discorso di ottimizzazione delle procedure di ricerca delle stringhe all'interno del testo...anche xkè sennò col classico confronto lettera per lettera ci vogliono anni, e bisogna introdurre i concetti di programmazione dinamica, o le euristiche di ricerca come allineamenti globali o locali.
Cmq sarebbe interessante vedere qual'è il livello prestazionale, magari potresti valutare il tempo di mappaggio rispetto alle dimensioni del genoma ed al numero di stringhe che hai valutato 
|
...e vicinissimo al mio sapere si accampa la tenebra della mia ignoranza... |
|
|
|
HnACP
Nuovo Arrivato
Prov.: MI
Città: Milano
78 Messaggi |
Inserito il - 22 maggio 2009 : 22:53:27
|
ho gia postato
scrivimi il codice per vedere il tempo di mappaggio, sembra interessante! |
|
|
|
HnACP
Nuovo Arrivato
Prov.: MI
Città: Milano
78 Messaggi |
Inserito il - 25 maggio 2009 : 19:17:11
|
| l'ho testato ma ci mette troppo tempo! necessita di una ottimizzazione! |
|
|
|
Ra1D
Utente Junior
Prov.: Bergamo
Città: Treviglio
174 Messaggi |
Inserito il - 26 maggio 2009 : 00:32:40
|
Per esempio, a livello di ottimizzazione concettuale, a prescindere dal codice ti faccio presente questo punto:
while len(sequence) > 15:
matches =re.search(sequence,self.filerawstr,re.IGNORECASE)
if matches is None:
sequence = sequence[:-1]
Ignorando il re.search che magari non è la funzione di ricerca più adeguata, tu hai impostato il ciclo chiedendo di partire dalla lunghezza della sequenza, e poi scalando via via.
Se per esempio la tua sequenza è lunga 20 ma c'è un match solo a 15, tu farai scorrere l'intero genoma 5 volte inutilmente, e ogni volta ripartirai da capo, ignorando che per 5 volte verrà trovato un match di 15...ma non verrà evidenziato xkè la tua sequenza durante il ciclo è lunga rispettivamente 20, 19, 18, 17, 16 nucleotidi.
Un'idea se non ricordo male è appunto quella di partire dal match più corto, poi chiedere in locale di "allungare" il match vedendo se c'è la tua seq più lunga oppure no.
|
...e vicinissimo al mio sapere si accampa la tenebra della mia ignoranza... |
|
|
|
dallolio_gm
Moderatore
Prov.: Bo!
Città: Barcelona/Bologna
2445 Messaggi |
Inserito il - 27 maggio 2009 : 19:25:55
|
Credo che sia la dimostrazione del fatto che i metodi di allineamento non euristici sono troppo lenti per essere utilizzati :-(
Secondo me, potresti risolvere il problema creando delle matrici tra il genoma e le sequence da allineare, in questo modo ti sarebbe piu' facile tornare indietro nel caso di mismatch.
Ovvero, prova con dei dotplot: 
|
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 - 27 maggio 2009 : 20:19:00
|
ho provato a partire dalla quindicesima base e aggiungerle man mano.. ma la realtà è che il tempo nn cambia di molto, considerando che stiamo parlando di 50 giorni con un quad. Questo credo perchè molte sequenze mappano gia dopo aver tolto la prima base e nonostante questo va lento.
ho visto implementazioni con psyco.. devo ancora provare!! potrebbe portare ad una ottimizzazione uguale a quella del C! però nn so..
come faccio a creare delle matrici? questo mi porterà ad una maggiore velocità? |
|
|
|
Ra1D
Utente Junior
Prov.: Bergamo
Città: Treviglio
174 Messaggi |
Inserito il - 28 maggio 2009 : 15:22:58
|
In teoria esistono vari algoritmi per il pattern matching che riducono l'ordine di grandezza dei confronti, per esempio c'è l'algoritmo di Knuth-Morris-Pratt (KMP) (http://code.activestate.com/recipes/117214/ implementazione in python)
però è utile quando si debba mappare un pattern P su un testo T, trovando ogni posizione in cui il pattern compare. http://it.wikipedia.org/wiki/Algoritmo_di_Knuth-Morris-Pratt
E a quanto ho capito a te basta invece "mappare" nel senso di identificare le sequenze che sono presenti nel genoma, senza tener conto della loro posizione e di quante volte compaiono, giusto??
Oltretutto, piccola curiosità forse pure stupida, xkè devi partire dalla sequenza e togliere via via una base fino a lunghezza 15??
Se per dire la tua sequenza fosse:
C AAAAA AAAAA AAAAA AAAAA
Mentre nel tuo genoma hai
G AAAAA AAAAA AAAAA AAAAA
Tu toglieresti via via le ultime 5 "A" senza trovare mai alcun match (a causa della prima C diversa da G) mentre se partissi a diminuire la sequenza dall'inizio avresti al 2° giro un match con le 20 A.
Cmq il tempo di calcolo è enorme, anche con gli algoritmi migliori, per questo poi si sceglie di usare delle erustiche tipo fasta o blast.
E 50 giorni su un quad core è ancora buono eheheheheheh
|
...e vicinissimo al mio sapere si accampa la tenebra della mia ignoranza... |
|
|
|
dallolio_gm
Moderatore
Prov.: Bo!
Città: Barcelona/Bologna
2445 Messaggi |
Inserito il - 28 maggio 2009 : 15:49:22
|
Dai una occhiata a exonerate (piuttosto che riscrivere blast da capo).
E' giá compilato, quindi non avrai problemi di dipendenze; funziona da linea di comando, e non devi far altro che lanciare exonerate -q file_sequenze.fasta -t file_genome.fasta per ottenere risultati; e infine, avrai risultati piu' accurati, con la possibilità di inserire gaps e penalties. |
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 - 28 maggio 2009 : 17:50:20
|
Citazione:
E a quanto ho capito a te basta invece "mappare" nel senso di identificare le sequenze che sono presenti nel genoma, senza tener conto della loro posizione e di quante volte compaiono, giusto??
si, esatto.
Citazione:
Oltretutto, piccola curiosità forse pure stupida, xkè devi partire dalla sequenza e togliere via via una base fino a lunghezza 15??
ehh.. devo fare un controllo di questo tipo. Lasciando andare un po il script lento ho trovate tante sequenze che mappavano...
Citazione:
Dai una occhiata a exonerate (piuttosto che riscrivere blast da capo).
E' giá compilato, quindi non avrai problemi di dipendenze; funziona da linea di comando, e non devi far altro che lanciare exonerate -q file_sequenze.fasta -t file_genome.fasta per ottenere risultati; e infine, avrai risultati piu' accurati, con la possibilità di inserire gaps e penalties.
l'ho provato ma le mie sequenze non sono in formato fasta e quindi nn mi trova niente. Tra l'altro neanche il file genome è fasta ma .raw |
|
|
|
Ra1D
Utente Junior
Prov.: Bergamo
Città: Treviglio
174 Messaggi |
Inserito il - 28 maggio 2009 : 19:35:20
|
Ok allora direi che non ci sono molti modi per ottimizzare...
io farei così:
parto dalle sequenze di lunghezza minima (15), uso l'implementazione di KMP in python modificando il ciclo in modo che al primo hit si blocchi e prenda in considerazione la sequenza successiva.
Una volta trovate le sequenze da 15 che combaciano, provo partendo da dove era matchata ad allungare di 1 in locale partendo dal punto di match sul genoma, finchè non c'è il mismatch o la sequenza finisce, altrimenti salto alla sequenza successiva.
Questo penso sia uno dei pochi modi per evitare che il numero di confronti sia esponenziale rispetto alla dimensione testo e dimensione pattern.
Unico particolare, ti perdi eventuali match più lunghi se la sequenza è ripetuta. Tipo se trovi la tua sequenza da 15, mentre più avanti nel genoma avresti anche la tua sequenza chessò da 20...
Se vuoi beccare tutto, ti tocca scorrere tutto il genoma mi sa, una marea di volte x ogni sequenza, cioè tempo di calcolo esagerato!
Potresti chiedere in prestito qualche supercomputer in alternativa |
...e vicinissimo al mio sapere si accampa la tenebra della mia ignoranza... |
|
|
|
chick80
Moderatore
Città: Edinburgh
11491 Messaggi |
Inserito il - 28 maggio 2009 : 22:51:38
|
Citazione:
l'ho provato ma le mie sequenze non sono in formato fasta e quindi nn mi trova niente.
Non è molto più veloce riscrivere le sequenze in FASTA? Non sono un esperto ma non mi sembra ci sia poi così tanta differenza con quelle che hai postato, si farebbe abbastanza velocemente con delle regexp. |
Sei un nuovo arrivato?
Leggi il regolamento del forum e presentati qui
My photo portfolio (now on G+!) |
|
|
|
Discussione |
|
mappatura con python
Guide tools online
