Pitanje:
Koji je najbrži način za izračunavanje broja nepoznatih nukleotida u FASTA / FASTQ datotekama?
Kamil S Jaron
2017-06-01 21:47:03 UTC
view on stackexchange narkive permalink

Prije sam radio s javno dostupnim genomskim referencama, gdje su obično dostupne osnovne statistike, a ako nisu, morate ih izračunati samo jednom, tako da nema razloga za brigu o izvedbi.

Nedavno Počeo sam s projektom sekvenciranja nekoliko različitih vrsta sa srednjim genomima (~ Gbp) i tijekom ispitivanja različitih montažnih cjevovoda mnogo puta sam izračunao broj nepoznatih nukleotida i u neobrađenim očitanjima (u fastq) i u montažnim skelama (u fasta), stoga sam mislio da bih želio optimizirati izračunavanje.

  • Za mene je razumno očekivati ​​4-linijske formatirane datoteke fastq, ali opće rješenje je i dalje poželjno
  • To bilo bi lijepo kada bi rješenje djelovalo i na gzipirane datoteke

P: Koji je najbrži način (s obzirom na izvedbu) za izračunavanje broja nepoznatih nukleotida (Ns) u datotekama fasta i fastq ?

Vjerojatno biste mogli htjeti ukupan broj nukleotida kao i broj N (kako biste dobili udio Ns). Ovo je prilično brza implementacija koja pruža ovo plus nekoliko drugih mjernih podataka za fasta: https://github.com/ENCODE-DCC/kentUtils/blob/master/src/utils/faSize/faSize.c
Deset odgovori:
#1
+22
user172818
2017-06-02 02:46:23 UTC
view on stackexchange narkive permalink

Za FASTQ:

  seqtk fqchk in.fq | head -2  

Daje vam postotak "N" baza, a ne točan broj.

Za FASTA:

  seqtk comp in.fa | awk '{x + = $ 9} END {print x}'  

Ovaj redak za naredbe također radi s FASTQ-om, ali bit će sporiji jer je awk spor.

UREDI : ok, na osnovu podsjetnika @ BaCH-a, idemo (trebate kseq.h za kompajliranje):

  / / za kompajliranje: gcc -O2 -o count-N this-prog.c -lz # include <zlib.h> # include <stdio.h> # include <stdint.h> # include "kseqnaFd, gzqnaFd, gZQQQF, gzqnaFn, gseqnaF. [256] = {0, 1, 2, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4 , 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 4, 4 , 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 0, 4, 1, 4, 4, 4, 2, 4 , 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 0, 0 , 4, 1, 4, 4, 4, 2, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 3, 4, 4, 4, 4, 4, 4 , 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4}; int glavni (int argc, char * argv []) {long i, n_n = 0, n_acgt = 0, n_gap = 0; gzFile fp; kseq_t * seq; if (argc == 1) {fprintf (stderr, "Upotreba: count-N <in.fa> \ n"); povratak 1; } if ((fp = gzopen (argv [1], "r")) == 0) {fprintf (stderr, "POGREŠKA: nije uspjelo otvoriti ulaznu datoteku \ n"); povratak 1; } seq = kseq_init (fp);
while (kseq_read (seq) > = 0) {for (i = 0; i < seq->seq.l; ++ i) {int c = dna5tbl [(unsigned char) seq->seq.s [i]]; if (c < 4) ++ n_acgt; inače ako (c == 4) ++ n_n; ostalo ++ n_gap; }} kseq_destroy (seq); gzclose (fp); printf ("% ld \ t% ld \ t% ld \ n", n_acgt, n_n, n_prozor); return 0;}  

Radi i za FASTA / Q i za gzip'ed FASTA / Q. Sljedeće koristi SeqAn:

  #include <seqan / seq_io.h>using namespace seqan; int main (int argc, char * argv []) {if ( argc == 1) {std :: cerr << "Upotreba: count-N <in.fastq>" << std :: endl; povratak 1; } std :: ios :: sync_with_stdio (netačno); CharString id; Dna5String seq; SeqFileIn seqFileIn (argv [1]); dugo i, n_n = 0, n_acgt = 0; while (! atEnd (seqFileIn)) {readRecord (id, seq, seqFileIn); for (i = beginPosition (seq); i < endPosition (seq); ++ i) if (seq [i] < 4) ++ n_acgt; ostalo ++ n_n; } std :: cout << n_acgt << '\ t' << n_n << std :: endl; return 0;}  

