Pitanje:
Očitajte distribuciju dužine iz datoteke FASTA
Scott Gigante
2017-05-17 09:38:52 UTC
view on stackexchange narkive permalink

Imam ~ 10 GB FASTA datoteke generirane iz pokretanja Miniona tvrtke Oxford Nanopore Technologies. Kako mogu brzo i učinkovito izračunati raspodjelu duljina čitanja?

Naivan pristup bio bi čitanje FASTA datoteke u Biopythonu, provjeravanje duljine svake sekvence, pohranjivanje duljina u numpirani niz i crtanje rezultati pomoću matplotlib-a, ali čini se da je to pronalazak kotača.

vidi https://www.biostars.org/p/72433/; https://www.biostars.org/p/45951/
@Pierre ova rješenja nisu dovoljna za dugo čitanje. Ako sam hej, izbacite jedan (tekst) redak po 1/10 baza, što bi dovelo do izlaza teksta više od 10.000 redaka (i potencijalno više od 10 puta) za dugo čitani fasta.
@ScottGigante jesu li to fasta datoteke od 10 GB ili niz datoteka čija ukupna količina iznosi 10 GB?
@amblina Imam jednu fasta datoteku kao izlaz iz poretools-a ili slično, s> 1M čitanja srednje dužine ~ 8Kb
Sedam odgovori:
#1
+17
Sam Nicholls
2017-05-17 15:43:50 UTC
view on stackexchange narkive permalink

Ako želite nešto brzo i prljavo, možete brzo indeksirati FASTA pomoću samtools faidx , a zatim u naredbeni redak staviti stupac dužina kroz R (dostupni su i drugi jezici).

pre> samtools faidx $ fastacut -f2 $ fasta.fai | Rscript -e 'podaci <- as.numeric (readLines ("stdin")); sažetak (podaci); hist (podaci) '

Ovo daje statistički sažetak i stvara PDF u trenutnom direktoriju zvanom Rplots.pdf, koji sadrži histogram.

