Frage:
Was ist der schnellste Weg, um die Anzahl unbekannter Nukleotide in FASTA / FASTQ-Dateien zu berechnen?
Kamil S Jaron
2017-06-01 21:47:03 UTC
view on stackexchange narkive permalink

Früher habe ich mit öffentlich verfügbaren Genomreferenzen gearbeitet, bei denen normalerweise grundlegende Statistiken verfügbar sind. Wenn dies nicht der Fall ist, müssen Sie diese nur einmal berechnen, damit Sie sich keine Sorgen um die Leistung machen müssen.

In letzter Zeit Ich begann mit dem Sequenzierungsprojekt einiger verschiedener Arten mit mittelgroßen Genomen (~ Gbp) und hatte beim Testen verschiedener Assemblierungspipelines viele Male die Anzahl unbekannter Nukleotide sowohl in Rohdaten (in Fastq) als auch in Assemblierungsgerüsten (in Fasta) berechnet. Daher dachte ich, dass ich die Berechnung optimieren möchte.

  • Für mich ist es vernünftig, 4-zeilig formatierte Fastq-Dateien zu erwarten, aber eine allgemeine Lösung wird immer noch bevorzugt
  • wäre schön, wenn die Lösung auch für komprimierte Dateien funktionieren würde

F: Was ist der schnellste Weg (in Bezug auf die Leistung), um die Anzahl unbekannter Nukleotide (Ns) in Fasta- und Fastq-Dateien zu berechnen ?

Vermutlich möchten Sie möglicherweise die Gesamtzahl der Nukleotide sowie die Anzahl der Ns (um den Anteil der Ns anzugeben). Dies ist eine ziemlich schnelle Implementierung, die dies und einige andere Metriken für fasta bereitstellt: https://github.com/ENCODE-DCC/kentUtils/blob/master/src/utils/faSize/faSize.c
Zehn antworten:
#1
+22
user172818
2017-06-02 02:46:23 UTC
view on stackexchange narkive permalink

Für FASTQ:

  seqtk fqchk in.fq | head -2  

Gibt den Prozentsatz der "N" -Basen an, jedoch nicht die genaue Anzahl.

Für FASTA:

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

Diese Befehlszeile funktioniert auch mit FASTQ, ist jedoch langsamer, da awk langsam ist.

BEARBEITEN : ok, basierend auf der Erinnerung von @ BaCH, los geht's (zum Kompilieren benötigen Sie kseq.h):

  / / zum Kompilieren: gcc -O2 -o count-N this-prog.c -lz # include <zlib.h> # include <stdio.h> # include <stdint.h> # include "kseq.h" [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, 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, 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, 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 main (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, "Verwendung: count-N <in.fa> \ n"); return 1; } if ((fp = gzopen (argv [1], "r")) == 0) {fprintf (stderr, "FEHLER: Eingabedatei kann nicht geöffnet werden \ n"); return 1; } seq = kseq_init (fp);
while (kseq_read (seq) > = 0) {für (i = 0; i < seq->seq.l; ++ i) {int c = dna5tbl [(vorzeichenloses Zeichen) seq->seq.s [i]]; if (c < 4) ++ n_acgt; sonst wenn (c == 4) ++ n_n; sonst ++ n_gap; }} kseq_destroy (seq); gzclose (fp); printf ("% ld \ t% ld \ t% ld \ n", n_acgt, n_n, n_gap); return 0;}  

Es funktioniert sowohl für FASTA / Q als auch für gzip'ed FASTA / Q. Im Folgenden wird SeqAn verwendet:

  #include <seqan / seq_io.h>using Namespace seqan; int main (int argc, char * argv []) {if ( argc == 1) {std :: cerr << "Verwendung: count-N <in.fastq>" << std :: endl; return 1; } std :: ios :: sync_with_stdio (false); CharString id; Dna5String seq; SeqFileIn seqFileIn (argv [1]); lang i, n_n = 0, n_acgt = 0; while (! atEnd (seqFileIn)) {readRecord (id, seq, seqFileIn); für (i = beginPosition (seq); i < endPosition (seq); ++ i) if (seq [i] < 4) ++ n_acgt; sonst ++ n_n; } std :: cout << n_acgt << '\ t' << n_n << std :: endl; return 0;}  

Auf einem FASTQ mit 4 Millionen 150 bp lautet:

  • Die C-Version: ~ 0,74 Sek.
  • Die C ++ - Version: ~ 2,15 Sek.
  • Eine ältere C-Version ohne Nachschlagetabelle (siehe die vorherige Bearbeitung): ~ 2,65 Sek.