Na FASTQ-u s 4 milijuna 150 bp glasi:

  • Verzija C: ~ 0,74 sek
  • The Verzija C ++: ~ 2,15 sek
  • Starija verzija C bez tablice pretraživanja (vidi prethodnu izmjenu): ~ 2,65 sek
na temelju benčmarkinga Devona Ryana najbrži ste do sada ...
Možete li dodati `#include ` na svoj prvi primjer C?
Ne mogu ne podržati SeqAn, ali ozbiljno sam zabrinut što je vaš C ++ kôd toliko sporiji od vašeg C koda. Trebali bi biti gotovo identični. To znači da niste naveli `std :: ios :: sync_with_stdio (false)`, što bi moglo objasniti neku razliku; a petlja bi se mogla pojednostaviti korištenjem iteratora, premda bih se iznenadio kad bi to utjecalo na brzinu bilo koje postavke iznad `-O1`.
@KonradRudolph `sync_with_stdio` malo pomaže pri čitanju iz normalne datoteke. Ne znam internu riječ seqan, pa ne mogu objasniti zašto je sporiji. To je reklo, ne bih se previše brinuo o izvedbi seqana. Gzip dekompresija je puno sporija od raščlanjivanja fastq. Dodan je @DevonRyan.
Razlika je u readfq. Ta zvijer je * brza *. Znam, pokušao sam zamijeniti Hengov kôd nečim manje ružnim u C ++-u i nisam se mogao približiti njemu u jednostavnim izvedbama čitanja. Najbolje što sam mogao smisliti u čistom C ++-u bio je 2x brži faktor brzine.
trebate dodati `-lz` u zastavice kompajliranja (barem za trenutni Ubuntu LTS)` gcc -O2 -o count-N this-prog.c -lz`
Temelje li se vaši brojevi na ispiranju predmemorije između uzoraka? Inače će ti brojevi gotovo sigurno biti netočni.
#2
+18
Devon Ryan
2017-06-02 03:20:21 UTC
view on stackexchange narkive permalink

5 sati i nema objavljenih mjerila? Jako sam razočaran.

Ograničit ću usporedbu na samo fasta datoteke, jer će fastq na kraju biti isti. Do sada su kandidati:

  1. R s paketom ShortRead (čak i ako ne najbrži, svakako super prikladan način).
  2. Cjevovod od grep -v "^ >" | tr -cd A | wc -c
  3. Cjevovod grep -v "^ >" | grep -io A | wc -l
  4. Cjevovod seqtk comp | awk
  5. Nešto prilagođeno u C / C ++

R kôd iz nekog razloga zapravo ne može izbrojati N zapisa, pa sam izbrojao A za sve.

Izuzet ću sva python rješenja, jer ne postoji zamislivi svemir u kojem je python brz. Za rješenje C / C ++ upravo sam upotrijebio ono što je objavio @ user172818 (sve što sam pokušao napisati trajalo je isto toliko vremena).

Ponavljajući 100 puta i uzimajući prosjek (da, treba uzeti srednju vrijednost):

  1. 1,7 sekundi
  2. 0,65 sekundi
  3. 15 sekundi
  4. 1,2 sekunde
  5. 0,48 sekundi

Nije iznenađujuće da će sve u pravom C ili C ++ pobijediti. grep s tr je iznenađujuće dobar, što me iznenađuje jer iako je grep vrlo optimiziran, ipak sam očekivao da će imati više troškova. Cjevovod do grep -io ubojica je performansi. R rješenje je iznenađujuće dobro s obzirom na tipične R. troškove.

