Domanda:
Leggere la distribuzione della lunghezza dal file FASTA
Scott Gigante
2017-05-17 09:38:52 UTC
view on stackexchange narkive permalink

Ho un file FASTA da ~ 10 GB generato da un'esecuzione MinION di Oxford Nanopore Technologies. Come posso calcolare in modo rapido ed efficiente la distribuzione delle lunghezze di lettura?

Un approccio ingenuo sarebbe leggere il file FASTA in Biopython, controllare la lunghezza di ogni sequenza, memorizzare le lunghezze in un array numpy e tracciare il risultati usando matplotlib, ma questo sembra reinventare la ruota.

vedere https://www.biostars.org/p/72433/; https://www.biostars.org/p/45951/
@Pierre queste soluzioni sono inadeguate per lunghe letture. Se eseguo l'output di una singola riga (di testo) per basi 1/10, ciò porterebbe a un output di testo superiore a 10.000 righe (e potenzialmente più di 10 volte) per una lunga lettura fasta.
@ScottGigante sono questi fasta in un file da 10 GB o in una serie di file per un totale di 10 GB?
@amblina Ho un singolo file fasta come output da poretools o simili, con> 1M di letture di lunghezza media ~ 8Kb
Sette risposte:
#1
+17
Sam Nicholls
2017-05-17 15:43:50 UTC
view on stackexchange narkive permalink

Se vuoi qualcosa di veloce e sporco, potresti indicizzare rapidamente FASTA con samtools faidx e quindi inserire la colonna delle lunghezze attraverso R (sono disponibili altre lingue) sulla riga di comando.

  samtools faidx $ fastacut -f2 $ fasta.fai | Rscript -e 'data <- as.numeric (readLines ("stdin")); riepilogo (dati); hist (data) ' 

Questo produce un riepilogo statistico e crea un PDF nella directory corrente chiamata Rplots.pdf, contenente un istogramma.