Basierend auf Devon Ryans Benchmarking sind Sie der bisher schnellste ...
Können Sie Ihrem ersten C-Beispiel `#include ` hinzufügen?
Ich kann SeqAn nicht positiv bewerten, aber ich bin ernsthaft besorgt, dass Ihr C ++ - Code so viel langsamer ist als Ihr C-Code. Sie sollten praktisch identisch sein. Sie haben jedoch nicht "std :: ios :: sync_with_stdio (false)" angegeben, was möglicherweise einen Teil des Unterschieds ausmacht. und die Schleife könnte durch Verwendung von Iteratoren vereinfacht werden, obwohl ich überrascht wäre, wenn dies die Geschwindigkeit bei einer Einstellung über "-O1" beeinflusst.
@KonradRudolph `sync_with_stdio` hilft wenig beim Lesen aus einer normalen Datei. Ich kenne das Innere von Seqan nicht, daher kann ich nicht erklären, warum es langsamer ist. Trotzdem würde ich mir keine Sorgen um die Leistung von Seqan machen. Die Gzip-Dekomprimierung ist viel langsamer als das Parsen von Fastq. @DevonRyan hinzugefügt.
Der Unterschied liegt in readfq. Das Biest ist * schnell *. Ich weiß, ich habe versucht, Hengs Code in C ++ durch etwas weniger Hässliches zu ersetzen, und konnte ihn bei einfacher Leseleistung nicht annähern. Das Beste, was ich in sauberem C ++ finden konnte, war ein Geschwindigkeitsfaktor von 2x langsamer.
Sie müssen `-lz` zu Ihren Kompilierungsflags hinzufügen (zumindest für das aktuelle Ubuntu LTS)` gcc -O2 -o count-N this-prog.c -lz`
Basieren Ihre Zahlen darauf, den Cache zwischen den Samples zu leeren? Andernfalls sind diese Zahlen mit ziemlicher Sicherheit falsch.
#2
+18
Devon Ryan
2017-06-02 03:20:21 UTC
view on stackexchange narkive permalink

5 Stunden und keine Benchmarks veröffentlicht? Ich bin zutiefst enttäuscht.

Ich werde den Vergleich auf Fasta-Dateien beschränken, da Fastq am Ende derselbe sein wird. Bisher sind die Konkurrenten:

  1. R mit dem ShortRead -Paket (auch wenn es nicht die schnellste ist, sicherlich eine sehr bequeme Methode).
  2. Eine Pipeline von grep -v "^ >" | tr-cd A | wc -c
  3. Eine Pipeline von grep -v "^ >" | grep -io A | wc -l
  4. Eine Pipeline von seqtk comp | awk
  5. Etwas Benutzerdefiniertes in C / C ++
  6. ol>

    Der R-Code kann aus irgendeinem Grund N Datensätze nicht zählen, wenn ich es versuche, also habe ich gezählt Ein für alles.

    Ich werde alle Python-Lösungen ausschließen, da es kein denkbares Universum gibt, in dem Python schnell ist. Für die C / C ++ - Lösung habe ich nur das verwendet, was @ user172818 gepostet hat (alles, was ich versucht habe zu schreiben, hat die gleiche Zeit in Anspruch genommen).

    100-mal wiederholen und den Durchschnitt nehmen (ja, man sollte den Median nehmen):

    1. 1,7 Sekunden
    2. 0,65 Sekunden
    3. 15 Sekunden
    4. 1,2 Sekunden
    5. 0,48 Sekunden
    6. ol>

      Es überrascht nicht, dass alles in Straight C oder C ++ gewinnen wird. grep mit tr ist überraschend gut, was mich überrascht, da ich, obwohl grep sehr, sehr optimiert ist, immer noch mehr Overhead erwartet habe. Das Piping zu grep -io ist ein Leistungskiller. Die R-Lösung ist angesichts des typischen R-Overheads überraschend gut.

      Update1: Wie in den Kommentaren

      1. Summe (x.count ('A') für x in open ('seqs.fa', 'r'), wenn x [0]! = '>') in Python
      2. ol >

        Dies ergibt eine Zeit von

        1. 1,6 Sekunden (ich bin jetzt auf der Arbeit, daher habe ich die Zeit so gut wie möglich angepasst, um die Unterschiede zu berücksichtigen Computer)
        2. ol>

          Update2: Weitere Benchmarks aus den Kommentaren

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

            Diese Fließzeiten betragen:

            1. 5,2 Sekunden
            2. 1,18 Sekunden
            3. ol>

              Die awk-Version ist der beste Hack, den ich mir vorstellen kann. Es gibt wahrscheinlich bessere Möglichkeiten.

              Update 3 : Richtig, in Bezug auf Fastq-Dateien führe ich dies auf einem aus 3,3 GB komprimierte Datei mit 10 Wiederholungen (die Ausführung dauert etwas), daher beschränke ich Tests zunächst nicht auf Befehle, die ich trivial ändern kann, um komprimierte Dateien zu verarbeiten (unkomprimierte Fastq-Dateien sind schließlich ein Gräuel). .

              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. Das C-Beispiel von user1728181 (der Vergleich mit C ++ ist dort bereits vorhanden, muss also nicht eingeschlossen werden).
              5. bioawk -c fastx '{n + = gsub (/ A /, "", $ seq)} END {print n}' foo.fastq.gz
              6. ol>

                Die durchschnittlichen Zeiten sind:

                1. 1 Minute 44 Sekunden
                2. 2 Minuten 6 Sekunden
                3. 1 Minute 52 Sekunden
                4. 1 Minute 15 Sekunden
                5. 3 Minuten 8 Sekunden
                6. ol>

                  Hinweis: Diese verarbeiten häufig keine Fastq-Dateien mit Einträgen, die mehr als 4 Zeilen umfassen. Solche Dateien sind jedoch die Dodo-Vögel der Bioinformatik, daher ist das praktisch nicht wichtig. Außerdem konnte ich ShortRead nicht dazu bringen, mit der komprimierten fastq-Datei zu arbeiten. Ich vermute, es gibt eine Möglichkeit, das zu beheben.