Ažuriranje1: Kao što je predloženo u komentarima

  1. zbroj (x.broj ('A') za x u otvorenom ('seqs.fa', 'r') ako je x [0]! = '>') u pythonu
  2. ol >

    To rezultira vremenom od

    1. 1,6 sekundi (sada sam na poslu, pa sam prilagodio vrijeme najbolje što mogu kako bih uzeo u obzir različite računala)

    Ažuriranje2: Još mjerila iz komentara

    1. awk -FA '! / ^ > "/ {cnt + = NF-1} END {print cnt} 'foo.fa
    2. perl -ne 'if (! / ^ > /) {$ count + = tr / A //} END {print "$ count \ n"}' foo.fa

    Ova vremena prinosa od:

    1. 5,2 sekunde
    2. 1,18 sekundi

    Verzija awk je najbolji hak koji sam mogao smisliti, vjerojatno postoje bolji načini.

    Ažuriranje 3 : Tačno, dakle, u vezi s datotekama fastq, ovo pokrećem na 3.3 gzipirana datoteka s 10 ponavljanja (za pokretanje je potrebno malo vremena), tako da neću u početku ograničavati testove na naredbe koje trivijalno mogu izmijeniti za obradu komprimiranih datoteka (uostalom, nekomprimirane datoteke fastq grozne su) .

    1. seqtk fqchk
    2. sed -n '2 ~ 4p' < (zcat foo.fastq.gz) | grep -io A | wc -l
    3. sed -n '1d; N; N; N; P; d' < (zcat Undetermined_S0_R1_001.fastq.gz) | tr -cd A | wc -c
    4. Primjer C iz korisnika user1728181 (tamo je već pružena usporedba sa C ++, pa ga nema potrebe uključivati).
    5. bioawk -c fastx '{n + = gsub (/ A /, "", $ seq)} END {print n}' foo.fastq.gz

    Prosječna vremena su:

    1. 1 minuta 44 sekunde
    2. 2 minute 6 sekundi
    3. 1 minuta 52 sekunde
    4. 1 minuta 15 sekundi
    5. 3 minute 8 sekundi

    NAPOMENA, one često neće obrađivati ​​fastq datoteke s unosima duljim od 4 retka. Međutim, takve su datoteke dodo ptice bioinformatike, pa to praktički nije važno. Također, nisam uspio dobiti ShortRead za rad sa komprimiranom datotekom fastq. Pretpostavljam da postoji način da se to popravi.

Lijepo mjerilo, koja je veličina vaše testne fasta datoteke? Također, zašto ne biste dodali Python ako imate vremena, da neće biti najbrži, ali zanimljiv za usporedbu, pokušao bih: `sum (x.count ('A') for x in open ('seqs.fa' , 'r') ako x [0]! = '>') `, vjerojatno će` BioPython` biti sporiji za ovo
@Chris_Rands Koristio sam dm6 (nekomprimirani), tako da nije užasno velik. Pretpostavio sam da će python biti toliko očigledno sporiji da se nisam htio gnjaviti, ali pretpostavljam da ne bi škodilo dodavanju.
@Chris_Rands Ažurirao sam s python referentnom vrijednošću. Dolazi kao i R, što je bolje nego što bih pretpostavio.
Hvala! Nisam toliko iznenađen jer je pythonov `str.count` implementiran u C. Također smatram da bi trebalo postojati samo brzo` awk` rješenje, ali nisam u stanju to shvatiti
Perl: `` `perl -ne 'if (! / ^> /) {$ Count + = tr / A //} END {print" $ count \ n "}' dm6.fa```
U redu, sada kada je pitanje ažurirano, možete ih izmijeniti tako da rade za opći slučaj FASTQ (višeredni slijed + niz kvalitete, pri čemu se u redovima kvalitete mogu pojaviti '@').
Moj awk-fu dovoljan je da smislim djelotvoran primjer, ali nedovoljno dobar da bih smislio jedan koji nije bolno spor.
@gringer Eto, danas sam produktivan ...
Možete li isprobati moj `bioawk -c fastx '{n + = gsub (/ A /," ", $ seq)} END {print n}'`? Radoznao sam vidjeti kako se uspoređuje s ostalim rješenjima vaših podataka.
@bli: Da, trenutno imam 5 drugih ...
R kod može brojati Ns nakon što ga pravilno napišete. ;-) (Moja greška)
Inačice awk i perl ne daju iste rezultate na mojoj testnoj fasta datoteci: inačica awk pronalazi još jedan "A" u genomu C. elegans od verzije perl. Pretpostavljam da ovisi hoće li jedan od slijedova započeti ili završiti s "A".
Dodao sam fastq usporedbu s nekim metodama (uključujući bioawk).
@bli U awk kodu trebao bi biti samo `-1`. Među mnogim razlozima da se umjesto toga koristi `bioawk`!
Mislim da zcat nije potreban s bioawkom. Poboljšava li brzinu?
@bli Nije potreban i to je moja pogreška (nisam zapravo koristio zcat u referentnoj vrijednosti). Hvala što ste to ulovili!
Zašto računate `A`, a ne` N`? Ili mi nešto nedostaje?
Budući da originalna varijanta ShortRead nije mogla podnijeti N i usporedbe performansi su iste bez obzira. "R kod zapravo ne može iz nekog razloga izbrojiti N zapisa kad ga probam, pa sam za sve izbrojao A."
Da, oprostite, to mi je nedostajalo kad sam prvi put pročitao odgovor.
Kažete: "Ponavljanje 100 puta i uzimanje prosjeka (da, treba uzeti srednju vrijednost):" Ali prosjek je srednja vrijednost, a ne medijana - pa koju ste koristili? Cijenim da je medijan manje osjetljiv na iskrivljenje, ali u ovom se slučaju zasigurno vrijeme izvođenja normalno raspoređuje (u tom je slučaju srednja vrijednost u redu), osim ako ne koristite sustav s vrlo visokim trenutnim opterećenjem CPU-a ili IO-a od drugih korisnika, što stvara sukob .
Srednja vrijednost, a ne srednja vrijednost, dakle upozorenje. Da, vrijeme izvođenja treba biti normalno raspoređeno, ali povremeno se događaju neobični IO skokovi i takvi se događaju.
Jeste li ispraznili predmemoriju sustava između svih ovih testova?
Ne, u praksi su datoteke bile u predmemoriji za sva izvođenja.
Trebali biste isprazniti predmemoriju između testova. Vaši će brojevi zavarati, u suprotnom (dobro, pogrešno, stvarno) bez obzira na to crtate li srednju vrijednost ili medijanu iz rezultata ispitivanja.
#3
+9
Karel Brinda
2017-06-01 22:22:20 UTC
view on stackexchange narkive permalink