Možete i jednostavnije s nečim poput: `faidx --transform chromsizes | izrezati -f2 | Rscript -e 'podaci <- as.numeric (readLines ("stdin")); sažetak (podaci); hist (podaci) ''. To zahtijeva pyfaidx: `pip install pyfaidx`.
#2
+12
gringer
2017-05-17 22:33:20 UTC
view on stackexchange narkive permalink

Statistika očitavanja nanopore je nezgodna zbog ogromnog raspona dužina čitanja koje mogu biti prisutne u jednom izvođenju. Otkrio sam da je najbolji način prikaza duljina korištenje skale dnevnika i na x osi (duljina) i na y osi (sekvencirane baze ili brojevi, ovisno o preferencijama).

Napisao sam moje vlastite skripte za to: jedno za generiranje duljina čitanja i drugo za crtanje raspodjele duljina na razne načine. Skripta koja generira duljine čitanja također ispljuva statistiku sažetka osnovne duljine na standardnu ​​pogrešku:

  $ ~ / scripts / fastx-length.pl > lengths_mtDNA_called.txtUkupne sekvence: 2110Ukupna duljina: 5.106649 MbNajduži slijed: 107.414 kbNajkraća sekvenca: 219 bSrednja dužina: 2,42 kbMedijalna dužina: 1,504 kbN50: 336 sekvenci; L50: 3.644 kbN90: 1359 sekvenci; L90: 1.103 kb $ ~ / scripts / length_plot.r lengths_mtDNA_called.txtlengths_mtDNA_called.txt ... gotovoBroj sekvenci: 2110 Duljina kvantila: 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 

Evo nekoliko izrađenih grafikona:

Digital electrophoresis plot

Read length distribution plot

Skripte za njihovo generiranje možete pronaći ovdje:

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

Korištenje Biopythona i matplotliba doista bi izgledalo kao put. Zapravo se svodi na tri retka koda da bi se dobio taj graf:

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

Naravno da biste mogli napraviti dužu skriptu koja se može pozvati iz naredbenog retka, s nekoliko opcija. Dobrodošli ste da koristite moj:

  #! / Usr / bin / env python2 "" "Prilagođena skripta za crtanje distribucije dužina u fasta datoteci. Napisao Lucas Sinclair.Kopimi. Ovu skriptu iz ljuske možete koristiti ovako: $ ./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) iz matplotli # python #################################################### ############################## desc = "fasta_length_hist v1.0" parser = argparse.ArgumentParser (description = desc, formatter_class = RawTextHelpFormatter ) # Svi potrebni argumenti # parser.add_argument ("- input", help = "Fasta datoteka za obradu", type = str) parser.add_argument ("- out", type = str) # Svi neobavezni argumenti # parser.add_argument ("- x_log", zadana vrijednost = Tačno, tip e = bool) parser.add_argument ("- y_log", default = True, type = bool) # Analizirajte ga #args = parser.parse_args () input_path = args.inputoutput_path = args.outx_log = bool (args.x_log) y_log = bool (args.y_log) ############################################# ####################################### Pročitajte #lengths = map (len, SeqIO.parse ( input_path, 'fasta')) # Report # sys.stderr.write ("Pročitajte sve dužine (% i sekvence) \ n"% len (lengths)) sys.stderr.write ("Najdulji niz:% i bp \ n" % max (duljine)) sys.stderr.write ("Najkraći niz:% i bp \ n"% min (duljine))
sys.stderr.write ("Izrada grafa ... \ n") # Podaci #values ​​= pandas.Series (duljine) # Plot #fig = pyplot.figure () osi = values.hist (boja = 'siva', kante = 1000) fig = pyplot.gcf () title = 'Raspodjela duljina sekvenci' osi.set_title (naslov) axes.set_xlabel ('Broj nukleotida u nizu') osi.set_ylabel ('Broj sekvenci s ovom duljinom') osi .xaxis.grid (False) # Log #if x_log: axes.set_yscale ('symlog') if y_log: axes.set_xscale ('symlog') # Adjust # width = 18.0; visina = 10,0; dno = 0,1; vrh = 0,93; lijevo = 0,07; desno = 0,98fig.set_figwidth (širina) fig.set_figheight (visina) fig.subplots_adjust (hspace = 0,0, dno = dno, gore = gore, lijevo = lijevo, desno = desno) # Podaci i izvor # fig.text (0,99, 0,98, time.asctime (), horizontalalignment = 'desno') fig.text (0,01, 0,98, 'korisnik:' + getpass.getuser (), horizontalalignment = 'left') # Lijepo grupiranje znamenki #sep = ('x' , 'y') ako je 'x' u sep: locale.setlocale (locale.LC_ALL, '') seperate = lambda x, pos: locale.format ("% d", x, grouping = True) axes.xaxis.set_major_formatter (matplotlib.ticker.FuncFormatter (odvojeno)) ako je 'y' u rupi: locale.setlocale (locale.LC_ALL, '') seperate = lambda x, pos: locale.format ("% d", x, grupiranje = True) axes.yaxis.set_major_formatter (matplotlib.ticker.FuncFormatter (odvojeno)) # Spremi # fig.savefig (output_path, format = 'pdf')  

UREDI - primjer rezultata:

Sequence length distribution

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

Postoji nekoliko potencijalnih pristupa. Na primjer:

Koji su od njih " brzo i učinkovito "pomoću datoteke od 10 GB ... teško je reći unaprijed. Možda ćete morati pokušati usporediti neke od njih.

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

bioawk mogao bi biti relativno učinkovit za ovu vrstu zadatka.

  $ bioawk -c fastx '{histo [length ($ seq)] ++} END {for (l in histo) print l, histo [l]} '\ | vrsta -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 govori programu da raščlani podatke kao fastq ili fasta. To daje pristup različitim dijelovima zapisa kao $ name , $ seq (i $ qual u slučaju fastq formata) u awk kôd (bioawk se temelji na awk, tako da možete koristiti sve jezične značajke koje želite od awk).

Između pojedinačnih navodnika dolazi niz <condition> {<action> } blokova.

Prvi nema <condition> dijela, što znači da se izvršava za svaki zapis. Ovdje ažurira brojanje duljina u tablici koju sam nazvao "histo". length unaprijed je definirana funkcija u awk-u.

U drugom bloku, uvjet END znači da želimo da se izvrši nakon što se izvrši sav unos obrađena. Akcijski dio sastoji se u petlji preko zabilježenih vrijednosti duljine i ispisu ih zajedno s pripadajućim brojanjem.

Izlaz se cijevi u sort -n kako bi se rezultati numerički sortirali.

Na mojoj radnoj stanici, gore navedenom kodu trebalo je 20 sekundi da se izvrši za 1.2G fasta datoteku.

Shvaćam da izlaz neće biti prikladan kada se radi s vrijednostima rijetke duljine, jer nema spajanja (ili, što je ekvivalentno tome, spremnici imaju širinu 1).
#6
+3
Mark Ebbert
2017-06-24 22:29:15 UTC
view on stackexchange narkive permalink

Izričito ste pitali za FASTA datoteke, ali važno je uvijek uzeti u obzir dužinu i kvalitetu čitanja prilikom procjene dugo pročitanih podataka s velikom pogreškom. FASTA datoteke ne pružaju kvalitetu. Ove će vam informacije pomoći da utvrdite koliko je uspješno trčanje bilo, koliko je čitanja bilo "visokokvalitetno" itd.

Ovdje sam izvorno objavio cjelovit odgovor, predlažući pauvre, ali zaključio sam da je to malo zaobišlo temu jer izgleda da imate samo FASTA datoteke.

Preporučujem generiranje FASTQ datoteka, ali nisam siguran imate li izvornu bazu -pozvane fast5 datoteke. Ako je tako, generirajte FASTQ-ove koristeći poretools kako slijedi ( poretools doc za generiranje FASTQ datoteka):

  poretools fastq fast5 /  

Tada preporučujem generiranje toplotne karte i grafikona margine histograma pomoću pauvre s dužinom čitanja i kvalitetom kao što je prikazano u nastavku.

My description

[Napomena: Ja nisam originalni autor za pauvre, ali sada doprinosim ovom projektu]

Bio bih zahvalan na bilo kakvim povratnim informacijama o tome zašto moj odgovor može biti nezadovoljan.
Dragi Mark, imam problema s pauvreom. Možeš li mi pomoći?
#7
+1
H. Gourlé
2017-05-17 11:38:50 UTC
view on stackexchange narkive permalink

Nije točno ono što ste tražili, ali možete generirati histogram distribucije duljine čitanja vaših podataka o nanoporima izravno iz HDF5 datoteka pomoću poretools

Zbog ograničenja prostora na disku, zadano ponašanje osnovnih pozivača nanopore više ne stvara datoteke FAST5 (HDF5) kao izlaz. Iako ih zapravo imam, većina korisnika to neće, a uz to se ova metoda ne bi generalizirala na druge sekvencere poput PacBio.


Ova pitanja su automatski prevedena s engleskog jezika.Izvorni sadržaj dostupan je na stackexchange-u, što zahvaljujemo na cc by-sa 3.0 licenci pod kojom se distribuira.
Loading...