Netter Benchmark, wie groß ist Ihre Test-Fasta-Datei? Warum nicht Python hinzufügen, wenn Sie Zeit haben? Ja, es wird nicht am schnellsten sein, aber zum Vergleich interessant. Ich würde versuchen: `sum (x.count ('A') für x in open ('seqs.fa') , 'r') wenn x [0]! = '>') `, wird` BioPython` dafür wahrscheinlich langsamer sein
@Chris_Rands Ich habe dm6 (unkomprimiert) verwendet, es ist also nicht besonders groß. Ich dachte mir, dass Python so offensichtlich langsamer sein würde, dass ich mich nicht darum kümmern wollte, aber ich nehme an, es würde nicht schaden, etwas hinzuzufügen.
@Chris_Rands Ich habe mit einem Python-Benchmark aktualisiert. Es funktioniert genauso gut wie R, was besser ist, als ich gedacht hätte.
Vielen Dank! Ich bin nicht so überrascht, da Pythons "str.count" in C implementiert ist. Ich denke auch, dass es eine schnelle "awk "-Lösung geben sollte, aber ich kann es nicht herausfinden
Perl: `` `perl -ne 'if (! / ^> /) {$ Count + = tr / A //} END {print" $ count \ n "}' dm6.fa```
Okay, jetzt, da die Frage aktualisiert wurde, können Sie diese so ändern, dass sie für den allgemeinen FASTQ-Fall funktionieren (mehrzeilige Sequenz + Qualitätszeichenfolge, wobei '@' möglicherweise in den Qualitätszeilen angezeigt wird).
Mein awk-fu reicht aus, um ein funktionierendes Beispiel zu finden, aber nicht gut genug, um eines zu finden, das nicht schmerzhaft langsam ist.
@gringer Da bin ich heute produktiv ...
Können Sie mein `bioawk -c fastx '{n + = gsub (/ A /," ", $ seq)} END {print n}'` ausprobieren? Ich bin gespannt, wie es mit den anderen Lösungen für Ihre Daten verglichen wird.
@bli: Ja, ich habe momentan 5 andere am Laufen ...
Der R-Code kann Ns zählen, sobald Sie ihn richtig geschrieben haben. ;-) (Mein Fehler)
Die awk- und perl-Versionen liefern in meiner Test-Fasta-Datei nicht die gleichen Ergebnisse: Die awk-Version enthält im Genom von C. elegans ein "A" mehr als die perl-Version. Ich denke, es hängt davon ab, ob eine der Sequenzen mit "A" beginnt oder endet.
Ich habe einen Fastq-Vergleich mit einigen Methoden (einschließlich Bioawk) hinzugefügt.
@bli Der awk-Code sollte nur ein "-1" enthalten. Unter den vielen Gründen, stattdessen "Bioawk" zu verwenden!
Ich denke, der Zcat ist mit Bioawk nicht notwendig. Verbessert es die Geschwindigkeit?
@bli Es wird nicht benötigt und das ist ein Fehler von meiner Seite (ich habe zcat im Benchmark nicht wirklich verwendet). Danke, dass du das verstanden hast!
Warum zählst du die "A" und nicht die "N"? Oder fehlt mir etwas?
Da die ursprüngliche ShortRead-Variante N nicht verarbeiten konnte und die Leistungsvergleiche unabhängig davon gleich sind. "Der R-Code kann aus irgendeinem Grund nicht N Datensätze zählen, wenn ich es versuche, also habe ich A für alles gezählt."
Ja, tut mir leid, das habe ich verpasst, als ich die Antwort zum ersten Mal gelesen habe.
Sie sagen: "100-mal wiederholen und den Durchschnitt nehmen (ja, man sollte den Median nehmen):" Aber der Durchschnitt ist der Mittelwert, nicht der Median - also welchen haben Sie verwendet? Ich schätze, dass der Median weniger anfällig für Verzerrungen ist, aber in diesem Fall ist die Laufzeit sicherlich normal verteilt (in diesem Fall ist der Mittelwert in Ordnung), es sei denn, Sie arbeiten auf einem System mit sehr hoher augenblicklicher CPU- oder E / A-Last von anderen Benutzern, was zu Konflikten führt .
Der Mittelwert, nicht der Median, also die Einschränkung. Ja, die Laufzeit sollte normal verteilt sein, aber gelegentlich treten seltsame E / A-Spitzen und dergleichen auf.
Haben Sie den Systemcache zwischen all diesen Tests geleert?
Nein, in der Praxis befanden sich die Dateien für alle Läufe im Cache.
Sie sollten den Cache zwischen den Tests leeren. Ihre Zahlen sind ansonsten irreführend (nun ja, wirklich falsch), unabhängig davon, ob Sie aus den Testergebnissen einen Mittelwert oder einen Median ziehen.
#3
+9
Karel Brinda
2017-06-01 22:22:20 UTC
view on stackexchange narkive permalink