Mislim da bi ovo trebalo biti prilično brzo:

FASTA:

  grep -v "^ >" seqs.fa | tr -cd N | wc -c  

FASTQ:

  sed -n '1d; N; N; N; P; d' seqs.fq | tr -cd N | wc -c  

Pogledajte ovaj odgovor na SO o tome kako brojati znakove u BASH-u koristeći različite pristupe.

Vaš FASTQ odgovor nažalost u osnovi ima isti problem kao i Iakov. Odnosno, FASTQ zapisi imaju promjenjiv broj linija.
Znam da izvorni standard dopušta FASTQ cijepanje linija. Međutim, jeste li ikada vidjeli takvu datoteku u praksi?
Nisam, ali nikada osobno nisam radio s Sanger sekvenciranjem / podacima PacBio. Mislim da biste ih tamo mogli susresti. Za kratko pročitane podatke pretpostavljam da je pretpostavka "zapis = 4 retka" sigurna.
Dugo čitani podaci ne koriste fastq. Trenutno PacBio koristi bam (da, stvarno bam) i ONT fast5 (h5 arhiva). Mislim da je vrlo razumno pretpostaviti fastq od 4 retka.
#4
+8
Iakov Davydov
2017-06-01 22:00:26 UTC
view on stackexchange narkive permalink

Kao što je istaknuto, fastq može biti kompliciran. Ali u jednostavnom slučaju kada imate četiri retka po zapisu, jedno od mogućih rješenja u bashu je:

  sed -n '2 ~ 4p' seqs.fastq | grep -io N | wc -l  
  • sed -n '2 ~ 4p' ispisat će svaki četvrti redak
  • grep -o N će ispisati redak s N za svaki odgovarajući simbol
  • wc -l će brojati redove

Sumnjam da će ovaj pristup pythona raditi brže:

  cat seqs.fastq | python3 -c "import sys; print (sum (line.upper (). count ('N') for line in sys.stdin))"  

FASTA

Coreutils:

  grep -v "^ >" seqs.fasta | grep -io N | wc -l  

Python:

  mačka seqs.fasta | python3 -c "import sys; print (sum (line.upper (). count ('N') za redak u sys.stdin ako nije line.startswith ('>')))"  
