Domanda:
Come trovo sequenze identiche in un file FASTA?
Cleb
2017-11-10 21:28:12 UTC
view on stackexchange narkive permalink

Voglio creare un database per uno studio di proteomica. Pertanto, la mappatura da una data sequenza a un ID proteina deve essere univoca. Mi chiedo se esiste già una funzione incorporata in Biopython per questo, ma non sono riuscito a trovarne alcuna. Le sequenze duplicati devono essere unite in una singola voce e là rispettivi ID dovrebbero essere semplicemente concatenati (vedi risultato atteso sotto).

Con il seguente esempio di input

  >prot_id_1MAWIGLDISKLFVENSRDAAA>prot_id_2MENVSRIWRYQQRRVLSRAFTHWYLMGLTKHNHPS>prot_id_3MRTRPSRDFHYIVGRITRESPEEKGIFYVFH>prot_id_4MEMVLSSANPSTTFADSYVV>prot_id_5MAWIGLDISKLFVENSRDAAA>prot_id_6MENVSRIWRYQQRRVLSRAFTHWYLMGLTKHNHPS>prot_id_7MAWIGLDISKLFVENSRDAAA  

Attualmente lo faccio come segue:

  from Bio import SeqIOimport numpy as npfrom Bio.Seq import Seqfrom Bio.SeqRecord import SeqRecordfrom Bio.Alphabet import IUPAC # Leggi voci nel dizionario, gi: sequencerecord_dict = SeqIO.to_dict (SeqIO.parse ("test.fasta", "fasta")) # Estrai la sequenza come stringrecord_dict = {gi: str (record_dict [gi] .seq) per gi in record_dict} # Verifica la presenza di sequenze duplicate: q_id_dict = {} per k, v in record_dict.items (): seq_id_dict.setdefault (v, []). append (k) # Filtra le sequenze che appaiono più di una voltafilter_dict = {seqi: li per seqi, li in seq_id_dict.items () if len (li) > 1} # Unisci gli identificatori in un single onemerged_ids = {"_". join (ordinato (li)): seqi per seqi, li in filter_dict.items ()} # Rimuove tutte le voci associate a ID duplicati per gi in set (np.concatenate (filter_dict.values ​​( ))): record_dict.pop (gi, None) # Aggiungi di nuovo gli ID unitirecord_dict.update (merged_ids) # Crea un elenco di oggetti Seqfinal_seq = [SeqRecord (Seq (seqi, IUPAC.protein), id = gi, name = '' , description = '') for gi, seqi in record_dict.items ()] # Scrivi fileSeqIO.write (final_seq, 'test_output.fasta', 'fasta')  

che mi dà il output desiderato:

  >prot_id_4MEMVLSSANPSTTFADSYVV>prot_id_1_prot_id_5_prot_id_7MAWIGLDISKLFVENSRDAAA>prot_id_3MRTRPSRDFHYIVGRITRESPEEKGIFYVFH>prot_id_2_prot_id_6MENVSRIWRYQQRRVLSRAFTHWYLMGLTKHNHPS  

C'è un modo più intelligente di realizzare questo?

Sei risposte:
#1
+8
Bioathlete
2017-11-10 21:45:20 UTC
view on stackexchange narkive permalink

Il trucco sarebbe scambiare la chiave nel dizionario con la sequenza stessa.

Inoltre consiglierei di usare un separatore diverso, "_" poiché è quello che hanno gli ID correnti in modo che tu può facilmente separare i singoli ID dall'ID concatenato. Ho usato una pipa "|" in questo esempio. Inoltre ho scritto manualmente l'output FASTA invece di utilizzare Biopython per farlo.

  da Bio import SeqIOfrom collections import defaultdictdedup_records = defaultdict (list) per record in SeqIO.parse ("test.fasta" , "fasta"): # Usa la sequenza come chiave e poi ottieni una lista di id come valore dedup_records [str (record.seq)]. append (record.id) with open ("Output.fasta", 'w ') come output: per seq, ids in dedup_records.items (): # Unisciti agli id ​​e scrivili come output di fasta.write ("> {} \ n" .format (' | '.join (ids)) ) output.write (seq + "\ n")  