Ich denke, das sollte ziemlich schnell gehen:

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  

In dieser Antwort auf SO erfahren Sie, wie Sie Zeichen in BASH mit verschiedenen Ansätzen zählen.

Ihre FASTQ-Antwort hat leider grundsätzlich das gleiche Problem wie die von Iakov. Das heißt, FASTQ-Datensätze haben eine variable Anzahl von Zeilen.
Ich weiß, dass der ursprüngliche Standard die FASTQ-Linienaufteilung zulässt. Haben Sie jemals eine solche Datei in der Praxis gesehen?
Ich habe es nicht getan, aber ich habe selbst nie mit Sanger-Sequenzierungs- / PacBio-Daten gearbeitet. Ich denke, Sie könnten ihnen dort begegnen. Für kurze Lesedaten denke ich, dass die Annahme „Datensatz = 4 Zeilen“ sicher ist.
Lang gelesene Daten verwenden kein Fastq. Derzeit verwendet PacBio bam (ja, wirklich bam) und ONT fast5 (h5-Archiv). Ich denke, es ist sehr vernünftig, 4 Zeilen Fastq anzunehmen.
#4
+8
Iakov Davydov
2017-06-01 22:00:26 UTC
view on stackexchange narkive permalink

FASTQ

Wie bereits erwähnt, kann fastq kompliziert sein. In einem einfachen Fall, wenn Sie vier Zeilen pro Datensatz haben, ist eine mögliche Lösung in bash:

  sed -n '2 ~ 4p' seqs.fastq | grep -io N | wc -l  
  • sed -n '2 ~ 4p' druckt jede vierte Zeile
  • grep -o N gibt für jedes übereinstimmende Symbol eine Zeile mit N aus.
  • wc -l zählt die Zeilen

Ich vermute, dass dieser Python-Ansatz schneller funktioniert:

  cat seqs.fastq | python3 -c "sys importieren; print (sum (line.upper (). count ('N') für Zeile in sys.stdin))"  

FASTA

Coreutils:

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

Python:

  cat seqs.fasta | python3 -c "sys importieren; print (sum (line.upper (). count ('N') für Zeile in sys.stdin, wenn nicht line.startswith ('>'))"  