Sada brojite retke koji sadrže Ns, a ne Ns. Također, FASTQ zapisi mogu obuhvaćati više od četiri retka. :-(
@KonradRudolph s -o grep daje pojedinačne simbole, razmislit ću kako ispraviti fastq kod
FASTQ je lukav jer `@` također može započeti kvalitetnu liniju niza; pogledajte http://chat.stackexchange.com/transcript/message/37809200#37809200 i sljedeće redove. Za referencu, ovdje je Perl skripta koja raščlanjuje FASTQ i uklanja podatke o spoju s identifikatora. Vjerojatno može poslužiti kao polazna točka; Mislim da je ovo vrlo blizu minimalnom kodu: https://gist.github.com/klmr/81a2b86cd93c706d28611f40bdfe2702
@KonradRudolph, Ažurirao sam odgovor ukazujući da to djeluje samo u jednostavnom slučaju
#5
+6
bli
2017-06-02 14:28:46 UTC
view on stackexchange narkive permalink

Korištenje bioawk

S bioawk (računajući "A" u genomu C. elegans , jer izgleda da u njima nema "N" ovu datoteku), na mom računalu:

   $ time bioawk -c fastx '{n + = gsub (/ A /, "", $  span> seq)} END {print n} 'genome.fa 32371810real 0m1.645suser 0m1.548ssys 0m0.088s  

bioawk je proširenje awka s prikladnim mogućnostima raščlanjivanja. Na primjer, -c fastx čini sekvencu dostupnom kao $ seq u formatu fasta i fastq ( uključujući gzipirane verzije ), pa je gornja naredba također bi trebao snažno rukovati formatom fastq .

Naredba gsub normalna je awk funkcija. Vraća broj zamijenjenih znakova (dobio ga je s https://unix.stackexchange.com/a/169575/55127). Pozdravljam komentare o tome kako brojati unutar awka na manje hakiran način.

Na ovom primjeru fasta gotovo je 3 puta sporiji od grep , tr , wc cjevovod:

  $ time grep -v "^ >" genome.fa | tr -cd "A" | wc -c32371810real 0m0.609suser 0m0.744ssys 0m0.184s  

(Ponovio sam mjerenja vremena, a gornje vrijednosti izgledaju reprezentativno)

Korištenje pythona

readfq

Isprobao sam python readfq preuzet iz spremišta spomenutog u ovom odgovoru: https://bioinformatics.stackexchange.com/a/365/292

  $ time python -c "iz sys import stdin; iz readfq import readfq; ispis zbroja ((seq.count ('A') for _, seq, _ in readfq (stdin)))" < genom .fa32371810real 0m0.976suser 0m0.876ssys 0m0.100s  

"čisti python"

Usporedba s nečim prilagođenim iz ovog odgovora: https: // bioinformatics. stackexchange.com/a/363/292

  $ time python -c "iz sys import stdin; ispisati zbroj ((line.count ('A') za redak u stdin ako ne i line.startswith ('>'))) "Genom <.fa32371810real 0m1.143user 0m1.100s
sys 0m0.040s  

(Sporši je s pythonom 3.6 nego s pythonom 2.7 na mom računalu)

pyGATB

Upravo sam naučio (08 / 06/2017) da GATB uključuje fasta / fastq parser i nedavno je objavio python API. Pokušao sam ju jučer iskoristiti za testiranje još jednog odgovora na sadašnje pitanje i pronašao sam grešku. Ova je greška sada ispravljena, pa evo odgovora na bazi pyGATB:

  $ time python3 -c "iz gatb import banke; print (sum ((seq. sequence.count (b \ "A \") za seq u Banci (\ "genome.fa \")))) "32371810real 0m0.663user 0m0.568ssys 0m0.092s  

( Možete i sequence.decode ("utf-8"). Count ("A") , ali to se čini malo sporije.)

Iako sam ovdje koristio python3.6 (pyGATB se čini samo python3), to je brže od druga dva python pristupa (za koja su izvještena vremena dobivena s python 2.7). Ovo je čak gotovo jednako brzo kao cjevovod grep , tr , wc .

Biopython

I, da bismo imali još više usporedbi, evo rješenja pomoću SeqIO.parse iz Biopythona (s python2.7):

  $ time python -c "iz Bio uvoz SeqIO; ispis (zbroj ((rec.seq.count (\ "A \") za rec u SeqIO.parse (\ "genome.fa \", \ "fasta \")))) "32371810real 0m1.632suser 0m1.532ssys 0m0.096s  

Ovo je nešto sporije od rješenja "čistog pythona", ali možda i robusnije.

Čini se da postoji malo poboljšanje sa Prijedlog @ peterjca za upotrebu niže razine SimpleFastaParser:

  $ time python -c "iz Bio.SeqIO.FastaIO import SimpleFastaParser; print (sum (seq.count ( 'A') za naslov, seq u SimpleFastaParser (open ('genome.fa'))))) "32371810real 0m1.618suser 0m1.500ssys 0m0.116s  

(Napravio sam niz vremena i pokušao uzeti jedan koji se činio reprezentativnim, ali ima puno preklapanja th vremena tajmera parsera više razine.)

scikit-bio

Danas sam pročitao o scikit-bio, koji ima čitač sekvenci.

  $ time python3 -c "import skbio; ispis (zbroj (seq.count ('A') za seq u skbio.io.read ('genome.fa', format = 'fasta', verify = False))) "32371810real 0m3.643user 0m3.440ssys 0m1.228s 

pyfastx

Danas sam pročitao o pyfastxu (27.8.2020.) koji ima čitač sekvenci.

  $ time python3 -c "iz pyfastx uvoza Fasta; ispis (zbroj (seq.count ('A') for (_, seq) u Fasta ('genome.fa', build_index = False)))" 32371810real 0m0.781suser 0m0.740ssys 0m0.040s  

(Vrijeme se vrši pomoću pythona 3.6, na istom stroju kao i prethodna mjerenja od prije 3 godine, pa bi trebao biti usporediv.)

fastx_read from mappy

Shvaćam (27.8.2020.) da nisam dodao mjerilo za fastx_read iz mappy koju redovito koristim već dosta dugo:

  $ time python3 -c "iz mappy import fastx_read; print (zbroj (seq.count ('A') for (_, seq, _) in fastx_read ('genome.fa'))) "32371810real 0m0.731user 0m0.648ssys 0m0.080s  

Mjerila na datoteci fastq.gz

Testirao sam neka od gornjih rješenja (ili njihove prilagodbe), računajući "N" u istoj datoteci koja je ovdje korištena: https: // bioinformatics .stackexchange.com / a / 400/292

bioawk

   $ time bioawk -c fastx '{ n + = gsub (/ N /, "", $  seq)} END {print n} 'SRR077487_2.filt.fastq.gz306072real 1m9.686suser 1m9.376ssys 0m0.304s  

pigz + readfq python modul

readfq se ne žali i vrlo je brz kad prođem izravno komprimirani fastq, ali vrati nešto pogrešno, zato ne zaboravite ručno se pobrinuti za dekompresiju.

Ovdje sam pokušao s pigz :

  $ time pigz -dc SRR077487_2.filt.fastq.gz | python -c "iz sys import stdin; iz readfq import readfq; ispis zbroja ((seq.count ('N') for _, seq, _ in readfq (stdin)))" 306072real 0m52.347suser 1m40.716ssys 0m8.604s 

Računalo ima 16 jezgri, ali pretpostavljam da je ograničavajući faktor za pigz čitanje s diska: procesori su daleko od pune brzine.

I s gzip:

  $ time gzip -dc SRR077487_2.filt.fastq.gz | python -c "iz sys import stdin; iz readfq uvozi readfq; ispiši zbroj ((seq.count ('N') for _, seq, _ in readfq (stdin)))" 306072real 0m49.448user 1m31.984ssys 0m2.312s 

Ovdje su i gzip i python koristili puni procesor. Bilans korištenja resursa bio je nešto bolji. Radim na stolnom računalu. Pretpostavljam da bi računalni poslužitelj bolje iskoristio pigz .

pyGATB

  $ time python3 -c "iz gatb import banke; print ( zbroj ((seq.sequence.count (b \ "N \") za seq u banci (\ "SRR077487_2.filt.fastq.gz \")))) "306072real 0m46.784suser 0m46.404ssys 0m0.296s  

biopython

  $ time python -c "iz gzip uvoza otvoren kao gzopen; iz Bio.SeqIO.QualityIO import FastqGeneralIterator; print (sum (seq.count (' N ') za (_, seq, _) u FastqGeneralIterator (gzopen (' SRR077487_2.filt.fastq.gz ')))) "" 306072real 3m18.103user 3m17.676ssys 0m0.428s  

pyfastx (dodano 27.08.2020.)

  $ time python3 -c "iz pyfastx import Fastq; print (sum (seq.count ('N') for (_, seq, _ ) u Fastqu ('SRR077487_2.filt.fastq.gz', build_index = False))) "306072real 0m51.140user 0m50.784ssys 0m0.352s  

fastx_read iz mappy (dodano 27/08/2020)

  $ time python3 -c "iz mappy import fastx_read; pri nt (zbroj (sek.broj ('N') za (_, sek, _) u fastx_read ('SRR077487_2.filt.fastq.gz'))) "" 306072real 0m52.336user 0m52.024ssys 0m0.300s  

Zaključak

pyGATB i bioawk obrađuju transparentno kompresiju (gzipirano ili ne) i razlike u formatu (fasta ili fastq). pyGATB prilično je nov alat, ali čini se učinkovitijim u usporedbi s ostalim python modulima koje sam testirao.

(Ažuriranje 27.8.2020.) pyfastx i mappy također se čine prilično prikladnima i samo su malo sporiji od pyGATB-a .

Možete li dodati brojeve FASTA parsera niskog nivoa Biopython? `` vrijeme python -c "iz Bio.SeqIO.FastaIO uvoz SimpleFastaParser; ispis (zbroj (sek.broj ('A') za naslov, sek u SimpleFastaParser (otvoren ('genome.fa'))))" "`
Zašto svi vaši primjeri broje slovo A, a ne N (i možda n) prema pitanju?
Brojim slovo "A" jer prva velika datoteka o kojoj sam razmišljao za svoje testiranje nije imala "N".
Mislim da će svi primjeri Pythona postići lošu ocjenu na ovoj maloj datoteci primjera zbog opterećenja učitavanja Pythona i uvoza biblioteka. Ali ako korisnik želi brojati samo jedan genom, to je svejedno fer test.
#6
+6
Matt Bashton
2017-06-03 01:28:27 UTC
view on stackexchange narkive permalink

Dekompresija gzipped FASTQ-a je glavno pitanje

Ako uzmemo gzipped FASTQ iz stvarnog svijeta (što bi, kako je predložio OP, bilo korisno), a ne trivijalnu FASTA kao polaznu točku, tada se stvarni problem zapravo dekomprimira datoteka koja ne broji N i u ovom slučaju program C count-N više nije najbrže rješenje.

Uz to bi bilo dobro koristiti određenu datoteku za usporedbu koja zapravo ima Ns, jer ćete dobiti neke prilično zanimljive razlike u vremenu izvršenja s nekim metodama računajući one koji se češće javljaju umjesto Ns.

Jedna od takvih datoteka je:

http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/HG00096/sequence_read/SRR077487_2.filt.fastq.gz

Također vrijedi provjeriti razne rješenja vraćaju točan odgovor, u gornjoj datoteci trebalo bi biti 306072 Ns.

Sljedeća napomena da je dekompresija ove datoteke preusmjerena na / dev / null sporija s zcat i gzip (koji su i gzip 1.6 na mom sustavu) nego recimo paralelna implementacija gzipa poput Marka Adlera pigz , za koji se čini da koristi 4 niti za dekompresiju. Sva mjerenja vremena predstavljaju prosječno 10 izvođenja koja izvještavaju stvarno (vrijeme zidnog sata).

  vrijeme pigz -dc SRR077487_2.filt.fastq.gz > / dev / nullreal 0m29.0132stime gzip -dc SRR077487_2. filt.fastq.gz > / dev / nullreal 0m40.6996s  

Između njih postoji ~ 11,7 sekundi razlike. Dalje, ako nakon toga pokušam usporediti referentnu liniju koja radi na FASTQ-u i daje točan odgovor (Napomena: još nisam naišao na FASTQ koji nije četveroredan, i ozbiljno tko generira ove datoteke!)

  vrijeme pigz -dc SRR077487_2.filt.fastq.gz | awk 'NR% 4 == 2 {print $ 1}' | tr -cd N | wc -c306072real 0m34.793s  

Kao što vidite, brojanje dodaje ~ 5,8 sekunde ukupnom vremenu izvođenja u odnosu na dekompresiju temeljenu na pigz . Uz to, ovaj put delta je veća kada se koristi gzip ~ 6,7 sekundi iznad same dekompresije gzip.

  vrijeme gzip -dc SRR077487_2.filt.fastq.gz | awk 'NR% 4 == 2 {print $ 1}' | tr -cd N | wc -c 306072real 0m44.399s  

Rješenje temeljeno na pigz awk tr wc je ~ 4,5 sekunde brže od count-N zasnovano na rješenju C koda:

  count count-N SRR077487_2.filt.fastq.gz 2385855128 306072 0real 0m39.266s  

Ova je razlika robusna za ponovno pokretanje koliko god puta želim. Pretpostavljam da biste mogli koristiti pthread u rješenju na bazi C ili ga izmijeniti da izvadite standard iz pigz-a, to bi također pokazalo povećanje performansi.

Benchmarking još jedna alternativna pigz grep čini se da treba manje-više isto vrijeme kao i varijanta temeljena na tr :

  vrijeme pigz -dc SRR077487_2.filt.fastq.gz | awk 'NR% 4 == 2 {print $ 1}' | grep -o N | wc -l306072real 0m34.869s  

Imajte na umu da je prethodno navedeno rješenje seqtk osjetno sporije,

  vrijeme seqtk comp SRR077487_2 .filt.fastq.gz | awk '{x + = $ 9} END {print x}' 306072real 1m42.062s 

Međutim, vrijedi napomenuti da seqtk comp radi malo više od ostalih rješenja .

#7
+5
Konrad Rudolph
2017-06-01 22:10:03 UTC
view on stackexchange narkive permalink

Iskreno, najlakši način (posebno za FASTQ) je vjerojatno korištenje posebne biblioteke za raščlanjivanje, poput R / Bioconductor:

  suppressMessages (library (ShortRead)) seq = readFasta (commandArgs (TRUE) [1]) # ili readFastqcat (colSums (alphabetFrequency (sread (seq))) [, 'N', drop = FALSE]), '\ n')  

Ovo može ne biti najbrži, ali je prilično brz , jer su relevantne metode visoko optimizirane. Najveći troškovi zapravo su vjerojatno od učitavanja ShortRead , što je neopravdani slogan.

#8
+5
BaCh
2017-06-01 22:15:23 UTC
view on stackexchange narkive permalink

Ako se bavite sirovom brzinom, tada je vjerojatno potrebno pisanje vlastitog malog C / C ++ programa.

Srećom, najgori dio (brz i pouzdan parser) ima već rješavan: readfq od Heng Lija vjerojatno je najbrži FASTA / FASTQ parser okolo.

I jednostavan je za upotrebu, primjer na GitHub-u lako se može proširiti kako bi se napravilo ono trebaš. Samo dodajte parametar parametara (za naziv datoteke koji želite analizirati), programirajte jednostavnu petlju 'N'-brojač i dajte rezultate ispisanim u stdout. Gotovo.

Inače, ako se radi o jednostavnosti, pogledajte Konradov odgovor za rješenje temeljeno na R. Ili vrlo jednostavna skripta koja koristi Biopython.

Bože, to je apsolutno grozni kod. ☹️ Mislim da je za neobrađenu brzinu [SeqAn] (http://docs.seqan.de/seqan/master/?p=SeqFileIn) (C ++) vjerojatno jednak i puno kvalitetniji.
Nije najljepše okolo. Pa opet, radi ono što treba i ne dolazi s cijelom knjižničnom ovisnošću. Imate li mjerilo za SeqAn i readfq?
Zanima me vidjeti kako se readfq i SeqAn uspoređuju s fasta / fastq parserom iz GATB-a: http://gatb-core.gforge.inria.fr/doc/api/classgatb_1_1core_1_1bank_1_1impl_1_1BankFasta.html#detailsNisam u mogućnosti kodirati u C / C ++ ipak.
#9
+3
Matt Shirley
2017-06-15 20:47:20 UTC
view on stackexchange narkive permalink

Za FASTA datoteke implementirao sam relativno učinkovitu metodu u pyfaidx v0.4.9.1. Ovaj me post natjerao da shvatim da je moj prethodni kôd bio prilično spor i jednostavan za zamjenu:

  $ pip install pyfaidx $ time faidx -i nukleotid ~ / Downloads / hg38.faname start end ATCGN otherschr1 1 248956422 67070277 67244164 48055043 48111528 18475410 chr10 1 133797422 38875926 39027555 27639505 27719976 534460 chr11 1 135086622 39286730 39361954 27903257 27981801 552880 chr11_KI270721v1_random 1 100316 18375 19887 31042 31012 0 ... chrX 1 156040895 46754807 46916701 30523780 30697741 1147866 CHRy 1 57227415 7886192 7956168 5285789 5286894 30812372 chrY_KI270740v1_random 1 37240 8274 12978 7232 8756 0 stvarni 2m28.251user 2m15.548ssys 0m9.951s  
#10
  0
Edward Kirton
2018-03-01 11:55:22 UTC
view on stackexchange narkive permalink

Ovdje su perl-lineri za datoteke fasta i fastq koji su dovoljno brzi ako želite samo ukupne vrijednosti:

$ perl -ne 'BEGIN {my $ n = 0}; / ^ > / ili $ n + = s / N // gi; END {print "N = $ n \ n"} 'contigs.fasta

Za svaki redak FASTA datoteke, ako nije zaglavlje (/ ^> /), to je red niza. Zamjena Ns, koja vraća broj izvršenih zamjena, koji se dodaje ukupnom iznosu.

$ zcat reads.fastq | perl -ne 'POČINI {moj $ n = 0}; / ^ [ATCGN] + $ / i i $ n + = s / N // gi; END {print "N = $ n \ n"} '

Za svaki redak sekvence FASTQ datoteke, podudarne sa / ^ [ATCGN] + $ / i (tj. Svi NT-ovi), zatim prebrojite zamjene izvršene na isti način.

Možete koristiti pigz umjesto zcat (ili gunzip -c) ako imate instaliran pigz.

Nekomprimirana datoteka fastq od 14M 250bp Illumina čitanja dovršena za 40 sek.

Ako biste mogli objasniti što čini svaki perl redak, to bi poboljšalo odgovor, posebno komentirajući razlike između oba pristupa


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