Sì, funziona bene ed è più condensato della mia soluzione (il tuo `dedup_records` è lo stesso di` seq_id_dict`). Ma è vero, tutti i filtri aggiuntivi non sono necessari; Venerdì pomeriggio... :)
(Aspetto con accettazione per un po 'solo per vedere se c'è già una sorta di built-in per questo)
Non esiste una funzione incorporata che io sappia che fa la concatenazione di ID. Se esistesse, probabilmente farebbe la stessa logica.
Sembra un buon approccio se inserire il file in memoria va bene, ma perché non usare biopython anche per scrivere? Questa asimmetria sembra strana, inoltre in questo modo non scaricherai tutta la sequenza su una singola riga
Soprattutto perché trovo che l'interfaccia BioPython sia ingombrante per scrivere sequenze in un formato fasta quando i dati non sono già un oggetto Seq. @cleb può usare Biopython se vuole che le modifiche principali siano nell'uso del defaultdict con la chiave scambiata.
Vedo, beh, considera ancora l'utilizzo di record.description e non record.id per assicurarti di ottenere tutte le intestazioni. Inoltre è possibile risparmiare memoria archiviando le sequenze come hash anziché le sequenze complete
OP si riferiva all'ID della proteina, quindi è quello che ho usato. Sono d'accordo che gli hash potrebbero essere utili per gestire grandi quantità di input.
@Chris_Rands: Mind aggiungendo la soluzione hash !? Nel mio caso ho solo gli ID, nessuna descrizione ulteriore ma se potessi aggiungere una risposta più dettagliata sarò più che felice di votarla.
@Cleb Volevo dire, il motivo per cui non ho aggiunto una soluzione hash è perché sarebbe difficile risolvere le collisioni hash (anche se sarebbero molto rare)
@Chris_Rands: Nessun problema, le soluzioni pubblicate sono attualmente sufficienti; sarebbe stato solo un tipo di cosa 'piacevole da vedere' ...
#2
+4
Cleb
2017-11-13 14:14:38 UTC
view on stackexchange narkive permalink

Proprio come una piccola variazione alla risposta di @ Bioathlete nel caso in cui desideri scrivere il fasta usando Biopython (ad esempio per aggiungere nomi e descrizione):

  da Bio import SeqIOfrom Bio.Seq import Seqfrom Bio.SeqRecord import SeqRecordfrom Bio.Alphabet import IUPACfrom collections import defaultdictdedup_records = defaultdict (list) per record in SeqIO.parse ("test.fasta", "fasta"): # Usa la sequenza come chiave e poi avere una lista di id come valore dedup_records [str (record.seq)]. append (record.id) # questo crea un generatore; se hai bisogno di un elenco fisico, sostituisci "(", ")" esterno con "[" e "]", rispettivamentefinal_seq = (SeqRecord (Seq (seqi, IUPAC.protein), id = "|" .join (gi) , name = '', description = '') for seqi, gi in dedup_records.items ()) # write fileSeqIO.write (final_seq, 'test_output_final.fasta', 'fasta')  

Questo darà il risultato desiderato:

  >prot_id_1 | prot_id_5 | prot_id_7MAWIGLDISKLFVENSRDAAA>prot_id_2 | prot_id_6MENVSRIWRYQQRRVLSRAFTHWYLMGLTKHNHPS>prot_id_4MEMVLSSANPSTTFADSYVV>prot_id_3MRTRPSRDFHYIVGRITRESPEEKGIFYVFH  
Penso che tu voglia dire che crea un generatore non un iteratore.
@Bioathlete: Ovviamente risolto, grazie per averlo individuato.
#3
+2
Pierre
2017-12-05 19:44:02 UTC
view on stackexchange narkive permalink