Jetzt zählen Sie Zeilen mit Ns, nicht Ns. Außerdem können FASTQ-Datensätze mehr als vier Zeilen umfassen. :-(
@KonradRudolph mit -o grep gibt einzelne Symbole aus, ich werde überlegen, wie man Fastq-Code korrigiert
FASTQ ist schwierig, da `@` auch die Qualitätszeichenfolge starten kann. Siehe http://chat.stackexchange.com/transcript/message/37809200#37809200 und die folgenden Zeilen. Als Referenz finden Sie hier ein Perl-Skript, das FASTQ analysiert und die Partnerinformationen von der Kennung abschneidet. Es kann wahrscheinlich als Ausgangspunkt dienen; Ich denke, dies kommt dem minimalen Code sehr nahe: https://gist.github.com/klmr/81a2b86cd93c706d28611f40bdfe2702
@KonradRudolph, Ich habe die Antwort aktualisiert, um anzuzeigen, dass dies nur im einfachen Fall funktioniert
#5
+6
bli
2017-06-02 14:28:46 UTC
view on stackexchange narkive permalink

Verwenden von Bioawk

Mit Bioawk (Zählen von "A" im Genom von C. elegans , da anscheinend kein "N" vorhanden ist diese Datei) auf meinem Computer:

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

bioawk ist eine Erweiterung von awk mit praktischen Analyseoptionen. Zum Beispiel stellt -c fastx die Sequenz als $ seq im Fasta- und Fastq-Format zur Verfügung ( einschließlich komprimierter Versionen ), sodass der obige Befehl sollte auch das Fastq-Format robust handhaben .

Der Befehl gsub ist eine normale awk-Funktion. Es gibt die Anzahl der ersetzten Zeichen zurück (von https://unix.stackexchange.com/a/169575/55127). Ich begrüße Kommentare darüber, wie man in awk weniger hackig zählt.

In diesem Fasta-Beispiel ist es fast dreimal langsamer als das grep , tr , wc Pipeline:

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

(Ich habe die Timings wiederholt und die obigen Werte scheinen repräsentativ zu sein)

Verwenden von Python

readfq

Ich habe versucht, die Python-Readfq aus dem in dieser Antwort genannten Repository abzurufen: https://bioinformatics.stackexchange.com/a/365/292

  $ time python -c "von sys import stdin; von readfq import readfq; Drucksumme ((seq.count ('A') für _, seq, _ in readfq (stdin)))" < Genom .fa32371810real 0m0.976suser 0m0.876ssys 0m0.100s  

"pure python"

Vergleich mit etwas, das aus dieser Antwort angepasst wurde: https: // bioinformatics. stackexchange.com/a/363/292

  $ time python -c "von sys import stdin; Drucksumme ((line.count ('A') für Zeile in stdin) wenn nicht line.startswith ('>'))) "< Genom.fa32371810real 0m1.143suser 0m1.100s
sys 0m0.040s  

(Es ist langsamer mit Python 3.6 als mit Python 2.7 auf meinem Computer)

pyGATB

Ich habe gerade gelernt (08 / 06/2017), dass GATB einen Fasta / Fastq-Parser enthält und kürzlich eine Python-API veröffentlicht hat. Ich habe gestern versucht, damit eine andere Antwort auf die vorliegende Frage zu testen, und habe einen Fehler gefunden. Dieser Fehler ist jetzt behoben. Hier ist eine pyGATB -basierte Antwort:

  $ time python3 -c "von der gatb import Bank; print (sum ((seq. sequence.count (b "A") für seq in Bank ("Genom.fa"))) 32371810real 0m0.663suser 0m0.568ssys 0m0.092s  

( Sie können auch sequence.decode ("utf-8"). Count ("A") ausführen, dies scheint jedoch etwas langsamer zu sein.)

Obwohl ich hier python3.6 verwendet habe (pyGATB scheint nur python3 zu sein), ist dies schneller als die beiden anderen Python-Ansätze (für die die angegebenen Timings mit Python 2.7 erhalten werden). Dies ist sogar fast so schnell wie die Pipeline grep , tr , wc .

Biopython

Und um noch mehr Vergleiche zu haben, hier eine Lösung mit SeqIO.parse von Biopython (mit python2.7):

  $ time python -c "von Bioimport SeqIO; print (sum ((rec.seq.count ("A") für rec in SeqIO.parse ("Genom.fa", "fasta"))) 32371810real 0m1.632suser 0m1.532ssys 0m0.096s  

Dies ist etwas langsamer als die "pure python" -Lösung, aber möglicherweise robuster.

Es scheint eine leichte Verbesserung mit zu geben @ peterjcs Vorschlag, die niedrigere Ebene zu verwenden SimpleFastaParser :

  $ time python -c "von Bio.SeqIO.FastaIO importieren SimpleFastaParser; print (sum (seq.count ( 'A') für Titel, seq in SimpleFastaParser (offen ('Genom.fa'))) 32371810real 0m1.618suser 0m1.500ssys 0m0.116s  

(Ich habe eine Reihe von Timings und versuchte, eine zu nehmen, die repräsentativ schien, aber es gibt viele Überlappungen wi th das Timing des übergeordneten Parsers.)

scikit-bio

Ich habe heute (23.06.2017) über scikit-bio gelesen, das über einen Sequenzleser verfügt.

  $ time python3 -c "import skbio; print (sum (seq.count ('A') für seq in skbio.io.read ('Genom.fa', format = 'fasta', verify = False)) "32371810real 0m3.643suser 0m3.440ssys 0m1.228s  

pyfastx

Ich habe heute (27.08.2020) über pyfastx gelesen, das über einen Sequenzleser verfügt.

  $ time python3 -c "aus pyfastx import Fasta; print (Summe (seq.count ('A') für (_, seq) in Fasta ('Genom.fa', build_index = False))" 32371810real 0m0.781suser 0m0.740ssys 0m0.040s  

(Das Timing erfolgt mit Python 3.6 auf demselben Computer wie die vorherigen Timings von vor 3 Jahren, daher sollte es vergleichbar sein.)

fastx_read von mappy

Mir ist klar (27.08.2020), dass ich keinen Benchmark für fastx_read von hinzugefügt habe mappy, das ich seit geraumer Zeit regelmäßig benutze:

  $ time python3 -c "von mappy import fastx_read; print (sum (seq.count ('A') für (_, seq, _) in fastx_read ('Genom.fa')) "32371810real 0m0.731suser 0m0.648ssys 0m0.080s  

Benchmarks für eine fastq.gz-Datei

Ich habe einige der oben genannten Lösungen (oder Anpassungen davon) getestet und "N" in derselben Datei gezählt, die hier verwendet wurde: https: // bioinformatics .stackexchange.com / a / 400/292

bioawk

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

pigz + readfq Python-Modul

readfq beschwert sich nicht und ist sehr schnell, wenn ich das komprimierte Fastq direkt übergebe, aber etwas Falsches zurückgibt. Vergessen Sie also nicht, die Dekomprimierung manuell durchzuführen.

Hier habe ich es mit pigz versucht:

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

Der Computer verfügt über 16 Kerne, aber ich vermute, dass der begrenzende Faktor für pigz das Lesen von der Festplatte ist: Die Prozessoren sind weit davon entfernt, die volle Geschwindigkeit zu erreichen.

Und mit gzip :

  $ time gzip -dc SRR077487_2.filt.fastq.gz | python -c "von sys import stdin; von readfq import readfq; Drucksumme ((seq.count ('N') für _, seq, _ in readfq (stdin)))" 306072real 0m49.448suser 1m31.984ssys 0m2.312s  

Hier verwendeten gzip und python beide einen vollen Prozessor. Das Gleichgewicht der Ressourcennutzung war etwas besser. Ich arbeite auf einem Desktop-Computer. Ich nehme an, ein Computerserver würde pigz besser nutzen.

pyGATB

  $ time python3 -c "von der gatb import Bank; print ( sum ((seq.sequence.count (b "N") für seq in Bank ("SRR077487_2.filt.fastq.gz"))) 306072real 0m46.784suser 0m46.404ssys 0m0.296s  

Biopython

  $ time python -c "vom gzip-Import als gzopen öffnen; vom Bio.SeqIO.QualityIO-Import FastqGeneralIterator; print (sum (seq.count (' N ') für (_, seq, _) in FastqGeneralIterator (gzopen (' SRR077487_2.filt.fastq.gz '))) 306072real 3m18.103suser 3m17.676ssys 0m0.428s  

pyfastx (hinzugefügt am 27.08.2020)

  $ time python3 -c "aus pyfastx importiert Fastq; print (sum (seq.count ('N') für (_, seq, _) ) in Fastq ('SRR077487_2.filt.fastq.gz', build_index = False)) 306072real 0m51.140suser 0m50.784ssys 0m0.352s  

fastx_read von mappy (hinzugefügt am 27.08.2020)

  $ time python3 -c "von mappy import fastx_read; pri nt (sum (seq.count ('N') für (_, seq, _) in fastx_read ('SRR077487_2.filt.fastq.gz')) 306072real 0m52.336suser 0m52.024ssys 0m0.300s  

Schlussfolgerung

pyGATB und bioawk behandeln beide transparent Komprimierungs- (gzipped oder nicht) und Formatunterschiede (fasta oder fastq). pyGATB ist ein ziemlich neues Tool, scheint aber im Vergleich zu den anderen von mir getesteten Python-Modulen effizienter zu sein.

(Update 27/08/2020) pyfastx und mappy scheinen ebenfalls recht praktisch zu sein und sind nur geringfügig langsamer als pyGATB .

Könnten Sie die FASTA-Parser-Nummern für Biopython-Zeichenfolgen auf niedriger Ebene hinzufügen? `` time python -c "von Bio.SeqIO.FastaIO importieren SimpleFastaParser; print (sum (seq.count ('A') für Titel, seq in SimpleFastaParser (open ('Genom.fa')))" ``
Warum zählen alle Ihre Beispiele den Buchstaben A anstelle von N (und vielleicht n) gemäß der Frage?
Ich zähle den Buchstaben "A", weil die erste große Datei, an die ich bei meinen Tests gedacht habe, kein "N" hatte.
Ich würde denken, dass alle Python-Beispiele in dieser kleinen Beispieldatei aufgrund des Overheads beim Laden von Python und beim Importieren von Bibliotheken schlecht abschneiden. Wenn der Benutzer jedoch nur ein Genom zählen möchte, ist dies dennoch ein fairer Test.
#6
+6
Matt Bashton
2017-06-03 01:28:27 UTC
view on stackexchange narkive permalink

Die Dekomprimierung von gezipptem FASTQ ist das Hauptproblem.

Wenn wir realisiertes komprimiertes FASTQ (was, wie das OP vorgeschlagen hat, von Vorteil wäre) anstelle von trivialem FASTA als Ausgangspunkt nehmen, dann ist das eigentliche Problem tatsächlich dekomprimiert Die Datei zählt nicht die Ns und in diesem Fall ist das C-Programm count-N nicht mehr die schnellste Lösung.

Außerdem wäre es gut, eine bestimmte Datei für das Benchmarking zu verwenden hat tatsächlich Ns, weil Sie einige interessante Unterschiede in der Ausführungszeit erhalten, wenn einige Methoden die häufiger auftretenden As anstelle von Ns zählen.

Eine solche Datei ist:

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

Es lohnt sich auch, die verschiedenen zu überprüfen Lösungen geben die richtige Antwort zurück. Die obige Datei sollte 306072 Ns enthalten.

Beachten Sie als nächstes, dass die Dekomprimierung dieser Datei, die nach / dev / null umgeleitet wird, mit zcat und gzip (die beide auf meinem System gzip 1.6 sind) als eine parallele Implementierung von gzip wie Mark Adlers pigz , die anscheinend 4 Threads zur Dekomprimierung verwendet. Alle Timings stellen einen Durchschnitt von 10 Läufen dar, die die reale Zeit (Wanduhrzeit) melden.

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

Zwischen beiden besteht ein Unterschied von ~ 11,7 Sekunden. Als nächstes, wenn ich dann versuche, einen Einzeiler zu bewerten, der mit FASTQ arbeitet und die richtige Antwort gibt (Hinweis: Ich habe noch keinen FASTQ gefunden, der nicht aus vier Zeilen besteht, und ernsthaft, wer diese Dateien generiert!)

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

Wie Sie sehen können, erhöht die Zählung die Gesamtlaufzeit um ~ 5,8 Sekunden gegenüber der auf pigz basierenden Dekomprimierung. Zusätzlich ist dieses Zeitdelta höher, wenn gzip ~ 6,7 Sekunden über der gzip-Dekomprimierung allein verwendet wird.

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

Die auf pigz awk tr wc basierende Lösung ist jedoch ~ 4,5 Sekunden schneller als die count-N basierte C-Code-Lösung:

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

Dieser Unterschied scheint robust zu sein so oft ich will wieder laufen. Ich gehe davon aus, dass wenn Sie pthread in der C-basierten Lösung verwenden oder ändern könnten, um den Standard von pigz zu entfernen, dies auch zu einer Leistungssteigerung führen würde.

Benchmarking einer anderen Alternative pigz grep Die Variante scheint mehr oder weniger dieselbe Zeit zu dauern wie die auf tr basierende Variante:

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

Beachten Sie, dass die oben diskutierte seqtk -basierte Lösung merklich langsamer ist,

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

Es ist jedoch erwähnenswert, dass seqtk comp etwas mehr leistet als die anderen Lösungen .

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

Ehrlich gesagt besteht der einfachste Weg (insbesondere für FASTQ) wahrscheinlich darin, eine dedizierte Parsing-Bibliothek wie R / Bioconductor zu verwenden:

  unterdrückenMessagen (Bibliothek (ShortRead)) seq = readFasta (commandArgs (TRUE) [1]) # oder readFastqcat (colSums (alphabetFrequency (sread (seq)) [, 'N', drop = FALSE]), '\ n')  

Dies kann sein nicht die schnellste sein, aber es ist ziemlich schnell , da die relevanten Methoden stark optimiert sind. Der größte Aufwand entsteht wahrscheinlich durch das Laden von ShortRead , was ein nicht zu rechtfertigender Slog ist.

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

Wenn Sie nach roher Geschwindigkeit suchen, müssen Sie wahrscheinlich ein eigenes kleines C / C ++ - Programm schreiben.

Glücklicherweise hat das Schlimmste (ein schneller und zuverlässiger Parser) bereits angegangen: Das readfq von Heng Li ist wahrscheinlich der schnellste FASTA / FASTQ-Parser.

Und es ist einfach zu bedienen, das Beispiel auf GitHub kann einfach erweitert werden, um was zu tun du brauchst. Fügen Sie einfach einen Parameter-Parser hinzu (für den Dateinamen, den Sie analysieren möchten), programmieren Sie eine einfache N-Zähler-Schleife und lassen Sie die Ergebnisse auf stdout drucken. Fertig.

Wenn Sie sonst nach Einfachheit suchen, lesen Sie Konrads Antwort auf eine R-basierte Lösung. Oder ein sehr einfaches Skript mit Biopython.

Gott, das ist absolut grausamer Code. ☹️ Ich denke, dass [SeqAn] (http://docs.seqan.de/seqan/master/?p=SeqFileIn) (C ++) für rohe Geschwindigkeit wahrscheinlich gleichwertig ist und eine viel höhere Qualität aufweist.
Es ist nicht das schönste. Andererseits macht es das, was es tun muss und kommt nicht mit einer ganzen Bibliotheksabhängigkeit. Haben Sie einen Benchmark mit SeqAn & readfq?
Ich bin gespannt, wie readfq ans SeqAn mit dem Fasta / Fastq-Parser von GATB verglichen werden kann: http://gatb-core.gforge.inria.fr/doc/api/classgatb_1_1core_1_1bank_1_1impl_1_1BankFasta.html#detailsIch kann nicht wirklich codieren in C / C ++ allerdings.
#9
+3
Matt Shirley
2017-06-15 20:47:20 UTC
view on stackexchange narkive permalink

Für FASTA-Dateien habe ich in pyfaidx v0.4.9.1 eine relativ effiziente Methode implementiert. Dieser Beitrag hat mir klar gemacht, dass mein vorheriger Code ziemlich langsam und einfach zu ersetzen war:

  $ pip install pyfaidx $ time faidx -i Nucleotid ~ / 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 real 2m28.251suser 2m15.548ssys 0m9.951s  
#10
  0
Edward Kirton
2018-03-01 11:55:22 UTC
view on stackexchange narkive permalink

Hier sind Perl-1-Liner für Fasta- und Fastq-Dateien, die schnell genug sind, wenn Sie nur Summen möchten:

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

Für jede Zeile einer FASTA-Datei ist es eine Sequenzzeile, wenn nicht der Header (/ ^> /). Ersetzen Sie Ns, wodurch die Anzahl der vorgenommenen Ersetzungen zurückgegeben wird, die zur Gesamtsumme addiert werden.

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

Für jede Sequenzzeile einer FASTQ-Datei, die mit / ^ [ATCGN] + $ / i (dh allen NTs) übereinstimmt, Zählen Sie dann die auf die gleiche Weise vorgenommenen Ersetzungen.

Sie können pigz anstelle von zcat (oder gunzip -c) verwenden, wenn pigz installiert ist.

Eine unkomprimierte fastq-Datei mit 14 MB 250 bp Illumina Lesevorgänge in 40 Sekunden abgeschlossen.

Wenn Sie erklären könnten, was jede Perl-Zeile tut, würde dies die Antwort verbessern und insbesondere die Unterschiede zwischen beiden Ansätzen kommentieren


Diese Fragen und Antworten wurden automatisch aus der englischen Sprache übersetzt.Der ursprüngliche Inhalt ist auf stackexchange verfügbar. Wir danken ihm für die cc by-sa 3.0-Lizenz, unter der er vertrieben wird.
Loading...