Puoi andare ancora più semplice con qualcosa come: `faidx --transform chromsizes | cut -f2 | Rscript -e 'data <- as.numeric (readLines ("stdin")); riepilogo (dati); hist (dati) ''. Ciò richiede pyfaidx: `pip install pyfaidx`.
#2
+12
gringer
2017-05-17 22:33:20 UTC
view on stackexchange narkive permalink

Le statistiche per le letture dei nanopori sono complicate a causa dell'enorme gamma di lunghezze di lettura che possono essere presenti in una singola esecuzione. Ho scoperto che il modo migliore per visualizzare le lunghezze è usare una scala logaritmica sia sull'asse x (lunghezza) che sull'asse y (basi sequenziate o conteggi, a seconda delle preferenze).

Ho scritto i miei script per fare questo: uno per generare le lunghezze di lettura e un altro per tracciare la distribuzione della lunghezza in vari modi. Lo script che genera le lunghezze di lettura restituisce anche statistiche di riepilogo della lunghezza di base all'errore standard:

  $ ~ / scripts / fastx-length.pl > lengths_mtDNA_called.txt Sequenze totali: 2110 Lunghezza totale: 5.106649 Mb Sequenza più lunga: 107,414 kb Sequenza più breve: 219 bLunghezza media: 2,42 kbLunghezza media: 1,504 kbN50: 336 sequenze; L50: 3,644 kb N90: 1359 sequenze; L90: 1.103 kb $ ~ / scripts / length_plot.r lengths_mtDNA_called.txtlengths_mtDNA_called.txt ... doneNumero di sequenze: 2110 Quantili di lunghezza: 0% 10% 20% 30% 40% 50% 60% 70% 219,0 506,9 724,4 953,0 1196,2 1503,0 1859,2 2347,3 80% 90% 100% 3128,2 4804,7 107414,0 

Ecco un paio di grafici prodotti:

Digital electrophoresis plot

Read length distribution plot

Gli script per generarli possono essere trovati qui:

#3
+7
xApple
2017-05-17 12:32:21 UTC
view on stackexchange narkive permalink

Usare Biopython e matplotlib sembrerebbe la strada da percorrere, anzi, si riduce a tre righe di codice per ottenere quel grafico:

  import Bio, pandaslengths = map (len , Bio.SeqIO.parse ('/ path / to / the / seqs.fasta', 'fasta')) pandas.Series (lengths) .hist (color = 'gray', bins = 1000)  

Ovviamente potresti voler creare uno script più lungo richiamabile dalla riga di comando, con un paio di opzioni. Puoi usare il mio:

  #! / Usr / bin / env python2 "" "Uno script personalizzato per tracciare la distribuzione delle lunghezze in un file fasta. Scritto da Lucas Sinclair.Kopimi. Puoi utilizzare questo script dalla shell in questo modo: $ ./fastq_length_hist --input seqs.fasta --out seqs.pdf "" "################### ################################################ ########## Moduli #import argparse, sys, time, getpass, localefrom argparse import RawTextHelpFormatterfrom Bio import SeqIOimport pandas # Matplotlib #import matplotlibmatplotlib.use ('Agg', warn = False) da matplotlib import pyplot # ################################################ ############################# desc = "fasta_length_hist v1.0" parser = argparse.ArgumentParser (description = desc, formatter_class = RawTextHelpFormatter ) # Tutti gli argomenti richiesti # parser.add_argument ("- input", help = "Il file fasta da elaborare", type = str) parser.add_argument ("- out", type = str) # Tutti gli argomenti opzionali # parser.add_argument ("- x_log", default = True, typ e = bool) parser.add_argument ("- y_log", default = True, type = bool) # Parse it #args = parser.parse_args () input_path = args.inputoutput_path = args.outx_log = bool (args.x_log) y_log = bool (args.y_log) ######################################### ##################################### Leggi #lengths = map (len, SeqIO.parse ( input_path, 'fasta')) # Report # sys.stderr.write ("Leggi tutte le lunghezze (% i sequenze) \ n"% len (lengths)) sys.stderr.write ("Sequenza più lunga:% i bp \ n" % max (lunghezze)) sys.stderr.write ("Sequenza più breve:% i bp \ n"% min (lunghezze))
sys.stderr.write ("Making graph ... \ n") # Data #values ​​= pandas.Series (lengths) # Plot #fig = pyplot.figure () axes = values.hist (color = 'gray', bins = 1000) fig = pyplot.gcf () title = 'Distribuzione delle lunghezze di sequenza'axes.set_title (titolo) axes.set_xlabel (' Numero di nucleotidi in sequenza ') axes.set_ylabel (' Numero di sequenze con questa lunghezza ') assi .xaxis.grid (False) # Log #if x_log: axes.set_yscale ('symlog') if y_log: axes.set_xscale ('symlog') # Regola # larghezza = 18.0; altezza = 10,0; in basso = 0,1; in alto = 0,93; sinistra = 0,07; right = 0.98fig.set_figwidth (width) fig.set_figheight (height) fig.subplots_adjust (hspace = 0.0, bottom = bottom, top = top, left = left = left, right = right) # Dati e sorgente # fig.text (0.99, 0.98, time.asctime (), horizontalalignment = 'right') fig.text (0.01, 0.98, 'user:' + getpass.getuser (), horizontalalignment = 'left') # Nice digit grouping #sep = ('x' , 'y') se 'x' in settembre: locale.setlocale (locale.LC_ALL, '') seperate = lambda x, pos: locale.format ("% d", x, grouping = True) axes.xaxis.set_major_formatter (matplotlib.ticker.FuncFormatter (seperate)) if 'y' in sep: locale.setlocale (locale.LC_ALL, '') seperate = lambda x, pos: locale.format ("% d", x, grouping = True) axes.yaxis.set_major_formatter (matplotlib.ticker.FuncFormatter (seperate)) # Save it # fig.savefig (output_path, format = 'pdf')  

EDIT - un output di esempio:

Sequence length distribution

#4
+5
neilfws
2017-05-17 10:09:24 UTC
view on stackexchange narkive permalink

Esistono diversi potenziali approcci. Ad esempio:

su quali di queste sono " veloce ed efficiente "utilizzando un file da 10 GB ... è difficile dirlo in anticipo. Potrebbe essere necessario provare a eseguire il benchmarking di alcuni di essi.

#5
+4
bli
2017-05-18 14:23:18 UTC
view on stackexchange narkive permalink

bioawk potrebbe essere ragionevolmente efficiente per questo tipo di attività.

  $ bioawk -c fastx '{histo [length ($ seq)] ++} END {for (l in histo) print l, histo [l]} '\ | sorta -n0 332701 15422 11323 33974 87765 118846 124747 143418 131659 1546710 2108911 3046912 4520413 6231114 8874415 11576716 14077017 19181018 31308819 51811120 109786721 472971522 657555723 273406224 101547625 49332326 32382727 16441928 10712029 7248730 4012031 2453832 2229533 1312134 938235 484736 385837 316138 285239 238840 163941 96142 68643 37744 19945 11446 7847 5948 5049 5250 4851 4252 3953 2854 4955 5956 5557 5158 5559 4360 5261 5662 4863 6764 9565 488  

-c fastx dice al programma di analizzare i dati come fastq o fasta. Questo dà accesso alle diverse parti dei record come $ name , $ seq (e $ qual nel caso del formato fastq) nel codice awk (bioawk è basato su awk, quindi puoi usare qualsiasi caratteristica del linguaggio desideri da awk).

Tra le virgolette singole c'è una serie di <condition> {<action> } blocchi.

Il primo non ha una parte <condition> , il che significa che viene eseguito per ogni record. Qui aggiorna i conteggi delle lunghezze in una tabella che ho chiamato "histo". length è una funzione predefinita in awk.

Nel secondo blocco, la condizione END significa che vogliamo che venga eseguita dopo che tutto l'input è stato elaborato. La parte dell'azione consiste nel scorrere in loop i valori di lunghezza registrati e stamparli insieme al conteggio associato.

L'output viene reindirizzato a sort -n per ordinare i risultati numericamente.

Sulla mia workstation, il codice sopra ha richiesto 20 secondi per essere eseguito per un file fasta 1.2G.

Mi rendo conto che l'output non sarà conveniente quando si tratta di valori di lunghezza sparsi, perché non c'è raggruppamento (o, equivalentemente, i contenitori hanno larghezza 1).
#6
+3
Mark Ebbert
2017-06-24 22:29:15 UTC
view on stackexchange narkive permalink

Hai chiesto specificamente i file FASTA, ma è importante considerare sempre la lunghezza e la qualità della lettura congiuntamente quando si valutano i dati di lettura lunga ad alto errore. I file FASTA non forniscono la qualità. Queste informazioni ti aiuteranno a determinare il successo della corsa, il numero di letture "di alta qualità", ecc.

All'inizio avevo pubblicato una risposta completa qui, suggerendo pauvre, ma ho deciso che era un po 'fuori tema visto che sembra che tu abbia solo i file FASTA.

Consiglio di generare file FASTQ, ma non sono sicuro che tu abbia la base originale -chiamato file fast5. In tal caso, genera FASTQ utilizzando poretools come segue ( poretools doc per la generazione di file FASTQ):

  poretools fastq fast5 /  

Quindi consiglio di generare una mappa termica e un grafico del margine dell'istogramma utilizzando pauvre con lunghezza e qualità di lettura come mostrato di seguito.

My description

[Nota: non sono l'autore originale di pauvre, ma ora contribuisco a questo progetto]

Apprezzerei qualsiasi feedback sul motivo per cui la mia risposta potrebbe essere sgradevole.
Caro Mark, ho problemi con il pauvre. Mi potete aiutare?
#7
+1
H. Gourlé
2017-05-17 11:38:50 UTC
view on stackexchange narkive permalink

Non è esattamente quello che hai chiesto, ma puoi generare un istogramma della distribuzione della lunghezza di lettura dei tuoi dati nanopori direttamente dai file HDF5 utilizzando poretools

A causa dei limiti di spazio su disco, il comportamento predefinito per i richiamo di base nanopori non crea più i file FAST5 (HDF5) come output. Anche se in realtà li ho, la maggior parte degli utenti no, e inoltre questo metodo non sarebbe generalizzato ad altri sequencer come PacBio.


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...