utilizzando datamash e alcuni awk

  $ cat input.fa | incolla - - | \ cut -c2- | sort -t $ '\ t' -k2,2 | \ datamash collapse 1 -g 2 | TR "" "_" | \ awk '{printf ( ">% s \ n% s \ n", $ 2, $ 1);}' >prot_id_1_prot_id_5_prot_id_7MAWIGLDISKLFVENSRDAAA>prot_id_4MEMVLSSANPSTTFADSYVV>prot_id_2_prot_id_6MENVSRIWRYQQRRVLSRAFTHWYLMGLTKHNHPS>prot_id_3MRTRPSRDFHYIVGRITRESPEEKGIFYVFH  
#4
+1
Tom
2017-12-05 19:14:46 UTC
view on stackexchange narkive permalink

Questo può anche essere fatto abbastanza facilmente con awk:

awk '{if (NR% 2) {printf ("% s \ t", $ 1)} else {print}} 'nomefile | sort -k2,2 | awk '{if ($ 2 == SEQ) {gsub (">", "", $ 1); ID = ID "|" $ 1} else {if (SEQ! = "") {printf ("% s \ n% s \ n ", ID, SEQ);} SEQ = $ 2; ID = $ 1;}} END {printf ("% s \ n% s \ n ", ID, SEQ)} '

Il primo awk converte il file fasta in un file separato da tabulazioni con formato ID \ tSequence , che viene quindi ordinato in sequenza da sort . L'ultimo awk passa attraverso il file ordinato guardando le sequenze: se la sequenza nella riga corrente è la stessa di quella nella riga precedente, aggiunge il suo id alla variabile ID. Altrimenti (cioè abbiamo incontrato una nuova sequenza), stampa i precedenti ID concatenati e la sequenza corrispondente.

#5
+1
Eslam Ibrahim
2019-04-18 18:06:16 UTC
view on stackexchange narkive permalink

Ecco il mio programma gratuito su Github Sequence database curator ( https://github.com/Eslam-Samir-Ragab/Sequence-database-curator)

È un programma molto veloce e può gestire:

  1. Sequenze nucleotidiche
  2. Sequenze proteiche

può funzionare con i sistemi operativi:

  1. Windows
  2. Mac
  3. Linux

Funziona anche per:

  1. Formato Fasta
  2. Formato Fastq

Cordiali saluti

Ciao Eslam, grazie per il tuo contributo. Sarebbe bello se potessi includere anche come utilizzare il tuo programma per ottenere ciò che il poster originale ha chiesto ... Qui cerchiamo di essere il più concreti possibile sia quando chiedi che rispondi :-).
Sono d'accordo con @KamilSJaron: Si prega di mostrare come utilizzare il pacchetto per risolvere il problema dall'alto (idealmente con lo stesso input). Ciò aiuterà gli altri a iniziare. Grazie!
#6
  0
Paulo Sergio Schlogl
2020-01-20 05:17:46 UTC
view on stackexchange narkive permalink

Se hai un file fasta molto grande (60 GB) perché usare SeqIO.parse e scrivere il file di output usando biopython?

Perché non usare invece SimpleFastaParser (ottieni una tupla (nome, seq) ) e non devi aspettare finché usi il parser e seqrecords.

Cosa ne pensate ragazzi?

Grazie

Paulo

Tipo di:

  from Bio.SeqIO.FastaIO import SimpleFastaParserfrom collections import defaultdictdedup_records = defaultdict (list) for (name, seq) in SimpleFastaParser (filename): # Usa la sequenza come chiave e quindi avere un elenco di ID come valore dedup_records [seq)]. append (name) with open (filename, 'w') as output: for seq, ids in dedup_records.items (): # Join the ids e scrivili come output di fasta.write ("> {} \ n" .format ('|' .join (ids))) output.write (seq + "\ n")  


Questa domanda e risposta è stata tradotta automaticamente dalla lingua inglese. Il contenuto originale è disponibile su stackexchange, che ringraziamo per la licenza cc by-sa 3.0 con cui è distribuito.
Loading...